Page 1 of 1

Rounding error for function in derivative block

Posted: Fri Mar 13, 2015 7:56 pm
by Tom Close
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)
    (nF) = (nanofarad)
}

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)
    (nF) = (nanofarad)
}

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

Re: Rounding error for function in derivative block

Posted: Mon Mar 16, 2015 1:43 pm
by hines

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.

Re: Rounding error for function in derivative block

Posted: Mon Mar 16, 2015 6:57 pm
by ted
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.

Re: Rounding error for function in derivative block

Posted: Mon Mar 16, 2015 9:10 pm
by Tom Close
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

Re: Rounding error for function in derivative block

Posted: Mon Mar 16, 2015 9:38 pm
by ted
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).

Re: Rounding error for function in derivative block

Posted: Thu Aug 31, 2017 5:54 am
by borismarin
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)
}

Re: Rounding error for function in derivative block

Posted: Sun Sep 03, 2017 12:00 am
by ted
The DE for z involves both z and cai, so derivimplicit is the correct integration method.