Intracellular resistance and potentials

[BoshuoW]

Hello all!

My questions are:

1. How does NEURON handle axial/intracellular resistance in its internal calculation, and how would setting intracellular resistivity/resistance to very small values (Ra<<default) affect the accuracy and stability of the simulation? Could very short pulse duration (<10 us) interact with such settings to produce unexpected behaviors (see descriptions below)?
2. Is it possible or are there better ways to specify equal intracellular potential for multiple compartments (connected sections each with nseg =1) than setting a very small intracellular resistivity/resistance for them?

The context of my questions is as follows:

I was exploring a modeling method described by Sergeev, E. N., Meffin, H., Tahayori, B., Grayden, D. B., and Burkitt, A. N. in "Effect of soma polarization on electrical stimulation thresholds of retinal ganglion cells." (2013 6th International IEEE/EMBS Conference on Neural Engineering (NER), pp. 1135––1138, 2013. http://ieeexplore.ieee.org/abstract/document/6696138/)

The main idea is to replace the single compartment soma with multiple compartments to simulate extracellular stimulation of transverse field (Fig.2).



As the author stated, "the multicompartmental soma was represented by a discretized sphere, with all the somatic compartments connected to the root vertex at the center of the sphere, by negligible axial resistances". This approximates the condition of equal intracellular potentials of the compartments while their extracellular potentials are different, thus representing patches of the somatic membrane with different transmembrane potentials.

I built a simpler model - a soma consisting of 15 compartments with HH membrane and without any dendritic and axonal connections. I tested stimulation thresholds under uniform extracellular fields for a range of pulse duration and compared results with a my own simulation in MATLAB (following Boinagrov, D., Loudin, J., and Palanker, D. "Strength–Duration Relationship for Extracellular Neural Stimulation: Numerical and Analytical Models." Journal of Neurophysiology, 104(4), pp. 2236–2248, 2010.) The difference for the MATLAB simulation is that the equal intracellular potential condition was explicitly enforced and not approximated. The multicompartment model in NEURON produced thresholds that agreed with the MATLAB simulation for most of the pulse widths, however, behaved differently for very short pulses.

Specifically, thresholds obtained from MATLAB and NEURON agreed within ~2% difference for pulses between 10 us to 10 ms duration, and the sub-chronaxie part of the strength-duration (S-D) curve (< hundreds us) plotted on log-log axes followed a linear relationship with negative slope. However, thresholds from NEURON started to deviate from this relationship when pulse widths decreased shorter than 10 us: the threshold values first increased to about 10%-30% larger compared to MATLAB, and then abruptly jumped to 2-3 orders of magnitude larger for pulses shorter than 3 us. I have suspected that the low intracellular resistivity was the problem, and tried to vary Ra. The behavior described above seems to be insensitive to Ra values between 1e-7 up to the default of 35.4 (which is not a negligible value at all). Less than 1e-8, simulation would not run, and the error message seems to indicate division by zero. Increasing Ra, from 35.4*10 up to 35.4 *1000, the thresholds for those short pulses decreased and came closer to the S-D curve (on the same order of magnitude); however, at this point the assumptions and behavior of the simulation have completely changed from the intended model.

Any insights and suggestions are appreciated.

[ted]

Remember, you're dealing with a computational representation of physical system which is continuous in space and time (when viewed at the scale of microns and microseconds). The continuous system is described by a PDE, but you're using numerical solution of a set of ODEs to approximate the solution to the PDE, and that requires discretizing both space and time. Spatial discretization approximates the PDE by a family of ODEs, where each ODE describes the time course of state variables at a different point in space. Accuracy of the entire solution depends not only on dt, the interval at which you sample time, but also on dx, the interval at which you sample space. Make dt very small, and you will expose spatial errors, which become predominant at short times.

So you either have to make dx very small (more, smaller compartments) or you have to "cheat" by eliminating space (which is what you did with Matlab by forcing the intracellular potential to be the same in all compartments). With NEURON you may be able to achieve the same effect by implementing your 15 compartments as completely independent sections ("cells"), each with nseg=1, then using the Linear Circuit Builder (or, equivalently, LinearMechanism) to connect their internal nodes at 0.5 to a common node. The Linear Circuit Builder's GUI might be adequate for the task, or if you can figure out the documentation of the LinearMechanism class you might just be able to write out the necessary code without having to resort to the GUI (maybe helped along by using the Builder to set up a toy model with just 2 or 3 sections, and using the Builder's ability to export code, which might give you enough of a hint how to generalize that to 15 sections).

[BoshuoW]

Hello Ted, thanks for your clarification and suggestion!

I'll explore both the Linear Circuit Builder and Linear Mechanism and see whether it works.

[ted]

A simpler/easier hack that I'd try first: put off dealing with LinearCircuit for now. Instead, connect the 0 ends of the 14 smaller sections to the 0.5 node of the largest section, and (1) see how small an Ra you can get away with, or (2) make all sections very "thin" (i.e. short--think blini). Just make sure that each section has the correct surface area and its vext is what you want. Doesn't force all sections to have identical internal potentials, but it may allow a closer approximation than what results from discretizing a single section into 15 segments.

[BoshuoW]

I have actually tried this "spoke-hub" topology as well. The results were almost identical as the other topology, with the same issue for short pulses. Although, I have not experimented so much with changing Ra (set only to 1e-3) or the section's size (set to cylinders of equal length and diameter first, then another attempt with thin/flat shapes by decreasing/increasing them 10 times) yet for this topology. I'll see whether smaller Ra or thinner sections improves the situation.

[BoshuoW]

I implemented the model with LinearMechanism, and unfortunately the S-D curve was pretty much identical compared to before. Therefore, I don't think that the very low axial resistance was the problem, but something else contributed to the behavior with short pulse duration and small time steps.

[ted]

Another possible explanation for peculiar results is that there may be some problem with NMODL-specified ion channel mechanisms that are used in the model. Opportunities for bugs abound--NMODL syntax is a bit obscure, and the equations for voltage-dependent rate constants (or steady state and time constant values) can themselves be a problem, e.g. singularities may occur (situations in which numerator & denominator fall to 0 at a particular voltage). The only way to rule that out is for an experienced NMODL programmer to examine the mod files. I'd be glad to do that for you, if you zip them up and email them to
ted dot carnevale at yale dot edu

[BoshuoW]

Ted, thank you for the suggestion!! I have doubted the HH.mod file as well, but wasn't very confident to change it. After all, this is the first time I worked with hoc code. I have adjusted the .mod file and the problem has been resolved even without using LinearMechanism. Thresholds could be obtained for all pulse duration down to 1 us and the results agreed with my MATLAB simulation!

Here are the changes I made:

Since the extracellular potentials could well exceed the original range of the table, I changed the table size and resolution for the voltage, which also significantly sped up the simulations

Code: Select all

TABLE minf, mtau, hinf, htau, ninf, ntau DEPEND celsius FROM -1000 TO 1000 WITH 5001

The vtrap function was changed to remove the -x/y/2 factor. Based on my own calculation using the L'Hôpital's rule, the limit of the vtrap with x->0 should be just y, which is the way I implemented the channels in MATLAB . The if statement wasn't changed, but probably could also be adjusted to (fabs(x) < 1e-3).

Code: Select all

FUNCTION vtrap(x,y) {  :Traps for 0 in denominator of rate eqns.
        if (fabs(x/y) < 1e-6) {
                vtrap = y :*(1 - x/y/2)
        }else{
                vtrap = x/(exp(x/y) - 1)
        }
}

[ted]

Not sure what HH.mod file you're referring to. Don't expect file names to be unique identifiers of anything. If you don't say where you found HH.mod, or present the contents of that file, nobody but you will know what's in it. The only exceptions might be the files that are distributed with NEURON, e.g. the particular hh.mod that is included--but even then you have to tell people that you're referring to the hh.mod that is distributed with NEURON.

With regard to the modifications you made, changing the range or resolution of a lookup table wouldn't be my first guess about how to fix a stability problem--but then I probably wouldn't be running simulations in which extracellular potential was so large in the first place.

The vtrap function was correct in its original form, which approximated x/(exp(x/y) - 1) by a straight line segment with the appropriate slope in the near neighborhood of |x/y|==0. You just replaced that approximation with something that will cause a discontinuity of a voltage-dependent gating parameter when v lands within a narrow band centered on x/y==0. Will that cause problems for you or anyone who reuses your code in the future? Maybe, if adaptive integration is used--the adaptive integrators require voltage-dependent parameters to be continuous functions of membrane potential, but the occasional violation of continuity may not always cause a noticeable problem.

[BoshuoW]

Sorry for not being clear. Since I was using the default HH channels in NEURON, I indeed meant the HH.mod file distributed with NEURON, which I located under the folder: \src\nrnoc.

Thanks for clarifying the vtrap function and continuity requirement. I thought that the 1e-6 band was narrow enough that the single value would be OK. I've change it back. And based on more testing, I think the table could be the main reason for the issue I had before. When I adjusted the model so that the threshold for a given pulse duration would increase (by changing temperature or gNabar), similar issue with short pulses occurred again, which could also be fixed by further increasing the range of the table. I don't think using large extracellular potentials by itself is a problem, since a large range of V_ext might be necessary and totally fine if a neuron spans a large distance. Having the extracellular potentials vary over very short distances could be problem.

[ted]

Since I was using the default HH channels in NEURON, I indeed meant the HH.mod file distributed with NEURON, which I located under the folder: \src\nrnoc.

The name of the file distributed with NEURON is hh.mod, not HH.mod. Case is important.