Optic Nerve Microcirculation: Fluid Flow and Electrodiffusion
11 Optic Nerve Microcirculation: Fluid Flow and Electro-Diffusion
Yi Zhu , Shixin Xu
2, * , Robert S. Eisenberg , and Huaxiong Huang [1] Department of Mathematics and Statistics, York University, Toronto, ON, M3J 1P3, Canada [2] Duke Kunshan University, 8 Duke Ave., Kunshan, Jiangsu, China [3] Department of Applied Mathematics, Illinois Institute of Technology, Chicago IL 60616 USA [4] Department of Physiology & Biophysics, Rush University Chicago IL 60612, USA [5] Research Centre for Mathematics, Advanced Institute of Natural Sciences, Beijing Normal University (Zhuhai), China [6] BNU- HKBU United International College, Zhuhai, China
ABSTRACT:
Flows of water and various organic and inorganic molecules in the central nervous system are important in a wide range of biological and medical processes, as has recently become apparent (Nicholson et al. 2017 [1]). However, the exact mechanisms that drive these flows are often not known. Here we investigate the forces and flows in a tridomain model of the central nervous system. We use the experimental system of the optic nerve, investigated by the Harvard group Orkand et al [3, 20] as a protype of the central nervous system in general. We construct a model and choose its parameters to fit the experimental data. Asymptotic analysis and numerical computation show the significant role of water in convective ion transport. The full model (including water) and the electro-diffusion model (excluding water) are compared in detail to show the main physiological consequences of the complex structure as viewed in our model. In the full model, convection due to water flow dominates inside the glial domain. This water flow in the glia contributes significantly to spatial buffering of potassium in the extracellular space. Convection in the extracellular domain does not contribute significantly to spatial buffering. Electrodiffusion is the dominant mechanism for flows confined to the extracellular domain. Introduction
Recent experimental studies [2] suggest that transport in the central nervous system during sleep plays a critical role in maintaining the health of brain tissue. Since the nervous system is densely packed with neurons communicating with each other, question arises: how is homeostasis maintained. A few action potentials are known to significantly alter ion concentration in the immediate vicinity of peripheral and optic nerve cells [3, 4] and that change in concentration acts on more than one axon, producing “ cross talk ” . The question is then how does the central nervous system deal with changes in ion concentration produced by hundreds or thousands of action potentials and maintain a healthy environment? So how does the central nervous system maintain concentrations in its narrow extracellular space? What are the roles played by of glial cells and extracellular space? Here we examine a question that occurs to a classical physiologist familiar with the circulatory system of mammals: concentrations throughout mammals are maintained by convection [5-8]. The circulatory system of animals provides the convection that evolution uses to control ion concentrations. The convection is driven by a pump, the heart. But the extracellular space of the central nervous system does not have a mechanical pump to drive that convection. If convection is important in the central nervous system, what pushes the flow? The central nervous system contains nerve fibers and glia, separated by a narrow extracellular space. We use three domains to describe the flow and diffusion of ions and water in the optic nerve bundle of the central nervous system, hoping to glimpse general properties by which the central nervous system controls the concentration of ions in such narrow confines. Accumulation of potassium ions in particular is of great interest because of its likely importance in phenomena slower than 10 msec involving many action potentials. Accumulation of potassium on these time scales is thought to be involved in diseases [9-13] that cripple a substantial number of people. The convective properties we investigate are clearly involved in glaucoma and other diseases of the eye and our model may be of use in extending and understanding computational models of glaucoma [14, 15]. We investigate the optic nerve bundle as a paradigm hoping to create a paradigm useful for the central nervous system in channel. The optic nerve bundle contains paired cranial nerve bundled with cell bodies in the retina. It reaches from the eye through the optic chiasma to the cortex. The optic nerve transfers visual information from the retina to the vision centers of the brain using digital (actually binary) electrical signals (action potentials). Classical physiologists, who discovered all or none (i.e., binary) nature of the nerve signal would no doubt be amused by the present usage of the word ‘digital’ (whi ch originally referred to the ten fingers of human hands) instead of binary. The designers of binary arithmetic used more or less universally in our ‘digital’ circuitry, would share those feelings, we suspect. The optic nerve is customarily separated into four main regions. These are well described in reference [16] Fig. 1 in Chapter 2 page 13 and in reference [17] Fig. 1 on page 2. The four regions are (1) intraocular nerve head, (2) intraorbital region, (3) intracanalicular and (4) intracranial [16, 18]. In this paper, we mainly focus on the intraorbital region, which occupies more than half of the optic nerve.
There are about one million optic nerve fibers in the optic nerve bundle. The ganglion cells that are the cell bodies of the axons are scattered on the retina and form into a bundle at the optic disc. The bundle passes through the mesh-like lamina cribrosa region into the intraorbital region. Like almost all nerve cells, optic nerve fibers are functionally isolated, nearly insulated one from another. Current flow down one axon cannot flow into the adjacent axon or glia [19, 20]. The ‘ephaptic communication’ of concern to pioneers in electrophysiology rare occurs.
Glial cells wrap the nerve fiber bundles producing a narrow cleft of extracellular space between nerve fiber and glia, not so different from the extracellular space between squid axon and Schwann cell studied in the classical work [4] (see its Fig. 1bc). Sometimes, central retinal blood vessels (CRV, arterioles in fact) are found in the center of the optic nerve bundle in the intraorbital region. Here we consider the case where the blood vessel is not present, as in the optic nerve of the mud puppy, the amphibian salamander Necturus used in the experiments of Orkand et al. [3, 20].
Glial cells are connected to each other through connexin proteins, called ` gap junctions’ in the older literature (because of their appearance as a black ‘gap’ between outer cell membranes in osmium-stained sections in the electron microscope . Glial cells form an electrical syncytium (as do so many other cells, e.g., epithelia, cardiac muscle, lens of the eye, liver, etc.) in which current flow in one cell spreads into another with little extra resistance. In syncytia like this, inorganic ions, and many organic molecules (typically less than 2 nanometer diameter) can diffuse from cell to cell with hardly any restriction and thus with mobility and ionic conductance similar to that in cytoplasm. Speaking crudely, only polymers like nucleic acids, polypeptides like proteins, and polysaccharides like glycogen are unable to move through a connexin [21-24]. Glial cells are thought to play an important role in accelerating K + clearance from the extracellular space [25, 26]. The membranes of glial cells [20, 27] have protein channels that connect their cytoplasm to the narrow extracellular space and are selective for potassium ions. These are mostly the K ir channel proteins [25, 28] that conduct more current when membrane potentials are made more negative (i.e., ’hyperpolarized’), and conduct less current when the membrane potential is made more positive (i.e., ’depolarized’) [29]. Such current voltage characteristics have been called inward rectifying for a long time [27], much to the confusion of present day scientists who no longer understand what is being rectified, as the original authors did (when they made crystal radios based on the rectification of current flowing from a crystals to a cat whisker). Note the connexin proteins do not produce a gap between cells. Quite the opposite. They produce a path between the interiors of cells. The gap is a visual reality only in cells prepared for microscopy in a particular way. The gap is a functional artifact.
Figure 1.1: Optic nerve structure. a: Longitudinal section of the optic nerve; b: cross section of the optic nerve.
The optic nerve bundles are surrounded by the meningeal sheath which consists of dura mater, arachnoid mater and pia mater, and cerebrospinal fluid (CSF) in the subarachnoid space (SAS) [18, 30]. Also see Fig. 1.1a. The pia matter and dura matter are thin deformable shells, with mechanical properties important in glaucoma [30-33]. Andrew et. al [34] and Killer et. al [31, 35] show that the dura mater contains lymphatic vessels that drain CSF out of SAS (also see [33, 36]. Pia matter forms a macroscopic semipermeable membrane made of many cells, not just one lipid bilayer [37]. Many layered epithelia have been characterized as "semipermeable membranes" in low resolution studies of epithelia for more than a century. Filipidis et. al. [38] have written a most helpful review that identifies analogous leptomeningeal structures important in the physiology of "like pleura [39-44], peritoneum [45-50], pericardium [38], fetal membranes [51, 52], and leptomeninges [53]," We imagine that a general tridomain model may help understand many of these tissues.
Various mathematical models have been proposed to model the optic nerve. Most of them focus on the optic nerve head [33, 54-56]. Band et. al [57] proposed a mechanical model for the whole nerve fiber bundle, where only the water flow in the axon is considered. A cable type model [58] focusing on electrical properties of optic nerve is presented in [59]. However, no model considers the glial cells, as far as we know, although they are a prominent feature of the nerve bundle occupying a substantial part, perhaps the majority of its volume. The interactions between neuronal cells and glial cells have been included in models of the important phenomenon of spreading depression [9-13] widely thought to be related to epilepsy and migraine. Some two-compartment models for potassium clearance (or spatial buffering) include interactions between neuron cells and extracellular space [60] or interactions between glial cells and extracellular space [61-63]. A tri-compartment model using ordinary differential equations (ODEs) was introduced by Sibille [28] to study the role of K ir4 channels. It shows that the flows play an important role in the central nerve system[25] via influx and efflux routes to help waste clearance, which has been called (with understandable enthusiasm if not hyperbole) a final frontier of neuroscience [1]. Some models, including flow but not electrodiffusion, were introduced to study the pressure effect on the flow direction [36, 57, 64, 65]. Mori [66] proposed a multidomain model for cortical spreading depression, where ionic electrodiffusion and osmosis between different compartment are considered. We extend those results and present some general conclusions based on the analysis of a specific set of experiments Orkand et al [3, 20] using a model of ionic electrodiffusion and osmosis that has proven quite helpful in dealing with the substantial literature on electrodiffusion of the lens of the eye [67]. The rest of the paper is organized as follows. In Section 2, we present the mathematical model in detail. To simplify the analysis and computation, we assume that the optical nerve is cylindrical with three compartments laying on top of each other: axon, glial and extracellular ones. The transport of ions is governed by time dependent convective-electro-diffusion-reaction equations and the flow of water follows the Darcy’s law. The three compartments are coupled via trans -membrane fluxes for three major ions, namely sodium, potassium and chloride, treated as reaction terms. Model calibration is discussed in Section 3 by matching extracellular potassium concentration accumulation after the optical nerve is stimulated by a train of electric current pulses. In Section 4, we present estimates using order of magnitude analysis of transport of ionic and water fluxes cross membranes. They provide useful insight into the mechanisms for potassium clearance. In Section 5, we present numerical simulations and investigate the role of water flow (convection) in ionic transport during and after stimulus of the optical nerve. Our analysis shows that convection is very important within the glia, but it is not important for flows in the extracellular space itself. Water flow in glia has an indirect but significant effect in clearing potassium from the narrow extracellular space. This may be an important role for glia wherever they are found in the central nervous system, and even in structures of the peripheral nervous system. In Section 6, we provide concluding remarks on the limitation of our study and directions for future research. Figure 2.1: Domain of axial symmetry model. The optic nerve Ω 𝑂𝑃 consist of axon compartment Ω 𝑎𝑥 , glial compartment Ω 𝑔𝑙 and extracellular space Ω 𝑒𝑥𝑂𝑃 . The subarachnoid space only has extracellular space Ω 𝑒𝑥𝑆𝐴𝑆 . Mathematical Model
The model is first proposed in [67]. Here in order to make this paper self-contained, we summarize the model. The model deals with two types of flow: the circulation of water (hydrodynamics) and the circulation of ions (electrodynamics) in the glial compartment Ω 𝑔𝑙 , axon compartment Ω 𝑎𝑥 and extracellular space Ω 𝑒𝑥 . The glial compartment and axon compartment are limited to the optic nerve bundle, while extracellular space exists both in the optic nerve bundle Ω 𝑒𝑥𝑂𝑃 and in the subarachnoid space Ω 𝑒𝑥𝑆𝐴𝑆 , (See Fig. 2.1) Ω 𝑂𝑃 = Ω 𝑎𝑥 ∪ Ω 𝑔𝑙 ∪ Ω 𝑒𝑥𝑂𝑃 , Ω 𝑆𝐴𝑆 = Ω 𝑒𝑥𝑆𝐴𝑆 . The model is mainly derived from laws of conservation of ions and water for flow through membranes between intracellular compartments and extracellular space [68]. We first introduce the following notations used in the paper, where 𝑖 = Na + , K + , Cl − for ion spices, 𝑙 = 𝑒𝑥, 𝑔𝑙, 𝑎𝑥 for extracellular space, glial compartment and axon compartment, and 𝑘 =𝑔𝑙, 𝑎𝑥 for glial or axon membrane in the optic nerve. The summary of notations is listed in Appendix A1. Water Circulation
If we use 𝜂 𝑙 to describe the volume fraction of 𝑙 domain (𝑙 = 𝑔𝑙, 𝑎𝑥, 𝑒𝑥) , then conservation of mass in each domain yields 𝜕𝜂 𝑔𝑙 𝜕𝑡 + ℳ 𝑔𝑙 𝐿 𝑔𝑙𝑚 𝑈 𝑔𝑙𝑚 +▽⋅ (𝜂 𝑔𝑙 𝒖 𝑔𝑙 ) = 0, 𝑖𝑛 Ω 𝑂𝑃 , (2.1) 𝜕𝜂 𝑎𝑥 𝜕𝑡 + ℳ 𝑎𝑥 𝐿 𝑎𝑥𝑚 𝑈 𝑎𝑥𝑚 + 𝜕𝜕𝑧 (𝜂 𝑎𝑥 𝑢 𝑎𝑥𝑧 ) = 0, 𝑖𝑛 Ω 𝑂𝑃 , (2.2) ▽⋅ (𝜂 𝑔𝑙 𝒖 𝑔𝑙 ) +▽⋅ (𝜂 𝑒𝑥 𝒖 𝑒𝑥 ) + 𝜕𝜕𝑧 (𝜂 𝑎𝑥 𝑢 𝑎𝑥𝑧 ) = 0, 𝑖𝑛 Ω 𝑂𝑃 , (2.3) 𝜂 𝑔𝑙 + 𝜂 𝑎𝑥 + 𝜂 𝑒𝑥 = 1, 𝑖𝑛 Ω . (2.4) where the transmembrane water flux is proportional to the intracellular/extracellular hydrostatic pressure and osmotic pressure differences, i.e., Starling’s law on the membrane , 𝑈 𝑔𝑙𝑚 = 𝑝 𝑔𝑙 − 𝑝 𝑒𝑥 − 𝛾 𝑔𝑙 𝑘 𝐵 𝑇(𝑂 𝑔𝑙 − 𝑂 𝑒𝑥 ) , 𝑈 𝑎𝑥𝑚 = 𝑝 𝑎𝑥 − 𝑝 𝑒𝑥 − 𝛾 𝑎𝑥 𝑘 𝐵 𝑇(𝑂 𝑎𝑥 − 𝑂 𝑒𝑥 ) . We assume that glial cells are isotropic, and axons are anisotropic. Here 𝒖 𝑙 and 𝑝 𝑙 with 𝑙 =𝑔𝑙, 𝑎𝑥, 𝑒𝑥 are the velocity and pressure in the glial cells and axons and extracellular space, respectively. And 𝑘 𝐵 𝑇𝑂 𝑙 , is the osmotic [69, 70] defined by 𝑂 𝑒𝑥 = ∑ 𝑖 𝑐 𝑒𝑥𝑖 , 𝑂 𝑙 = ∑ 𝑖 𝑐 𝑙𝑖 + 𝐴 𝑙 𝜂 𝑗𝑟𝑒 𝜂 𝑙 , 𝑙 = 𝑔𝑙, 𝑎𝑥, where 𝐴 𝑙 𝜂 𝑙𝑟𝑒 𝜂 𝑙 > 0 ( 𝑙 = 𝑔𝑙, 𝑎𝑥) is the density of the permanent negatively charged protein in glial cell and axons that varies with the volume (fraction) of the region. The relation between the hydrostatic pressure 𝑝 𝑙 and volume fraction 𝜂 𝑙 (𝑙 = 𝑒𝑥, 𝑔𝑙, 𝑎𝑥) is connected by the force balance on the membrane 𝑘 (= 𝑔𝑙, 𝑎𝑥) [66, 69] 𝐾 𝑔𝑙 (𝜂 𝑔𝑙 − 𝜂 𝑔𝑙𝑟𝑒 ) = 𝑝 𝑔𝑙 − 𝑝 𝑒𝑥 − (𝑝 𝑔𝑙𝑟𝑒 − 𝑝 𝑒𝑥𝑟𝑒 ), 𝑖𝑛 Ω 𝑂𝑃 , (2.5) 𝐾 𝑎𝑥 (𝜂 𝑎𝑥 − 𝜂 𝑎𝑥𝑟𝑒 ) = 𝑝 𝑎𝑥 − 𝑝 𝑒𝑥 − (𝑝 𝑎𝑥𝑟𝑒 − 𝑝 𝑒𝑥𝑟𝑒 ), 𝑖𝑛 Ω 𝑂𝑃 , (2.6) where 𝐾 𝑔𝑙,𝑎𝑥 is the stiffness constant and 𝜂 𝑙𝑟𝑒 and 𝑝 𝑙𝑟𝑒 ( 𝑙 = 𝑒𝑥, 𝑔𝑙, 𝑎𝑥 ) are the resting state volume fraction and hydrostatic pressure. Fluid Velocity in the Glial Compartment.
As we mentioned before, the glial space is a connected space, where water can flow from cell to cell through connexin proteins joining membranes of neighboring cells. The velocity of fluid in glial syncytium 𝒖 𝑔𝑙 depends on the gradients of hydrostatic pressure and osmotic pressure: 𝑢 𝑔𝑙𝑟 = − 𝜅 𝑔𝑙 𝜏 𝑔𝑙 𝜇 ( 𝜕𝑝 𝑔𝑙 𝜕𝑟 − 𝛾 𝑔𝑙 𝑘 𝐵 𝑇 𝜕𝑂 𝑔𝑙 𝜕𝑟 ), (2.7) 𝑢 𝑔𝑙𝜃 = − 𝜅 𝑔𝑙 𝜏 𝑔𝑙 𝜇 (
1𝑟 𝜕𝑝 𝑔𝑙 𝜕𝜃 − 𝛾 𝑔𝑙 𝑘 𝐵 𝑇
1𝑟 𝜕𝑂 𝑔𝑙 𝜕𝜃 ), (2.8) 𝑢 𝑔𝑙𝑧 = − 𝜅 𝑔𝑙 𝜏 𝑔𝑙 𝜇 ( 𝜕𝑝 𝑔𝑙 𝜕𝑧 − 𝛾 𝑔𝑙 𝑘 𝐵 𝑇 𝜕𝑂 𝑔𝑙 𝜕𝑧 ). (2.9) The boundary conditions of fluid in the glial syncytium are as follows { 𝒖 𝑔𝑙 ⋅ 𝒓̂ = 0 , 𝑜𝑛 Γ ▽ 𝑝 𝑔𝑙 ⋅ 𝒛̂ = 0, 𝑜𝑛 Γ ▽ 𝑝 𝑔𝑙 ⋅ 𝒛̂ = 0, 𝑜𝑛 Γ 𝒖 𝑔𝑙 ⋅ 𝒓̂ = 0, 𝑜𝑛 Γ (2.10) Fluid Velocity in the Axon Compartment.
Since the axons are only connected in the longitudinal direction, the fluid velocity in axons region is defined along 𝑧 direction as 𝑢 𝑎𝑥𝑟 = 0, (2.11) 𝑢 𝑎𝑥𝜃 = 0, (2.12) 𝑢 𝑎𝑥𝑧 = − 𝜅 𝑎𝑥 𝜇 𝜕𝑝 𝑎𝑥 𝜕𝑧 , (2.13) Dirichlet boundary conditions are used to the fluid velocity in axons {𝒖 𝑎𝑥 ⋅ 𝒛̂ = 𝑢 , 𝑜𝑛 Γ 𝒖 𝑎𝑥 ⋅ 𝒛̂ = 𝑢 , 𝑜𝑛 Γ (2.14) where 𝑢 and 𝑢 are inflow and outflow velocities at Γ and Γ , respectively. Fluid Velocity in the Extracellular Space.
The extracellular space is narrow, and the extracellular velocity is determined by the gradients of hydrostatic pressure and electric potential 𝑢 𝑒𝑥𝑟 = − 𝜅 𝑒𝑥 𝜏 𝑒𝑥 𝜇 𝜕𝑝 𝑒𝑥 𝜕𝑟 − 𝑘 𝑒 𝜏 𝑒𝑥 𝜕𝜙 𝑒𝑥 𝜕𝑟 , (2.15) 𝑢 𝑒𝑥𝜃 = − 𝜅 𝑒𝑥 𝜏 𝑒𝑥 𝜇 1𝑟 𝜕𝑝 𝑒𝑥 𝜕𝜃 − 𝑘 𝑒 𝜏 𝑒𝑥 1𝑟 𝜕𝜙 𝑒𝑥 𝜕𝜃 , (2.16) 𝑢 𝑒𝑥𝑧 = − 𝜅 𝑒𝑥 𝜏 𝑒𝑥 𝜇 𝜕𝑝 𝑒𝑥 𝜕𝑧 − 𝑘 𝑒 𝜏 𝑒𝑥 𝜕𝜙 𝑒𝑥 𝜕𝑧 , (2.17) where 𝜙 𝑒𝑥 is the electric potential in the extracellular space, 𝜏 𝑒𝑥 is the tortuosity of extracellular region [68, 71] and 𝜇 is the viscosity of water, 𝑘 𝑒 is introduced to describe the effect of electro-osmotic flow [72-74], 𝜅 𝑒𝑥 is the permeability of extracellular space. Here the hydro permeability 𝜅 𝑒𝑥 , electric osmotic 𝑘 𝑒 and tortuosity 𝜏 𝑒𝑥 have two distinguished values in the region Ω 𝑒𝑥𝑂𝑃 and Ω 𝑒𝑥𝑆𝐴𝑆 . 𝜅 𝑒𝑥 = { 𝜅 𝑒𝑥𝑂𝑃 , 𝑖𝑛 Ω 𝑂𝑃 ,𝜅 𝑒𝑥𝑆𝐴𝑆 , 𝑖𝑛 Ω 𝑆𝐴𝑆 , 𝜏 𝑒𝑥 = { 𝜏 𝑒𝑥𝑂𝑃 , 𝑖𝑛 Ω 𝑂𝑃 ,𝜏 𝑒𝑥𝑆𝐴𝑆 , 𝑖𝑛 Ω 𝑆𝐴𝑆 , 𝑘 𝑒𝑥 = { 𝑘 𝑒𝑥𝑂𝑃 , 𝑖𝑛 Ω 𝑂𝑃 ,𝑘 𝑒𝑥𝑆𝐴𝑆 , 𝑖𝑛 Ω 𝑆𝐴𝑆 (2.18)
Since Γ ∪ Γ are the far end of optic nerve away from eyeball and next to connect to optic canal, we assume the hydrostatic pressure of extracellular is equal to the cerebrospinal fluid pressure. On the other hand, the intraocular pressure (IOP) is imposed at Γ where the extracellular space is connected to the retina. At boundary Γ , we assume a non-permeable boundary. We are aware of the significance of the pressures and flows at these boundaries for clinical phenomena including glaucoma [14, 32, 57, 75] and will return to that subject in later publications. The water flow across the semi-permeable membrane Γ is produced by the lymphatic drainage on the dura membrane, which depends on the difference between extracellular pressure and orbital pressure (OBP). We assume the velocity across the pia membrane Γ , is continuous and determined by the combination of hydrostatic and osmotic pressures. To summarize, the boundary conditions of the extracellular fluid are { 𝒖 𝑒𝑥 ⋅ 𝒓̂ = 0,𝑝 𝑒𝑥 = 𝑝 𝐶𝑆𝐹 ,𝒖 𝑒𝑥𝑆𝐴𝑆 ⋅ 𝒓̂ = 𝐿 𝑑𝑟𝑚 (𝑝 𝑒𝑥𝑆𝐴𝑆 − 𝑝 𝑂𝐵𝑃 ),𝒖 𝑒𝑥 ⋅ 𝒓̂ = 0,𝑝 𝑒𝑥 = 𝑝 𝐼𝐶𝑃 ,𝒖 𝑒𝑥𝑂𝑃 ⋅ 𝒓̂ = 𝒖 𝑒𝑥𝑆𝐴𝑆 ⋅ 𝒓̂ = 𝐿 𝑝𝑖𝑎𝑚 (𝑝 𝑒𝑥𝑂𝑃 − 𝑝 𝑒𝑥𝑆𝐴𝑆 − 𝛾 𝑝𝑖𝑎 𝑘 𝐵 𝑇(𝑂 𝑒𝑥𝑂𝑃 − 𝑂 𝑒𝑥𝑆𝐴𝑆 ) ), (2.19) where 𝑝 𝐶𝑆𝐹 is the cerebrospinal fluid pressure [57] and 𝑝 𝐼𝐶𝑃 is the pressure in the eye and 𝑝 𝑂𝐵𝑃 is the orbital pressure on the dura mater.
Ion Transport
The conservation of chemical species implies the following system of partial differential equations to describe the dynamics of ions in each region, for 𝑖 = Na + , K + , Cl − 𝜕(𝜂 𝑔𝑙 𝑐 𝑔𝑙𝑖 )𝜕𝑡 + ℳ 𝑔𝑙 𝐽 𝑔𝑙𝑚,𝑖 +▽⋅ (𝜂 𝑔𝑙 𝒋 𝑔𝑙𝑖 ) = 0, 𝑖𝑛 Ω 𝑂𝑃 , (2.20) 𝜕(𝜂 𝑎𝑥 𝑐 𝑎𝑥𝑖 )𝜕𝑡 + ℳ 𝑎𝑥 𝐽 𝑎𝑥𝑚,𝑖 + 𝜕𝜕𝑧 (𝜂 𝑎𝑥 𝑗 𝑎𝑥,𝑧𝑖 ) = 0, 𝑖𝑛 Ω 𝑂𝑃 , (2.21) 𝜕(𝜂 𝑒𝑥 𝑐 𝑒𝑥𝑖 )𝜕𝑡 − ℳ 𝑎𝑥 𝐽 𝑎𝑥𝑚,𝑖 − ℳ 𝑔𝑙 𝐽 𝑔𝑙𝑚,𝑖 +▽⋅ (𝜂 𝑒𝑥 𝒋 𝑒𝑥𝑖 ) = 0, 𝑖𝑛 Ω 𝑂𝑃 , (2.22) where the last equation reduces to the following in the Ω 𝑆𝐴𝑆 region, 𝜕𝑐 𝑒𝑥𝑖,𝑆𝐴𝑆 𝜕𝑡 +▽⋅ 𝒋 𝑒𝑥𝑖,𝑆𝐴𝑆 = 0. (2.23) The transmembrane ion flux 𝐽 𝑘𝑚,𝑖 (𝑘 = 𝑎𝑥, 𝑔𝑙) consists of active ion pump source 𝐽 𝑝,𝑘𝑖 and passive ion channel source 𝐽 𝑐,𝑘𝑖 , on the 𝑘 membrane, 𝐽 𝑘𝑚,𝑖 = 𝐽 𝑝,𝑘𝑖 + 𝐽 𝑐,𝑘𝑖 , 𝑘 = 𝑎𝑥, 𝑔𝑙, 𝑖 = Na + , K + , Cl − . Here 𝐽 𝑝,𝑘𝑖 is the active Na/K pump source and 𝐽 𝑐,𝑘𝑖 is passive source for ion 𝑖 (usually) through channels in the axons ( 𝑘 = 𝑎𝑥 ) or glial cells membranes ( 𝑘 = 𝑔𝑙 ). In the glial cell membranes. The 𝐽 𝑐,𝑔𝑙𝑖 is defined as 𝐽 𝑐,𝑔𝑙𝑖 = 𝑔 𝑔𝑙𝑖 𝑧 𝑖 𝑒 (𝜙 𝑔𝑙 − 𝜙 𝑒𝑥 − 𝐸 𝑔𝑙𝑖 ), 𝑖 = Na + , K + , Cl − . (2.24) where the Nernst potential is used to describe the gradient of chemical potential 𝐸 𝑔𝑙𝑖 = 𝑘 𝐵 𝑇𝑒𝑧 𝑖 log ( 𝑐 𝑒𝑥𝑖 𝑐 𝑔𝑙𝑖 ) and the conductance 𝑔 𝑔𝑙𝑖 for ion 𝑖 in the glial membrane is a fixed constant, independent of voltage and time. Note that in models like ours the concentrations are functions of time, unlike in the original Hodgkin Huxley model [76] . On the ax on’s membrane, 𝐽 𝑐,𝑎𝑥𝑖 is defined as 𝐽 𝑐,𝑎𝑥𝑖 = 𝑔 𝑎𝑥𝑖 𝑧 𝑖 𝑒 (𝜙 𝑎𝑥 − 𝜙 𝑒𝑥 − 𝐸 𝑎𝑥𝑖 ), 𝑖 = Na + , K + , Cl − , (2.25) where 𝑔 𝑎𝑥𝑁𝑎 = 𝑔̅ 𝑁𝑎 𝑚 ℎ + 𝑔 𝑙𝑒𝑎𝑘𝑁𝑎 , 𝑔 𝑎𝑥𝐾 = 𝑔̅ 𝐾 𝑛 + 𝑔 𝑙𝑒𝑎𝑘𝐾 , 𝑔 𝑎𝑥𝐶𝑙 = 𝑔 𝑙𝑒𝑎𝑘𝐶𝑙 . The time dependent dynamic of open probability, often loosely called ‘gating’ is governed by the
Hodgkin-Huxley model [76, 77]. 𝑑𝑛𝑑𝑡 = 𝛼 𝑛 (1 − 𝑛) − 𝛽 𝑛 𝑛, 𝑑𝑚𝑑𝑡 = 𝛼 𝑚 (1 − 𝑚) − 𝛽 𝑚 𝑚, 𝑑ℎ𝑑𝑡 = 𝛼 ℎ (1 − ℎ) − 𝛽 ℎ ℎ, (2.26) We assume that the only pump is the
Na/K active transporter. We are more than aware that other active transport systems can and likely do move ions and thus water in this system. They will be included as experimental information becomes available. In the case of the Na / K pump 𝐽 𝑝,𝑘𝑖 (𝑘 = 𝑎𝑥, 𝑔𝑙), the strength of the pump 𝐼 𝑘 depends on the concentration in the intracellular and extracellular space [77, 78], i.e. 𝐽 𝑝,𝑘𝑁𝑎 = 3 𝐼 𝑘 𝑒 , 𝐽 𝑝,𝑘𝐾 = −2 𝐼 𝑘 𝑒 , 𝐽 𝑝,𝑘𝐶𝑙 = 0, 𝑘 = 𝑎𝑥, 𝑔𝑙, (2.27) where 𝐼 𝑔𝑙 = 𝐼 𝑔𝑙,1 ( 𝑐 𝑔𝑙𝑁𝑎 𝑐 𝑔𝑙𝑁𝑎 +𝐾 𝑁𝑎1 ) ( 𝑐 𝑒𝑥𝐾 𝑐 𝑒𝑥𝐾 +𝐾 𝐾1 ) + 𝐼 𝑔𝑙,2 ( 𝑐 𝑔𝑙𝑁𝑎 𝑐 𝑔𝑙𝑁𝑎 +𝐾 𝑁𝑎2 ) ( 𝑐 𝑒𝑥𝐾 𝑐 𝑒𝑥𝐾 +𝐾 𝐾2 ) , (2.28) 𝐼 𝑎𝑥 = 𝐼 𝑎𝑥,1 ( 𝑐 𝑎𝑥𝑁𝑎 𝑐 𝑎𝑥𝑁𝑎 +𝐾 𝑁𝑎1 ) ( 𝑐 𝑒𝑥𝐾 𝑐 𝑒𝑥𝐾 +𝐾 𝐾1 ) + 𝐼 𝑎𝑥,2 ( 𝑐 𝑎𝑥𝑁𝑎 𝑐 𝑎𝑥𝑁𝑎 +𝐾 𝑁𝑎2 ) ( 𝑐 𝑒𝑥𝐾 𝑐 𝑒𝑥𝐾 +𝐾 𝐾2 ) , 𝐼 𝑘,1 and 𝐼 𝑘,2 are related to 𝛼 − and 𝛼 − isoform of Na/K pump on 𝑘 membrane. The definitions of ion flux in each domain are as follows, for 𝑖 = Na + , K + , Cl − , 0 𝒋 𝑙𝑖 = 𝑐 𝑙𝑖 𝒖 𝑙 − 𝐷 𝑔𝑙𝑖 𝜏 𝑙 (▽ 𝑐 𝑙𝑖 + 𝑧 𝑖 𝑒𝑘 𝐵 𝑇 𝑐 𝑙𝑖 ▽ 𝜙 𝑙 ) , 𝑙 = 𝑔𝑙, 𝑒𝑥, (2.29) 𝑗 𝑎𝑥,𝑧𝑖 = 𝑐 𝑎𝑥𝑖 𝑢 𝑎𝑥𝑧 − 𝐷 𝑎𝑥𝑖 ( 𝜕𝑐 𝑎𝑥𝑖 𝜕𝑧 + 𝑧 𝑖 𝑒𝑘 𝐵 𝑇 𝑐 𝑎𝑥𝑖 𝜕𝜙 𝑎𝑥 𝜕𝑧 ). (2.30) For the axon compartment and glial compartment boundary condition, we have 𝑐 𝑎𝑥𝑖 = 𝑐 𝑎𝑥𝑖,𝑟𝑒 , 𝑜𝑛 Γ ∪ Γ , (2.31) and { 𝒋 𝑔𝑙𝑖 ⋅ 𝒓̂ = 0, 𝑜𝑛 Γ , 𝑐 𝑔𝑙𝑖 = 𝑐 𝑔𝑙𝑖,𝑟𝑒 , 𝑜𝑛 Γ ∪ Γ ,𝒋 𝑔𝑙𝑖 ⋅ 𝒓̂ = 0, 𝑜𝑛 Γ , (2.32) where the Dirichlet boundary conditions are used at locations Γ ∪ Γ for axons and glial cell membranes, and a non-flux boundary condition is used for glial cells ions flux on pia mater Γ . For the extracellular space boundary condition, similar boundary conditions are imposed except on the pia mater Γ . The flux across the pia mater is assumed continuous and Ohm’s law is used [70]. Additionally, a non-permeable boundary condition is used at location Γ and a homogeneous Neumann boundary condition is applied at the location of the dura mater Γ , { 𝒋 𝑒𝑥𝑖 ⋅ 𝒓̂ = 0, 𝑜𝑛 Γ ,𝑐 𝑒𝑥𝑖 = 𝑐 𝑐𝑠𝑓𝑖 , 𝑜𝑛 Γ ∪ Γ ,▽ 𝑐 𝑒𝑥𝑖 ⋅ 𝒓̂ = 0, 𝑜𝑛 Γ ,𝒋 𝑒𝑥𝑖 ⋅ 𝒛̂ = 0, 𝑜𝑛 Γ ,𝑐 𝑒𝑥𝑖 = 𝑐 𝑒𝑦𝑒𝑖 , 𝑜𝑛 Γ ,𝒋 𝑒𝑥𝑖,𝑂𝑃 ⋅ 𝒓̂ = 𝒋 𝑒𝑥𝑖,𝑆𝐴𝑆 ⋅ 𝒓̂ = 𝐺 𝑝𝑖𝑎𝑖 𝑧 𝑖 𝑒 (𝜙 𝑒𝑥𝑂𝑃 − 𝜙 𝑒𝑥𝑆𝐴𝑆 − 𝐸 𝑝𝑖𝑎𝑖 ), 𝑜𝑛 Γ , (2.33) Multiplying equations in (2.20-2.22) with 𝑧 𝑖 𝑒 respectively, summing up, and using the charge neutrality condition, we have the following system for the electric fields in 𝑎𝑥, 𝑔𝑙, 𝑒𝑥 , ∑ 𝑖 𝑧 𝑖 𝑒ℳ 𝑔𝑙 (𝐽 𝑝,𝑔𝑙𝑖 + 𝐽 𝑐,𝑔𝑙𝑖 ) + ∑ 𝑖 𝑧 𝑖 𝑒 ▽⋅ (𝜂 𝑔𝑙 𝒋 𝑔𝑙𝑖 ) = 0, (2.34) ∑ 𝑖 𝑧 𝑖 𝑒ℳ 𝑎𝑥 (𝐽 𝑝,𝑎𝑥𝑖 + 𝐽 𝑐,𝑎𝑥𝑖 ) + ∑ 𝑖 𝑧 𝑖 𝑒 𝜕𝜕𝑧 (𝜂 𝑎𝑥 𝑗 𝑎𝑥,𝑧𝑖 ) = 0, (2.35) ∑ 𝑖 𝑧 𝑖 𝑒 ▽⋅ (𝜂 𝑔𝑙 𝒋 𝑔𝑙𝑖 ) + ∑ 𝑖 𝑧 𝑖 𝑒 𝜕𝜕𝑧 (𝜂 𝑎𝑥 𝑗 𝑎𝑥,𝑧𝑖 ) + ∑ 𝑖 𝑧 𝑖 𝑒 ▽⋅ (𝜂 𝑒𝑥 𝒋 𝑒𝑥𝑖 ) = 0. (2.36) In the subarachnoid space Ω 𝑆𝐴𝑆 , the extracellular equations reduce to ∑ 𝑖 𝑧 𝑖 𝑒 ▽⋅ (∑ 𝑖 𝒋 𝑒𝑥𝑖,𝑆𝐴𝑆 ) = 0. (2.37) The boundary conditions for electric fields 𝜙 𝑎𝑥 , 𝜙 𝑔𝑙 and 𝜙 𝑒𝑥 are given below. In the axon compartment: 1 {▽ 𝜙 𝑎𝑥 ⋅ 𝒛̂ = 0, 𝑜𝑛 Γ ▽ 𝜙 𝑎𝑥 ⋅ 𝒛̂ = 0, 𝑜𝑛 Γ . (2.38) In the glial compartment: { 𝜙 𝑔𝑙 ⋅ 𝒓̂ = 0, 𝑜𝑛 Γ ,▽ 𝜙 𝑔𝑙 ⋅ 𝒛̂ = 0, 𝑜𝑛 Γ ,▽ 𝜙 𝑔𝑙 ⋅ 𝒛̂ = 0, 𝑜𝑛 Γ ,∑ 𝑖 𝑧 𝑖 𝑒𝒋 𝑔𝑙𝑖 ⋅ 𝒓̂ = 0, 𝑜𝑛 Γ , (2.39) and in the extracellular space: { ∑ 𝑖 𝑧 𝑖 𝑒𝒋 𝑒𝑥𝑖 ⋅ 𝒓̂ = 0, 𝑜𝑛 Γ ,▽ 𝜙 𝑒𝑥 ⋅ 𝒛̂ = 0, 𝑜𝑛 Γ ∪ Γ ,▽ 𝜙 𝑒𝑥 ⋅ 𝒓̂ = 0, 𝑜𝑛 Γ ,∑ 𝑖 𝑧 𝑖 𝑒𝒋 𝑒𝑥𝑖 ⋅ 𝒛̂ = 0, 𝑜𝑛 Γ ,▽ 𝜙 𝑒𝑥 ⋅ 𝒛̂ = 0, 𝑜𝑛 Γ ,∑ 𝑖 𝑧 𝑖 𝑒𝒋 𝑒𝑥𝑖,𝑂𝑃 ⋅ 𝒓̂ = ∑ 𝑖 𝑧 𝑖 𝑒𝒋 𝑒𝑥𝑖,𝑆𝐴𝑆 ⋅ 𝒓̂ = ∑ 𝑖 𝐺 𝑝𝑖𝑎𝑖 (𝜙 𝑒𝑥𝑂𝑃 − 𝜙 𝑒𝑥𝑆𝐴𝑆 − 𝐸 𝑝𝑖𝑎𝑖 ), 𝑜𝑛 Γ . (2.40) In the rest of this paper, the full electric-diffusion-convection model is defined by equations (2.1) through (2.40). The electric-diffusion model is defined by equations (2.20)-(2.40). The electric diffusion model is a reduced version of the full model in which water is neglected.
3. Model Calibration
We calibrate our model with the Orkand et al experimental data [3] that measures the change in potential across the glial membrane produced by a train of action potentials. In the experiment, optic nerve has been put in bathing solutions with three different K + concentration (1.5 mM, 3 mM, 4.5 mM) and the resting potential across the glia membrane was measured. Then the axon was stimulated to give a train of action potentials. The action potentials increased K + in extracellular space (ECS). The accumulated K + then made the glia membrane potential more positive. In the simulation, we applied a train of stimuli with frequency for 1 s to the axon membrane at 𝑧 = 1.5 mm, 13.5 mm , 𝑎 = 48 μm . Each individual stimulus in the train lasted 3 ms (as Orkand’s paper indicated) and had strength 3 mA/m . The stimulus was large enough to exceed threshold and generate action potentials. We set the ECS K + to be 1.5 mM , 3 mM , or 4.5 mM and record the largest absolute value of the change in glial membrane potential in each case as in the Fig. 3.2. This number is loosely called ‘the depolarization’ in most laboratories.
The blue symbols show experimental data, red ones are the simulations results of electrodiffusion model and the green ones are the full model. Fig. 3.2 shows that both the full model and electrodiffusion model could match the experimental resting potentials (solid symbols) and depolarizations (open symbols) very well for the different ECS K + concentrations. 2 Figure 3.1: a: axon membrane potential profile when eye-end axon stimulated. b: axon membrane potential profile when two-end axon simulated.
Fig. 3.1 shows the propagation of the axon action potential. The membrane potential from axons at the center of the optic nerve bundle is shown when different locations of the axon had been stimulated. In both eye-end and two-end cases, the stimulus current was applied at 𝑡 = 1 ms . In Fig. 3.1a, the stimulus was applied near to the optic nerve near the eye-end ( 𝑧 = 1.5 mm ). At 𝑡 =10 ms , the action potential completely has propagated and left the location near far-eye-end ( ). The axon in the optic nerve of the mud puppy is unmyelinated. This speed of action potential propagation in the model lies in the range of the action potential speeds typical of unmyelinated axons, i.e., between and [85]. In the Fig. 3.1b, when the two-ends of the axon stimulated, the axon membrane potential has is more uniform spatially at each time point in compare to the single side stimulus case. Presumably, Orkand et al used the dual stimulation to more closely approximate a ‘space clamp’.
Figure 3.2 : The comparison between the experiment ( [3] Figure 5) and simulation on the effect of nerve impulses on the membrane potential of glial cells. The solid symbols are resting potentials and the open symbols are depolarization potentials with different ECS K + concentrations. The comparison between the experimental data and full model is presented in detail in 67. Zhu, Y., et al.,
A Tridomain Model for Potassium Clearance in Optic Nerve. arXiv preprint arXiv:2012.03303, 2020.
4. Effects of Water Flow
In this section, we estimate ion fluxes in regions during the time when the nerve is stimulated. We examine water flow and the pattern of water circulation (driven by osmosis) and their effects on ion concentrations and flow, within and between each compartment. Numerical simulations are compared to asymptotic estimations. We compare the electrodiffusion model and the full convection-electrodiffusion model to understand how the nervous (neuron-glia) system interacts with the extracellular space to create microcirculation.
In this section, when part of the nerve is stimulated, we estimate the transmembrane fluxes and the resulting accumulation of ions in the extracellular space and glial cells. The estimations show that the variation of osmotic pressure between extracellular space and glial cells is the dominant mechanism that drives water flow. The circulation pattern and strength of water flow in optic nerve are also presented.
To simplify our discussions, we focus our analyses on an idealized setting where the stimulus is applied at an inner part of the axon compartment. As shown in Figure 4.1, the stimulus was applied at 𝑠𝑡𝑖 at a given location 𝑧 = 𝑧 . This stimulus is within the optic nerve, so 𝑟 𝑠𝑡𝑖 <𝑅 𝑎 = 𝑟 ∗ ( 𝑟 ∗ is introduced for asymptotic analysis) shown in Fig. 4.1. We distinguish the stimulated region and the non-stimulated region in the optic nerve Ω 𝑂𝑃 shown in the Fig. 4.1, since the electrical signal propagates in the 𝑧 direction in the axon compartment. We do not put the stimulus everywhere in this region, rather we only apply the stimulus at the location (𝑧 ) within a radial. Figure 4.1: Stimulated region and non-stimulated region in the optic nerve ( Ω 𝑂𝑃 ). The stimulus is applied in the axon compartment where 𝑠𝑡𝑖 at a given location 𝑧 = 𝑧 . To understand the mechanism inducing the water circulation, we first estimate the variations of ion concentrations from axon to the extracellular space during a single action potential. Then we analyze the different transmembrane current on the glial cells and identify the dominant K + current. Finally, we study osmotic pressure change after a train of action potentials on axon. We first estimate the amount of ion exchange between axon and extracellular space during a single action potential. We assume that during the single action potential, the volume fraction 𝜂 𝑙 , 𝑙 = 𝑎𝑥, 𝑔𝑙, 𝑒𝑥, does not differ from their resting state. We find then that the variation of Na + and K + in the stimulated extracellular region is the same to leading order, and that agrees with experimental observations [64,83,84]. Although our estimation is based on the classic Hodgkin-Huxley model, the methods are general and can be applied to systems with other channels and transporters. When an action potential occurs in the nerve, the equilibrium (or steady state) balance between the ions and electric fields is lost and resting state changes. We introduce notation to separate the resting state (with superscript ‘re’) before the action potential from that during the action potentials (with superscript ‘dy’). We introduce the total current of the 𝑖 th ionic species through on axon and glial membrane as 𝐼 𝑘𝑖,𝑗 = 𝑧 𝑖 𝑒 𝐽 𝑘𝑚,𝑖,𝑗 = 𝑧 𝑖 𝑒𝐽 𝑐,𝑘𝑖,𝑗 + 𝑧 𝑖 𝑒𝐽 𝑝,𝑘𝑖,𝑗 𝑖 = Na + , K + , Cl − , 𝑗 = 𝑟𝑒, 𝑑𝑦, 𝑘 = 𝑎𝑥, 𝑔𝑙 . where 𝐽 𝑝,𝑘𝑖,𝑗 is the active Na/K pump effect and 𝐽 𝑐,𝑘𝑖,𝑗 is passive source for ion 𝑖 on the axons ( 𝑘 =𝑎𝑥 ) or glial cells membranes ( 𝑘 = 𝑔𝑙 ) at resting state ( 𝑗 = 𝑟𝑒 ) before the action potential or during the action potentials ( 𝑗 = 𝑑𝑦 ). In the resting state, Na / K pump 𝐽 𝑝,𝑘𝑖,𝑟𝑒 and flux through the ion channels 𝐽 𝑐,𝑘𝑖,𝑟𝑒 in the axon membrane 𝑘 = 𝑎𝑥) and glial membrane ( 𝑘 = 𝑔𝑙) satisfy 𝐽 𝑝,𝑘𝑁𝑎,𝑟𝑒 = 3 𝐼 𝑘𝑟𝑒 𝑒 , 𝐽 𝑝,𝑘𝐾,𝑟𝑒 = −2 𝐼 𝑘𝑟𝑒 𝑒 , 𝐽 𝑝,𝑘𝐶𝑙,𝑟𝑒 = 0, 𝐽 𝑐,𝑘𝑖,𝑟𝑒 = 𝑔 𝑘𝑖,𝑟𝑒 𝑧 𝑖 𝑒 (𝑉 𝑘𝑟𝑒 − 𝐸 𝑘𝑖,𝑟𝑒 ) , 𝑘 = 𝑎𝑥, 𝑔𝑙 , 𝑖 = Na + , K + , Cl − . where the membrane potential is 𝑉 𝑘𝑟𝑒 = 𝜙 𝑘𝑟𝑒 − 𝜙 𝑒𝑥𝑟𝑒 , 𝑘 = 𝑎𝑥, 𝑔𝑙. The ion channel conductance on the axon membrane is defined as in the classical Hodgkin-Huxley model and on the glial membrane is a fixed constant, 𝑔 𝑎𝑥𝑁𝑎,𝑟𝑒 = 𝑔̅ 𝑁𝑎 (𝑚 𝑟𝑒 ) ℎ 𝑟𝑒 + 𝑔 𝑙𝑒𝑎𝑘𝑁𝑎 , 𝑔 𝑎𝑥𝐾,𝑟𝑒 = 𝑔̅ 𝐾 (𝑛 𝑟𝑒 ) + 𝑔 𝑙𝑒𝑎𝑘𝐾 , 𝑔 𝑎𝑥𝐶𝑙,𝑟𝑒 = 𝑔 𝑙𝑒𝑎𝑘𝐶𝑙 , 𝑔 𝑔𝑙𝑖,𝑟𝑒 = 𝑔 𝑔𝑙𝑖 . The kinetic variables 𝑚 𝑟𝑒 , ℎ 𝑟𝑒 and 𝑛 𝑟𝑒 are measures of the resting state open probability for the voltage-gated Na + and K + channel in the axon membrane. In addition, in the resting state, we have for ion flux through ion channel 𝐽 𝑐,𝑘𝑖,𝑟𝑒 and ion flux through pump 𝐽 𝑝,𝑘𝑖,𝑟𝑒 on glial membrane ( 𝑘 = 𝑔𝑙 ) or axon membrane ( 𝑘 = 𝑎𝑥 ) are at same magnitude, such that 5 𝑂(|𝐽 𝑝,𝑘𝑖,𝑟𝑒 |) = 𝑂(|𝐽 𝑐,𝑘𝑖,𝑟𝑒 |), 𝑖 = Na + , K + , Cl − , 𝑘 = 𝑎𝑥, 𝑔𝑙 . During a single action potential, the ion fluxes through active Na / K pump is presented as 𝐽 𝑝,𝑘𝑁𝑎,𝑑𝑦 = 3 𝐼 𝑘𝑟𝑒 +𝛿𝐼 𝑘 𝑒 , 𝐽 𝑝,𝑘𝐾,𝑑𝑦 = −2 𝐼 𝑘𝑟𝑒 +𝛿𝐼 𝑘 𝑒 , 𝑘 = 𝑎𝑥, 𝑔𝑙 , and through ion channels are 𝐽 𝑐,𝑘𝑖,𝑑𝑦 = 𝑔 𝑘𝑖,𝑑𝑦 𝑧 𝑖 𝑒 (𝑉 𝑘𝑟𝑒 − 𝐸 𝑘𝑖,𝑟𝑒 ) + 𝑔 𝑘𝑖,𝑑𝑦 𝑧 𝑖 𝑒 (𝛿𝑉 𝑘 − 𝛿𝐸 𝑘𝑖 ) , 𝑘 = 𝑎𝑥, 𝑔𝑙 , 𝑖 = Na + , K + , Cl − , where 𝛿X 𝑘 = X 𝑘𝑑𝑦 − X 𝑘𝑟𝑒 is the deviation of X away from the resting state value with X = 𝑉, 𝐸, 𝐼 on the membrane 𝑘 and 𝑔 𝑎𝑥𝑁𝑎,𝑑𝑦 = 𝑔̅ 𝑁𝑎 (𝑚 𝑑𝑦 ) ℎ 𝑑𝑦 + 𝑔 𝑙𝑒𝑎𝑘𝑁𝑎 , 𝑔 𝑎𝑥𝐾,𝑑𝑦 = 𝑔̅ 𝐾 (𝑛 𝑑𝑦 ) + 𝑔 𝑙𝑒𝑎𝑘𝐾 , 𝑔 𝑎𝑥𝐶𝑙,𝑑𝑦 = 𝑔 𝑎𝑥𝐶𝑙,𝑟𝑒 , 𝑔 𝑔𝑙𝑖,𝑑𝑦 = 𝑔 𝑔𝑙𝑖,𝑟𝑒 , 𝑖 = Na + , K + , Cl − . We claim that 𝛿𝐸 𝑎𝑥𝑖 ≪ 𝛿𝑉 𝑎𝑥 , 𝑖 = Na + , K + , Cl − (see Appendix A2) during the time the axon is firing and 𝐽 𝑝,𝑎𝑥𝑖,𝑑𝑦 ≪ 𝑔 𝑎𝑥𝑖,𝑑𝑦 𝑧 𝑖 𝑒 (𝑉 𝑎𝑥𝑟𝑒 − 𝐸 𝑎𝑥𝑖,𝑟𝑒 ), 𝑖 = Na + , K + . During the time the axon is firing, the ion conductance of the axon membrane satisfies 𝑔 𝑎𝑥𝑖,𝑟𝑒 ≪𝑔 𝑎𝑥𝑖,𝑑𝑦 , and 𝐽 𝑝,𝑎𝑥𝑖,𝑑𝑦 = 𝑂(𝐽 𝑝,𝑎𝑥𝑖,𝑟𝑒 ) = 𝑂(−𝐽 𝑐,𝑎𝑥𝑖,𝑟𝑒 ) because the variation of the concentration during a single action potential is small and the ion fluxes through the Na / K pump is controlled by its maximum current. Therefore, we can approximate the transmembrane ion currents through the axon membrane as 𝐼 𝑎𝑥𝑖,𝑑𝑦 ≈ 𝑔 𝑎𝑥𝑖,𝑑𝑦 (𝑉 ax𝑟𝑒 − 𝐸 𝑎𝑥𝑖,𝑟𝑒 ) + 𝑔 𝑎𝑥𝑖,𝑑𝑦 𝛿V ax , 𝑖 = Na + , K + , Cl − . (4.1) In the next paragraphs, eq. (4.1) is used to estimate the total axon transmembrane Na + and K + flux during a single action potential. This estimation helps us estimate the concentration changes in the extracellular space due to the axon firing. The governing equation of 𝑚 𝑑𝑦 in the Hodgkin-Huxley model is 𝑑𝑚 𝑑𝑦 𝑑𝑡 = 𝛼 𝑚 (1 − 𝑚 𝑑𝑦 ) − 𝛽 𝑚 𝑚 𝑑𝑦 , (4.2) where 𝛼 𝑚 =
110 25 − 𝛿𝑉 𝑎𝑥 exp( )−1 , 𝛽 𝑚 = 4exp (− 𝛿𝑉 𝑎𝑥 ) , (4.3) and 𝛿𝑉 𝑎𝑥 = 𝑉 𝑎𝑥𝑑𝑦 − 𝑉 𝑎𝑥𝑟𝑒 . The integral solution for eq. (4.2) is 6 𝑚 𝑑𝑦 (𝑡) = 𝑚 exp (− ∫ 𝛼 𝑚 (𝑠) + 𝛽 𝑚 (𝑠)𝑑𝑠 𝑡0 ) + ∫ 𝛼 𝑚 (𝑠) exp (− ∫ 𝛼 𝑚 (𝑢) + 𝛽 𝑚 (𝑢)𝑑𝑢 𝑡𝑠 ) 𝑑𝑠 𝑡0 . (4.4) During a single action potential period [0, 𝑇 𝑎𝑥∗ ] , we define two distinguished time intervals based on Na + channel as showed in the Fig. 4.2 below. The first period [0, 𝑡 𝑚1 ] is when the Na + channel becomes fully open status and the action membrane potential moves positive from its resting value to its most positive value, at the peak of its overshoot, as it was classically named. The second period[ 𝑡 𝑚1 , 𝑇 𝑎𝑥∗ = 𝑡 𝑚1 + 𝑡 𝑚2 ] occurs when the Na + channel closes and the action potential recovers from peak value to the hyperpolarization value. Figure 4.2: Two distinguished time intervals used in the estimation during a single action potential. The blue line is the 𝛿𝑉 𝑎𝑥 (= 𝑉 𝑎𝑥𝑑𝑦 − 𝑉 𝑎𝑥𝑟𝑒 ) during a single action potential. The dark dash line is the linear approximation of the 𝛿𝑉 𝑎𝑥 . 𝑡 𝑚1 and 𝑡 𝑚2 are the two time intervals we estimate by eq. (4.7) and eq. (4.9). In the first timer interval [0, 𝑡 𝑚1 ] , we assume that 𝛿𝑉 𝑎𝑥 increases monotonically from to 𝐸 𝑎𝑥𝑁𝑎,𝑟𝑒 − 𝑉 𝑎𝑥𝑟𝑒 , where we approximate the peak value of action potential by the Nernst potential of Na + in the resting state such that 𝛿𝑉 𝑎𝑥 (𝑡) = 𝐸 𝑎𝑥𝑁𝑎,𝑟𝑒 − 𝑉 𝑎𝑥𝑟𝑒 𝑡 𝑚1 𝑡 , 𝑡 ∈ [0, 𝑡 𝑚1 ] , (4.5) where 𝐸 𝑎𝑥𝑁𝑎,𝑟𝑒 − 𝑉 𝑎𝑥𝑟𝑒 ≈ 1.4 × 10 mV and the slope is the unknown variable in eq.(4.5). The initial value 𝑚 is chosen when 𝛿𝑉 𝑎𝑥 = 0 mV, such that 𝑚 = 𝑚 𝑟𝑒 = 𝑚 𝑒𝑞 (0) , where 𝑚 𝑒𝑞 is the equilibrium state of eq. (4.2) depending on 𝛿𝑉 𝑎𝑥 , 𝑚 𝑒𝑞 (𝛿𝑉 𝑎𝑥 ) = 𝛼 𝑚 (𝛿𝑉 𝑎𝑥 )𝛼 𝑚 (𝛿𝑉 𝑎𝑥 )+𝛽 𝑚 (𝛿𝑉 𝑎𝑥 ) . (4.6) Then by substituting eq. (4.3), eq. (4.5) into eq. (4.4) , we get 7 m dy (t m1 ) = m exp ( m1 (exp (− ) − 1 ) + t m1 [Li (exp(x)) + xln(1 − exp(x)) − x ] | ) − t m1 ∫ sexp(s)−1 exp ( m1 (exp (− ) − exp (− )) + −11.52.5 t m1 [Li (exp(x)) + xln(1 − exp(x)) − x ] | s−11.5 ) ds . (4.7) Without loss of generality, we assume the voltage-gated Na + channel is almost fully open when 𝑡 = 𝑡 𝑚1 . Therefore, at 𝑡 = 𝑡 𝑚1 , we assume that 𝑚 𝑑𝑦 (𝑡 𝑚1 ) =0.95 < 𝑚 𝑒𝑞 (𝛿𝑉 𝑎𝑥 (𝑡 𝑚1 )) . The estimation from above eq. (4.7) gives 𝑡 𝑚1 ≈ ms . In the second timer interval, we use the homogeneous property of eq. (4.2) and move the time interval [ 𝑡 𝑚1 , 𝑇 𝑎𝑥∗ = 𝑡 𝑚1 + 𝑡 𝑚2 ] to [0, 𝑡 𝑚2 ] to simplify the notation. We assume that 𝛿𝑉 𝑎𝑥 decreases monotonically from 𝐸 𝑎𝑥𝑁𝑎,𝑟𝑒 − 𝑉 𝑎𝑥𝑟𝑒 to 𝐸 𝑎𝑥𝐾,𝑟𝑒 − 𝑉 𝑎𝑥𝑟𝑒 at time 𝑡 𝑚2 such that 𝛿𝑉 𝑎𝑥 (𝑡) = 𝐸 𝑎𝑥𝑁𝑎,𝑟𝑒 − 𝑉 𝑎𝑥𝑟𝑒 − 𝐸 𝑎𝑥𝑁𝑎,𝑟𝑒 −𝐸 𝑎𝑥𝐾,𝑟𝑒 𝑡 𝑚2 𝑡 , 𝑡 ∈ [0, 𝑡 𝑚2 ] , (4.8) where 𝐸 𝑎𝑥𝑁𝑎,𝑟𝑒 − 𝐸 𝑎𝑥𝐾,𝑟𝑒 ≈ 1.5 × 10 mV . We assume that the initial value 𝑚 in the second timer period is 𝑚 = 𝑚 𝑑𝑦 (𝑡 𝑚1 ) , and Na + channel is in a nearly closed state when the 𝛿𝑉 𝑎𝑥 approaching 𝐸 𝑎𝑥𝐾,𝑟𝑒 − 𝑉 𝑎𝑥𝑟𝑒 , namely, 𝑚 𝑒𝑞 (𝛿𝑉 𝑎𝑥 (𝑡 𝑚2 )) < 𝑚(𝑡 𝑚2 ) . We let 𝑚 𝑑𝑦 (𝑡 𝑚2 ) =0.1 , so that (𝑚 𝑑𝑦 (𝑡 𝑚2 )) = 10 −3 . The estimation of 𝑡 𝑚2 comes from the following equation by substituting eq. (4.3), eq. (4.8) into eq. (4.4) m dy (t m2 ) = m exp ( m2 (exp (− ) − exp ( )) + t m2 [Li (exp(x)) + xln(1 − exp(x)) − x ] | ) + t m2 ∫ sexp(s)−1 exp ( m2 (exp (− ) − exp ( ) ) + t m2 [Li (exp(x)) + xln(1 − exp(x)) − x ] | ) ds (4.9) where the above estimation gives 𝑡 𝑚2 ≈ ms . We provide the robust analysis of estimation using eq. (4.8) and eq. (4.9) in the Appendix A3. In sum, based on estimated 𝑡 𝑚1 and 𝑡 𝑚2 , 𝛿𝑉 𝑎𝑥 and ℎ have the explicit expressions as follow, for 𝑡 ∈ [0, 𝑇 𝑎𝑥∗ = 𝑡 𝑚1 + 𝑡 𝑚2 ] 𝛿𝑉 𝑎𝑥 = { 𝐸 𝑎𝑥𝑁𝑎,𝑟𝑒 − 𝑉 𝑎𝑥𝑟𝑒 𝑡 𝑚1 t , 𝑡 ∈ [0, 𝑡 𝑚1 ] , 𝐸 𝑎𝑥𝑁𝑎,𝑟𝑒 − 𝑉 𝑎𝑥𝑟𝑒 − 𝐸 𝑎𝑥𝑁𝑎,𝑟𝑒 − 𝐸 𝑎𝑥𝐾,𝑟𝑒 𝑡 𝑚2 (𝑡 − 𝑡 𝑚1 ) , 𝑡 ∈ [𝑡 𝑚1 , 𝑇 𝑎𝑥∗ ] . and ℎ 𝑑𝑦 (𝑡) = ℎ exp (− ∫ 𝛼 ℎ (𝑠) + 𝛽 ℎ (𝑠)𝑑𝑠 𝑡0 ) + ∫ 𝛼 ℎ (𝑠) exp (− ∫ 𝛼 ℎ (𝑢) + 𝛽 ℎ (𝑢)𝑑𝑢 𝑡𝑠 ) 𝑑𝑠 𝑡0 , where 𝛼 ℎ = exp (− 𝛿𝑉 𝑎𝑥 ) , 𝛽 ℎ = )+1 , and the initial value ℎ = ℎ 𝑟𝑒 = 𝛼 ℎ (0)𝛼 ℎ (0)+𝛽 ℎ (0) . Therefore, by using eq. (4.1), we estimate the total Na + flux during a single action potential [0, 𝑇 𝑎𝑥∗ ] by ∫ 𝐼 𝑎𝑥𝑁𝑎,𝑑𝑦 𝑧 𝑁𝑎 𝑒 𝑇 𝑎𝑥∗ 𝑑𝑡 ≈ ∫ 𝑔̅ 𝑁𝑎 ℎ 𝑑𝑦 (𝑚 𝑑𝑦 ) 𝑧 𝑁𝑎 𝑒 (𝑉 𝑎𝑥𝑟𝑒 − 𝐸 𝑎𝑥𝑁𝑎,𝑟𝑒 )+ 𝑔̅ 𝑁𝑎 ℎ 𝑑𝑦 (𝑚 𝑑𝑦 ) 𝑧 𝑁𝑎 𝑒 𝛿𝑉 𝑎𝑥𝑡 𝑚1 +𝑡 𝑚2 𝑑𝑡 ≈ − × 10 −9 mol/m , (4.10) where 𝑡 𝑚1 = 0.67 ms and 𝑡 𝑚2 = 3.00 ms have been used. In the next, we have the total Cl − flux though axon membrane during a single action potential is ∫ 𝐼 𝑎𝑥𝐶𝑙,𝑑𝑦 𝑧 𝐶𝑙 𝑒𝑇 𝑎𝑥∗ 𝑑𝑡 ≈ ∫ 𝑔 𝑎𝑥𝐶𝑙 𝛿𝑉 𝑎𝑥 𝑧 𝐶𝑙 𝑒 𝑑𝑡 𝑡 𝑚1 +𝑡 𝑚2 ≈ −3.7 × 10 −10 mol/m , (4.11) since the Cl − current through axon membrane during an action potential can be approximated by 𝐼 𝑎𝑥𝐶𝑙,𝑑𝑦 = 𝑔 𝑎𝑥𝐶𝑙 (𝑉 𝑎𝑥𝑟𝑒 − 𝐸 𝑎𝑥𝐶𝑙,𝑟𝑒 ) +𝑔 𝑎𝑥𝐶𝑙 (𝛿𝑉 𝑎𝑥 − 𝛿𝐸 𝑎𝑥𝐶𝑙 ) ≈ 𝑔 𝑎𝑥𝐶𝑙 𝛿𝑉 𝑎𝑥 . where we used 𝑉 𝑎𝑥𝑟𝑒 − 𝐸 𝑎𝑥𝐶𝑙,𝑟𝑒 , 𝛿𝐸 𝑎𝑥𝐶𝑙 ≪ 𝛿V ax . We still need the estimation of the total K + flux through axon membrane during a single action potential. The governing equation of 𝜙 𝑎𝑥 yields ∑ 𝑖 𝑧 𝑖 𝑒 𝜕𝜕𝑧 (𝜂 𝑎𝑥 𝑗 𝑎𝑥𝑖 ) = −ℳ 𝑎𝑥 ( 𝐼 𝑎𝑥𝐾,𝑑𝑦 + 𝐼 𝑎𝑥𝑁𝑎,𝑑𝑦 + 𝐼 𝑎𝑥𝐶𝑙,𝑑𝑦 ), (4.12) At every location of the stimulated region, the duration of a single action potential is 𝑇 𝑎𝑥∗ . However, there also a signal transportation time 𝑇 𝑎𝑙𝑙∗ during which the electric signal propagates from the axon near the head of the optic nerve ( 𝑧 = 0 ) to far-eye-side ( 𝑧 = 𝐿 ) axon, which 9 connects the brain. By integrating right-hand side of eq. ( ) over time [0, 𝑇 𝑎𝑙𝑙∗ ] , we have −ℳ 𝑎𝑥 ∫ ∫ 𝐼 𝑎𝑥𝐾,𝑑𝑦 + 𝐼 𝑎𝑥𝑁𝑎,𝑑𝑦 + 𝐼 𝑎𝑥𝐶𝑙,𝑑𝑦𝐿0 𝑑𝑧𝑑𝑡 𝑇 𝑎𝑙𝑙∗ ≈ −ℳ 𝑎𝑥 𝐿 ∫ 𝐼 𝑎𝑥𝐾,𝑑𝑦 + 𝐼 𝑎𝑥𝑁𝑎,𝑑𝑦 + 𝐼 𝑎𝑥𝐶𝑙,𝑑𝑦𝑇 𝑎𝑥∗ 𝑑𝑡, (4.13) where we use the propagation property of the action potential along 𝑧 direction and only the firing period is taken into consideration. By integrating the left-hand side of eq. (4.12), we have ∫ ∫ ∑ 𝑖 𝑧 𝑖 𝑒 𝜕𝜕𝑧 (𝜂 𝑎𝑥 𝑗 𝑎𝑥𝑖 ) 𝐿0𝑇 𝑎𝑙𝑙∗ 𝑑𝑧𝑑𝑡 = 𝑂(𝑇 𝑎𝑙𝑙∗ 𝑒𝜂 𝑎𝑥 𝑗 𝑎𝑥𝑏𝑑 ) . (4.14) We assume that the characteristic time scale of 𝑇 𝑎𝑙𝑙∗ = O (10 −3 ). The scale of ion flux at the boundary ( 𝑧 = 0, 𝐿 ) is dominated by the diffusion term and we estimate its scale as 𝑗 𝑎𝑥𝑏𝑑 = 𝑂 (𝐷 𝑎𝑥∗ 𝛿𝑐 𝑎𝑥∗ 𝑧 ∗ ) , since the boundary conditions are 𝜕𝜙 𝑎𝑥 𝜕𝑧 | 𝑧=0,𝐿 = 0 and 𝑢 𝑎𝑥 (0) = 𝑢 𝑎𝑥 (𝐿) = 0. The 𝛿𝑐 𝑎𝑥∗ is the characteristic variation between the ion concentration at its boundary value and the ion concentration inside the axon, after a single action potential. Based on the Na + flux estimation in Eq. (4.10), we estimate 𝛿𝑐 𝑎𝑥∗ = 𝑂(10 −1 ) and we have 𝑂(T 𝑎𝑙𝑙∗ 𝜂 𝑎𝑥 𝑗 𝑎𝑥𝑏𝑑 ) ≪ 𝑂 (ℳ 𝑎𝑥 𝐿 |∫ 𝐽 𝑎𝑥𝑚,𝐶𝑙,𝑑𝑦𝑇 𝑎𝑥∗ 𝑑𝑡|) ≪ 𝑂 (ℳ 𝑎𝑥 𝐿 |∫ 𝐽 𝑎𝑥𝑚,𝑁𝑎,𝑑𝑦𝑇 𝑎𝑥∗ 𝑑𝑡 |) . (4.15) In other words, based on eq. (4.13) , eq. (4.14) and eq. (4.15), we have 𝑂 (|∫ 𝐽 𝑎𝑥𝑚,𝐾,𝑑𝑦𝑇 𝑎𝑥∗ 𝑑𝑡|) = 𝑂 (|∫ 𝐽 𝑎𝑥𝑚,𝑁𝑎,𝑑𝑦𝑇 𝑎𝑥∗ 𝑑𝑡|) , (4.16) Based on eq. (4.10), the total axon transmembrane K + flux should be ∫ 𝐽 𝑎𝑥𝑚,𝐾,𝑑𝑦𝑇 𝑎𝑥∗ 𝑑𝑡 ≈ − ∫ 𝐽 𝑎𝑥𝑚,𝑁𝑎,𝑑𝑦𝑇 𝑎𝑥∗ 𝑑𝑡 ≈ 2 × 10 −9 mol/m . (4.17) where [0, 𝑇 𝑎𝑥∗ ] is the time interval enclosing a single action potential. Remark 4.1 : Eq. (4.16) shows that for a single action potential, the leading order of the potassium flux out of the axon to the extracellular space equals the leading order of the sodium flux into the axon from the extracellular space. This estimation is consistent with observations in the literature (64,83,84). The time scale of 𝑇 𝑎𝑥∗ is milliseconds. In this short period, the cumulative ion fluxes through the axon membrane change the ion concentration in the extracellular space, namely, 𝜂 𝑒𝑥 𝛿𝑐 𝑒𝑥𝑖 = ℳ 𝑎𝑥 ∫ 𝐽 𝑎𝑥𝑚,𝑖,𝑑𝑦𝑇 𝑎𝑥∗ 𝑑𝑡, 𝑖 = Na + , K + . (4.18)
0 Based on eq. (4.17) and (4.18), the variations of Na + and K + concentrations in the extracellular space with action potentials, can be written as 𝛿𝑐 𝑠𝑡𝑖 = 𝑂 ( ℳ 𝑎𝑥 𝜂 𝑒𝑥 |∫ 𝐽 𝑎𝑥𝑚,𝑖,𝑑𝑦𝑇 𝑎𝑥∗ 𝑑𝑡|) , 𝑖 = Na + , K + . (4.19) In the following discussion, we use 𝛿𝐶 𝑠𝑡𝑖 describes the concentration changes in the stimulated extracellular space after a single action potential, 𝛿𝑐 𝑠𝑡𝑖 = 0.12 mM . In this section, we estimate the changes in the glial transmembrane current when the concentration varies by 𝛿𝑐 𝑠𝑡𝑖 in the stimulated extracellular region. We also find that the electric field 𝜙 𝑔𝑙 responds immediately to changes of the glial potassium Nernst potential. The variation of extracellular electric potential 𝛿𝜙 𝑒𝑥 is small in the stimulated region compared to the variation of glial electric potential 𝛿𝜙 𝑔𝑙 ,. We conclude that the dominant current through the glial membrane in the stimulated region is through the passive K + channel, rather than the Na/K pump. At the same time, in the non- stimulated extracellular region, almost the same amount of K + moves from the glial compartment to extracellular space. In other words, both the glial cells and extracellular space in the non-stimulated region participate in the spatial buffering process to help potassium clearance [79, 80]. In the stimulated region, the Nernst potential of K + across the glial membrane changes because of additional potassium 𝛿𝑐 𝑒𝑥𝐾 in the extracellular space 𝛿𝐸 𝑔𝑙𝐾 = 𝑘 𝐵 𝑇𝑧 𝐾 𝑒 ( log (1 + 𝛿𝑐 𝑒𝑥𝐾 𝑐 𝑒𝑥𝐾,𝑟𝑒 ) − log (1 + 𝛿𝑐 𝑔𝑙𝐾 𝑐 𝑔𝑙𝐾,𝑟𝑒 ) ) , (4.20) where 𝛿𝑐 𝑙𝐾 , 𝑙 = 𝑒𝑥, 𝑔𝑙 are the variations of concentrations in the 𝑙 compartment. The variation of glial potassium concentration 𝛿𝑐 𝑔𝑙𝐾 is a result of the 𝛿𝑐 𝑒𝑥𝐾 produced by the glial transmembrane flux. Recall that that the volume fraction ( 𝜂 𝑔𝑙 ) of the glial compartment is much larger than the extracellular space ( 𝜂 𝑒𝑥 ). During a single action potential, we have then 𝛿𝑐 𝑔𝑙𝐾 < 𝛿𝑐 𝑒𝑥𝐾 ≪ 𝑐 𝑒𝑥𝐾,𝑟𝑒 ≪ 𝑐 𝑔𝑙𝐾,𝑟𝑒 . The 𝛿𝐸 𝑔𝑙𝐾 in eq. (4.20) can be approximated by a Taylor expansion. 𝛿𝐸 𝑔𝑙𝐾 ≈ 𝑘 𝐵 𝑇𝑧 𝐾 𝑒 𝛿𝑐 𝑒𝑥𝐾 𝑐 𝑒𝑥𝐾,𝑟𝑒 . (4.21)
1 At the same time, the variation in the potassium Nernst potential eq. (4.21) in the stimulated region, induces the changes of the glial membrane potential 𝑉 𝑔𝑙 and glial compartment electric potential 𝜙 𝑔𝑙 . We move on now to estimate the variations of potential in the extracellular space and glial compartment. The governing equation for 𝜙 𝑒𝑥 yields ∑ 𝑖 𝑧 𝑖 𝑒 ▽⋅ (𝜂 𝑒𝑥 𝒋 𝑒𝑥𝑖 ) = ∑ 𝑖 𝑧 𝑖 𝑒ℳ 𝑔𝑙 (𝐽 𝑝,𝑔𝑙𝑖 + 𝐽 𝑐,𝑔𝑙𝑖 ) + ∑ 𝑖 𝑧 𝑖 𝑒ℳ 𝑎𝑥 (𝐽 𝑝,𝑎𝑥𝑖 + 𝐽 𝑐,𝑎𝑥𝑖 ), (4.22) where 𝒋 𝑒𝑥𝑖 = 𝑐 𝑒𝑥𝑖 𝒖 𝑒𝑥 − 𝐷 𝑒𝑥𝑖 𝜏 𝑒𝑥 ( ▽ 𝑐 𝑒𝑥𝑖 + 𝑧 𝑖 𝑒𝑘 𝐵 𝑇 𝑐 𝑒𝑥𝑖 ▽ 𝜙 𝑒𝑥 ). We claim that after the axon stops firing, the major glial transmembrane current is through glial membrane K + channels, such that ∑ 𝑖 𝑧 𝑖 𝑒ℳ 𝑔𝑙 (𝐽 𝑝,𝑔𝑙𝑖 + 𝐽 𝑐,𝑔𝑙𝑖 ) + ∑ 𝑖 𝑧 𝑖 𝑒ℳ 𝑎𝑥 (𝐽 𝑝,𝑎𝑥𝑖 + 𝐽 𝑐,𝑎𝑥𝑖 ) ≈ ℳ 𝑔𝑙 𝑔 𝑔𝑙𝐾 (𝛿𝑉 𝑔𝑙 − 𝛿𝐸 𝑔𝑙𝐾 ) . (4.23) We provide detail estimation in the Appendix A4. Next, we integrate eq. ( ) over the stimulated region 𝑉 𝑆 = {(𝑟, 𝑧, 𝜃)| 𝑟 ∈ [0, 𝑟 𝑠𝑡𝑖 ], 𝑧 ∈[0, 𝐿], θ ∈ [0,2π]} , through which the action potential propagates shown in Fig. 4.1. By eq. ( ) , we have the total current approximation as ∫ ℳ 𝑔𝑙 𝑔 𝑔𝑙𝐾 ( 𝛿𝑉 𝑔𝑙 − 𝛿𝐸 𝑔𝑙𝐾 ) 𝑉 𝑆 𝑑𝑣 ≈ 𝜋𝑟 𝑠𝑡𝑖2 𝐿ℳ 𝑔𝑙 𝑔 𝑔𝑙𝐾 (𝛿𝑉 𝑔𝑙 − 𝛿𝐸 𝑔𝑙𝐾 ). (4.24) For the left-hand side of eq. ( ) , by the charge neutrality assumption in the extracellular space, we naturally have ∑ 𝑧 𝑖 𝑒 𝑖 𝑐 𝑒𝑥𝑖 𝒖 𝑒𝑥 = 0 . Based on eq. (4.11) and eq. (4.17), we know after a single action potential that the extracellular concentration changes as O( 𝛿𝑐 𝑒𝑥𝑁𝑎 ) = −𝛿𝑐 𝑠𝑡𝑖 , 𝑂(𝛿𝑐 𝑒𝑥𝐾 ) = 𝛿𝑐 𝑠𝑡𝑖 , 𝑂(𝛿𝑐 𝑒𝑥𝐶𝑙 ) ≪ 𝛿𝑐 𝑠𝑡𝑖 . (4.25) Using eq. (4.25) and boundary condition eq. (2.48), the diffusion term in eq. ( ) can be approximated as − ∫ ∑ 𝑧 𝑖𝑖 𝑒 ▽⋅ (𝜂 𝑒𝑥 𝐷 𝑒𝑥𝑖 𝜏 𝑒𝑥 ▽ 𝑐 𝑒𝑥𝑖 ) 𝑑𝑣 𝑉 𝑆 ≈ 2𝜋𝑟 𝑠𝑡𝑖 𝑒𝐿𝜂 𝑒𝑥 𝐷 𝑒𝑥diff 𝜏 𝑒𝑥 𝛿𝑐 𝑠𝑡𝑖 𝑟 ∗ , (4.26) where 𝐷 𝑒𝑥diff = 𝐷 𝑒𝑥𝐾 − 𝐷 𝑒𝑥𝑁𝑎 . In eq. (4.26), we have assumed that the total current through the left 2 boundary ( 𝑧 = 0 ) and right boundary ( 𝑧 = 𝐿 ) (of the stimulated region 𝑉 𝑆 ) is much smaller compared to the current in the radial direction because of the different length scaling in the 𝑧 and 𝑟 direction. The transition region 𝑆 𝑡 = {(𝑟, 𝑧, 𝜃)|𝑟 = 𝑟 𝑠𝑡𝑖 , 𝑧 ∈ [0, 𝐿], 𝜃 ∈ [0,2𝜋]} area is much larger than the left and right boundary area of 𝑉 𝑆 . Similarly, the integration of the electric drift term in eq. ( ) becomes − ∫ ∑ 𝑧 𝑖𝑖 𝑒 ▽⋅ (𝜂 𝑒𝑥 𝐷 𝑒𝑥𝑖 𝜏 𝑒𝑥 𝑧 𝑖 𝑒𝑘 𝐵 𝑇 𝑐 𝑒𝑥𝑖 ▽ 𝜙 𝑒𝑥 ) 𝑑𝑣 𝑉 𝑆 = 2𝜋𝑟 𝑠𝑡𝑖 𝐿 𝜂 𝑒𝑥 𝜏 𝑒𝑥 𝑒 𝑘 𝐵 𝑇 𝜎 𝑒𝑥 𝛿𝜙 𝑒𝑥 𝑟 ∗ , (4.27) where 𝜎 𝑒𝑥 = ∑ (𝑧 𝑖 ) 𝐷 𝑒𝑥𝑖 𝑐 𝑒𝑥𝑖 . Combining equations eq. ( ) , eq. ( ) and eq. ( ) , we have 𝑠𝑡𝑖 ( 𝜂 𝑒𝑥 𝜏 𝑒𝑥 D exdiff 𝑀 𝑔𝑙 𝛿𝑐 𝑠𝑡𝑖 𝑟 ∗ + 𝜂 𝑒𝑥 𝜏 𝑒𝑥 𝜎 𝑒𝑥 𝑒𝑀 𝑔𝑙 𝑘 𝐵 𝑇 𝛿𝜙 𝑒𝑥 𝑟 ∗ ) ≈ 𝑔 𝑔𝑙𝐾 𝑒 ( 𝛿𝑉 𝑔𝑙 − 𝛿𝐸 𝑔𝑙𝐾 ) (4.28) In the same way, from the governing equation of 𝜙 𝑔𝑙 , ∑ 𝑖 𝑧 𝑖 𝑒 ▽⋅ (𝜂 𝑔𝑙 𝒋 𝑔𝑙𝑖 ) = − ∑ 𝑖 𝑧 𝑖 𝑒ℳ 𝑔𝑙 (𝐽 𝑝,𝑔𝑙𝑖 + 𝐽 𝑐,𝑔𝑙𝑖 ) , where 𝒋 𝑔𝑙𝑖 = 𝑐 𝑔𝑙𝑖 𝒖 𝑔𝑙 − 𝐷 𝑔𝑙𝑖 𝜏 𝑔𝑙 (▽ 𝑐 𝑔𝑙𝑖 + 𝑧 𝑖 𝑒𝑘 𝐵 𝑇 𝑐 𝑔𝑙𝑖 ∇𝜙 𝑔𝑙 ) . We obtain the estimation as − 𝑠𝑡𝑖 𝜂 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝑒𝑀 𝑔𝑙 𝑘 𝐵 𝑇 𝛿𝜙 𝑔𝑙 𝑟 ∗ ≈ 𝑔 𝑔𝑙𝐾 𝑒 ( 𝛿𝑉 𝑔𝑙 − 𝛿𝐸 𝑔𝑙𝐾 ) , (4.29) where 𝜎 𝑔𝑙 = ∑ (𝑧 𝑖 ) 𝐷 𝑔𝑙𝑖 𝑐 𝑔𝑙𝑖 . We neglect the diffusion and convection terms in the glial compartment because these require much longer time to respond to the extracellular concentration change. Therefore, based on eq. (4.28) and eq. (4.29) , we have 𝛿𝜙 𝑒𝑥 = − 𝜂 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝜂 𝑒𝑥 𝜏 𝑒𝑥 𝜎 𝑒𝑥 𝛿𝜙 𝑔𝑙 − D exdiff 𝑘 𝐵 𝑇𝛿𝑐 𝑠𝑡𝑖 𝜎 𝑒𝑥 𝑒 , (4.30) where 𝑂 ( 𝜂 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝜂 𝑒𝑥 𝜏 𝑒𝑥 𝜎 𝑒𝑥 ) = 10 −2 and 𝑂 ( D exdiff 𝑘 𝐵 𝑇𝛿𝑐 𝑠𝑡𝑖 𝜎 𝑒𝑥 𝑒 ) = 10 −6 . By matching orders, we claim that |𝛿𝜙 𝑔𝑙 | ≫ |𝛿𝜙 𝑒𝑥 | (see Appendix A5) and 𝑂(𝛿𝑉 𝑔𝑙 ) = 𝑂(𝛿𝜙 𝑔𝑙 − 𝛿𝜙 𝑒𝑥 ) = 𝑂(𝛿𝜙 𝑔𝑙 ) . (4.31) In the next step, we estimate the leaking potassium current through the glial membrane. Based on eq. ( ) , we have 3 𝑔 𝑔𝑙𝐾 𝑒 ( 𝛿𝜙 𝑔𝑙 − 𝛿𝐸 𝑔𝑙𝐾 ) ≈ 𝑔 𝑔𝑙𝐾 𝑒 ( 𝛿𝑉 𝑔𝑙 − 𝛿𝐸 𝑔𝑙𝐾 ) ≈ − 2𝑟 𝑠𝑡𝑖 𝜂 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝑒ℳ 𝑔𝑙 𝑘 𝐵 𝑇 𝛿𝜙 𝑔𝑙 𝑟 ∗ . Hence, we have the relation between 𝛿𝐸 𝑔𝑙𝐾 and 𝛿𝜙 𝑔𝑙 𝛿𝐸 𝑔𝑙𝐾 ≈ (1 + ℎ 𝜖 ) 𝛿𝜙 𝑔𝑙 , (4.32) where ℎ 𝜖 = 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝑒 𝑟 𝑠𝑡𝑖 𝑀 𝑔𝑙 𝑘 𝐵 𝑇𝑔 𝑔𝑙𝐾 𝑟 ∗ . Therefore, we have 𝑔 𝑔𝑙𝐾 (𝛿𝑉 𝑔𝑙 − 𝛿𝐸 𝑔𝑙𝐾 ) = − 𝑔 𝑔𝑙𝐾 ℎ 𝜖 𝜖 𝛿𝐸 𝑔𝑙𝐾 . (4.33) Furthermore, from eq. (4.30), eq. (4.32) and eq. (4.21), we have 𝛿𝜙 𝑒𝑥 ≈ − 𝜂 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝜂 𝑒𝑥 𝜏 𝑒𝑥 𝜎 𝑒𝑥 𝜖 ) 𝑘 𝐵 𝑇𝑧 𝐾 𝑒 𝛿𝑐 𝑒𝑥𝐾 𝑐 𝑒𝑥𝐾,𝑟𝑒 − D exdiff 𝑘 𝐵 𝑇𝛿𝑐 𝑠𝑡𝑖 𝜎 𝑒𝑥 𝑒 ≈ − 𝜂 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝜂 𝑒𝑥 𝜏 𝑒𝑥 𝜎 𝑒𝑥 𝜖 ) 𝑘 𝐵 𝑇𝑧 𝐾 𝑒 𝛿𝑐 𝑒𝑥𝐾 𝑐 𝑒𝑥𝐾,𝑟𝑒 , (4.34) since we have D exdiff ≪ 𝜂 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝜂 𝑒𝑥 𝜏 𝑒𝑥 𝑐 𝑒𝑥𝐾,𝑟𝑒 (1+ℎ 𝜖 ) . The variation of electric field 𝛿𝜙 𝑔𝑙 in both stimulated and non-stimulated regions is produced (without delay) by the change in 𝛿𝐸 𝑔𝑙𝐾 in the stimulated region, as described by the governing equation for 𝜙 𝑔𝑙 in eq. (2.42). The potassium leak current is the major current through the glial membrane in the non-stimulated region as it is in the stimulated region, because the flux through the ion channel is voltage 𝜙 𝑔𝑙 dependent and 𝑔 𝑔𝑙𝐶𝑙 , 𝑔 𝑔𝑙𝑁𝑎 << 𝑔 𝑔𝑙𝐾 . In the next pages, we introduce the superscript notation ‘ S ’ for the stimulated region variables and superscript ‘ NS ’ for non -stimulated region ones. For the glial transmembrane fluxes, we have ∑ 𝑖 𝑧 𝑖 𝑒ℳ 𝑔𝑙 (𝐽 𝑝,𝑔𝑙𝑆,𝑖 + 𝐽 𝑐,𝑔𝑙𝑆,𝑖 ) ≈ ℳ 𝑔𝑙 𝑔 𝑔𝑙𝐾 ( 𝛿𝑉 𝑔𝑙𝑆 − 𝛿𝐸 𝑔𝑙𝑆,𝐾 ) , ∑ 𝑖 𝑧 𝑖 𝑒ℳ 𝑔𝑙 (𝐽 𝑝,𝑔𝑙𝑁𝑆,𝑖 + 𝐽 𝑐,𝑔𝑙𝑁𝑆,𝑖 ) ≈ ℳ 𝑔𝑙 𝑔 𝑔𝑙𝐾 ( 𝛿𝑉 𝑔𝑙𝑁𝑆 − 𝛿𝐸 𝑔𝑙𝑁𝑆,𝐾 ). By integration of the 𝜙 𝑔𝑙 eq. (2.42) over the stimulated region 𝑉 𝑆 and integration over the non- stimulated region 𝑉 𝑁𝑆 (respectively), we have 4 { ∫ ∑ 𝑖 𝑧 𝑖 e ▽⋅ (𝜂 𝑔𝑙 𝒋 𝑔𝑙𝑆,𝑖 )𝑑𝑣 𝑉 𝑆 ≈ ∫ ℳ 𝑔𝑙 𝑔 𝑔𝑙𝐾 ( 𝛿𝑉 𝑔𝑙𝑆 − 𝛿𝐸 𝑔𝑙𝑆,𝐾 ) 𝑉 𝑆 ,∫ ∑ 𝑖 𝑧 𝑖 e ▽⋅ (𝜂 𝑔𝑙 𝒋 𝑔𝑙𝑁𝑆,𝑖 )𝑑𝑣 𝑉 𝑁𝑆 ≈ ∫ ℳ 𝑔𝑙 𝑔 𝑔𝑙𝐾 ( 𝛿𝑉 𝑔𝑙𝑁𝑆𝑉 𝑁𝑆 − 𝛿𝐸 𝑔𝑙𝑁𝑆,𝐾 ), (4.35) Most of the current flows between region 𝑉 𝑆 and region 𝑉 𝑁𝑆 are through the transition region 𝑆 𝑡 = {(𝑟, 𝑧, 𝜃)|𝑟 = 𝑟 𝑠𝑡𝑖 , 𝑧 ∈ [0, 𝐿], 𝜃 ∈ [0,2𝜋]} . The flows through the left and right boundaries of 𝑉 𝑆 and 𝑉 𝑁𝑆 are much smaller compared to that through the transition region 𝑆 𝑡 . Then by eq. (4.35) and boundary conditions for 𝜙 𝑔𝑙 we have ∫ ℳ 𝑔𝑙 𝑔 𝑔𝑙𝐾 ( 𝛿𝑉 𝑔𝑙𝑆 − 𝛿𝐸 𝑔𝑙𝑆,𝐾 ) 𝑉 𝑆 𝑑𝑣 ≈ − ∫ ℳ 𝑔𝑙 𝑔 𝑔𝑙𝐾 (𝛿𝑉 𝑔𝑙𝑁𝑆 − 𝛿𝐸 𝑔𝑙𝑁𝑆,𝐾 ) 𝑉 𝑁𝑆 𝑑𝑣 . (4.36) The average K + flux through the glial membrane in the non-stimulated region is defined by 𝑔 𝑔𝑙𝐾 𝑧 𝐾 𝑒 (𝛿𝑉 𝑔𝑙𝑁𝑆 − 𝛿𝐸 𝑔𝑙𝑁𝑆,𝐾 ) = − 𝑟 𝑠𝑡𝑖2 𝑟 ∗2 −𝑟 𝑠𝑡𝑖2 𝑔 𝑔𝑙𝐾 𝑧 𝐾 𝑒 (𝛿𝑉 𝑔𝑙𝑆 − 𝛿𝐸 𝑔𝑙𝑆,𝐾 ) , (4.37) In summary, eq. (4.36) and eq. (4.37), show how the glial compartment in the non-stimulated region and the extracellular space serve as spatial buffers and help clear potassium from the extracellular space outside the stimulated axons [81]. Remark 4.2: The glial compartment serves as an important and quick potassium transport device to remove accumulated potassium during the axon firing as shown in Fig. 4.3a. In the stimulated region, the change in the potassium Nernst potential change makes the glial membrane potential more positive and moves potassium through ion channels into the glial compartment. In the non-stimulated region, since glia is an electrical syncytium, the glial membrane potential simultaneously increases as it does in the stimulated region. However, the glia potassium Nernst potential in the non-stimulated region is not very different from that in the resting state. These potentials produce an outward potassium flux from the glial compartment in the non-stimulated region. Interacting regions of this sort depend on spatial variables and the properties of the glia as a syncytium. It is difficult to capture these effects in models that do not include space as an independent variable. Even if such compartment models capture these effects correctly in one set of conditions (because parameters are chose to make the description correct), they are unlikely to describe the effects of changes in conditions consistently, including membrane potential.
In this section, we discuss the water circulation between the stimulated region and non- stimulated region. As extra K + moves from axon to extracellular space, accumulated ions produce an osmotic pressure difference between the intra- and inter- domain. The osmotic pressure difference drives the transmembrane water flow and circulation.
5 Now we consider a train of stimulus stimulated at frequency 𝑓 𝑚 in the axon region ( 𝑠𝑡𝑖 , 𝑧 = 𝑧 ) during time [0, 𝑇 𝑠𝑡𝑖 ]. The estimation depends on the potassium and sodium variation of concentration in the extracellular space and charge neutrality condition. The clearance of the extra amount of potassium 𝛿𝐶 𝑒𝑥𝐾 in the stimulated extracellular space mostly goes through glial membrane and extracellular pathway (see Appendix A6) as 𝑑(𝜂 𝑒𝑥 𝛿𝑐 𝑒𝑥𝐾 )𝑑𝑡 = −(𝜆 𝑔𝑙𝑚,𝐾 + 𝜆 𝑒𝑥𝐾 ) 𝛿𝑐 𝑒𝑥𝐾 , (4.38) where 𝜆 𝑔𝑙𝑚,𝐾 ≈ ℳ 𝑔𝑙 𝑔 𝑔𝑙𝐾 ℎ 𝜖 𝑘 𝐵 𝑇𝑧 𝐾 (1+ℎ 𝜖 )𝑒 𝑐 𝑒𝑥𝐾,𝑟𝑒 and 𝜆 𝑒𝑥𝐾 = 𝑒𝑥𝐾 𝜏 𝑒𝑥 𝜂 𝑒𝑥 𝑟 ∗ 𝑟 𝑠𝑡𝑖 . (4.39) The 𝜆 𝑒𝑥𝐾 describes the spatial effect of extracellular communication between stimulated region and non- stimulated region. This spatial effect is not negligible since 𝜆 𝑒𝑥𝐾 is in a comparable level to the effect of glial transmembrane flux 𝜆 𝑔𝑙𝑚,𝐾 and depends on the stimulus location. The initial value of eq. (4.38) starts with the first stimulus on axon as 𝛿𝑐 𝑒𝑥𝐾 (0) = 𝛿𝑐 𝑠𝑡𝑖 , and at the beginning of each stimulus, there is an additional 𝛿𝑐 𝑠𝑡𝑖 amount of K + accumulated in the extracellular space due to the axon firing as 𝛿𝑐 𝑒𝑥𝐾 (𝑖𝑇) = 𝛿𝑐 𝑒𝑥𝐾 (𝑖𝑇) + 𝛿𝑐 𝑠𝑡𝑖 , 𝑖 = 1 ⋯ 𝑛 − 1. where 𝑛 ( = 𝑇 𝑠𝑡𝑖 𝑓 𝑚 ) is a total number of periodic periods. In the above, we view the extracellular K + concentration changes due axon firing as a source term 𝛿𝑐 𝑠𝑡𝑖 . Remark 4.3:
The concentration in the stimulated extracellular region changes rapidly because of the transmembrane action potentials, as well as the extracellular electric potential 𝜙 𝑒𝑥 . The effect of fluid circulation is the cumulative result of the above 𝛿𝑂 𝑒𝑥 in eq. (4.46). The fluid flows from the non-stimulated region to the stimulated region are dominated by the trans-glia-membrane flow. So, the convection in the extracellular reduces (i.e., flattens) the variation of osmotic pressure. Remark 4.4:
These effects make our spatially inhomogeneous model quite different from existing ODE models [63, 82], since those ODE models either take the extracellular ion concentration as constant or they do not consider the ion exchange between the extracellular space at all. In a recent work, Marte J. et al [83] introduce a compartment model similar to eq. (4.38) by considering ion flux between neuron, glia and extracellular regions in both the dendrite and soma region. It is always possible to take a field theory and approximate its 𝑥 dependence into compartments. But it is quite difficult to know how to describe the parameter dependence, and compartment inter-dependence in such models consistently. Field theories show the interdependence as outputs of the analysis. Because field models are 6 consistent, and their solutions are unique, parameter dependence and compartmental interdependence is unique. In compartment models, different assumptions are possible and difficult to compare. Analysis with assumed compartments is likely then to give different results in the hands of different investigators, slowing progress. Field models have many fewer assumptions and are more productive. However, they involve considerably more mathematical analysis [69, 70] and numerical difficulties. The time course of variation of sodium ( 𝛿𝑐 𝑒𝑥𝑁𝑎 ) in the extracellular space is 𝑑(𝜂 𝑒𝑥 𝛿𝑐 𝑒𝑥𝑁𝑎 )𝑑𝑡 = − 𝜆 𝑒𝑥𝑁𝑎,1 𝛿𝑐 𝑒𝑥𝑁𝑎 + 𝜆 𝑒𝑥𝑁𝑎,2 𝛿𝑐 𝑒𝑥𝐾 , (4.40) and 𝛿𝑐 𝑒𝑥𝑁𝑎 (0) = −𝛿𝑐 𝑠𝑡𝑖 , 𝛿𝑐 𝑒𝑥𝑁𝑎 (𝑖𝑇) = 𝛿𝑐 𝑒𝑥𝑁𝑎 (𝑖𝑇) − 𝛿𝑐 𝑠𝑡𝑖 , 𝑖 = 1 ⋯ 𝑛 − 1 . In eq. (4.40), 𝜆 𝑒𝑥𝑁𝑎,1 is the effect of extracellular diffusion and 𝜆 𝑒𝑥𝑁𝑎,2 is the effect of extracellular electric drift between stimulated and non-stimulated regions as (see Appendix A6) 𝜆 𝑒𝑥𝑁𝑎,1 = 𝑒𝑥𝑁𝑎 𝜏 𝑒𝑥 𝜂 𝑒𝑥 𝑟 ∗ 𝑟 𝑠𝑡𝑖 , 𝜆 𝑒𝑥𝑁𝑎,2 = 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝑐 𝑒𝑥𝑁𝑎,𝑟𝑒 𝐷 𝑒𝑥𝑁𝑎 𝜎 𝑒𝑥 (1+ℎ 𝜖 )𝑐 𝑒𝑥𝐾,𝑟𝑒 𝑟 ∗ 𝑟 𝑠𝑡𝑖 . (4.41) In addition, for K + flux in the extracellular space, the electric drift has a much smaller magnitude compared to diffusion flux, since the ratio 𝑅 𝑒𝑥𝐾 between the scale of electric drift and scale of diffusion for K + is (see Appendix A6) 𝑅 𝑒𝑥𝐾 = 𝜂 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝜂 𝑒𝑥 𝜏 𝑒𝑥 𝜎 𝑒𝑥 (1+ℎ 𝜖 ) = 𝑜(1) . (4.42) However, for Na + flux in the extracellular space, the magnitude of electric drift is comparable to diffusion flux as (see Appendix A6). K + and Na + are so different because of the dramatically different concentrations in extracellular solutions. 𝑅 𝑒𝑥𝑁𝑎 = 𝜂 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝜂 𝑒𝑥 𝜏 𝑒𝑥 𝜎 𝑒𝑥 (1+ℎ 𝜖 ) 𝑐 𝑒𝑥𝑁𝑎 𝑐 𝑒𝑥𝐾 = 𝑂(1) . (4.43) We solve the coupled system (4.38) and (4.40) at the time 𝑇 𝑠𝑡𝑖 = 𝑛𝑇 , yields (see Appendix A6) 𝛿𝑐 𝑒𝑥𝐾 (𝑛𝑇) = 𝛿𝑐 𝑠𝑡𝑖 ∑ exp(𝑖𝐴 𝑇) 𝑛𝑖=1 = 𝛿𝑐 𝑠𝑡𝑖 exp(𝐴 𝑇)−exp((𝑛+1)𝐴 𝑇)1−exp(𝐴 𝑇) , (4.44) and 𝛿𝑐 𝑒𝑥𝑁𝑎 (𝑛𝑇) = ∑ 𝐴 𝛿𝑐 𝑒𝑥𝐾 ((𝑖−1)𝑇)𝐴 −𝐴 (exp(𝐴 𝑇) − exp(𝐴 𝑇)) 𝑛𝑖=1 exp( (𝑛 − 𝑖)𝐴 𝑇 ) − 𝛿𝑐 𝑠𝑡𝑖 ∑ exp(𝑖𝐴 𝑇) 𝑛𝑖=1 . (4.45) where 7 𝛿𝑐 𝑒𝑥𝐾 (𝑗𝑇) = 𝛿𝑐 𝑠𝑡𝑖 1−exp((𝑗+1)𝐴 𝑇)1−exp(𝐴 𝑇) , 𝑗 = 0,1, … 𝑛 − 1. By the charge neutrality condition, the osmotic variation in the stimulated extracellular space 𝛿𝑂 𝑒𝑥 is 𝛿𝑂 𝑒𝑥 = 2(𝛿𝑐 𝑒𝑥𝐾 + 𝛿𝑐 𝑒𝑥𝑁𝑎 ), (4.46) where 𝛿𝑐 𝑒𝑥𝐾 and 𝛿𝑐 𝑒𝑥𝑁𝑎 is written in eq. (4.44) and eq. (4.45). In the following, we provide the estimations on the scales of the glial transmembrane velocity, glial compartment radial velocity, and extracellular space radial velocity. The variation in osmotic pressure in the stimulated region is the driving force for the water flow. Our estimation is based on the equations governing fluid flow and the spatial variation of osmotic pressure.
From the conservation of mass in glial compartment, we have 𝜕𝜂 𝑔𝑙 𝜕𝑡 + ℳ 𝑔𝑙 𝐿 𝑔𝑙𝑚 (𝑝 𝑔𝑙 − 𝑝 𝑒𝑥 − 𝛾 𝑔𝑙 𝑘 𝐵 𝑇(𝑂 𝑔𝑙 − 𝑂 𝑒𝑥 )) +▽⋅ (𝜂 𝑔𝑙 𝒖 𝑔𝑙 ) = 0. (4.47) By (4.46), we know that at 𝑡 = 𝑇 𝑠𝑡𝑖 , there is cumulative variation 𝛿𝑂 𝑒𝑥 (𝑇 𝑠𝑡𝑖 ) in the stimulated region. Since the glial compartment volume fraction is larger than the extracellular volume fraction, we have |𝛿𝑂 𝑔𝑙 (𝑇 𝑠𝑡𝑖 )| < |𝛿𝑂 𝑒𝑥 (𝑇 𝑠𝑡𝑖 )|. Therefore, we assume that the 𝛿𝑂 𝑒𝑥 is the driving force for hydrostatic pressure variation. Since at the resting state, eq. (4.47) yields ℳ 𝑔𝑙 𝐿 𝑔𝑙𝑚 (𝑝 𝑔𝑙𝑟𝑒 − 𝑝 𝑒𝑥𝑟𝑒 − 𝛾 𝑔𝑙 𝑘 𝐵 𝑇(𝑂 𝑔𝑙𝑟𝑒 − 𝑂 𝑒𝑥𝑟𝑒 )) +▽⋅ (𝜂 𝑔𝑙𝑟𝑒 𝒖 𝑔𝑙𝑟𝑒 ) = 0, we have 𝜕𝛿𝜂 𝑔𝑙 𝜕𝑡 + ℳ 𝑔𝑙 𝐿 𝑔𝑙𝑚 (𝛿𝑝 𝑔𝑙 − 𝛿𝑝 𝑒𝑥 − 𝛾 𝑔𝑙 𝑘 𝐵 𝑇(𝛿𝑂 𝑔𝑙 − 𝛿𝑂 𝑒𝑥 )) +▽⋅ (𝛿(𝜂 𝑔𝑙 𝒖 𝑔𝑙 )) = 0. (4.48) From the scales in the second term and third term in eq. (4.48), we have ℳ 𝑔𝑙 𝐿 𝑔𝑙𝑚 ≫ 𝜂 𝑔𝑙 𝜅 𝑔𝑙 𝜏 𝑔𝑙 𝜇(𝑟 ∗ ) . (4.49) Therefore, eq. (4.48), in the glial stimulated region of the optic nerve can be approximated by 𝜕(𝛿𝑝 𝑔𝑙 −𝛿𝑝 𝑒𝑥 )𝐾 𝑔𝑙 𝜕𝑡 + ℳ 𝑔𝑙 𝐿 𝑔𝑙𝑚 (𝛿𝑝 𝑔𝑙 − 𝛿𝑝 𝑒𝑥 ) + ℳ 𝑔𝑙 𝐿 𝑔𝑙𝑚 𝛾 𝑔𝑙 𝑘 𝐵 𝑇𝛿𝑂 𝑒𝑥 = 0, (4.50) where we use the hydraulic pressure and volume fraction relation 𝐾 𝑔𝑙 𝛿𝜂 𝑔𝑙 = 𝛿𝑝 𝑔𝑙 − 𝛿𝑝 𝑒𝑥 .
8 We approximate 𝛿𝑂 𝑒𝑥 by a linear function as 𝛿𝑂 𝑒𝑥 (𝑡) = 𝛿𝑂 𝑒𝑥 (𝑇 𝑠𝑡𝑖 )𝑇 𝑠𝑡𝑖 𝑡, 𝑡 ∈ [0, 𝑇 𝑠𝑡𝑖 ], So, the solution of 𝛿𝑝 𝑔𝑙 − 𝛿𝑝 𝑒𝑥 in eq. (4.50) can be written as 𝛿𝑝 𝑔𝑙 (𝑡) − 𝛿𝑝 𝑒𝑥 (𝑡) = ( 𝐵𝑡𝐴 exp(𝐴𝑡) − 𝐵𝐴 (exp(𝐴𝑡) − 1)) exp (−𝐴𝑡) , (4.51) where 𝐴 = ℳ 𝑔𝑙 𝐿 𝑔𝑙𝑚 𝐾 𝑔𝑙 , 𝐵 = −𝐾 𝑔𝑙 ℳ 𝑔𝑙 𝐿 𝑔𝑙𝑚 𝛾 𝑔𝑙 𝑘 𝐵 𝑇 𝛿𝑂 𝑒𝑥 (𝑇 𝑠𝑡𝑖 )𝑇 𝑠𝑡𝑖 . with initial value as 𝛿𝑝 𝑔𝑙 (0) − 𝛿𝑝 𝑒𝑥 (0) = 0. So the average transmembrane water flow velocity through the glial membrane in the stimulated region can be approximated as 𝑈 𝑔𝑙𝑚 (𝑡) = 𝐿 𝑔𝑙𝑚 (𝛿𝑝 𝑔𝑙 (𝑡) − 𝛿𝑝 𝑒𝑥 (𝑡) + 𝛾 𝑔𝑙 𝑘 𝐵 𝑇𝛿𝑂 𝑒𝑥 (𝑡)) . (4.52) and the scale of glial transmembrane velocity in stimulated region as 𝑈 𝑔𝑙𝑚,∗ = |𝐿 𝑔𝑙𝑚 (𝛿𝑝 𝑔𝑙 (𝑡) − 𝛿𝑝 𝑒𝑥 (𝑡) + 𝛾 𝑔𝑙 𝑘 𝐵 𝑇𝛿𝑂 𝑒𝑥 (𝑡))| 𝑡=𝑇 𝑠𝑡𝑖 . (4.53) In the next step, we estimate the radial velocity scale of the glial compartment 𝑢 𝑔𝑙𝑟∗ and radial velocity scale of the extracellular space 𝑢 𝑒𝑥𝑟∗ . By the incompressible condition, we have ▽⋅ (𝜂 𝑔𝑙 𝒖 𝑔𝑙 ) +▽⋅ (𝜂 𝑒𝑥 𝒖 𝑒𝑥 ) + ∂(𝜂 𝑎𝑥 𝑢 𝑎𝑥𝑧 )𝜕𝑧 = 0. It can be approximated by 𝜕(𝜂 𝑔𝑙 𝑢 𝑔𝑙𝑟 )𝜕𝑟 + 𝜕(𝜂 𝑒𝑥 𝑢 𝑒𝑥𝑟 )𝜕𝑟 = 0 , (4.54) because the propagation of the axon signal is only in the longitudinal direction and the water flow and osmotic difference are in the radial direction. Thanks to the boundary conditions at 𝑟 =0 , 𝑢 𝑔𝑙𝑟 = 𝑢 𝑒𝑥𝑟 = 0 , we have 𝜂 𝑔𝑙 𝑢 𝑔𝑙𝑟 + 𝜂 𝑒𝑥 𝑢 𝑒𝑥𝑟 = 0 . (4.55) Then the glial compartment velocity 𝑢 𝑔𝑙𝑟 can be rewritten as 𝑢 𝑔𝑙𝑟 = (1 − 𝛼)𝑢 𝑔𝑙𝑟 − 𝛼 𝜂 𝑒𝑥 𝜂 𝑔𝑙 𝑢 𝑒𝑥𝑟 , (4.56) where 𝛼 are defined as 9 𝛼 = 𝜅 𝑔𝑙 𝜏 𝑔𝑙𝜂𝑒𝑥𝜂𝑔𝑙 𝜅 𝑒𝑥 𝜏 𝑒𝑥 +𝜅 𝑔𝑙 𝜏 𝑔𝑙 . Based on eq. (4.56), eq. (2.15) and eq. (2.23), we estimate the scale of water velocity in the glial compartment along radius direction at 𝑡 = 𝑇 𝑠𝑡𝑖 as 𝑢 𝑔𝑙𝑟∗ = |(1 − 𝛼) 𝜅 𝑔𝑙 𝜏 𝑔𝑙 𝜇 𝛿𝑝 𝑔𝑙 −𝛿𝑝 𝑒𝑥 𝑟 ∗ − (1 − 𝛼) 𝜅 𝑔𝑙 𝜏 𝑔𝑙 𝜇 𝛾 𝑔𝑙 𝑘 𝐵 𝑇 𝛿𝑂 𝑔𝑙 𝑟 ∗ − 𝛼 𝜂 𝑒𝑥 𝜂 𝑔𝑙 𝑘 𝑒 𝜏 𝑒𝑥 𝛿𝜙 𝑒𝑥 𝑟 ∗ | 𝑡=𝑇 𝑠𝑡𝑖 . (4.57) We claim that the 𝛿𝑂 𝑔𝑙 in the stimulated region is due to change in the volume fraction of glial compartment 𝛿𝜂 𝑔𝑙 (see Remark 4.5) estimated as 𝛿𝑂 𝑔𝑙 ≈ 𝜂 𝑔𝑙𝑟𝑒 𝑂 𝑔𝑙𝑟𝑒 𝜂 𝑔𝑙𝑟𝑒 +𝛿𝜂 𝑔𝑙 − 𝑂 𝑔𝑙𝑟𝑒 = − 𝛿𝜂 𝑔𝑙 𝜂 𝑔𝑙𝑟𝑒 +𝛿𝜂 𝑔𝑙 𝑂 𝑔𝑙𝑟𝑒 . where 𝛿𝜂 𝑔𝑙 = 𝛿𝑝 𝑔𝑙 −𝛿𝑝 𝑒𝑥 𝐾 𝑔𝑙 . In eq. (4.57), the 𝛿𝜙 𝑒𝑥 estimation is given by eq. (4.34). Furthermore, by eq. (4.54), the scale of radial direction extracellular region velocity scale given by 𝑢 𝑒𝑥𝑟∗ = 𝜂 𝑔𝑙 𝜂 𝑒𝑥 𝑢 𝑔𝑙𝑟∗ . (4.58) In the eq. (4.52), the hydrostatic pressure is a passive reaction to the reduced osmotic concentration in the extracellular space. Fig. 4.3b shows that the water flow exhibits circulation patterns between the extracellular space and glial compartment. The water flow in the glial compartment is from the stimulated region to the non-stimulated region in the radial direction. In extracellular space, the water flow in the radial direction is from the non-stimulated region to stimulated region. Remark 4.5:
We assume the average total number of molecules (not concentration) in the stimulated glial region does not change since the major glial transmembrane ion flux in the stimulated region is K + flux and this K + flux from the stimulated extracellular space moves through the glial transition 𝑆 𝑡 to the non-stimulated extracellular space as eq. (4.36). 0 Figure 4.3: a: Schematic graph of the potassium flux when inner part axon stimulated. In the stimulated region, the potassium takes the way of extracellular pathway and through the glial compartment via glial membrane. In the non- stimulated region, the potassium leaks out to the extracellular space through the glial membrane. b: Schematic graph of the water circulation when inner part axon stimulated. In the stimulated region, the glial transmembrane water flow goes from extracellular space into glial compartment as the effect of osmosis difference. In the extracellular space, water goes from non-stimulated region to stimulated region in radial direction. In the glia compartment goes in the opposite direction. This compartment drawing is given only to aid qualitative understanding. Our model is in three dimensions. We do not compute with a compartment model, only with the full model described by partial differential equations in space and time. relative importance of flux components
In this section, we discuss the relative importance of the components of ion flux, diffusion, convection, and electric drift, in the glial and extracellular regions, respectively. The major flux is in the radial direction in each region and we restrict our discussion to that. In the extracellular space, the scales of convection term and electric drift term for ion 𝑖 =Na + , K + in the radial direction are described by 𝑐 𝑒𝑥𝑖 𝑢 𝑒𝑥𝑟,∗ and 𝑐 𝑒𝑥𝑖 | 𝛿𝜙 𝑒𝑥 𝑟 ∗ | , where 𝛿𝜙 𝑒𝑥 is estimated from eq. (4.34). The comparable radial diffusion flux term in the extracellular space is |𝐷 𝑒𝑥∗ 𝜏 𝑒𝑥 𝛿𝑐 𝑒𝑥𝑖 𝑟 ∗ | , 𝑖 = Na + , K + . We characterize the relative importance of electric drift and diffusion (of potassium and sodium) in the extracellular space by the ratios 𝑅 𝑒𝑥𝐾 and 𝑅 𝑒𝑥𝑁𝑎 analyzed in eq (4.42) and eq. (4.43) 𝑅 𝑒𝑥𝐾 = | 𝜂 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝜂 𝑒𝑥 𝜏 𝑒𝑥 𝜎 𝑒𝑥 (1+ℎ 𝜖 ) | , 𝑅 𝑒𝑥𝑁𝑎 = | 𝜂 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝜂 𝑒𝑥 𝜏 𝑒𝑥 𝜎 𝑒𝑥 (1+ℎ 𝜖 ) 𝑐 𝑒𝑥𝑁𝑎 𝑐 𝑒𝑥𝐾 | . For radial direction flux, the scale ratio between convection and diffusion in the extracellular space is estimated by the Peclet number 1 𝑃𝑒 𝑒𝑥𝑖 = | 𝑐 𝑒𝑥𝑖 𝑢 𝑒𝑥𝑟∗ 𝑟 ∗ 𝐷 𝑒𝑥𝑖 𝜏 𝑒𝑥 𝛿𝑐 𝑒𝑥𝑖 | , 𝑖 = Na + , K + . (4.59) Note that the Peclet number for Na + and K + are significantly different since their concentrations are different. In a similar way, we can estimate the Peclet number of the glial compartment. 𝑃𝑒 𝑔𝑙𝑖 = | 𝑐 𝑔𝑙𝑖 𝑢 𝑔𝑙𝑟∗ 𝑟 ∗ 𝐷 𝑔𝑙∗ 𝜏 𝑔𝑙 𝛿𝑐 𝑔𝑙𝑖 | , 𝑖 = Na + , K + . (4.60) The scale ratio between electric drift and diffusion in the glial compartment is 𝑅 𝑔𝑙𝐾 = | 𝜖 𝑐 𝑔𝑙𝐾 𝛿𝑐 𝑒𝑥𝐾 𝑐 𝑒𝑥𝐾 𝛿𝑐 𝑔𝑙𝐾 |, 𝑅 𝑔𝑙𝑁𝑎 = | 𝜖 𝑐 𝑔𝑙𝑁𝑎 𝛿𝑐 𝑒𝑥𝐾 𝑐 𝑒𝑥𝐾 𝛿𝑐 𝑔𝑙𝑁𝑎 | , (4.61) where we have used eq. (4.32) and eq. (4.21). We estimate the potassium concentration changes after the stimulus 𝑡 = 𝑇 𝑠𝑡𝑖 in the stimulated glial compartment is 𝛿𝑐 𝑔𝑙𝐾 , 𝛿𝑐 𝑔𝑙𝐾 ≈ (𝑛𝑐 𝑠𝑡𝑖 − 𝛿𝑐 𝑒𝑥𝐾 ) 𝜆 𝑔𝑙𝑚,𝐾 𝜆 𝑔𝑙𝑚,𝐾 +𝜆 𝑒𝑥𝐾 𝜂 𝑒𝑥 𝜂 𝑔𝑙 , (4.62) We estimate 𝛿𝑐 𝑔𝑙𝑁𝑎 in the stimulated glial compartment as 𝛿𝑐 𝑔𝑙𝑁𝑎 ≈ − 𝑎𝑥 𝑔 𝑔𝑙𝐾 (𝛿𝑉 𝑔𝑙 −𝛿𝐸 𝑔𝑙𝐾 ) 𝛿𝑐 𝑔𝑙𝐾 , (4.63) where 𝛿𝐼 𝑎𝑥 and 𝑔 𝑔𝑙𝐾 (𝛿𝑉 𝑔𝑙 − 𝛿𝐸 𝑔𝑙𝐾 ) are derived as 𝛿𝐼 𝑎𝑥 ≈ 2 ( K K1 𝐼 𝑎𝑥𝑟𝑒,1 𝑐 𝑒𝑥𝐾,𝑒𝑞 (𝑐 𝑒𝑥𝐾,𝑟𝑒 +𝐾 𝐾1 ) + K K2 𝐼 𝑎𝑥𝑟𝑒,2 𝑐 𝑒𝑥𝐾,𝑒𝑞 (𝐶 𝑒𝑥𝐾,𝑟𝑒 +𝐾 𝐾2 ) ) 𝛿𝑐 𝑒𝑥𝐾 , 𝑔 𝑔𝑙𝐾 (𝛿𝑉 𝑔𝑙 − 𝛿𝐸 𝑔𝑙𝐾 ) = − 𝑔 𝑔𝑙𝐾 ℎ 𝜖 𝜖 𝛿𝐸 𝑔𝑙𝐾 . In the next section, we carry out a numeric simulation as mentioned previously. Furthermore, we compare the results between the electrodiffusion model with the convection-electrodiffusion (full) model.
In this section, we applied a train of current pulses to stimulate the axon membrane near the left boundary ( {(𝑧 , 𝑟)|𝑧 = 1.5 mm , 𝑟 < 𝑟 𝑠𝑡𝑖 = 𝑟 ∗ = 24 μm}) . Each single stimulus has current strength 𝐼 𝑠𝑡𝑖 = 3 × 10 −3 A/m and duration . The frequency is 50 Hz (𝑇 = 0.02 s) and total stimulus time 𝑇 𝑠𝑡𝑖 = 0.2 s . We first estimate how large are the fluid velocities of extracellular space and glial compartment can be generated by a train of stimuli. From eq. (4.44) and eq. (4.45), we estimate the concentration variation in the stimulated extracellular region at 𝑡 = 𝑇 𝑠𝑡𝑖 as 2 𝛿𝑐 𝑒𝑥𝑁𝑎 ≈ −1.06 mM, 𝛿𝑐 𝑒𝑥𝐾 ≈ 0.89 mM, 𝛿𝑂 𝑒𝑥 ≈ − 0.34 mM . The estimated glial transmembrane velocity in eq. (4.53) is 𝑈 𝑔𝑙𝑚∗ ≈ 9.78 × 10 −2 nm/s . From (4.57) and (4.58), the estimated scale of radial velocities of glial compartment and extracellular space are 𝑢 𝑒𝑥𝑟,∗ ≈ 1.56 × 10 nm/s, 𝑢 𝑔𝑙𝑟,∗ ≈ 3.90 nm/s . In Fig. 4.4a-c, we plot the computed average variation of concentrations in the stimulated extracellular region. These computed concentration changes are consistent with the estimates presented previously. The change of concentration reaches its peak at the end of the train of stimulus (𝑡 = 𝑇 𝑠𝑡𝑖 ) and quickly returns to its previous equilibrium value. In Fig. 4.4f, we plot the computed average transmembrane water flow through the glial membrane in the stimulated region. We see Fig. 4.3b that water flows into the glial compartment from the extracellular space in the stimulated region. This transmembrane water flow generates the water circulation between the stimulated region and non-stimulated region in the radial direction. As in the Fig. 4.3b, in the extracellular compartment, the water flow goes from the non-stimulated region to the stimulated region and in the glial compartment, water flows in the opposite (radial) direction. In the Fig. 4.4d-e, we plot the computed average water velocity in the radial direction in the glial compartment and in the extracellular space. The computations are consistent with our estimation above. Figure 4.4: a-c: average variation of concentration in the stimulated extracellular space . d-e: the average intradomain radial velocity (radial direction is the normal direction). f: average glial transmembrane velocity in the stimulated region (the normal direction points from glial compartment to extracellular space.) In the Fig. 4.5a, we show the transmembrane water flow through the glial membrane in the non-stimulated region as in the Schematic Fig. 4.3b. This water flow to the extracellular space produces widening of the extracellular space volume in the non-stimulated region, as shown in 3 Fig. 4.5b. At the same time, the extracellular space volume shrinks (in the stimulated region) as shown in Fig. 4.5c. The shrinkage is produced by the inward water flow through the glial membrane in stimulated region, as in Fig. 4.5f. Our simulation is consistent with the experiments in references [84, 85], where the extracellular space become smaller in the middle cortical layers (where the stimulus is applied) but widens in the most superficial and deep cortical layers (where no stimulus is applied).
Figure 4.5: a: average glial transmembrane velocity in the non-stimulated region (the normal direction points from glial compartment to extracellular space.). b-c: average variation of the extracellular volume fraction in non- stimulated region and stimulated regions.
In this section, we try to investigate the importance of the fluid convection during the potassium clearance. We examine the estimated Peclet number for Na + and K + in the extracellular and glia compartments. By eq. (4.59), the Peclet number in the extracellular space (for radial flux) is 𝑃𝑒 𝑒𝑥𝐾 = | 𝑐 𝑒𝑥𝐾 𝑢 𝑒𝑥𝑟∗ 𝑟 ∗ 𝐷 𝑒𝑥𝐾 𝜏 𝑒𝑥 𝛿𝑐 𝑒𝑥𝐾 | ≈ 1.0 × 10 −2 , 𝑃𝑒 𝑒𝑥𝑁𝑎 = | 𝑐 𝑒𝑥𝑁𝑎 𝑢 𝑒𝑥𝑟∗ 𝑟 ∗ 𝐷 𝑒𝑥𝑁𝑎 𝜏 𝑒𝑥 𝛿𝑐 𝑒𝑥𝑁𝑎 | ≈ 3.5 × 10 −1 . By eq. (4.42) and eq. (4.43), the scale ratio between electric drift and diffusion in the extracellular space (for radial flux) are 𝑅 𝑒𝑥𝐾 ≈ 6.2 × 10 −2 , 𝑅 𝑒𝑥𝑁𝑎 ≈ 2.3. In the glial compartment, based on eq. (4.59), eq. (4.61) and eq. (4.62), we have 𝑃𝑒 𝑔𝑙𝐾 = | 𝑐 𝑔𝑙𝐾 𝑢 𝑔𝑙𝑟∗ 𝑟 ∗ 𝐷 𝑔𝑙𝐾 𝜏 𝑔𝑙 𝛿𝑐 𝑔𝑙𝐾 | ≈ 2.9 × 10 , 𝑃𝑒 𝑔𝑙𝑁𝑎 = | 𝑐 𝑔𝑙𝑁𝑎 𝑢 𝑔𝑙𝑟∗ 𝑟 ∗ 𝐷 𝑔𝑙𝑁𝑎 𝜏 𝑔𝑙 𝛿𝑐 𝑔𝑙𝑁𝑎 | ≈ 1.7 × 10 . By eq. (4.61), the scale ratio between electric drift and diffusion in the glial space in the radial direction is 𝑅 𝑔𝑙𝐾 = | 𝜖 𝑐 𝑔𝑙𝐾 𝛿𝑐 𝑒𝑥𝐾 𝑐 𝑒𝑥𝐾 𝛿𝑐 𝑔𝑙𝐾 | ≈ 4.3 × 10 , 𝑅 𝑔𝑙𝑁𝑎 = | 𝜖 𝑐 𝑔𝑙𝑁𝑎 𝛿𝑐 𝑒𝑥𝐾 𝑐 𝑒𝑥𝐾 𝛿𝑐 𝑔𝑙𝑁𝑎 | ≈ 1.7 × 10 . 4 In Fig. 4.6, we plot the computed potassium and sodium fluxes in the extracellular and glial compartments (in the radial direction) Figure 4.6: a: average radial direction fluxes components in the extracellular space. b: average radial direction fluxes components in the glial compartment (radial direction as normal direction).
In the extracellular space, estimations of the importance of different fluxes are more complicated because they depend on the species of ion. Diffusion flux of potassium (Fig. 4.6 a upper panel), dominant, but not for the sodium (Fig. 4.6 a lower panel). The three fluxes, diffusion, convection and electric drift, are comparable to the larger electric drift flux. These calculations agree with our estimations above. In the extracellular space, potassium’s
Pelect number 𝑃𝑒 𝑒𝑥𝐾 and the ratio 𝑅 𝑒𝑥𝐾 are in O(10 −2 ) , while sodium’s Pelect number 𝑃𝑒 𝑒𝑥𝑁𝑎 is order of O(10 −1 ) and the ratio 𝑅 𝑒𝑥𝑁𝑎 is in O(1) . In the glial compartments (Fig. 4.6b), the situation is different form the extracellular space. The electric drift is the prime dominant one and convection flux come as second dominant for both sodium and potassium (see Fig. 4.6). The water flow has more important effect for potassium in the glial compartment than in the extracellular space. The maximum of the convection flux occurs after the stimuli since it takes that long for osmotic pressure to accumulate. Also, it lasts longer when the effect of electric drift diminished. In the Fig. 4.7a and 4.7b, the potassium and sodium flux through the glial membrane are presented and the results are consistent with our estimates. The major current through the glial membrane is through the potassium channel in both stimulated region and non-stimulated region. Fig. 4.7c compares the stimulated and non-stimulated region. Fig. 4.7c shows the total potassium flux through potassium channels (integrated over all the glial membrane). The total potassium flux has different direction in the stimulated region and non-stimulated region, as shown in our estimation in eq. (4.36). The strength is the same, but the direction is different. 5
Figure 4.7: a: potassium and sodium flux variation through Na/K pump and ion channels on the glial membrane in the stimulated region. b: potassium and sodium flux variation through Na/K pump and ion channels on the glial membrane in the non-stimulated region. c: the total potassium flux through potassium channel on the glial membrane.
Fig. 4.8 compares the potassium flux in the electrodiffusion (ED) model and convection-electrodiffusion (full) model. In the full model, the water circulation between the stimulated and non-stimulated region in both extracellular and glial compartments have an important role in the circulation of potassium. The water circulation has an important role in buffering potassium in the optic nerve bundle. The water circulation increases the potassium flow through the glial compartment. Fig. 4.8f show how water flow increases the potassium flux through the glia in the transition region between the stimulated and non-stimulated region. The potassium flux moves back to the stimulated extracellular region from non-stimulated extracellular region through the extracellular pathway, as shown in Fig. 4.8b. The time rate of change of the cumulative K + flux through the extracellular transition region decreases after stimulus. Figure 4.8: a-d: comparison between electrodiffusion model and full model on average K + flux, cumulative K + flux and K + flux components on extracellular transition region 𝑆 𝑡 . e-h: the same comparison on the glial transition region 𝑆 𝑡 (radial direction as normal direction). Multiple trains of action potentials strengthen the effect of water flow on the transport through the glial compartment. In the Fig. 4.9, three trains of action potentials occur with 0.2 s resting 6 period between each. Fig. 4.9f shows that water flow increases the amount of cumulative potassium flux through the transition region in the glial compartment, beyond the potassium flow in the electro-diffusion model. Consequently, the amount of cumulative potassium flux through the transition region in the extracellular space is less than in the electro-diffusion model see Fig. 4.9b. Figure 4.9: comparison between electrodiffusion model and full model when three trains of the stimulus current applied on the axon (radial direction as normal direction.).
5. Conclusion
This work provides a comprehensive set of estimates and computations, showing the water circulation in the optic nerve. The water flow is generated by the osmotic difference between the glial compartment and extracellular space. Through the estimation, we show that in the stimulated region, the extracellular osmotic changes are not induced by ion fluxes from the axon compartment when the axon is firing. Indeed, based on the analysis, we found that the leading order of potassium flux out and sodium flux into axon is the same during the action potential, which is consistent with the literature (64,83,84). The osmotic difference is generated due to the sodium and potassium conductance difference in the glial membrane. In other words, more potassium leaks into the glial compartment, and less sodium leaks out. As a result of this glial transmembrane water flow in the stimulated region, it forms a water circulation in the radial direction between the stimulated region and the non-stimulated region. Our estimation of the velocity scales in the glial compartment and extracellular space shows that this water flow has a considerable effect on potassium flux in the glial compartment. By comparing the full model (including water) with the electrodiffusion model (exclude water), we validate that water circulation through the glial pathway helps clearance of potassium in the extracellular space and enhance the glial buffering effect. With additional numerical simulations, we show that the repetitive activity of the nerve fibers further increases the importance of water flow, and the water flow contribution to glia buffering, which is likely to dramatically dominate pathological situations of repetitive activity. 7 Besides, through our analysis, we show that the electrical syncytium property of the glial cells is critical for clearing potassium (from the extracellular space) when the neuron fires. Based on the governing equation of glial electric potential, we explain why the inward glial transmembrane potassium flux in the stimulated region is almost the same as the outward potassium flux out to the extracellular space in the non-stimulated region when axon firing. This is because the electric potential spreads through the connected cells in the glial compartment. The glial electric potential in the non-stimulated region becomes more positive in response to the depolarization of the glial electric potential in the stimulated region. This electric property for the glial compartment is always exist as long as there exists two distinguish stimulated region and non-stimulated region. The glial wrap the axon like a faster potassium transporter, which quickly remove the extra potassium (in the extracellular space) from the stimulated region to the non-stimulated region.
Finally, we’d like to point out that the coupling of ionic and water flows is not unique to optical nerve. It is ubiquitous in many parts of the mammalian body and other biological tissues. Our analysis of the model for the optical nerve is just a first small step towards the understanding of the mechanisms of various transport processes and the consequences of a disrupted process under pathological conditions.
Reference
1. Nicholson , C. and S. Hrabětová,
Brain Extracellular Space: The Final Frontier of Neuroscience.
Biophysical Journal, 2017. (10): p. 2133-2142. 2. Nedergaard, M. and S.A. Goldman,
Glymphatic failure as a final common pathway to dementia.
Science, 2020. (6512): p. 50-56. 3. Orkand, R.K., J.G. Nicholls, and S.W. Kuffler,
Effect of nerve impulses on the membrane potential of glial cells in the central nervous system of amphibia.
J Neurophysiol, 1966. (4): p. 788-806. 4. Frankenhaeuser, B. and A.L. Hodgkin, The after-effects of impulses in the giant nerve fibres of Loligo.
Journal of Physiology (London), 1956. : p. 341-376. 5. Boron, W. and E. Boulpaep,
Medical Physiology . 2008, New York: Saunders. 1352. 6. Davson, H.,
A textbook of general physiology, 4th Edition . 4th ed. 1970, New York: Churchill. 1694. 7. Ruch, T.C. and H.D. Patton,
Physiology and Biophysics, Volume 2: Circulation, Respiration and Balance . Vol. 2. 1973, Philadelphia: W.B. Saunders Company. 558. 8. Ray, L., J.J. Iliff, and J.J. Heys,
Analysis of convective and diffusive transport in the brain interstitium.
Fluids and Barriers of the CNS, 2019. (1): p. 6. 9. Tuckwell, H.C. and R.M. Miura, A mathematical model for spreading cortical depression.
Biophysical Journal, 1978. (2): p. 257-276. 10. Postnov, D.E., L.S. Ryazanova, and O.V. Sosnovtseva, Functional modeling of neural – glial interaction. BioSystems, 2007. (1-3): p. 84-91. 11. Fu, K.Y., et al., Sciatic nerve regeneration by microporous nerve conduits seeded with glial cell line ‐ derived neurotrophic factor or brain ‐ derived neurotrophic factor gene transfected neural stem cells. Artificial organs, 2011. (4): p. 363-372. 12. Chang, J.C., et al., A mathematical model of the metabolic and perfusion effects on cortical spreading depression.
PLoS One, 2013. (8): p. e70469. 13. Yao, W., H. Huang, and R.M. Miura, A continuum neuronal model for the instigation and propagation of cortical spreading depression.
Bulletin of mathematical biology, 2011. (11): p. 2773-2790. 14. Gardiner, B.S., et al., Computational modeling of fluid flow and intra-ocular pressure following glaucoma surgery.
PloS one, 2010. (10). 15. Norman, R.E., et al., Finite element modeling of the human sclera: influence on optic nerve head biomechanics and connections with glaucoma.
Experimental eye research, 2011. (1): p. 4-12. 16. Salazar, J.J., et al., Anatomy of the Human Optic Nerve: Structure and Function , in
Optic Nerve . 2018, IntechOpen. 17. Selhorst, J.B. and Y. Chen.
The optic nerve . in
Seminars in neurology . 2009. © Thieme Medical Publishers. 18. Hayreh, S.S.,
Ischemic optic neuropathy.
Progress in retinal and eye research, 2009. (1): p. 34-62. 19. Atchison, D. and G. Smith, Optics of the human eye . 2000: Butterworth Heinemann. 20. Kuffler, S.W., J.G. Nicholls, and R.K. Orkand,
Physiological properties of glial cells in the central nervous system of amphibia.
J Neurophysiol, 1966. (4): p. 768-87. 21. Luo, Y., Angelo R. Rossi, and Andrew L. Harris, Computational Studies of Molecular Permeation through Connexin26 Channels.
Biophysical Journal. (3): p. 584-599. 22. Harris, Andrew L. and I.D. Locke, eds.
Connexins: a guide . 2008, Springer. 573. 23. Harris, A.L.,
Connexin channel permeability to cytoplasmic molecules.
Progress in biophysics and molecular biology, 2007. (1-2): p. 120-143. 24. Simpson, I., B. Rose, and W.R. Loewenstein, Size limit of molecules permeating the junctional membrane channels.
Science, 1977. (4275): p. 294-296. 25. Bellot-Saez, A., et al.,
Astrocytic modulation of neuronal excitability through K+ spatial buffering.
Neuroscience & Biobehavioral Reviews, 2017. : p. 87-97. 26. Wallraff, A., et al., The impact of astrocytic gap junctional coupling on potassium buffering in the hippocampus.
Journal of Neuroscience, 2006. (20): p. 5438-5447. 27. Katz, B., Les constantes electrique de la membrane du muscle.
Arch.Sci.physiol., 1949. ( ): p. 285-300. 28. Sibille, J., et al., The neuroglial potassium cycle during neurotransmission: role of Kir4. 1 channels.
PLoS computational biology, 2015. (3): p. e1004137. 29. Kubo, Y., et al., Primary structure and functional expression of a mouse inward rectifier potassium channel.
Nature, 1993. (6416): p. 127-133. 30. Hayreh, S.S.,
The sheath of the optic nerve.
Ophthalmologica, 1984. (1-2): p. 54-63. 31. Killer, H., et al.,
Architecture of arachnoid trabeculae, pillars, and septa in the subarachnoid space of the human optic nerve: anatomy and clinical considerations.
British Journal of Ophthalmology, 2003. (6): p. 777-781. 9 32. Pache, M. and P. Meyer, Morphological changes of the retrobulbar optic nerve and its meningeal sheaths in glaucoma.
Ophthalmologica, 2006. (6): p. 393-396. 33. Hua, Y., A.P. Voorhees, and I.A. Sigal,
Cerebrospinal fluid pressure: revisiting factors influencing optic nerve head biomechanics.
Investigative ophthalmology & visual science, 2018. (1): p. 154-165. 34. Andres, K., et al., Nerve fibres and their terminals of the dura mater encephali of the rat.
Anatomy and embryology, 1987. (3): p. 289-301. 35. Killer, H.E., H.R. Laeng, and P. Groscurth,
Lymphatic capillaries in the meninges of the human optic nerve.
Journal of Neuro Ophthalmology, 1999. (4): p. 222-228. 36. Morgan, W.H., et al., Cerebrospinal fluid pressure and the eye.
British Journal of Ophthalmology, 2016. (1): p. 71-77. 37. Filippidis, A.S., et al.,
Permeability of the arachnoid and pia mater. The role of ion channels in the leptomeningeal physiology.
Child's Nervous System, 2012. (4): p. 533-540. 38. Vogiatzidis, K., et al., μ -Opioid influence on transmesothelial resistance of isolated sheep pleura and parietal pericardium. European journal of pharmacology, 2006. (3): p. 276-280. 39. Hatzoglou, C., K. Gourgoulianis, and P. Molyvdas,
Effects of SNP, ouabain, and amiloride on electrical potential profile of isolated sheep pleura.
Journal of Applied Physiology, 2001. (4): p. 1565-1569. 40. Payne, D.K., G.T. Kinasewitz, and E. Gonzalez, Comparative permeability of canine visceral and parietal pleura.
Journal of Applied Physiology, 1988. (6): p. 2558-2564. 41. Sarkos, S., et al., Effect of amiloride in human and sheep parietal pleura.
Respiratory physiology & neurobiology, 2002. (2): p. 233-237. 42. Zarogiannis, S., et al.,
Comparison of the electrophysiological properties of the sheep isolated costal and diaphragmatic parietal pleura.
Clinical and experimental pharmacology & physiology, 2007. (1-2): p. 129-131. 43. Zarogiannis, S., et al., Adrenergic influence on the permeability of sheep diaphragmatic parietal pleura.
Respiration, 2007. (1): p. 118-120. 44. Zarogiannis, S., et al., Dexamethasone decreases the transmesothelial electrical resistance of the parietal and visceral pleura.
The Journal of Physiological Sciences, 2009. (4): p. 335-339. 45. Li, F.K., et al., Electrophysiology and glucose transport of human peritoneal mesothelial cells: implications for peritoneal dialysis.
Peritoneal dialysis international, 2001. (2): p. 115-121. 46. Simon, M., Peritoneal mesothelium in vitro: an electrophysiologic study.
Peritoneal dialysis international, 1996. (4): p. 393-397. 47. Stefanidis, I., et al., Amiloride-sensitive sodium channels on the parietal human peritoneum: evidence by ussing-type chamber experiments.
Asaio Journal, 2007. (3): p. 335-338. 48. Zarogiannis, S., et al. Influence of the sodium transport inhibition by amiloride on the transmesothelial resistance of isolated visceral sheep peritoneum . in
Advances in peritoneal dialysis. Conference on Peritoneal Dialysis . 2005. 49. Zarogiannis, S., et al.,
Effect of sodium-potassium pump inhibition by ouabain on the permeability of isolated visceral sheep peritoneum.
Adv Perit Dial, 2007. : p. 43-47. 0 50. Zarogiannis, S., et al., µ-opioid stimulation of isolated parietal sheep peritoneum decreases peritoneal permeability in vitro. Adv. Perit. Dial, 2007. : p. 34-37. 51. Verikouki, C., et al., Rapid effect of progesterone on transepithelial resistance of human fetal membranes: evidence for non ‐ genomic action. Clinical and Experimental Pharmacology and Physiology, 2008. (2): p. 174-179. 52. Adams, E.A., et al., Comparison of amniotic and intramembranous unidirectional permeabilities in late-gestation sheep.
American journal of obstetrics and gynecology, 2005. (1): p. 247-255. 53. Filippidis, A., et al.,
Transmembrane resistance and histology of isolated sheep leptomeninges.
Neurological research, 2010. (2): p. 205-208. 54. Causin, P., et al., A poroelastic model for the perfusion of the lamina cribrosa in the optic nerve head.
Mathematical biosciences, 2014. : p. 33-41. 55. Jansonius, N.M., et al.,
A mathematical model for describing the retinal nerve fiber bundle trajectories in the human eye: average course, variability, and influence of refraction, optic disc size and optic disc position.
Experimental eye research, 2012. : p. 70-78. 56. Sigal, I.A., et al.,
Finite element modeling of optic nerve head biomechanics.
Investigative ophthalmology & visual science, 2004. (12): p. 4378-4387. 57. Band, L.R., et al., Intracellular flow in optic nerve axons: a mechanism for cell death in glaucoma.
Investigative ophthalmology & visual science, 2009. (8): p. 3750-3758. 58. McNeal, D.R., Analysis of a model for excitation of myelinated nerve.
IEEE Transactions on Biomedical Engineering, 1976(4): p. 329-337. 59. Oozeer, M., et al.,
Simulation of intra-orbital optic nerve electrical stimulation.
Medical and Biological Engineering and Computing, 2005. (5): p. 608-617. 60. Bellinger, S., G. Miyazawa, and P. Steinmetz, Submyelin potassium accumulation may functionally block subsets of local axons during deep brain stimulation: a modeling study.
Journal of neural engineering, 2008. (3): p. 263. 61. Chen, K.C. and C. Nicholson, Spatial buffering of potassium ions in brain extracellular space.
Biophysical journal, 2000. (6): p. 2776-2797. 62. Reichenbach, A., et al., What do retinal Müller (glial) cells do for their neuronal ‘small siblings’?
Journal of chemical neuroanatomy, 1993. (4): p. 201-213. 63. Ø stby, I., et al., Astrocytic mechanisms explaining neural-activity-induced shrinkage of extraneuronal space.
PLoS computational biology, 2009. (1): p. e1000272. 64. Hou, R., et al., Intracranial pressure (ICP) and optic nerve subarachnoid space pressure (ONSP) correlation in the optic nerve chamber: the Beijing Intracranial and Intraocular Pressure (iCOP) study. brain research, 2016. : p. 201-208. 65. Killer, H., et al.,
Cerebrospinal fluid dynamics between the intracranial and the subarachnoid space of the optic nerve. Is it always bidirectional?
Brain, 2007. (2): p. 514-520. 66. Mori, Y.,
A multidomain model for ionic electrodiffusion and osmosis with an application to cortical spreading depression.
Physica D: Nonlinear Phenomena, 2015. (15): p. 94-108. 67. Zhu, Y., et al.,
A Tridomain Model for Potassium Clearance in Optic Nerve. arXiv preprint arXiv:2012.03303, 2020. 68. Nicholson, C.,
Diffusion and related transport mechanisms in brain tissue.
Reports on 1 progress in Physics, 2001. (7): p. 815. 69. Xu, S., et al., Osmosis through a Semi-permeable Membrane: a Consistent Approach to Interactions. arXiv preprint arXiv:1806.00646, 2018. 70. Zhu, Y., et al.,
A Bidomain Model for Lens Microcirculation
Biophysical Journal, 2019. (6): p. 1171-1184 Preprint available at https://arxiv.org/abs/1810.04162. 71. Pérez-Pinzón, M., L. Tao, and C. Nicholson,
Extracellular potassium, volume fraction, and tortuosity in rat hippocampal CA1, CA3, and cortical slices during ischemia.
Journal of Neurophysiology, 1995. (2): p. 565-573. 72. McLaughlin, S. and R.T. Mathias, Electro-osmosis and the reabsorption of fluid in renal proximal tubules.
J Gen Physiol, 1985. (5): p. 699-728. 73. Vaghefi, E., et al., Development of a 3D finite element model of lens microcirculation.
Biomed Eng Online, 2012. : p. 69. 74. Wan, L., et al., Self-Consistent Approach to Global Charge Neutrality in Electrokinetics: A Surface Potential Trap Model.
Physical Review X, 2014. (1): p. 011042. 75. Ethier, C., et al., Finite Element Modeling of the Human Sclera: Influence on ONH Biomechanics and Connections With Glaucoma.
Investigative Ophthalmology & Visual Science, 2009. (13): p. 4889-4889. 76. Hodgkin, A.L. and A.F. Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve.
J. Physiol., 1952. : p. 500-544. 77. Fitzhugh, R.,
Thresholds and plateaus in the Hodgkin-Huxley nerve equations.
The Journal of general physiology, 1960. (5): p. 867-896. 78. Gao, J., et al., Isoform-specific function and distribution of Na/K pumps in the frog lens epithelium.
J Membr Biol, 2000. (2): p. 89-101. 79. Busija, D.W., et al.,
Mechanisms involved in the cerebrovascular dilator effects of cortical spreading depression.
Progress in neurobiology, 2008. (4): p. 417-433. 80. Smith, J.M., et al., Physiological studies of cortical spreading depression.
Biological Reviews, 2006. (4): p. 457-481. 81. Dreier, J.P. and C. Reiffurth, The stroke-migraine depolarization continuum.
Neuron, 2015. (4): p. 902-922. 82. Murakami, S. and Y. Kurachi, Mechanisms of astrocytic K+ clearance and swelling under high extracellular K+ concentrations.
The Journal of Physiological Sciences, 2016. (2): p. 127-142. 83. Sæ tra, M.J., G. Halnes, and G.T. Einevoll, An electrodiffusive neuron-extracellular-glia model with somatodendritic interactions. bioRxiv, 2020. 84. Holthoff, K. and O.W. Witte,
Directed spatial potassium redistribution in rat neocortex.
Glia, 2000. (3): p. 288-292. 85. Kofuji, P. and E.A. Newman, Potassium buffering in the central nervous system.
Neuroscience, 2004. (4): p. 1043-1054. 83. Dietzel, I., Heinemann, U., Hofmeier, G. and Lux, H.D., 1982. Stimulus-induced changes in extracellular Na+ and Cl− concentration in relation to changes in the size of the extracellular space. Experimental brain research, 46(1), pp.73-84. 84. Keynes, R.D., 1951. The ionic movements during nervous activity. The Journal of physiology, 114(1-2), p.119. 2 85. Susuki, K. (2010) Myelin: A Specialized Membrane for Cell Communication. Nature Education 3(9):59
Appendix: Supporting Information
A1. Notations 𝑐 𝑙𝑖 : Ion 𝑖 concentration in the 𝑙 region, 𝐾 𝑘 : Stiffness constant of k membrane, 𝜙 𝑙 : Electric potential in 𝑙 region, 𝜏 𝑙 : Tortuosity of 𝑙 region, 𝑝 𝑙 : Hydrostatic pressure in 𝑙 region, 𝑧 𝑖 : Valence of the ion 𝑖 , 𝒖 𝑙 : Fuild velocity inside of the 𝑙 region, 𝐴 𝑙 : Negative charged protein density in 𝑙 region, 𝜂 𝑙 : Volume fraction of 𝑙 region, 𝐽 𝑝,k𝑖 : Active ATP based ion 𝑖 pump on k membrane, 𝑂 𝑙 : Osmotic concentration in 𝑙 region, 𝐽 𝑐,k𝑖 : Passive transmembrane source of k membrane, ℳ 𝑘 : Membrane area k in per unit control volume Δ𝑉 , 𝑔 𝑘𝑖 : Conductance of 𝑘 membrane for ion i , 𝜅 𝑙 : Water permeability of 𝑙 region, 𝑔 𝑖 ̅ : Maximum conductance of axon membrane for ion i , 𝐿 𝑘𝑚 : Membrane hydrostatic permeability of k membrane, 𝑔 𝑙𝑒𝑎𝑘𝑖 : Leak conductance of axon membrane for ion 𝑖, 𝜇 : Fluid viscosity. A2. Comparison between membrane potential and Nernst potential on axon membrane
The classical Hodgkin Huxley analysis of a single action potential [77] assumes that changes in concentration of ions are much less important than current flow in determining the shape of the action potential. In other words, the change in the Nernst (i.e., equilibrium) potential is much less than the change in the membrane potential. In this part, we show that in our model the variation of the Nernst potential for Na + , K + and Cl − is much smaller than the variation of the axon membrane potential during axon firing, namely, 𝛿𝐸 𝑎𝑥𝑖 ≪ 𝛿V ax , 𝑖 = Na + , K + , Cl − . During a single action potential can be approximated by the difference between the resting state Na + and K + Nernst potential, 𝛿𝑉 𝑎𝑥∗ = 𝑂(𝐸 𝑎𝑥𝑁𝑎,𝑟𝑒 − 𝐸 𝑎𝑥𝐾,𝑟𝑒 ) ≫ 𝑂(𝑉 ∗ ), (A-1) where 𝑉 ∗ = 𝑘 𝐵 𝑇𝑒 . Taking the Cl − for example, by the charge neutrality condition, we have 𝛿𝑐 𝑎𝑥𝐶𝑙 ≈ − 𝜂 𝑒𝑥 𝜂 𝑎𝑥 𝛿𝑐 𝑒𝑥𝐶𝑙 . ( A-2)
By eq. (A-2), the Cl − Nernst potential on axon membrane yields 𝛿𝐸 𝑎𝑥𝐶𝑙 = 𝑉 ∗ (log ( 𝑐 𝑒𝑥𝐶𝑙,𝑟𝑒 +𝛿𝑐 𝑒𝑥𝐶𝑙 𝑐 𝑎𝑥𝐶𝑙,𝑟𝑒 + 𝛿𝑐 𝑎𝑥𝐶𝑙 ) − log ( 𝑐 𝑒𝑥𝐶𝑙,𝑟𝑒 𝑐 𝑎𝑥𝐶𝑙,𝑟𝑒 )) ≈ 𝑉 ∗ (log (1 + 𝛿𝑐 𝑒𝑥𝐶𝑙 𝑐 𝑒𝑥𝐶𝑙,𝑟𝑒 ) − log (1 − 𝜂 𝑒𝑥 𝛿𝑐 𝑒𝑥𝐶𝑙 𝜂 𝑎𝑥 𝑐 𝑎𝑥𝐶𝑙,𝑟𝑒 )) . ( A-3)
In eq. (A-3), we know 𝑐 𝑒𝑥𝐶𝑙,𝑟𝑒 ~10 mM and 𝜂 𝑎𝑥 𝑐 𝑎𝑥𝐶𝑙,𝑟𝑒 𝜂 𝑒𝑥 ~10 mM . We know the characteristic time for a single action potential 𝑇 𝑎𝑥∗ = 𝑂(10 −3 ) (millisecond) so the scale of 𝛿𝑐 𝑒𝑥𝐶𝑙 is 3 𝛿𝑐 𝑒𝑥𝐶𝑙,∗ = 𝑂(𝛿𝑐 𝑒𝑥𝑁𝑎 + 𝛿𝑐 𝑒𝑥𝐾 ) < 𝑂 ( 𝑇 𝑎𝑥∗ ℳ𝑎𝑥 𝑔̅ 𝑁𝑎 𝑉 𝑎𝑥∗ 𝑒𝜂 𝑒𝑥 ) = O (1). ( A-4) where we use charge neutrality condition and assume the voltage-gated ion channel fully open to its maximum conductance. Therefore, eq. (
A-3) yields 𝛿𝐸 𝑎𝑥𝐶𝑙 ≈ 𝑉 ∗ ( 𝑒𝑥𝐶𝑙,𝑟𝑒 + 𝜂 𝑒𝑥 𝜂 𝑎𝑥 𝑐 𝑎𝑥𝐶𝑙,𝑟𝑒 ) 𝛿𝑐 𝑒𝑥𝐶𝑙 . ( A-5)
Therefore, based on eq. (A-1) and eq. (A-5) , we have 𝛿𝐸 𝑎𝑥𝐶𝑙 ≪ 𝛿V ax . A similar analysis of Na + and K + gives 𝛿𝐸 𝑎𝑥𝐾,𝑁𝑎 ≪ 𝛿𝑉 𝑎𝑥 . A3. Robust analysis for the HH-model estimation
In the following Table 1-2, we show that the estimation of 𝑡 𝑚1 and 𝑡 𝑚2 are robust if different open probability of activation of Na + channel, 𝑚 𝑑𝑦 (𝑡 𝑚1 ) and 𝑚 𝑑𝑦 (𝑡 𝑚2 ), are chosen in eq. (4.7) and eq. (4.9) . In Table 1 below, the estimation of the 𝑡 𝑚1 is consistent when we choose different probabilities for 𝑚 𝑑𝑦 at time 𝑡 = 𝑡 𝑚1 . Table 1: Estimation of 𝒕 𝒎𝟏 with 𝒎 𝒅𝒚 (𝒕 𝒎𝟏 ) Table 2 uses 𝑚 𝑑𝑦 (𝑡 𝑚1 ) = 0.95 as the initial value 𝑚 in eq. (4.9) for the second period 𝑡 𝑚2 estimation. In a similar way, we choose different open probability of activation of Na + channel 𝑚 𝑑𝑦 at time 𝑡 = 𝑡 𝑚2 . In Table 2, it shows that the estimation of 𝑡 𝑚2 is consistent Table 2: Estimation of 𝒕 𝒎𝟐 with 𝒎 𝒅𝒚 (𝒕 𝒎𝟐 ) In sum, based on the results in Table 1-2, we confirm that by using eq. (4.7) and eq. (4.9) to estimate the slope of the linear approximation 𝛿𝑉 𝑎𝑥 during a single action potential has the consistent result. A4. Estimation magnitude of current through membrane
We assume that voltage-gated channels on axon membrane have almost returned to their resting state after axon firing, namely, 𝑔 𝑎𝑥𝑁𝑎,𝑑𝑦 ≈ 𝑔 𝑎𝑥𝑁𝑎,𝑟𝑒 , 𝑔 𝑎𝑥𝐾,𝑑𝑦 ≈ 𝑔 𝑎𝑥𝐾,𝑟𝑒 . We have ion channel conductance on the glial and axon membrane as 𝑔 𝑎𝑥𝑁𝑎,𝑟𝑒 , 𝑔 𝑎𝑥𝐾,𝑟𝑒 , 𝑔 𝑎𝑥𝐶𝑙 , 𝑔 𝑔𝑙𝐶𝑙 , 𝑔 𝑔𝑙𝑁𝑎 << 𝑔 𝑔𝑙𝐾 , 𝑚 𝑑𝑦 (𝑡 𝑚1 ) 𝑡 𝑚1 ms ms ms 𝑚 𝑑𝑦 (𝑡 𝑚2 ) 𝑡 𝑚2 ms ms ms
4 Similar to eq. (4.21), we claim 𝛿𝐸 𝑔𝑙,𝑎𝑥𝑁𝑎,𝐶𝑙 ≪ 𝛿𝐸 𝑔𝑙𝐾 , since O( 𝛿𝑐 𝑒𝑥𝑁𝑎 ) = −𝛿𝑐 𝑠𝑡𝑖 , 𝑂(𝛿𝑐 𝑒𝑥𝐾 ) = 𝛿𝑐 𝑠𝑡𝑖 , 𝑂(𝛿𝑐 𝑒𝑥𝐶𝑙 ) ≪ 𝛿𝑐 𝑠𝑡𝑖 and 𝑐 𝑒𝑥𝐾,𝑟𝑒 ≪ 𝑐 𝑒𝑥𝑁𝑎,𝐶𝑙,𝑟𝑒 . In addition, for the Na / K pump 𝛿𝐽 𝑝,𝑙𝑁𝑎 + 𝛿𝐽 𝑝,𝑙𝐾 = 𝛿𝐼 𝑙 𝑒 , 𝑙 = 𝑎𝑥, 𝑔𝑙. A Taylor expansion approximates the increase of the pump current due to the extracellular K + concentration variation . 𝛿𝐼 𝑙 ≈ 2 ( K K1 𝐼 𝑙𝑟𝑒,1 𝑐 𝑒𝑥𝐾,𝑟𝑒 (𝑐 𝑒𝑥𝐾,𝑟𝑒 +𝐾 𝐾1 ) + K K2 𝐼 𝑙𝑟𝑒,2 𝑐 𝑒𝑥𝐾,𝑟𝑒 (𝑐 𝑒𝑥𝐾,𝑟𝑒 +𝐾 𝐾2 ) ) 𝛿𝑐 𝑒𝑥𝐾 . (A - From the eq. (4.21) and eq. ( A-6 ), we have 𝛿𝐼 𝑙 << 𝑔 𝑔𝑙𝐾 𝛿𝐸 𝑔𝑙𝐾 , 𝑙 = 𝑎𝑥, 𝑔𝑙. In all, we claim the dominated current in the right hand side of eq. (4.23) is through glial membrane K + channel, such that ∑ 𝑖 𝑧 𝑖 𝑒ℳ 𝑔𝑙 (𝐽 𝑝,𝑔𝑙𝑖 + 𝐽 𝑐,𝑔𝑙𝑖 ) + ∑ 𝑖 𝑧 𝑖 𝑒ℳ 𝑎𝑥 (𝐽 𝑝,𝑎𝑥𝑖 + 𝐽 𝑐,𝑎𝑥𝑖 ) ≈ ℳ 𝑔𝑙 𝑔 𝑔𝑙𝐾 (𝛿𝑉 𝑔𝑙 − 𝛿𝐸 𝑔𝑙𝐾 ) . where we use the fact that at the resting state, the transmembrane currents in both axon membrane and glial membrane are negligible in compare to the source term 𝑔 𝑔𝑙𝐾 𝛿𝐸 𝑔𝑙𝐾 . A5. Magnitude comparison between 𝜹𝝓 𝒈𝒍 and 𝜹𝝓 𝒆𝒙 In this section, we claim that |𝛿𝜙 𝑔𝑙 | ≫ |𝛿𝜙 𝑒𝑥 | based on eq. (4.30) where we know 𝑂 ( 𝜂 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝜂 𝑒𝑥 𝜏 𝑒𝑥 𝜎 𝑒𝑥 ) = 10 −2 and 𝑂 ( D exdiff 𝑘 𝐵 𝑇𝛿𝑐 𝑠𝑡𝑖 𝜎 𝑒𝑥 𝑒 ) = 10 −6 . If the claim |𝛿𝜙 𝑔𝑙 | ≫ |𝛿𝜙 𝑒𝑥 | is not true, then from the expression of 𝛿𝜙 𝑒𝑥 in (4.30), the dominated term would be 𝑂 ( D exdiff 𝑘 𝐵 𝑇𝛿𝑐 𝑠𝑡𝑖 𝜎 𝑒𝑥 𝑒 ) = 10 −6 and 𝑂(𝛿𝜙 𝑔𝑙 ) ≤ 10 −5 . Then the right-hand side of (4.29) gives |𝑔 𝑔𝑙𝐾 𝑒 ( 𝛿𝑉 𝑔𝑙 − 𝛿𝐸 𝑔𝑙𝐾 )| ≈ |𝑔 𝑔𝑙𝐾 𝑒 𝛿𝐸 𝑔𝑙𝐾 | ~ 10 −8 mol/m s, where we use the estimation of 𝛿𝐸 𝑔𝑙𝐾 (= 𝑂(10 −3 )) in eq. (4.21) and 𝑂(𝛿𝑉 𝑔𝑙 ) = 𝑂(𝛿𝜙 𝑔𝑙 − 𝛿𝜙 𝑒𝑥 ) ≤ 10 −5 . While the left-hand side of eq. (4.29) gives | 𝑠𝑡𝑖 𝜂 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝑒 ℳ𝑔𝑙 𝑘 𝐵 𝑇 𝛿𝜙 𝑔𝑙 𝑟 ∗ | ~ 10 −11 mol/m s. Therefore, in eq. ( ), the order of right-hand side does match the order of the left-hand side. Therefore, we conclude that |𝛿𝜙 𝑔𝑙 | ≫ |𝛿𝜙 𝑒𝑥 |. A6. Estimation on the transport of 𝐍𝐚 + and 𝐊 + in extracellular space For the K + in the extracellular space in eq. (4.38), by using eq. (4.21) and eq. (4.33), we have the glial transmembrane effect of transport K + 𝜆 𝑔𝑙𝑚,𝐾 ≈ ℳ𝑔𝑙 𝑔 𝑔𝑙𝐾 ℎ 𝜖 𝑘 𝐵 𝑇𝑧 𝐾 (1+ℎ 𝜖 )𝑒 𝑐 𝑒𝑥𝐾,𝑟𝑒 . ( A - ) For K + flux through the extracellular pathway, we only consider the diffusion and electric drift effect in the radial direction because the length scale difference between longitudinal direction (𝑧 ∗ ) and radial direction (𝑟 ∗ ) . And the effect of extracellular osmosis variation produces water circulation in the optic nerve has been discussed it in Remark 4.3. The scale of the diffusion K + flux term in the radial direction in the extracellular space can be approximated as 𝑂 (−𝐷 𝑒𝑥𝐾 𝜏 𝑒𝑥 𝑑𝑐 𝑒𝑥𝐾 𝑑𝑟 ) ≈ 𝐷 𝑒𝑥𝐾 𝜏 𝑒𝑥 𝛿𝑐 𝑒𝑥𝐾 𝑟 ∗ . ( A - ) The scale of the electric drift of K + flux in the radial direction in the extracellular space is 𝑂 (− 𝐷 𝑒𝑥𝐾 𝜏 𝑒𝑥 𝑒𝑘 𝐵 𝑇 𝑐 𝑒𝑥𝐾 𝑑𝜙 𝑒𝑥 𝑑𝑟 ) ≈ 𝐷 𝑒𝑥𝐾 𝜏 𝑒𝑥 𝑒𝑘 𝐵 𝑇 𝑐 𝑒𝑥𝐾 𝛿𝜙 𝑒𝑥 𝑟 ∗ , ≈ − 𝜂 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝜂 𝑒𝑥 𝜏 𝑒𝑥 𝜎 𝑒𝑥 (1+ℎ 𝜖 ) 𝐷 𝑒𝑥𝐾 𝜏 𝑒𝑥 𝛿𝑐 𝑒𝑥𝐾 𝑟 ∗ ( A - ) where 𝛿𝜙 𝑒𝑥 used the estimation from eq. (4.34). Based on eq. (A-8) and eq. (A-9), in the extracellular space, we note that the electric drift for K + flux is in the opposite direction to the diffusion K + flux. The electric drift has a much smaller magnitude than the diffusion flux, because the ratio 𝑅 𝑒𝑥𝐾 between the scale of electric drift and scale of diffusion is 𝑅 𝑒𝑥𝐾 = 𝜂 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝜂 𝑒𝑥 𝜏 𝑒𝑥 𝜎 𝑒𝑥 (1+ℎ 𝜖 ) = 𝑜(1) . ( A - ) Therefore, the effect of K + flux in extracellular space in eq. (4.38) can be approximated as 𝜆 𝑒𝑥𝐾 = 𝑒𝑥𝐾 𝜏 𝑒𝑥 𝜂 𝑒𝑥 𝑟 ∗ 𝑟 𝑠𝑡𝑖 , ( A - ) 6 where we use the ratio of volume and effective surface, 𝜂 𝑒𝑥 𝑆 𝑡 𝑉 𝑆 = 𝑒𝑥 𝑟 𝑠𝑡𝑖 . In eq. ( ), we first look for the effect of Na + fluxes in the extracellular space. Similar to eq. (A-8) the (scale of) diffusion Na + flux term in the radial direction in the extracellular space is 𝑂 (−𝐷 𝑒𝑥𝑁𝑎 𝜏 𝑒𝑥 𝑑𝑐 𝑒𝑥𝑁𝑎 𝑑𝑟 ) ≈ 𝐷 𝑒𝑥𝑁𝑎 𝜏 𝑒𝑥 𝛿𝑐 𝑒𝑥𝑁𝑎 𝑟 ∗ . (A - The scale of the electric drift of Na + in the radial direction in the extracellular space is 𝑂 (− 𝐷 𝑒𝑥𝑁𝑎 𝜏 𝑒𝑥 𝑒𝑘 𝐵 𝑇 𝑐 𝑒𝑥𝑁𝑎 𝑑𝜙 𝑒𝑥 𝑑𝑟 ) ≈ 𝐷 𝑒𝑥𝑁𝑎 𝜏 𝑒𝑥 𝑒𝑘 𝐵 𝑇 𝑐 𝑒𝑥𝑁𝑎 𝛿𝜙 𝑒𝑥 𝑟 ∗ , ≈ − 𝜂 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝜂 𝑒𝑥 𝜏 𝑒𝑥 𝜎 𝑒𝑥 (1+ℎ 𝜖 ) 𝑐 𝑒𝑥𝑁𝑎 𝑐 𝑒𝑥𝐾 𝐷 𝑒𝑥𝑁𝑎 𝜏 𝑒𝑥 𝛿𝑐 𝑒𝑥𝐾 𝑟 ∗ . ( A - ) For Na + in the extracellular space, the electric drift flux for Na + is in the same direction as the diffusion K + flux since in the stimulated region 𝛿𝑐 𝑒𝑥𝑁𝑎 = −𝑂(𝛿𝑐 𝑒𝑥𝐾 ) . ( A - ) The Na + diffusion flux in the extracellular is about as big as the electric drift component, namely from eq. (A-12) and eq. (A-13) 𝑅 𝑒𝑥𝑁𝑎 = 𝜂 𝑔𝑙 𝜏 𝑔𝑙 𝜎 𝑔𝑙 𝜂 𝑒𝑥 𝜏 𝑒𝑥 𝜎 𝑒𝑥 (1+ℎ 𝜖 ) 𝑐 𝑒𝑥𝑁𝑎 𝑐 𝑒𝑥𝐾 = 𝑂(1) . (A-15) The Na + flux through glial transmembrane is much smaller than K + flux such that 𝜆 𝑔𝑙𝑚,𝑁𝑎 = 𝑜(𝜆 𝑔𝑙𝑚,𝐾 ) . (A-16) since the conductance on the glial membrane 𝑔 𝑔𝑙𝑁𝑎 ≪ 𝑔 𝑔𝑙𝐾 . For Na + , the effect of Na + flux through glial transmembrane can be neglected in eq. ( ), since eq. (A-12) and eq. (A-8) are in the same magnitude and eq. (A-16) show the for Na + flux is much smaller than K + on the glial membrane. In the end of this section, we consider the solution for the coupled dynamic system for eq. (4.38) and eq. (4.40), 𝑑𝑑𝑡 ( 𝛿𝑐 𝑒𝑥𝐾 𝛿𝑐 𝑒𝑥𝑁𝑎 ) = 𝐴 ( 𝛿𝑐 𝑒𝑥𝐾 𝛿𝑐 𝑒𝑥𝑁𝑎 ) , (A-17) where 𝐴 = [𝐴 𝐴 ] = [−(𝜆 𝑔𝑙𝑚,𝐾 + 𝜆 𝑒𝑥𝐾 )/𝜂 𝑒𝑥𝑟𝑒 𝑒𝑥𝑁𝑎,2 /𝜂 𝑒𝑥𝑟𝑒 −𝜆 𝑒𝑥𝑁𝑎,1 /𝜂 𝑒𝑥𝑟𝑒 ] . In the system (A-17), we assume that 𝜂 𝑒𝑥 keep at its resting state as 𝜂 𝑒𝑥𝑟𝑒 . The initial value of the system (A-17) is 7 ( 𝛿𝑐 𝑒𝑥𝐾,0 𝛿𝑐 𝑒𝑥𝑁𝑎,0 ) = ( 𝛿𝑐 𝑠𝑡𝑖 −𝛿𝑐 𝑠𝑡𝑖 ). When 𝑡 ∈ [0, 𝑇] , we have the following solution for the variations of potassium and solidum in the extracellular space as { 𝛿𝑐 𝑒𝑥𝐾 (𝑡) = 𝛿𝑐 𝑠𝑡𝑖 exp(𝐴 𝑡), 𝛿𝑐 𝑒𝑥𝑁𝑎 (𝑡) = 𝐴 𝛿𝑐 𝑠𝑡𝑖 𝐴 −𝐴 (exp(𝐴 𝑡) − exp(𝐴 𝑡)) − 𝛿𝑐 𝑠𝑡𝑖 exp(𝐴 𝑡 ) . (A-18) Based on the solution in (A-18), we provide solution when there is a frequency of 𝑛 stimuli in the time interval [0, 𝑇 𝑠𝑡𝑖 = 𝑛𝑇] . We view axon firing in the stimulated region, as a source term for 𝑐 𝑒𝑥𝐾 and a sink term for 𝑐 𝑒𝑥𝑁𝑎 . We have 𝑐 𝑒𝑥𝐾 and 𝑐 𝑒𝑥𝑁𝑎 change at each time point 𝑡 =𝑇, 2𝑇, … (𝑛 − 1)𝑇 as 𝛿𝑐 𝑒𝑥𝐾 (𝑖𝑇) = 𝛿𝑐 𝑒𝑥𝐾 (𝑖𝑇) + 𝛿𝑐 𝑠𝑡𝑖 , 𝛿𝑐 𝑒𝑥𝑁𝑎 (𝑖𝑇) = 𝛿𝑐 𝑒𝑥𝑁𝑎 (𝑖𝑇) − 𝛿𝑐 𝑠𝑡𝑖 , for 𝑖 = 1, ⋯ 𝑛 − 1. By using eq. (A-18), we have at time 𝑡 = 𝑛𝑇 𝛿𝑐 𝑒𝑥𝐾 (𝑛𝑇) = 𝛿𝑐 𝑠𝑡𝑖 ∑ exp(𝑖𝐴 𝑇) 𝑛𝑖=1 = 𝛿𝑐 𝑠𝑡𝑖 exp(𝐴 𝑇)−exp((𝑛+1)𝐴 𝑇)1−exp(𝐴 𝑇) , and 𝛿𝑐 𝑒𝑥𝑁𝑎 (𝑛𝑇) = ∑ 𝐴 𝛿𝑐 𝑒𝑥𝐾 ((𝑖−1)𝑇)𝐴 −𝐴 (exp(𝐴 𝑇) − exp(𝐴 𝑇)) 𝑛𝑖=1 exp( (𝑛 − 𝑖)𝐴 𝑇 ) − 𝛿𝑐 𝑠𝑡𝑖 ∑ exp(𝑖𝐴 𝑇) 𝑛𝑖=1 . where 𝛿𝑐 𝑒𝑥𝐾 (𝑗𝑇) = 𝛿𝑐 𝑠𝑡𝑖 1−exp((𝑗+1)𝐴 𝑇)1−exp(𝐴 𝑇) , 𝑗 = 0,1, … 𝑛 − 1., 𝑗 = 0,1, … 𝑛 − 1.