mod file CVODE compatibility

[Brad]

We are looking at the persistent sodium channel model (specified in "nap.mod") published in

Poirazi, P. Brannon, T. & Mel, B.W.
"Arithmetic of Subthreshold Synaptic Summation in a Model CA1 Pyramidal Cell."
Neuron 977-987, February 2003.
(file from the ModelDB is http://tinyurl.com/ctor8)

If one makes some minor modifications (removing explicit references to t and dt), we obtain something like the code snippet as follows:




TITLE Na persistent channel

NEURON {
SUFFIX nap2
USEION na READ ena WRITE ina
RANGE gnabar,vhalf, K

}

UNITS {
(mA) = (milliamp)
(mV) = (millivolt)
}

PARAMETER {
v (mV)
ena = 50 (mV) : Na reversal potential (reset in cell-setup.hoc)
K = 4.5 (1) : slope of steady state variable
gnabar = 0 : initialized conductance
vhalf = -50.4 (mV) : half potential
}

STATE { n }

ASSIGNED {
ina (mA/cm2)
}

BREAKPOINT {
SOLVE states
ina = gnabar*n^3*(v-ena)
}

PROCEDURE states() { : exact when v held constant; integrates over dt step
n = 1 / (1 + (exp(vhalf - v)/K)) : steady state value
}




Putting aside all considerations about the accuracy/correctness of the model itself...

1. When we translate this file, the following line appears among the output

"Notice: This mechanism cannot be used with CVODE"

What exactly about this file is incompatible with CVODE? It appears to conform to the most modern NMDOL idioms (i.e. no explicit references to t or dt).

2. In this environment, how big are the computational savings from using multiple multiplication statements ("n*n*n") versus using exponential syntax ("n^3"). The later seems so much more aesthetically and intellectually pleasing.

Brad

[ted]

STATE { n }
. . .
BREAKPOINT {
SOLVE states
ina = gnabar*n^3*(v-ena)
}

PROCEDURE states() { : exact when v held constant; integrates over dt step
n = 1 / (1 + (exp(vhalf - v)/K)) : steady state value
}

. . .

1. When we translate this file, the following line appears among the output
"Notice: This mechanism cannot be used with CVODE"
What exactly about this file is incompatible with CVODE? It appears to conform to the most modern NMDOL idioms (i.e. no explicit references to t or dt).

Very good question. See
http://www.neuron.yale.edu/neuron/stati ... tionIssues
for a concise discussion of what is needed for CVODE compatibility. Here the
problem is that there's no KINETIC or DERIVATIVE block, so there's no differential
equations to be solved. And the BREAKPOINT block doesn't have a SOLVE
statement, so this NMODL code isn't telling the compiler which integration algorithms
to use.

The clues as to what the equations should look like are in the t & dt dependent
stuff you removed.

2. In this environment, how big are the computational savings from using multiple multiplication statements ("n*n*n") versus using exponential syntax ("n^3"). The later seems so much more aesthetically and intellectually pleasing.

Probably no difference, or at least not one that matters.
BTW I haven't seen tremendous savings from using precomputed tables of rates and
steady-state values, either--some models are maybe 10-20% faster, if that much.
Most users would get a bigger bang from spending a few bucks to increase their
PC's RAM.

[Brad]

Ted, the mechanism doesn't actually have any gating variables that are specified by differential equations, just the one algebraic equation that you see in STATES. So I don't believe I deleted anything of consequence to the actual functioning of the channel model. I removed exactly the following lines:

1. "INDEPENDENT {t FROM 0 TO 1 WITH 100 (ms)}"
2. "dt (ms)" from the PARAMETER block (which is probably just left over when the original author cut and pasted from another file)
3. " VERBATIM
return 0;
ENDVERBATIM " from the STATES block.

Which raises another question: One reference I have ("User Defined Membrane Mechanisms" http://neuron.duke.edu/userman/10/modl.html#VERBATIM) indicates that a PROCEDURE solved by a SOLVE statement may return an error code. Is this still the case with the current version? I don't notice any ill effects if those lines are commented out. [The documentation does say "may" return an error code, so perhaps there are other situations where it is necessary.]

At any rate, I edited a few blocks to read as follows




STATE { m } : dummy gate variable so integrator has something to chew on

ASSIGNED {
ina (mA/cm2)
n
}

INITIAL { : dummy block so the integrator has something to "solve"
m = 0
}

BREAKPOINT {
SOLVE states METHOD cnexp
ina = gnabar*n^3*(v-ena)
}

DERIVATIVE states {
m' = 0
n = 1 / (1 + (exp(vhalf - v)/K)) : steady state value
}




The modification makes the mechanism compatible with CVODE and the tested behavior is identical to the original model.

I am also glad you made a comment on the use of TABLES, which always seemed to me to be of questionable advantage. I was hard-pressed to tell a difference.

***What is the format code for tabs?

[ted]

Sorry, I misinterpreted your code. In too much of a hurry to wrap up
preparations for the NEURON course.

First a word about CVODE and DASPK. These are names for adaptive
integrators that were incorporated into NEURON that have subsequently
been replaced by CVODES and IDA, respectively (parts of the SUNDIALS
package). However, the names CVODE and DASPK are what these
things are called in NEURON for "historical" reasons.

CVODES (and its predecessor CVODE) cannot deal with a system in
which one or more STATEs is described by an algebraic equation.
Your original mechanism has a STATE whose voltage-dependence is
described by an algebraic equation, not a differential equation. That's
why you got the msg that it wasn't compatible with CVODE. However, it
is compatible with IDA (and its predecessor DASPK). I think that
NEURON's adaptive integrator is smart enough to automatically use
DASPK when your original mechanism is present. The way to tell is
to just use your original mechanism in a model, then click on the
Use variable dt button in the VariableTimeStep GUI tool, and finally
click on that tool's Details button, which brings up the Numerical Value
Selection window in which you will see several radio buttons, one of which
is labeled while another is labeled Daspk. The Daspk button should be
checked.

Creating a dummy state that does nothing is an interesting ploy, but I
think it may potentially lead to instability or inaccuracy because all it
does is trick the compiler into not complaining. You really want to use
IDA (called "DASPK" in NEURON) for this, not CVODES.

[ted]

***What is the format code for tabs?

In NEURON or phpBB?

In NEURON it's \t, the same as in C.
In the Programmer's Reference look in the alphabetical listing for printf-IO
http://www.neuron.yale.edu/neuron/stati ... tml#printf

phpBB has limited capabilities for preserving text formatting.
You can use HTML's pre & /pre (inside angle brackets, of course)
but the results are ugly. You _can_ select text and then click on
the Code button. phpBB can be enhanced by allowing additional
text formatting capabilities but we haven't begun to really tinker
with it yet. Check out the FAQ link near the top of the main Forum
page.

[Brad]

In NEURON or phpBB?

The formatting code of the forum. My code snippets looked ugly. Guess I'll just have to enter a few spaces by hand. Posting to forums is new to me and I didn't see the "CODE" button. What I was asking was how to do something like this...


DERIVATIVE states {
\tm' = ...
\tn' = ...
}


...and preserve the indentations that I had in the original text file.

Thanks for the help on the bigger issue, though.

[ramcdougal]

I was looking at the same file earlier today. It seems like I can do without using any solve statements or state variables, as in:

Code: Select all

NEURON {
    SUFFIX nap
    USEION na READ ena WRITE ina
    RANGE  gnabar,vhalf, K
}

UNITS {
    (mA) = (milliamp)
    (mV) = (millivolt)
}

PARAMETER {                         : parameters that can be entered when function is called in cell-setup 
    v (mV)
    ena = 50           (mV)     : Na reversal potential  (reset in cell-setup.hoc)
    K = 4.5            (1)      : slope of steady state variable
    gnabar = 0        (mho/cm2) : initialized conductance
    vhalf  = -50.4    (mV)      : half potential
}	

ASSIGNED {
    ina (mA/cm2)
    n   (1)
}

INITIAL {
    : nothing to do
}

BREAKPOINT {
    n = 1 / (1 + (exp((vhalf - v)/(1 (mV)))/K))  : steady state value
    ina = gnabar*n*n*n*(v-ena)
}

Is this legit? It compiles, passes modlunit, and cvode doesn't complain*.

Should the K be inside the exponential and have units of mV? That would eliminate the weird divided by 1 mV and it seems right when I compare this file to, say, car.mod from the same source.

* I haven't gotten the model containing the file to advance past t=0, but this file alone is not the problem. I just meant that cvode doesn't complain about me including this mechanism.

[ted]

Time to raise the issue of "correctness of the implementation," which was prematurely nixed when this thread was first opened.

The code may "work" in the sense that it may contain no syntax errors and may generate no run time errors, and that would be fine, but for one problem: there appears to be a mismatch between the implementations and the intent of their authors. In the code as written, K is in the wrong place to have any effect on dn/dv when n == 0.5 -- it merely shifts the value of v at which n == 0.5. It needs to be in the denominator of the argument to exp(). If K is placed in that location, and its declaration in the PARAMETER block asserted units of (mV), units consistency would require deleting the (1 (mV)).

[ramcdougal]

Thanks for clarifying about the K.

Can I safely include mod files without SOLVE statements without affecting the stability of the rest of the system?

[ted]

Thanks for clarifying about the K.

I should have remarked on it four years ago, but my brain must have turned off the instant I read "Putting aside all considerations about the accuracy/correctness . . . "

I wonder whatever happened to Brad and his project. I hope he wasn't going to use those simulation results for anything that MicroSoft's EULAs always warn against, like live control of jet planes or nuclear reactors, life-and-death situations in the operating room, computerized hedge fund management etc..

Can I safely . . .

My intuition is that problems are most likely to occur with mechanisms that introduce a phenomenological negative resistance that is large and activates instantaneously--which is exactly what this mechanism will do if its maximum conductance is large compared to other conductances in a model. But that remains to be seen; try it and find out.

[Darshan]

Hi,

I was trying to make this mod file work with variable time step.

The mod file uses dt to calculate the mexp and hexp

Code: Select all

PROCEDURE trates(v) {  
                      
        LOCAL tinc
        TABLE minf, mexp, hinf, hexp
	DEPEND dt	
	FROM vmin TO vmax WITH 199

	rates(v): not consistently executed from here if usetable == 1

        tinc = -dt 

        mexp = 1 - exp(tinc/mtau)
        hexp = 1 - exp(tinc/htau)
}

which is then used to calculate m and h:

Code: Select all

        m = m + mexp*(minf-m)
        h = h + hexp*(hinf-h)

Also, there is no METHOD used with SOLVE statement

Code: Select all

SOLVE states

Could you suggest ways to make this model work with variable time step.

Thanks,
Darshan

[ted]

Could you suggest ways to make this model work with variable time step.

In the Programmer's Reference read the general comments under the heading CVode at https://www.neuron.yale.edu/neuron/stat ... html#cvode
In particular note the fourth paragraph which starts with

The two major energy barriers to using the method are the requirement that hh type models be expressed in a DERIVATIVE block (instead of the explicit exponential integration step commonly implemented in a PROCEDURE)

and ends with

These issues are discussed in Channels and Events.

Click on Channels and there's the answer.

[Darshan]

Thank you for the quick response!