Calculating and plotting variables

[AnotherNewUser]

Hi,

A completely new user and new to programming in general.

I quite happily implemented a simple network described in a paper and, to my surprise, it reproduced the results.

I'd like to determine the average postsynaptic current in one polulation of the neurons and have become completely vexed in the process. I would have thought something like

meani=0
for k=0,ncells-1{
for j=0,ncells-1{
meani = meani+(pypyG[k][j].i+inpyA[k][j].i+inpyB[k][j].i+tcpyG[k][j].i)/ncells
}
}

was the way to start, but apparently not. Further, I don't see how to plot a quantity like this using the GUI, although this is probably straightforward too. Any guidence would be appreciated.

[ted]

I quite happily implemented a simple network described in a paper and, to my surprise, it reproduced the results.

A bigger surprise is that the paper contained enough detail so that it was reproducible.
But the truly astonishing fact is that the journal printed those details, free of typos.

I'd like to determine the average postsynaptic current in one polulation of the neurons

Please be more specific. Do you mean the instantaneous average, or the average over
the course of the entire simulation? Does "postsynaptic current" mean the current that
enters through synapses, or all (? ionic) transmembrane currents in the postsynaptic
cells?

[AnotherNewUser]

Hi,

Thanks for the very rapid reply. My hope is to keep track of the instantaneous average current through the synapses of one cell type in the network, not the totality of their transmembrane currents and not the average current over the course of the simulation.

Thanks Again

[ted]

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()
}

[AnotherNewUser]

Thank you. It's very helpful to see two ways of working out the problem.