Finding middle of a section

When Python is the interpreter, what is a good
design for the interface to the basic NEURON
concepts.

Moderator: hines

Post Reply
ssrihari
Posts: 4
Joined: Mon Dec 30, 2019 11:05 pm

Finding middle of a section

Post by ssrihari »

Hi! I am currently working on an algorithm for finding the min. current needed to produce an action potential using extracellular stimulation and the MRG axon as my model. I used h.APCount(<insert section name>(0.5)) to monitor whether the voltage in the middle of my sections crosses a certain threshold. However, I noticed I was getting varying results depending on the method I use to determine the coordinates of the middle of my sections (which I need in order to determine their distance from the current source and thus the voltage to be assigned to that section's e_extracellular). For example, when I use this method:

Code: Select all

if (sec.n3d() % 2 == 0):
        midUp = sec.n3d() / 2
        midDown = midUp - 1
        secCoor = Point((sec.x3d(midUp) + sec.x3d(midDown)) / 2e3, (sec.y3d(midUp) + sec.y3d(midDown)) / 2e3, (sec.z3d(midUp) + sec.z3d(midDown)) / 2e3)
    else:
        mid = sec.n3d() / 2
        secCoor = Point(sec.x3d(mid) / 1e3, sec.y3d(mid) / 1e3, sec.z3d(mid) / 1e3)

    r = calcDistance(source, secCoor)
I get 2.1 mA for my current threshold (min. current necessary for an AP). Whereas when I use this method:

Code: Select all

xSum = 0
    ySum = 0
    zSum = 0

    for i in range(sec.n3d()):
        xSum += sec.x3d(i)
        ySum += sec.y3d(i)
        zSum += sec.z3d(i)
    
    return Point(float(xSum) / sec.n3d(), float(ySum) / sec.n3d(), float(zSum) / sec.n3d())
I get 3.1 mA for my current threshold even though the distances are extremely similar. Finally, when I try using x_xtra, y_xtra, and z_xtra to find the midpoints, my apc's (APCount) find no action potentials occured even for large current ranges (e.g. 0 - 100 mA). My question then is which of these corresponds to the true middle of a section (sec(0.5))?
ramcdougal
Posts: 267
Joined: Fri Nov 28, 2008 3:38 pm
Location: Yale School of Public Health

Re: Finding middle of a section

Post by ramcdougal »

Neither of those strategies necessarily identifies the coordinates of the center of the section. The second solution doesn't necessarily lead to a point on the section (consider a circular section with evenly spaced 3D points; it would find the center of the circle).

The explanation is simple: 3D points don't necessarily have any connection to the numerical discretization. In particular, you could have all but one 3D point near the 0 end; their purpose is only to define the shape.

The function coordinate in the below returns the (x, y, z) value of the section(x) that is passed in, as demonstrated in the below by plotting evenly along a non-straight path:

Code: Select all

from collections import namedtuple
import numpy as np
from matplotlib import pyplot as plt
from neuron import h

Point = namedtuple('Point', ['x', 'y', 'z'])

def coordinate(seg):
    sec = seg.sec
    n3d = sec.n3d()
    d = seg.x * sec.L
    xs = [sec.x3d(i) for i in range(n3d)]
    ys = [sec.y3d(i) for i in range(n3d)]
    zs = [sec.z3d(i) for i in range(n3d)]
    arcs = [sec.arc3d(i) for i in range(n3d)]
    return Point(np.interp(d, arcs, xs), np.interp(d, arcs, ys), np.interp(d, arcs, zs))

soma = h.Section(name='soma')
soma.pt3dadd(0, 0, 0, 1)
soma.pt3dadd(5, 5, 0, 2)
soma.pt3dadd(10, 5, 0, 2)

for i in range(11):
    pt = coordinate(soma(i * 0.1))
    plt.scatter(pt.x, pt.y)

plt.axis('equal')
plt.show()
(Importantly, this is not necessarily the center of the corresponding segment; seg.x returns the x value used to create the segment not the normalized poisition at the center of the segment... this is why one segment sufficed for the above code.)

Ignoring all that: calling coordinate(my_sec(0.5)) will return the (x, y, z) coordinate of the center of my_sec.
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Finding middle of a section

Post by ted »

I am currently working on an algorithm for finding the min. current needed to produce an action potential using extracellular stimulation and the MRG axon
You say nothing about the geometry of the stimulating electrode(s) (point electrodes vs. electrodes with finite surface area, monopolar vs. bipolar vs. multipolar stimulation), the time course of stimulating current), or the location of the axon relative to the location of the stimulating electrode(s). The following comments assume monopolar stimulation with a short duration current pulse (about 0.1 ms).

APCount isn't smart. It can't tell whether a threshold crossing is caused by an action potential or is simply stimulus artifact. This is why it is best to attach the APCount to a node of Ranvier that is far away from the stimulating electrode.
I was getting varying results depending on the method I use to determine the coordinates of the middle of my sections (which I need in order to determine their distance from the current source and thus the voltage to be assigned to that section's e_extracellular).
Middle of your sections? That should be middles of the model's segments, even if you are certain that no section in the model has nseg>1. Even if the model is "designed so that each section has nseg == 1", you will want to test for adequacy of spatial discretization--which requires multiplying nseg of at least all "long" sections (internodes in this case) by an odd factor and running a new simulation to see how changing nseg affects simulation results.

This code

Code: Select all

if (sec.n3d() % 2 == 0):
        midUp = sec.n3d() / 2
        midDown = midUp - 1
        secCoor = Point((sec.x3d(midUp) + . . .
    else:
        mid = sec.n3d() / 2
        . . .
appears to have been written under the assumption that the extracellular field will affect only those segments that contain a section's midpoint, or, in the case of a section with even nseg, the segments that are adjacent to the section's midpoint. That assumption is incorrect. Every segment that is in the conductive medium is exposed to the extracellular field, and its e_extracellular must be controlled.

This code

Code: Select all

xSum = 0
    ySum = 0
    zSum = 0

    for i in range(sec.n3d()):
        xSum += sec.x3d(i)
        . . .
is guaranteed to produce incorrect values for the xyz coordinates of segment centers. Why? For one thing, it will fail dramatically in a model that has branching.

"Well, the MRG axon isn't branched, so . . . "

Even in an unbranched model it risks producing completely incorrect values for the xyz coordinates of segment centers. Read the documentation of define_shape() in the Programmer's Reference, especially this sentence:
Sections that already have pt3d info are translated to ensure that their first point is at the same location as the parent.
If a sections' geometry is defined by pt3dadd statements in which the first xyz coordinate is 0,0,0, then as soon as h.define_shape() is called (or if a PlotShape is created, either via the GUI or a hoc or Python statement), that section's pt3d data will all be changed automatically. That's a "feature" in the sense that it makes shape plots render the shape of the cell correctly, but it breaks your algorithm for calculating the locations of xyz points.

So it's best to call h.define_shape() preemptively, right after all sections have been created, connected, assigned their geometries via L and diam or pt3dadd statements, and their nseg parameters have been set. And after you do that, every section will have its own pt3d data that will include the correct xyz coordinates of its 0 and 1 ends and of the centers of each of its segments. Ready to be used. No need to loop over its pt3d data and add anything up.

What happens if you then change any section's nseg or its length? Will that section's pt3d data be revised automatically? I'm not sure--but why take a chance? Just call h.define_shape() again and let NEURON figure it out.
ramcdougal
Posts: 267
Joined: Fri Nov 28, 2008 3:38 pm
Location: Yale School of Public Health

Re: Finding middle of a section

Post by ramcdougal »

So it's best to call h.define_shape() preemptively, right after all sections have been created, connected, assigned their geometries via L and diam or pt3dadd statements, and their nseg parameters have been set. And after you do that, every section will have its own pt3d data that will include the correct xyz coordinates of its 0 and 1 ends and of the centers of each of its segments.
Note: h.define_shape() does not create new points for sections whose morphology was defined by calling the pt3dadd method.

For example,

Code: Select all

from neuron import h

soma = h.Section(name='soma')
soma.pt3dadd(0, 0, 0, 1)
soma.pt3dadd(5, 0, 0, 1)
soma.nseg = 5

h.define_shape()

print(soma.n3d())  # prints: 2
shows "2" despite there being 5 segments; there are no points at the center of the segments.
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Finding middle of a section

Post by ted »

h.define_shape() does not create new points for sections whose morphology was defined by calling the pt3dadd method.
True. My mistake.

The way to discover the xyz locations of segment centers is to treat a section's pt3d data as piecewise linear functions of path length from the section's 0 end, and apply linear interpolation to those functions at path distances from the 0 end that correspond to
L*(0.5 + i/nseg)
where i ranges from 0 to nseg-1.
Post Reply