classes ModelDescriptionIssues event record states active jacobian record_remove statistics atol maxorder rtol stiff atolscale netconlist solve dstates re_init statename
objref cvodecvode = new CVode()cvode = new CVode(0 - 1)This class interfaces a variant of the CVODE Software Package by Scott D. Cohen and Alan C. Hindmarsh obtained from http://www.netlib.org. The changes to this package (located in $NEURONHOME/src/cvode/cvodemlh.[ch]) have to do with handling abrupt changes in parameters (eg current stimulus pulses, synaptic conductance onset, etc).
When this class is active the finitialize/fadvance calls use the CVode integrator. The argument to the CVode constructor defines whether the instance starts out active (1) or inactive (0). With no argument the instance is inactive. Only one instance of cvode can exist at one time at present. In the default variable step context, the integrator chooses the time step and fadvance returns after one step. Local accuracy is determined by atol and rtol .
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 that my variant of the CVode solver be explicitly notified of the exact time of any discontinuity such as step changes, pulses, and synaptic conductance onset's. These issues are discussed in Channels and Events .
After your mod files are generalized it will probably be convenient to compare the default method with CVode by toggling the UseCVode checkbox in the RunControl .
Only the SEClamp works with CVode. VClamp cannot be used with this method.
Also .mod authors must take measures to handle step changes in parameters ( Events )
CVode
cvode.solve()cvode.solve(tout)For backward compatibility with finitialize/fadvance it is better to use the active method instead of calling solve directly.
CVode
cvode.statistics()
CVode
x = cvode.rtol()x = cvode.rtol(relative)The solver attempts to use a step size so that the local error for each state is less than
rtol*|state| + atol*atolscale_for_stateThe error test passes if the error in each state, e[i], is such that e[i]/state[i] < rtol OR e[i] < atol*atolscale_for_state (the default atolscale_for_state is 1, see atolscale )
CVode
x = cvode.atol()x = cvode.atol(absolute)The solver attempts to use a step size so that the local error for each state is less than
rtol*|state| + atol*atolscale_for_stateThe error test passes if the error in each state, e[i], is such that e[i]/state[i] < rtol OR e[i] < atol*atolscale_for_state
Therefore states should be scaled (or the absolute tolerance reduced) so that when the value is close to 0, the error is not too large.
(See atolscale for how to set distinct absolute multiplier tolerances for different states.)
Either rtol or atol may be set to 0 but not both. (pure absolute tolerance or pure relative tolerance respectively).
CVode
tol = cvode.atolscale(&var, toleranceMultiplier)tol = cvode.atolscale(&var)cvode.atolscale(&soma.v(.5), 1e-8) sets
the absolute tolerance multiplier for all membrane potentials everywhere.
(The syntax for merely specifying a name is admittedly cumbersome but
the function is not often needed and it avoids the necessity of
explicitly having to parse strings such as "TrigKSyn.G".)
The currently specified multiplier for that state name
is returned by the function call.
Specification of a particular STATEs absolute tolerance multiplier is only needed if its scale is extremely small or large and is best indicated within the model description file itself using the STATE declaration syntax:n
See nrn/demo/release/cabpump.mod for an example of a model which needs a specific scaling of absolute tolerances (ie, calcium concentration and pump density).state (units) <tolerance>
CVode
cvode.re_init()
CVode
x = cvode.stiff()x = cvode.stiff(0-2)1 only membrane potential computed implicitly. 0 Adams-Bashforth integration.
CVode
x = cvode.active()x = cvode.active(0)x = cvode.active(1)following two not yet implementedx = cvode.active(1, dt)x = cvode.active(tvec)This function allows one to toggle between the normal integration method and the CVode method with no changes to existing intepreter code. The return value is whether CVode is active.
With only a single 1 arg, the fadvance calls CVode to do a single variable time step.
With the dt arg, fadvance returns at t+dt.
With a Vector tvec argument, CVode is made active and a sequence of calls to fadvance returns at the times given by the elements of tvec. After the last tvec element, fadvance returns after each step.
CVode
x = cvode.maxorder()x = cvode.maxorder(0 - 12)
CVode
x = cvode.jacobian()x = cvode.jacobian(0 - 2)
CVode
objref dest_vectordest_vector = new Vector()cvode.states(dest_vector)
CVode
cvode.dstates(dest_vector) Fill the destination Vector with the values of d(state)/dt.
CVode
strdef dest_stringcvode.statename(i, dest_string)print dest_string
CVode
List = cvode.netconlist(precell, postcell, target)List = cvode.netconlist(precell, postcell, target, list)
iterator ltr() {local i, cnt
for i = 0, $o2.count - 1 {
$o1 = $o2.object(i)
iterator_statement
}
}
objref xo
for ltr(xo, cvode.netconlist(precell, "", "")) {
print xo.postcell
}
CVode
cvode.record(&rangevar, yvec, tvec)cvode.record(&rangevar, yvec, tvec, 1)During a run, record the stream of values in the specified range variable into the yvec Vector along with time values into the tvec Vector. Note that each recorded range variable must have a separate tvec which will be different for different cells. On initialization the yvec and tvec Vectors are resized to 1 and the initial value of the range variable and time is stored in the Vectors.
To stop recording into a particular vector, remove all the references either to tvec or yvec or call record_remove .
If the third argument is present and equal to 1, the yvec is recorded only at the t values in tvec.
CVode
cvode.record_remove(yvec)
CVode
cvode.event(t)
CVode Channels Concentrations EventsThe following aspects of model descriptions (.mod files) are relevant to their use with CVode.
KINETIC block - No changes required.
DERIVATIVE block - No changes required. The Jacobian is approximated as a diagonal matrix. If the states are linear in state' = f(state) the diagonal elements are calculated analytically, otherwise the diagonal elements are calculated using the numerical derivative (f(s+.01) - f(s))/.001 .
LINEAR, NONLINEAR blocks - No changes required. However, at this time they can only be SOLVED from a PROCEDURE or FUNCTION, not from the BREAKPOINT block. The nrn/src/nrnoc/vclmp.mod file gives an example of correct usage in which the function icur is called from the BREAKPOINT block and in turn SOLVE's a LINEAR block. If desired, it will be a simple matter to allow these blocks to be solved from the BREAKPOINT block.
SOLVE PROCEDURE within a BREAKPOINT block - Changes probably required. Such a procedure is called once after each return from CVode.solve().
ModelDescriptionIssuesThe SOLVE PROCEDURE form was often used to implement the exponential integration method for HH like states and was very efficient in the context of the Crank-Nicholson like staggered time step approach historically used by NEURON. Furthermore the exponential integration often used tables of rates which were calculated under the assumption of a fixed time step, dt. Although it can still be used under some circumstances, the usage to integrate states should be considered obsolete and converted to a DERIVATIVE form. To do this, 1) replace the PROCEDURE block with a DERIVATIVE block, eg. DERIVATIVE states { m' = (minf - m)/mtau ... } 2) replace the SOLVE statement in the BREAKPOINT block with SOLVE states METHOD cnexp 3) if using tables, store mtau instead of (1 - exp(-dt/mtau)) The nmodl translator will emit c code for both the staggered time step and high order variable time step methods. The only downside is slightly less efficiency with the staggered time step method since the exp(-dt...) is calculated instead of looked up in tables.
In summary, no model should anymore depend on dt.
ModelDescriptionIssues
ModelDescriptionIssues
How does one handle events? This is really the only serious difficulty in writing models that work properly in the context of a variable time step method. All models which involve discontinuous functions of time, eg steps, pulses, synaptic onset, require special provision to notify the integrator that an event has occurred within this time step, ie between t-dt and t. If this is not done, the time step may be so large that it completely misses a pulse or synaptic event. And if it does see the effect of the event, there is a huge inefficiency involved in the variable step method's search for the location of the event and the concomitant tremendous reduction in size of dt.
So, if you change any variable discontinuously in the model at some time tevent, call call
at_time(tevent)
The user may check the return value of this function to decide
if something needs changing. Examples of the two styles of usage are:
1) Just notify and do the logic separately.
at_time(del)
at_time(del + dur)
if (t >= del && t <= dur) {
istim = on_value
}else{
istim = 0
}
2) Use the at_time return value to do the logic.
INITIAL {
istim = 0
}
...
if (at_time(del)) {
istim = on_value
}
if (at_time(del + dur)) {
istim = 0
}
Notice the requirement of initialization or else if the previous
run was stopped before del + dur the value of istim would be on_value
at the beginning of the next run.
What happens internally when at_time(tevent) is called?
The interesting case (t-dt < tevent <= t) --- First, at_time returns 0. Then CVode changes its step size to (tevent - (t-dt) - epsilon) and redoes the step starting at t-dt. Note that this should be safely prior to the event (so at_time still returns 0), but if not then the above process will repeat until a step size is found for which there is no event. CVode then re-initializes it's internal state and restarts from a new initial condition at tevent+epsilon. Now when at_time is called, it returns 1. Note that in its single step mode, CVode.solve() will return at t = tevent-epsilon, the subsequent call will start the initial condition at t = tevent + epsilon and return after a normal step (usually quite small).
The case (tevent <= t - dt) --- at_time returns 0.
The case (tevent > t) --- at_time returns 0.
Note that an action potential model with axonal delay delivering a "message" to a synaptic model may or may not think it worthwhile to call at_time at the time of threshold (I would just do my own interpolation to set t_threshold) but will certainly call at_time(t_threshold + delay) (and possibly not allow t_threshold to change again until it returns a 1);
I am sorry that the variable time step method requires that the model author take careful account of events but I see no way to have them automatically taken care of.