Vector play problems

[carl]

I am having trouble getting my own complex SEClamp to work. I load a waveform from a file int o a vector and create an appropriate tvec, and then set SEClamp to play it. But my program bombs with a segmentation fault. Here is what the code looks like:

-----------------------------------------------------------------------------------------------------
// make the vector of data by loading a file

volts_clamp_vec = new Vector()
sprint(tmp_str,"%s/%s", param_dir, volts_clamp_file_name)
volts_clamp_file_ref = new File(tmp_str)
volts_clamp_file_ref.ropen()
volts_clamp_vec.scanf(volts_clamp_file_ref)
volts_clamp_file_ref.close()


// make an appropriate time vector: the file contained data sampled at the rate
// of 100 kHz so that should be the time step

volts_clamp_tvec = new Vector(volts_clamp_vec.size())
volts_clamp_tvec.indgen(0.01)


// this loop confirms that the time and data vectors contain what I expect
// (this works!)

for i=0, volts_clamp_vec.size()-1 {
print "wav(", i, ") t = ", volts_clamp_tvec.x, " / V = ", volts_clamp_vec.x
}

// now construct the voltage SEClamp - as advised in previous posts, just give a duration
// that is much too long
soma volts_clamp_ref = new SEClamp(.5)
volts_clamp_ref.dur1 = 500
volts_clamp_ref.rs=.01

// finally, tell the voltage vector to play into the SEClamp
volts_clamp_vec.play(&volts_clamp_ref.amp1, volts_clamp_tvec)
--------------------------------------------------------------------------------------------------------

But when I run my simulation I get

/usr/local/bin/nrniv: Segmentation violation

What might be the problem? Thanks!

[ted]

You're just in a hurry to get results, and accidentally left something out. The code you wrote
is fine; the problem is the code that you didn't write. Variable names must be declared as
objrefs before they can be used as objrefs.

[carl]

Actually, I keep all my object refs in one "header" file and didn't mention it in my post. Sorry for not being clear on all the details. So this bit of code is loaded before what I mention above:

// for a voltage clamp, if used
objref volts_clamp_ref
objref volts_clamp_vec, volts_clamp_tvec
strdef volts_clamp_file_name
objref volts_clamp_file_ref

As I mentioned, the loop displaying the vector contents (after scanf) runs and gives the expected result. So I have the vector created. And the code creating the SEClamp and calling the vector 'play' function runs. And then some other generic setup for my simulation takes place. And then my simulation runs for like 0.01 ms (perhaps the first CVODE step, but I'm not actually checking how many steps it goes for) BUT then it dies with the segmentation fault. Can you think of anything else that could be the problem?

I few other details are that I'm using CVODE and have secondorder set to 2. Could there be any problems associated with this? I read in the previous posts the VClamp is not CVODE compataible but SEClamp is.

Another question is how do you pass in the parameter for "continuous" is that meant to be written as a striing, as in simply writing that text as a third parameter:

volts_clamp_vec.play(&volts_clamp_ref.amp1, volts_clamp_tvec, continuous)

or is that a variable, like

continous = 1
volts_clamp_vec.play(&volts_clamp_ref.amp1, volts_clamp_tvec, continuous)

I guess that sounds like a stupid question but its not clear from the documentation. I had decided to load my vector with very high frequency sampled data under the assumption that if CVODE takes bigger steps it will just skip over the data for the very short dt I provide. But if this continuous option does what I think it does that would be a better way, i.e. provide data at 20 kHz or so and let it interpolate the missing points. Right?

[ted]

Caught in a logical inconsistency again. If I keep that up, eventually I'll pass the Turing
test . . .

The problem lies not in the code you show, but in the code you don't. Do you get
segmentation faults when you run a simulation with fixed dt?

I'm using CVODE and have secondorder set to 2. Could there be any problems associated with this?

Other than the fact that it's like trying to mount oars on a bicycle, nothing in particular. Have
you read chapter 4 in The NEURON Book? Just the last 3 pages will do, i.e. 86-88.
Either you're using a fixed dt method (the default implicit Euler or one of the
Crank-Nicholson variants) or you're using adaptive integration (aka "variable time step"
or CVODE, both of which are misnomers because (1) it's really variable order AND
variable time step, and (2) it's really CVODES or IDA). If cvode is active, the value of
secondorder is irrelevant.

how do you pass in the parameter for "continuous"

To quote from the documentation

When continuous is 1 then linear interpolation is used to define the values between time points.

To quote from your code,

Code: Select all

volts_clamp_vec.play(&volts_clamp_ref.amp1, volts_clamp_tvec, 1)

I had decided to load my vector with very high frequency sampled data under the assumption that if CVODE takes bigger steps it will just skip over the data for the very short dt I provide.

The integrator doesn't skip over anything. Whether you use a time vector or Vector
play's dt syntax, the integrator will meticulously stop at every specified time and grind
out a new solution point. But, as per the documentation of Vector play,

for variable step methods, unless continuity is specifically requested, the function is a step function.

That is, the transferred variable will be treated as constant during each time step unless
you explicitly request "continuous." If you want the simulation to behave as if the driven
variable is a smooth function of time, interpolation is necessary.

[carl]

Well, when I switch to a fixed step (no CVODE) I can actually run the code without crashing. I set the fixed step to the same as my input vector (i.e. 0.01)

I still have some problems, but I now have a workaround...

1) the Vector.play call I described in previous posts simply does not work - it is ignored and the SEClamp stays put at whatever voltage I initialized it to. What does work is if i explicitly reset the SEClamp's amp1 in my main loop. i.e.

step = 0
while (t < tstop) {
volts_clamp_ref.amp1 = volts_clamp_vec.x[step]
fadvance()
step = step +1
}

The above is not very elegant, but it gets the job done. I don't have any idea why the Vector play method is not working. The call to play is accepted by the interpreter, but nothing ever results.

2) In my experiments with trying to get SEClamp to work I also tried putting in fixed voltges and durations. In this process I found that amp2, dur2, amp3, dur3 in the SEClamp had no effect : only amp1 and dur1 did anything.

3) Not so much a problem as an observation: the simulation is VERY unstable with secondorder=2. If anyone else tries something similar (SEClamping to a high frequency waveform with a large, complex cell) I would say you'd better stick to secondorder=0.

Maybe I should mention that I am running:
NEURON -- Version 5.7 2005-4-21 11:7:7 Main (159)

Perhaps a later version would have better luck? (I mean with the Vector.play and the amp2,dur2 stuf) If so I can bug my sysadmin about upgrading us when he gets back from vacation - for now I think I can run my experiments with the workaround described above.

[ted]

Well, when I switch to a fixed step (no CVODE) I can actually run the code without crashing.

Yep.

1) the Vector.play call I described in previous posts simply does not work - it is ignored and the SEClamp stays put at whatever voltage I initialized it to. What does work is if i explicitly reset the SEClamp's amp1 in my main loop. i.e.

step = 0
while (t < tstop) {
volts_clamp_ref.amp1 = volts_clamp_vec.x[step]
fadvance()
step = step +1
}

The above is not very elegant, but it gets the job done. I don't have any idea why the Vector play method is not working.

Time for me to ask a few questions. Your answers will help me decide how best to
answer your questions. They will also be very helpful in terms of guiding the
development of new documentation.

Before you started wrestling with Vector play, were you using NEURON's standard run
system (executing simulations by simply calling run(), which initializes the model and then
launches a simulation that continues until t=tstop)? Or were you using your own main
computational loop even then?
Have you ever tried Vector play before, with any model, and gotten it to work?
Does your model involve anything other than NEURON's built-in mechanisms, i.e. pas,
hh, extracellular? If so, are those mechanisms specified with NMODL code in mod files,
or with the Channel Builder? Are there any ion accumulation mechanisms, pumps, or
buffers? I ask because running with fixed dt worked, more or less, but not with cvode.

3) Not so much a problem as an observation: the simulation is VERY unstable with secondorder=2. If anyone else tries something similar (SEClamping to a high frequency waveform with a large, complex cell) I would say you'd better stick to secondorder=0.

If secondorder is 0, the integrator is a fully implicit Euler method which has excellent
stability. When secondorder is 1 or 2, the integrator is a Crank-Nicholson method. CN
is susceptible to spurious (albeit damped) oscillations when system equations are stiff
(involve time constants that span a wide range). Voltage clamping a cell forces the
effective membrane time constant in the vicinity of the clamp to be extremely short.
Consequently dt must also be very short, or else the simulation will be marred by
oscillations. Bottom line: if you're using any kind of voltage clamp--whether SEClamp or
VClamp--don't use CN. Use implicit Euler (secondorder = 0) or cvode.

Perhaps a later version would have better luck? (I mean with the Vector.play and the amp2,dur2 stuf). If so I can bug my sysadmin about upgrading us

Has nothing to do with it. However, there are many other reasons why you should get
the latest version (bug fixes and improved functionality of many tools). If you are using
MSWin, be sure to get the latest alpha version from
http://www.neuron.yale.edu/ftp/neuron/versions/alpha/
because it contains a critical fix for modlunit.

[carl]

were you using NEURON's standard run system (executing simulations by simply calling run(), which initializes the model and then launches a simulation that continues until t=tstop)? Or were you using your own main computational loop even then?

My own main computational loop - I have some special steps that must be taken on every fadvance (exporting details of all the membrane currents.) Is there a better way to accomplish this?

Have you ever tried Vector play before, with any model, and gotten it to work?

No

Does your model involve anything other than NEURON's built-in mechanisms, i.e. pas,
hh, extracellular? If so, are those mechanisms specified with NMODL code in mod files,
or with the Channel Builder?

Several mechanisms specified with NMODL code.

Are there any ion accumulation mechanisms, pumps, or
buffers?

Yes: calcium accumulation with the standard cadif.mod

[ted]

[quote="carl"]My own main computational loop{/quote]
That might be why Vector play isn't working. play and record take advantage of the
standard run system's internals to make sure that the vectors are properly connected
to the variables, and to ensure that values are transferred correctly as the simulation
advances.
[quote]I have some special steps that must be taken on every fadvance (exporting details of all the membrane currents.) Is there a better way to accomplish this?{/quote]
Do you just need to capture them all to Vectors? That should be doable with Vector
record or CVode record--see
http://www.neuron.yale.edu/neuron/stati ... tml#record
--without any modification to the standard run system. Special stuff that must be done on
a per-step basis can be accomplished by customizing proc advance(). The standard
run system's advance() is

Code: Select all

proc advance() {
  fadvance()
}

If you have to chew on results, or stuff them away anywhere other than in Vectors, just
put your special code in a proc of its own. Then, somewhere near the end of your other
code, put this so that it will replace the proc advance() that was defined when nrngui.hoc
or noload.hoc was loaded:

Code: Select all

proc advance() {
  fadvance()
  myspecialproc()
}

Some complex models run nicely with fixed dt but crash and burn with cvode. This tends
to happen when there are mechanisms specified with NMODL that have peculiar
nonlinearities. It can also happen when calcium accumulation is involved, unless special
care is taken when setting error tolerances. The Atol Scale Tool (part of the Variable
Step Control) can do a test run and then, at the click of a button, customize the error
tolerances of individual variables for you. This can make a lot of otherwise recalcitrant
models work nicely with cvode. That tool also makes it easy to manually adjust individual
tolerances, should it be necessary.