Solving 8 interdependent differential equations

NMODL and the Channel Builder.
Post Reply
bdut
Posts: 5
Joined: Thu Sep 05, 2019 11:53 am

Solving 8 interdependent differential equations

Post by bdut » Thu Sep 19, 2019 12:53 am

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.

ted
Site Admin
Posts: 5632
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Solving 8 interdependent differential equations

Post by ted » Thu Sep 19, 2019 10:09 am

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?

bdut
Posts: 5
Joined: Thu Sep 05, 2019 11:53 am

Re: Solving 8 interdependent differential equations

Post by bdut » Thu Sep 19, 2019 10:37 am

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.

ted
Site Admin
Posts: 5632
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Solving 8 interdependent differential equations

Post by ted » Thu Sep 19, 2019 10:58 am

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.

bdut
Posts: 5
Joined: Thu Sep 05, 2019 11:53 am

Re: Solving 8 interdependent differential equations

Post by bdut » Fri Sep 20, 2019 1:20 am

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

ted
Site Admin
Posts: 5632
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Solving 8 interdependent differential equations

Post by ted » Fri Sep 20, 2019 10:47 am

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.

bdut
Posts: 5
Joined: Thu Sep 05, 2019 11:53 am

Re: Solving 8 interdependent differential equations

Post by bdut » Fri Sep 20, 2019 11:36 am

I have mailed you the zip of the code. It can be used with a single compartment model. Thank you for your timely help.

Post Reply