current is zero in some segments

[mflynn]

I noticed a problem when I looked at the axial currents in my model. There are segments where the current is zero. The problem is not with the model, since I found the same thing in a simple, unbranched passive cylindrical neuron. When I inject current into one end of the neuron I found that the voltage doesn't decay smoothly along the the length of the neuron even though there is no change in the axial resistance. Even though the segments are all of equal length, the voltage in some of the segments adjacent to each other are identical. I can't figure out why the voltage is not strictly decreasing.

[ted]

Can you provide a very brief example? One section should do it, two if you absolutely must,
please keep nseg as small as possible, pas mechanism only should be sufficient.

[mflynn]

Hi Ted,
Here's the simplest example of where I see the problem:

Code: Select all

xopen("$(NEURONHOME)/lib/hoc/noload.hoc")
tstop = 15
steps_per_ms = 1
dt = 1e-2

create axon
access axon

nseg = 5
diam = 10
L = 500

insert pas
ra        = 100

rm        = 30000

objref stim
stim = new IClamp(0)
stim.del = 5
stim.amp = 10
stim.dur = 0.5

The voltage difference between segments 3 and 4 is zero, so there's no current flowing between them, which would seem to violate Kirchoff's current law. Is it an issue with the integration routine? Perhaps a higher order integration method would be more accurate and not produce these anomolies?

Thanks,
Mark

[csh]

Code: Select all

steps_per_ms = 1
dt = 1e-2

This will result in a very low resolution of your plots. Try

Code: Select all

dt = 1e-2
steps_per_ms = 1/dt

to plot 1 point per time step.

Code: Select all

ra        = 100

hoc is case-sensitive. Axial resistance is Ra.

Code: Select all

rm        = 30000

rm is not a valid hoc keyword. If you use the pas mechanism for passive (leak) conductance, you should use g_pas, which is in (S/cm2). If you meant 30000 Ohm cm2, you could simply write

Code: Select all

g_pas = 1/rm

. In addition, it's helpful to start simulation at rest so that membrane potential is constant at the beginning. You could achieve this with

Code: Select all

e_pas = -70
v_init = e_pas

The voltage difference between segments 3 and 4 is zero,

It is zero at t == 0, because voltage is uniform across the whole section. However, as soon as your IClamp starts at t == 5 ms, voltage will be different in the two segments. You could check this with the following code:

Code: Select all

xopen("$(NEURONHOME)/lib/hoc/stdrun.hoc")
tstop = 5.5 // end of current injection
dt = 1e-2
steps_per_ms = 1/dt

create axon
access axon

nseg = 5
diam = 10
L = 500
Ra        = 100

insert pas
rm        = 30000
g_pas = 1/rm
e_pas = -70
v_init = -70

objref stim
stim = new IClamp(0)
stim.del = 5
stim.amp = 10
stim.dur = 0.5

run()

for (x) if (x != 0 && x != 1) print "Voltage at ", x," = ", axon.v(x)

Hope that helps,
Christoph

[mflynn]

I tried changing steps_per_ms and it didn't change the qualitative results. I'm recording the voltage at the ends of the segments. Is that not the correct way to get an accurate spatial distribution of voltages? I use the voltage difference between the beginning and end of a segment and the resistance of that segment to calculate the axial current.

Mark

[ted]

Here's the simplest example of where I see the problem:

In addition to the technical issues that Christoph noted, consider the following:

Because steps_per_ms is 1, you're not observing membrane potential during the
current pulse. The first measurement isn't until 0.5 ms after the end of the current.
This matters because the physical dimensions and biophysical properties of the model
are such that, in the absence of injected current, longitudinal voltage gradients dissipate
extremely quickly--this is a fat leaky cable with membrane time constant of 1 ms, and its
higher order time constants (which govern relaxation of longitudinal voltage gradients)
are _much_ shorter than 1 ms.

To see this in excruciating detail, temporarily abandon writing hoc code and instead use
the GUI to set up the model and create a space plot of membrane potential along the
length of axon. Be sure to set Init to -70, so you don't have to wait 5 ms after the start of
the simulation for v to settle to its resting level. It's OK to leave dt = 0.025 and
Points plotted/ms = 40. Use the Movie Run tool, and turn on Keep Lines in the space plot.
Result: plenty of longitudinal gradient during the current pulse, vanishes quickly afterwards.

[ted]

I use the voltage difference between the beginning and end of a segment and the resistance of that segment to calculate the axial current.

NEURON uses the central difference approximation to discretize the cable equation. This
means that it computes v at the middle of each segment, not at the ends (see chapter 5 of
The NEURON Book, or
Hines, M.L. and Carnevale, N.T. The NEURON simulation environment.
Neural Computation 9:1179-1209, 1997
preprint available from http://www.neuron.yale.edu/neuron/bib/nrnpubs.html).
Consequently v(x) returns v at the middle of the segment that contains x. If x lies on the
boundary between segments, the returned value will depend on roundoff error (i.e. whether
you get the value from the "lower" node or the "upper" node).

[mflynn]

Because steps_per_ms is 1, you're not observing membrane potential during the
current pulse. The first measurement isn't until 0.5 ms after the end of the current.
This matters because the physical dimensions and biophysical properties of the model
are such that, in the absence of injected current, longitudinal voltage gradients dissipate
extremely quickly--this is a fat leaky cable with membrane time constant of 1 ms, and its
higher order time constants (which govern relaxation of longitudinal voltage gradients)
are _much_ shorter than 1 ms.

Yes, I thought of that after I had posted my question and used a higher steps_per_ms. I still see this problem. If I picked voltage points in the middle of each segment, then the voltage between adjacent segments would be different and hence non-zero currents. From reading the User's Forum posts, I thought the correct way to calculate currents was to measure voltages at each end of a segment and use the axial resistance in that segment. Is that incorrect?

Thanks,
Mark

[csh]

From reading the User's Forum posts, I thought the correct way to calculate currents was to measure voltages at each end of a segment and use the axial resistance in that segment. Is that incorrect?

If you literally mean the beginning and end of one and the same segment, this is wrong, because voltage is defined at a node in the middle of a segment and considered to be uniform throughout this segment (neglecting, for a moment, the slightly more complicated situation at the terminal ends of a section). You will have to use the voltage difference between the middle of adjacent segments to calculate the axial current. The for (x) idiom will take care of that by returning the middle of each segment in relative units of total section length.

[mflynn]

Yes, that worked. I thought I had read in the User's Forum that I should use the voltage at the ends. So, then I should also use the axial resistance in the middle as well?

Mark

[csh]

So, then I should also use the axial resistance in the middle as well?

The specific axial resistance Ra is a section variable that is constant throughout a section. There is a helper function called ri(x) that makes it easier to calculate the axial resistance between two adjacent nodes (http://www.neuron.yale.edu/neuron/stati ... ry.html#ri). You are of course free to implement your own version from ri = Ra * l / A, where l is the length and A is the cross-sectional area between the two adjacent segments you are interested in.

[ted]

voltage is defined at a node in the middle of a segment and considered to be uniform throughout this segment

I hesitated to interject these pearls for fear of further muddying the already turbid waters,
and running the risk of fomenting new misunderstandings. But here they are anyway.

A fine point: the central difference method produces a numerical result that is second
order correct in space. This means that the value of v at any location between adjacent
internal nodes can be found by linear interpolation. For example, in a section with
nseg = 5, NEURON comuptes v to second order precision at 0.1, 0.3, 0.5, 0.7, and 0.9.
If one wishes to know v(0.4), it can be calculated as v(0.3) + 0.5*(v(0.5) - v(0.3)).

Another fine point: current density is treated as if it were constant along the length of any segment.

Enough.

[ted]

I thought I had read in the User's Forum that I should use the voltage at the ends.

If such a statement is somewhere in this Forum, it is untrue or taken out of context.

[mflynn]

So, then I should also use the axial resistance in the middle as well?

The specific axial resistance Ra is a section variable that is constant throughout a section. There is a helper function called ri(x) that makes it easier to calculate the axial resistance between two adjacent nodes (http://www.neuron.yale.edu/neuron/stati ... ry.html#ri). You are of course free to implement your own version from ri = Ra * l / A, where l is the length and A is the cross-sectional area between the two adjacent segments you are interested in.

Yes, but is the x in ri(x) in the middle of the segment where I'm measuring the voltage or is it at the ends?

Thanks,
Mark

[ted]

is the x in ri(x) in the middle of the segment where I'm measuring the voltage or is it at the ends?

RTM
http://www.neuron.yale.edu/neuron/stati ... ry.html#ri