Multiple runs in NEURON result in different results

Anything that doesn't fit elsewhere.
Post Reply
Skyfire79

Multiple runs in NEURON result in different results

Post by Skyfire79 » Fri Jan 25, 2013 4:36 pm

Hi, Ted,

I want to study parameter variation of my NEURON model and hope it can run automatically for multiple times (each time with new values of certain parameters) and save the data in different directories for later analysis. The problem now I have is that even with identical parameters (nothing change in the code), the data I got from two different runs are different. What I did is simple:

proc run() {
running_ = 1
stdinit()
continuerun(tstop)
}

proc multi_run() {
for mi = 1, 2 {
// set_parameter( ) // New parameter will set in this function
run( )
save_data(mi)
}
}

multi_run()

I also run the simulation twice using the "Init & Run" button in the NEURON Control Menu, and the results I got from these two runs are still different: same as using the multi_run() function I wrote. When I compared the voltages I recorded from these two runs, some neurons have the identical voltages, some don't. Here is an example for a neuron with different voltage:

1st run:
-70
-69.9829
-69.8064
-69.5509
-69.2699
-68.9842
-68.7011
-68.4232
-68.1517
-67.8904
-67.6496
-67.4391
-67.2598
-67.1078
-66.978


2nd run:
-70
-69.9829
-69.8064
-69.5509
-69.2694
-68.9393
-68.5515
-68.1485
-67.7518
-67.3786
-67.0432
-66.7491
-66.4922
-66.2665
-66.0666


That is quite surprising to me since it means the "Init & Run" button in the Control Menu could produce different results when pressed for multiple times, even if nothing in the code change. I am also very confused about this since the run() should initialize all the variables back to the same initial conditions. Please advise and thank you very much.

Best,
Guoshi

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

Re: Multiple runs in NEURON result in different results

Post by ted » Fri Jan 25, 2013 5:10 pm

This is usually a symptom of incorrect initialization, typically (but not always) caused by a poorly written mod file. It is a good thing that you discovered this; many modelers (not just NEURON users!) don't bother to check for run-to-run reproducibility, and aren't even aware that initialization is important. To tell you how to fix the problem, I will have to be able to reproduce it. If you zip up just enough code for me to run your simulation--hoc files, ses files, mod files, but no .o, .c, or .dll please--and email it to
ted dot carnevale at yale dot edu
I will tell you what needs to be done.

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

Re: Multiple runs in NEURON result in different results

Post by ted » Wed Jan 30, 2013 4:48 pm

All of the mechanisms used in the model are properly initialized--as long as initial conditions are unchanged, and model parameters remain constant, they generate the same results on every run.

So that leaves the other common cause of lack of reproducibility: failure to reinitialize pseudorandom number generators that are used in the course of a simulation. This model uses 150 NetStims to represent various kinds of afferent activity, and each NetStim's random is 1. However, the pseudorandom generator that supplies their random values is not reseeded before each run. Consequently if a series of runs is executed without exiting and restarting NEURON, each run has a different spatiotemporal pattern of synaptic activation.

How to see this for yourself: Your program uses this code to launch multiple runs:

Code: Select all

proc multi_run() {
   for mi = 1, 2  {
   print mi
   run()
   save_data(mi)
   }
}
multi_run()
Copy Net.hoc to testNet.hoc, then edit testNet.hoc by commenting out the "multi_run()" line, and inserting the following below it:

Code: Select all

///// record "afferent drive" (events generated by NetStims)
// there are 150 NetStims, each with random = 1

NSNUM = 150

objref nsstvec, nsidvec, tnc, nil
nsstvec = new Vector() // to hold all NetStim's event times
nsidvec = new Vector() // to hold the IDs of the active NetStims

for ii=0,NSNUM-1 {
  tnc = new NetCon(NetStim[ii],nil)
  tnc.record(nsstvec,nsidvec,ii)
}

objref tnc // don't need this any more
This uses the NetCon class's record() method to capture each NetStim's "spike times" to Vector nsstvec; in addition, for each spike time in nsstvec, the "index" of the NetStim that produced it will be captured to Vector nsidvec.

Now we need a way to run two simulations and plot the NetStims' spike rasters. Append this code after the "objref tnc" statement:

Code: Select all

///// run two simulations and show the events generated by NetStims 
NSSEED = 31 // I just picked this number out of a hat

proc onerun() {
//  NetStim[0].seed(NSSEED)
  mi = $1
  print mi
  run()
  save_data(mi)
}

objref g
g = new Graph(0)
g.size(0,3,0,150)
g.view(0, 0, 3, 150, 329, 25, 588.48, 472.96)

onerun(1)
print "Run 1: black marks"
nsidvec.mark(g, nsstvec, "|", 6, 1, 2) // thick black line
onerun(2)
print "Run 1: red marks"
nsidvec.mark(g, nsstvec, "|", 6, 2, 1) // thin red line
This runs two simulations and plots the times at which each of the NetStims spiked.

Use NEURON to execute testNet.hoc and you'll see that the first and second runs have very different afferent activity. That's why you're getting different voltage responses in your biophysical model cells.

Exit NEURON, uncomment the
// NetStim[0].seed(NSSEED)
line in testNet.hoc, save that file, and use NEURON to execute it again. Note that the red marks (second run's spike times) lie right on top of the black marks (first run's spike times). Compare the voltages recorded from your biophysical model cells and note that they are identical.

How did I know that NetStim had a seed() method? I figured it had to have a seed method, so I read the source code in netstim.mod, which is located in nrn-x.x/src/nrnoc/netstim.mod (download the source code for NEURON and expand it, if you don't already have it on your machine).

Why did I only have to call NetStim.seed() for just one NetStim? Because by default all NetStims use a shared pseudorandom sequence generator. Seed one and you have seeded them all.

If you want each NetStim to have its own sequence of spike times, independent of all other NetStims, you have to associate it with its own pseudorandom sequence generator. This is important for models that will be run on parallel hardware, since it ensures that noisy spike trains will be independent of the number of processors or how cells are distributed over those processors (see
Randomness in NEURON models
on the Documentation page of NEURON's WWW site
http://www.neuron.yale.edu/neuron/docs
).

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

Re: Multiple runs in NEURON result in different results

Post by ted » Wed Jan 30, 2013 4:51 pm

After you have had the opportunity to read this thread, I will move it to the Modeling networks area of the Forum, where it will probably do the most good.

Skyfire79

Re: Multiple runs in NEURON result in different results

Post by Skyfire79 » Fri Feb 01, 2013 11:26 am

Thank you very much for figuring this out, Ted! I really appreicate it.

I have tried the testNet.hoc and it works exactly as you described. A few questions I have are:

1) I did not know we can directly use NetStim[ii], instead of the objref associated with the NetStim object.
2) If I reseed all the NetStim objects instead of the first one at the beginning of each run, the results would be the same, right?
3) You mentioned that if I want each NetStim to have its own sequence of spike times, independent of all other NetStims, I have to associate it with its own pseudorandom sequence generator. But right now, even if all the NetStim objects have the same pseudorandom sequence generator, they still generate different spike trains. Does it mean they are independent of each other, or not?

Thanks again!
Guoshi

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

Re: Multiple runs in NEURON result in different results

Post by ted » Fri Feb 01, 2013 12:44 pm

Very good questions.
Skyfire79 wrote:I did not know we can directly use NetStim[ii], instead of the objref associated with the NetStim object.
That was a quick and dirty hack, done for the purpose of diagnosis. It is not a good programming practice. You'll probably find it very useful to read Object Oriented Programming in OC http://www.neuron.yale.edu/neuron/stati ... n/obj.html. Quoting from that document:
An object does have a unique name that can be used anywhere a reference to the object is used. But that feature is for the use of the graphical interface and should never be used in programming.
. . .
Object names or id's are not guaranteed to be the same between different sessions of NEURON unless the sequence of creation and destruction of objects of the same type are identical. The object name is defined as classname[index] where the index is automatically incremented everytime one is created. Indexes are not reused after objects are deleted except when there are no existing objects of that type; then the index starts over again at 0.
Skyfire79 wrote:If I reseed all the NetStim objects instead of the first one at the beginning of each run, the results would be the same, right?
Yes, if you used the same seed value for each NetStim. If you used different seeds, the last one you used would govern all the NetStims.
Skyfire79 wrote:You mentioned that if I want each NetStim to have its own sequence of spike times, independent of all other NetStims, I have to associate it with its own pseudorandom sequence generator.
True.
Skyfire79 wrote:But right now, even if all the NetStim objects have the same pseudorandom sequence generator, they still generate different spike trains. Does it mean they are independent of each other, or not?
It does not mean that they are independent. Consider this simple example:

Case 1: a model with one NetStim ns0 that draws from a pseudorandom number generator (rng). Suppose the rng produces this sequence of values
r0, r1, r2, r3 . . .
At the beginning of the simulation, ns0 draws r0 from the rng. This is the time ns0 must wait until it fires its first spike.
Each time ns0 spikes, it draws a new number from the rng, and that number is the latency ns0 must wait until it fires its next spike. So r0 is the time of ns0's first spike, and the rj for j>0 are the interspike intervals of ns0's spike train.

Now we change the model by adding another NetStim ns1. If ns0 and ns1 have independent spike trains, this should not affect the times at which ns0 spikes, and the times at which ns0 spikes shouldn't affect ns1's spike times. But what really happens?

Case 2: a model with two NetStims ns0 and ns1 that draw from the same pseudorandom number generator (rng). Again suppose the rng produces the same sequence of values
r0, r1, r2, r3 . . .
r0 is the time of ns0's first spike, and r1 is the time of ns1's first spike.
If r0 < r1, the first spike in the model is produced by ns0, so r2 becomes ns0's first interspike interval.
But if r1 < r0, ns1 produces the first spike in the model, and r2 becomes ns1's first interspike interval.

So adding another NetStim causes a huge difference in ns0's spike train. Adding a third NetStim will have a huge effect on ns0's and ns1's spike trains. Similarly, eliminating any NetStim will have major effects on the spike trains produced by all the remaining NetStims. Their spike trains are not independent of each other, because they are all drawing values from the same pseudorandom sequence. It is as if the toss of dice by a gambler at one craps table affected what happened to the dice thrown by a gambler at a different table.

Skyfire79

Re: Multiple runs in NEURON result in different results

Post by Skyfire79 » Fri Feb 01, 2013 1:46 pm

Excellent response and it makes perfect sense to me!
Thank you very much again, Ted!

Best,
GL

Post Reply