A space-fractional cable equation for the propagation of action potentials in myelinated neurons
A non-local model of the propagation of action potentials in myelinated neurons
Corina S. Drapaca a , Sahin Ozdemir a , Elizabeth A. Proctor b,c a Department of Engineering Science and Mechanics, Pennsylvania State University, University Park PA 16802, USA b Departments of Neurosurgery and Pharmacology, Pennsylvania State College of Medicine, Hershey PA 17033, USA c Department of Biomedical Engineering, Pennsylvania State University, University Park PA 16802, USA
Abstract
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.
Keywords:
Non-locality; fractional calculus; action potentials; Hodgkin-Huxley model.
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 long-range 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 voltage-gated 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. Figure 1. (a) Schematic representation of the interactions among myelinated neurons (N), astrocytes (A) and myelinating oligodendrocytes (O). Astrocytes have numerous specialized connections (endfeet) to other astrocytes, the soma, dendrites and axons of neurons, the soma and processes of oligodendrocytes, and the brain fluid barriers. Neurons and astrocytes control the production and assembly of the myelin sheath around the axons by the oligodendrocytes. One oligodendrocyte can myelinate many axonal segments belonging to different neurons. Interactions among astrocytes and oligodendrocytes are controlled by gap junctions, while neuronal interactions happen through electrical and chemical synapses. (adapted from [16]) (b) Schematic representation of the diffusion barriers of the ECS: geometry, dead space microdomains, obstruction by the extracellular matrix, biding sites, and fixed negative charges on the extracellular matrix (adapted from [12]). * 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. A O O N N Axon
Terminals
Dendrites
Myelinated
Axon
Astrocyte
Endfeet - - Shape
Fixed Charge
Binding
Obstruction
Dead Space
ECS (a) (b)
3 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 one-dimensional 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.
Figure 2. Schematic of a myelinated neuron: the information travels from the cell body to the terminal through the axon made of myelinated regions separated by the nodes of Ranvier. The information is transmitted through action potentials which happen only at the nodes of Ranvier and are represented by vertical double arrows. Here we consider neurons which do not have backpropagating action potentials, and thus the only direction of propagation of the action potentials is from the cell body to the axon terminals.
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-34].
Definition 2.1:
1. Let 𝑓 ∈ 𝐿 (𝑎, 𝑏) and 𝛼 > 0 . The left-sided Riemann-Liouville fractional integral of order 𝛼 is: 𝐼 𝑎+𝛼 𝑓(𝑥) = 1Γ(𝛼) ∫ 𝑓(𝑦)(𝑥 − 𝑦) 𝑑𝑦 𝑥𝑎 (1) Cell
Body
Axon
Terminals
Nodes of Ranvier
Myelin
Dendrites
Direction of Propagation of Action Potentials
Axon and the right-sided Riemann-Liouville fractional integral of order 𝛼 is: 𝐼 𝑏−𝛼 𝑓(𝑥) = 1Γ(𝛼) ∫ 𝑓(𝑦)(𝑦 − 𝑥) 𝑑𝑦 𝑏𝑥 (2) where Γ(𝑠) = ∫ 𝑡 𝑠−1 𝑒 −𝑡 𝑑𝑡 ∞0 is the gamma function. In particular, 𝐼 𝑎+0 𝑓(𝑥) = 𝐼 𝑏−0 𝑓(𝑥) = 𝑓(𝑥).
2. If 𝑓: [𝑎, 𝑏] → 𝑹 is 𝑛 -times differentiable on (𝑎, 𝑏) ⊂ 𝑹 with 𝑓 (𝑛) ∈ 𝐿 (𝑎, 𝑏) and 𝑛 − 1 < 𝛼 < 𝑛, 𝑛 = 1,2,3 … , then the left-sided Caputo fractional derivative of order 𝛼 is: 𝐷 𝑎+𝛼 𝑓(𝑥) = 1Γ(𝑛 − 𝛼) ∫ 𝑓 (𝑛) (𝑦)(𝑥 − 𝑦) 𝛼+1−𝑛 𝑑𝑦 𝑥𝑎 = 𝐼 𝑎+𝑛−𝛼 𝑓 (𝑛) (𝑥) (3) and the right-sided Caputo fractional derivative of order 𝛼 is: 𝐷 𝑏−𝛼 𝑓(𝑥) = (−1) 𝑛 Γ(𝑛 − 𝛼) ∫ 𝑓 (𝑛) (𝑦)(𝑦 − 𝑥) 𝛼−𝑛+1 𝑑𝑦 𝑏𝑥 = (−1) 𝑛 𝐼 𝑏−𝑛−𝛼 𝑓 (𝑛) (𝑥) (4) where 𝑓 (𝑛) is the 𝑛 th -order derivative of 𝑓 . In particular, 𝐷 𝑎+𝑛 𝑓(𝑥) = 𝑓 (𝑛) (𝑥) , and 𝐷 𝑏−𝑛 𝑓(𝑥) = (−1) 𝑛 𝑓 (𝑛) (𝑥) . 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: 𝐷 𝑎+𝑛𝛼 𝑓(𝑥) = 𝐷 𝑎+𝛼 (𝐷 𝑎+(𝑛−1)𝛼 ) 𝑓(𝑥), 𝐷 𝑏−𝑛𝛼 𝑓(𝑥) = 𝐷 𝑏−𝛼 (𝐷 𝑏−(𝑛−1)𝛼 ) 𝑓(𝑥) (5) 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.
Proposition 2.1:
1. If 𝛼 > 0, 𝛾 > −1, and 𝑥 > 0 then: 𝐼 𝑥 𝛾 = Γ(𝛾 + 1)Γ(𝛾 + 𝛼 + 1) 𝑥 𝛾+𝛼 (6) 𝐷 𝑥 𝛾 = 𝛤(𝛾 + 1)𝛤(𝛾 − 𝛼 + 1) 𝑥 𝛾−𝛼 (7) 𝐷 𝑐 = 0 , c constant (8)
2. If 𝑛 − 1 < 𝛼 ≤ 𝑛 for 𝑛 = 1, 2, 3 … and the function 𝑓 has the property that there exist 𝑝 > 𝜇 ≥ −1 and a continuous function 𝑔 such that 𝑓 (𝑛) (𝑥) = 𝑥 𝑝 𝑔(𝑥) for 𝑥 > 0 , then: 𝐷 𝐼 𝑓(𝑥) = 𝑓(𝑥), 𝐼 𝐷 𝑓(𝑥) = 𝑓(𝑥) − ∑ 𝑓 (𝑗) (0+) 𝑥 𝑗 𝑗! 𝑛−1𝑗=0 (9)
3. If and the function 𝑓 has continuous sequential Caputo derivatives 𝐷 𝑘𝛼 𝑓 for 𝑘 = 0, 1, … , 𝑛 + 1, then the following generalized Taylor’s formula holds for 𝑥 > 0 : 𝑓(𝑥) = ∑ 𝐷 𝑓(0)𝛤(𝑘𝛼 + 1) 𝑛𝑘=0 𝑥 𝑘𝛼 + 𝐷 𝑓(𝜉) 𝛤((𝑛 + 1)𝛼 + 1) 𝑥 (𝑛+1)𝛼 (10)
5 with
The following representation is also valid: 𝑓(𝑥) = 𝑓(𝑥 ) + 𝐷 𝑓(𝑥 )Γ(𝛼 + 1) (𝑥 𝛼 − 𝑥 ) + 𝑅 (11) where 𝑅 is a reminder term.
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.
Definition 3.1:
The non-local voltage of order 𝜶(𝑡) = (𝛼 (𝑡), 𝛼 (𝑡), 𝛼 (𝑡)), 𝛼 𝑖 (𝑡) > 0, 𝑖 = 1,2,3 between two points in 𝑹 of position vectors and 𝒙 = (𝑥 , 𝑥 , 𝑥 ) is: 𝑉(𝒙, 𝑡) = − 1𝐿 (𝑡)−1 Γ(𝛼 (𝑡)) ∫ (𝑥 − 𝑦 ) 𝛼 (𝑡)−1 𝐸 (𝒚, 𝑡)𝑑𝑦 − 1𝐿 (𝑡)−1 Γ(𝛼 (𝑡)) ∫ (𝑥 − 𝑦 ) 𝛼 (𝑡)−1 𝐸 (𝒚, 𝑡)𝑑𝑦 − 1𝐿 (𝑡)−1 Γ(𝛼 (𝑡)) ∫ (𝑥 − 𝑦 ) 𝛼 (𝑡)−1 𝐸 (𝒚, 𝑡)𝑑𝑦 (12) where 𝑬(𝒚, 𝑡) = (𝐸 (𝒚, 𝑡), 𝐸 (𝒚, 𝑡), 𝐸 (𝒚, 𝑡)) is the electric field vector, 𝑡 ≥ 0 is a non-dimensional time, and 𝐿 𝑖 , 𝑖 =1, 2, 3 are characteristic lengths. As in the classic approach, we take 𝑉(𝟎) = 0 for mathematical convenience. 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: 𝛼 (𝑡) = 𝛼 (𝑡) = 𝛼 (𝑡) = 𝛼 = 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡, 0 < 𝛼 ≤ 1 (13) In this case, the element of path along the integration of 𝐸⃗ (a generalization of a measure proposed in [35]): 𝑑𝒚 𝜶(𝑡) = ( 1Γ(𝛼 (𝑡)) 𝑦 (𝑡)−1 𝑑𝑦 , 1Γ(𝛼 (𝑡)) 𝑦 (𝑡)−1 𝑑𝑦 , 1Γ(𝛼 (𝑡)) 𝑦 (𝑡)−1 𝑑𝑦 ) reduces to: 𝑑𝒚 𝜶 = 𝒚 𝛼−1 Γ(𝛼) 𝑑𝒚 . We assume further that the electric field is 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 formula (12) becomes:
𝑉(𝑥, 𝑡) = − 𝐸𝐿 𝛼−1 ∫ 𝑑(𝑥 − 𝑦) 𝛼 = − 𝑥0 𝐸𝐿 𝛼−1 ∫ 𝑑 [(𝑥 − 𝑦) 𝛼 Γ(𝛼 + 1)] 𝑥0 = 𝐸𝐿 𝛼−1 𝑥 𝛼 Γ(𝛼 + 1) (14)
By comparing formulas (14) and (6) we notice that:
𝑉(𝑥, 𝑡) = 𝐼 ( 𝐸𝐿 𝛼−1 ) From formula (14) we further get:
Δ𝑉(𝑥, 𝑡) = 𝑉(𝑥 + Δ𝑥, 𝑡) − 𝑉(𝑥, 𝑡) = 𝐸𝐿 𝛼−1
Γ(𝛼 + 1) ((𝑥 + Δ𝑥) 𝛼 − 𝑥 𝛼 ) (15) On the other hand, formula (11) gives:
Δ𝑉(𝑥, 𝑡) = 𝑉(𝑥 + Δ𝑥, 𝑡) − 𝑉(𝑥, 𝑡) ≈ 𝐷
𝑉(𝑥, 𝑡)Γ(𝛼 + 1) ((𝑥 + Δ𝑥) 𝛼 − 𝑥 𝛼 ) (16) where the Caputo fractional derivative is taken with respect to the spatial variable 𝑥 . Thus, as Δ𝑥 → 0 , formulas (15) and (16) yield: 𝐷 𝑉(𝑥, 𝑡) = 𝐸𝐿 𝛼−1 (17)
The uniformity of the current density in every cross-sectional area of the cylindrical neuron gives [27]:
𝐸 = 𝑟 𝐿 𝐼𝜋𝑟 (18) where 𝑟 𝐿 is the intracellular resistance and 𝐼 is the electric current. By combining formulas (17) and (18) we get the expression of the longitudinal current: 𝐼 𝐿 = − 𝜋𝑟 𝐿 𝛼−1 𝑟 𝐿 𝐷 𝑉(𝑥, 𝑡) (19) 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: 𝑅 𝛼 = 𝑟 𝐿 𝜋𝑟 𝐿 𝛼−1 Γ(𝛼 + 1) ((𝑥 + Δ𝑥) 𝛼 − 𝑥 𝛼 ) (20) Further we denote by: 𝑟̃ 𝐿 = 𝑟 𝐿 𝐿 𝛼−1 As in the classic case [1], we introduce the following currents: 𝐼 𝑚 = 2𝜋𝑟 Δ𝑥 𝑖 𝑚 𝐼 𝑒 = 2𝜋𝑟 Δ𝑥 𝑖 𝑒 𝐼 𝑐 = 2𝜋𝑟 Δ𝑥 𝑐 𝑚 𝜕𝑉𝜕𝑡 (21) 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: 𝐼 𝑐 + 𝐼 𝑚 − 𝐼 𝐿|𝑙𝑒𝑓𝑡 + 𝐼
𝐿|𝑟𝑖𝑔ℎ𝑡 − 𝐼 𝑒 = 0
7 we obtain: 𝑐 𝑚 𝜕𝑉𝜕𝑡 = 𝑟2𝑟̃ 𝐿 𝐷 𝑉 |𝑟𝑖𝑔ℎ𝑡 − 𝐷 𝑉 |𝑙𝑒𝑓𝑡 Δ𝑥 + 𝑖 𝑒 − 𝑖 𝑚 (22) As Δ𝑥 → 0 : 𝐷 𝑉 |𝑟𝑖𝑔ℎ𝑡 − 𝐷 𝑉 |𝑙𝑒𝑓𝑡 Δ𝑥 ≈ 𝜕𝜕𝑥 (𝐷 𝑉) which replaced in equation (22) and letting Δ𝑥 → 0 give the expression of the non-local cable equation : 𝑐 𝑚 𝜕𝑉𝜕𝑡 = 𝑟2𝑟̃ 𝐿 𝜕𝜕𝑥 (𝐷 𝑉) + 𝑖 𝑒 − 𝑖 𝑚 (23) Equation (23) is a space-fractional reaction-diffusion equation. For 𝛼 = 1 , equation (23) reduces to the classic cable equation [27]: 𝑐 𝑚 𝜕𝑉𝜕𝑡 = 𝑟2𝑟 𝐿 𝜕 𝑉𝜕𝑥 + 𝑖 𝑒 − 𝑖 𝑚 (24) 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]: 𝑖 𝑚 = 𝑉 − 𝑉 𝑟𝑒𝑠𝑡 𝑟 𝑚 ≔ 𝑣𝑟 𝑚 where 𝑣 = 𝑉 − 𝑉 𝑟𝑒𝑠𝑡 is the voltage relative to the resting potential, and 𝑟 𝑚 is the specific membrane resistance. Since 𝑉 𝑟𝑒𝑠𝑡 is constant it follows from formula (8) that: 𝐷 𝑉 = 𝐷 𝑣 Also: 𝜕𝑉𝜕𝑡 = 𝜕𝑣𝜕𝑡
Thus, equation (23) becomes: 𝜏 𝑚 𝜕𝑣𝜕𝑡 = 𝜆 𝛼+1 𝜕𝜕𝑥 (𝐷 𝑣) + 𝑟 𝑚 𝑖 𝑒 − 𝑣 (25) where we denoted by: 𝜏 𝑚 = 𝑟 𝑚 𝑐 𝑚 , 𝜆 𝛼+1 = 𝑟𝑟 𝑚 𝐿 𝐿 𝛼−1 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-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). Figure 3. (a) Schematic of a part of a neuron made of an internodal (myelinated) region of length 𝑳 and a node of Ranvier. The radius of the neuron is denoted by 𝒓 . (b) The electric currents entering (the longitudinal current 𝑰 𝑳|𝒍𝒆𝒇𝒕 and external current 𝑰 𝒆 ) and exiting ( the longitudinal current 𝑰 𝑳|𝒓𝒊𝒈𝒉𝒕 , the membrane current 𝑰 𝒎 , and the capacitor current 𝑰 𝒄 ) a neuronal segment of length 𝚫𝒙 . 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: 𝜏 𝑚 𝜕𝑣𝜕𝑡 = 𝜆 𝛼+1 𝜕𝜕𝑥 (𝐷 𝑣) − 𝑣, 𝑥 ∈ (0, 𝐿) 𝑣(𝑥, 0) = 0, 𝑣(0, 𝑡) = 0, 𝑣(𝐿, 𝑡) = 𝑓(𝑡) (26) where 𝑓(𝑡) = 𝑉(𝑡) − 𝑉 𝑟𝑒𝑠𝑡 and 𝑉(𝑡) , the membrane potential of the node, is solution of the modified Hodgkin-Huxley model [40] reproduced here for completeness. L x Node of
Ranvier
Myelin r Δ x 𝑰 𝑳 | 𝒍𝒆𝒇𝒕 𝑰 𝑳 | 𝒓𝒊𝒈𝒉𝒕 𝑰 𝒆 𝑰 𝒎 𝑰 𝒄 (a) (b) 𝑐 𝑚 𝑑𝑉𝑑𝑡 = 𝐼 − (𝐺 𝑁𝑎 𝑚 ℎ + 𝐺 𝑁𝑎𝐿 )(𝑉 − 𝐸 𝑁𝑎 ) − (𝐺 𝐾 𝑛 + 𝐺 𝐾𝐿 )(𝑉 − 𝐸 𝐾 )− 𝐺 𝐶𝑙𝐿 (𝑉 − 𝐸 𝐶𝑙 ) 𝑑𝑚𝑑𝑡 = 𝛼 𝑚 (1 − 𝑚) − 𝛽 𝑚 𝑚 𝑑𝑛𝑑𝑡 = 𝛼 𝑛 (1 − 𝑛) − 𝛽 𝑛 𝑛 𝑑ℎ𝑑𝑡 = 𝛼 ℎ (1 − ℎ) − 𝛽 ℎ ℎ 𝑉(0) = 𝑉 𝑟𝑒𝑠𝑡 𝑚(0) = 𝛼 𝑚 (𝑉(0))𝛼 𝑚 (𝑉(0)) + 𝛽 𝑚 (𝑉(0)) 𝑛(0) = 𝛼 𝑛 (𝑉(0))𝛼 𝑛 (𝑉(0)) + 𝛽 𝑛 (𝑉(0)) ℎ(0) = 𝛼 ℎ (𝑉(0))𝛼 ℎ (𝑉(0)) + 𝛽 ℎ (𝑉(0)) 𝛼 𝑚 = 0.32(𝑉 + 54)1 − exp(−0.25(𝑉 + 54)) , 𝛽 𝑚 = 0.28(𝑉 + 27)exp(0.2(𝑉 + 27)) − 1 𝛼 𝑛 = 0.032(𝑉 + 52)1 − exp(−0.2(𝑉 + 52)) , 𝛽 𝑛 = 0.5 exp (− 𝑉 + 5740 ) 𝛼 ℎ = 0.128 exp (− 𝑉 + 5018 ) , 𝛽 ℎ = 41 + exp(−0.2(𝑉 + 27)) (27) 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: 𝑤(𝑥, 𝑡) = 𝑣(𝑥, 𝑡) − 𝑣̃(𝑥, 𝑡) (28) where 𝑣̃(𝑥, 𝑡) is solution to the following boundary value problem: 𝜕𝜕𝑥 (𝐷 𝑣̃) = 0 𝑣̃(0, 𝑡) = 0 = 𝑣(0, 𝑡) 𝑣̃(𝐿, 𝑡) = 𝑓(𝑡) = 𝑣(𝐿, 𝑡) (29) Then the following identities hold: 𝜕𝜕𝑥 (𝐷 𝑤) = 𝜕𝜕𝑥 (𝐷 𝑣), 𝑤(0, 𝑡) = 𝑤(𝐿, 𝑡) = 0 (30)
The solution to problem (29) can be found as follows. Integrating one the differential equation gives: 𝐷 𝑣̃ = 𝑐 = constant The above identity can be transformed further into: 0 𝑣̃(𝑥, 𝑡) − 𝑣̃(0+, 𝑡) = 𝐼 (𝐷 𝑣̃) = 𝐼 𝑐 = 𝑐𝑥 𝛼 Γ(𝛼 + 1) (31) by using formulas (9) and (6). Finally, using the boundary conditions of problem (29) in formula (31) we find the solution to problem (29): 𝑣̃(𝑥, 𝑡) = 𝑓(𝑡) (𝑥𝐿) 𝛼 (32) 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: 𝜏 𝑚 𝜕𝑤𝜕𝑡 = 𝜆 𝛼+1 𝜕𝜕𝑥 (𝐷 𝑤) − 𝑤 − (𝑓(𝑡) + 𝜏 𝑚 𝑑𝑓(𝑡)𝑑𝑡 ) (𝑥𝐿) 𝛼 𝑤(0, 𝑡) = 𝑤(𝐿, 𝑡) = 0 𝑤(𝑥, 0) = 0 (33) 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 < 𝑥 < ⋯ < 𝑥 𝑁−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: 𝜕𝑤(𝑥 𝑗 , 𝑡)𝜕𝑡 = 𝜆 𝛼+1 𝜏 𝑚 Δ𝑥 𝛼+1 ∑ 𝑏 𝑖𝑗 𝑤(𝑥 𝑗 , 𝑡) 𝑁𝑖=0 − 1𝜏 𝑚 𝑤(𝑥 𝑗 , 𝑡) − 1𝜏 𝑚 (𝑓(𝑡) + 𝜏 𝑚 𝑑𝑓(𝑡)𝑑𝑡 ) (𝑥 𝑗 𝐿 ) 𝛼 , (34) with 𝑏 𝑖𝑗 = { (−1) 𝑗−𝑖+1 Γ(𝛼 + 2)Γ(j − i + 2)Γ(α + 1 − j + i) ,0, otherwise 0 < 𝑗 < 𝑁, 𝑖 ≤ 𝑗 + 1 (35)
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: (𝛼 + 1) 𝜆 𝛼+1 𝜏 𝑚 Δ𝑥 𝛼+1 Δ𝑡 ≤ 1 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 for the internodal length. In all the numerical simulations we used a time step ∆𝑡 = 0.01. Table 1. List of parameters with corresponding values and physical units.
Parameters Values and Units Reference 𝑉 𝑟𝑒𝑠𝑡 -65 𝑚𝑉 [40] 𝐸 𝑁𝑎 𝑚𝑉 [40] 𝐸 𝐾 -88 𝑚𝑉 [40] 𝐸 𝐶𝑙 -61 𝑚𝑉 [40] 𝐺 𝑁𝑎 𝑚𝑆/𝑚𝑚 [40] 𝐺 𝐾 𝑚𝑆/𝑚𝑚 [40] 𝐺 𝑁𝑎𝐿 𝑚𝑆/𝑚𝑚 [40] 𝐺 𝐾𝐿 𝑚𝑆/𝑚𝑚 [40] 𝐺 𝐶𝑙𝐿 𝑚𝑆/𝑚𝑚 [40] 𝑐 𝑚 𝜇𝐹/𝑚𝑚 [40] 𝐼 𝜇𝐴/𝑚𝑚 𝑟 𝑚𝑚 [27] 𝑟 𝑚 𝑘Ω ∙ 𝑚𝑚 [27] 𝑟 𝐿 𝑘Ω ∙ 𝑚𝑚 [27] 𝐿 𝑚𝑚 [27] 𝛼 Figure 4 shows the well-known solution of problem (27), while figure 5 highligths the differences in the spatio-temporal 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)). 2 Figure 4. Time variations of (a) voltage
𝑽(𝒕) , and (b) gating variables 𝒏(𝒕), 𝒎(𝒕), and 𝒉(𝒕) , obtained by solving the modified Hodgkin-Huxley equations (27).
Figure 5. Spatio-temporal variations of voltage 𝒗(𝒙, 𝒕) for (a) 𝜶 = 𝟏 and (b) 𝜶 = 𝟎. 𝟔𝟓 . Figure 6. Profiles of 𝒗(𝒙, 𝒕) at (a) 𝒕 = 𝟎. 𝟏 𝒎𝒔 and (b) 𝒙 = 𝟓𝟎𝟎 𝝁𝒎 . The fixed time in part (a) is almost the time when
𝑽(𝒕) reaches its first pick (see figure 4(a)).
Figure 7. Profiles of 𝒗(𝒙, 𝒕) at (a) 𝒕 = 𝟎. 𝟏 𝒎𝒔 and (b) 𝒙 = 𝟕𝟎𝟎 𝝁𝒎 for various values of parameter 𝜶 . Figure 8. Profiles of 𝒗(𝒙, 𝒕) for (a) 𝜶 = 𝟎. 𝟔𝟓 and (b) 𝜶 = 𝟏 (bottom row) at 𝒕 = 𝟎. 𝟏 𝒎𝒔 and (b) 𝒙 = 𝟕𝟎𝟎 𝝁𝒎 for various values of the axon’s radius 𝒓 . Figure 9. Profiles of 𝒗(𝒙, 𝒕) for (a) 𝜶 = 𝟎. 𝟔𝟓 and (b) 𝜶 = 𝟏 at 𝒙 = 𝟕𝟎𝟎 𝝁𝒎 for various values of the axon’s radius 𝒓 . Figure 10. Profiles of 𝒗(𝒙, 𝒕) for 𝜶 = 𝟎. 𝟔𝟓 (black) and 𝜶 = 𝟏 (red) at 𝒕 = 𝟎. 𝟏 𝒎𝒔 and for various values of the internodal length 𝑳 . The shapes for the two values of 𝜶 are identical for shorter internodal lengths. 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 non-locality 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.
6- References [1] Kwapieri, J., and S. Drozdz. “Physical Approach to Complex Systems” Physical Reports 515 (2012): 115- 226. [2] Ellis, G.F.R. “Physics, Complexity and Causality” Nature 435 (2005): 743. [3] Tuckwell, H.C. “Introduction to Theoretical Neurobiology: Nonlinear and Stochastic Theories” vol. 2. London, UK: Cambridge University Press (1988). [4] Dinariev, O.Y. “Equivalence between Classical Statistical Mechanics and Nonlocal Hydrodynamics for a Certain Class of External Forces” Russian Physics Journal 41(3) (1998): 211-216. [5] Dinariev, O.Y. “Dynamic Theory of Thermal Fluctuations with Allowance for the Spatiotemporal Nonlocality” Russian Physics Journal 43(4) (2000): 279-282. [6] Eisenberg, B. “Crowded Charges in Ion Channels” Advances in Chemical Physics S.A. Rice, Ed., John Wiley & Sons, Inc. 148 (2011): 77-223. [7] Liu, J.L., D. Xie, and B. Eisenberg. “Poisson-Fermi Formulation of Nonlocal Electrostatics in Electrolyte Solutions” Mol. Based Math. Biol. 5 (2017): 116-124. [8] Berthoumieux, H., and F. Paillusson. “Dielectric Response in the Vicinity of an Ion: A Nonlocal and Nonlinear Model of the Dielectric Properties of Water” J. Chem. Phys. 150 (2019): 094507. [9] Bakker, H.J. “Structural Dynamics of Aqueous Salt Solutions” Chem.Rev. 108 (2008): 1456-1473. [10] Ai, B., Z. Shao, and W. Zhong. “Rectified Brownian Transport in Corrugated Channels: Fractional Brownian Motion and Lévy Flights” J. Chem. Phys. 137 (2012): 174101. [11] Nicholson, C., P. Kamali-Zare, and L. Tao. “Brain Extracellular Space as a Diffusion Barrier” Comput.Vis.Sci. 14(7) (2011): 309-325. [12] Nicholson, C., and S. Hrab ě tov á . “Brain Extracellular Space: The Final Frontier in Neuroscience” Biophysical Journal 113 (2017): 2133-2142. [13] Morawski, M., T. Reinert, W. Meyer-Klaucke, F.E. Wagner, W. Troger, A. Reinert, C. Jager, G. Bruckner, and T. Arendt. “Ion Exchanger in the Brain: Quantitative Analysis of Perineuronally Fixed Anionic Binding Sites Suggests Diffusion Barriers with Ion Sorting Properties” Scientific Reports 5 (2015): 16471. [14] Min, R., and M.S. van der Knaap. “Genetic Defects Disrupting Glial Ion and Water Homeostasis in the Brain” Brain Pathology 28 (2018): 372-387. [15] Simard, M., and M. Nedergaard. “The Neurobiology of Glia in the Context of Water and Ion Homeostasis” Neuroscience 129(4) (2004): 877-896. [16] Baumann, N., and D. Pham-Dinh. “Biology of Oligodendrocyte and Myelin in the Mammalian Central Nervous System” Physiological Reviews 81(2) (2001): 871-927. [17] Bradl, M., and H. Lassmann. “Oligodendrocytes: Biology and Pathology” Acta Neuropathol. 119(1) (2010): 37-53. [18] Brazhe, A.R., G.V. Maksimov, E. Mosekilde, and O.V. Sosnovtseva. “Excitation Block in a Nerve Fibre Model Owing to Potassium-Dependent Changes in Myelin Resistance” Interface Focus 1 (2011): 86- 100. [19] Freeman, S.A., A. Desmazieres, D. Fricker, C. Lubetzki, and N. Sol-Foulon. “Mechanisms of Sodium Channel Clustering and its Influence on Axonal Impulse Conduction” Cell. Lol. Life Sci. 73 (2016): 723-735. [20] Weigel, A.V., B. Simon, M.M. Tamkun, and D. Krapf. “Ergodic and Nonergodic Processes Coexist in the Plasma Membrane as Observed by Single-Molecule Tracking” PNAS 108(16) (2011): 6438-6443. [21] West, B.J. “Fractional Calculus View of Complexity Tomorrow’s Science” Boca Raton, FL: CRC Press (2016). [22] West, B.J. “Nature’s Patterns and the Fractional Calculus” Series Fractional Calculus in Applied Sciences and Engineering 2. Berlin, Germany: De Gruyter (2017). [23] Metzler, R., and J. Klafter. “The Restaurant at the End of the Random Walk: Recent Developments in the Description of Anomalous Transport by Fractional Dynamics” J. Phys. A: Math. Gen. 37 (2004): R161-R208. [24] Stuart, G., N. Spruston, B. Sakmann, M. Hausser. “Action Potential Initiation and Backpropagation in Neurons of the Mammalian CNS” TINS 20(3) (1997): 125-131. [25] Kjellander, R. “Oscillatory and Long-Range Monotonic Exponential Decays of Electrostatic Interactions in Ionic Liquids and Other Electrolytes: The Significance of Dielectric Permittivity and Renormalized Charges” The Journal of Chemical Physics 148 (2018): 193701. [26] Bochud, T., and D. Challet. “Optimal Approximations of Power Laws with Exponentials: Application to Volatility Models with Long Memory” Quantitative Finance 7(6) (2007): 585-589. [27] Dayan, P, and L.F. Abbott. “Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems” Cambridge, MA: MIT Press (2001). [28] Henry, B.I, T.A.M. Langlands, and S.L. Wearne. “Fractional Cable Models for Spiny Neuronal Dendrites” Physical Review Letters 100 (2008): 128103. [29] Langlands, T.A.M., B.I. Henry, and S.L. Wearne. “Fractional Cable Equation Models for Anomalous Electrodiffusion in Nerve Cells: Infinite Domain Solutions” Journal of Mathematical Biology 59(6) (2009): 761-808. [30] Sambandham. B., and A.S. Vatsala. “Basic Results for Sequential Caputo Fractional Differential Equations” Mathematics 3 (2015): 76-91. [31] Samko, S., A.A. Kilbas, and O.I. Marichev. “Fractional Integrals and Derivatives: Theory and Applications” London, UK: Gordon and Breach Science Publishers (2000).17