Neural activity induces strongly coupled electro-chemo-mechanical interactions and fluid flow in astrocyte networks and extracellular space—A computational study

The complex interplay between chemical, electrical, and mechanical factors is fundamental to the function and homeostasis of the brain, but the effect of electrochemical gradients on brain interstitial fluid flow, solute transport, and clearance remains poorly quantified. Here, via in-silico experiments based on biophysical modeling, we estimate water movement across astrocyte cell membranes, within astrocyte networks, and within the extracellular space (ECS) induced by neuronal activity, and quantify the relative role of different forces (osmotic, hydrostatic, and electrical) on transport and fluid flow under such conditions. We find that neuronal activity alone may induce intracellular fluid velocities in astrocyte networks of up to 14μm/min, and fluid velocities in the ECS of similar magnitude. These velocities are dominated by an osmotic contribution in the intracellular compartment; without it, the estimated fluid velocities drop by a factor of ×34–45. Further, the compartmental fluid flow has a pronounced effect on transport: advection accelerates ionic transport within astrocytic networks by a factor of ×1–5 compared to diffusion alone.

both ionic electrodiffusion and fluid movement in an arbitrary number of compartments is presented by Mori [38]. In this framework, hydrostatic and osmotic pressure gradients are assumed to drive fluid flow in both intra-and extracellular spaces. Zhu et al. [39] later extend this framework by including electro-osmosis as a driving force for interstitial fluid movement. They apply their framework to study the role of fluid flow on ionic transport in the optical nerve. However, little attention has been paid to quantifying the contributions of different driving forces for interstitial fluid flow during neuronal activity in the cortex.
Here, our target is two-fold: we aim to estimate the water movement induced by neuronal activity across astrocyte cell membranes, within astrocyte networks, and within the extracellular space, and to determine the relative role of different forces (osmotic, hydrostatic, and electrical) on ionic transport and fluid flow under such conditions. To estimate this electrochemo-mechanical response, we introduce a high-fidelity computational model describing the spatial and temporal dynamics at the micro/milliscale of volume fractions, electrical potentials, ion concentrations, and hydrostatic pressures in an intracellular space (ICS) representing different astrocyte configurations and the extracellular space (ECS) (Fig 1). The model is embedded in the electrodiffusive Kirchhoff-Nernst-Planck framework and builds on previous work [40] incorporating ionic electrodiffusion [34,38], fluid dynamics [39], and astrocyte modeling [34].
Our findings show that neuronal activity in the form of extracellular ionic input fluxes induces complex and strongly-coupled chemical-electrical-mechanical interactions in the astrocytic ICS and ECS. The response is characterized by membrane (electric) potential depolarization on the order of tens of millivolts and ECS water potentials ranging from a few to a hundred kilo-pascals, spatial differences in osmolarity on the order of several tens of millimolars, and fluid velocities ranging from a fraction-to tens of micrometers per minute. The fluid dynamics are crucially coupled to the spatial organization of the intracellular network. We observe intracellular fluid velocities in astrocyte networks of up to 14 μm/min, and fluid velocities in the ECS of similar magnitude. These velocities are dominated by an osmotic contribution in the intracellular compartment; without it, the estimated fluid velocities drop by a factor of ×34-45. Furthermore, the compartmental fluid flow has a pronounced effect on transport: advection accelerates ionic transport within astrocytic networks by a factor of ×1-5 compared to diffusion alone.

Results
In order to quantify the relative role of osmotic, hydrostatic, and electrical forces on transport and flow in cortical tissue, we ask the following questions. How do astrocyte and extracellular ion concentrations, electric potentials, pressures, and interstitial fluid velocities respond to changes in extracellular ion concentrations mirroring neural activity on the time scale of seconds? Moreover, to what extent do the mechanical responses (cellular swelling, fluid flow) contribute to alleviating ionic and mechanical ECS distress? To address these questions, we introduce a set of biophysical models for these quantities of interest, governed by the balance of mass, momentum, and charge, in combination with astrocyte membrane mechanisms in a representative volume.

A model for electrodiffusive, osmotic, and hydrostatic interplay in astrocyte networks
Ion-and fluid movement in an astrocyte network (ICS) and the extracellular space (ECS) is modeled via coupled partial differential equations (PDEs) in a homogenized model domain (Fig 1, Methods). Specifically, we consider a 1D domain of length 300 μm representing brain tissue between two blood vessels, e.g., an arteriole and a venule. The model predicts the evolution in time and distribution in space of the volume fraction α r , the ion concentrations [Na + ] r , [K + ] r , and [Cl − ] r , the electrical potential ϕ r , and the hydrostatic pressure p r in both ICS (r = i) and ECS (r = e). Ionic transport is driven by diffusion, electric drift, and advection. To model fluid movement in each compartment, we consider three different model scenarios:

M1
The intra-and extracellular fluid flow is driven by hydrostatic pressure gradients. The astrocytic compartment can be interpreted as either a single closed cell or as a syncytium of cells without intercellular osmotic flow.

M2
The intracellular fluid flow is driven by osmotic and hydrostatic pressure gradients, and the extracellular fluid flow is driven by the same mechanism as in M1. The astrocytic compartment can be interpreted as a syncytium of cells where osmosis acts as a driving force for fluid flow.

M3
The intracellular fluid flow model is the same as in M2, and the extracellular fluid flow is driven by electro-osmosis in addition to hydrostatic pressure gradients. To include electroosmosis as a driving force for ECS fluid flow is motivated by the narrowness of the ECS [39,41].
For comparison, we also consider a zero-flow scenario [34]: M0 The compartmental fluid velocities and the transmembrane fluid flow (and thus cellular swelling) are assumed to be zero.
Transmembrane fluid flow in model scenarios M1-M3 is driven by hydrostatic and osmotic pressure differences. At the boundaries, we assume that no fluid and no ionic fluxes enter or leave the system. To account for transmembrane ionic movement, we include an Illustration of brain tissue between two blood vessels with astrocytes (purple), neurons (grey), and ECS with neural activity in the center (A). The tissue is represented as a 1D domain of length 300 μm including ICS (astrocytes) and the ECS (B). Within each compartment, the model describes the dynamics of the volume fraction (α), the Na + , K + , and Cl − concentrations ([Na + ], [K + ], [Cl − ]), the electrical potential (ϕ), and the hydrostatic pressure (p). Neuronal activity is implicitly represented by K + and Na + input currents (j K input and j Na input ) in the input zone (of length 30 μm) and decay currents (j K decay and j Na decay ) across the whole domain. Transmembrane currents include an inward rectifying K + current (j Kir ), Na + and Cl − leak currents (j K leak and j Cl leak ), and a Na + /K + pump current (j pump ). Intra-and extracellular currents (j k i and j k e ) are driven by electrodiffusion and advection. Fluid can travel across the membrane (w m ) and compartmentally in the intra-and extracellular space (u i and u e ).
https://doi.org/10.1371/journal.pcbi.1010996.g001 inward rectifying K + channel, passive Na + and Cl − channels, and a Na + /K + pump. To ensure electroneutrality of the system, we include a set of immobile anions a r . The immobile anions contribute to the osmotic pressures.
We mimic a scenario of high local neuronal activity by injecting a constant K + current into the ECS and simultaneously removing Na + ions in a stimulus zone in the middle of the computational domain. We set the strength of the input currents such that the extracellular K + concentration in the input zone reaches a maximum value of approximately 10 mM during the simulations. At this concentration level, we expect the K + buffering process to play a critical role; still, the concentration is below the level we observe in pathological conditions such as spreading depression (see Halnes et al. [34] and references therein). To maintain electroneutrality of the system, the K + and Na + input currents are of the same magnitude. Neuronal pumps and cotransporters are accounted for by removing K + ions from the ECS at a given decay rate and adding the same amount of Na + ions. The decay is proportional to the extracellular K + concentration and defined across the whole domain. We also study the effect of timevarying input by injecting pulsatile K + currents into the ECS. Note that the stimulus does not induce any osmotic pressure changes, as it does not affect the total osmotic concentration of the ECS.

Neuronal activity induces complex chemical-electrical-mechanical interplay
In order to understand and quantify the baseline electrical and mechanical response to chemical alterations, we first consider the model scenario where the compartmental fluid flow is only driven by hydrostatic pressure gradients (M1). Turning on the input currents at t = 10 s leads to changes in the ion concentrations, cellular swelling, depolarization of the glial membrane, and an increase in the transmembrane hydrostatic pressure difference. After about 40 seconds, the system reaches a new steady state, before all fields return to baseline levels after input offset at t = 210 s (Fig 2).
The input currents (Fig 2A) and the subsequent astrocytic activity lead to an increase of 6.68 mM in [K + ] e , a decrease of 18.7 mM in [Na + ] e , and a decrease of 16.9 mM in [Cl − ] e , measured at the center of the input zone ( Fig 2C). The increase in [K + ] e activates the K + -and Na +decay currents (Fig 2B), which eventually lead the system back to baseline. Intracellularly, we observe an initial peak in Δ[K + ] i of 3.19 mM before it settles on 0.428 mM. [Na + ] i and [Cl − ] i decrease by 6.44 mM and increase by 6.55 mM, respectively ( Fig 2D). In response to the ionic shifts, the astrocytic compartment swells: the ICS volume increases by 13% ( Fig 2F) and the ECS shrinks correspondingly by 26% ( Fig 2E). As the initial size of the ECS is half that of the ICS, a change in ECS volume twice that of the ICS volume is as expected. Further, these volume changes affect the hydrostatic pressures: the transmembrane hydrostatic pressure difference increases by 118 Pa (Fig 2G). Finally, the glial membrane potential depolarizes from −86 mV to −61 mV (Fig 2H).

Transmembrane dynamics induce hydrostatic pressure gradients and compartmental fluid flow
Osmotically driven transport of fluid through AQP4, and possibly other membrane mechanisms, play an important role in cellular swelling and volume control of the ECS [42][43][44][45]. Whether cellular swelling induces hydrostatic pressure gradients driving compartmental fluid flow in the ICS and ECS is, however, far from settled [34]. We therefore also assess to what extent osmotic pressures induce hydrostatic pressures and fluid flow, still in the model scenario with hydrostatic-pressure-driven compartmental fluid flow (M1).
The concentration shifts following astrocytic activity result in altered ICS and ECS osmolarities. Notably, the change in osmolarities peak in the input zone, where the intra-and extracellular osmolarities decrease by maximum 20.5 mM and 20.7 mM, respectively ( Fig 3A). Consequently, the osmotic pressure across the astrocytic membrane decreases by a maximum of 713 Pa (Fig 3B). The osmotic pressure drives fluid across the astrocytic membrane, with a maximum velocity of 0.003 μm/min (Fig 3C), resulting in cellular swelling. Following swelling, the ICS hydrostatic pressure increases by at most 21.2 Pa in the input zone, while the ECS hydrostatic pressure drops by at most 97.1 Pa (Fig 3D). Note that the changes in compartmental hydrostatic pressures lead to a change in the transmembrane hydrostatic pressure gradient (Fig 3B), which affects the transmembrane fluid flow. The ICS and ECS hydrostatic pressure gradients drive compartmental fluid flow forming two circulation zones (Fig 3E and 3F). The intra-and extracellular superficial fluid velocities, α r u r , peak at 0.31 μm/min ( Fig 3E). Note that the two superficial fluid velocities are opposite in direction but have the same magnitude-this is a direct consequence of the incompressibility condition and no-flux boundary conditions in one dimension. During steady-state and cellular swelling in the input zone, fluid flows across the membrane into the ICS (Fig 3C). In the ICS, the fluid consequently flows away from the input zone, with positive flow to the right of the swelling and negative flow to the left of the swelling (Fig 3E). The ECS fluid flow is in the opposite direction, towards the input zone ( Fig 3E).

Extracellular and intracellular fluid flow alleviate osmotic pressure build-up
In a previous study [46], we predicted the strength of osmotic pressure build-up across an astrocyte membrane using a classical electrodiffusive model not accounting for fluid flow. However, to what extent will swelling and compartmental fluid flow affect osmotic pressure build-up across the membrane? Here, we investigate this question by comparing predictions of the zero-flow model (M0), the flow model without intercellular osmotic flow (M1), and the flow model with osmotic intercellular flow (M2).
The ICS and ECS osmolarities are altered by astrocytic activity in all model scenarios, notably peaking in the input zone (Fig 4A and 4B). In the ECS, the osmolarity decrease by 36.5 mM, 20.7 mM, and 6.74 mM for M0, M1, and M2, respectively ( Fig 4A). Interestingly, the ICS osmolarity increases by 10.5 mM and 1.09 mM for M0 and M2, respectively, whereas it decreases for model scenario M1 by 20.45 mM (Fig 4B). The decrease in ICS osmolarity in model scenario M1 results from cellular swelling: the osmolarity is defined as the amount of ions per unit volume. An increase in cell volume may thus cause the ion concentration to drop even if the number of ions increases. Furthermore, we can convert the intra-and extracellular osmolarities to intra-and extracellular solute potentials, P i and P e , respectively (see Section 4.2 for further details). Taking the difference in solute potential across the membrane gives us the osmotic pressure, which differs substantially between the models: M0, M1, and M2 predict a maximum drop in osmotic pressure of respectively 121 kPa, 0.713 kPa, and 20.2 kPa (Fig 4C). Allowing for cellular swelling and compartmental fluid flow thus reduces the osmotic pressure across the membrane by 99.4% (M1) and 83.3% (M2). These findings suggest that model scenario M0, or generally any model for electrodiffusion not taking into account fluid dynamics, highly overestimates the osmotic pressure building up across the membrane.

Volume dynamics is essential for ECS homeostasis
Water uptake in astrocytes via e.g. AQP4 has been hypothesized to contribute to stabilizing ECS ion concentrations, and thus prevent severe neuronal swelling. Swelling is driven by the difference in intra-and extracellular water potential, which is given by the solute potential P r plus the hydrostatic pressure p r (see Methods for details). An increase in extracellular water potential will result in neuronal swelling as water flows along the potential gradient.
For all model scenarios (M0, M1, M2), the ECS solute potential increases ( Fig 4D). The increase is most severe for the zero-flow model (M0), which predicts a maximum increase of 94.1 kPa. When taking swelling and compartmental hydrostatic-pressure-driven fluid flow into account (M1), the ECS solute potential increases by maximum 53.5 kPa, whereas adding ICS osmotic forces (M2) leads to an increase of maximum 17.4 kPa. In M0, the change in p e is zero by definition, whereas M1 predicts a negligible maximal drop (0.0971 kPa, Fig 4E). Conversely, model scenario M2 predicts a pronounced hydrostatic pressure drop (5.67 kPa, Fig  4E). Consequently, the maximal change in the ECS water potential is substantially smaller in M2 (11.7 kPa, Fig 4F) than in M0 and M1 (respectively 94.1 kPa and 53.4 kPa, Fig 4F). The maximal change in the ECS solute potential is less severe for M2 than for M0 and M1 ( Fig 4D). Additionally, while the ECS solute potentials increase (Fig 4D), the ECS hydrostatic pressures decrease (Fig 4E), and thus drive the water potential in opposite directions. Consequently, the ECS hydrostatic pressure change in M2 reduces the contribution from the less severe change in ECS solute potential, resulting in a lower ECS water potential. Thus, cellular swelling and osmotic transport within the astrocytic network (M2) contribute to prevent water potential build-up in the ECS.

Astrocyte osmotics strengthens compartmental fluid flow
The presence of driving forces for fluid flow-both at the brain microscale and organ-levelremain an open question. To assess the potential contributions from osmosis in the astrocytic network (M2) and electro-osmosis in the ECS (M3), we here compare the compartmental fluid velocities, cellular swelling, and hydrostatic pressures predicted by model scenarios M1, M2, and M3.
The maximum superficial fluid velocity predicted by M2 is 14 μm/min-about 45 times larger than for M1 ( Fig 5A and 5B and Table 1). For M3, the maximum superficial fluid velocity is 13 μm/min-slightly smaller than for M2 ( Fig 5C and 5D and Table 1). For both M2 and M3, the fluid velocities are dominated by osmosis in the ICS (Fig 5A and 5C) and hydrostatic forces in the ECS (Fig 5B and 5D). Interestingly, the intracellular hydrostatic pressure gradient drives fluid towards the input zone in M2 and M3 ( Fig 5A and 5C), as opposed to the intracellular hydrostatic pressure in M1 which drives fluid away from the input zone ( Fig 3C). The difference in ICS flow direction predicted by M1, M2, and M3 arises from the coupling of the hydrostatic-, osmotic-, and electro-osmotic forces in the mathematical model: the osmoticand electro-osmotic forces are given by the ion concentration-and electrical potential gradients, respectively, whereas the hydrostatic pressure result from the incompressibility of the interstitial fluid (cf. Eq (7)). Less cellular swelling is predicted by M2 and M3 than by M1: the astrocyte swells by respectively 12.9%, 3.74%, and 4.55% in M1, M2, and M3 (Table 1). Finally, we observe notable differences in the intra-and extracellular hydrostatic pressures: the maximum ICS hydrostatic pressure is 1.02 kPa, −4.64 kPa, and −10.3 kPa in M1, M2, and M3, respectively ( Table 1). The maximum ECS hydrostatic pressure is -0.0971 kPa, -5.67 kPa, and -11.3 kPa in M1, M2, and M3, respectively ( Table 1). The transmembrane hydrostatic pressures are similar in M2 and M3 even if the intra-and extracellular pressures are different: the maximum transmembrane hydrostatic pressure is respectively 1.03 kPa and 1.04 kPa in M2 and M3 (Table 1). The maximum transmembrane hydrostatic pressure in M1 is 1.12 kPa (Table 1).

Astrocyte osmotics accelerates ionic transport and alters role of advection
In model scenarios M1-M3, the compartmental fluid velocities will contribute to ionic transport via advection. To assess the role of advection in compartmental ionic transport, we compare model scenarios M1 and M3. Specifically, we decompose the intra-and extracellular ionic fluxes and calculate the advection/diffusion fraction (F diff ) and the advection/drift fraction (F drift ) for each of the ionic species (see Methods for further details).
We observe that F diff and F drift range from 0.002 to 0.062 in model scenario M1, indicating that advection plays a negligible role in ionic transport (Fig 6A-6F). In model scenario M3, however, we observe a larger variability in the advection/diffusion-and advection/drift rates: F diff ranges from 0.058 to 4.519, and F drift ranges from 0.325 to 0.949 (Fig 6G-6L). For K + transport in the M3 model, the advective flux is on the same order of magnitude as the intraand extracellular electric drift and about 4.5 times stronger than the intracellular diffusion (F diff = 4.519, Fig 6G and 6J). Advection plays the most important role intracellularly and accelerates the K + transport; For M1, the intracellular K + flux has a maximum value of 58 μmol/(m 2 s) (Fig 6D), whereas for M3, the maximum intracellular K + flux is 70 μmol/(m 2 s) ( Fig 6J). For Na + and Cl − transport in the M3 model, advection is on the same order of magnitude as diffusion and electric drift (Fig 6H-6I and 6K-6L). The advection even dominates intracellular diffusion of Na + (F diff = 1.286, Fig 6K) and extracellular diffusion of Cl − (F diff = 4.079, Fig 6I). Overall, advection accelerates the transport of total charge in the system: For M1, the charge flux (defined here as z K j K + z Na j Na + z Cl j Cl ) is maximum 59 μmol/(m 2 s), whereas for M3, the charge flux is maximum 71 μmol/(m 2 s).

Fluid flow patterns induced by low-frequency activity waves
Till now, we have studied the quasi-static response resulting from a step-wise change in (constant) input fluxes. To study how low-frequency waves of ionic changes affect the osmotic, hydrostatic, and electrical response of the system, we compare three different stimulus protocols via model scenario M3: (I) (slow) input fluxes varying at a low (1Hz) frequency, (II) (ultraslow) input fluxes varying at an ultralow (0.05 Hz) frequency, and (III) constant input fluxes, with on-and offset times t = 10 s and t = 210 s for all protocols and same amplitude for all input currents (Fig 7A-7C, see Methods for further details).
The ECS potential, ECS K + concentration, and ECS volume fraction reach periodic steady states for both the slow and ultraslow stimuli (Fig 7D, 7G, 7J, 7E, 7H and 7K, respectively). Interestingly, we observe phase shifts in ϕ e , α e , and [K + ] e compared to the input flux. The phase shifts are more pronounced for the slow than the ultraslow stimuli. As expected, when j K input increases, ϕ e and α e decrease, while [K + ] e increases on average. The phase of [K + ] e is shifted by 10% and 5.0% (relative to the cycle period and the input wave) for the slow and ultraslow stimuli, respectively, with the corresponding numbers being 60% and 53% for ϕ e and 74% and 63% for α e (Fig 7A-7B, 7D-7E, 7G-7H and 7J-7K). The amplitudes of ϕ e , [K + ] e , and α e are smaller for the slow stimulus (Fig 7D, 7G and 7J) than the ultraslow (Fig 7E, 7H and 7K) and constant (Fig 7F, 7I and 7L) stimuli. For the ultraslow stimulus, the amplitude of ϕ e is the same as for the constant stimulus (Fig 7E  and 7F), whereas the amplitudes of [K + ] e and α e are smaller compared to the constant stimulus (Fig 7H-7I and 7K-7L). These reduced amplitudes are a consequence of the lower total input flux for the pulsatile stimuli, with a more substantial reduction for the slow stimulus.
The slow and ultraslow stimuli also cause periodic variations in the superficial fluid velocities, as opposed to the constant input, which causes the maximum superficial fluid velocity to plateau at 13 μm/min (Fig 8). The slow stimulus induces maximum superficial fluid velocities varying between 4.1 μm/min and 9.1 μm/min (Fig 8D), whereas the ultraslow stimulus causes the maximum superficial fluid velocity to vary between 0.39 μm/min and 13 μm/min (Fig 8E). The maximum superficial fluid velocity peaks near nadir input for the slow stimulus (Fig 8D), and near peak input for the ultraslow stimulus (Fig 8E).

Flow sensitivity to changes in permeabilities, stiffness, and input flux density
The magnitude of the compartmental fluid velocities α � u is likely to depend on the choice of parameters. We thus perform a sensitivity analysis where we measure the maximum superficial fluid velocity for different values of the compartmental permeability κ, the membrane stiffness K m , and the membrane water permeability η m (Fig 9). Specifically, we compare the sensitivity min for M1 and from 0.0 μm/min to 15 μm/min for M3 (Fig 9A). The superficial fluid velocity in M3 converges around κ = 1.4 � 10 −13 m 2 /(Pas), while the superficial fluid velocity in M1 continues to increase through κ = 1 � 10 −12 m 2 /(Pas). The difference between M1 and M3 in superficial fluid velocities decreases for κ values above the default parameter choice, but we still observe a notable difference of 215% at κ = 1 � 10 −12 m 2 /(Pas). For κ values below the default parameter choice, the absolute difference between M1 and M3 fluid velocities decreases, but the relative difference increases. Increasing K m from 0.0 Pa to 3000 Pa (with the default value set to 2294 Pa) increases the superficial fluid velocity linearly from 0.0 μm/min to 0.40 μm/min for M1 (Fig 9B). For M3, the superficial fluid velocity is 13 μm/min for all values of K m . The small changes in M1 fluid velocities and constant M3 fluid velocity lead to a close-to-constant absolute difference between the predictions made by the two modeling setups when we vary K m . By increasing the membrane water permeability η m from 0.0m/(Pas) to 1.628 � 10 −13 m/ (Pas) (with the default value set to 8.14 � 10 −14 m/(Pas)), the superficial fluid velocity increases from 0.0 μm/min to 0.32 μm/min for M1 and from 0.0 μm/min to 14 μm/min for M3 (Fig 9C). The difference in superficial fluid velocity between M1 and M3 decreases for η m values below the default choice, but we still observe that M3 predicts a 16 times higher superficial fluid velocity than M1 (4.2 μm/min vs. 0.27 μm/min) for η m as low as 5.43 � 10 −15 m/(Pas). To study how changes in the input flux density j K input affect the superficial fluid velocities, we run the simulations using 20%, 40%, 60%, 80%, and 100% of the default input flux value (8.0 � 10 −7 mol/ (m 2 s) for M1 and 9.05 � 10 −7 mol/(m 2 s) for M3). Both modeling setups show a linear relationship between j K input and the superficial fluid velocities (Fig 9D). The absolute difference between M1 and M3 decreases when j K input decreases, but M3 still predicts a superficial fluid velocity that is 47 times higher than the prediction made by M1 when the input flux density is set to 20% of the default value.

Discussion
Our results demonstrate that localized extracellular K + influx in conjunction with Na + efflux, reflecting a zone of high neuronal activity, induces a strongly coupled and complex chemicalelectrical-mechanical response in astrocyte ICS and the ECS with spatial and temporal changes in osmolarities, swelling, electrical potentials, pressures, and fluid flow. Cellular swelling and, importantly, osmotically driven fluid flow within the astrocytic network contribute to preventing high levels of extracellular water potential, effectively protecting against neuronal swelling. Fluid flow within each compartment may reach tens of μm/min and, as such, substantially contributes to the overall dynamics. Compartmental fluid flow, in concert with cellular swelling, alleviates osmotic pressure build-up and accelerates ionic transport within astrocytic networks by a factor of ×1-5 compared to diffusion alone.
The shifts in ECS ion concentrations are in line with experimental observations and reports from comparable modeling studies. Experimentally, ECS K + is measured to increase by 6-12 mM during sensory stimulation and strong electrical stimulation [2]. Halnes et al. [34,47], Østby et al. [35], and Saetra et al. [48] report an increase in ECS K + in the range of 5-10 mM via in-silico studies. We observe an increase in ECS K + of 6.68 mM in the input zone during stimuli. Further, we observe a decrease in both ECS Na + and Cl − concentrations, which is in agreement with previous modeling studies [34,47]. Interestingly, we observe a decrease in the ICS Na + concentration: although the number of Na + ions increases, the increase in ICS volume results in a total decrease in concentration. Our findings suggest that both ICS Cl − and K + concentrations increase during neuronal activity. This is in agreement with previous experimental reports [49] and with K + uptake in astrocytes facilitating clearance of excess interstitial K + following neuronal activity [2].
Astrocytes swell in response to K + influx, reducing the interstitial space volume by up to 30% [42,50,51]. The local interplay between astrocytic K + uptake and ECS shrinkage has previously been studied computationally using single neuron-glia-ECS unit models, not considering spatial buffering [35,36]. Østby et al. [35] report that inclusion of the co-transporters NBC and NKCC1, together with NKCC1-dependent water transport, is necessary to obtain ECS shrinkage that matches experimentally observed values. In their model scenario including only the major basic membrane processes (i.e., Na + /K + pump, passive ion transport, and osmotically driven water transport), they obtain an ECS shrinkage of 10.8 ± 4.0%. In contrast, Jin et al. [36] report that their model accounts for experimental observations without including non-AQP4 water transport pathways. We observe a 25% shrinkage of the ECS with only passive transmembrane water transport during stimuli comparable to that of Østby et al. [35] and Jin et al. [36] While the importance of osmotic effects has been widely recognized in the context of interstitial fluid flow and production, it has remained and remains hard to quantify [3,16]. When prescribing a hydrostatic pressure difference of 1 mmHg/mm but ignoring electrochemical contributions and interactions between the ECS and ICS, Holter et al. [27] arrived at a superficial ECS fluid velocity estimate of less than 0.  [56], as reported by Nicholson [57], the average interstitial bulk velocities in humans of 1-10 μm/min as quantified by Vinje et al. [32], and in the lower range of the 7-50 μm/min bulk flow velocities identified by Ray et al [58]. Comparing with the pioneering modeling study by Asgari et al. [37], they report baseline flow estimates resulting from a hydrostatic pressure difference (of unknown origin) alone of �1-3 � 10 −2 μm 3 /s. Interestingly, the intra-astrocytic and extracellular fluid velocities induced in our study result directly from the chemomechanical interactions following extracellular K + influx.
Part of our interest in interstitial fluid flow and the interplay between electrochemical activity and fluid dynamics relate to the intriguing effects of sleep on brain signaling, brain solute transport, and clearance [1,[59][60][61]. Both sleep spindles (11-16 Hz) and slow oscillations, i.e., low frequency (<1Hz) waves of neural activity, may support memory formation and learning [1,62]. Fultz et al. [1] used multi-modal imaging to reveal coherent patterns of slow neural-, hemodynamical-and cerebrospinal fluid waves and suggest that these waves contribute to enhance brain clearance during sleep. Our results indicate that slow and ultraslow local waves of extracellular K + and Na + flux can induce pulsatile flow of interstitial fluid in the extracellular space as well as in astrocyte networks. The slow and ultraslow regimes express different electrochemical and fluid flow wave characteristics, including differences in wave amplitudes and phase shifts.
It is more challenging to compare the absolute and relative hydrostatic pressures obtained here with clinical or experimental measurements. The coupling between electro-chemical and mechanical effects leads to an inherent cascade in which ionic concentration differences induce an osmotic pressure difference across the membrane, transmembrane water flux and cellular swelling, and intracompartmental fluid flow. The difference p i − p e between the hydrostatic pressures p i and p e is regulated by the intracellular volume changes resulting from transmembrane water movement, modulated by the elastic stiffness K m of the cell membrane. On the other hand, the absolute value of these hydrostatic pressures are determined by the incompressibility of the fluid environment and the permeability of each compartment. Within this paradigm, our simulations suggest that the ion dynamics can induce localized differences in extracellular hydrostatic pressures of several kPa over a distance of 150 μm, which would correspond to an average spatial gradient of tens of MPa/m. These values are out of range when compared with e.g. the intracranial pressure (ICP), which pulsates with the cardiac and respiratory cycles with (supine) mean ICP values of *7-15 mmHg relative to baseline atmospheric pressure (1 mmHg = *133 Pa) in healthy subjects [63,64], mean ICP wave amplitudes of *1.5-7 mmHg [65,66], and spatial differences of less than 1-3 mmHg/m [67,68]; but more comparable with normal cerebral perfusion pressure (representing the difference between the mean arterial pressure and the ICP) of 50-150 mmHg (6.7-20 kPa). In comparison, an osmotic pressure of 1 kPa (7.5 mmHg) across a cellular membrane corresponds to an osmotic concentration difference of only *0.4 mM. Future modeling work is required to couple e.g. vascular or perivascular pressure pulsations with the interstitial dynamics presented here.
Fluid flow may enhance ion and solute transport by advection in addition to the intrinsically present diffusion. We here discuss the advective contribution to K + transport through astrocytic networks. Previous computational studies of spatial K + buffering has utilized either models based on cable theory, where contributions from diffusive currents are neglected [69][70][71][72], or applied electrodiffusive frameworks accounting for the coupling of spatiotemporal variations in ion concentrations and electrical potentials [34,73]. In contrast to the model presented here, neither of these previous models account for cellular swelling and advective transport. We find that K + is mainly transported through the ICS, and notably that electrical drift dominates both diffusive and advective transport. Interestingly, we observe a net transport of K + away from the input zone even in model scenario M1, where diffusion drives K + in the opposite direction (i.e., towards the input zone). The strong dominance of electrical drift in ICS K + transport is in accordance with the findings reported in Halnes et al. [34] and Zhu et al. [39]. Further, our findings indicate that osmotically driven flow through the astrocytic syncytium facilitates spatial buffering via advection: our model scenario M3 predicts a 21.5% higher K + transport (away from the input zone) than Halnes et al. [34].
Although the exact mechanisms for water transport across astrocytic membranes are debated [3], it is well established that astrocytes have higher water permeability than neurons in part since neurons do not express AQP4 [74]. The high astrocytic water permeability has been hypothesized to stabilize extracellular ion concentrations, shielding neurons from severe swelling [34]. Our findings indeed indicate that cellular swelling, and importantly osmotic transport within the astrocytic network, facilitate in preventing water potential build-up in the ECS and thus neuronal swelling. Further, we find that models for ionic electrodiffusion must account for fluid dynamics and cellular swelling to estimate osmotic pressure across glial-ECS membranes adequately.
The computational model considered here is complex, with numerous model parameters, giving rise to considerable uncertainty. Our sensitivity analysis shows that the maximum superficial fluid velocity varies substantially under variations in a selection of the parameters related to mechanics (compartmental permeability, membrane stiffness, and membrane water permeability). Still, the differences we observe between M1 and M3 are robust to the choice of these model parameters. As such, we deem the current model useful for pointing at a mechanistic understanding of how astrocytic response to neuronal activity may impact fluid movement in the brain.
In order to maintain a reasonable level of complexity, our model distinctly includes electrodiffusion, osmosis, and hydrostatic pressures, while the representation of a number of other mechanisms are substantially simplified. In particular, we let the transmembrane water transport be parameterized by a single transmembrane water permeability (η m ), assuming that fluid is carried by passive transporters only. The parameter η m captures the permeability of all (passive) membrane fluid transporters lumped together. However, through which channels fluid flows and in what direction is still an unresolved issue [3] that could be studied further using the here proposed model as a starting point. Furthermore, we have not allowed fluid to enter or leave the system (i.e., closed boundary conditions). A natural next step for this modeling work would be to open up the boundaries to investigate the coupling between the brain tissue/neuropil and e.g. perivascular spaces. By allowing for fluid and ionic flow across the domain boundaries, one may model and explore the interplay between vascular or perivascular pressure pulsations at different frequencies and amplitudes, electrochemical activity, and interstitial fluid flow. In particular, this approach could be used to study the polarization of AQP4 channels, which are known to be densely packed at the astrocyte endfeet [75].
We conclude that the framework presented here is a promising tool for predicting complex phenomena related to electro-chemical-mechanical interplay in brain tissue. We point at a mechanistic understanding of how astrocytic response to neuronal activity and permeabilities may impact fluid movement in the brain. Our sensitivity analysis supports the idea that reduced glial water permeability may reduce (ICS and) ECS fluid velocities. Further in-silico studies with more physiological transmembrane water mechanisms and boundary conditions representing vascular or perivascular pressure pulsations could elucidate open questions related to the role of AQP4 and astrocytes in brain water clearance and homeostasis.

Methods
The homogenized tissue is represented by a one-dimensional domain O, with length L and outer boundary Γ = @O. We assume that the tissue consists of two compartments representing the ICS (denoted by subscript r = i) and the ECS (denoted by r = e). We predict the evolution in time and distribution in space of the volume fractions α r , the ion concentrations [Na + ] r , [K + ] r , and [Cl − ] r , the electrical potentials ϕ r , and the hydrostatic pressures p r . The model is embedded in the electrodiffusive Kirchhoff-Nernst-Planck framework and builds on previous work on ionic electrodiffusion [34,38], fluid dynamics [39], and astrocyte modeling [34].

Governing equations
We consider the following system of coupled, time-dependent, nonlinear partial differential equations. Find the ICS volume fraction α i : O × (0, T] ! [0, 1) such that for each t 2 (0, T]: where u i : O � ð0; T� ! R (m/s) is the ICS fluid velocity field. The transmembrane water flux w m is driven by osmotic and hydrostatic pressure, and will be discussed further in Section 4.3. The coefficient γ m (1/m) represents the area of cell membrane per unit volume of tissue. By definition, the total volume fractions sum to 1, and we assume that neurons occupy 40% of the total tissue volume [34,76]. We thus have that: Further, for each ion species k 2 {Na + , K + , Cl − } and for r = {i, e}, find the ion concentration ½k� r : O � ð0; T� ! R and the electrical potential � r : O � ð0; T� ! R such that for each t 2 (0, T]: where j k i ¼ j k i ðx; tÞ and j k e ¼ j k e ðx; tÞ (mol/(m 2 s)) are the compartmental ionic flux densities for each ion species k. Modeling of the transmembrane ion flux density j k m , the input ion flux density j k input , and the decay ion flux density j k decay will be discussed further in Sections 4.4-4.5. Note that (3) follows from first principles and express conservation of ion concentrations in each region. Moreover, we assume that the ion flux densities satisfy: where z k (unitless) is the valence of ion species k. Note that (4) arises from assuming electroneutrality, i.e., that the sum of all charge inside a compartment is zero. Given the smallness of the capacitance and that we do not model action potentials, this is a well-established approximation to use in place of the charge-capacitor relation that is commonly used when modeling electrodiffusion [38]. We further assume that the compartmental ionic flux densities j k r : O � ð0; T� ! R are driven by diffusion, electric drift, and advection: Here, D k (m 2 /s) denotes the diffusion coefficient of ion species k and λ r (unitless) denotes the tortuosity of compartment r. The constant ψ = RT/F combines Faraday's constant F (C/ mol), the absolute temperature T (K), and the gas constant R (J/(molK)).
We now turn to the dynamics of the fluid velocity fields u r and the hydrostatic pressures p r : O � ð0; T� ! R (Pa). We will consider three different models (M1, M2, and M3) for the compartmental fluid velocities u i and u e that are detailed in Section 4.2. The relation between the intra-and extracellular hydrostatic pressure p i and p e is given by the force balance on the membrane [38,39]: here K m (Pa) denotes the membrane stiffness, α i,init (unitless) denotes the initial intracellular volume fraction, and p m,init (Pa) is the initial hydrostatic pressure difference across the membrane (see Table 2). Furthermore, we assume that the volume-fraction weighted fluid velocity is divergence free; that is: By inserting (6) and the relevant expressions for the compartmental fluid velocities u i and u e (see Section 4.2 for further details) into (7), we obtain an equation for the extracellular hydrostatic pressure.
The combination of (1), (3), (4), and (7) with insertion of (5) and (6), and the expressions for u r (c.f. Section 4.2) define a system of 10 differential equations for the 10 unknown fields (α i , [k] r , ϕ r , and p e ). Note that the extracellular volume fraction α e and the intracellular hydrostatic pressure p i can be calculated using respectively (2) and (6). Appropriate initial conditions, boundary conditions, and importantly membrane mechanisms close the system.

Expressions for fluid velocities (model scenarios M1, M2, and M3)
To model compartmental fluid flow, we consider three different modeling setups:

M1
We assume that the compartmental fluid flow is driven by hydrostatic pressure gradients.
Specifically, the fluid velocities are given by: where κ r (m 2 /(Pa s)) denotes the mobility of compartment r. We will refer to κ r as the compartmental permeability.

M2
We assume that fluid flow in the glial network is driven by hydrostatic and osmotic pressure gradients. Since Na + , K + , and Cl − can move through the ICS via gap junctions, we assume that they do not contribute to osmosis. Thus, only the immobile ions drive the osmotic flow. As osmotic forces only act across membranes, we assume that fluid flow in the ECS is only driven by hydrostatic pressure gradients. The fluid velocities are given by: where i (unitless) is the Van't Hoff factor, and a i (mol/m 3 ) is the concentration of immobile ions.

M3
We follow Zhu et al. 2021 [39] and assume that the intracellular fluid flow is driven by hydrostatic and osmotic forces, while the extracellular fluid flow is driven by hydrostatic and electro-osmotic forces. To model the electro-osmotic flow, we use the Helmholtz-Smoluchowski approximation [41]. The fluid velocities are given by: Here, � r (unitless) is the relative permittivity of the extracellular solution, � 0 (F/m) is the vacuum permittivity, z (V) is the zeta-potential, and μ (Pas) is the viscosity of water.

Transmembrane fluid flow
The transmembrane fluid flow is driven by hydrostatic and osmotic pressure gradients, and the fluid velocity, w m (m/s), is expressed as: Here, η m (m/(Pas)) is the membrane water permeability and O r is the osmolarity of compartment r. We assume that all ion species contribute to the osmolarity, which is given by Multiplying O r by −iRT gives us the solute potential in compartment r, P r .

Membrane mechanisms
We adopt the ionic membrane mechanisms from Halnes et al. 2013 [34]. The mechanisms include a Na + leak channel, a Cl − leak channel, an inward-rectifying K + channel, and a Na + /K + pump. The membrane flux densities (mol/(m 2 s)) are given by: where � g k (S/m 2 ) is the membrane conductance for ion species k, ϕ m (V) is the membrane potential (defined as ϕ i − ϕ e ), and E k (V) is the reversal potential. The reversal potentials are given by the Nernst equation: Further, the Kir-function, f Kir , which modifies the inward-rectifying K + current, is given by: Here, Δϕ = ϕ m − E k , and E K,init is the reversal potential for K + at initial ion concentrations. Finally, the pump flux density (mol/(m 2 s)) is given by: where ρ pump (mol/(m 2 s)) is the maximum pump rate, P Nai (mol/m 3 ) is the [Na + ] i threshold, and P Ke (mol/m 3 ) is the [K + ] e threshold.

Input/decay fluxes
To stimulate the system, we follow the same procedure as in Halnes et al. 2013 [34]. We assume that there is a group of highly active neurons within the input zone, defined to be the interval [L 1 , L 2 ] with L 1 = 1.35 � 10 −4 m and L 2 = 1.65 � 10 −4 m, during the time interval [10 s, 210 s]. The neurons are not modeled explicitly, but we mimic their activity by injecting a K + current into the ECS and removing the same amount of Na + ions simultaneously. The input flux j k input : ½L 1 ; L 2 � � ½10 s; 210 s� ! R is given by three different protocols: where j in (mol/(m 2 s)) is constant. The constant stimulus is used for all simulations shown throughout this paper, except for the simulations illustrated by dark green and blue lines in Figs 7 and 8, which are run using the slow (1 Hz) and ultraslow (0.05 Hz) stimuli.
To mimic neuronal pumps and cotransporters, we remove K + ions from the extracellular space at a given decay rate and add the same amount of Na + ions. The decay is proportional to the extracellular K + concentration and defined across the whole domain. Specifically, the decay current j k decay : O � ð0; T� ! R (mol/(m 2 s)) is given by: where k dec denotes the [K + ] e decay factor, and [K + ] e,init the initial extracellular K + concentration.

Model parameters
All parameters of the system are as listed in Table 2 Fig 1D, Fig S6] report the Young's modulus E of individual astrocytes in the hippocampus for different deformation frequencies with E � 300 Pa for 30 Hz, E � 420 Pa for 100 Hz, and E � 520 Pa for 200 Hz, and a Poisson's ratio ν = 0.47. We use the mean Young's modulus (413 Pa) to calculate the bulk modulus (stiffness constant) K m as: ¼ 2:294 � 10 3 Pa:

Boundary conditions
We apply sealed-end boundary conditions to the system, that is, no ions and no fluid are allowed to enter or leave the system on the boundary Γ: where n Γ is the outward pointing normal vector. The extracellular electrical potential, ϕ e , and the extracellular hydrostatic pressure, p e , are only determined up to constants. To constrain the electrical potential, we require that: We enforce this zero-average constraint by introducing an additional unknown (a Lagrange multiplier) c e . For the extracellular hydrostatic pressure, we set

Initial conditions
We obtain initial conditions for the system through a two-step procedure. First, we specify a set of pre-calibrated initial values (Table 3, Pre-calibrated column). Second, we calibrate the model by running a simulation for 1 � 10 6 s. For the calibration, we set the transmembrane water permeability to zero and use N = 100 and Δt = 10 −2 . The final values from the calibration are listed in Table 3 (Post-calibrated column).
To ensure fluid equilibrium at t = 0 s and electroneutrality of the system, we define a set of immobile macromolecules, a i and a e (mol/m 3 ), with charge number z 0 (unitless) based on the initial ion concentrations. These are defined as constant concentrations with respect to the total volume of the system. Requiring an electroneutral system and zero transmembrane fluid flux at t = 0 s gives: Solving (21) gives the following expressions for z 0 , a i , and a e : To ensure strict electroneutrality and fluid equilibrium, we calculate the values at the beginning of each simulation. By using the post-calibrated initial conditions listed in Table 3, the values are approx. z 0 = −0.6, a e = 4.7 mM, and a i = 73.5 mM. Note that z 0 is interpreted as the average charge number of the macromolecules, and may thus be a decimal number (e.g., if not all the immobile macromolecules added to the system are charged).

Numerical implementation and verification
We discretize the system using a finite element method in space with characteristic mesh size Δx and a first order implicit finite difference scheme in time with time step Δt. Our numerical scheme and implementation builds on previous work presented in Ellingsrud et al. 2021 [40]. The numerical scheme is implemented via the FEniCS finite element library [79] (Python 3.8), and the code is openly available at https://github.com/martejulie/fluid-flow-in-astrocytenetworks. To determine the time step Δt and the mesh size Δx, we perform a numerical convergence study. Specifically, we apply the numerical scheme outlined above for model scenario M1 for different mesh resolutions and time steps: Δx = L/N for N = 25, 50, 100, 200, 400, 800 and Δt = 1, 10 −1 , 10 −2 , 10 −3 , 10 −4 s. For each simulation, we calculate the peak extracellular superficial fluid velocity and the peak concentration of ECS K + at t = 20 s. We find that the peak extracellular superficial fluid velocity and the ECS K + concentration converge towards 0.271 μm/min and 9.186 mM, respectively, as the temporal and spatial resolution increase (Table 4). Based on our findings, we choose Δt = 10 −3 s and N = 400 for all simulations presented within this work.

Calculation of advection/diffusion and advection/drift fractions
We calculate the advection/diffusion fraction F diff and the advection/electric-drift fraction F drift for ion species k as follows: where j k adv (mol/(m 2 s)), j k diff (mol/(m 2 s)), and j k drift (mol/(m 2 s)) are the advective, diffusive, and electric-drift components of the peak total ionic flux at t = 200 s, respectively.