Fig. 1 shows the structure we are modeling. This is a subset of all the structures of interest in the optic nerve, but it is a good place to start. Computational domain $\mathrm{\Omega }$consists of the subarachnoid space (SAS) region ${\mathrm{\Omega }}_{SAS}$ and optic nerve region ${\mathrm{\Omega }}_{Op}$, i.e.,

$$\mathrm{\Omega }={\mathrm{\Omega }}_{OP}\cup {\mathrm{\Omega }}_{SAS},{\mathrm{\Omega }}_{OP}\cap {\mathrm{\Omega }}_{SAS}={\mathrm{\Gamma }}_7,$$

where the domain ${\mathrm{\Omega }}_{SAS}$ is filled with CSF, enclosed by dura mater ${\mathrm{\Gamma }}_7$ and pia mater ${\mathrm{\Gamma }}_4$.

For the optic nerve region, based on the model first proposed in Ref. [22], we introduce a six-domain model for microcirculation of the optic nerve, which is summarized in Table 1. We must remind the reader that the optic nerve bundler contains many nerve fibers, at least one blood vessel, and a glial syncytium as shown in Fig. 1. In addition to the axon (${\mathrm{\Omega }}_{ax}$), glial (${\mathrm{\Omega }}_{gl}$), and extracellular (${\mathrm{\Omega }}_{ex}$) compartments, the model incorporates three perivascular spaces: perivascular space A (surrounding the artery, ${\mathrm{\Omega }}_{pa}$), perivascular space V (surrounding the vein, ${\mathrm{\Omega }}_{pv}$), and perivascular space C (surrounding the capillary, ${\mathrm{\Omega }}_{pc}$), see Fig. 2.

$${\mathrm{\Omega }}_{OP}={\mathrm{\Omega }}_{ax}\cup {\mathrm{\Omega }}_{gl}\cup {\mathrm{\Omega }}_{ex}\cup {\mathrm{\Omega }}_{pa}\cup {\mathrm{\Omega }}_{pv}\cup {\mathrm{\Omega }}_{pc}.$$

We apply compartmental mass conservation laws to describe the dynamics of ions and water within and across cellular and extracellular domains. These laws ensure that the model respects fundamental physical principles—such as the conservation of mass and charge—and establish a robust foundation for physiologically realistic and numerically stable simulations. Specifically, the model is derived from ion and fluid conservation laws governing transmembrane transport between intracellular compartments and the ECS, following the framework introduced in [24]. These conservation principles are applied in each of the six domains: ${\mathrm{\Omega }}_l$, $l=ax$, $gl$, $ex$, $pa$, $pv$, $pc$.

$$\frac{\partial }{\partial t}({\eta }_l{f}_l)+\nabla \cdot ({\eta }_l{𝑱}_l)+S=0$$

where ${\eta }_l$ is the volume fraction of l compartment, $f_l$ is the concentration of a given substance, ${𝑱}_l$ is the flux inside the compartment, and S is the source term induced by trans-membrane communications due to the active pumps and passive leak channels on the membranes.

For the boundaries, ${\mathrm{\Gamma }}_1$ is the central retinal blood wall of the optic nerve; ${\mathrm{\Gamma }}_2$ and ${\mathrm{\Gamma }}_3$ are the far end (away from the eyeball) of the optic nerve which is connected to optic canal region [25]. ${\mathrm{\Gamma }}_5$ is used to model the dura mater connected to the sclera (the white matter of the eye) and assumed to be non-permeable [12]. ${\mathrm{\Gamma }}_6$ is used to denote the lamina cribrosa where the optic nerve head exits the eye posteriorly through pores of the lamina cribrosa [26]. We adopt membrane-type boundary conditions to model transmembrane fluxes and compartmental interactions. This choice is motivated by the anatomical context: the intraorbital segment of the optic nerve is embedded in and interacts with surrounding tissues rather than existing in isolation. The membrane boundary conditions allow us to incorporate physiological coupling with adjacent structures, and any region-specific properties (e.g., diffusivity, permeability) can be tuned via coefficients without altering the core model. Additional technical details regarding the boundary condition implementation are provided in the Supplementary material.

In this model, the electric potential ${\varphi }_l$​ within each compartment ${\mathrm{\Omega }}_l$​ is not externally imposed but emerges from the redistribution of mobile ions during neuronal activation. When potassium ions efflux from stimulated axons into the ECS, they cause a local accumulation of positive charge. To restore ionic balance, these ions are subsequently taken up by neighboring glial cells, leading to regional shifts in intracellular electric potential. This spatial potential difference generates an electric field that further drives ionic electrodiffusion within and across compartments.

Rather than solving the Poisson equation for ${\varphi }_l$, we adopt a macroscopic electroneutrality assumption, which is valid at tissue scale due to the extremely short Debye length in biological media. Under this assumption, charge separation is negligible compared to the domain size, and the electric potential is determined implicitly by enforcing charge neutrality in each compartment, combined with ionic flux dynamics. This approach is widely used in tissue-level models of electrodiffusion and has the advantage of avoiding the numerical stiffness associated with solving elliptic equations for the electric field. Mathematically, the formulation of electric potential under electroneutrality constraints is provided in Supplementary material, where ${\varphi }_l$ is recovered from ionic concentrations via an algebraic closure condition derived from the ion flux equations. To ensure consistency, we apply homogeneous Neumann boundary conditions for the electric potential on the domain boundaries implying no external electric field enters or leaves the system. This setup isolates the modeled domain electrically and ensures that all electric fields arise from internal biophysical processes, such as ionic imbalance and transmembrane fluxes, consistent with volume conductor theory in neuroscience.

Based on the structure of the optic nerve, we have the following global assumptions for the model (Table 1).

$\bullet$ Axial symmetry: For simplicity, axial symmetry is assumed. To clarify, the model is fundamentally three-dimensional. When axial symmetry is assumed—as in our current implementation—it simplifies the numerical computation without altering the underlying physics. Importantly, the model formulation itself does not rely on symmetry assumptions and can be extended to fully three-dimensional, non-axisymmetric geometries to accommodate spatial heterogeneity or more anatomically realistic structures.

$\bullet$ Ions type: To focus on the primary mechanisms of ionic homeostasis and fluid–ion interaction, the model includes only three major ion species: sodium (Na⁺), potassium (K⁺), and chloride (Cl^–). These ions are the principal contributors to transmembrane potential, osmotic gradients, and electrochemical fluxes under typical physiological conditions. We acknowledge that other ions and metabolites—such as calcium (Ca²⁺), protons (H⁺), lactate, and glutamate—also play important roles in central nervous system function, particularly in regulating excitability, pH balance, and neuro-glial signaling. However, these species are not included in the current model in order to maintain computational tractability and focus on the dominant components involved in potassium buffering and clearance.

$\bullet$ Charge neutrality: In each domain, we assume electroneutrality [27],

$${\eta }_{gl}\sum _i{z}^i{C}_{gl}^i+z^{gl}{\eta }_{gl}^{re}{A}_{gl}=0$$

$${\eta }_{ax}\sum _i{z}^i{C}_{ax}^i+z^{ax}{\eta }_{ax}^{re}{A}_{ax}=0$$

$$\sum _i{z}^i{C}_l^i=0,l=ex,pa,pv,pc,csf,$$

where $A_l>0$ with $l=ax$, $gl$ is the density of proteins in axons or glial cells. Most CNS proteins—such as albumin—have isoelectric points (pI) well below the physiological pH of CSF, which is approximately 7.3–7.4. At this pH, these proteins exist predominantly in their deprotonated form and therefore carry a net negative charge. This property is a general consequence of the amino acid composition of these proteins, where acidic residues (e.g., glutamate and aspartate) contribute to the lower pI values. In our model, we incorporate this biochemical reality by assuming that these negatively charged proteins are permanently distributed within axonal and glial compartments, with an effective valence $z^l=-1,l=ax,gl$. The ${\eta }_{ax}$ and ${\eta }_{gl}$ are the volume fraction of axon and glial compartments in the optic nerve and ${\eta }_{gl}^{re}$ and ${\eta }_{ax}^{re}$ are the resting state volume fractions. Under physiological conditions, the extracellular fluid—such as CSF and interstitial fluid—maintains macroscopic electroneutrality, meaning the total concentrations of positive and negative charges are balanced. While transient microscopic charge separations occur (e.g., across the membrane during action potentials), the tissue-scale charge density remains effectively zero, minimizing the risk of unbalanced electrostatic fields. This assumption not only reflects physiological reality but also circumvents the numerical complexity of solving the full Poisson equation, thus enhancing computational efficiency without sacrificing model fidelity.

$\bullet$ Anisotropy of axon compartment and isotropy of other compartments: The axons are separated, more or less parallel cylindrical cells that form separate compartments. They are electrically isolated and molecules cannot diffuse directly from one to another. As a result, the intra-compartment flux of ions and fluid are only along the axis direction. In contrast, other compartments are fully connected. For example, the glial cells are connected by connexins and form a syncytium, allowing the intra-compartment fluxes flow in both axial and radial directions.

$\bullet$ Communications between compartments: The communications among different compartments are illustrated in Fig. 3. Especially, there is no direct interaction between glial (or perivascular space) and axon compartments. Interactions occur only through changes in concentration, electrical potential, and flows in ECS [28].

The details of the mathematical model can be found in the SI.

Under a pressure gradient of 0.0083 mmHg/mm, CSF flows into the subarachnoid space (SAS) from the intracranial region, with an average velocity of approximately 250 µm/s in the z-direction, consistent with previous findings [31]. CSF then leaks into perivascular space A. The flow direction in perivascular space A (or V) matches the blood flow in the central retinal artery (or vein). The fluid flow in perivascular space C adjusts dynamically through exchanges with other compartments, primarily influenced by perivascular space V due to a larger pressure gradient at boundary $\mathrm{\Gamma }$₁, with fluid exiting the optic nerve at boundary $\mathrm{\Gamma }$₂ under membrane boundary conditions.

In perivascular space A, the spatially averaged fluid velocity under a 0.007 mmHg/mm pressure gradient in the z-direction is 5 µm/s from the intracanalicular space to the intraorbital region. In perivascular space V, the spatially averaged fluid velocity in the z-direction under a –0.012 mmHg/mm pressure gradient is 4 µm/s from the intraorbital region to the intracanalicular space. These results align with the observations in [17, 32].

Simultaneously, fluid flows from pvsA into the ECS and the glial compartment through gaps and Aquaporin channels located near or on the glial cell feet [14], and then into perivascular space V. Fluid also moves from the ECS into the glial compartment, in agreement with the glymphatic system’s role [5, 13, 33].

This section examines the glymphatic system’s role (including glial cells and perivascular spaces) in metabolic waste clearance during stimulation. As depicted in Fig. 4, the stimulus is applied to the axon membrane within the region at a specified location $z=z_0$, with $\theta \in [0,\mathrm{\hspace{0.25em}2}\pi ]$. We differentiate between the stimulated and non-stimulated regions in the optic nerve ${\mathrm{\Omega }}_{Op}$, where electrical signals propagate in the z-direction within the axon compartment. The stimulus frequency is 50 Hz (T = 0.02 s) with a duration of 0.2 s. Each stimulus has a current strength of I_s⁢t⁢i = 3 $\times$ 10^–⁢3 A/m² and lasts for 3 ms.

In this section, we investigate the role of the glial compartment in extracellular potassium buffering by varying the transmembrane ionic conductance ($g_{gl}^i$), hydraulic permeability, and modifying connexin connectivity, which reduces diffusion coefficients and permeability within the glial compartment.