Simulating intermittent barrage of network activity

[pascal]

I have inherited some code in which a network of 100 neurons receives a barrage of external stimulation at randomly-selected intervals (this is meant to simulate epileptic activity). The problem is that the code currently stores the time of every single externally-supplied synaptic event in a vector before the numerical integration, which results in the program running out of memory very quickly. I need to modify the code to eliminate these memory issues. I've read the post about using vecevent (http://www.neuron.yale.edu/phpBB/viewtopic.php?f=16&t=1454&p=5088&hilit=netstim+net_receive#p5088), but I don't think this is useful in my case because my code already stores just event times (as opposed to a vector of 0's and 1's for each time step), and it still has memory issues.

So it seems that the best thing to do is generate noise events "on the fly." I'm new to NEURON, so I coded up a very simple network of just one biophysical neuron driven by one NetStim "neuron," and I tried to see if I could get the NetStim neuron to activate at two different times which I specified in the course of integration, so that my run() procedure looks like this (cells.object(1) is the NetStim cell):

Code: Select all

proc run(){
   stdinit()
   continuerun(3000)
   cells.object(1).pp.start = 3001
   continuerun(tstop)
   showraster()
}


The NetStim template looked like this:

Code: Select all

begintemplate NetStim_NetStim
public pp, connect2target, x, y, z, position, is_art
objref pp
proc init() {
  pp = new NetStim()
    pp.start = 1000
    pp.number = 100
    pp.noise = 0.5
}
func is_art() { return 1 }
obfunc connect2target() { localobj nc
  nc = new NetCon(pp, $o1)
  if (numarg() == 2) { $o2 = nc }
  return nc
}
proc position(){x=$1  y=$2  z=$3}
endtemplate NetStim_NetStim


so you can see that originally the NetStim object started stimulating the biophysical neuron at t=1000 ms.

Here is the raster plot for the biophysical neuron:



Obviously, the NetStim did not turn back on at t=3001 ms like I wanted it to. It seems like there should be a way to do this, but I've read the NetStim documentation and still don't see how. I saw that NetStim can also have a NET_RECEIVE block, which would allow other neurons in the network to initiate further NetStim activity, but this is not helpful to me because I want to directly control when the NetStim object is active. I'd also like to be able to modify how long it is active within a given burst. Is this possible? And if not, what is the best approach to administer bursts of noisy stimulation "on the fly"? Thanks.

[ted]

I have inherited some code in which a network of 100 neurons receives a barrage of external stimulation at randomly-selected intervals (this is meant to simulate epileptic activity). The problem is that the code currently stores the time of every single externally-supplied synaptic event in a vector before the numerical integration, which results in the program running out of memory very quickly.

There must be one heck of a lot of "externally supplied events" if the event times can't be read into a Vector without causing memory overflow.

I've read the post about using vecevent

Perhaps that thread is too confusing. 0s and 1s have nothing to do with the VecStim class that is defined by vecevent.mod--see the vecevent* files in nrn/examples/nrniv/netcon. But that's irrelevant, because if there are so many "externally supplied events" that they can't be stored in Vectors, VecStim isn't going to help.

Code: Select all

So it seems that the best thing to do is generate noise events "on the fly." I'm new to NEURON, so I coded up a very simple network of just one biophysical neuron driven by one NetStim "neuron," and I tried to see if I could get the NetStim neuron to activate at two different times which I specified in the course of integration

There is no need to write a custom run() procedure to accomplish this. The documentation of the NetStim class http://www.neuron.yale.edu/neuron/static/docs/help/neuron/neuron/mech.html#NetStim tells you what to do: drive it with events.

"Where will those events come from?"

From anything that can serve as an event source--a NetStim, a VecStim, some other artificial spiking cell, or a biophysical model spiking cell, which is connected to the target NetStim by a NetCon.

[pascal]

Thanks, Ted, that is really helpful. Following your advice, I was able to induce two separate barrages of input from NetStim using the following code:

Code: Select all

proc run(){ local i
   stdinit()
   //solve(tstop)
   //added the following 2 lines
   continuerun(3000)
   for(i=3001;i<4000;i=i+5){
      nclist.object(0).event(i)
   }
   continuerun(5000)
   showraster()
}


I am curious to read about how NetCon events work. I looked through The NEURON Book (including pp. 288-290) and the Programmer's Reference, but they didn't answer all my questions. For example, is it possible to delete events from the event queue? Also, are there separate queues for NetCon events and for CVODE events? (I have seen some hoc code that calls cvode.event and netcon.event in succession.) Thanks for the help.

[ted]

Interesting. It works, but stuffing a bunch of events into the event queue is not generally a good idea (for one thing, it increases the overhead of maintaining the queue). It is better to do this: after the network setup code has executed, but before calling run(), fill a Vector with the event times e.g. for your particular example

Code: Select all

objref stimes
stimes = new Vector(1000/5)
stimes.indgen.mul(5).add(3001)

(this does the following things, in sequence, at machine language speeds:
fills the Vector with 0,1,2...
multiplies each element by 5
adds 3001 to each element).
Then use an instance of the VecStim class to deliver the events in the course of simulation. You could then run your simulation and display results by calling myrun()

Code: Select all

tstop = 5000
proc myrun() {
  run()
  showraster()
}
myrun()

This has the added advantage of leaving the standard run system intact, which can be important if you ever have the need to further develop or debug your code.

You'll find the NMODL source code for VecStim, plus a couple of other files that demonstrate its usage, in nrn/share/nrn/examples/nrniv/netcon -- look for the vecevent* files.

Good questions. A lot of very useful things can be done with the event delivery system, including sophisticated control of program flow, e.g. for implementing "virtual experimental protocols" in which model parameters are adjusted on the fly, or complex stimuli are generated dynamically.

The event delivery system had just been added to NEURON when The NEURON Book was written. You'll want to read the Programmer's Reference entries on the CVode class's event() method, and the NetCon class's event() method. You'll also want to read about the FInitializeHandler class.

is it possible to delete events from the event queue?

No. What did you have in mind that would require that?

are there separate queues for NetCon events and for CVODE events?

No. The event delivery system has only one queue.

I have seen some hoc code that calls cvode.event and netcon.event in succession.

Really? Where? I wonder what the programmers were trying to achieve.

[pascal]

Sorry I'm slow with this...I'm confused because I thought you advised me to drive the NetStim with events. Also, I don't understand why it's not a good idea to stuff a bunch of events into the event queue--isn't that what NetStim typically does before a simulation starts? I thought netstim.start, netstim.interval, and netstim.number helped to define (with the help of a random number generator) a bunch of event times that are placed in the event queue before the simulation begins. Why does it matter whether events are added before or in the middle of a simulation?

The reason I want to add events in the middle of the simulation is because the code I've inherited is running out of memory when a TON of events are defined before the simulation starts. (To be clear, I don't think the code is using NetStim. That's just something I tried in the toy model above). So currently, the code looks like this:

[define a ton of events]
psolve(Tstop)

So I thought, why not break up the definition of the events into smaller pieces? Then we could do something like this:

[define a number of events to take us through the first 1/10 of simulation]

for i=1,10{
psolve(i/10*Tstop)
[delete the previous events in the queue, and add events for the next tenth of the simulation]
}

This is why I am interested in whether you can delete events in the event queue, so that they don't add up and consume all available memory. I think events are deleted after they are realized, but I'm not sure.

[ted]

The discussion drifting away from the problem that motivated your first post: how to implement repetitive bursts of afferent events in a way that avoids storing event times in a Vector. Let me return to that before I try to address your most recent questions.

The simplest implementation would take advantage of this property of the NetStim class (quoting from the Programmer's Reference entry):

If the stimulator is in the on=0 state and receives a positive weight event, then the stimulator changes to the on=1 state and goes through its burst sequence before changing to the on=0 state. During that time it ignores any positive weight events. If, in the on=1 state, the stimulator receives a negative weight event, the stimulator will change to the off state. In the off state, it will ignore negative weight events. A change to the on state immediately causes the first spike.

This could be done with pair of NetStims connected like so

Code: Select all

NS1         NS2
   o-------------< o-------------<
      NC1

NS2 produces the events that you deliver to whatever targets you want.
NS2.start = 1e9 so it never turns on by itself.
NS2.interval, number, and noise specify the statistics of the spikes in any one of its bursts.
NS1's parameters govern the times at which NS2's bursts will start.
NC1.weight = 1, delay = 0.

If you prefer to use precomputed event times that are read from a text file, you could replace NC1 with a VecStim.

Sorry I'm slow with this...I'm confused because I thought you advised me to drive the NetStim with events.

Now you know what I meant to say.

Also, I don't understand why it's not a good idea to stuff a bunch of events into the event queue

It doesn't scale well.

isn't that what NetStim typically does before a simulation starts?

No. When a NetStim generates an output event, it also determines if it should produce another output event at some future time. If the answer is "yes", it sends itself a self-event that returns at that time. Arrival of that self-event triggers production of an output event and another determination of if & when to produce the next output event. So a NetStim never puts more than two events on the queue--an output event, and a self-event--regardless of the values of its "number," "interval," or "noise" parameters.

[pascal]

Thanks so much, Ted. This toy model is officially working. I used a VecStim artificial cell to stimulate a NetStim artificial cell, which in turn stimulates a biophysical neuron. Here is the code for the network setup:

Code: Select all

begintemplate Cell_Cell
public is_art
public init, topol, basic_shape, subsets, geom, biophys, geom_nseg, biophys_inhomo
public synlist, x, y, z, position, connect2target

public soma
public all

objref synlist

proc init() {
  topol()
  subsets()
  geom()
  biophys()
  geom_nseg()
  synlist = new List()
  synapses()
  x = y = z = 0 // only change via position
}

create soma

proc topol() { local i
  basic_shape()
}
proc basic_shape() {
  soma {pt3dclear() pt3dadd(0, 0, 0, 1) pt3dadd(15, 0, 0, 1)}
}

objref all
proc subsets() { local i
  objref all
  all = new SectionList()
    soma all.append()

}
proc geom() {
}
external lambda_f
proc geom_nseg() {
}
proc biophys() {
  soma {
    insert hh
      gnabar_hh = 0.12
      gkbar_hh = 0.036
      gl_hh = 0.0003
      el_hh = -54.3
  }
}
proc biophys_inhomo(){}
proc position() { local i
  soma for i = 0, n3d()-1 {
    pt3dchange(i, $1-x+x3d(i), $2-y+y3d(i), $3-z+z3d(i), diam3d(i))
  }
  x = $1  y = $2  z = $3
}
obfunc connect2target() { localobj nc //$o1 target point process, optional $o2 returned NetCon
  soma nc = new NetCon(&v(1), $o1)
  nc.threshold = 10
  if (numarg() == 2) { $o2 = nc } // for backward compatibility
  return nc
}
objref syn_
proc synapses() { //this is the result of my adding one synapse in the GUI
  /* ExpSyn0 */   soma syn_ = new ExpSyn(0.6)  synlist.append(syn_)
}
func is_art() { return 0 }

endtemplate Cell_Cell



begintemplate NetStim_NetStim
public pp, connect2target, x, y, z, position, is_art
objref pp
proc init() {
  pp = new NetStim()
    pp.start = -1
    pp.number = 25
    pp.noise = 0.5
}
func is_art() { return 1 }
obfunc connect2target() { localobj nc
  nc = new NetCon(pp, $o1)
  if (numarg() == 2) { $o2 = nc }
  return nc
}
proc position(){x=$1  y=$2  z=$3}
endtemplate NetStim_NetStim

//VecStim template
begintemplate vs_VecStim
public pp, connect2target, x, y, z, position, is_art
objref pp
proc init() {
  pp = new VecStim()
}
func is_art() { return 1 }
obfunc connect2target() { localobj nc
  nc = new NetCon(pp, $o1)
  if (numarg() == 2) { $o2 = nc }
  return nc
}
proc position(){x=$1  y=$2  z=$3}
endtemplate vs_VecStim

//Network specification interface

objref cells, nclist, netcon
{cells = new List()  nclist = new List()}

func cell_append() {cells.append($o1)  $o1.position($2,$3,$4)
   return cells.count - 1
}

//this function creates a connection between two cells, and adds it to the list of netcons
func nc_append() {//srcindex, tarcelindex, synindex
  if ($3 >= 0) {
    //cells.object($1) picks out the $1 element of the 'cells' list; then you connect this cell to the target specified in the argument to 'connect2target,' in this case the $2 element of the 'cells' list, this cell might have multiple synapses, so you connect it to the $3 synapse on that cell
    netcon = cells.object($1).connect2target(cells.object($2).synlist.object($3))
    //this just specifies the weight and delay associated with this particular netcon
    netcon.weight = $4   netcon.delay = $5
  }else{
    netcon = cells.object($1).connect2target(cells.object($2).pp)
    netcon.weight = $4   netcon.delay = $5
  }
  nclist.append(netcon)
  return nclist.count - 1
}

//Network instantiation

  //so this just adds the first (and only) HH cell to the network, and specifies its location
  /* Cell0 */  cell_append(new Cell_Cell(),   -35,    -48, 0)
 
//this adds a netstim cell to the network
  /* NetStim1 */  cell_append(new NetStim_NetStim(),   5,    -8, 0)
  //now add a vecstim
  cell_append(new vs_VecStim(),0,0,0)

  //this adds a connection from Cell #1 (netstim) to Cell #0 (HH cell), to the 0th synapse on Cell #0, with a weight of 0 and a delay of 1 (ms?)
  /* NetStim1 -> Cell0.ExpSyn0 */  nc_append(1,   0, 0,  0.5,1)

  //add connection from Cell #2 (VecStim) to Cell #1 (netstim)
  //(important to make third entry '-1' (see pg. 328 NEURON Book)
   nc_append(2,1,-1,1,0)


I got the template for the VecStim by loading "vecevent.ses" in nrn\examples\nrniv\netcon and then having the network builder generate code. And here is the code for the simulation and plotting:

Code: Select all

v_init = -65
tstop = 5000
dt=0.025

objref rect,recv
rect = new Vector()
recv = new Vector()

//make *sure* to include '.object' when it's a list; otherwise NEURON will just crash
recv.record(&cells.object[0].soma.v(0.5))
rect.record(&t)

//prepare to record and display spike trains (from p. 338 in The NEURON Book)
objref netcon,vec,spikes,nil,graster

proc preprasterplot(){
   spikes = new List()
   //for i=0,cells.count()-1{
      vec = new Vector()
      netcon = new NetCon(&cells.object(0).soma.v,nil)
      netcon.record(vec)
      spikes.append(vec)
   //}
   objref netcon, vec
 
   graster = new Graph(0)
   graster.view(0,0,tstop,cells.count(), 300,105,300.48,200.32)
}

preprasterplot()
objref spikey

proc showraster(){
   graster.erase_all()
   //for i = 0,cells.count()-1 {
      spikey = spikes.object(0).c
      spikey.fill(1)
      spikey.mark(graster,spikes.object(0), "|", 6)
   //}
   objref spikey
}

objref vec
vec = new Vector(4)
//define the times at which you want to stimulate a barrage of activity
vec.x[0]=1000
vec.x[1]=2000
vec.x[2]=3000
vec.x[3]=4000
cells.object(2).pp.play(vec)

proc go(){
   stdinit()
   run()
   showraster()
}


The output of the biophysical neuron of the biophysical neuron shows that the simulation is working as intended, inducing bursts of activity at 1000, 2000, 3000, and 4000 ms:



One thing to note is that when I set netstim.start=1e9, I could not elicit any spiking. Apparently netstim will not respond to any input from vecstim until after the time specified by netstim.start. Therefore setting netstim.start=-1 solved the problem. Thanks for all the help, Ted. I think VecStim should be able to solve the memory problems in the code I mentioned at the start of this post...

[ted]

One thing to note is that when I set netstim.start=1e9, I could not elicit any spiking.

My mistake. I was trying to prevent the NetStim from producing events on its own, but making start == 1e9 prevented from doing anything until after t == 1 billion ms. It's good that setting start < 0 did the job.

Looking at the NMODL source code, I see why this worked. start is used in the INITIAL block, where this conditional clause determines whether or not to launch the self-event that triggers the first "spike" when the NetStim is "allowed" to start spiking on its own:

Code: Select all

        if (start >= 0 && number > 0) {
                on = 1
                : randomize the first spike so on average it occurs at
                : start + noise*interval
                event = start + invl(interval) - interval*(1. - noise)
                : but not earlier than 0
                if (event < 0) {
                        event = 0
                }
                net_send(event, 3)
        }

If start<0, that first self-event is not generated, so the NetStim stays stumm until an input event with weight>0 kicks it into action.

[rth]

Hi Ted,
Maybe a little bit off-top.
I use NetStim->NetCon->Exp2Syn system to get a barrage of synaptic inputs on distal dendrites of 8 compartments pyramidal cell model. (The code is too long and too simple to post it here) It worked very nice and did what I want until I tried to run simulation with ParallelContext using threads. NEURON returns an error:
/usr/local/nrn/i686/bin/nrniv: multithread random in NetStim only via hoc Random
So I'm looking for the way to fix it. I can generate Poisson's process by random generator, but how create events for synapses isn't clear to me. Maybe you can give some simple example?
Thanks,
Ruben

[ted]

This post is for the benefit of future readers of this thread.

At first I didn't understand the question, partly because of the cryptic nature of the error message. However, my confusion was dispelled when rth showed me what was happening. Each NetStim was producing an event stream whose timings were generated by a shared pseudorandom number generator. However, this wouldn't work during multithreaded execution; instead, each NetStim would have to be paired with its own instance of the Random class, and each of those would have to be seeded in such a way as to ensure that its stream doesn't overlap with that produced by any other stream. This can be done with the same basic strategy that is used when randomness is needed in a distributed network model--see
Randomness in NEURON models
on the Documentation page http://www.neuron.yale.edu/neuron/docs