Variable dt methods alter output; at_time() missing?

[tcleland]

Hi all --

I have an issue that may or may not be a problem, but is certainly interesting, and one peculiar associated problem.

I have a model cell exhibiting subthreshold oscillations and bursts of spikes, in which burst length depends on the level of (steady) current injection. Started from default parameters, simulations become stationary after 1000-2000 ms simulation time, which I now avoid by initializing all the free variables based on a set of values taken from a timepoint during stationarity. Even if I change parameters slightly, stationarity comes much more quickly; if I change them a lot, then I make a new initial-conditions variable set. Works fine.

I am trying to get variable dt to work, to take advantage of its speed bonus. The NMODL files are all OK, but the code as a whole has been giving me trouble tied to the insertion of an IClamp. I have solved this, I think (no error upon initialization, runs fine in multiple variable-dt modes as well as in fixed-timestep mode). However, the results I get are different depending upon whether I use fixed-timestep, CVODE, or DASPK, and also on the error level I use (default .001 or more stringent 0.0001). Specifically, I'm just counting the pattern of spikes-per-burst... fixed-timestep at a certain input level is stationary at 2 spikes/3 spikeless oscillations/repeat. Solving this with the different variable-dt methods gives me quasi-chaotic burst patterns (3, four x 2, 4, four x 2, 3, three x 2, 3, five x 2, 4, four x 2, 3, four x 2 is one example), and that's even neglecting to count the number of inter-burst oscillations. Reducing the tolerable error tenfold changes the burst pattern, and the pattern still also differs at this error level between CVODE and DASPK.

So, I would be very excited about this, except that I assume it's a numerical error of some sort based on my inappropriate expectations of the variable-dt solvers. But I can't really fathom what this would be. There is constant current injection (IClamp), but it doesn't vary over the course of the simulation (except at time=0 when it steps on). Plus, I'm simulating out to 10000 ms and the differences in spike patterns depending on solver mechanism persist (There are also clear differences based on initial conditions than dissipate after a second or three of sim time -- I'm not worried about those).

So that's the big maybe-a-problem.

One of the things I wanted to try was to use at_time() to tell NEURON about the 'step' at time=0, just in case. Probably wouldn't be the solution. But NEURON doesn't recognize the existence of at_time() in any form.

oc>at_time(0.1)
nrniv: at_time undefined function
near line 240
at_time(0.1)
^
at_time(0.1 )
oc>

Has at_time() been deprecated? (It's still in the documentation). I've loaded nrngui.hoc -- is there some other /lib/hoc library not in the nrngui.hoc cascade that I'm somehow not loading? Hopefully this is a simple problem.

I'm using NEURON -- Release 5.9.10 (1601) 2006-11-01. But at_time() is also missing in my NEURON 5.8 installation.

I'd appreciate any thoughts or solutions at all... why might I be getting these persistent differences among the different solvers? Why isn't at_time() working?

Thanks a lot, folks --

Thom

[ted]

Interesting questions, Thom.

Remember the slide from the NEURON course that shows an HH spike computed
under different conditions? If the stimulus current pulse is big (well above threshold),
the simulations were nearly identical, regardless of whether computed with the
default implicit Euler method (first order accuracy), Crank-Nicholson (second
order accuracy) with the same dt (0.025 ms), Crank-Nicholson with dt = 0.001 ms,
or adaptive integration (cvode, atol = 0.01). However, if the stimulus current
pulse is just a bit above threshold, the default implicit integrator's simulation
generated a spike that was smaller and about 0.3 ms later than the spikes computed
with the other integrators. The point of that slide was to demonstrate the sensitivity
of simulation results to simulation parameters (integration method, time step) when
solution trajectories pass close to spike threshold.

I have a model cell exhibiting subthreshold oscillations and bursts of spikes

Subthreshold oscillations can be particularly susceptible to numerical errors.
During the oscillations, the trajectory repeatedly passes close to spike threshold
(near each successive peak). At each of these close approaches to threshold,
solutions accumulate more and more error. Furthermore, the error that accrues in
one cycle of oscillation becomes part of the "initial condition"from which the next
cycle of oscillation is computed, so error doesn't just accumulate, it is amplified,
growing exponentially with time (i.e. according to the Lyapunov exponent of the
system of equations that is being simulated).

"So why does this affect burst length, number of spikes etc.?"

Here's my educated guess.
Burst duration and interburst interval (during which the slow oscillation occurs)
are governed by the balance between currents that favor regenerative depolarization/
hyperpolarization, and currents that oppose it. Burst duration is controlled by the
cumulative inactivation of regenerative inward currents (usually sodium and/or calcium
currents) plus cumulative activation of nonregenerative outward currents (various
voltage- and/or calcium-gated potassium currents). In the interburst interval, where the
subthreshold oscillation occurs, the inward current mechanisms slowy recover from
inactivation while the outward currents gradually deactivate (and cai falls because of e.g.
diffusion, buffering, accumulation in organelles, extrusion by active transport or
exchange). Anything that affects the interburst interval will affect the degree of recovery
of inward currents--and that in turn will affect the starting point from which the next burst
is launched, therefore changing burst duration and number of spikes etc..

the results I get are different depending upon whether I use fixed-timestep, CVODE, or DASPK, and also on the error level I use (default .001 or more stringent 0.0001).

Right. I'd bet a nickel that your model includes calcium accumulation and a pump. Even if
it doesn't, I'd bet that the same absolute error tolerance is not appropriate for all
mechanisms. Some mechanisms may need tighter tolerances.

Try the VariableTimeStep tool's "Atol Scale Tool". Make sure that the "Use variable dt"
checkbox shows a red check mark, then click on Atol Scale Tool. In the panel that pops
up, click on the "Analysis Run" button. Then click on the Rescale button. To discover what
this tool is doing, click on its Help button and read the description. You can close the
Atol Scale Tool and save the VariableTimeStep tool to a session file for future reuse; it
will "remember" the nondefault scale factors.

I would be very excited about this, except that I assume it's a numerical error of some sort

True. Any oddities that occur in simulations that involve close approaches to firing
threshold should be regarded with a skeptical eye.

at_time() has nothing to do with any of this. at_time() is a feature of NMODL, isn't callable
from hoc, and is used in the IClamp's own NMODL source code (see stim.mod in
c:\nrnxx\src\nrnoc\stim.mod (MSWin) or nrn-x.x/src/nrnoc/stim.mod (Linux)). If at_time()
wasn't working, neither would IClamp.

[tcleland]

Thanks, Ted. I couldn't figure out why the differences between methods persisted when I narrowed atol (even to the point where variable dt was running slower than fixed-step). I guess that my (uniform) atol settings may have been overkill for most variables and insufficiently precise for the calcium mechanisms (good call there). I'll do as you suggest.

thanks again --
Thom