I can think of two straightforward ways to do this. Both take advantage of NEURON's
standard run system. Method 1 is fast but may require lots of storage, and the result
is available only after the end of a run. Method 2 uses little storage, and the result is
available for plotting during a run (no advantage for large scale simulations, which
are never run in interactive mode) but may be much slower (potentially a big drawback--
your mileage may vary).
Both require global time steps. If you must employ local variable time steps, method 1
may still be usable, but at the cost of almost twice as much storage, plus added time on
your behalf (because of having to write and debug more complex code) and that of the
computer (for some additional postprocessing, but who cares, human time is what
counts).
Method 1.
This uses a List and its append(), count() and object() methods. It also uses Vectors
and the Vector class's record(), add(), and div() methods. You will also need a custom
proc run() if you want the postprocessing to happen automatically.
When setting up your network, each time you create a synapse of interest, also create
a Vector and make it record that synapse's current. Also append each of these Vectors
to a List. It would be convenient to do all of these things with a proc.
Code: Select all
objref tvec // you'll want to keep track of time too
tvec = new Vector()
tvec.record(&t)
objref veclist, synlist
veclist = new List()
synlist = new List()
// call in the context of a currently referenced section
// arg is the location on the section where the synapse is to be placed
proc setsyn() { localobj tempsyn, tempvec
tempsyn = new ExpSyn($1) // or Exp2Syn if you prefer
. . . specify taus and reversal potential . . .
synlist.append(tempsyn)
tempvec = new Vector()
tempvec.record(&tempsyn.i)
veclist.append(tempvec)
}
After each simulation, use add() to create a new Vector whose elements equal the sum
of the corresponding elements in the Vectors in the List. Use div() to divide the elements
of this new Vector by the number of Vectors in the List.
Code: Select all
objref avgisyn
proc postprocess() { local i
avgisyn = new Vector()
avgisyn = veclist.object(0).c
for i = 0, veclist.count()-1 avgisyn.add(veclist.object(i))
avgisyn.div(veclist.count())
}
proc run() {
stdinit()
continuerun(tstop)
postprocess()
}
At the end of a run, avgisyn will contain the sequence of average synaptic currents,
and tvec will hold the corresponding times.
To use this approach with local variable time step integration, you'll have to use
the CVode class's record() method instead of Vector record(), and be sure to record
a time vector for each cell. Postprocessing will involve interpolating current values at
regular intervals for each current, time Vector pair, adding up these interpolated values,
etc..
Method 2.
Simpler but slower. Uses just a List.
Each time you set up a synapse of interest, append it to the List. You can do this with
a proc setsyn() that just leaves out the Vector stuff.
At each time step during your solution, iterate over the objects in the List to add up all
of the synaptic currents. Divide this sum by the number of synapses. Easily done with
a custom proc advance(). You will also want a custom proc init() to force your
"average current variable" to 0 on every initialization.
Code: Select all
avgisyn = 0
proc init() {
finitialize(v_init)
avgisyn = 0
if (cvode.active()) {
cvode.re_init()
} else {
fcurrent()
}
frecord_init()
}
proc advance() { local i
fadvance()
for i = 0, synlist.count()-1 avgisyn+=synlist.object(i).i
avgisyn/=synlist.count()
}