The math behind extracellular mechanisms
The math behind extracellular mechanisms
Hello, we're trying to port a NEURON model over to another simulator system in order to take advantage of some nonstandard hardware and accelerate our simulations. The model we're working with makes use of the extracellular mechanism, and I'm running into a problem getting the results to line up between NEURON and the other simulator (the simulators line up when the extracellular mechanisms are not included). I think that the problem comes down to understanding what the underlying math used by the extracellular mechanism is, but I'm having a great deal of trouble understanding that from the extcell.c source.
The model is using the extracellular mechanism to model the effects of a myelin sheath on an axon, so there is only one extracellular layer for which xraxial, xc, xg, and vext are defined. From the mechanism documentation at http://www.neuron.yale.edu/neuron/stati ... racellular, my understanding is that there should a differential equation associated with vext, something along the lines of d_vext/dt = (1/xc)(I_membrane + I_myelin + I_periaxonal) where I_membrane is the current flowing through the cell membrane, I_myelin is the current flowing through the myelin sheath, which should be equal to xg*vext and I_periaxonal is the current flowing along the periaxonal space, which should be equal to xraxial*(vext  vext_neighbor) where vext_neighbor is the neighboring section.
My questions are:
1. Are the equations I've outlined above correct, or is there something else going on in the extracellular mechanism that I'm missing.
2. For converting the "densities" of xraxial, xc, and xg to absolute values (i.e. Mhos/cm^2 to Mhos) does the extracellular mechanism use the standard length and diameter defined for the section? Because we're modeling a myelin sheath the diameter of the sheath is different from the diameter of the section, but I don't see anywhere that this would be defined in the extracellular mechanism.
The model is using the extracellular mechanism to model the effects of a myelin sheath on an axon, so there is only one extracellular layer for which xraxial, xc, xg, and vext are defined. From the mechanism documentation at http://www.neuron.yale.edu/neuron/stati ... racellular, my understanding is that there should a differential equation associated with vext, something along the lines of d_vext/dt = (1/xc)(I_membrane + I_myelin + I_periaxonal) where I_membrane is the current flowing through the cell membrane, I_myelin is the current flowing through the myelin sheath, which should be equal to xg*vext and I_periaxonal is the current flowing along the periaxonal space, which should be equal to xraxial*(vext  vext_neighbor) where vext_neighbor is the neighboring section.
My questions are:
1. Are the equations I've outlined above correct, or is there something else going on in the extracellular mechanism that I'm missing.
2. For converting the "densities" of xraxial, xc, and xg to absolute values (i.e. Mhos/cm^2 to Mhos) does the extracellular mechanism use the standard length and diameter defined for the section? Because we're modeling a myelin sheath the diameter of the sheath is different from the diameter of the section, but I don't see anywhere that this would be defined in the extracellular mechanism.

 Site Admin
 Posts: 5885
 Joined: Wed May 18, 2005 4:50 pm
 Location: Yale University School of Medicine
 Contact:
Re: The math behind extracellular mechanisms
Yes.mesorensen wrote:Are the equations I've outlined above correct
Yes. Don't forget about nseg. If diameter changes along the length of the section, be sure to use diam(x) and area(x) where x is the normalized distance from the 0 end of a given segment.For converting the "densities" of xraxial, xc, and xg to absolute values (i.e. Mhos/cm^2 to Mhos) does the extracellular mechanism use the standard length and diameter defined for the section?
With NEURON xc and xg must be specified relative to the surface area of the neurite itself, not the myelin sheath. For example, given a (cylindrical) neurite with diameter 1 um, and a sheath with external diameter 2 um, the numeric value of the sheath's specific capacitance specified relative to the area of the neurite would be 2 times larger than if it were specified relative to the external surface area of the sheath.Because we're modeling a myelin sheath the diameter of the sheath is different from the diameter of the section, but I don't see anywhere that this would be defined in the extracellular mechanism.
Re: The math behind extracellular mechanisms
Ted,
Thank you very much for your reply. Something still does not appear to line up, and I'm afraid I must be missing something that's going on behind the scenes in NEURON's extracellular mechanism. The time course of Vm for the "internal" section is correct, but the time course of the extracellular potential is much too long for the same value of xc; interestingly, the steadystate value of the extracellular potential is the same as neuron computes. Do you see something I'm missing from the equations below?
Here's the .hoc description of the section I'm working with (I'm disregarding the contribution of neighboring sections right now):
axon{
nseg=1
diam = 1.9
L = 100
Ra = 642.8
cm = 0.66
insert pas
g_pas = 0.00033
e_pas = 80
insert extracellular xraxial=337397 xg=4.16667e06 xc=0.000416667
}
And here are the equations/code I'm trying to use to replicate NEURON's behavior:
diameter = 1.9 //um
length = 100.0 //um
Ra = 642.8 //Ohm*cm
cm = 0.66 //uF/cm^2
g_pas = 0.00033 //S/cm^2
e_pas = 80 //mV
xg = 4.16667e6 //S/cm^2
xc = 0.000416667 //uF/cm^2
//Calculate surface area of section
SA = (3.1415926*diameter*length)/1e8 //cm^2
Cm = cm*SA*1e6 //pF
Cp = xc*SA*1e6 //pF
GPas = g_pas*SA*1e9 //nS
Gmyelin = xg*SA*1e9 //nS
//Calculate Passive currents
IPas = GPas*(Vm  e_pas)
Imyelin = Gmyelin*Vm_ext
//Membrane potential differential equation
Vm' = (1/Cm)*(Istim  IPas)
Vm_ext' = (1/Cp)*(IPas  Imyelin)
Thank you very much for your reply. Something still does not appear to line up, and I'm afraid I must be missing something that's going on behind the scenes in NEURON's extracellular mechanism. The time course of Vm for the "internal" section is correct, but the time course of the extracellular potential is much too long for the same value of xc; interestingly, the steadystate value of the extracellular potential is the same as neuron computes. Do you see something I'm missing from the equations below?
Here's the .hoc description of the section I'm working with (I'm disregarding the contribution of neighboring sections right now):
axon{
nseg=1
diam = 1.9
L = 100
Ra = 642.8
cm = 0.66
insert pas
g_pas = 0.00033
e_pas = 80
insert extracellular xraxial=337397 xg=4.16667e06 xc=0.000416667
}
And here are the equations/code I'm trying to use to replicate NEURON's behavior:
diameter = 1.9 //um
length = 100.0 //um
Ra = 642.8 //Ohm*cm
cm = 0.66 //uF/cm^2
g_pas = 0.00033 //S/cm^2
e_pas = 80 //mV
xg = 4.16667e6 //S/cm^2
xc = 0.000416667 //uF/cm^2
//Calculate surface area of section
SA = (3.1415926*diameter*length)/1e8 //cm^2
Cm = cm*SA*1e6 //pF
Cp = xc*SA*1e6 //pF
GPas = g_pas*SA*1e9 //nS
Gmyelin = xg*SA*1e9 //nS
//Calculate Passive currents
IPas = GPas*(Vm  e_pas)
Imyelin = Gmyelin*Vm_ext
//Membrane potential differential equation
Vm' = (1/Cm)*(Istim  IPas)
Vm_ext' = (1/Cp)*(IPas  Imyelin)

 Site Admin
 Posts: 5885
 Joined: Wed May 18, 2005 4:50 pm
 Location: Yale University School of Medicine
 Contact:
Re: The math behind extracellular mechanisms
I'd like to see this. Send just enough hoc, ses, and mod files to reproduce your findings (in a zipped archive) to ted dot carnevale at yale dot edumesorensen wrote:The time course of Vm for the "internal" section is correct, but the time course of the extracellular potential is much too long for the same value of xc; interestingly, the steadystate value of the extracellular potential is the same as neuron computes.
Re: The math behind extracellular mechanisms
Ted,
Thanks for your help. We've resolved this issue in our code to correspond with NEURON's results. One final question: if you have connected together two sections, each with nseg = 1, and the same diameter, like so:
connect sectionA(1), sectionB(0)
and are injecting a current at sectionA, location 0.5, is NEURON using the passive cable equation to calculate Vm at location 1 on section A, and then using that Vm to calculate the axonal or extracellular current exchange between the two sections? Or with nseg=1 is each section an isopotential compartment?
Thanks again.
Thanks for your help. We've resolved this issue in our code to correspond with NEURON's results. One final question: if you have connected together two sections, each with nseg = 1, and the same diameter, like so:
connect sectionA(1), sectionB(0)
and are injecting a current at sectionA, location 0.5, is NEURON using the passive cable equation to calculate Vm at location 1 on section A, and then using that Vm to calculate the axonal or extracellular current exchange between the two sections? Or with nseg=1 is each section an isopotential compartment?
Thanks again.

 Site Admin
 Posts: 5885
 Joined: Wed May 18, 2005 4:50 pm
 Location: Yale University School of Medicine
 Contact:
Re: The math behind extracellular mechanisms
Considering the simplest case, which omits extracellular
The equivalent circuit of a section with nseg==1 is a T circuit in which all membrane properties (capacitance and ion channels) are lumped into the vertical stem of the T. Read this
https://neuron.yale.edu/neuron/static/p ... p2.htm#3.2
The axial resistance from the 0 to the 1 end of the section constitutes the total resistance of the two arms (top bar) of the T. The ungrounded end of the T's stem is attached to the junction of the two arms. If the section's geometry is specified with the stylized method (i.e. L and diam), the resistances of the two arms of the T will be identical (as if diameter were constant). If geometry is specified with the pt3d method (i.e. a sequence of pt3dadd statements), the arms may have different resistances because the 3d diameters and distances between adjacent pt3d points are taken into consideration.
For the case that you describe, the equivalent circuit is a PI with four nodes. Two of the nodes are the unattached ends of the two sections (0 in A and 1 in B). The other two nodes are the sites of attachment of the lumped membrane properties of the two sections (0.5 in A and B). The resistance between these two nodes is the sum of the axial resistance between 0.5 and 1 in A and the axial resistance between 0 and 0.5 in B.
A comment on this particular use of the connect statement: it makes B the parent of A, and neither the 1 end of A nor the 0 end of B has a parent node. This can lead to user confusion, and may produce unexpected resultsfor example, read about ri in the Programmer's Reference, and then see what you get by executing
forall for (x) print secname(), " ", x, ri(x)
This is avoided by treating parent sections as proximal, child sections as distal,
and
0 ends as proximal, 1 ends as distal. Many of NEURON's GUI tools, and much hoc code, is written with this convention in mind.
Guessing that you intended A to be the parent, the corresponding connect statement would be
connect B(0), A(1)
The equivalent circuit of a section with nseg==1 is a T circuit in which all membrane properties (capacitance and ion channels) are lumped into the vertical stem of the T. Read this
https://neuron.yale.edu/neuron/static/p ... p2.htm#3.2
The axial resistance from the 0 to the 1 end of the section constitutes the total resistance of the two arms (top bar) of the T. The ungrounded end of the T's stem is attached to the junction of the two arms. If the section's geometry is specified with the stylized method (i.e. L and diam), the resistances of the two arms of the T will be identical (as if diameter were constant). If geometry is specified with the pt3d method (i.e. a sequence of pt3dadd statements), the arms may have different resistances because the 3d diameters and distances between adjacent pt3d points are taken into consideration.
For the case that you describe, the equivalent circuit is a PI with four nodes. Two of the nodes are the unattached ends of the two sections (0 in A and 1 in B). The other two nodes are the sites of attachment of the lumped membrane properties of the two sections (0.5 in A and B). The resistance between these two nodes is the sum of the axial resistance between 0.5 and 1 in A and the axial resistance between 0 and 0.5 in B.
A comment on this particular use of the connect statement: it makes B the parent of A, and neither the 1 end of A nor the 0 end of B has a parent node. This can lead to user confusion, and may produce unexpected resultsfor example, read about ri in the Programmer's Reference, and then see what you get by executing
forall for (x) print secname(), " ", x, ri(x)
This is avoided by treating parent sections as proximal, child sections as distal,
and
0 ends as proximal, 1 ends as distal. Many of NEURON's GUI tools, and much hoc code, is written with this convention in mind.
Guessing that you intended A to be the parent, the corresponding connect statement would be
connect B(0), A(1)
Re: The math behind extracellular mechanisms
Thanks, Ted. We should be able to figure this out now.
Re: The math behind extracellular mechanisms
Ted
Thanks again for all your help. On further examination, we're still having trouble getting the extracellular mechanisms working correctly. I have several questions right now, but I'll try taking them one at a time to see if I can troubleshoot this without taking too much of your time.
The first question is how does the inclusion of an extracellular layer in NEURON effect the evolution of a section's internal membrane potential? That is, for a normal passive section with currents Iaxonal (the currents from neighboring sections), Istim (externally driven stimulus current), and Ipas (passive current with g_pas and e_pas), the membrane capacitor current should be:
ICmem = Istim  Ipas  Iaxonal
and dVm/dt = (1/Cmem)*ICmem
where
Ipas = g_pas*(Vme_pas)
and Iaxonal is:
sum((Vm  Vm_neighbor)/Raxonal)
When we include an extracellular layer around the section, how do the equations governing ICmem and dVm/dt change? For example, would Ipas change to:
Ipas = g_pas*(Vm  Vm_ext  e_pas)
to account that the driving force is reduced/augmented by the voltage of the external section?
Thanks again for all your help. On further examination, we're still having trouble getting the extracellular mechanisms working correctly. I have several questions right now, but I'll try taking them one at a time to see if I can troubleshoot this without taking too much of your time.
The first question is how does the inclusion of an extracellular layer in NEURON effect the evolution of a section's internal membrane potential? That is, for a normal passive section with currents Iaxonal (the currents from neighboring sections), Istim (externally driven stimulus current), and Ipas (passive current with g_pas and e_pas), the membrane capacitor current should be:
ICmem = Istim  Ipas  Iaxonal
and dVm/dt = (1/Cmem)*ICmem
where
Ipas = g_pas*(Vme_pas)
and Iaxonal is:
sum((Vm  Vm_neighbor)/Raxonal)
When we include an extracellular layer around the section, how do the equations governing ICmem and dVm/dt change? For example, would Ipas change to:
Ipas = g_pas*(Vm  Vm_ext  e_pas)
to account that the driving force is reduced/augmented by the voltage of the external section?

 Site Admin
 Posts: 5885
 Joined: Wed May 18, 2005 4:50 pm
 Location: Yale University School of Medicine
 Contact:
Re: The math behind extracellular mechanisms
First, an important item of terminology: there is no "external section."mesorensen wrote:When we include an extracellular layer around the section, how do the equations governing ICmem and dVm/dt change? For example, would Ipas change to:
Ipas = g_pas*(Vm  Vm_ext  e_pas)
to account that the driving force is reduced/augmented by the voltage of the external section?
Second, an important fact: the intracellular potential at any location is the sum of the transmembrane potential and the extracellular potential immediately adjacent to the external surface of the cell at that location. It is possible to write the node equations so as to exploit this fact.
NEURON solves the node equations, which can be derived by referring to the equivalent circuit and applying Kirchhoff's current law at each node. For the standard case (extracellular adds two layerssee
http://www.neuron.yale.edu/neuron/stati ... racellular
) each intracellular node has two extracellular counterparts. Proceeding from the outside in, the potentials at the ith nodes are vext[1]_i, vext_i, and v_i + vext_i.
For a passive model, assume that:
indices increase from left to right
the ith axial resistors are to the left of the ith nodes
all resistances, conductaces, and capacitances are in absolute units (not density units).
Adopt the following conventions:
use D as the derivative operator
don't bother enclosing incremented or decremented indices in parentheses
omit the node indices of cm, g_pas, e_pas, xc, xg, xc[1], xg[1], and e_extracellular (which are all understood to be i)
take advantage of the fact that, for the ith node triplet, the absolute potential at the innermost node is vext_i + v_i
The equations for the ith node triplet, with capacitive terms on the left hand side of each equation, can then be written "by inspection":
v_i * cm*D
= v_i1 * (1/Ra_i)  v_i * ((1/Ra_i) + (1/Ra_i+1)) + v_i+1 * (1/Ra_i+1)  g_pas * (v_i  e_pas)
 v_i * cm*D + vext_i * (cm + xc)*D  vext[1]_i * xc*D
= vext_i1 * (1/xraxial_i)  vext_i * ((1/xraxial_i) + (1/xraxial_i+1)) + vext_i+1 * (1/xraxial_i+1) + g_pas * (v_i  e_pas)  xg * (vext_i  vext[1]_i)
 vext_i * xc*D + vext[1]_i * (xc + xc[1])*D
= vext[1]_i1 * (1/xraxial[1]_i)  vext[1]_i * ((1/xraxial[1]_i) + (1/xraxial[1]_i+1)) + vext[1]_i+1 * (1/xraxial[1]_i+1) + xg * vext_i  xg[1] * (vext[1]_i  e_extracellular)
Alternatively, the first two equations could be written first in terms of the absolute potential at the intracellular nodes, after which one would substitute vext + v for the absolute potential and rearrange to get the form shown above.
Although I have proofread this several times, it would not surprise me if these equations contain one or more errors (especially in the indices).
Re: The math behind extracellular mechanisms
Ted
That got it for us; we're now only getting slight discrepencies that are probably due to the use of a different integration method. Thanks again for all your help.
Michael
That got it for us; we're now only getting slight discrepencies that are probably due to the use of a different integration method. Thanks again for all your help.
Michael
Re: The math behind extracellular mechanisms
Hello,
I've been working with extracellular for some weeks now to simulate electrical stimulation and I assumed I had understood it. However I just thought of a question I could not answer. I believe the answer should easy, especially for you, but apparently I'm just a little too slow today :)
Referring to viewtopic.php?f=28&t=168&p=378 and estim.zip from viewtopic.php?f=13&t=212&p=4447&hilit=ecstim.zip#p4447:
You mention Rattay who described that the extracellular stimulation can be simulated by a current source proportional to the second spatial derivative of the extracellular field. Could you please explain why for the extracellular mechanism a linear e_extracellular gradient along the axon still works?
Thank you!
I've been working with extracellular for some weeks now to simulate electrical stimulation and I assumed I had understood it. However I just thought of a question I could not answer. I believe the answer should easy, especially for you, but apparently I'm just a little too slow today :)
Referring to viewtopic.php?f=28&t=168&p=378 and estim.zip from viewtopic.php?f=13&t=212&p=4447&hilit=ecstim.zip#p4447:
You mention Rattay who described that the extracellular stimulation can be simulated by a current source proportional to the second spatial derivative of the extracellular field. Could you please explain why for the extracellular mechanism a linear e_extracellular gradient along the axon still works?
Thank you!

 Site Admin
 Posts: 5885
 Joined: Wed May 18, 2005 4:50 pm
 Location: Yale University School of Medicine
 Contact:
Re: The math behind extracellular mechanisms
A very good question. The short answer is that a linear extracelluar gradient also works as an effective stimulus in the real world, so the activating function approach can't be valid under all circumstances, can it?Could you please explain why for the extracellular mechanism a linear e_extracellular gradient along the axon still works?
If you can lay your hands on a copy of Rattay's book
Electrical Nerve Stimulation (Springer, New York, 1990)
you'll find his derivation of the activating function on pages 122124, but near the bottom of page 124 comes
By first approach he really means first approximation. Also note that on page 131 the last sentence of paragraph 2 saysIn unmyelinated fibers, a first approach of the influence of extracellular electrodes is given by the activating function with is the second derivative of the extracellular potential along the fiber.
where "these points of interest" means the points in space at which the discretized cable equation is integrated in time. So in settings with complex extracellular geometries he abandons the notion of the activating function and falls back on . . . the same approach that I used with NEURON, i.e. represent extracellular potential explicitly.In the case of inhomogeneous media . . . or for complicated shapes of electrodes, we will only calculate the extracellular voltage at these points of interest.
"But what's complex about a 'uniform extracellular field' (i.e. a field with constant potential gradient)?"
Nothing, really. The complexity lies in the cell. The activation function derivation made several tacit assumptions. Chief among these is that it assumes the axon is infinitely long. This allowed Rattay to completely ignore boundary conditions at the ends of the axon. And, in theory, it is true that an infinitely long axon exposed to a constant extracellular field along its entire length will remain at resting potential.
However, if you look at what happens to a finite axon that is exposed to a uniform extracellular field, you will see that its membrane potential is at rest in the middle, but it diverges from rest as you approach its sealed ends. For long axons, the slope of v vs. distance is far from constant, approaching 0 at the middle but increasing rapidly in the near vicinity of the ends.
Re: The math behind extracellular mechanisms
Thank you very much for the detailed answer! That clarifies a lot!
Re: The math behind extracellular mechanisms
Hi Ted,
I'am trying to derive the equations from one of your previous posts myself, in order to understand them, but I can't figure out what exactly the variables v, vext, vext[1] are.
v+vext is the potential difference between inneraxonal and the first layer, vext the difference between first and second layer and so forth. then v should be the potential difference between the inneraxonal and the outermost layer.
but then I can't derive the same equations then you did.
Am I right with my assumptions? and what is it that I am doing wrong?
Thanks,
Simon
I'am trying to derive the equations from one of your previous posts myself, in order to understand them, but I can't figure out what exactly the variables v, vext, vext[1] are.
v+vext is the potential difference between inneraxonal and the first layer, vext the difference between first and second layer and so forth. then v should be the potential difference between the inneraxonal and the outermost layer.
but then I can't derive the same equations then you did.
Am I right with my assumptions? and what is it that I am doing wrong?
Thanks,
Simon

 Site Admin
 Posts: 5885
 Joined: Wed May 18, 2005 4:50 pm
 Location: Yale University School of Medicine
 Contact:
Re: The math behind extracellular mechanisms
The Programmer's Reference documentation of extracellular spells this out and includes a picture.sdanner wrote:I'am trying to derive the equations from one of your previous posts myself, in order to understand them, but I can't figure out what exactly the variables v, vext, vext[1] are.
No. It is the intracellular potential relative to the bath (AKA "ground" in the Programmer's Reference figure). Please see the documentation.v+vext is the potential difference between inneraxonal and the first layer
No. It is the extracellular potential adjacent to the cell membrane relative to the bath. See the documentation.vext the difference between first and second layer and so forth
No. It is the membrane potential, i.e. the potential difference between the internal and external surfaces of the neuronal membrane. This is the potential difference that is responsible for the electrical field that bends channel proteins. Channel proteins aren't affected by the potential difference between the inside of the neuronal membrane and the bath, unless the bath is right up against the outside surface of the neuronal membrane.then v should be the potential difference between the inneraxonal and the outermost layer.