CVODE: record, play, events

The basics of how to develop, test, and use models.
nikkip
Posts: 28
Joined: Mon Jun 24, 2013 6:09 pm

CVODE: record, play, events

Post by nikkip »

I'm using CVODE for the first time, and I have a few implementation questions. Sorry for the long post; I'm just missing some basics despite digging through the documentation/forum/book.

1) I'm running my simulation multiple times (parameter sweep). Should I call cvode.record(...) just once (outside my simulation procedure) or for every simulation (within the simulation procedure)?

2) I want to use cvode.record(...) to record the current of my IClamp point process. Is it possible to use the cvode class's record function for a point process? If so, what should the first argument be? If not, how else could I record the current time course?

3) I would also like to use cvode.record(...) to record Vm. However, I'm not clear on the syntax for the first argument (pointer to the voltage) given that my nodes are a vector:

Code: Select all

create node[axonnodes]
Let’s say I want to record Vm for node[0]. I’ve tried two ways:

Code: Select all

 objref t_vec
t_vec = new Vector()
objref Vm_vec
Vm_vec = new Vector()
cvode.record(&node[0].v(0.5),Vm_vec,t_vec)

Code: Select all

 objref t_vec
t_vec = new Vector()
objref Vm_vec
Vm_vec = new Vector()
node[0] {
	cvode.record(&v(0.5),Vm_vec,t_vec)
}
But in both cases, after running my simulation, Vm_vec and t_vec are empty (size = 0).

4) I’m using extracellular current point source stimulation (phi = I/(4*pi*sigma*r) at each node, where I is a current pulse from an IClamp). I want to set up the extracellular potential time course for each node and use myvec.play(…) to apply them. I’ve tried looking at the folder "extracellular_stim_and_rec", but I’m having trouble grasping the pieces that I need. From calcrxc.hoc, I understand that I can create a vector (rx_xtra) that contains the transfer impedance from the extracellular point source to each node of my axon. Then, using rx_xtra, I could just loop through all the nodes and all time points with a fixed time step…

Code: Select all

for K = 0, axonnodes-1 {
	for J = 0, ntsteps – 1 {
		// compute Ve at each node at each time point
		if ((J*dt > delay) && (J*dt < delay + dur)) {
			Ve[K].x[J] = I/(4*pi*sigma*r.x[K])
		}
	}
}
Is there a more efficient way, in which I’d actually use the IClamp rather than figuring out if the stimulation is on or not manually as I did above?

Then I think I’d have to apply (i.e. play) the potentials with something like…

Code: Select all

for K = 0, axonnodes-1 {
	Ve[K].play(&node[K].e_extracellular(0.5), dt)
}
I will also use cvode.event(t) to add events signaling the step up and step down of my current pulse.

Thank you!
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: CVODE: record, play, events

Post by ted »

nikkip wrote:I'm using CVODE for the first time
First a question for you: do you need cvode.record? No need for that unless you're using the local variable time step method, and you won't be using that unless (1) your model has multiple cells, and (2) you have found that simulations run quicker or are significantly more accurate when executed with the local variable time step method.

The following answers assume that you do not need cvode.record.
I'm running my simulation multiple times (parameter sweep). Should I call cvode.record(...) just once (outside my simulation procedure) or for every simulation (within the simulation procedure)?
I don't know what you mean by "simulation procedure." In the context of computational modeling it is most useful to think of the following kinds of code:
model setup code--specifies the properties of the physical system that will be represented in the computer
instrumentation code--specifies signal sources or inputs to the physical system, and which signals or variables are to be observed or recorded
simulation control code--specifies how simulations will be initialized, and what happens in the course of simulation execution (duration of simulation, integration method, dt or error tolerance, conditions or times at which perturbations will occur (think "virtual experimental protocol"))

In general, properly written model setup and instrumentation code only needs to be executed once. This is certainly the case for statements that specify which variables are to be recorded. That said, if you are going to be writing code that contains unfamiliar features of NEURON, it is a very good idea to first try those features in the context of a toy problem (a simple task involving a very simple model for a very short run time) just to make sure that there will be a close match between what you want and what your code actually does.
2) I want to use cvode.record(...) to record the current of my IClamp point process.
Unless you're using local variable timestep integration, use the Vector class's record method. In either case, make a toy model (e.g. a single compartment with an IClamp attached to it), and see if the values recorded are what you expected based on the IClamp's parameters. The documentation of Vector.record provides several examples, the first of which is
vdest.record(&var)
so for anything you might want to record, first ask "what is the name of the variable in hoc, and what is the name of the Vector to which I want to record it?" and that should give you a good idea of what the necessary statement is. Examples:
Suppose you have an IClamp called stim and you want to record its time course to a vector called ivec. The current that stim delivers is called stim.____ so the statement is ivec.record(&_______).
Suppose you also want to record time at each step in the simulation to a vector called tvec (a good idea, because you can then use the Vector class's plot method to show the time course of stim.i in a graph. Time is called t, so the statement is tvec.record(&______).
Now make a toy model and try out these statements.
3) I would also like to use cvode.record(...) to record Vm.
The documentation of Vector.record tells what to do.
my nodes are a vector:

Code: Select all

create node[axonnodes]
Actually, that's not a vector--it's an example of the use of array notation to manage multiple things (sections in this case) that have the same base name. The complete name of any one of your nodes is node where i is either a whole number or a variable whose value lies in the range 0..axonnodes-1. So if the way to record v at the midpoint of a section called terminal to a vector called dv is
dv.record(&terminal.v(0.5))
what do you think would be the way to record terminal.v(0.5) if the section's name is node[3] instead of terminal? How would you test this idea?

With regard to the examples you provide, the problem is a very understandable misinterpretation of the documentation. Referring to the documentation for Vector.record, here are the key points:
vdest.record(&var) means that the values of var will be captured to vdest after finitialize() and each fadvance().
vdest.record(&var, Dt) means that values will be captured to vdest at the first simulation times that are >= i*Dt where i = 0, 1, 2 . . .. If Dt is a multiple of dt, this will capture every "Dt/dt"th value generated in the simulation, allowing for roundoff error here or there.
vdest.record(&var, tvec), where tvec is monotonically increasing, means that values are captured to vdest at the first simulation times that are >= the values of the elements in tvec. In other words, tvec is a user-supplied vector of times at which you want the value of var to be sampled, again allowing for roundoff error here and there.

So your code examples don't capture anything because the vector of sample times is empty. But why bother with that syntax, unless you are using adaptive integration and really want to capture values at specific times? And even then you only need to use Vector record unless the model has multiple cells and the simulation is being executed with local variable time steps.

4) I’m using extracellular current point source stimulation (phi = I/(4*pi*sigma*r) at each node, where I is a current pulse from an IClamp). I want to set up the extracellular potential time course for each node and use myvec.play(…) to apply them. I’ve tried looking at the folder "extracellular_stim_and_rec", but I’m having trouble grasping the pieces that I need. From calcrxc.hoc, I understand that I can create a vector (rx_xtra) that contains the transfer impedance from the extracellular point source to each node of my axon. Then, using rx_xtra, I could just loop through all the nodes and all time points with a fixed time step…

Code: Select all

for K = 0, axonnodes-1 {
	for J = 0, ntsteps – 1 {
		// compute Ve at each node at each time point
		if ((J*dt > delay) && (J*dt < delay + dur)) {
			Ve[K].x[J] = I/(4*pi*sigma*r.x[K])
		}
	}
}
Is there a more efficient way, in which I’d actually use the IClamp rather than figuring out if the stimulation is on or not manually as I did above?

Then I think I’d have to apply (i.e. play) the potentials with something like…

Code: Select all

for K = 0, axonnodes-1 {
	Ve[K].play(&node[K].e_extracellular(0.5), dt)
}
I will also use cvode.event(t) to add events signaling the step up and step down of my current pulse.
This is about a very different topic, so I will split it off into a different thread and answer it there at a later time.
nikkip
Posts: 28
Joined: Mon Jun 24, 2013 6:09 pm

Re: CVODE: record, play, events

Post by nikkip »

Thank you Ted. So let me start with addressing the issue of the possible need for CVODE. I have a model for which I'm extracting the conduction speed and strength-duration curve for validation against published experimental data. This model is a modification of the MRG model. Based upon this thread...

http://www.neuron.yale.edu/phpbb/viewto ... f=8&t=2271

...it seemed that the accuracy of my conduction speed estimate would be significantly compromised without CVODE (or an unfeasibly small fixed time step). Unfortunately, it is not easy for me to try a local variable time step with my current code structure. But perhaps I've been over-concerned and should simply use my current code without CVODE.
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: CVODE: record, play, events

Post by ted »

nikkip wrote:it seemed that the accuracy of my conduction speed estimate would be significantly compromised without CVODE
Did I say not to use adaptive integration (aka "CVODE" aka "CVode" aka "cvode")?

Time for some disambiguation. While the CVode class is used to manage several different adaptive integration methods, this does not mean that all adaptive integration methods are the same. In particular note that some methods use global time steps (state variables of all cells are always at the same point in time) with the local variable timestep method (each cell advances at its own pace, so the ith and jth cells will be at times ti and tj and there is no guarantee that ti = tj except at the start and end of a simulation). This point may not be sufficiently explicit among the welter of details that are covered so concisely in the Programmer's Reference, so look back at the materials from the NEURON Summer Course to find two figures that contrast simulations of a network of three cells using global and local variable time step integration, and a sequence of figures that detail the step-by-step advance in time of a local variable time step simulation of a two cell network.

I hope that clarifies the difference between global and local variable time step methods. Now you know what I meant when I wrote that you don't need local variable time steps in a model that has only one cell, and perhaps my other comments will be more understandable. This would be a good time to review the programmer's reference documentation of the CVode class and the Vector record method.
nikkip
Posts: 28
Joined: Mon Jun 24, 2013 6:09 pm

Re: CVODE: record, play, events

Post by nikkip »

Very helpful clarifications. I've sorted out most of my issues (and I'll try to make a post this weekend about my solutions for future readers). I think I have the myvec.play(...) figured out as well.

I'm stumped on one bug though, which is causing my code to crash with the following error...
t: 1.975
t: 1.9796648
t: 1.98
t: 1.985
t: 1.99
t: 1.995
t: 2
t: 2
t: 2
t: 2
err=-6
)
, , )
, , , , , )
)
IDASolve-- at t = 2 and h = 6.11934e-16, the error test
failed repeatedly or with |h| = hmin.

1 nrniv: variable step integrator error
1 in MRGaxon_unmyelinated_StrengthDuration_v5.hoc near line 588
1 {pc.runworker()}
^
1 fadvance(1 stimul_extracellular(2.01-2e+080.51 find_thresh(12.01200.501 ParallelContext[0].runworker(application called MPI_Abort(MPI_COMM_WORLD, -1) - process 1
[0]0:Return code = 0, signaled with Quit
[0]1:Return code = 255
So I'm printing out the time points within my while(t<tstop) loop, and it gets stuck at t=2ms, which is when my stimulus turns on, even though I've added a cvode.event at that time.

Here's why I'm very perplexed. I've been running 4 different types of runs: with and without CVODE, and with intracellular or extracellular stimulation. They all work, except for CVODE with extracellular stimulation, which produces the above error (i.e. integration issues at the time when my stimulus turns on). But my extracellular stimulation code works with fixed time step. And analogous code with an intracellular stimulus using CVODE (with a cvode.event at 2ms as well) works.

I realize code snippets may be required to elucidate the issue, but I'm quite unclear as to what the cause is. Further, in case the issue has an obvious solution, I thought I'd start by stating the problem and error.

Thank you!
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: CVODE: record, play, events

Post by ted »

Abrupt discontinuities in a driving function can cause problems. Example: Vector play with this pair of vectors

Code: Select all

t V
0 0
1 0
1 1
2 1
2 0
3 0
describes a square pulse that starts abruptly at 1 ms and ends abruptly at 2 ms, even if played with interpolation. For such cases it may be useful to replace the jumps by ramps, e.g.

Code: Select all

t V
0 0
0.9 0
1.1 1
1.9 1
2.1 0
3 0
nikkip
Posts: 28
Joined: Mon Jun 24, 2013 6:09 pm

Re: CVODE: record, play, events

Post by nikkip »

Quick follow-up questions:

1) How is your last example a ramp...? Doesn't it still turn on and off abruptly?

2) Did I just get lucky then that I did not get an error when I used an intracellular current pulse?

3) Should I still use cvode.event of some form when using a ramp?

4) How do I know how gentle of a slope is required to avoid the error?

5) Is there an easy way (e.g. an existing function) to implement a ramp..? All I'm trying to do is apply a pulse of a certain amplitude and duration to get conduction speed and strength-duration data.

Thanks!
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: CVODE: record, play, events

Post by ted »

nikkip wrote:1) How is your last example a ramp...? Doesn't it still turn on and off abruptly?
Excellent question. I meant to stipulate that the Vector should be played "with interpolation," which means the t, V pairs are treated as breakpoints of a piecewise linear function and that V values at times that lie between any of the t, V pairs are calculated by linear interpolation. Using fixed dt simulation with dt = 0.025 ms, the transition from 0 at 0.9 ms to 1.1 at 1 ms would take place over 8 time steps. Using adaptive integration, there would be one or two very short time steps after each breakpoint at which slope changes abruptly, and after that step length would be governed entirely by the dynamics of the model's system equations until the next slope change is encountered.
2) Did I just get lucky then that I did not get an error when I used an intracellular current pulse?
Not necessarily. An intracellular current pulse with abrupt onset does not make local v change instantaneously--instead, v changes very gradually ("the high frequency components of the injected current are filtered out by membrane capacitance"). A change of extracellular potential that turns on abruptly affects multiple extracellular nodes suddenly, and it must also affect the potential at multiple intracellular nodes just as abruptly. This may be difficult for the adaptive integration algorithms to deal with.
3) Should I still use cvode.event of some form when using a ramp?
I don't see why it is needed at all, unless you are using it to force an abrupt change in some parameter. If you do that, you must follow it with a cvode.re_init because what you have done is equivalent to starting the solution of a new initial value problem. Ramps or pulses or other waveforms generated by Vector play or code in a mod file do not need cvode.event or cvode.re_init.
4) How do I know how gentle of a slope is required to avoid the error?
It's always an empirical question.
5) Is there an easy way (e.g. an existing function) to implement a ramp..?
Vector play with interpolation allows you to do it entirely in hoc. Could also be done with NMODL.
nikkip
Posts: 28
Joined: Mon Jun 24, 2013 6:09 pm

Re: CVODE: record, play, events

Post by nikkip »

So, I'm still getting an integration error (below) at the time of my stimulus onset, even though I'm using dt/1000 for 5*dt before and after the time point when the stimulus turns on. Any other suggestions? All I'm trying to do is apply a pulse of extracellular current (i.e. "pulse" of Ve) to get conduction speed and extracellular strength-duration data using CVODE.

Code: Select all

IDA initialization failure, weighted norm of residual=95.2975
err=-6
        )
      , 1 nrniv: variable step integrator error
1  in MRGaxon_unmyelinated_StrengthDuration_v7.hoc near line 688
1  {pc.runworker()}
                 ^
1 fadvance(1 stimul_extracellular(2.01-2e+080.51 find_thresh(12.01200.501 ParallelContext[0].runworker(application called MPI_Abort(MPI_COMM_WORLD, -1) - process 1
, )
    , , , , , )
  )
[0]0:Return code = 0, signaled with Quit
[0]1:Return code = 255

Thanks!

Nikki
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: CVODE: record, play, events

Post by ted »

I guess I'll have to see this myself in order to find out what is happening and what might be done about it. Can you zip up just the hoc, ses, and mod files that are necessary to reproduce the problem and email that to
ted dot carnevale at yale dot edu
?
nikkip
Posts: 28
Joined: Mon Jun 24, 2013 6:09 pm

Re: CVODE: record, play, events

Post by nikkip »

A few things I figured out, in case they're helpful for future readers.

myvec.record(…) - Vector class record

1) I’m doing a parameter sweep to get strength duration data. I’m sweeping the pulse width and using a binary search algorithm to find threshold at each pulse width. For some of my vectors, I call myvec.record() once outside my “run” procedure, so it doesn’t get recalled for each parameter value. Thus, even though I haven’t emptied my recording vector, nor have I called myvec.record() at the start of each sim, it seems to work.

2) Here are a few examples of how to use myvec.record(…) :

Record the current of an IClamp point process object:

Code: Select all

create node[axonnodes]
…
dt_fixed = 0.001		// [ms]
delay = 2			// [ms]
myamp = 20			// [nA]
myPW = 0.1			// [ms]
objref istim_vec
istim_vec = new Vector()
node[x_elec_ind]{
	stim = new IClamp()
	stim.loc(.5)
	stim.del = delay
	stim.amp = myamp
	stim.dur = myPW
}
istim_vec.record(&stim.i,dt_fixed)
Record time vector:

Code: Select all

dt_fixed = 0.001		// [ms]
objref t_vec
t_vec = new Vector()
t_vec.record(&t,dt_fixed)
Record transmembrane voltage:

Code: Select all

create node[axonnodes]
…
dt_fixed = 0.001		// [ms]
objref Vm_vec[axonnodes]
for i = 0, axonnodes - 1 {
	Vm_vec[i] = new Vector()
	Vm_vec[i].record(&node[i].v(0.5),dt_fixed)
}
myvec.play(…) - Vector class play
In case this example is useful for anyone else. The following code sets up the extracellular voltages for an extracellular current point source that turns on at t=delay and turns off at t=delay+myPW:

Code: Select all

objref Ve_vec[axonnodes]
for i = 0, axonnodes - 1 {
	Ve_vec[i] = new Vector(t_vec_play.size(),0)
}

for i = 0, axonnodes-1 {
	for J = 0, t_vec_play.size()-1 {
		// compute Ve at each node at each time point
		if ((t_vec_play.x[J] > delay) && (t_vec_play.x[J] < delay + myPW)) {
			Ve_vec[i].x[J] = (10^-5)*rhoe*1/(4*PI*sqrt((x_axon.x[i]-x_elec)^2 + (myzelec)^2))
		}
	}
	Ve_vec[i].play(&node[i].e_extracellular(0.5), t_vec_play)
}
Adaptive integration getting stuck
You can allow vector.play() to “interpolate” in time rather than just playing a vector with a fixed time step by setting the third argument to be non-zero:

Code: Select all

sourcevec.play(&targetvar,timevec,1)
That being said, it seems that with some models, it is not possible to apply an extracellular stimulation pulse with CVODE (i.e. adaptive time step) because the integrator gets stuck at the pulse onset. You can diagnose this issue (at least partially) by printing out the time after every fadvance(). So in my case, my stimulus pulse started at t=2ms, so it would repeatedly print 2 before it crashed, and gave a cvode-related error.

Hope that some of that is helpful to someone!

EDIT (June 16, 2014): Vector class record, not CVODE class record.
Last edited by nikkip on Mon Jun 16, 2014 10:19 am, edited 2 times in total.
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: CVODE: record, play, events

Post by ted »

The examples are of the Vector class's record method. Did you mean to include an example of the use of cvode.record?
nikkip
Posts: 28
Joined: Mon Jun 24, 2013 6:09 pm

Re: CVODE: record, play, events

Post by nikkip »

ted wrote:The examples are of the Vector class's record method. Did you mean to include an example of the use of cvode.record?

Oops! No, I never ended up using cvode.record(). Thanks for catching that.
MBeining

Re: CVODE: record, play, events

Post by MBeining »

I have a question related to that.
Since I use a multistep IClamp I play the amps into my IClamp with a play vector.
However that causes my time vector in cvode mode to have two values which are the same.

E.g. if a step is at t = 55 ms, my record time vector has two entries which are exactly 55 ms. And when I record the current of my IClamp, it has at those two indices the entries:
[i0, i1]
where i0 is the amplitude pre-55ms and i1 is the amplitude post-55ms.

Is there any possibility to let cvode.record only record the second entry? Since the amplitude changes exactly at 55ms to i1, it should only record i1 at 55ms and not i0 and i1 ;-)

Thank you in advance.
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: CVODE: record, play, events

Post by ted »

First, it's not cvode.record--we're discussing Vector.record. Second, the answer to your question will be obvious if you think about it. Or if you read the Programmer's reference entries on adaptive integration. If that doesn't make clarity immediately crystallize in your mind, try answering this question: how would your suggestion allow one to disambiguate a sawtooth wave from a series of rectangular pulses?
Post Reply