Page 1 of 1

Solving 8 interdependent differential equations

Posted: Thu Sep 19, 2019 12:53 am
by bdut
In my .mod file, I have 8 DEs to be solved, and all are interdependent. DERIVATIVE block is as follows:

Code: Select all

DERIVATIVE blrrate {
	calcParameters(blr, aG, cAMP, aCaMK, Vcilia, CaCaM, Ca, IX)
	blr' = k1 * od * (Rtot - blr) - r1 * blr
	aG' = kG * (Gtot - aG) - rG
	cAMP' = syncAMP - pd * cAMP
	Ca' = inf * Icng - JNCX - (cc1 - cc2 * CaCaM)
	CaCaM' = cc1lin * Ca - cc2 * CaCaM
	aCaMK' = ck1lin * CaCaM - ck2 * aCaMK
	IX' = cx1lin * Ca - cx2 * IX
	Vcilia' = (Icng+IClCa+IL)/ CapCilia
}
calcParameters(...) is a procedure that evaluates the formulas.

When I debugged with printf statement, only the first 2 equations are getting solved. The rest of the variables do not have values.
t: 4.8, blr: 7.25551e-011, aG: 7.26423e-007, cAMP: -1.#IND, Ca: -1.#IND, CaCaM: -1.#IND, aCaMK: -1.#IND, IX: -1.#IND, Vcilia: -1.#IND
t: 4.825, blr: 6.10429e-011, aG: 6.52792e-007, cAMP: -1.#IND, Ca: -1.#IND, CaCaM: -1.#IND, aCaMK: -1.#IND, IX: -1.#IND, Vcilia: -1.#IND
t: 4.85, blr: 5.13573e-011, aG: 5.86624e-007, cAMP: -1.#IND, Ca: -1.#IND, CaCaM: -1.#IND, aCaMK: -1.#IND, IX: -1.#IND, Vcilia: -1.#IND
t: 4.875, blr: 4.32085e-011, aG: 5.27162e-007, cAMP: -1.#IND, Ca: -1.#IND, CaCaM: -1.#IND, aCaMK: -1.#IND, IX: -1.#IND, Vcilia: -1.#IND
t: 4.9, blr: 3.63526e-011, aG: 4.73726e-007, cAMP: -1.#IND, Ca: -1.#IND, CaCaM: -1.#IND, aCaMK: -1.#IND, IX: -1.#IND, Vcilia: -1.#IND

Could you please provide some pointers as to why this is happening? The same equations work fine in Matlab.

Re: Solving 8 interdependent differential equations

Posted: Thu Sep 19, 2019 10:09 am
by ted
Does the STATE block declare all of the variables that appear on the LHS of the ODEs?
What are the contents of the BREAKPOINT block?
Is initialization of the state variables correct?
Is PROCEDURE calcparameters() working properly?

Re: Solving 8 interdependent differential equations

Posted: Thu Sep 19, 2019 10:37 am
by bdut
Thank you Ted for your reply.

STATE block declares all of the variables on LHS of the ODE's.

BREAKPOINT contains below code:

Code: Select all

BREAKPOINT {
	hv1 = 1/(1+(exp(-(t-1)/sharpness)))
	hv2 = 1/(1+(exp(-(t-2)/sharpness)))
	ipulse = hv1 - hv2 
	od = ostim * ipulse
	SOLVE blrrate METHOD cnexp

	printf("t: %g, blr: %g, aG: %g, cAMP: %g, Ca: %g, CaCaM: %g, aCaMK: %g, IX: %g, Vcilia: %g \n", t, blr, aG, cAMP, Ca, CaCaM, aCaMK, IX, Vcilia)	
}
I have initialized the state variables as follows:

Code: Select all

INITIAL {
	blr=1.e-8 
	aG=1.e-8
	cAMP=1.e-8
	Ca=1.e-8
	CaCaM=1.e-8
	aCaMK=1.e-8
	IX=1.e-8
	Vcilia=-69.6650	
}
I will check the calcparameters procedure once if it calculates all the variables properly with a printf() method.

Re: Solving 8 interdependent differential equations

Posted: Thu Sep 19, 2019 10:58 am
by ted
I'd suggest changing
SOLVE blrrate METHOD cnexp
to
SOLVE blrrate METHOD derivimplicit
because many of the state variables have very small values and should not become negative.

Re: Solving 8 interdependent differential equations

Posted: Fri Sep 20, 2019 1:20 am
by bdut
I have changed to derivimplicit, but there is no change.
I printed the values for the PROCEDURE formulas and found they are not being calculated.

kG: -1.#IND, rG: -1.#IND
t: 5, syncAMP: -1.#IND, Icng: -1.#IND, inhcng: -1.#IND, JNCX: -1.#IND, cc1: -1.#IND, ck1: -1.#IND, cx1: -1.#IND, IClCa: -1.#IND, IL: -1.#IND
kG: -1.#IND, rG: -1.#IND

My procedure is as follows:

Code: Select all

PROCEDURE calcParameters(blr, aG, cAMP, aCaMK, Vcilia, CaCaM, Ca, IX) {
	kG = k2 * blr
	rG = r2 * aG
	syncAMP = aG * (smax/(1+(aCaMK/kinh)^ninh))
	Icng = cnmax * (cAMP^n1/(cAMP^n1 + (inhcng * hmc1)^n1)) * (vcng - Vcilia)
	inhcng = 1 + ((inhmax-1) * CaCaM^ninhcng)/(CaCaM^ninhcng + kinhcng^ninhcng)
	JNCX = ef * Ca / (1+(IX/k1)^nI)
	cc1 = cc1lin * Ca
	ck1 = ck1lin * CaCaM
	cx1 = cx1lin * Ca
	IClCa = clmax * (Ca^n2 / (Ca^n2 + hmc2^n2)) * (vcl - Vcilia)
	IL = gl * (vl - Vcilia)
}
A few values of Icng, IClCa , IL, from Matlab proram are
-9.690379e-38 -9.543601e-38 -9.374277e-38
-9.069165e-38 -5.428998e-38 -5.216164e-38
-4.923153e-38 -4.245737e-38 -2.283475e-38
-5.370498e-03 -4.374749e-03 -3.263507e-03
-2.848488e-03 -1.624059e-03 -1.051240e-03
-9.310657e-04 -9.305656e-04 -6.166836e-04

Do these values pose a problem?

Thank you

Re: Solving 8 interdependent differential equations

Posted: Fri Sep 20, 2019 10:47 am
by ted
This is not a problem of "magic numbers." There's a bug, either a subtle syntax error in your NMODL code or in NEURON itself. I'll need to work with your NMODL file in a toy model. If the NMODL file can be used by itself in a single compartment model, that's all I'll need. If it requires other mechanisms or a particular cellular anatomy, I'll need whatever additional code is appropriate. Zip it up and email it to
ted dot carnevale at yale dot edu
and I'll see what I can find out.

Re: Solving 8 interdependent differential equations

Posted: Fri Sep 20, 2019 11:36 am
by bdut
I have mailed you the zip of the code. It can be used with a single compartment model. Thank you for your timely help.