recording currents

Anything that doesn't fit elsewhere.
maria.kesa
Posts: 36
Joined: Thu Nov 12, 2015 10:45 am

Re: recording currents

Post by maria.kesa » Thu Jan 21, 2016 2:56 pm

I set the recording resolution to 1 ms, maybe this will help. The simulation is running right now.

I did as you said. Summed the vectors together at the end of the program.

maria.kesa
Posts: 36
Joined: Thu Nov 12, 2015 10:45 am

Re: recording currents

Post by maria.kesa » Fri Jan 22, 2016 8:16 am

Thanks for your help, Ted. Recording at 1 ms precision makes the memory manageable and the NSG cluster runs the simulation.

And you were right, only the pyramidal cell inputs should be counted in the approximation to the LFP.

maria.kesa
Posts: 36
Joined: Thu Nov 12, 2015 10:45 am

Re: recording currents

Post by maria.kesa » Fri Jan 22, 2016 10:10 am

Dear Ted,

I still have memory issues at NSG even after setting the timestep for recording at 1ms.
Right now i tried removing the currents to inhibitory neurons-- you were right, they don't contribute so much to the lfp. But I think that memory will still be an issue?

In terms of the memory requirements, I think the only option left is to do recording at every times step. Could you use fadvance() somehow for that? I tried, but failed. I would appreciate your help.

Thanks,
Maria

maria.kesa
Posts: 36
Joined: Thu Nov 12, 2015 10:45 am

Re: recording currents

Post by maria.kesa » Sun Jan 24, 2016 1:38 pm

I tried this, but at the end sumlfp.count()
0

Code: Select all


proc record_timestep()  { local i  localobj tmpobj, intermsumlfp, tmp
  intermsumlfp= new Vector()
  for i=0,$o1.count()-1 {
    tmpobj = $o1.o(i)
    if (tmpobj.isampa == 1) {
      tmp=abs(tmpobj.iampa)+abs(tmpobj.inmda)
      intermsumlfp.append(tmp)
    } else {
      tmp=abs(tmpobj.igaba)
      intermsumlfp.append(tmp)
    }
  } 
}

proc advance() {
  fadvance()
  record_timestep($o1)
  sumlfp.append(intermsumlfp.sum())
}
advance(synapselist)
Could you help me to modify the code so that it would work?

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

Re: recording currents

Post by ted » Sun Jan 24, 2016 9:49 pm

maria.kesa wrote:I think the only option left is to do recording at every times step.
You don't want to do that--performance would be bad when executed serially, and terrible if executed in parallel. A far better strategy is to break the simulation into intervals just big enough that the recorded data will fit into memory with a bit of room to spare, and at the end of each interval add up the absolute values of recorded currents and append that result to a single Vector that finally contains the complete time course of the summed absolute currents. I have just now implemented and tested that, but will need to clean up and comment the code a bit before giving it to you tomorrow.

maria.kesa
Posts: 36
Joined: Thu Nov 12, 2015 10:45 am

Re: recording currents

Post by maria.kesa » Mon Jan 25, 2016 10:16 am

Thank you, Ted! I'll be waiting.

Maria

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

Re: recording currents

Post by ted » Tue Jan 26, 2016 12:48 am

The basic idea is to break a long simulation into a sequence of intervals or segments short enough that the segment's recorded data, and the cumulative record of time and sum of absolute synaptic currents, will fit into memory. Here is a minimal generic proc segrun() with comments that indicate where application-specific statements would be inserted. Note that the localobj declarations in the opening statement of proc segrun are all application-specific.

Code: Select all

proc errexit() {
  print $s1
  stop
}

// executes a long simulation as a series of shorter segments
// args 1-3 are inputs to segrun, 4 and 5 are outputs from segrun
// $1 total duration of simulation
// $2 number of time segments
// application-specific comments

// note:  the localobj declarations are all application-specific
proc segrun() { local i, interval, count, t2, lastpass // localobj ...
  ///// bulletproofing /////
  if ($1<=0) errexit("stop time must be > 0")
  if ($2<1) errexit("number of time intervals must be >= 1")

  ///// start of task-specific statements that must be executed before launching simulation
  // . . .
  ///// end of task-specific statements that must be executed before launching simulation

  interval = $1/$2 // segment length in ms (final segment may be shorter)
  count = 0 // index of pass through the main computational loop
  lastpass = 0 // == 1 if current simulation segment is last pass through loop
  stdinit() // assumes v_init has been specified elsewhere

  ///// main computational loop
  while (t<$1) {
    t2 = t+interval
    if (t2>=$1) {
      t2=$1
      lastpass = 1 // this is the last pass through the loop
    }

    continuerun(t2)

    ///// start of task-specific statements done after completion of a simulation segment
    // . . .
    ///// end of task-specific statements done after completion of a simulation segment

    count+=1
  }

  ///// start of any other task-specific statements that must be done before exiting proc segrun()
  // . . .
  ///// end of any other task-specific statements that must be done before exiting proc segrun()
}
The advantage of this approach is that all operations on data are carried out by compiled vector operations that are much faster than iterating over all synaptic instances to add up the absolute values of their currents at each time step.

Here are the application-specific statements for your particular problem. This assumes that you already have a List called synlist whose elements point to the synaptic mechanism instances.

First, proc segrun's initial and "declaration statement" become

Code: Select all

// executes a long simulation as a series of shorter segments
// args 1-3 are inputs to segrun, 4 and 5 are outputs from segrun
// $1 total duration of simulation
// $2 number of time segments
// $o3 List of objrefs that point to all synaptic mechanisms whose currents are to be recorded
// $o4 Vector that will contain entire record of time
// $o5 Vector that will contain entire record of sum of absolute synaptic currents

proc segrun() { local i, interval, count, t2, lastpass \
                localobj segtvec, tvec, tmpvec, isynvecs, segsumabsisyn, sumabsisyn
Next, the task-specific statements that must be executed before launching a simulation are

Code: Select all

  segtvec = new Vector() // records t in current segment
  segtvec.record(&t)
  tvec = new Vector() // will contain entire record of time

  tmpvec = new Vector()
  isynvecs = new List() // elements will be Vectors that record individual synaptic currents
  for i=0,$o3.count()-1 {
    tmpvec.record(&$o3.o(i).i)
    isynvecs.append(tmpvec)
    tmpvec = new Vector()
  }

  segsumabsisyn = new Vector() // sum of absolute synaptic currents for current segment
  sumabsisyn = new Vector() // will contain entire record of sum of absolute synaptic currents
The task-specific statements done after completion of a simulation segment are

Code: Select all

    tvec.append(segtvec) // save results from this interval and prepare for next interval

    segsumabsisyn = sumabsvecs(isynvecs)
    sumabsisyn.append(segsumabsisyn)

    printf("interval # %d ended at %f\n", count, t)
    printf("%d data points in this interval\n", segtvec.size())
    printf("%d data points captured so far\n", tvec.size())

    if (lastpass==0) { // prevent duplication at end of internal segments
      tvec.remove(tvec.size()-1)
      sumabsisyn.remove(sumabsisyn.size()-1)
    }
    frecord_init() // reset Vector recording counters
And the "other task-specific statements" done before exiting proc segrun() are

Code: Select all

  $o4 = tvec.c
  $o5 = sumabsisyn.c
Of course this depends on the following obfunc

Code: Select all

// $o1 List of Vectors
// returns a Vector whose elements are the sum of the abs values
// of the corresponding elements in $o1's Vectors

obfunc sumabsvecs() { local i  localobj tmpvec, tmpsum
  tmpsum = new Vector($o1.o(0).size) // same size as vecs in $o1, initialized to 0s
  for i=0,$o1.count()-1 {
    tmpvec = $o1.o(i).c
    tmpvec.abs()
    tmpsum.add(tmpvec)
  }
  return tmpsum
}
which can be conveniently inserted right after the declaration of proc errexit().

https://www.neuron.yale.edu/ftp/ted/neu ... xample.zip
contains a complete working example. Use NEURON to execute init_segrun.hoc. All of the interesting stuff starts below the comment

Code: Select all

///// "model specification" is now complete

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

Re: recording currents

Post by ted » Tue Jan 26, 2016 2:25 pm

A couple more comments are appropriate.

First, the example that I uploaded involves synaptic mechanisms that generate currents called i. The synaptic mechanisms you're using don't have a RANGE variable called i, but you probably already know how to deal with that. For the benefit of others who may read this thread, the steps are
1. Insert
RANGE i
as the last statement of the NEURON block.
2. Insert the declaration
i (nA)
as the last statement of the ASSIGNED block.
3. If the mechanism is GABAergic, insert
i = igaba
in the BREAKPOINT block right after the statement that calculates the value of igaba--the latter is probably something that looks like

Code: Select all

igaba = W*g_gaba*(v - Erev_gaba)
If instead the mechanism is AMPAergic, find the statements that calculate the values of inmda and iampa, and after both have been calculated, insert this statement
i = iampa + inmda

Second, the code I provided will work only if you execute the simulations serially. It is possible to revise the code to allow parallel simulation, but that involves other complications (transferring Vectors from the various hosts to host 0 at the end of each segment, then adding up the currents on host 0 and appending the results etc.), but that would take more time to develop. Furthermore, I think the effort would be misplaced--it would be better to try to revise the synaptic mechanisms so that each instance could handle converging inputs from multiple presynaptic spike sources. That would significantly reduce the number of Vectors that must be used to record synaptic currents, and decrease memory requirements substantially.

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

Re: recording currents

Post by ted » Tue Jan 26, 2016 3:58 pm

One more caveat about this particular model: I just noticed that at least two of the synaptic mechanisms calculate the synaptic currents in the DERIVATIVE block, not in the BREAKPOINT block. I suspect that all of the point processes that generate currents do the same thing. This introduces an error into simulations, because it means that the dynamics of the model are completely unaffected by the conductance changes caused by synaptic activation. You see, at each fadvance the effects of membrane conductance changes on model dynamics are taken into account by executing statements in the BREAKPOINT block twice--once with an argument of v+0.001 and again with an argument of v--in order to calculate a numerical estimate of the derivative di/dv. This conductance is added to the corresponding diagonal element of the Jacobian matrix that had been set up during initialization. This means that if a density mechanism or point process generates a current, the statement in which that current is calculated should be executed in the BREAKPOINT block.

For example, look at pyrD2interD_STFD.mod and you'll see that the BREAKPOINT block is simply

Code: Select all

BREAKPOINT {
  SOLVE release METHOD cnexp
}
and that the synaptic currents are computed in the DERIVATIVE block

Code: Select all

DERIVATIVE release {
  . . .
  inmda = W_nmda*g_nmda*(v - Erev_nmda)*sfunc(v)
  . . .
  iampa = W*g_ampa*(v - Erev_ampa)
  . . .
}
This mechanism can be fixed by simply copying them to the BREAKPOINT block after the SOLVE statement, like so

Code: Select all

BREAKPOINT {
  SOLVE release METHOD cnexp
  inmda = W_nmda*g_nmda*(v - Erev_nmda)*sfunc(v)
  iampa = W*g_ampa*(v - Erev_ampa)
}
Of course, to prevent confusion among those who may subsequently look at this file and wonder "why are these two variables calculated in two places?", it would be a very good idea to comment them out in the DERIVATIVE block (that is, insert a colon character i.e. ":" as the first character in those two lines in the DERIVATIVE block) and maybe also a short comment that says why, e.g.

Code: Select all

: must calculate inmda and iampa in the BREAKPOINT block
And, of course, in the context of my previous email, the final revised BREAKPOINT block would look like this:

Code: Select all

BREAKPOINT {
  SOLVE release METHOD cnexp
  inmda = W_nmda*g_nmda*(v - Erev_nmda)*sfunc(v)
  iampa = W*g_ampa*(v - Erev_ampa)
  i = inmda + iampa : a convenience for calculating the LFP
}
Now we consider a potentially difficult issue. Someone who reads this might well ask "Wait a minute--if these synaptic mechanisms have to be revised to eliminate an error, doesn't this invalidate the entire model and the conclusions that the original authors drew from their simulations?"

Not necessarily. First realize that all numerical integration methods produce approximate results, so there is a sense in which all numerical simulations are "incorrect." Second, it remains to be determined whether "fixing" these synaptic mechanisms will have any significant effect on simulation results. Yes, in a toy model (single compartment with hh, with a single instance of pyrD2interD_STFD attached) there is a slight difference in the time course of membrane potential during a spike elicited by a single input event. That does shift the time of the spike by a time step or two. But possibly more important would be what happens to the dynamics of the model cells when membrane potential is fluctuating in the subthreshold region; this is the critical zone, because it's where cells spend most of their time, and what they're doing is "thinking" about whether or not to spike. Anything that affects membrane conductance while a cell is in the subthreshold region could possibly have a big effect on network operation.

So here's what needs to be done. Someone (the model's original authors, preferably) has to make sure that all of the synaptic mechanisms (and all of the density mechanisms, too) compute their currents correctly--with statements in the BREAKPOINT block. And that someone then needs to repeat the simulations that served as the basis for the results reported by the original model's authors in order to verify whether and how the reported findings are affected. If the reported findings are qualitative observations, there is a good chance that the paper's conclusions remain valid. Semiquantitative findings are also likely to be robust--that is, essentially unaffected. Strongly quantitative findings are more likely to need to be revised.

Post Reply