CICR and CVODE issues with Kinetic Scheme

[jjnaylor]

I'm attempting to model ER/cytosol calcium oscillations with a synaptic calcium current. The ER oscillation bit is a reduction of the Keizer-Levine model. My state variable 'l' represents the fraction of non-inactivated IP3 channels.

Code: Select all

STATE {
	cai 	(mM)
	caE		(mM)
	l		(1)
	q		(1)
}
BREAKPOINT {		
		SOLVE state METHOD sparse
}
KINETIC state {
	rates()
	COMPARTMENT PI*diam*diam/4 {cai}
	COMPARTMENT sigma*PI*diam*diam/4 {caE}
	~ caE <-> cai (fi*Jin, fi*(Jin+Jout))
	~ l <-> q (loff, lon)
:	~cai<< (-ica*PI*diam/(2*FARADAY))       Include synaptic input after debug ER oscil.
	CONSERVE cai + caE = caT*PI*diam*diam/4
	CONSERVE l + q = 1
}

PROCEDURE rates() {
	Jin  = (Lip3 + Pip3 * ((ip3*cai*l)/( (ip3 + ki) * (cai+ka) ))^3)
	Jout = serca(cai)
	loff = A*cai
	lon = A*Kd
}
FUNCTION serca(x (mM)) {
	serca	=	vser * ( x / (x^2 + kser^2) )		
}	

While this yields ER emptying and refilling periodically with a fixed timestep, it empties and refills only once with a variable timestep, and depending on certain values of A can become unstable. Additionally, when I compile, I see a statement that says NEURON's CVODE doesn't recognize CONSERVE statements. Although I don't intuitively see how, should I try to rework this scheme to not have reaction rates dependent on state variables (the NEURON book says that could complicate things and make the variable step integrator compensate by reducing dt and thus slowing simulations)?

Regarding units/dimensions: should the COMPARTMENT statement take a volume or volume per length term? do I want Jin/Jout to be in units of inverse time or volume per unit time (since the compartment statement will divide the rate by the compartment volume)? Thanks for your help!

[ted]

The first question is whether the "reduction of the Keizer-Levine model" that you want to implement captures the dynamics you are interested in, i.e. has the necessary nonlinear dynamical analysis been performed and the results found to be appropriate?

The second question is whether or not the NMODL code is a correct implementation of the reduced K-L model. Not knowing the reduced K-L model's equations (the original had 4 states), I can't say more.

While this yields ER emptying and refilling periodically with a fixed timestep, it empties and refills only once with a variable timestep, and depending on certain values of A can become unstable.

Not knowing anything about the reduced K-L model I can't tell whether this is to be expected or whether it indicates that {the "reduction" is incorrect} OR {the "reduction" is correct, and this is just how these equations behave when solved with fixed dt or adaptive integration}.

Additionally, when I compile, I see a statement that says NEURON's CVODE doesn't recognize CONSERVE statements.

CONSERVE is merely a hint to the NMODL translator that the referenced STATEs are not independent i.e. one of them can be calculated from the others by simple algebraic manipulation.

Although I don't intuitively see how, should I try to rework this scheme to not have reaction rates dependent on state variables (the NEURON book says that could complicate things and make the variable step integrator compensate by reducing dt and thus slowing simulations)?

Before trying to recast the equations in a different form, ask whether slowing of the adaptive integrator is the problem--doesn't seem so. How are you using adaptive integration--entirely from hoc, in which case I would ask how you selected the values for the error tolerances are appropriate for the model's state variables, or are you using an "analysis run" launched from the VariableStepControl's "Atol Scale Tool," then clicking on that tool's Rescale button?

Regarding units/dimensions: should the COMPARTMENT statement take a volume or volume per length term? do I want Jin/Jout to be in units of inverse time or volume per unit time (since the compartment statement will divide the rate by the compartment volume)?

COMPARTMENT statements require (volume/unit length) for STATEs that have units of (/volume) e.g. stuff in solution, and (area/unit length) for STATEs that have units of (/area) e.g. pumps and other stuff that is bound to membranes. For anything other than typical v-gated HH-style channel models, it is usually best to employ modlunit to detect instances of unit inconsistency, because it is very fastidious and provides useful hints for how to resolve such inconsistencies. modlunit also detects many syntax errors.

If you're nearing your wit's end with NMODL (very easy to do) and would like additional help with this task, zip up just enough source code (hoc, ses, mod files) to run a test simulation, and email to me
ted dot carnevale at yale dot edu
and I'll tell you what I find. Also, please point me to the source of the reduced equations so I can verify the translation from published source to NMODL.