## Stability of steady state approximation

NMODL and the Channel Builder.
Raj
Posts: 220
Joined: Thu Jun 09, 2005 1:09 pm
Location: Hanze University of Applied Sciences
Contact:

### Stability of steady state approximation

Dear Forum,

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

### Re: Stability of steady state approximation

The problem almost certainly arises in the expression for minf, which contains
subexpressions that become indeterminate when v == 35 (as does the expression
for malpha). Try rewriting these using L'Hospital's rule and see if that doesn't fix
the problem. Since lim x / (exp(k*x) - 1) as x -> 0 is 1/(k*exp(k*x)), instead of

Code: Select all

``malpha = (-v-35)/(10*(exp((-v-35)/10)-1) )``
write this

Code: Select all

``````: malpha = (-v-35)/(10*(exp((-v-35)/10)-1) )
malpha = vtrap((-v-35), 1/10) / 10
. . .

FUNCTION vtrap(x, k) {
if (fabs(k*x) < 1e-6) {
vtrap = 1/ (k * exp(k*x))
} else {
vtrap = x / (exp(k*x) - 1)
}
}``````
While you're at it, why not make your code more readable and easier to maintain
by calculating malpha and mbeta first, so you can then just write

Code: Select all

``minf = malpha/(malpha + mbeta)``
As always, test everything--I may have made a typographical error in the
above code, and even if I didn't, you might when you try to use it. Of course,
if we're both careless but real lucky, your error might compensate for mine . . .
Raj
Posts: 220
Joined: Thu Jun 09, 2005 1:09 pm
Location: Hanze University of Applied Sciences
Contact:

### Re: Stability of steady state approximation

Unfortunately the corrected code, although better and less prone to future bugs, does not solve the stability problem. I made a test set http://www.vanelburg.net/neuronforum/Te ... bility.zip which comes without the complications of my complete simulation.

The dynamical mechanism seems stable. The stability of the steady state mechanism seems to be highly dependent on the decay time of the synapse, both with the expsyn synapse and the exp2syn synapse. For decay times of 5 ms and higher both mechanisms seem to work fine but for shorter decaytimes like 2 ms the steady state mechanism becomes instable.

Reading in the NEURON documentation on the numerical methods I read about stiff dynamics, i.e. dynamics with widely varying time constants. Now I'm wondering whether the steady state approximation is introducing excessive stiffness (in the form of a time constant equal to 0) into the dynamics, leading to instability of the steady state approximation?
ted