I'm using the mod-file included below in a neuron model, depending on whether I use the steady state activation approximation in the breakpoint or the dynamic activation I get completely different stability properties. If I use the steady state I have to go to dt=0.001 ms for stable numerics but the dynamic approximation works fine with dt=0.025. Is there a explanation for this?
More importantly can it be and is it worthwhile mending so one can use the steady state approximations with large timesteps?
Regards,
Ronald van Elburg
(I'm aware that the present mod-code looks odd because it contains calculations both for the dynamic and steady state approximation, but ... I am experimenting!)
Code: Select all
TITLE Sodium Current
COMMENT
Author: Ronald van Elburg
Taken from model for fast spiking neuron in TCW_2002
Modifications:
ID Date Authors Email Description
M_001
ENDCOMMENT
INDEPENDENT { t FROM 0 TO 1 WITH 1 (ms) }
UNITS {
: (Abbreviation)= (Unit)
(mV) = (millivolt)
(mA) = (milliamp)
(pS) = (picosiemens)
(um) = (micron)
: Abbreviation = (Constant) (Unit)
}
NEURON {
SUFFIX NaPyr
USEION na READ ena WRITE ina
RANGE gbar, ina, minf,malpha, mbeta
GLOBAL v_table_min, v_table_max
}
PARAMETER {
: Parameter =Initial Value (Units) Description
gbar = 350 (pS/um2)
v (mV)
ena (mV)
phih = 5 (1)
phim =5 (1)
: Table settings
: Do not change these without looking at the table, make sure that he value
: of m for v=-35 is obtained through interpolation. If you fail to do this
: it is likely you will obtain division by zero errors.
:
v_table_min = -120 (mV)
v_table_max = 100 (mV)
}
ASSIGNED {
: Parameter Units Description
ina (mA/cm2)
minf (1)
halpha (1/ms)
hbeta (1/ms)
malpha (1/ms) :Dynamic activation
mbeta (1/ms) :Dynamic (de)activation
}
STATE {
h (1)
m (1)
}
BREAKPOINT {
settables(v)
SOLVE states METHOD cnexp
:ina =(1e-4)*gbar * minf * minf * minf * h * (v - ena) :Steady state activation approximation
ina =(1e-4)*gbar * m * m* m * h * (v - ena) :Dynamic activation
}
INITIAL {
settables(v)
h=halpha/(halpha+hbeta)
m=minf
}
DERIVATIVE states {
h' =phih* ( halpha*(1-h) -hbeta*h)
m' =phim* ( malpha*(1-m) - mbeta*m) :Dynamic activation
}
UNITSOFF
PROCEDURE settables(v (mV)) {
TABLE minf, halpha, hbeta, malpha, mbeta FROM v_table_min TO v_table_max WITH 961 :For dynamic activation add malpha, mbeta
halpha = 0.07 * ( exp( ( - v - 58 ) / 20) )
hbeta = 1/( exp((-v -28 )/10 ) +1)
minf = (-v-35)/(10*(exp((-v-35)/10)-1) )/((-v-35)/(10*(exp((-v-35)/10)-1) )+4*exp((-v-60)/18)) :Steady state activation approximation
malpha = (-v-35)/(10*(exp((-v-35)/10)-1) ) :Dynamic activation
mbeta = 4*exp((-v-60)/18) :Dynamic activation
}
UNITSON