The Impact of Anomalous Diffusion on Action Potentials in Myelinated Neurons

: Action potentials in myelinated neurons happen only at specialized locations of the axons known as the nodes of Ranvier. The shapes, timings, and propagation speeds of these action potentials are controlled by biochemical interactions among neurons, glial cells, and the extracellular space. The complexity of brain structure and processes suggests that anomalous diffusion could affect the propagation of action potentials. In this paper, a spatio-temporal fractional cable equation for action potentials propagation in myelinated neurons is proposed. The impact of the ionic anomalous diffusion on the distribution of the membrane potential is investigated using numerical simulations. The results show spatially narrower action potentials at the nodes of Ranvier when using spatial derivatives of the fractional order only and delayed or lack of action potentials when adding a temporal derivative of the fractional order. These ﬁndings could reveal the pathological patterns of brain diseases such as epilepsy, multiple sclerosis, and Alzheimer’s disease, which have become more prevalent in the latest years. the in the shapes of the the of caused by the combined with could the presence of conditions. Fu- ture work on a


Introduction
Myelination of neurons is a critical growth process taking place during the first two years of life, which facilitates the fast transmission of information throughout the body. The process involves the wrapping of the neuronal axons by myelin sheaths made of multiple layers of glial plasma membrane. The myelinated regions are separated by the nodes of Ranvier ( Figure 1). Some specialized glial cells called oligodendrocytes produce myelin and create the neuro-glial-controlled clustering of voltage-gated sodium channels at the nodes of Ranvier and the aggregation of fast voltage-gated potassium channels in the myelin sheath [1,2]. The existence of some slow potassium channels in the nodes of Ranvier was reported in [3]. Experiments performed in cell cultures and in vivo reveal that nodal-like groups are created just before the myelination starts [4]. This indicates that myelination is highly dependent on timing and neuronal surroundings, especially the amounts of certain chemicals and water. An impediment to the myelination process or damage to the myelin sheath, which may or may not lead to demyelination, could cause serious, possibly life-threatening, neurologic disorders. For instance, multiple sclerosis is the most common demyelinating disorder affecting young adults that currently has no cure [5]. Mathematical models can provide essential insights into myelination/demyelination processes that can further inspire better diagnostic and therapeutic procedures for various neurological disorders.
Neurons transmit information throughout the body via action potentials. In a myelinated neuron an action potential happens when sodium flows inside the neuron at the node of Ranvier, and potassium is released into the extracellular space between the axon and the myelin (the periaxonal space). To prevent 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 astrocyte endfeet that form gap junctions with the outermost myelin layer [6]. Thus, in
3. [27,30,31]: 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 It was observed experimentally that the voltage-gated potassium channel Kv2.1 in the plasma membrane of a cultured human embryotic kidney cell displays anomalous diffusion [7]. Because the Kv2.1 channel shows similar clustering patterns in the membranes of cultured human embryonic kidney cells and native neurons [7], it is possible that anomalous diffusion of potassium can happen in the myelin sheath of neurons as well. The crowdedness of the very narrow ion channels [8] could be the cause of the anomalous diffusion of ions through the membrane's channels. Anomalous diffusion of ions that could affect the propagation of action potentials does not happen only through the ion channels of the neuronal membrane. According to [4], clusters of cytoskeletal scaffolding proteins, cell-adhesion molecules, and parts of the extracellular matrix (ECM) can be found near the nodes of Ranvier, which could influence ion flow. Furthermore, anomalous diffusion of ions through the extracellular space (ECS) is also possible due to (1) the ECS geometry; (2) dead spaces of the ECS, voids, or glial wrapping around cells; (3) flow obstruction by the ECM; (4) binding sites on the cell membranes of the ECM; (5) the presence of fixed negative charges on the ECM [9] that control ion mobility [10]; and (6) the tight control of the ion and water movement in the ECS by astrocytes [6,11]. In particular, anomalous diffusion might be involved in the regulation of the extracellular potassium concentration by astrocytes, which is essential for the proper creation and propagation of action potentials (for example, an unregulated elevated extracellular potassium concentration can cause an epileptic seizure [12]). Thus, anomalous diffusion of ions and water in the vicinity of the nodes of Ranvier is caused by a combination of structural geometry and highly dynamic and possibly emerging biochemical processes.
Anomalous diffusion through various materials has been successfully modeled using fractional calculus [13]. A few generalizations of the classic cable equation that models the spatio-temporal propagation of action potentials in neurons that involve fractional derivatives have been proposed in the literature. In [14,15], a fractional cable equation was proposed that uses a temporal Riemann-Liouville fractional derivative acting on the classic Laplace operator to model anomalous electro-diffusion in neurons. In [16], a space-fractional cable equation was derived to model the non-local forward (from the neuronal body to the axon terminals) propagation of action potentials in mature myelinated neurons. In this paper, a spatio-temporal fractional cable equation is proposed that involves both spatial and temporal fractional derivatives. The equation models the bidirectional (forward and backward) propagation of the action potentials in the presence of anomalous diffusion. The proposed work is a straightforward generalization of the mathematical model in [16]. In the case when the order of the temporal derivative is 1 and only the forward propagation of the action potentials is considered, the results presented in [16] are recovered. The proposed equation is coupled with the fractional Hodgkin-Huxley equations (a slight variation of those proposed in [17]) and solved numerically using the methods in [18,19]. Computer simulations show spatially narrower action potentials at the nodes of Ranvier at a fixed time when using spatial derivatives of the fractional order only and delayed or lack of action potentials when also using a temporal derivative of the fractional order. The spatially narrower shape of the action potentials resembles the experimentally observed action potentials [20] better than the shape predicted by the classic cable equation. Delayed or lack of action potentials bear the mark of neurological disorders [21][22][23][24][25][26]. It is important to notice that these signs of brain pathology were obtained without adding many more equations and parameters to the cable equation, which is the current mathematical modeling approach.
The structure of the paper is as follows. The mathematical preliminaries needed in the paper are given in Section 2. The derivation of the spatio-temporal fractional cable equation is presented in Section 3, together with a numerical solution to the equation coupled with the fractional Hodgkin-Huxley equations. Numerical simulations are shown in Section 4. The paper ends with a section of conclusions and further work.

Mathematical Model
In this section, the mathematical derivation of the spatio-temporal fractional cable equation is presented. The bidirectional propagation of action potentials between the cell body and the axon terminals is modeled using left-and right-sided spatial fractional order operators. For simplicity, in what follows the notations D α a+ , D α b− will be used for the spatial Caputo fractional derivatives and the notation ∂ α a+ will be used for the temporal Caputo fractional derivative.
For now, the time-dependent, spatial, non-local effects are neglected. Because the axonal branches are long and narrow, the action potentials depend only on one spatial variable [33] and thus the problem is one-dimensional. Thus, the axon can be modeled as a circular cylinder of constant radius r and one characteristic length L, the length of the internodal region, and assume that In this case, the element of the path along the integration of → E (a generalization of a measure proposed in [34]): . (14) since αΓ(α) = Γ(α + 1). The electric field is assumed to be uniform along the axon and the standard convention of current direction is valid: the membrane and synaptic current are positive when they are outward, and the electrode currents are positive when they are inward. Then, by using Expression (14), Formula (12) becomes Comparing Formulas (15) and (6) gives Consider an equally spaced discretization of [0, R] with step size ∆x and the number of nodes chosen such that, when α = 1, the following formula reduces to the classic forward approximation of the first-order derivative. Then Formula (15) yields On the other hand, using Formula (11) and the fact that p + q = 1 give where the Caputo fractional derivatives are taken with respect to the spatial variable x. Thus, as ∆x → 0 , Formulas (16) and (17) yield The uniformity of the current density in every cross-sectional area of the cylindrical neuron gives [33] where r L is the intracellular resistance and I is the electric current. Combining Formulas (18) and (19) gives the expression of the longitudinal current: where the negative sign comes from the current sign convention. Replacing Formula (19) into Formula (16) also yields Ohm's law: for the generalized electric resistance: As in [17,33,35], the following currents are introduced: where I m and I e are the membrane and external currents, respectively, and i m and i e are their corresponding currents per unit area. The capacitor current is I c and where c m is the specific membrane capacitance, and T is a characteristic time. The use of a temporal Caputo fractional derivative in the expression of I c was justified in [17] as giving a more accurate characterization of experimentally observed propagation of action potentials. The replacement of Formulas (20) and (22) in Kirchhoff's law in the element is shown in Figure 2b: gives where with initial conditions In Equations (35)- (38), is an externally applied current per unit area; , , and are the reverse potentials; and , , , , and are, respectively, the voltage-gated maximal conductances of and , and the leak conductances of , , and . The gating variables , , and ℎ are the activations of the and channels and the inactivation of the channel, respectively. Parameters , , , , , and were found by fitting the classic Hodgkin-Huxley model ( = 1 in Equations (35)- (38)) to experimental data (see, for instance, Section 2.2 of [38] and references within) and thus they should not be confused with the newly introduced non-local parameters and .

Numerical Discretization
A numerical discretization of the right-hand side of Equation (33) is given first. As ∆x → 0 : which replaced in Equation (25) and letting ∆x → 0 give the expression of the spatiotemporal fractional cable equation: If V(0, t) = V(R, t) = 0, then Equation (27) can be written in the equivalent form: For p = 1, q = 0 (or p = 0, q = 1) and α = β = 1, Equation (27) reduces to the classic cable equation [33]: A numerical discretization of the second-order spatial derivative in Equation (29) shows that in myelinated neurons the voltage at node n depends only on the voltages at the adjacent nodes of Ranvier, n − 1 and n + 1. However, a numerical discretization of the spatial Caputo fractional derivatives in Equation (27)  unwald-Letnikov formula can also be used to approximate numerically the temporal Caputo fractional derivative in Equation (27). (Although, in this paper, a different numerical scheme is applied to discretize the temporal Caputo fractional derivative (see later), it is somehow easier to mention the Gr .. unwald-Letnikov formula here to highlight temporal non-locality.) Thus, Equation (27) models the long-term memory and long-range interactions of the voltage in myelinated neurons.
Near the resting potential V rest the membrane current per unit area can be approximated as [33] where v = V − V rest is the voltage relative to the resting potential, and r m is the specific membrane resistance. Since V rest is constant, it follows from Formula (8) that Thus, Equation (27) becomes where The generalized electrotonic length λ has the same expression as in [16] and captures the effect of spatial non-locality via the parameter α ∈ (0, 1) on the relationship among the geometric parameters r, L and the resistances r m , r L of a neuron. Parameter β ∈ (0, 1) controls the time scale and accounts for long-term memory effects. Since Equation (27) remains valid for a time-varying α(t), the variations of α(t) and β, and the interplay between them, may reveal how aging and myelin assembly and remodeling affect the timing, speed, and shape of the action potentials.

Dynamic Problem
In this sub-section, a numerical solution to Equation (31) coupled with the fractional Hodgkin-Huxley equations is presented. The problem is stated and solved numerically in a leaky internodal region of length L with one isopotential node at x = L (Figure 2a).
The assumption that the length of the node of Ranvier can be neglected in the model is supported by the fact that in most neurons the length of the internodal region is of the order 1000 ÷ 2000 µm, which is much bigger than the usual nodal length of about 3 µm [36]. It is also assumed that the potential v(x, t) is zero at the middle of the internodal region and symmetric about the vertical line x = L/2. Thus, the spatial domain is [L/2, L]. The membrane potential of the internodal region is the solution to the following initial boundary value problem: v(x, 0) = 0, v(L/2, t) = 0, v(L, t) = f (t) (34) where f (t) = V(t) − V rest and V(t), the membrane potential of the node, is the solution of the following fractional Hodgkin-Huxley model, adapted from [37]: with initial conditions In Equations (35)- (38), I is an externally applied current per unit area; E Na , E K , and E Cl are the reverse potentials; and G Na , G K , G NaL , G KL , and G ClL are, respectively, the voltagegated maximal conductances of Na + and K + , and the leak conductances of Na + , K + , and Cl − . The gating variables m, n, and h are the activations of the Na + and K + channels and the inactivation of the Na + channel, respectively. Parameters α m , α n , α h , β m , β n , and β h were found by fitting the classic Hodgkin-Huxley model (β = 1 in Equations (35)-(38)) to experimental data (see, for instance, Section 2.2 of [38] and references within) and thus they should not be confused with the newly introduced non-local parameters α and β.

Numerical Discretization
A numerical discretization of the right-hand side of Equation (33) is given first. Let L/2 = x 0 < x 1 < · · · < x N−1 < x N = L be an equally spaced discretization of the interval [L/2, L] of constant step size ∆x. The spatial fractional derivatives are approximated using the right-and left-shifted Gr .. unwald-Letnikov formulas given in [18,39] and the boundary conditions in (34) (the link between the Caputo fractional derivatives and the Gr ..

unwald-Letnikov formulas is established through an integration by parts):
where Formula (41) can be re-written by changing the indexes of the terms in the sums: As in [18], Formula (43) introduces the following matrixL of dimension (N − 1) × (N − 1): Similarly, a matrixĽ of dimension (N − 1) × (N − 1) corresponding to Formula (44) can be introduced. However, using Formula (42), it is easy to show thatĽ =L T . The terms involving the boundary conditions in Formulas (43) and (44) are not included in the expressions of matricesL andĽ. The boundary terms are collected from the sums in Formulas (43) and (44) and stored in the following vectors: where the boundary conditions (34) have been used. Replacing Formulas (43)- (45) in Equation (33) and grouping the like terms give where and the last term of Formula (44) is stored in the following vector: The real part of (−1) α is used in the numerical simulations. Equations (46) and (35)-(38) form a system of coupled, non-linear, fractional ordinary differential equations that needs to be solved numerically using the initial conditions (34) and (39). Let Then, Equations (46) and (35)- (38) can be written in compact form: where F j (t), 1 ≤ j ≤ N + 3, are the corresponding right-hand sides of Equations (46) and (35)- (38). The initial conditions (34) and (39) are denoted by u j,0 , 1 ≤ j ≤ N + 3. A numerical discretization of system (47) is presented next. The fractional Hodgkin-Huxley equations in [17,35] were solved using either a nonstandard finite difference scheme [35] that combines the Gr .. unwald-Letnikov formula and a denominator function dependent on the step size satisfying certain conditions, or a semi-analytic approach [17] that involves the predator-corrector method given in [40]. The numerical scheme used in this paper is a variant of the one in [40]. More precisely, system (47) is solved numerically using MATLAB's function fde12 [19], which is an implementation of (1) a modified predictor-corrector method, which combines generalizations of the Adam-Bashforth (predictor) and Adam-Moulton (corrector) methods [41]; and (2) a fast Fourier transform (FFT) algorithm that reduces the computational cost [42]. For the sake of completeness, the main steps of the modified predictor-corrector method in [40,41] are presented further.
Applying the left-sided Riemann-Liouville fractional integral of order β ∈ (0, 1) to system (47) and using Formulas (9) and (6) give Let 0 = t 0 < t 1 < . . . < t M be an equally spaced discretization of a time interval [0, t] with constant step size ∆t, denoted by u j (t n ) = u n j . The generalization of the Adams-Bashforth method involves the use of the product rectangle rule to approximate: By replacing Formula (49) in Expression (48), the predator formula is obtained: The generalization of the Adams-Moulton method involves the use of the product trapezoidal quadrature formula to approximate: where a piecewise linear interpolant of F j is used on the left-hand side of Formula (51). Thus, where ϕ k,n+1 are the piecewise linear (hat) basis functions associated with nodes t k , 0 ≤ k ≤ n + 1. By now replacing Formula (51) in Expression (48), the corrector formula is obtained: Lastly, the sums in the predictor and corrector Formulas (50) and (52) can be written as discrete convolutions, which are efficiently calculated using the FFT algorithm given in [42].

Results
The parameters used in the numerical simulations are given in Table 1. The simulations use an internodal length of L = 1000 µm, which gives the maximum propagation speed of an action potential [43]. This value of L is within the range of values for the internodal length found in the literature ( [44] gives L = 200 ÷ 2000 µm where the lower (higher) value corresponds to the initial myelination (full maturity) stage [45]). In the simulations, the time step size is ∆t = 0.0001 and the space step size ∆x = 0.01. For simplicity, T = 1. The plots showing the spatial variations of v(x, t) (Figures 3a, 4a, 6a and 8a) represent the solutions v of problems (33)-(34) on the interval [L/2, L] and their mirror images about the vertical line at the node of Ranvier. This approach offers a better visualization of the action potentials at the node of Ranvier for a fixed time. Figures 3 and 4 show spatial and temporal variations of v(x, t) for β = 1 and various values of α when p = 1, q = 0 ( Figure 3) and p = 0, q = 1 ( Figure 4). As expected, when α = 1, there is no difference between the two cases since the classic solution of the cable equation is recovered. At fixed time t = 0.1 ms, the shape of the action potentials at the node of Ranvier becomes narrower as α decreases, with the narrower shapes obtained for the case p = 0, q = 1. (Figures 3a and 4a). These spatially narrower action potentials at the nodes look like the ones in ( [20], Figure 7D). In the case p = 0, q = 1, the numerical solution becomes unstable for α < 0.55. However, a solution in this scenario might not be physiologically meaningful since a too narrow action potential will lack the gradual spread of v within the internodal region seen in optical recordings of membrane potentials [20].
At fixed location x = 700 µm, the amplitude of v decreases and shifts to the right as α decreases (Figures 3b and 4b). The amplitude of v decreases faster with decreasing α in the case p = 0, q = 1 (Figure 4b) than in the case p = 1, q = 0 (Figure 3b). This effect of the bidirectional propagation on the temporal variations of the membrane potential within the internodal region has not been observed experimentally. Figure 5 shows the time variations of the gating variables n(t), m(t), and h(t) for β = 1 and various values of α when p = 1, q = 0. Since the gating variables are independent of the spatial variable x, the values of α, p, and q do not influence the variations of n(t), m(t), and h(t). Thus, the plots in Figure 5 are the solutions to the classic Hodgkin-Huxley equations.  Figures 3 and 4 show spatial and temporal variations of ( , ) for = 1 and various values of when = 1, = 0 ( Figure 3) and = 0, = 1 ( Figure 4). As expected, when = 1, there is no difference between the two cases since the classic solution of the cable equation is recovered. At fixed time = 0.1 ms, the shape of the action potentials at the node of Ranvier becomes narrower as decreases, with the narrower shapes obtained for the case = 0, = 1. (Figures 3a and 4a). These spatially narrower action potentials at the nodes look like the ones in ( [20], Figure 7D). In the case = 0, = 1, the numerical solution becomes unstable for < 0.55. However, a solution in this scenario might not be physiologically meaningful since a too narrow action potential will lack the gradual spread of within the internodal region seen in optical recordings of membrane potentials [20]. At fixed location = 700 μm, the amplitude of decreases and shifts to the right as decreases (Figures 3b and 4b). The amplitude of decreases faster with decreasing in the case = 0, = 1 (Figure 4b) than in the case = 1, = 0 (Figure 3b). This effect of the bidirectional propagation on the temporal variations of the membrane potential within the internodal region has not been observed experimentally. Figure 5 shows the time variations of the gating variables ( ), ( ), and ℎ( ) for = 1 and various values of when = 1, = 0. Since the gating variables are independent of the spatial variable , the values of , , and do not influence the variations of ( ), ( ), and ℎ( ). Thus, the plots in Figure 5 are the solutions to the classic Hodgkin-Huxley equations.       (Figure 6a), the amplitude of the potential v at fixed location x = 700 µm decreases as p decreases (q increases) (Figure 6b). The long-term memory represented by β = 0.7 < 1 causes delays in the opening and closing of the ion channels, which can be seen in the time variations of the gating variables n(t), m(t), and h(t) (Figure 7), which further produce the delay in the membrane potential seen in Figure 6b. Again, the values of α, p and q do not influence the variations in n(t), m(t), and h(t) (Figure 7).       Figures 6-10 suggest that the bidirectional propagation has little to no effect on the shapes of the action potentials at the nodes of Ranvier. At the fixed time t = 0.1 ms, the amplitude of the action potential at the node of Ranvier increases as the long-term memory increases (corresponding to β decreasing) and its shape becomes wider (Figure 8a), while at fixed location x = 700 µm the potentials flatten with decreasing β (Figure 8b). As the long-term memory increases, the delay in the opening and closing of the ion channels increases (Figure 9a-c) until they stop operating at β = 0.66 (Figure 9d). Figure 8b shows that the potential reaches a steady state value at β = 0.66. The spatio-temporal variations of v(x, t) in Figure 10 show how possible defects in the membrane potentials start to develop as β decreases. The increasing delay in the action potentials at the node of Ranvier with decreasing β (Figure 10b) looks like the shapes of the membrane potentials observed experimentally in transgenic Na v 1.8-null knockout DRG neurons ( [21], Figures 8b and 9c). Abnormal Na v 1.8 expression in some neurons may be linked to multiple sclerosis in humans [22]. Mutations of some sodium channels can cause persistent sodium currents that, combined with the blockade of potassium and calcium currents, lead to prolonged plateau potentials in the neurons of mice [25] (also [24], Figure 1B). These potentials resemble the ones in Figure 10b. According to [24], persistent sodium currents are linked to epilepsy. The same shape of the membrane potential shown in Figure 10b was also observed experimentally in neurons under spreading depolarization (SD) and a blockade of astrocytic glutamate transporters ( [26], Figure 8a-c). SD occurs in traumatic brain injuries and ischemic stroke. Lastly, cerebral anoxic depolarization that happen during a stroke or brain ischemia causes a potential like the one in Figure 10c [23]. Without immediate treatment, this condition leads to brain death. Therefore, there is no physiological justification to further decrease the value of β. Nevertheless, it was observed that the solution becomes unstable for β < 0.63.

Conclusions
In this paper, it was assumed that the anomalous diffusion of ions and water through the myelin sheath, nodes of Ranvier and ECS affect the propagation of action potentials in myelinated neurons. A spatio-temporal fractional cable equation is proposed to model In conclusion, the spatial non-locality controls the width of action potentials at the nodes of Ranvier, the narrower shapes being closer to the action potentials observed experimentally. The spatial parameters p and q control the strength of the forward and backward propagations of the action potentials; the weaker the forward propagation is, the lower the amplitude of the action potentials becomes. Lastly, the long-term memory in the presence of spatial non-locality could delay or stop the propagation of action potentials and widen the shapes of the action potentials at the nodes of Ranvier.

Conclusions
In this paper, it was assumed that the anomalous diffusion of ions and water through the myelin sheath, nodes of Ranvier and ECS affect the propagation of action potentials in myelinated neurons. A spatio-temporal fractional cable equation is proposed to model the effects of the long-range interactions and long-term memory of ions on the bidirectional propagation of action potentials in myelinated neurons. Caputo fractional derivatives and their properties were used to derive the equation. Numerical simulations showed the spatio-temporal distributions of the membrane potential in a leaky internodal region with one isopotential node described by the fractional Hodgkin-Huxley equations. The main finding is that the spatially wider shapes of the action potentials at the nodes of Ranvier and lack or delayed action potentials caused by the long-term memory combined with spatial non-locality could signal the presence of certain pathological conditions. Future work will explore possible effects on the propagation of action potentials of a time-dependent, non-local parameter α(t), and the interplay between α(t) and β.