Rounding error for function in derivative block

NMODL and the Channel Builder.

Rounding error for function in derivative block

Postby Tom Close » Fri Mar 13, 2015 7:56 pm

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
Tom Close
 
Posts: 18
Joined: Wed Sep 05, 2012 11:59 pm

Re: Rounding error for function in derivative block

Postby hines » Mon Mar 16, 2015 1:43 pm

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

Re: Rounding error for function in derivative block

Postby ted » Mon Mar 16, 2015 6:57 pm

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

Re: Rounding error for function in derivative block

Postby Tom Close » Mon Mar 16, 2015 9:10 pm

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
Tom Close
 
Posts: 18
Joined: Wed Sep 05, 2012 11:59 pm

Re: Rounding error for function in derivative block

Postby ted » Mon Mar 16, 2015 9:38 pm

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
https://www.neuron.yale.edu/phpBB/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).
ted
Site Admin
 
Posts: 5040
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine

Re: Rounding error for function in derivative block

Postby borismarin » Thu Aug 31, 2017 5:54 am

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

Re: Rounding error for function in derivative block

Postby ted » Sun Sep 03, 2017 12:00 am

The DE for z involves both z and cai, so derivimplicit is the correct integration method.
ted
Site Admin
 
Posts: 5040
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine


Return to Adding new mechanisms and functions to NEURON

Who is online

Users browsing this forum: No registered users and 1 guest