A non-local model of the propagation of action potentials in myelinated neurons

Myelinated neurons are characterized by the presence of myelin, a multilaminated wrapping around the axons formed by specialized neuroglial cells. Myelin acts as an electrical insulator and therefore, in myelinated neurons, the action potentials do not propagate within the axons but happen only at the nodes of Ranvier which are gaps in the axonal myelination. Recent advancements in brain science have shown that the shapes, timings, and propagation speeds of these so-called saltatory action potentials are controlled by various biochemical interactions among neurons, glial cells and the extracellular space. Given the complexity of brain's structure and processes, the work hypothesis made in this paper is that non-local effects are involved in the optimal propagation of action potentials. A non-local model of the action potentials propagation in myelinated neurons is proposed that involves spatial derivatives of fractional order. The effects of non-locality on the distribution of the membrane potential are investigated using numerical simulations.


1-Introduction
The human brain is made of billions of neurons and glial cells separated by the extracellular space (ECS) that facilitates a huge number of cellular interconnections needed for the proper functionality of brain. Neurons are the most important brain cells for the cognitive functions of the brain since they provide brain's functional abilities by transmitting information throughout the body via action potentials. Understanding the fundamental science that governs neuronal processes is crucial for the development of reliable diagnoses and treatments of brain diseases.
While a lot of progress has been made in neuroscience especially in the last few decades when the advancements in technology allowed for much improved observations of cellular and molecular phenomena in the living brain, important knowledge about neuronal processes that could translate into successful treatments of some brain diseases with increasing incidence such as epilepsy, multiple sclerosis and Alzheimer's disease continues to be missing. The main difficulty in achieving this much needed knowledge is the complex nature of neurons and their interactions. Neuronal complexity comes from the fractal-like branched structure of neurons and the highly non-linear processes taking place inside and outside the neurons [1]. Neurons have the following characteristics specific to complex systems: flexibility, adaptability, emergence and noise [1,2]. Of particular interest are emergence and noise. Emergence is a spontaneous event that cannot be predicted from knowledge of the structure and interactions of neurons [1]. According to [3], (membrane) noise is the collection of random fluctuations in the potential of the neuronal membrane which is responsible for action potentials. These fluctuations are caused by the Brownian motion of ions and electrons due to thermal variations (Gaussian white noise), discrete nature of the electric currents through the membrane (shot noise), and conductance changes due to the random opening and closing of ion channels. The noise, emergence and the connections between them represent strong non-linear dynamics which are non-local in space and time [2]. Fluctuation effects at the pico-and nanometric length scales of the ion channels and the ECS cannot be neglected. Dinariev [4,5] proved mathematically that thermal fluctuations can lead to non-local hydrodynamics. Also, at these length scales the longrange interactions among ions and between ions and water molecules are significant [6,7]. Thanks to water's network of intermolecular hydrogen bonds, nanometric dipolar correlations among the water molecules develop at nanoscales which is a non-local dielectric feature of water [8]. Adding ions to water produces solvation shells and experiments and molecular dynamics simulations suggest that the ions and their first hydration shells behave as colloidal suspensions in pure liquid water [9] which can have enhanced long-range interactions at nanoscales.
Another noise-inducing process that could lead to spatio-temporal non-locality is anomalous diffusion of ions described as fractional Brownian motion (fractional Gaussian noise) or Lévy flights (Lévy-stable noise) [10]. For instance, non-locality in neuronal dynamics may be facilitated by the ECS that surrounds the extensively interconnected brain cells ( figure 1(a)). The ECS is a collection of narrow channels (20-60 nm width) with convoluted geometries that are filled with interstitial fluid, free and fixed long-chain macromolecules that form the extracellular matrix, and ions involved in cellular resting, action potentials and release of transmitter molecules [11,12]. This space is a constantly fluctuating environment since action potentials and some biochemical interactions among neurons, glial cells and blood capillaries are mediated by the ECS. Diffusion of ions through the ECS might be anomalous due to various factors such as 1) the geometry of the ECS, 2) dead spaces caused by local enlargements of the ECS, voids or glial wrapping around cells, 3) obstruction by the extracellular matrix, 4) biding sites on cell membranes of the extracellular matrix, and 5) the presence of fixed negative charges on the extracellular matrix [12] ( figure 1(b)). The fixed negative charges regulate ion mobility [13]. Glial cells, especially astrocytes, also control the ion and water flow in the ECS [14,15]. In particular, the stellate-shaped astrocytes have numerous specialized endfeet linking with other glial cells, neurons, and brain-fluid barriers that are endowed with ion and water channels for the regulation of ion and water homeostasis. For example, the regulation of the extracellular potassium concentration * by astrocytes is achieved through two separate processes: 1) an active uptake, accumulation and release of potassium, and 2) passive spatial buffering where the uptake at one location is coupled to the release at a different location. Since both processes cause astrocyte swelling, an intricate process of volume regulation is activated in the astrocyte that will release ion and water in the ECS [14]. We suspect that the tight control of the ion and water movement by astrocytes may also contribute to the ECS non-locality.
An additional source of spatio-temporal non-locality could be the anomalous diffusion of ions through the myelin sheath of neurons. In a myelinated neuron the axon is made of myelinated regions separated by the nodes of Ranvier ( figure 2). Myelin is a multilaminated wrapping around the axons made of layers of glial plasma membrane. Oligodendrocytes are specialized glial cells that produce myelin and induce (through processes regulated by neurons and astrocytes) the clustering of voltage-gated sodium channels at the nodes of Ranvier and the aggregation of fast voltagegated potassium channels in the myelin sheath [16,17]. A few slow potassium channels might also exist in the nodes of Ranvier [18]. Some experiments in cell cultures and in vivo have shown that nodal-like clusters appear just before the deposition of myelin along axons begins [19]. During an action potential, sodium flows inside the neuron at the node of Ranvier, and potassium is released into the periaxonal space which is the extracellular space between the axon and the myelin. To avoid the accumulation of potassium and osmotically driven water in the periaxonal space, potassium diffuses through the myelin's potassium channels and the gap junctions connecting the myelin layers and is removed by the astrocytes endfeet that form gap junctions with the outer most myelin layer [14]. Thus, in myelinated neurons, the action potentials do not propagate within the axons but happen only at the nodes of Ranvier. Experiments performed on the plasma membrane of a cultured human embryonic kidney cell showed that the voltage-gated potassium channel Kv2.1 displays anomalous diffusion [20]. Given that the Kv2.1 channel shows similar clustering patterns in the membranes of cultured human embryonic kidney cells and native neurons [20], it is possible that anomalous diffusion of potassium may also happen in the myelin sheath of neurons. * A disruption in the potassium regulation is a sign of neuronal network dysfunction [20]. For instance, during an epileptic seizure, the extracellular potassium concentration increases. In this paper we assume that the spatio-temporal non-locality due to, among other factors, long-range interactions of ions and water molecules at nanoscales and the anomalous diffusion of ions through the myelin sheath and the ECS, has a non-negligible effect on the propagation of action potentials in myelinated neurons. We propose a mathematical model of the action potential propagation that accounts for spatial non-local effects. To model spatial non-locality, we use fractional calculus which has been successfully employed to describe spatio-temporal non-locality in various complex systems (see for instance [21,22] and references within, and [23] for anomalous diffusion). We generalize the classic cable equation that models the spatio-temporal propagation of action potentials in neurons. For now, we derive the model for mature myelinated neurons which do not have backpropagating action potentials [24], and thus the unidirectional propagation of the action potentials will be represented by asymmetric operators. We start by introducing the non-local voltage as a convolution whose kernel is a decaying power law of fractional order. This kernel function can model the Lévy flights/fractional Brownian motion characterizing anomalous diffusion [10,23] and the superposition of exponentially decaying screened Coulomb potentials representing long-range electrostatic interactions of ions [25]. An optimal approximation of a decaying power law by exponentials can be found in [26]. We further derive a non-local cable equation using fractional order derivatives and the same steps as in [27] where the classic cable equation is presented. By replacing the spatial derivatives of integer order in the cable equation (which is a one-dimensional parabolic partial differential equation) with spatial fractional Caputo derivatives of fractional order we obtain a onedimensional integro-differential equation. We notice that such an approach has not been proposed in the literature until now. The fractional cable equation proposed in [28,29] uses a temporal fractional order derivative to introduce a diffusion operator in the classic cable equation. This equation models only temporal anomalous electro-diffusion in neurons, not spatial non-locality.
The structure of the paper is as follows. Some mathematical preliminaries needed in the paper are given in section 2. The proposed mathematical model is presented in section 3, together with a numerical solution to the non-local cable equation coupled with the Hodgkin-Huxley equations. Numerical simulations are shown in section 4. Lastly, the paper ends with a section of conclusions and further work.

2-Preliminaries
In this section we present some fractional calculus results that will be used in the paper. The mathematical concepts we need in this paper are taken from [30][31][32][33][34].

Direction of Propagation of Action Potentials
Axon and the right-sided Riemann-Liouville fractional integral of order is: is the gamma function. In particular, + 0 ( ) = − 0 ( ) = ( ).
3. The left (right)-sided Caputo fractional derivative of order , with − 1 < < , = 2,3, …, is said to be a left(right)-sided sequential Caputo fractional derivative if: Since in our model we will use only the left-sided fractional operators to describe the unidirectional propagation of action potentials, some useful properties of these operators are listed in the following proposition.

3-Mathematical Model
In this section we introduce a generalization of the cable equation using spatial derivatives of fractional order. The unidirectional propagation of action potentials from the cell body to the axon terminals is modeled using left-sided fractional order operators.
For now, we do not consider time-dependent spatial non-local effects. In addition, since the axonal branches are long and narrow, the action potentials depend only on one spatial variable [27] and thus the problem is one-dimensional. Therefore, we model the axon as a circular cylinder of constant radius and one characteristic length , the length of the internodal region, and assume that: On the other hand, formula (11) gives: where the Caputo fractional derivative is taken with respect to the spatial variable .
Thus, as Δ → 0, formulas (15) and (16) yield: The uniformity of the current density in every cross-sectional area of the cylindrical neuron gives [27]: where is the intracellular resistance and is the electric current. By combining formulas (17) and (18) we get the expression of the longitudinal current: where the negative sign comes from the current sign convention. Also, replacing formula (18) into formula (15) yields the Ohm's law:

Δ ( , ) =
for the generalized electric resistance: Further we denote by: As in the classic case [1], we introduce the following currents: where and are the membrane and external currents, respectively, and and are their corresponding currents per unit area. The capacitor current is and the specific membrane capacitance is .
If we replace now formulas (19) and (21) in the Kirchhoff's law in the element shown in figure 3b: As Δ → 0: which replaced in equation (22) A numerical discretization of the second order spatial derivative in equation (24) shows that in myelinated neurons the voltage at node depends only on the voltages at the adjacent nodes of Ranvier − 1 and + 1. However, a numerical discretization of the spatial fractional derivative in equation (23) based on the Grünwald-Letnikov formula (see later) shows that the voltage at node depends on the voltages at all the previous nodes 0, 1, … , − 1, and at node + 1.
Near the resting potential the membrane current per unit area can be approximated as [27]: is the voltage relative to the resting potential, and is the specific membrane resistance. Since is constant it follows from formula (8) where we denoted by: Here is a generalized electrotonic length that shows the role of parameter ∈ (0,1) as a non-local corrector of the relationship among the geometric parameters , and the resistances , of a neuron. We notice that equation (23) remains valid for a time-varying ( ), thus allowing the modeling of dynamic effects on , , , and through varying . This might provide relevant insights into how aging and neuro-glial processes control the myelin assembly so that the timing, speed and shape of the action potentials are optimized. For instance, by relating ( ) to aging we might be able to understand the variations in the lengths and diameters of the internodal regions and nodes along one axon or between two neurons [36][37][38]. It may also be possible to link parameter ( ) to a mechanical model of swelling to account for the change in the axon's diameter during an action potential [39], or to the diffusion of the extracellular potassium to account for changes in the membrane resistance [18].
In the following sub-section, we will present a numerical solution to equation (25) coupled with the Hodgkin-Huxley equations. The problem is stated and solved in the internodal region (figure 3a).

3-1-Dynamic Problem
We consider a leaky internodal region with one isopotential node described by the modified Hodgkin-Huxley model proposed in [40]. The reason why we chose the Hodgkin-Huxley model instead of the more commonly used (see for instance [41,42]) model by Frankenhaeuser and Huxley [43] is that the Frankenhaeuser-Huxley model provides an incorrect description of the sodium currents [3]. The membrane potential of the internodal region is solution to the following initial boundary value problem: Above is an externally applied current per unit area, , , and are the reverse potentials, , , , , and are respectively the voltage-gated maximal conductances of + and + , and the leak conductances of + , + , and − . Lastly, the gating variables , , and ℎ represent the activations of the + and + channels and the inactivation of the + channel, respectively. Let: where ̃( , ) is solution to the following boundary value problem: Then the following identities hold: The solution to problem (29) can be found as follows. Integrating one the differential equation gives: The above identity can be transformed further into: ̃( , ) −̃(0+, ) = 0+ ( 0+̃) = 0+ = Γ( + 1) by using formulas (9) and (6). Finally, using the boundary conditions of problem (29) in formula (31) we find the solution to problem (29): From the definition of ( ) it follows that ̃( , 0) = 0 which implies that ( , 0) = 0. Thus an initial boundary value problem equivalent to problem (29) is: We use now the numerical scheme proposed in [44] to discretize the spatial derivatives. The numerical scheme is based on a shifted Grünwald-Letnikov formula that is stable. Let 0 = 0 < 1 < ⋯ < −1 < = be an equally spaced discretization of the interval [0, ] of constant step size Δ . Then the spatially discretized equation (33) with zero Dirichlet boundary conditions is: System (34) with the zero initial condition of problem (33) together with problem (27) are further solved numerically using Matlab's built-in function ode15s. This function solves systems of stiff first order ordinary differential equations by combining a modified linear multistep backward difference formula of order up to 5 and an adaptive step size that changes according to an algorithm that calculates relative and absolute error tolerances [45]. In our simulations we kept the default values of ode15s for the relative error tolerance (10 −3 ) and the absolute error tolerance (10 −6 ). Lastly, the spatial step size Δ was chosen such that the following stability condition for system (34) is satisfied: where Δ is the time step.

4-Results
The parameters used in our numerical simulations are given in Table 1. For our simulations we chose = 1000 which is within the range of values for the internodal region found in the literature (for instance, [46] reports = 200 ÷ 2000 where the lower (higher) values correspond to initial myelination (full maturity) stage [36]). Also, according to the results in [37], a maximum propagation speed of an action potential is reached around a critical value of 1000 for the internodal length. In all the numerical simulations we used a time step ∆ = 0.01.  Figure 4 shows the well-known solution of problem (27), while figure 5 highligths the differences in the spatiotemporal distribution of the membrane potential ( , ), solution to problem (26), for = 1 (classic case; figure 5(a)) and for = 0.65 ( figure 5(b)). It is apparent that near = 0 the amplitude of oscillations of the membrane potential are higher for = 0.65 than for = 1, while near the node of Ranvier ( = ) the oscillations are higher for = 1. Profiles through the three-dimensional plots in figure 5 are shown in figure 6. At fixed time ( = 0.1 ), the potential increases sharper near the node of Ranvier for = 0.65 than for = 1 (figure 6(a)), while at the middle of the internodal region, the amplitude of potential's oscillations is slightly higher for = 1 than for = 0.65 ( figure 6(b)).
Figures 7-10 show profiles of ( , ) in the case when the potential is assumed to be zero at the middle of the internodal revion and symmetric at the node of Ranvier. As → 1, the action potential loses its sharpness at the node for the fixed time = 0.1 (figure 7(a)), while its time variations are almost identical at the fixed location = 700 ( figure 7(b)). The results shown in figures 8 and 9 are obtained for = 0.65 (figures 8(a) and 9(a)) and for = 1 (figures 8(b) and 9(b)). Figure 8(a) shows that at time = 0.1 the membrane potential becomes smoother at the node as the axon's diameter increases. Also, the long-range effects increase as the diamater increases: the potential is non-zero almost everywhere except in a very small neighbourhood of the middle of the internodal region. By comparison, figure 8(b) shows that in the classic case ( = 1) the potential decays to zero away from the node regardless of the size of the axon's diameter. This suggests a local behavior. The decay becomes slower as the diameter increases. Figure 9 shows that at fixed location = 700 , the amplitude of oscillations increases with the axon's diamater in both cases ( = 0.65 and = 1 ) which is in agreement with experimental observations. While the shapes of the oscillations are the same for the two values of , the amplitude is slightly higher for = 0.65 ( figure 9(a)).       Lastly, figure 10 shows that at time = 0.1 the shape of the membrane potential changes with the length of the internodal region for both = 0.65 and = 1, and, as expected, when = 0.65 the long range effects diminish as the internodal length decreases.

5-Conclusion
In this paper we assumed that spatial non-locality due to processes like long-range interactions of ions and water molecules at nanoscales and the anomalous diffusion of ions through the myelin sheath and the ECS affects the propagation of action potentials in myelinated neurons. We proposed a non-local cable equation to model spatial nonlocality of the action potentials propagation. We used spatial Caputo fractional derivatives and their properties to derive this equation. Further we used a numerical method to find the spatio-temporal distribution of the membrane potential in a leaky internodal region with one isopotential node described by a modified Hodgkin-Huxley model. Numerical results were shown for this dynamic problem. In our future work we plan to expand our model to neurons that also have backpropagating action potentials by using both left-and right-sided Caputo fractional spatial derivatives, and explore effects on the propagation of action potentials by introducing a link between a time-dependent non-local parameter ( ) and the diffusion of the extracellular potassium.