Resolved: calculating/recording LFP, actually was a psolve q

[neuromau]

Note: the problem I had here turned out to be relying on advance() while trying to use pc.solve(), which does not call advance(). But I'm posting my whole thought & debugging process anyway in case this is helpful for others:

I'm just starting to wrap my head around LFPs, the extracellular mechanism, and the role the extracellular mechanism can play in calculating a local field potential (LFP).

I looked through its entry in the old reference, but I can't find a proper entry for it in the new documentation, nor can I find the actual mechanism definition in my NEURON tree (searching for 'SUFFIX extracellular' or an extracellular.mod file). I was looking for the definition to get a better idea of how i_membrane works.

I also just downloaded the sample code from the Extracellular stimulation and recording post, which looks quite helpful. But for that, I see that it is calculating an LFP for a cell that is being stimulated by an intracellular current injection or an extracellular current injection. Now I'm trying to figure out what I would do differently for a detailed model of a single real cell, that is being stimulated by 100s or 1000s of netstims connected to various synapses on the cell. Note that even though I am using parallel code, for now I am just running this model for onecell on oneprocessor and haven't yet considered how I would need to modify things if I have multiple cells and/or multiple processors.

My first attempt at using the example code in my model resulted in all 0s for vrec. Looking through the forum*, I saw Ted's comment in this post that

A more important one is that the spatial distribution of current sources and sinks must also be taken into consideration. Simply adding up all transmembrane currents at any time in the simulation will produce a net result of 0.

But I *think* I got all of the necessary functionality for considering the spatial distribution of the model cell compartments from the example into my model code, namely the calcrx and grindaway functionality, and I call "define_shape()" on my model after specifying L and diam for each compartment, so it should have all the 3d point info necessary. I also checked my cell's position, which at (1.17,0,0) seems close enough to the recording electrode's location of (50,0,0). But just in case I tried moving the electrode to (5,0,0) and the vrec was still always 0.

Then, I wondered if the i_membrane actually knows about synaptic currents. From another post, I saw that Ionic and capacitive current are included:

Use extracellular's i_membrane, which includes all ionic currents and membrane capacitive current.

But then I found this post that says that
NONSPECIFIC_CURRENTs, such as those in my model synapses, are also included:

Then membrane current density will be available to you as i_membrane, and this value will include contributions from membrane capacitance plus all all ionic species and all NONSPECIFIC_CURRENTs.

I added print statements to my xtra mechanism to confirm that xtra's er was non-zero during the simulation and it was. Since I am not using an extracellular stimulating electrode, I guess I won't be playing anything into is, and therefore the is and ex within xtra should always be 0.

Then I realized that the problem was in redefining the advance() proc to include the calculation step:

proc advance() {
fadvance()
fieldrec()
}

Because I am running the simulation using pc.psolve(), which never calls advance(). The recommended alternative was to create a mod file that performs the necessary calculation after every SOLVE. In my case, looking through the fieldrec() function, I think I can accomplish its calculations by adding the following lines to my xtra mechanism:

Code: Select all

	GLOBAL lfp : added to the NEURON block
	lfp (microvolts) : added to the ASSIGNED block
        lfp = 0 : added to the BEFORE BREAKPOINT block
	lfp = lfp + er : added to the end of the INITIAL and AFTER SOLVE blocks

And then I set up a vector to record the lfp value at each time step:

Code: Select all

	myvrec.record(&lfp_xtra)

To make sure I didn't introduce any bugs with this, I did the following checks:
1. Set up an event handler that calls the fieldrec() function every so often during psolve(), and then, for several values of t, compared the values generated from the fieldrec() function with those recorded into my vector from lfp_xtra:

Code: Select all

objref fihw
fihw = new FInitializeHandler(2, "midbal()")
proc midbal() {local wt, thisstep
	print "at time t = ", t, ", the vrec was ", fieldrec()
	cvode.event(t+10, "midbal(0)")
}

Also, since I like to use cache_efficient, and now I've introduced some pointers to my code, I ran the code with cache_efficient on and off to make sure I got the same values for my cell's voltage trace and the lfp either way:

Code: Select all

cvode.cache_efficient(0) // 0 off, 1 on

So the lfp calculation that I added to the xtra mechanism, the output from the fieldrec(), and the voltage trace were consistent with each other and with cache_efficient on or off.


* the search function of the forum doesn't return anything if I search for "LFP" though I see a lot of posts with that abbreviation. I was still able to access a bunch of relevant ones with "local field potential" anyway.

[ted]

can't find a proper entry for [extracellular] in the new documentation

At
http://www.neuron.yale.edu/neuron/stati ... index.html
note the "Quick Link" called "Index". That takes you to
http://www.neuron.yale.edu/neuron/stati ... index.html
where a click on the "E" jumps to a point where it is easy to identify
"extracellular (mechanism)"
which of course brings up
http://www.neuron.yale.edu/neuron/stati ... ml#index-5

nor can I find the actual mechanism definition in my NEURON tree (searching for 'SUFFIX extracellular' or an extracellular.mod file

True, because it's not implemented as a mod file. "Inserting" extracelluar into a section means changing that section from an ordinary cable model to a cable model that has two additional layers of parallel R-C elements. Thiat's a pretty substantial change to the system equations that describe a model.

trying to figure out what I would do differently for a detailed model of a single real cell, that is being stimulated by 100s or 1000s of netstims connected to various synapses on the cell.

Nothing. Just insert extracelluar into all of its sections and simulate away.

My first attempt at using the example code in my model resulted in all 0s for vrec.

That will happen if your model cell has only one compartment, or if the transfer resistances from all of its compartments to the recording site are identical. The potential observed at some point by a monopolar electrode is the weighted sum of all of the model's membrane currents. The weights are the transfer resistances between the recording site and the membrane associated with each of the model cell's segments.

I call "define_shape()" on my model after specifying L and diam for each compartment, so it should have all the 3d point info necessary

Except for very simple stylized models, it is probably best to specify each section's geometry by the pt3d method. That's the only way to have absolute control over the orientation of each section. In a models specified with the stylized (L, diam) method, NEURON chooses the angles and branch orientations at branch points, which is OK as long as you don't really care about the shape of the model cell. But as soon as you say "I wonder what this cell's membrane currents will do to extracellular potentials" you must also be quite concerned about the cell's shape.

Then I realized that the problem was in redefining the advance() proc to include the calculation step:

Code: Select all

    proc advance() {
    fadvance()
    fieldrec()
    }

Definitely not what you want to do in a parallel simulation. You may be interested in this thread
viewtopic.php?f=31&t=3083
and in particular Michael Hines's posts about pybs.mod and discussion of allreduce.