Variable dt methods alter output; at_time() missing?
Posted: Sun Feb 04, 2007 5:17 pm
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
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