Rounding error for function in derivative block

NMODL and the Channel Builder.
Post Reply
Tom Close

Rounding error for function in derivative block

Post 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
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: Rounding error for function in derivative block

Post 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.
ted
Site Admin
Posts: 6286
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Rounding error for function in derivative block

Post 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.
Tom Close

Re: Rounding error for function in derivative block

Post 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
ted
Site Admin
Posts: 6286
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Rounding error for function in derivative block

Post 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).
borismarin
Posts: 1
Joined: Tue Aug 29, 2017 11:18 am

Re: Rounding error for function in derivative block

Post 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)
}
ted
Site Admin
Posts: 6286
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Rounding error for function in derivative block

Post by ted »

The DE for z involves both z and cai, so derivimplicit is the correct integration method.
Post Reply