recording currents
-
- Posts: 36
- Joined: Thu Nov 12, 2015 10:45 am
Re: recording currents
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.
I did as you said. Summed the vectors together at the end of the program.
-
- Posts: 36
- Joined: Thu Nov 12, 2015 10:45 am
Re: recording currents
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.
And you were right, only the pyramidal cell inputs should be counted in the approximation to the LFP.
-
- Posts: 36
- Joined: Thu Nov 12, 2015 10:45 am
Re: recording currents
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
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
-
- Posts: 36
- Joined: Thu Nov 12, 2015 10:45 am
Re: recording currents
I tried this, but at the end sumlfp.count()
0
Could you help me to modify the code so that it would work?
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)
-
- Site Admin
- Posts: 6299
- Joined: Wed May 18, 2005 4:50 pm
- Location: Yale University School of Medicine
- Contact:
Re: recording currents
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 wrote:I think the only option left is to do recording at every times step.
-
- Posts: 36
- Joined: Thu Nov 12, 2015 10:45 am
Re: recording currents
Thank you, Ted! I'll be waiting.
Maria
Maria
-
- Site Admin
- Posts: 6299
- Joined: Wed May 18, 2005 4:50 pm
- Location: Yale University School of Medicine
- Contact:
Re: recording currents
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.
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
Next, the task-specific statements that must be executed before launching a simulation are
The task-specific statements done after completion of a simulation segment are
And the "other task-specific statements" done before exiting proc segrun() are
Of course this depends on the following obfuncwhich 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
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()
}
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
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
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
Code: Select all
$o4 = tvec.c
$o5 = sumabsisyn.c
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
}
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
-
- Site Admin
- Posts: 6299
- Joined: Wed May 18, 2005 4:50 pm
- Location: Yale University School of Medicine
- Contact:
Re: recording currents
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
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.
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)
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.
-
- Site Admin
- Posts: 6299
- Joined: Wed May 18, 2005 4:50 pm
- Location: Yale University School of Medicine
- Contact:
Re: recording currents
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 simplyand that the synaptic currents are computed in the DERIVATIVE blockThis mechanism can be fixed by simply copying them to the BREAKPOINT block after the SOLVE statement, like soOf 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.
And, of course, in the context of my previous email, the final revised BREAKPOINT block would look like this:
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.
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
}
Code: Select all
DERIVATIVE release {
. . .
inmda = W_nmda*g_nmda*(v - Erev_nmda)*sfunc(v)
. . .
iampa = W*g_ampa*(v - Erev_ampa)
. . .
}
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)
}
Code: Select all
: must calculate inmda and iampa in the BREAKPOINT block
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
}
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.