Implement extracellular in pre-existing model

[ramanujan.srinath]

I downloaded the Mainen & Sejnowski, 1996 model from ModelDB and wish to insert the extracellular mechanism to record the extracellular potential. (For now, I have one cell but eventually I will implement a network.) After inserting extracellular into all sections, all I need is the transmembrane current to be stored into a file. Is this correct?

Code: Select all

i_net = 0

forall insert extracellular

objref rect, reci
rect = new Vector()
reci = new Vector()

rect.record(&t)
reci.record(&i_net)

finitialize(-65)

objref savdata
savdata = new File()
savdata.wopen("i_membrane.dat")

while (t < tstop) {
        i_net = 0
        fadvance()
        forall {
                i_net = i_net + i_membrane(0)
                }
        print i_net
		savdata.printf("%g %g\n", t, i_net)
}

savdata.close()

Also, i_membrane is the transmembrane current density. I need the current. So I multiply by the segment area?

(I went through the extracellular_stim_and_rec file and that seems a little extensive for my project. I just need the transmembrane current out of an existing model.)

[ted]

wish to insert the extracellular mechanism to record the extracellular potential.

Do you mean the potential at a particular location in a conductive medium? Or do you want to be able to compute the potential at any arbitrary point, or multiple points?

(For now, I have one cell but eventually I will implement a network.) After inserting extracellular into all sections, all I need is the transmembrane current to be stored into a file. Is this correct?

Partly, and only if you do each of the following:
1. figure out the location in space of the center of each compartment
2. record the time course of transmembrane current generated by each compartment (segment) of the model and save all of these data to one (or more) file(s)
3. perform offline calculations with some other software that combines the time course of each current and its location in space with a model of the conductive medium in which your network is embedded, in order to calculate the potential at one or more extracellular locations in that volume

Your code excerpt is a pretty good start, but it omits many important details. Also, if you're going to write data to an output file at each time point in the simulation, there's no need to use the Vector class's record() method, is there?

i_membrane is the transmembrane current density. I need the current. So I multiply by the segment area?

Yes, this is one of those details.

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. To understand why this is so, apply an inductive proof that starts with a single compartment model, then extend to a two compartment model, and finally to an arbitrary number of compartments. For a single compartment model cell, net membrane current is 0 (capacitive current must be equal and opposite to the sum of ionic currents).

For a two compartment model with compartments numbered 0 and 1, there are three currents to consider: net membrane currents im0 and im1 (positive for outward current), and axial current that spreads from compartment 0 to 1 called iax01, the current balance equations are
im0 + iax01 = 0
-iax01 + im1 = 0
from which im0 + im1 = 0, but if the compartment centers are at different locations in space, the model should generate a dipole field whose intensity depends on the product of im0 and the distance between centers.

For a three compartment model . . . the remainder of the derivation is "left as an exercise to the reader" (a delightful quote from an old calculus textbook by Johnson and Kiokemeister), but you get the point.

(I went through the extracellular_stim_and_rec file and that seems a little extensive for my project.

It calculates the location of the center of each compartment. You might find that part helpful. It also combines that information with a very naive conceptual model of the volume conductor in which the cell is embedded (purely resistive, homogenous, isotropic, effectively infinite in size) to calculate the potential at a particular location in space.

If all you need are the membrane currents so you can use them as driving functions for a volume conductor model that is simulated with some other software, you might find it helpful to use the following mechanism (simplified from xtra.mod), because it uses compiled code to calculate the product of compartment area with local i_membrane.

Code: Select all

: xtr.mod, based on xtra.mod

COMMENT
This mechanism is intended to be used in conjunction 
with the extracellular mechanism.  A pointer specified 
at the hoc level must be used to connect the 
extracellular mechanism's i_membrane to this mechanism's im.

xtr does two useful things:

1. Calculates net membrane current inet from the product 
of local i_membrane and segment area.

2. Allows local storage of xyz coordinates interpolated from 
the pt3d data.  Once stored, these coordinates may be read out
at any time e.g. for writing to a file.
ENDCOMMENT

NEURON {
	SUFFIX xtr
	RANGE x, y, z
    RANGE inet
	POINTER im
}

PARAMETER {
	x = 0 (1) : spatial coords
	y = 0 (1)
	z = 0 (1)
}

ASSIGNED {
	v (millivolts)
	inet (nanoamp)
	im (milliamp/cm2)
	area (micron2)
}

INITIAL {
	inet = (0.01)*im*area
: this demonstrates that area is known
: UNITSOFF
: printf("area = %f\n", area)
: UNITSON
}

: Use BREAKPOINT for NEURON 5.4 and earlier
: BREAKPOINT {
:	SOLVE f METHOD cvode_t
: }
:
: PROCEDURE f() {
:	: 1 mA * 1 megohm is 1000 volts
:	: but ex is in mV
:	inet = (0.01)*im*area
: }

: With NEURON 5.5 and later, abandon the BREAKPOINT block and PROCEDURE f(),
: and instead use AFTER SOLVE

AFTER SOLVE { : after each solution step
	inet = (0.01)*im*area
}

Code: Select all

i_net = 0

forall insert extracellular

objref rect, reci
rect = new Vector()
reci = new Vector()

rect.record(&t)
reci.record(&i_net)

finitialize(-65)

objref savdata
savdata = new File()
savdata.wopen("i_membrane.dat")

while (t < tstop) {
        i_net = 0
        fadvance()
        forall {
                i_net = i_net + i_membrane(0)
                }
        print i_net
		savdata.printf("%g %g\n", t, i_net)
}

savdata.close()

[ramanujan.srinath]

Thanks for the reply Ted. That helped me in setting up and starting off with the project and, in fact, clear a few questions I had about biophysics in general!

I am directly using the "extracellular_stim_and_rec" code now. I am trying to do the exact same thing it does except with a population.

1. I exported the cell as a class and created a population with the synapses in place. (I am quite comfortable using cell/networkBuilder. Is there a way of building several thousand neuron populations with networkBuilder?)
2. Now I have to call setpointers on each of them, right?
3. And after that, I can specify the distance between each cell and the recording electrode with setelec(x,y,z).

What are the potential threats here? How do I go about integrating the vrec vectors; I can just add them up right?

[ted]

Is there a way of building several thousand neuron populations with networkBuilder?

Only by a lot of clicking and dragging. Your best bet is to save a hoc file from the NetWork Builder, extract the cell class template from that, and use the template with a bit of your own code. Example (assuming that the cell class is called SomeCell()):

Code: Select all

NCELLS = 1000
load_file("class.hoc") // contains the template of a cell class extracted from a NetWork Builder's hoc output
objref cells
cells = new List()
for i=1,NCELLS cells.append(new SomeCell())

Now I have to call setpointers on each of them, right?

Yes--you have to ensure that each section of each cell has both extracellular and xtra, and then you have to link im_xtra and ex_xtra in each segment to the i_membrane and e_extracellular in the same segment.

after that, I can specify the distance between each cell and the recording electrode with setelec(x,y,z).

setelec() specifies the location of the stimulating electrode in space. You need to make sure that each cell has the proper location and orientation in space.

Cell location is easy to handle--every cell template exported from a CellBuilder or by the Network Builder will have a public method called position() that can be used to set the coordinates of the root of the cell (the 0 end of the root section, which is almost always soma). This will control the location of all other points in the cell--move the root, and the rest of the cell will move with it automatically.

Not quite so easy is cell orientation. If you want to control the orientation of the cell in space, you need to apply a rotation transformation to the x3d, y3d, and z3d coordinates of each section's pt3d points. I'd suggest this sequence:
1. use position() to translate the root node of the cell to 0,0,0
2. control the model cell's orientation by applying the desired rotation transformation to the x3d, y3d, and z3d coordinates of each sectinon's pt3d points
3. use the model cell's position() method to place the root node of the oriented cell at the desired location in space.

Cell orientation and location must be completed before calling the procedure that finds segment centers (proc grindaway() in interpxyz.hoc).

How do I go about integrating the vrec vectors; I can just add them up right?

Yes.

[ted]

How do I go about integrating the vrec vectors; I can just add them up right?

Yes.

Sorry, that was a careless and potentially misleading answer. For the sake of clarity (mostly for anybody else who might read this thread):

If your aim is to calculate the extracellular potential generated at a particular location (x,y,z) by the activity of multiple cells, the principle of superposition applies as long as the extracellular medium is linear. That is, the potential is the sum of all the contributions from all of the membrane currents of all cells. For any given cell in a purely resistive medium, that contribution is
SUMMA rx_j * im_j
where im_j is the total membrane current of the jth compartment in the cell and rx_j is the transfer resistance between the center of that compartment and the point at (x,y,z).

If you had multiple cells and calculated a different xrec vector for each cell, you could then just do an elementwise sum of all of the vectors to get the vector that describes the time course of v(x,y,z) produced by the activity of all cells.