## Rounding error for function in derivative block

NMODL and the Channel Builder.
Tom Close
Posts: 18
Joined: Wed Sep 05, 2012 11:59 pm

### Rounding error for function in derivative block

I appear to be getting a rounding error when encapsulating a time derivative in the derivative block within a function, which I don't understand.

Using the PyNN izhikevich mod file as an example, could you explain to me why

Code: Select all

``````: Izhikevich artificial neuron model from
: EM Izhikevich "Simple Model of Spiking Neurons"
: IEEE Transactions On Neural Networks, Vol. 14, No. 6, November 2003 pp 1569-1572
:
: This NMODL file may not work properly if used outside of PyNN.
: Ted Carnevale has written a more complete, general-purpose implementation - see http://www.neuron.yale.edu/ftp/ted/devel/izhdistrib.zip

NEURON {
POINT_PROCESS Izhikevich1
RANGE a, b, c, d, u, Cm, uinit, vthresh
NONSPECIFIC_CURRENT i
}

UNITS {
(mV) = (millivolt)
(nA) = (nanoamp)
}

INITIAL {
u = uinit
net_send(0, 1)
}

PARAMETER {
a       = 0.02 (/ms)
b       = 0.2  (/ms)
c       = -65  (mV)   : reset potential after a spike
d       = 2    (mV/ms)
vthresh = 30   (mV)   : spike threshold
Cm      = 0.001  (nF)
uinit   = -14  (mV/ms)
}

ASSIGNED {
v (mV)
i (nA)
}

STATE {
u (mV/ms)
}

BREAKPOINT {
SOLVE states METHOD cnexp  : derivimplicit
i = -Cm * (0.04*v*v + 5*v + 140 - u)
:printf("t=%f, v=%f u=%f, i=%f, dv=%f, du=%f\n", t, v, u, i, 0.04*v*v + 5*v + 140 - u, a*(b*v-u))
}

DERIVATIVE states {
u' = deriv_u(u)
}

FUNCTION deriv_u(u) {
deriv_u = a*(b*v - u)
}

NET_RECEIVE (weight (mV)) {
if (flag == 1) {
WATCH (v > vthresh) 2
} else if (flag == 2) {
net_event(t)
v = c
u = u + d
} else { : synaptic activation
v = v + weight
}
}
``````
produces different results (with 'u' diverging) to

Code: Select all

``````: Izhikevich artificial neuron model from
: EM Izhikevich "Simple Model of Spiking Neurons"
: IEEE Transactions On Neural Networks, Vol. 14, No. 6, November 2003 pp 1569-1572
:
: This NMODL file may not work properly if used outside of PyNN.
: Ted Carnevale has written a more complete, general-purpose implementation - see http://www.neuron.yale.edu/ftp/ted/devel/izhdistrib.zip

NEURON {
POINT_PROCESS Izhikevich2
RANGE a, b, c, d, u, Cm, uinit, vthresh
NONSPECIFIC_CURRENT i
}

UNITS {
(mV) = (millivolt)
(nA) = (nanoamp)
}

INITIAL {
u = uinit
net_send(0, 1)
}

PARAMETER {
a       = 0.02 (/ms)
b       = 0.2  (/ms)
c       = -65  (mV)   : reset potential after a spike
d       = 2    (mV/ms)
vthresh = 30   (mV)   : spike threshold
Cm      = 0.001  (nF)
uinit   = -14  (mV/ms)
}

ASSIGNED {
v (mV)
i (nA)
}

STATE {
u (mV/ms)
}

BREAKPOINT {
SOLVE states METHOD cnexp  : derivimplicit
i = -Cm * (0.04*v*v + 5*v + 140 - u)
:printf("t=%f, v=%f u=%f, i=%f, dv=%f, du=%f\n", t, v, u, i, 0.04*v*v + 5*v + 140 - u, a*(b*v-u))
}

DERIVATIVE states {
u' = a*(b*v - u)
}

NET_RECEIVE (weight (mV)) {
if (flag == 1) {
WATCH (v > vthresh) 2
} else if (flag == 2) {
net_event(t)
v = c
u = u + d
} else { : synaptic activation
v = v + weight
}
}``````
Is there a work around for this if you want to implement multiple regimes, i.e. solve different time derivatives depending on a flag set in the mod file.

Cheers,

Tom

hines
Posts: 1558
Joined: Wed May 18, 2005 3:32 pm

### Re: Rounding error for function in derivative block

Code: Select all

``````    SOLVE states METHOD cnexp
...
u' = deriv_u(u)
``````
The cnexp method is invalid for that derivative equation.
The cnexp method requires that each y' equation be linear and involve no other state than y.
That is, cnexp was intended for use only in HH style gating equations where the form looked like
m' = (minf(v) - m)/mtau(v)
or
m' = bm(v)*(1. - m) - am(v)*m

For general (nonlinear and coupled state ) equations, use derivimplicit. I need to change the nmodl translator to generate an error message if the equation forms do not
satisfy the cnexp requirements and will look into the possibility of doing that.
I don't think "Rounding error" is a good set of keywords for this topic. Perhaps nonsense numerical results would be better.
I should mention that even derivimplicit will not result in a simulation that uses the full jacobian in this case but only the diagonal terms dv'/dv and du'/du . That is usually sufficient for numerical stability in most practical circumstances.

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

### Re: Rounding error for function in derivative block

hines wrote:

Code: Select all

``````    SOLVE states METHOD cnexp
...
u' = deriv_u(u)
``````
The cnexp method is invalid for that derivative equation.
The cnexp method requires that each y' equation be linear and involve no other state than y.
That is, cnexp was intended for use only in HH style gating equations where the form looked like
m' = (minf(v) - m)/mtau(v)
or
m' = bm(v)*(1. - m) - am(v)*m

For general (nonlinear and coupled state ) equations, use derivimplicit.
In light of which, I have now updated
http://www.neuron.yale.edu/ftp/ted/devel/izhdistrib.zip
accordingly.

Tom Close
Posts: 18
Joined: Wed Sep 05, 2012 11:59 pm

### Re: Rounding error for function in derivative block

Hi Michael,

Thanks, I switched to derivimplicit and the results are the same now. I think it would be really helpful for nrnivmodl to spit out warnings if the integration method is not applicable/optimal for the provided mechanism.

I will create a ticket so that this also gets updated in the PyNN version, and will take a closer look into the different integration options so I am a bit more aware of their limitations/trade-offs.

Cheers,

Tom

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

### Re: Rounding error for function in derivative block

Tom Close wrote:I will create a ticket so that this also gets updated in the PyNN version, and will take a closer look into the different integration options so I am a bit more aware of their limitations/trade-offs.
There's a summary in
Integration methods for SOLVE statements
viewtopic.php?f=28&t=592
which I will update to mention using derivimplicit for general (nonlinear and coupled state ) equations, and to remind users that v is indeed a state variable (even though it never appears in a mod file's STATE block).

borismarin
Posts: 1
Joined: Tue Aug 29, 2017 11:18 am

### Re: Rounding error for function in derivative block

I am also facing numerical inconsistencies when switching from cnexp to derivimplicit in a mod file that appears more than once on modeldb.

In this case, would derivimplicit be the "correct" choice (or in which sense would cnexp be "incorrect"?), given that z dynamics depends on Ca - even if the kinetics is HH-like?

Code: Select all

``````DERIVATIVE states {
rates(cai)
z' = (zInf - z) / zTau
}

PROCEDURE rates(ca(mM)) {
if(ca < 1e-7){
ca = ca + 1e-07
}
zInf = 1/(1 + (0.00043 / ca)^4.8)
}``````

ted