Results Don't Converge with Decreasing Time Step

[ajpc6d]

Why would NEURON give different results when I vary the time step outside the stimulus period? Below, I contextualize this question in the abstract, but I can provide a code example if necessary.

My situation:
I have an extremely short-duration E-field pulse I'm applying to a neuron model to assess excitation threshold, as for a strength-duration curve. To make each simulation manageable (in terms of computation time), I've used CVode().event() objects to decrease the time step to an appropriate size during the pulse duration and return to baseline afterwards. Simultaneously, I build the array representing my stimulus in the time domain such that every step NEURON takes should find an associated value (and interpolate further, just in case). This all seems to work well, and the model fires an action potential at a certain threshold.

However, when I tested decreasing the 'baseline', or extra-stimulus, time step size, the excitation threshold increased significantly. And as far as I can tell, there is no limit to this inverse relationship (i.e., as small as I have the patience to reduce the time step, the threshold rises proportionately).

My expectation is that there ought to be a value for time step size below which further variation does not affect the results. But is this a valid expectation? If it is, what could be going wrong?

[ted]

Alas, Carnac the Magnificent has shuffled off his mortal coil, leaving behind those who can only scry bugs in software that they hold in hand, or at least run on their own machines to reproduce whatever it is. So if you zip up just enough code to reproduce what you describe, I'll see what I can discover and get back to you with informed advice.

That said, this
"I have an extremely short-duration E-field pulse I'm applying to a neuron model to assess excitation threshold" and "
and this
"I've used CVode().event() objects to decrease the time step to an appropriate size during the pulse duration and return to baseline afterwards"
bring a couple of thoughts to mind. There are three reasons why it is useful to employ a brief stimulus when trying to determine excitation threshold. First, keeping stimulus duration much shorter than membrane time constant means that charge is deposited on membrane capacitance so quickly that it doesn't have time to escape from the cell through ion channels. Second, for a cell with extended geometry, keeping stimulus duration short means that the stimulus ends before the spatial gradient of transmembrane potential that it elicited can be dissipated by longitudinal current flow away from where charge entered. And these two consequences have a third consequence that favors excitation: local membrane depolarization is established so quickly that spike sodium channels are driven to open before they deactivate, and delayed rectifier (spike potassium) channels lag behind.

How "fast" should the stimulus be? For most cells and most voltage gated spike channels, 0.1 ms is probably fast enough. How can you tell if your stimulus is short enough? The same way you discover if dt, your fixed time step, is short enough: run a computational experiment. And don't bother changing dt in the middle of a simulation. All you're doing is complicating your code. If you're really concerned about numerical accuracy, use secondorder and a dt of 0.001 ms, or use adaptive integration (but make sure that the various states' error tolerances have been properly scaled; let me know if you aren't sure how to do that).

"To make each simulation manageable (in terms of computation time)" keep run times short. Do this by
(1) making sure the model is properly initialized so it starts out in steady state and doesn't need to run for hundreds or even tens of ms before it settles down
and
(2) applying the stimulus as soon as reasonable (e.g. after 1 ms that demonstrates the model was in steady state)
and
(3) letting the simulation run only long enough to demonstrate that a spike was elicited. In the case of an axon model, you might want to know that the spike actually propagates, but at most that only requires waiting for a few ms of model time after the stimulus.
For reasons that are totally unclear, it is not uncommon to see code that produces simulations in which a model sits silent, doing nothing useful, for hundreds or thousands of ms, after which a single stimulus is applied that elicits a single spike that lasts about 1 ms, yet the simulation continues for many more hundreds or thousands of ms.

[ajpc6d]

zip up just enough code to reproduce what you describe

Will do, I thought it may have more to do with NEURON than any specific code chunk.

If you're really concerned about numerical accuracy, use secondorder and a dt of 0.001 ms, or use adaptive integration (but make sure that the various states' error tolerances have been properly scaled; let me know if you aren't sure how to do that).

Indeed, I'm a bit lost by this. Appreciate anything you can offer.

(1) making sure the model is properly initialized so it starts out in steady state and doesn't need to run for hundreds or even tens of ms before it settles down

I thought there must be a way to do this, but I could never get it to work. I tried to run a separate 'equilibration' simulation to get the necessary state variable values to initialize the model to a perfect resting potential. But when I ran the 'real' simulation with a time-matched stimulus vector played into distributed mechanisms, things went awry (I believe because I could not get h.t to set where I wanted it).
But maybe I made it more complicated than it really needs to be. Is there a 'standard' process to properly initialize a model?

[ted]

I thought it may have more to do with NEURON than any specific code chunk

NEURON is a "permissive environment" in the sense that there are lots of ways to write code that runs without error messages but actually does something other than what its author thought it was doing. In that sense, yes, it may have a lot to do with NEURON, even if the problem actually lies in user-written code. But sometimes user-written code exposes an obscure bug. In either case, interactive debugging is often needed to diagnose and fix the problem.