www.neuron.yale.edu

Posted: **Wed Sep 08, 2010 12:26 pm**

I would like to do several simulations of a single neuron in current clamp, for which the input is read from a file (stim1.dat, stim2.dat, etc; a different one for each turn). I want the simulations to be a follow-up of each other, so the initialization of the next turn should be the same as the last state of the previous turn. After defining the neuron and its geometry I was working with the following code:

Code: Select all

objref r,f, svstate
strdef filename, line
f = new File()
r = new Vector()
svstate = new SaveState() 

proc init() {
  finitialize(v_init)
  svstate.restore()
  t = 0 // t is one of the "states"
  if (cvode.active()) {
    cvode.re_init()
  } else {
    fcurrent()
  }
  frecord_init()
}


for(i = 1; i <= 2; i = i + 1) {
  sprint(filename, "stim%d.dat", i)
  f.ropen(filename)
  r.scanf(f)
  r.play(&IClamp[1].amp, dt)   
  init()
  run()  
  svstate.save()
  f.close()
}

However, this gives me the following error:

Code: Select all

opening file: stim1.dat
starting run 1
running custom init
SaveState warning: 3 sections exist but saved 0
nrniv: SaveState: Stored state inconsistent with current neuron structure
 in D:/Neuron/dendtc/tc3_cc_hcur_loop.hoc near line 394
 }
  ^
SaveState[0].restore(        )
init(      )

I can understand that for the first time, there was no SaveState yet, so this could be a problem. However, adding 'svstate.save()' after 'svstate = new SaveState()' makes the first simulation give wrong results and gives another error in the second run:

Code: Select all

'Failed assertion nprs_==prl->count() at line 574 of file/home/hines/neuron/releases/nrn-7.1/src/nrniv/savstate.cpp'

Does anyone know how to do several consecutive simulations like this?

Posted: **Mon Sep 13, 2010 3:21 pm**

You want the initialization for the first run to do what unmodified proc init() does; only initializations for the subsequent runs may use SaveState restore().

The "brute force" way to do this is with inline code, i.e. a program organized like so:
(presumes biology and instrumentation--stimulators, graphs, Vector recording etc.--have all been done prior to this point)

Code: Select all

objref r,f, svstate
strdef filename
svstate = new SaveState() 
f = new File()
r = new Vector()
r.play(&IClamp...) // do this only once,
  // and do it before any SaveState save or restore

filename = "stim1.dat"
f.ropen(filename)
r.scanf(f)
run() // first pass uses NEURON's built-in proc init()
svstate.save()

// now declare your custom proc init()
// just as in the example you provided

proc init() { . . .

// finally your for loop, but starting with i = 2

for (i = 2; . . .

That isn't elegant but it should work.

A more sophisticated strategy would employ a single custom init procedure that expects an argument that specifies whether or not to to use SaveState restore()--in other words, whether or not the first run in a series of runs has already been executed. Something like the following should work:

Code: Select all

objref r,f, svstate
strdef filename
f = new File()
r = new Vector()
r.play(&IClamp...) // do this only once
svstate = new SaveState() 

proc myinit() {
  finitialize(v_init)
  if (!$1) { // this is not the first run in a series
    svstate.restore()
    t = 0 // t is one of the "states"
    if (cvode.active()) {
      cvode.re_init()
    } else {
      fcurrent()
    }
    frecord_init()
  }
}

proc mystdinit() {
  realtime=0
  startsw()
  setdt()
  myinit($1)
  initPlot()
}

proc myrun() {
  mystdinit($1)
  continuerun(tstop)
}

proc batchrun() { local firstrun
  firstrun = 1
  for(i = 1; i <= $1; i = i + 1) {
    sprint(filename, "stim%d.dat", i)
    f.ropen(filename)
    r.scanf(f)
    f.close()
    myrun(firstrun)
    svstate.save()
    firstrun = 0
  }
}

Then you can call batchrun() with an integer argument that specifies how many passes to make.

Comments on this implementation:

batchrun() contains the main administrative loop--which launches a series of simulations. During the first pass, the firstrun is 1. The value of this variable controls whether or not states are restored from a SaveState object. This means the value of firstrun must be passed to the custom init procedure. I could have made firstrun be a variable with global scope, but it is safer to pass its value as an argument. In NEURON's standard run system the call chain from run to init is
run() -> stdinit() -> init()
so passing firstrun's value as an argument requires customization of not only init but also stdinit and run. Rather than overwrite the standard run time library's run, stdinit, and init, I decided to add my own procs with similar (but unique) names. This leaves the standard run system remains intact, so I can use the RunControl panel for exploratory/development/debugging purposes.

One more comment--note that it is only necessary to call Vector play() once. Vector play() remains in effect until the Vector is destroyed, which it is not (you're just reading new values from a series of files). Also, not knowing the implementational details of SaveState save(), I would hesitate to call Vector play() _between_ a SaveState save() and a subsequent restore().

Posted: **Tue Sep 14, 2010 9:03 am**

Dear Ted,

Thank you so much for your fast reply. I tried what you said, but unfortunately, something goes wrong with the vector-play and savestate combination:

For the first method I changed my code to

Code: Select all

// Variable current electrode
objectvar El2             
soma El2 = new IClamp(0.5) 
El2.del = 0
El2.dur = tstop-tstart   

objref r,f, svstate
strdef filename
svstate = new SaveState()
f = new File()
r = new Vector()
r.play(&El2.amp, dt)    // do this only once,
  // and do it before any SaveState save or restore

filename = "stim1.dat"
f.ropen(filename)
r.scanf(f)
run() // first pass uses NEURON's built-in proc init()
svstate.save()
f.close

// now declare your custom proc init()
// just as in the example you provided

proc init() {
  finitialize(v_init)
  svstate.restore()
  t = 0 // t is one of the "states"
  if (cvode.active()) {
    cvode.re_init()
  } else {
    fcurrent()
  }
  frecord_init()
}

// finally your for loop, but starting with i = 2
for(i = 2; i <= 2; i = i + 1) {
  sprint(filename, "stim%d.dat", i)
  f.ropen(filename)         
  r.scanf(f) 
  init()
  run() 
    
  svstate.save()
  f.close()
}

as you said.
This way the simulation was indeed started from the last state of the previous run each time, but only the first 'stim'-file was played into the current clamp, the second run had no time-dependent stimulation. The right file was read into the vector r, but not played into the electrode: only the last value was used. If I added 'r.play(&IClamp[1].amp, dt) ' in the for loop after r.scanf(f), or after init(), it did not change this, but now an error in savestate occured (the savestate.save()).

The batchrun has the same problem: savestate works, but only the last value of the 'second' vector r is used.

Somehow, savestate.save() and r.play(&El2.amp, dt) interfere with one another, but I cannot really see how.

Posted: **Tue Sep 14, 2010 10:00 pm**

I suppose it is possible that the SaveState object contains not only the states of the "represented biology" but also whatever r contained before SaveState save() was called. If that is the case, it is necessary to defer reading the stimulus waveform into r until after return from init().

Code: Select all

proc myrun() {
  mystdinit($1)
  f.ropen(filename)
  r.scanf(f)
  f.close()
  continuerun(tstop)
}

proc batchrun() { local firstrun
  firstrun = 1
  for (i = 1; i <= $1; i = i + 1) {
    sprint(filename, "stim%d.dat", i)
/*
    f.ropen(filename)
    r.scanf(f)
    f.close()
*/
    myrun(firstrun)
    svstate.save()
    firstrun = 0
  }
}

Posted: **Thu Sep 16, 2010 11:54 am**

Thank you again for the fast reply. Unfortunately, the simulations still have the same problem.

Now the first time I do a batchrun (for instance batchrun(2) ), also the first vector/file is not played into the stimulus, and during the second run only the first value of the first vector is used.

The second time I do a batchrun, he first vector/file is played correctly, but not the second one. Again, the second run only uses the first value of the first vector.

I think the files are read into the vector correctly (I checked it using r.x() ), but somehow they are not played into the stimulus.

When I take away

Code: Select all

svstate.save()
firstrun = 0

The first time I do a batchrun, the first vector is not played, the second is. The second time I do a batchrun they are both played correctly, but of course they are starting in the same state each time, instead of using the last state of the previous run.

So somehow both things can work (reading the files and playing them consecutively, using savestate), but not at the same time.

Posted: **Thu Sep 16, 2010 1:49 pm**

fleurzeldenrust wrote:So somehow both things can work (reading the files and playing them consecutively, using savestate), but not at the same time.

More correctly, if there is a way to do it, you and I haven't stumbled into it yet. I'll have more to post soon.

Posted: **Thu Sep 16, 2010 2:19 pm**

Well, I'm stuck. I can get the first pass to work, but subsequent passes result in no stimulus, even though the data are read from the file into r. Inserting
r.play(&stim.amp, DT)
right after the statement that reads the data into r, so that it is executed on each pass, doesn't help.

Posted: **Fri Sep 17, 2010 2:10 pm**

Michael Hines suggested calling SaveState restore() with an argument of 1, because it doesn't clear the event queue, and sure enough it works. Vector play works properly in all runs, not just the first, and each run's initial conditions are the conditions at the end of the previous run (except for t, which is always reset to 0).

Here's all of the code that I used to test this:

init.hoc

Code: Select all

load_file("nrngui.hoc")
load_file("rig.ses")

for ii=0,1 Graph[ii].exec_menu("Keep Lines")

load_file("newctl.hoc")

rig.ses
rig.ses is a session file that recreates the following items:

RunControl panel with Tstop 5 ms, dt 0.025 ms, Points plotted/ms 40.
SingleCompartment (NEURON Main Menu / Build / single compartment) with "pas" selected, so it creates a single compartment model with surface area 100 um2 and pas mechanism inserted.
PointProcessManager configured as an IClamp.
Two Graphs, both "Voltage axis" so x coordinates are integer multiples of dt (one graph plots v(0.5), the other shows IClamp[0].i)

newctll.hoc

Code: Select all

objref stim
stim = IClamp[0]
stim.del = 0
stim.dur = 1e9
stim.amp = 0

DT = 0.025 // stimulus sampling interval

objref r,f, svstate
strdef filename
f = new File()
r = new Vector()
r.play(&stim.amp, DT) // do this only once

svstate = new SaveState()

proc myinit() {
  finitialize(v_init)
  if (!$1) { // this is not the first run in a series
    svstate.restore(1) // don't clear the event queue
    t = 0 // t is one of the "states"
    if (cvode.active()) {
      cvode.re_init()
    } else {
      fcurrent()
    }
    frecord_init()
  }
}

proc mystdinit() {
  realtime=0
  startsw()
  setdt()
  myinit($1)
  initPlot()
}

proc myrun() {
  mystdinit($1)
  continuerun(tstop)
}

proc batchrun() { local firstrun  localobj mtrx // data were saved from a "picked" Vector
  mtrx = new Matrix(1,2)
  firstrun = 1
  for(i = 1; i <= $1; i = i + 1) {
    sprint(filename, "stim%d.dat", i)
    f.ropen(filename)
    mtrx.scanf(f)
    f.close()
    mtrx.getcol(1, r) // doesn't destroy the existing r object, unlike r = mtrx.getcol(1)
    myrun(firstrun)
    svstate.save()
    firstrun = 0
  }
}

Note that my data files were created by "picking" the time course of current pulse waveforms from a "Voltage axis graph"* (select the graph menu's Pick Vector, then click on the trace to be "picked") and saving them to files with
NEURON Main Menu / Vector / Save to File
This produced files that begin with a two-line header, which is followed by pairs of xy values, e.g.

Code: Select all

label:IClamp[0].i
201
0       0
0.025   0
0.05    0
 . . .

I edited these files, removing the first line and appending the number of columns (entries per line) to the second line, like so

Code: Select all

because mtrx.scanf(f) expects the first line to contain the number of "rows and columns" of data that are contained in a file.

*--in a "Voltage axis graph" the x coordinates are integer multiples of dt.

Posted: **Thu Sep 23, 2010 6:38 pm**

Sorry for the late reply, I have been travelling to the US and a meeting. Thank you for the effort, I will try the solution as soon as I can, I will post the results soon.

Posted: **Sat Sep 25, 2010 2:49 pm**

Works like a charm, thank you very much!

Posted: **Sun Sep 26, 2010 10:18 am**

Thank you for using NEURON in your research. If this work leads to publication(s), please be sure to cite NEURON (see How to cite NEURON viewtopic.php?f=22&t=73), and also let me know so we can add your publication(s) to the Bibliography page http://www.neuron.yale.edu/neuron/stati ... ednrn.html.

www.neuron.yale.edu

Simulation divided into short ones using SaveState

Simulation divided into short ones using SaveState

Re: Simulation divided into short ones using SaveState

Re: Simulation divided into short ones using SaveState

Re: Simulation divided into short ones using SaveState

Re: Simulation divided into short ones using SaveState

Re: Simulation divided into short ones using SaveState

Re: Simulation divided into short ones using SaveState

Re: Simulation divided into short ones using SaveState

Re: Simulation divided into short ones using SaveState

Re: Simulation divided into short ones using SaveState

Re: Simulation divided into short ones using SaveState