IClamp affects potential even though amp=0

[bosse89]

I am trying to do a phase-response curve (PRC) for my neuron. I did this before in a single-compartment neuron without any problems. This new neuron is multicompartment and a lot of different currents.
In a prc, first you have a spike-train, then you stimulate with a low current for a small amount of time at different points of one of the spikes phase.
When I set IClamp-stimulus current amplitud equal to 0 it seems I dont get the same voltage when I change the duration of the pulse.
The duration of the pulse should not matter since the amplitude is 0 anyway?
Why does these kind of things affect the spike train? The difference is small but noticable and changing my PRC calculation.

[ted]

First a couple of questions:
What version of NEURON are you using? (start NEURON and tell me what is in the first line that it prints to its terminal)
Are you using fixed time step integration, or variable time steps?

[bosse89]

NEURON -- VERSION 7.3 ansi (1078:2b0c984183df) 2014-04-04

Not sure I understand the question: "Are you using fixed time step integration, or variable time steps?"
I run time steps in neuron in this way:
"while (t < tstop) { fadvance() }"

I re-run my neuron-simulation using matlab and changing input parameters , telling neuron when to stimulate with IClamp.

If I did not answer your question, could you rephrase it?

[ted]

I'm going to have to reproduce the problem in order to be able to tell you how to fix it. Can you zip up just enough code that will allow me to see this for myself, and send it to
ted dot carnevale at yale dot edu
Please see
What to include in a zip file, and what to leave out
viewtopic.php?f=28&t=560

I re-run my neuron-simulation using matlab and changing input parameters , telling neuron when to stimulate with IClamp.

Please be sure to include the matlab code that you use to control NEURON.

[ted]

Your prc.hoc file uses variable time step integration (or, more correctly, adaptive integration)*, and that's what causes the effect you encountered**.

*--How do I know? It calls cvode.active with an argument of 1.

**--Why does this happen? If adaptive integration is used, a parameter discontinuity forces re-initialization of the integrator. An IClamp object forces the integrator to re-initialize when t==del and again when t==del+dur. This happens even if the IClamp's amp parameter is 0. This changes the times at which the solution is computed, and the trajectory of the solution; after all, numerical integration is only approximate. The resulting perturbation is typically small and doesn't have a significant effect on qualitative results, but you'll have to decide for yourself about your own particular application.

"Why can't NEURON test to see if the IClamp's amp parameter is 0, and if it's 0, just don't bother to re-initialize the integrator?"

That trades one problem for another. Why? Because if
|IClamp.amp| < float_epsilon (read about float_epsilon in the Programmer's Reference)
the simulation will proceed unperturbed, but if
|IClamp.amp|>=float_epsilon
then all of a sudden the simulation will be perturbed at del and del+dur. And don't forget that float_epsilon can be changed by the user. Sounds likely to cause difficulties.

"Well, what else can I do?"

You can reduce the perturbation by using a smaller absolute error tolerance with cvode.atol() (the default value is 1e-3), or by using cvode.atolscale() to apply well-chosen multipliers to the various states (see Programmer's Reference for documentation of these). Or forget about adaptive integration and stick with fixed time steps.

"But if I switch to fixed time step integration, my simulations take forever!"

Then design a simulation protocol that doesn't waste time. If, for example, you have a model that takes a long time to settle into a stable trajectory, run a single "warmup" simulation that allows the model to converge on that trajectory. After it settles down, use the SaveState class to capture the model's states at a point on that trajectory,
save them to a file, and then for future runs you can just read the states from that file and use them to put the model right back on that trajectory, at the exact point where you stopped it?

[bosse89]

I will try the fixed time steps. Thank you!

[ted]

The resulting perturbation is typically small and doesn't have a significant effect on qualitative results, but you'll have to decide for yourself about your own particular application.

In making this decision, first realize that all experimental and computational approaches for estimating phase response curves are approximations. Do you know the accuracy of the particular method that you are using? Since you're determining the prc of a computational model, do you know the accuracy of your model cell's anatomical and biophysical parameters? I won't go through the entire list of potential error sources, but before you abandon adaptive integration, you might want to ask first if the error that it introduces is small compared to other error sources.

You might also try this computational experiment:
Generate a series of prcs using fixed time steps with two or three different values of dt. Presumably the smallest dt will give the most accurate prc, but that's not the point here. The aim is to see how the prc estimates converge as dt becomes smaller. Then generate another prc with adaptive integration, and compare it with the prcs you generated with fixed dt. If you find that it is sufficiently close to the limit of the fixed dt prc as dt->0, then go ahead and use adaptive integration.