First run of simulation acting differently

[aluchko]

I'm having a lot of trouble understanding why the first run of this loop is producing different results than the subsequent runs, the issue seems to come up when there's a single fadvance() after the initialization and the SEClamp is somewhat configured. I can't really understand what I'm doing wrong.

Code: Select all

dt=0.1

create soma
access soma

soma {
        nseg = 1
        diatstop = i
        L = 20
        diam=20
        Ra = 150
        insert hh
        ek = -103
        insert pas
        g_pas = 0.00003334
}

objectvar stim
soma stim = new SEClamp(0.5)
tstop = 200

finitialize(-110)
stim.dur1 = 50
//stim.dur2 = 400 // Uncomment this line and the issue goes away
stim.amp1 = -110
stim.amp2 = -40
fadvance()

for i=0,2 {
  print "Run ", i
  finitialize(-110)
  stim.dur1 = 50
  stim.dur2 = 400
  stim.amp1 = -110
  stim.amp2 = -40
  while (t < tstop) {
    if (t > tstop - .2) {
      print soma.v, stim.amp2
    }     
    fadvance()
  }
}

Outputs

Code: Select all

Run 0 
-69.484068 -40 
-69.484068 -40 
Run 1 
-43.717401 -40 
-43.717401 -40 
Run 2 
-43.717401 -40 
-43.717401 -40

[ted]

It is generally not a good idea to interleave instrumentation code (especially statements that alter SEClamps, IClamps, or other signal sources) and simulation control code (statements that initialize or launch simulations). The problem you encountered happens because the code calls finitialize too soon--specifically, before one of the simulation parameters is changed (i.e. before the duration of one of the pulses that is to be generated by the SEClamp). The solution is to defer calling finitialize until just before the main computational loop (the "while" loop).

finitialize() does many things. It is best to avoid writing one's own simulation control code; instead, take advantage of NEURON's own run control system. For this particular example, it would be better to

Code: Select all

load_file("nrngui.hoc") // or load_file("stdgui.hoc") if you don't want to see any GUI stuff
// now all of NEURON's libraries are available
 . . . model specification code here . . .
 . . . instrumentation code here . . .
// finally the simulation control code
v_init = -110
tstop = 200
steps_per_ms = 10 // the standard run system forces dt to being an integer fraction of 1/steps_per_ms
dt = 1/steps_per_ms // so you are sure to get the dt value you want
run()

If there is some code that must be executed at a particular time, e.g. to report the value of some variable, use cvode.event() and an FInitializeHandler to make that happen. Example: after model setup is complete, but before calling run(), the following statements will make something special happen a particular number of times, beginning at a particular start time, and recurring at a particular interval--

Code: Select all

T_EVENT1 = 13 // the first time it should happen
EVENT_INTERVAL = 5 // how often it should repeat
NEVENTS = 2 // how many times it should be repeated
nevents = NEVENTS // nevents is the number of events that are yet to be done

func event_handler() {
  if ($1 > 0) { // act only if something remains to be done
    do_something() // whatever you need to happen
    $1 -= 1
    if ($1 > 0) cvode.event(t + EVENT_INTERVAL, "nevents = event_handler(nevents)") // prepare for the next one
  }
  return $1
}

proc launch_events() {
  nevents = NEVENTS
  cvode.event(T_EVENT1, "nevents = event_handler(nevents)")
}

objref fih
fih = new FInitializeHandler("launch_events()")