Axonal Conduction Velocity: A Computer Study

Arthur D. Snider¹, Aman Chawla², Salvatore D. Morgera¹

¹Department of Electrical Engineering, University of South Florida, Tampa, USA.
²Department of Computer Science and Engineering, Indian Institute of Technology, New Delhi, India.

DOI: 10.4236/jamp.2024.121007                

Abstract

This paper derives rigorous statements concerning the propagation velocity of action potentials in axons. The authors use the Green’s function approach to approximate the action potential and find a relation between conduction velocity and the impulse profile. Computer simulations are used to bolster the analysis.

Keywords

Neuron, Axon, Action Potential, Conduction Velocity, Internode

1. Introduction

The nerve in the CNS of vertebrates with nervous systems is composed of extended parts called axons. These specializations are responsible for carrying the information integrated in the cell body and delivering it to the post-synaptic cells. If one were to look at the nervous system as a network of neurons, with the cell bodies being the computing nodes and the axons being the communication links, then the bottleneck in network information throughput clearly is at these links. More importantly, it is limited by the conduction speed of the signals down the axons. Modeling and measuring the conduction velocity of axonal impulses correctly, therefore acquires a great significance. Furthermore, in [1] [2] [3] [4] [5] , it was shown that in the case of nerves with ephaptic coupling, precise adjustments to this speed take place. Thus, it becomes imperative to clearly understand the underlying single-axon impulse speed, especially in an analytical manner.

The Hodgkin-Huxley equation has no known analytic solution and has a complicated mechanism at the nodes. There have even been attempts in the literature to simplify the nodal equations for better understanding [6] [7] . The present paper takes the problem from the opposite angle, and studies the solutions to this equation analytically and numerically and by zooming in on time. This investigation provides the insight that a “wave” is a misnomer for the propagation solution and thereby calls for an alternate understanding of the propagating impulse; one such alternate view is also suggested.

Much work has been done on the analysis of the wave properties of the action potential as it propagates down an axon [8] - [19] , and the nomenclature has become so diverse that some semantic coincidences have spawned inappropriate generalizations. This paper is intended to restore a proper perspective and derive rigorous statements concerning the propagation velocity.

The paper is structured as follows. In Section 2, we investigate conduction velocity in an axon. We conclude in Section 3. An Appendix provides more details related to the simulations performed. Before we move to the next section, we provide a table (Table 1) summarizing the notation used in the paper succinctly.

2. Action Potential Conduction Velocity

Conduction velocity of nerve fibers has long been a subject of intense interest. See, for example, the review at the beginning of [20] . In the present paper, we go over the conduction velocity formulation in the context of myelinated nerve fibers required for computer simulation purposes and come to a conclusion as to the best method of computing it from simulation output.

2.1. Preliminaries

We first summarize the physical configuration under study. An axon is a long, slender projection of a nerve cell, or neuron, which conducts electrical impulses away from the neuron’s cell body (see Figure 1). It is composed of relatively long sections sheathed in myelin, interrupted by short nodes of Ranvier. Ideally, the axon is straight and cylindrically symmetric and the myelinated sections are longitudinally homogeneous. The electrical impulses are exhibited as voltage drops $V\left(x\mathrm{,}t\right)$ from the center line of the axon to the external medium, with x measuring distance along the axon and t denoting time.

The voltage $V\left(x\mathrm{,}t\right)$ satisfies the cable equation:

$C\frac{\partial V}{\partial t}=\frac{1}{R}\frac{{\partial }^{2}V}{\partial x^2}-GV+I$ (1)

where the transversal capacitance C, longitudinal resistance R, and leakage conductance G are uniform (but different) in each of the two regions (myelinated and Ranvier); I denotes any injected current, plus the ion currents generated solely in the nodes of Ranvier in accordance with the Hodgkin-Huxley model [12] (or later variations [9] ). The parameters in Equation (1) are per-unit-length.

2.2. Wave Interpretation of Simulation Results

Simulations were obtained by applying the Crank-Nicholson scheme to Equation (1). We provide numerical details in the Appendix. In Figures 2-5, we display simulated voltage profiles at selected times for a 4 cm length of axon with nine

Table 1. Symbols and their meanings as used in the main text and Appendix. Units and values are indicated in the third column, where appropriate.

nodes of Ranvier, each having length 2.5 × 10⁻⁴ cm and spaced 0.2 cm apart. To eliminate the effect of boundary conditions (Neumann on the left end and Dirichlet on the right end) we stimulated the central node of Ranvier and discontinued the simulation when the boundary voltages became appreciable. Thanks to symmetry, we only need to display the profiles to the right of the central node.

Figure 1. Neuron, with a focus on the axon cross-section. The nodes of Ranvier are marked with short red arrows and the myelin sheath with longer red ones.

Figure 2. Simulated action potential, spreading to the right. The horizontal axis is location and the vertical axis is voltage. Here, $t_1<t_2<t_3<t_4$ .

Figure 3. Simulated action potential, exciting the second node of Ranvier. Here, $t_5<t_6<t_7<t_8$ .

The figures display the stimulated action potential ( $t_1$ ), spreading out ( $t_2\mathrm{,}{t}_3\mathrm{,}{t}_4$ ) until it starts to stimulate the adjacent node of Ranvier ( $t_5$ ), eventually triggering an action potential at that node ( $t_6$ ). This action potential spreads out ( $t_7\mathrm{,}{t}_8$ ), triggering another action potential at the next downrange node ( $t_9\mathrm{,}{t}_{10}$ ). And in turn further downrange nodes are triggered. (Note that, in agreement with experiment, the subsequent action potentials only spawn further potentials downrange, because the preceding nodes of Ranvier need a recovery time before

Figure 4. Simulated action potential, exciting the third node of Ranvier. Here, $t_9<t_{10}<t_{11}<t_{12}$ .

Figure 5. Simulated action potential exciting the fourth node of Ranvier and then triggering the fifth node. Here, $t_{13}<t_{14}<t_{15}<t_{16}$ .

they can be triggered again.)

It is a little misleading to call this propagation a “wave,” since the profile does not retain its shape as it propagates down the myelinated segments. However, if we only display the profiles at the instants of triggering, we get the apparent wave profiles in Figure 6.

Possibly, the misleading Figure 6 has led some investigators to conjecture that the voltage $V\left(x\mathrm{,}t\right)$ satisfies the wave equation:

$\frac{{\partial }^{2}V}{\partial t^2}=v^2\frac{{\partial }^{2}V}{\partial x^2}$ (2)

whose solutions are superposed shape-invariant functions $V\left(x\mathrm{,}t\right)=f\left(x\pm vt\right)$ . Indeed, it has often been proposed to combine Equation (2) with the cable equation, Equation (1), in the hopes of deriving an expression for the propagation velocity v. The futility of this strategy is highlighted by the following observations:

1) The only simultaneous solutions of Equation (2) and Equation (1) (with $G=I=0$ ) are constant or unbounded:

$V\left(x\mathrm{,}t\right)=M+N{\text{e}}^{RCv\left(vt-x\right)}+P{\text{e}}^{RCv\left(vt+x\right)}\mathrm{<mo>\; ;\; </mo>}$ (3)

Figure 6. Action potential profiles at triggering instants. Note the equality in the heights of the impulses.

Figure 7. Comparison of spatial and temporal second derivatives of the transmembrane potential, $\frac{{\partial }^{2}V}{\partial t^2}$ and $\frac{{\partial }^{2}V}{\partial x^2}$ . Here, the horizontal axis represents position.

2) The graphs of $\frac{{\partial }^{2}V}{\partial t^2}$ and $\frac{{\partial }^{2}V}{\partial x^2}$ are definitely not proportional (Figure 7).

Indeed, examination of Figures 2-5 underscores the difficulty of assigning a “propagation velocity” to the phenomenon. The profiles do not evolve as “waves” within a myelinated segment at all; they resemble more a tidal “swell,” or diffusion spreading. In fact, the solution of the cable Equation (1) with $I=0$ , evolving from an initial delta-function shape supported at $x_0$ , is given by [17] .

$V_{Green}\left(x\mathrm{,}t\right)=\frac{{\text{e}}^{-Gt/C}{\text{e}}^{-{\left(x-x_0\right)}^{2}RC/4t}}{\sqrt{\pi t}}$ (4)

(This is, essentially, the Green’s function for the diffusion equation.) Figure 8 displays how closely the profiles of the action potential are tracked by a judiciously scaled replica of the Green’s function. For this calculation the nodes of Ranvier have been removed for simplicity.

Figure 8. Simulated voltage profiles (black) and Green’s function (blue). Note the close overlap between the two curves, in each sub-figure. The progression in time is from the left-most sub-figure to the right-most sub-figure.

2.3. The Conduction Velocity

Since the voltage profiles only occasionally appear as shifted replicas of themselves, it would appear that the wisest strategy for assigning a conduction velocity v to the spreading of the action potential is suggested by Figure 6; we take v to be the ratio of the distance L between two consecutive Ranvier nodes, divided by the time lapse $\tau$ between the triggering of action potentials at the nodes.

The distance L, of course, is well defined. To demonstrate the feasibility of specifying the time lapse $\tau$ , let us assume for the moment that the birth of an action potential at a node of Ranvier occurs when the voltage at that node $V\left(x\mathrm{,}t\right)$ achieves a critical value $V_c$ . And we assume that the action potential is well approximated by a scaled version of the Green’s function, Equation (4). Then the triggering event is described by the equation:

$V_{Green}\left(\tau \right)=\frac{K{\text{e}}^{-G\tau /C}{\text{e}}^{-L^{2}RC/4\tau }}{\sqrt{\pi \tau }}=V_c$ (5)

and the velocity is given formally by:

$v=\frac{L}{\tau }=\frac{L}{V_{Green}^{-1}\left(V_c\right)}$ (6)

A plot of $V_{Green}\left(\tau \right)$ versus $\tau$ demonstrates that Equation (5) produces a well-defined time lapse $\tau$ (Figure 9).

Obviously, (5) does not predict a constant $\frac{L}{\tau }$ ratio; the conduction velocity varies with the myelinated segment length L.

But the point is made: a triggering condition based on voltage $V\left(L\mathrm{,}\tau \right)=V_c$ , on longitudinal current $I\left(L,\tau \right)=-\frac{1}{R}\frac{\partial V\left(L,\tau \right)}{\partial x}=I_c$ , or on charge accumulation $Q={\displaystyle {\int }_0^{\tau }I\left(L,t\right)\text{d}t}=Q_c$ [14] specifies a well-defined time lapse $\tau$ and thus a conduction velocity. But this computation need not be based on crude approximations like (5) (which makes no attempt to allow for different values of the electric parameters in the nodes of Ranvier); it can be gleaned from simulations. One can visually identify the times at which the action potential spikes peak, and deduce $\tau$ by subtraction.

Figure 9. Plot of the Green’s function with respect to the second argument, $V_{Green}\left(L\mathrm{,}\tau \right)$ . The critical value of transmembrane voltage at which an action potential is triggered is indicated on the vertical axis, as $V_c$ .

2.4. Observations

In this subsection, we list a few observations based on our simulations.

1) For a given internode distance L, the time lapse between every consecutive pair of action potential triggerings appears to be the same.

2) However, the velocity v does vary with internode distance L, apparently possessing a maximum. See Figure 10, compiled from simulations.

3) Which parameter (V, I, or Q) triggers the action potential upon achieving a critical value? Some authors suggest it is accumulated charge Q (which may be equivalent to voltage, through capacitance). Action potential initiation is very difficult to detect precisely from simulations, because once an action potential has been initiated it dominates, and thus distorts, the voltage profile. Also, there is the issue of where the parameter should be measured. We circumvented the issue of when the action potential is triggered by simply noting the instants of peak voltage.

3. Conclusions

This article addresses a simplified model of a neuron, admitting the approximation of the action potential’s progression by a Green’s function for the linear, constant coefficient cable equation (Equation (1))—an approximation justified solely by visual inspection of Figure 8. Clearly, more sophisticated models can be evoked, but at the cost of sacrificing the transparency of the analytic nature of the formulation. Moreover, the issue of whether the birth of an action potential is triggered by critical values of the voltage, charge, or current should be addressed;

Figure 10. Plot of the conduction velocity of an action potential versus internodal length. The horizontal axis represents internodal length in centimeters and the vertical axis represents conduction velocity in centimeters per second.

but note that one can retain the idealized analytic formulation for investigating any of these criteria. However, the immediate focus of the study is to redirect the misleading paradigm of the potential as a wave.

To conclude, specialists recognize that the cable equation is of the parabolic type (diffusive), and not hyperbolic (wave). The propagation between nodes of Ranvier is more of a surge than a wave. The “chain reaction’’ nature of the action potentials is only suggestive, and does not justify attributing any of the features of the wave equation to the phenomenon. By discussing these features rigorously, we hope to have illuminated them in a way that is useful to computational neuroscientists and mathematicians alike.

Appendix: Details of the Simulations

The simulations were obtained by applying the Crank-Nicolson scheme to the cable equation simplified with an integrating factor:

$\frac{\partial \left[{\text{e}}^{Gt/C}V\right]}{\partial t}=\frac{1}{RC}\frac{{\partial }^2\left[{\text{e}}^{Gt/C}V\right]}{\partial x^2}+\frac{{\text{e}}^{Gt/C}I}{C}$ (7)

employing mesh sizes $\delta x=0.02$ cms, $\delta t=2\text{e}-6$ seconds (and confirmed with $\delta x=0.04$ cm, $\delta t=8\text{e}-6$ seconds). Overall dimensions and boundary conditions are described in what follows; relevant values used in the simulations reported in the Appendix are shown in Table 1.

The injected current was −4 × 10⁻⁵ A/cm for a duration of 4 × 10⁻⁶ seconds. No leakage current was presumed at the nodes of Ranvier. The sodium, potassium and “nonspecific’’ (p) ion currents at these nodes are expressed in terms of Hodgkin-Huxley activation/inactivation variables m, n, h and p (which are dimensionless) as follows: $j_{Na}=\pi d{P}_{Na}{m}^{2}hZ\left(Na\right)$ , $j_K=\pi d{P}_K{n}^{2}Z\left(K\right)$ , $j_p=\pi d{P}_p{p}^{2}Z\left(p\right)$ A/cm, where, with $Y=Na,K$ , we define:

$Z\left(Y\right)=\frac{F^2\left(V-0.07\right)}{\kappa T}\times \frac{{\left[Y\right]}_0-{\left[Y\right]}_i{\text{e}}^{F\left(V-0.07\right)/\kappa T}}{{\text{e}}^{F\left(V-0.07\right)/\kappa T}-1}$

in units of Coulombs/cubic centimeter. Here, the resting potential is −0.07 volts, T is the absolute temperature (25 K + 273.15 K), F is the Faraday constant (96,485 Coulombs/mole), and $\kappa$ is the gas constant (8.3145 J/K-mole). $P_{Na}\mathrm{,}{P}_K\mathrm{,}{P}_p$ denote the permeabilities. The outer and inner concentrations are taken as ${\left[Na\right]}_o=1.145\times {10}^{-4}$ , ${\left[Na\right]}_i=1.374\times {10}^{-5}$ , ${\left[K\right]}_o=2.5\times {10}^{-6}$ , ${\left[K\right]}_i=1.2\times {10}^{-4}$ moles/cm³. The activation/inactivation variables each satisfy [9] [12] :

$\frac{\text{d}y}{\text{d}t}=\left[{\alpha }_y\left(1-y\right)-{\beta }_{y}y\right]\rho \left(T\right)$ (8)

where $y=\left(m\mathrm{,}n\mathrm{,}h\mathrm{,}p\right)$ . Here, $\rho \left(T\right)$ is a temperature factor (with T = 298.15 K) and the coefficients ( ${\alpha }_y\mathrm{,}{\beta }_y$ ) are voltage-dependent [ $V=V\left(x\mathrm{,}t\right)$ ] functions whose tabulations have persistently been reported with incorrect dimensions in the literature, even in some of the references cited in this paper. In our simulations, we used:

${\alpha }_h=\frac{A_{\ast }\left(B_{\ast }-V\right)}{1-{\text{e}}^{\left(V-B_{\ast }\right)/C_{\ast }}}$ , ${\beta }_h=\frac{A_{\ast }}{1+{\text{e}}^{\left(B_{\ast }-V\right)/C_{\ast }}}$ ,

${\alpha }_y=\frac{A_{\ast }\left(V-B_{\ast }\right)}{1-{\text{e}}^{\left(B_{\ast }-V\right)/C_{\ast }}}$ , ${\beta }_y=\frac{A_{\ast }\left(B_{\ast }-V\right)}{1-{\text{e}}^{\left(V-B_{\ast }\right)/C_{\ast }}}$ , $y=m\mathrm{,}n\mathrm{,}p$ ;

$A_{{\alpha }_h}=1\text{e}5$ , $A_{{\alpha }_m}=3.6\text{e}5$ , $A_{{\alpha }_n}=2\text{e}4$ , $A_{{\alpha }_p}=6\text{e}3$ per volt-s

$A_{{\beta }_h}=4.5\text{e}3$ per s

$A_{{\beta }_m}=4\text{e}5$ , $A_{{\beta }_n}=5\text{e}4$ , $A_{{\beta }_p}=9\text{e}4$ per volt-s

$B_{{\alpha }_h}=-1\text{e}-2$ , $B_{{\alpha }_m}=2.2\text{e}-2$ , $B_{{\alpha }_n}=3.5\text{e}-2$ , $B_{{\alpha }_p}=4\text{e}-2$ volts

$B_{{\beta }_h}=4.5\text{e}-2$ , $B_{{\beta }_m}=1.3\text{e}-2$ , $B_{{\beta }_n}=1\text{e}-2$ , $B_{{\beta }_p}=-2.5\text{e}-2$ volts

$C_{{\alpha }_h}=6\text{e}-3$ , $C_{{\alpha }_m}=3\text{e}-3$ , $C_{{\alpha }_n}=1\text{e}-2$ , $C_{{\alpha }_p}=1\text{e}-2$ volts

$C_{{\beta }_h}=1\text{e}-2$ , $C_{{\beta }_m}=2\text{e}-2$ , $C_{{\beta }_n}=1\text{e}-2$ , $C_{{\beta }_p}=2\text{e}-2$ volts.

We adopted values for the permeabilities from [9] as follows: $P_{Na}=8\text{e}-3$ cm/s, $P_p=0.54\text{e}-3$ cm/s, $P_K=1.2\text{e}-3$ cm/s. If we regard $V\left(x\mathrm{,}t\right)$ as constant over time intervals of length $\delta t$ , the coefficients in Equation (8) are also constant and its solutions satisfy:

$y\left(t+\delta t\right)=y\left(t\right){\text{e}}^{-\left({\alpha }_y+{\beta }_y\right)\rho \delta t}+\left(1-{\text{e}}^{-\left({\alpha }_y+{\beta }_y\right)\rho \delta t}\right)\frac{{\alpha }_y}{{\alpha }_y+{\beta }_y}$ ; (9)

we use this formula in the simulations.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Chawla, A. and Morgera, S.D. (2014) Ephaptic Synchronization as a Mechanism for Selective Amplification of Stimuli. BMC Neuroscience, 15, Article No. P87. https://doi.org/10.1186/1471-2202-15-S1-P87
[2]	Scott, A.C. (2002) Neuroscience: A Mathematical Primer. Springer-Verlag, New York.
[3]	Chawla, A., Morgera, S. and Snider, A. (2019) On Axon Interaction and Its Role in Neurological Networks. IEEE/ACM Transactions on Computational Biology and Bioinformatics, 18, 790-796.
[4]	Chawla, A., Morgera, S.D. and Snider, A.D. (2021) Fields, Geometry, and Their Impact on Axon Interaction. Journal of Applied Mathematics and Physics, 9, 751-778. https://doi.org/10.4236/jamp.2021.94053
[5]	Chawla, A. and Morgera, S.D. (2022) On the Geometry of Interacting Axon Tracts in Relation to Their Functional Properties. https://doi.org/10.1101/2022.03.29.486317
[6]	FitzHugh, R. (1955) Mathematical Models of Threshold Phenomena in the Nerve Membrane. The Bulletin of Mathematical Biophysics, 17, 257-278. https://doi.org/10.1007/BF02477753
[7]	Markin, V.S. (1970) Electric Interaction of Parallel Nonmyelinized Nerve Fibers. I. Change in the Excitability of an Adjacent Fiber. Biofizika, 15, 120-129.
[8]	Debanne, D., Campanac, E., Bialowas, A., Carlier, E. and Alcaraz, G. (2011) Axon Physiology. Physiological Reviews, 91, 555-602. https://doi.org/10.1152/physrev.00048.2009
[9]	Frankenhaeuser, B. and Huxley, A.F. (1964) The Action Potential in the Myelinated Nerve Fibre of Xenopus laevis as Computed on the Basis of Voltage Clamp Data. The Journal of Physiology, 171, 302-315. https://doi.org/10.1113/jphysiol.1964.sp007378
[10]	Halter, A. (1989) A Distributed-Parameter Model for Myelinated Nerve Fibre. Ph.D. Thesis, Rice University, Texas, AAT9136095.
[11]	Goldman, L. and Albus, J.S. (1968) Computation of Impulse Conduction in Myelinated Fibers; Theoretical Basis of the Velocity-Diameter Relation. Biophysical Journal, 8, 596-607. https://doi.org/10.1016/S0006-3495(68)86510-5
[12]	Hodgkin, A.L. and Huxley, A.F. (1952) A Quantitative Description of Membrane Current and Its Application to Conduction and Excitation in Nerve. The Journal of Physiology, 117, 500-544. https://doi.org/10.1113/jphysiol.1952.sp004764
[13]	Landahl, H.D. and Podolsky, R.J. (1949) On the Velocity of Conduction in Nerve Fibre with Saltatory Transmission. The Bulletin of Mathematical Biophysics, 11, 19-27. https://doi.org/10.1007/BF02477911
[14]	Malmivuo, J. and Plonsey, R. (1995) Bioelectromagnetism: Principles and Applications of Bioelectric and Biomagnetic Fields. Oxford University Press, Oxford. https://doi.org/10.1093/acprof:oso/9780195058239.001.0001
[15]	Reutskiy, S., Rossoni, E. and Tirozzi, B. (2003) Conduction in Bundles of Demyelinated Nerve Fibers: Computer Simulation. Biological Cybernetics, 89, 439-448. https://doi.org/10.1007/s00422-003-0430-x
[16]	Rushton, W.A.H. (1951) A Theory of the Effects of Fibre Size in Medullated Nerve. The Journal of Physiology, 115, 101-122. https://doi.org/10.1113/jphysiol.1951.sp004655
[17]	Snider, A.D. (2006) Partial Differential Equations: Sources and Solutions. Dover Publications, Mineola.
[18]	Waxman, S.G. and Wood, S.L. (1984) Impulse Conduction in Inhomogeneous Axons: Effects of Variation in Voltage-Sensitive Ionic Conductances on Invasion of Demyelinated Axon Segments and Preterminal Fibers. Brain Research, 294, 111-122. https://doi.org/10.1016/0006-8993(84)91314-3
[19]	Waxman, S.G. (1988) Biophysical Mechanisms of Impulse Conduction in Demyelinated Axons. Advances in Neurology, 47, 185-213.
[20]	Matsumoto, G. and Tasaki, I. (1977) A Study of Conduction Velocity in Nonmyelinated Nerve Fibers. Biophysical Journal, 20, 1-13. https://doi.org/10.1016/S0006-3495(77)85532-X