Now some comments about stability.
The authors concocted the mathematical form of the instantanously activating regenerative current (second term on the right hand side of their Eq. 3) simply to produce threshold and spiking behavior. Instantaneously activating regenerative currents can pose special challenges to the stability of numerical integration. Attempts to deal with this may result in a NEURON implementation that is stable at the cost of generating simulation results that differ quantitatively from the authors' original implementation. If you're lucky, the quantitative differences will not be so large as to cause qualitative differences which might destroy the utility of the NEURON implementation.
One modification that might fix the stability problem is to introduce a simple first order delay; this seems most likely to avoid destructive qualitative differences. Another is to substitute a different mathematical form for the voltage-dependence of the regenerative current; to me, this seems more problematical.
Let me know whether the following changes, which introduce a first order delay between "membrane potential" and the regenerative current, improve stability sufficiently. If problems persist, then I will propose a specific alternative form for the voltage dependence of m.
PARAMETER block: declare
taum = 0.1 (millisecond)
ASSIGNED block: declare
minf (1)
STATE block: declare
m (1)
BREAKPOINT block: change
iion = gs*vs - ge*deltaT*exp((vs-vt)/deltaT)
to
iion = gs*vs - ge*deltaT*m
INITIAL block: insert
rates()
m = minf
iion = gs*vs - ge*deltaT*m
DERIVATIVE block: insert
rates()
m' = (minf - m)/taum
Also add the following PROCEDURE block:
Code: Select all
PROCEDURE rates() {
minf = exp((vs-vt)/deltaT)
}
