Question about proc batchrun()

[erhervert]

I'm working along with a cell model taken from Muller et al.,2009. I'm trying to execute multiple simulations varying the syn.onset parameter using the batch run procedure according with http://www.neuron.yale.edu/phpbb/viewto ... =28&t=2333. The problem is that when I run the simulation with the command go(), the model just runs once and the syn.onset parameter does not vary. I'd like to execute it 5 times with syn.onset increments of 10 ms each, could you tell me if I'm making mistakes please?

This is what I have:

Code: Select all

///////////////////////////
//* simulation control *//
///////////////////////////

dt = 0.025
tstop = 50
v_init = -65

proc batchrun() { local i, tmp
  tmp = syn.onset
  for i=0,$0-5 {
    syn.onset = tmp*i
    run()
  }
}

proc initialize() {
  finitialize(v_init)
  fcurrent()
}

proc integrate() {
  g.begin()
  while (t<tstop) {
    fadvance()
    g.plot(t)
  }
  g.flush()
}

proc go() {
  initialize()
  integrate()
}

Thanks in advance

[ted]

Here are some hints and questions that will help you answer your own question.

But first this question: why bother to reinvent NEURON's standard run system? And a small subset of it, at that. I refer to this excerpt from your code:

Code: Select all

proc initialize() {
  finitialize(v_init)
  fcurrent()
}

proc integrate() {
  g.begin()
  while (t<tstop) {
    fadvance()
    g.plot(t)
  }
  g.flush()
}

proc go() {
  initialize()
  integrate()
}

If you simply
load_file("nrngui.hoc")
at the start of your program (or if you don't want to see the NEURON Main Menu toolbar, just
load_file("stdgui.hoc")), you can launch a single run with just one command:
run()
and it will automatically take care of initialization and updating Graphs. In other words, run() will take care of everything that's in your code's initialize(), integrate(), and go() procedures--and its easy to type and doesn't clutter up your source code.

Are you following an example you have seen elsewhere? If so, please let me know so that I can find out if it is possible to do something about it (i.e. something that will help prevent others from spending time reinventing the standard run system).

"Well, I have cobbled up a Graph that isn't updating from one step to the next."

Fine. If you're using the standard run sytem, you can just append your non-updating Graph to a graphList. graphList[0] is a good choice if you want points to be plotted vs. t.
graphList[0].append(g)
will do the job. Now run() will launch a simulation and your special Graph will update properly.

when I run the simulation with the command go(), the model just runs once

So walk through your code to discover what happens when you execute go(). Pretend you're the computer, read each statement line by line, and imagine that you do what it says.
go() calls initialize()
initialize() calls finitialize(v_init), then calls fcurrent() and returns
go() next calls integrate()
integrate does what?
Did you see anything that calls batchrun()?

And about batchrun()--

Code: Select all

proc batchrun() { local i, tmp
  tmp = syn.onset
  for i=0,$0-5 {
    syn.onset = tmp*i
    run()
  }
}

have you tried calling this yourself from the oc> prompt? That is, at the oc> prompt type
batchrun()
and see what happens. Did it work as expected? Or were you going to call it with a numerical argument, and if so, what? Maybe it would be useful to comment out the run() line and insert a print statement so you can get a better idea of what will happen when batchrun() is called. That is, change the for loop to

Code: Select all

  for i=0,$0-5 {
    syn.onset = tmp*i
//    run()
print "i ", i, "   syn.onset ", syn.onset 
  }

[erhervert]

Hello Ted!

To give an answer to the first question I'm working with a model published in the following paper --> Hines ML, Davison AP, Muller E. NEURON and Python. Front Neuroinform. 2009;3:1. Epub 2009 Jan 28 and I really don't know why they used a modified version of the standard run system but I was executing it through the command go().

I added the line load_file("nrngui.hoc") at the start of my program and the graphList[0].append(g) as you said and it worked very good. Now the graph updates each time I type the command run() at the oc> prompt.

The next thing I do was to test what happen with the batchrun() process when I typed batchrun() directly into the oc> prompt and I got and error related to the arguments. So I commented out the run() line and insert the print statement into the loop just as you suggest and it returned the same error but inmediatly I added a numerical value to the batchrun() and run it again. Surprisingly, it returned a list of values from 0 to 5 (because I used the number 5 as numerical argument) that looks like the following:

i0 syn.onset 10
i1 syn.onset 10
i2 syn.onset 10
i3 syn.onset 10
i4 syn.onset 10

Then I modified the proc batchrun() to get something like this at the end of the program:

Code: Select all

proc batchrun() { local i, tmp
  tmp = syn.onset 
  for i=1, $1-1 {
    syn.onset = tmp*i
    run()
  }
}

And now when I execute my program using the batchrun() with a numerical argument it does exactly what I had in mind. I can run my model n times varying the syn.onset parameter in each iteration to get a beautiful plot that shows the results of the simulation, everything is perfect right now.

This is the full code:

Code: Select all

/////////////////////////
/* model specification */
/////////////////////////
/////// topology ////////
{load_file("nrngui.hoc")}
create soma, apical, basilar, axon
connect apical(0), soma(1)
connect basilar(0), soma(0)
connect axon(0), soma(0)
/////// geometry ////////
soma {
  L = 30
  diam = 30
  nseg = 1
}

apical {
  L = 600
  diam = 1
  nseg = 23
}

basilar {
  L = 200
  diam = 2
  nseg = 5
}

axon {
  L = 1000
  diam = 1
  nseg = 37
}
////// biophysics //////
forall {
  Ra = 100
  cm = 1
}

soma { 
  insert hh
}

apical {
  insert pas
  g_pas = 0.0002
  e_pas = -65
}

basilar {
  insert pas
  g_pas = 0.0002
  e_pas = -65
}

axon {
  insert hh
}

//////////////////////////
/*    Instrumentation   */
//////////////////////////

///// synaptic input /////

objref syn
soma syn = new AlphaSynapse(0.5)
syn.onset = 10
syn.tau= 0.1
syn.gmax = 0.05
syn.e = 0
/// grphical display ///
objref g
g = new Graph()
g.size(0,50,-80,40)
g.addvar("soma.v(0.5)", 1, 1, 0.6, 0.9, 2)
graphList[0].append(g)

///////////////////////////
//* simulation control *//
///////////////////////////
dt = 0.025
tstop = 50
v_init = -65

proc batchrun() { local i, tmp
  tmp = syn.onset 
  for i=1, $1-1 {
    syn.onset = tmp*i
    run()
  }
}

Thanks for your help Ted

[ted]

To give an answer to the first question I'm working with a model published in the following paper --> Hines ML, Davison AP, Muller E. NEURON and Python. Front Neuroinform. 2009;3:1. Epub 2009 Jan 28

Well, that's a fairly authoritative source.

I really don't know why they used a modified version of the standard run system

I'm not sure myself--maybe because they wanted a "minimal" run time system. But I wouldn't advise anyone to replicate that particular aspect of their code, on the principles of "less code is better code" and "don't reinvent the wheel".

I added the line load_file("nrngui.hoc") at the start of my program and the graphList[0].append(g) as you said and it worked very good. Now the graph updates each time I type the command run() at the oc> prompt.

That's good.

So you worked through those suggestions, ran into some problems, read some documentation and discovered the cause of the problems and how to fix them. Successful debugging. The next time a problem occurs--and problems have a way of happening during code development--you'll be more likely to be able to discover and fix the cause.

And now when I execute my program using the batchrun() with a numerical argument it does exactly what I had in mind.

Well, if I call batchrun(2) I only see one run, which corresponds to i == 1. Changing
for i=1, $1-1 {
to
for i=1, $1 {
makes batchrun(2) generate two runs, which seems to be an improvement.

[erhervert]

Well, if I call batchrun(2) I only see one run, which corresponds to i == 1. Changing
for i=1, $1-1 {
to
for i=1, $1 {
makes batchrun(2) generate two runs, which seems to be an improvement.

That's true. I also added a number generator to vary the value of syn.onset randomly. This is what I got:

Code: Select all

objref r
r = new Random()
r.normal(.5, .05)


proc batchrun() { local i, index
  for i=0, $1 {
    index = r.repick()
    syn.onset = (((10-5)*index)+5)
    run()
  }
}

The only thing is that when I load the file, there are some numbers (1, 1, 1, 0.62880253) below the NEURON header. Do you know why those numbers appear? May I remove them?

NEURON -- Release 7.2 (562:42a47463b504) 2011-12-21
Duke, Yale, and the BlueBrain Project -- Copyright 1984-2008
See http://www.neuron.yale.edu/credits.html

1
1
1
0.62880253

Thanks a lot!

[ted]

Some statements automatically return values that are printed to the terminal. Simple example:
oc>1+2
3
To suppress this, surround the statement with a pair of curly brackets. Example:
oc>{1+2}
(prints nothing)

Always remember that the normal distribution extends to +- infinity. For your particular use of random numbers, a large negative number would be meaningless. You might want to prevent that from happening. One way would be to replace

Code: Select all

    index = r.repick()
    syn.onset = (((10-5)*index)+5)

with
syn.onset = tsynon()
where tsynon is a func that returns values that are >= 0. A simple algorithm to use in such a function is

Code: Select all

ton = -1
while (ton<0) {
  ton = some expression involving r.repick()
}
return ton

The while loop exits when ton is >= 0.

[erhervert]

You are right. I implement this in my program and it works fine.

Code: Select all

objref r
r = new Random()
{r.normal(.5, .05)}

objref avg

func tsynon() {
  ton = -1
  while (ton<0) {
    ton = (((10-5)*r.repick())+5)
  }
  return ton
}

proc batchavg() { local i
  for i=0, $1 {
    syn.onset = tsynon()
    run()
  }
  avg = new Vector()
  if (i==1) {
    tmp = avg.c()
  } else {
      avg.add(tmp)
    }
  avg.div($1)
}

I also found that statements that were printing stuff below the NEURON header and I fixed them. Now they are not printing anything and the program is still working fine.

Code: Select all

{g.size(0,50,-80,40)}
{g.addvar("soma.v(0.5)", 1, 1, 0.6, 0.9, 2)}
{graphList[0].append(g)}

And

Code: Select all

{r.normal(.5, .05)}

Thank you very much!