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.