Stability of steady state approximation

[Raj]

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]

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]

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]

[quote=]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?[/quote]
Probably. Using your toy model "as is," reducing dt to 0.005 ms or changing to adaptive
integration makes the problem go away.