Validating a Computational Framework for Ionic Electrodiffusion with Cortical Spreading Depression as a Case Study

Abstract Cortical spreading depression (CSD) is a wave of pronounced depolarization of brain tissue accompanied by substantial shifts in ionic concentrations and cellular swelling. Here, we validate a computational framework for modeling electrical potentials, ionic movement, and cellular swelling in brain tissue during CSD. We consider different model variations representing wild-type (WT) or knock-out/knock-down mice and systematically compare the numerical results with reports from a selection of experimental studies. We find that the data for several CSD hallmarks obtained computationally, including wave propagation speed, direct current shift duration, peak in extracellular K+ concentration as well as a pronounced shrinkage of extracellular space (ECS) are well in line with what has previously been observed experimentally. Further, we assess how key model parameters including cellular diffusivity, structural ratios, membrane water and/or K+ permeabilities affect the set of CSD characteristics.


Introduction
Cortical spreading depression (CSD) is a slowly propagating wave of depolarization of brain cells followed by temporary silencing of electrical brain activity because of a complete collapse of cellular ion homeostasis (Pietrobon and Moskowitz, 2014). CSD is characterized by elevated levels of extracellular K 1 and glutamate (Somjen, 2001), swelling of neuronal somata and dendrites (Takano et al., 2007), swelling of astrocyte endfeet (Rosic et al., 2019), and pronounced shrinkage of the extracellular space (ECS; Mazel et al., 2002). Analyzing the sequence of events in CSD may provide new insight into physiological processes underlying both normal brain function and pathophysiological processes pertinent to a range of brain disorders (Enger et al., 2017).
Despite extensive research over the last decades, even basic questions relating to the mechanisms underlying CSD remain unanswered (Miura et al., 2007). Computational, or in silico, modeling may play an important role going forward, not least since simulation studies can address hypotheses and points of debate that are difficult to isolate or address experimentally. Here, we consider a comprehensive computational framework (the Mori framework), describing spatial and temporal dynamics of intracellular and extracellular ion concentrations, electric potentials, and volume fractions (Mori, 2015). The framework has previously been applied to study the roles of glial cells, NMDA receptors and glutamate propagation in CSD (O'Connell and Mori, 2016;Tuttle et al., 2019). However, as simulations have not to any significant extent been compared with experimental findings, the physiological validity of the computational framework remains an open question.
To systematically address this issue, we here simulate CSD in different model scenarios, and compare the computational predictions with values from the experimental literature. Our scenarios mimic different mouse models with varying structural and functional parameters, including varying intracellular diffusion, varying transmembrane water and K 1 permeabilities, and varying membrane characteristics. These choices of model scenarios are in part motivated by the incomplete or disparate findings on the role of AQP4 (Thrane et al., 2013;Yao et al., 2015;Enger et al., 2017;Rosic et al., 2019), and K ir 4.1 channels (Djukic et al., 2007;Enger et al., 2015;Ohno, 2018), both expressed in the glial cell membranes, in CSD.
Overall, we find that the range of wave speeds, direct current (DC) shift durations, peak in extracellular K 1 , neuronal changes in volume fraction, and ECS shrinkage obtained computationally all overlap with the experimentally reported ranges in wild-type (WT) mice. Further, the intracellular glial diffusivity strongly influences the DC shift, while the ratio of neuronal and glial membrane area-to-tissue volume strongly affects the CSD wave speed. Reducing the glial water permeability has a pronounced effect on cellular swelling, whereas the CSD depolarization wavefront speed and the other quantities of interest remain unaltered. In addition, we find that reducing the K ir 4.1 expression results in reduced glial swelling and depolarization of the glial membrane during CSD.

Mathematical and computational framework
We model ionic electrodiffusion and osmotic water flow in brain tissue via the Mori framework, as introduced by Mori (2015), studied numerically by Ellingsrud et al. (2021), and applied by O'Connell and Mori (2016) and Tuttle et al. (2019). This framework describes tissue dynamics in an arbitrary number of cellular compartments and the ECS via coupled ordinary differential equations and partial differential equations in a model domain. Here, we consider a version of the general framework presented by Mori (2015): we use the set of equations and physical parameters introduced by Tuttle et al. (2019). Specifically, we consider three compartments representing neurons, glial cells and the ECS (see Fig. 1), and induce CSD computationally in a one-dimensional model domain of length 10 mm. In particular, we consider the glial subtype astrocyte (for convenience, we will henceforth use the terms glia and astrocyte interchangeably). The model predicts the evolution in time and distribution in space of the volume fraction a, the electrical potential f , and the concentrations ½Na 1 ; ½K 1 ; ½Cl À ; ½Glu in each of these compartments. The membrane potential is defined as the difference between the intracellular and extracellular potentials, both for the neurons and the glial cells.
Interaction between the compartments is modelled by exchange across the neuron-extracellular and glia-extracellular membranes. In terms of osmotic water flow, the water permeability and area-to-volume ratio of the cellular membranes are key model parameters. To account for transmembrane ion movement across the neuron-extracellular membrane, we consider leak channels, voltagegated K 1 and Na 1 channels, NMDA receptors, and the Na 1 /K 1 -ATPase. Across the glia-extracellular interface, we model leak channels, the inward rectifying K 1 channel K ir 4.1, the Na/K/Cl co-transporter, and the Na 1 /K 1 -ATPase. The precise mathematical formulation is detailed in Tuttle et al. (2019), and for the sake of completeness included in Extended Data 1.

Triggers of CSD and spreading depolarization (SD)
CSD can be triggered in various ways experimentally: by electrical stimulation on the surface of the cerebral cortex (Leao, 1944), pinprick Rosic et al., 2019), or by topical application of K 1 (Yao et al., 2015;Enger et al., 2017). The consequent wave of SD is followed by a period of depressed spontaneous neuronal activity. Notably, such spontaneous activity may already be extinguished in metabolically compromised tissue (Dreier, 2011;Lauritzen et al., 2011). Therefore, the term SD is used to denote the phenomenon in the setting of hypoxia and ischemia (Somjen, 2001;Dreier, 2011). As opposed to in CSD, the SD associated conditions, such as traumatic brain injury, ischemic stroke and intracranial hemorrhage, are associated with ATP depletion and a consequent failure of the Na 1 /K 1 -ATPase (Lauritzen et al., 2011;Cozzolino et al., 2018). One notable difference between CSD and SD is the time course of recovery: the restoration of ionic gradients, repolarization of the tissue and recovery of brain function typically take longer in SD than in CSD, depending on the degree of local metabolic compromise and impairment of Na 1 /K 1 -ATPase activity (Pietrobon and Moskowitz, 2014). All these scenarios result in a local extracellular K 1 increase, which in turn causes opening of voltage gated cation channels. Here, we consider three different mechanisms for triggering CSD and SD. Specifically, the first two mechanisms trigger CSD, the latter SD: • Excitatory fluxes: a flux of Na 1 , K 1 , and Clis introduced over the neuronal membrane at the first (leftmost) 0.02 mm of the computational domain for the first 2 s of simulation time.
• Topical application of K 1 : the initial values for the extracellular K 1 and Clconcentrations are increased in the first 1 mm of the computational domain.
• Disabled Na 1 /K 1 -ATPase: Na 1 /K 1 -ATPase is disabled by setting the neuron and glia maximum pump rates to zero in the first 1 mm of the computational domain for the first 2 s of simulation time.
The precise expressions are detailed in Extended Data 1.

Quantities of interest
Experimental studies of CSD have reported on the speed of CSD waves, the duration and amplitude of extracellular K 1 and glutamate rises, the DC shift, neuronal and glial swelling, as well as extracellular shrinkage. To compare computational results to the experimental findings, we define the following quantities of interest.
• Mean wave propagation speed (mm/min): given the point x i at which the neuron potential peaks at time t i , we define the wave speed v i as v i ¼ ðx i À x iÀ1 Þ=ðt i À t iÀ1 Þ at all times t i for which the neuron potential f n ðx i Þ has passed a depolarization threshold of -20 mV and after the wave is fully initiated. We then set the mean wave propagation speed v as the average of the v i .
• Duration and amplitude of the DC shift: we define the DC shift in terms of the change in extracellular potential as follows. The amplitude of the DC shift is the maximal spatial difference in the extracellular potential (sampled at t = 60 s). The duration of the DC shift is the difference between the latest and earliest time for which the extracellular potential is below a threshold of -0.05 mV from baseline (sampled at x = 1 mm).
• Duration and amplitude of the neuronal, glial and extracellular swelling: these are defined analogously as the duration and amplitude of the DC shift. The lower threshold for neuronal and glial swelling and extracellular shrinkage is set at 0.5%. Shrinkage is defined as negative swelling.
• Duration and amplitude of elevated extracellular K 1 , Cl -, and glutamate: these are defined similarly as for the DC shift, with lower thresholds of 8 mM for K 1 and 002 mM for glutamate, and an upper threshold of 111 mM for Cl -.
• Duration and amplitude of neuronal and glial membrane depolarization: these are defined similarly as the above with a lower threshold of -66 mV for the neuronal membrane potential and -77 mV for the glial membrane potential.
In addition to these specific quantities of interest, we plot snapshots in time (i.e., plots of the respective fields vs the spatial coordinate x) and the time evolution of the computed fields evaluated at x = 1.0 mm.

Computational model and model variations
As a baseline, we define a WT model (Model A; see Table 1) with default model parameters (listed in Extended Data 1). Within each compartment, the model describes the dynamics of the volume fraction (a), the Na 1 , K 1 , Cland glutamate concentrations (½Na 1 ; ½K 1 ; ½Cl À ; ½Glu), and the potential (f ). Communication between the compartments occur via ionic and/or water membrane fluxes.
In addition, we consider four classes of model variations as described in the following (see also Table 1). Several parameters in the WT model (Model A) are associated with substantial uncertainty, in particular the strength of the glial gap junction coupling and the ratio of membrane area to tissue volume. As such, we present variations of the WT model (Model B and Model C), where we vary these two parameters gradually to explore the response in the CSD wave characteristics. Further, we introduce model versions where vary the glial transmembrane water and K 1 permeabilities, motivated by the incomplete or disparate findings on the role of AQP4 (Thrane et al., 2013;Yao et al., 2015;Enger et al., 2017;Rosic et al., 2019), and K ir 4.1 channels (Djukic et al., 2007;Enger et al., 2015;Ohno, 2018).

Reduced glial gap junctions
Interconnected astrocytes form syncytia (networks) by gap junctions. Intercellular transport through the astrocytic networks likely facilitate the removal of excess extracellular K 1 through spatial buffering in the hippocampus (Wallraff et al., 2006). In our computational models, the glial gap junction factor x g represents glial gap junctions and defines the effective intercellular diffusion through the astrocytic networks. To explore how diffusion through astrocytic networks affects the CSD wave, we consider a version of the WT model with no glial gap junctions (x g = 0, i.e., a 100% reduction) and as a consequence, a zero effective diffusion coefficient in the glial compartment (Model B). In addition, we examine the graded response of intercellular glial diffusion by reducing the glial gap junction factor x g by 25%, 50%, and 75%.

Reduced membrane area-to-volume
The ratio of membrane area to tissue volume (membrane area-to-volume) for the neuronal and glial compartments are averaged model parameters (g ne ; g ge ) that may be difficult to determine experimentally and that are associated with substantial uncertainty. Kager et al. (2000) report an estimated neuronal membrane area-tovolume based on measurements of a rat hippocampal CA1 neuron, and O'Connell and Mori (2016) assume the same value for the glial membrane area-to-volume. Halnes et al. (2019) independently estimate a higher neuronal and glial area-to-volume value. To explore the effect of reducing the neuronal and glial membrane area-to-volume, we consider model versions where we gradually reduce the area-to-volume parameters by 25%, 50%, and 75% (Model C).

Reduced glial membrane water permeability
To study the role and effect of water movement across the astrocytic membrane on CSD dynamics, we define a model variant (Model D) by setting the water permeability of the glial membrane h ge to 0 (a 100% reduction). This model thus stipulates that water cannot cross the glial membrane; neither via the lipid bilayer itself nor other membrane mechanisms such as AQP4 channels, VRAC, NKCC, or other co-transporters. As such, it can be viewed as an extreme case providing an upper bound on the effect of, e.g., reduced AQP4 expression on CSD wave characteristics, and allows for comparing computational predictions with the relatively large body of literature on AQP4 in CSD (Thrane et al., 2013;Yao et al., 2015;Enger et al., 2017;Rosic et al., 2019). Note that the water permeability is reduced while keeping all other model parameters constant, and as such, Model D does account for any potential physiological compensatory mechanisms. In addition to the extreme case where the glial water permeability h ge is set to zero, we let h ge be reduced by 25%, 50%, and 75%.

Reduced K ir 4.1 expression
To study the effect of potassium movement across the astrocytic membrane on the CSD dynamics, we define a model variant (Model E) by reducing the K ir 4.1 resting conductance of the glial membrane g K ir 4:1 . Specifically, we reduce g K ir 4:1 by 10%, 20%, 30% (Model E), 40%, 50%, 60%, 70%, and 80%. For K ir 4.1 resting conductance reductions by .30%, the computational (nonlinear) solver fails. Thus, Model E represent a partial K ir 4.1 knock-down. Changes in the membrane parameters affect the steady state of the models, and in response, we also consider new initial states of the system for these model variations (see Extended Data 1).

Numerical solution and verification
We apply the solution algorithm presented previously (Ellingsrud et al., 2021) and approximate the mathematical model numerically by a finite element method in space, a BDF2 scheme in time and a Strang splitting method. We use spatial and temporal discretization sizes of Dx = 1.25 mm and Dx = 3.125 ms, respectively. Numerical verification tests have been conducted to ensure convergence of solutions. The numerical error in the calculated mean wave propagation speed is estimated to be ,1.5%, whereas for the other quantities of interest we expect the numerical errors to be negligible (Ellingsrud et al., 2021). Table 1: Overview of the computational models with parameter values: vg: glial gap junction factor, cne: neuronal membrane area-to-volume, cge: glial membrane area-to-volume, gge: glial membrane water permeability, g Kir4:1 : glial Kir 4.1 resting conductance
Model A corresponds to the default parameters (given in Extended Data 1, only three significant digits included here). The dash (-) indicates no change from the default values.

Validation principles
We perform a quantitative and qualitative comparison of experimental and computational results via (a selection of the) quantities of interest listed above (Quantities of interest). We collected experimental findings from a set of recent experimental studies on CSD (Lauritzen and Hansen, 1992;Theis et al., 2003;Padmawar et al., 2005;Takano et al., 2007;Chang et al., 2010;Zhou et al., 2010;Thrane et al., 2013;Enger et al., 2015Enger et al., , 2017Yao et al., 2015;Kucharz et al., 2017). For the comparison, we define intervals of computational and experimental values for each of the quantities of interest reported in (a subset of) the experimental studies (wave speed, DC shift and duration, peak extracellular K 1 concentration, neuronal swelling and ECS shrinkage, elevated extracellular glutamate duration). The experimental ranges are defined by the 25th percentile (Q1, splits off the lowest 25% of the data from the highest 75%) and the 75th percentile (Q3, splits off the highest 25% of the data from the lowest 75%) of the mean values reported in the experimental studies. The computational ranges are defined by the 25th and 75th percentiles of the values from Models A, B, and C. To indicate variability outside the 25th and 75th percentiles, we define minimum and maximum values (the min/max-range) by the lowest and highest data points still within 1.5 IQR (where IQR = Q3 -Q1) of the lower and upper quartiles, respectively. Values outside the min/max-range are taken to be outliers. We qualitatively classify the match between computational and experimental results as follows: in good agreement if the computational and experimental intervals overlap, overlap in range if the min/max-ranges overlap but the intervals do not overlap, and not in agreement if the min/ max-ranges do not overlap.

Calculating changes in volume fractions relative to baseline
Figures and tables shown in Results display the change in volume fraction relative to baseline (Da r ¼ a r À a 0 r a 0 r Â 100%).
We remark that the volume fractions (a r , r = {n, g, e}) sum to 1 (cf. model equations in Extended Data 1), and thus that total volume is conserved in the computational model.

Code accessibility
The simulation software described in the paper is freely available online at https://bitbucket.org/adaje/supplementarycode-validating-a-computational-framework-for/src/master/. The code is available as Extended Data 2. The simulations were run in serial on a Lenovo ThinkPad Carbon X1 11th Gen 2.80 GHz Â 8 Intel Core i7-1165G7 CPU with Ubuntu 20.04 using FEniCS 2019.1.0, a software library enabling automated solution of (partial) differential equations.

Results
Excitatory fluxes trigger wave of depolarization, ionic changes, and swelling In the WT computational model (Model A), the excitatory fluxes trigger a wave of neuronal and glial depolarization, changes in ionic concentrations and cellular swelling Research Article: Methods/New Tools spreading through the tissue domain (Fig. 2). We observe a depolarization of the neuronal and glial potentials from -68.5 to -15.5 and -82.0 to -38.9 mV, respectively, accompanied by a DC shift with an amplitude of 11.02 mV, with a duration of 32 s (Fig. 2C,G). The neuronal depolarization wave is followed by an increase in the concentrations of extracellular K 1 of 76.4 mM (Fig. 2B,F), and glutamate of 1.38 mM ( Fig.  2A,E), and decreases in extracellular Na 1 and Clconcentrations (Fig. 2B,F). The increased levels of extracellular K 1 and glutamate persist for 26 and 20 s, respectively, whereas the drop in extracellular Cllasts for 84 s. In response to the ionic shifts, both the neurons and the glial cells swell: we observe a neuronal swelling of 11.7% and a glial swelling of 7.13%; the ECS shrinks correspondingly. We observe altered volume fractions for 103-135 s (Table 2). Finally, the wave front has reached 5.25 mm after 60 s, and we observe a mean wave propagation speed of 5.84 mm/ min. We note that the wave speeds v i used to calculate the mean wave propagation speed, varies between 5.78 and 5.85 mm/min.

Different computational CSD and SD triggers give comparable wave characteristics
CSD and SD can be triggered by different mechanisms experimentally, and we find that the same holds computationally. All three triggering mechanisms considered here, excitatory fluxes, topical application of K 1 and disabled Na 1 /K 1 -ATPase (see Materials and Methods), induce a propagating SD wave with nearly identical wave characteristics, including a mean wave propagation speed of 5.84, 5.82, and 5.83 mm/min, respectively.

Reduced glial intercellular diffusion reduces DC shift but maintains membrane potentials
The glial gap junction factor x g modulates the intercellular diffusion through the astrocytic networks. Reducing the effective intercellular diffusion (Model B) does not lead to substantial changes in the CSD wave speed, ionic concentration changes, or cellular swelling (Table 2, Model A vs B; Fig. 3C,G): Model B gives a 5% reduction in the mean wave propagation speed (to 5.52 mm/min), and ,2% in the other (ionic concentration or cellular swelling) quantities of interest. On the other hand, the DC shifts notably differ: the DC shift amplitude of Model B (3.80 mV) is 66% smaller than in Model A (11.02 mV). Yet, we observe that both the glial and neuronal potentials depolarize more in Model B compared with Model A, and thus the membrane potentials do not differ substantially between the two models. In the case of 25%, 50%, and 75% reductions in the effective intercellular diffusion, we observe a similar behavior: the DC shift amplitude gradually decreases as the glial gap junction factor is decreased, whereas the neuronal and glial depolarizations gradually increase (Fig. 4).

Reduced membrane area-to-volume reduces CSD wave speed and amplitudes
Different values for the ratios of cell membrane area to unit tissue volume g in brain tissue have been reported in the literature (Kager et al., 2000;Halnes et al., 2019). These parameters are thus uncertain and it is key to understand their effect on CSD wave characteristics.
Reducing the membrane area-to-volume for the neurons g ne and glial cells g ge by 75% (Model C) substantially alters the CSD wave characteristics (Fig. 5). In particular, the amplitudes of the ECS K 1 and glutamate elevations are reduced by 29% and 86%, respectively. We observe that both neurons and glial cells swell less, and correspondingly the ECS shrinks less. Further, the amplitudes of the neuronal and the glial membrane potentials are, respectively, 18% and 16% smaller in the Model C than in the Model A (Table 2). Remarkably, the reduced membrane area-to-volume substantially slows down the CSD wave: the wave speed is 45% slower in Model C (3.19 mm/min) compared with Model A (5.84 mm/min; Table 2). Gradually reducing the membrane area-to-volume for the neurons g ne and the glial cells g ge (by 25%, 50%, and 75%) results in a gradual reduction in the amplitudes of the ECS K 1 and glutamate elevations, the neuronal and glial swelling, the amplitudes of the neuronal and the glial membrane potentials, and the wave propagation speed (Fig. 6).

Computational versus experimental quantities of interest in WT mice
To evaluate the computational predictions, we compare our simulation results from the WT models (A, B, C) with findings from a selection of experimental studies (summarized in Table 3). For the comparison, we define intervals of computational and experimental values for each of the relevant quantities of interest (Fig. 7). The experimental ranges are defined by the 25th percentile and the 75th percentile of the mean values reported in the experimental    studies. The computational ranges are defined by the 25th and 75th percentiles of the values from Models A, B, and C. To indicate variability outside the 25th and 75th percentiles, we also define minimum and maximum values (the min/max-range; for details, see Materials and Methods, Validation principles). Values outside the min/ max-range are taken to be outliers. CSD wave propagation speeds in WT mice are reported in the interval [3.61, 4.4] mm/min, whereas the computational WT models give mean wave propagation speeds between 3.52 and 4.68 mm/min (Fig. 7A). For the wave speed, the computational and experimental results are thus in good agreement. As for the DC shift, the experimental values are in the min/max-range [12.44, 23] mV, whereas we observe DC shifts within a min/max-range of [3.80, 11.02] mV in the numerical simulations (Fig. 7B). We note that the lowest computational DC shift of 3.8 mV (Model B) likely is an underestimate, as we expect some diffusion through the astrocytic networks. We thus find that the computational DC shifts are not in agreement with the experimental values. Regarding the duration of the DC shift, the experimental interval ([66, 66.7] 7D). Experimentally, elevated glutamate levels are typically indicated via the relative change in mean fluorescence (DF/F) over time (Enger et al., 2015(Enger et al., , 2017. Comparing their duration, we find that elevation of extracellular glutamate in our computational models last for 19-25 s (min/max-range), whereas the duration of the increase in DF/F has been reported experimentally in a min/max-range of [18.6, 19.5] s (Fig. 7E). The ranges of elevated extracellular glutamate duration thus overlap.

A B C
Neuronal swelling is reported in the min/max-range [11, 37.1]% in the collection of experimental studies, with a duration of minimum 300 s and maximum 600 s. We observe a neuronal swelling in the min/max-range [5.69, 11.69]%, lasting for 109-183 s (Fig. 7F,G). The range of computational neuronal swelling amplitudes thus overlap with the experimental range, whereas the computational durations do not agree with experimental reports. We also note that neither Zhou et al. (2010) nor Takano et al. (2007) find that astrocytes swell significantly during CSD. In contrast, our numerical findings show glial (astrocytic) swelling within a range of [4.05, 7.13]%. Regarding extracellular shrinkage, Yao et al. (2015) report a 70.6% reduction of the ECS volume. The computational values are not in agreement: predicting a reduction in the extracellular volume fraction in a min/max-range of [19.32, 39.79]% (Fig. 7H).

Reducing the glial water permeability affects cellular swelling
When setting the glial water permeability to zero, we naturally observe no glial swelling, while the neuronal swelling increases by 27.7% compared with the baseline and the resulting shrinkage of the ECS is 6.2% smaller (Model D; Table 2). The remaining CSD wave characteristics, including the CSD wave propagation speed, do not change notably. Reducing the water permeability by 25%, 50%, or 75% results in unaltered CSD wave characteristics, notably including unaltered cellular swelling (Fig. 8). Experimentally, Thrane et al. (2013), Yao et al. (2015), and Enger et al. (2017) have studied CSD in AQP4 knockout mice. Enger et al. (2017) observe no difference in the CSD wave propagation speeds between the WT and AQP4 knock-outs (4.6 6 0.2 mm/min). This is in contrast to the findings of Yao et al. (2015), who report that the CSD wave propagation speed is reduced by 22% in AQP4 knock-outs. Yao et al. (2015) and Thrane et al. (2013) also report a small, but significant, reduction in the DC shift amplitude in AQP4 knock-outs compared with WT mice, whereas Enger et al. (2017) observe no significant difference in either the duration or the amplitude of the DC shift between AQP4 -/and WT. Enger et al. (2017) additionally study extracellular glutamate elevations during CSD. They report a 20% reduction in the duration of glutamate elevation for AQP4 knock-outs (18.6 6 1.7 s in WT vs 15.7 6 1.2 s in AQP4 knock-outs). Moreover, Thrane et al. (2013) report that the amplitude of elevated levels of extracellular K 1 is significantly lower in AQP4 -/than in WT mice. Our computational findings (modulo the duration of glutamate elevation) are thus in agreement with the experimental results by Enger et al. (2017).

Reduced K ir 4.1 expression changes glial and DC shift dynamics
Reducing the K ir 4.1 channel conductivity by 30% (Model E) alters the CSD wave characteristics, inducing changes in the glial potential, glial swelling and DC shift (Table 2; Fig. 9). In particular, we observe that the glial membrane potential amplitude drops from 55.14 to 48.23 when the K ir 4.1 expression is reduced. Moreover, we observe a 6% reduction in glial swelling amplitude, a 1% increase in neuronal swelling amplitude and a corresponding decrease in the ECS shrinkage amplitude. Reducing K ir 4.1 also results in a slightly higher ECS K 1 amplitude (0.2%), a slightly lower ECS glutamate amplitude (1%), and a reduced wave speed (9.8%). The computational (nonlinear) solver diverges when CSD is attempted induced in variants for which g K ir 4:1 is reduced by 40% or more. However, we observe that the new steady state value for the glial potential increases as we reduce g K ir 4:1 (Fig. 10A, no CSD induced). An 80% reduction in g K ir 4:1 results in a 63% increase in the glia resting potential (from -82 to -30 mV; Fig. 10A).

Discussion
In this study, we have simulated CSD in computational models including multiple triggering mechanisms under variations in morphologic properties, intercellular diffusivity, membrane water permeabilities, and channel conductance, and compared the computational findings with experimental reports. The range of wave speeds, DC shifts duration, peak extracellular K 1 concentration, neuronal changes in volume fraction, and extracellular shrinkage obtained computationally all overlap with the experimentally reported ranges. Our findings show that intercellular glial diffusion strongly affects the DC shift, and that the ratio of cellular membrane area to tissue volume strongly affects the CSD wave speed. The computational model predicts that a reduced K ir 4.1 expression will reduce the glial swelling and the depolarization of the glial membrane during CSD.
While several papers consider mathematical and computational modeling of CSD (Kager et al., 2000(Kager et al., , 2002Shapiro, 2001;Almeida et al., 2004;Florence et al., 2009;Wei et al., 2014;Mori, 2015;O'Connell and Mori, 2016;Tuttle et al., 2019), few have performed an systematic comparison with experimental literature. The effects of varying glial gap junction strength and the K ir 4.1 conductance on the DC shift and CSD wave speed has previously been computationally assessed by O'Connell and Mori (2016). They find that the glial gap junction coupling impacts the DC shift, which is in agreement with our computational results. Further, they find that the glial cell either buffer or broadcast K 1 , depending on the values of the gap junction and K ir 4.1 conductance. To our knowledge, there are no computational studies examining the effects of reduced glial water permeability or exploring the effect of membrane area-to-volume ratios during CSD.
There are several experimental studies on the role of AQP4 in CSD (Thrane et al., 2013;Yao et al., 2015;Enger et al., 2017;Rosic et al., 2019). Astroglial transmembrane water transport may also occur via other membrane mechanisms (e.g., EAAT, NKCC, VRAC;MacAulay et al., 2004;Pedersen et al., 2016), and the individual contributions to the total transmembrane water permeability of these different mechanisms are debated. Instead of targeting specific water mechanisms, we have here studied the effect of reducing the total glial transmembrane water permeability. When eliminating the glial membrane water transport entirely (Model D), we find, as expected, that the glial cells do not swell at all, whereas the neuronal swelling increases in comparison with the baseline model (Model A). We do not find any difference in CSD wave speed nor in any other computational quantity of interest between these two models. We note that the experimental literature is inconclusive with regard to this point. Our computational findings are largely in agreement with the experimental findings of Enger et al. (2017), knocking out AQP4 does not substantially alter the CSD wave characteristics. In contrast, Yao et al. (2015) and Thrane et al. (2013) report significant differences in their AQP4 -/mice. It is however clear that AQP4 knock-out mice may have additional morphologic differences, e.g., related to the ECS volume: Yao et al. (2008Yao et al. ( , 2015) report a ;30% larger ECS baseline volume in AQP4 -/than in WT mice. Our findings suggest that if a reduction in CSD wave speed in AQP4 -/is present, it originates from other factors than the water permeability of the glial membrane. In particular, the ratio between glial and neuronal membrane area and tissue volume was found to be an important factor for the CSD wave speed.
Diffusion of extracellular K 1 has been hypothesized to be an underlying mechanisms in CSD wave propagation (Grafstein, 1956;Obrenovitch and Zilkha, 1995;Pietrobon and Moskowitz, 2014). Local elevations in ½K 1 e will, via diffusion, increase the levels of ½K 1 e in neighboring regions, activating voltage-gated and/or ½K 1 e -dependent channels which leads to a further depolarization of the membrane and thus further release of K 1 into the ECS. Our models predict extracellular K 1 amplitudes and wave speeds that are in line with studies performed without the dampening effects of anesthesia; thus, we conjecture that these observations in our computational results may be connected. Indeed, in the model with reduced membrane area-to-volume, we observe lower extracellular K 1 amplitudes and lower wave speeds. Spatial K 1 buffering by astrocytes is hypothesized to be a key mechanism in controlling extracellular K 1 levels. Astrocytes buffer and redistribute extracellular K 1 via the astrocytic networks by transferring K 1 ions from regions with high K 1 concentration to regions with lower K 1 levels (Orkand et al., 1966;Kofuji and Newman, 2004). There are open questions related to the role of spatial K 1 buffering in CSD. Spatial buffering may prevent buildup of extracellular K 1 , or it could potentially increase the CSD wave speed through K 1 distribution in the direction of the wave propagation. The influx of K 1 during spatial buffering is mainly mediated by the inward rectifying K 1 channel K ir 4.1 (Ohno, 2018;Rasmussen et al., 2020). In general, studying the role of K ir 4.1 channels in vivo is challenging: the lack of K ir 4.1 has been reported to cause hyperexcitable neurons, epileptic seizures, and premature death (within three to four weeks) in mice (Djukic et al., 2007). Modelling K ir 4.1 knock-outs computationally could potentially aid in investigating the role of K ir 4.1. Our findings show that a reduced K ir 4.1 expression will reduce glial swelling and the depolarization of the glial membrane during CSD. We also observe a slightly reduced CSD wave propagation speed. Further, glial buffering currents have been suggested as an important mechanism behind DC potentials in brain tissue (Herreras and Makarova, 2020;Saetra et al., 2021). This is in line with our observation that a reduction in glial intercellular diffusion reduces the DC shift: reducing glial intercellular diffusion will effectively remove the glial buffering currents.
In terms of limitations, the model framework is founded on an homogenized representation of the tissue: the ECS, the cell membrane and the intracellular space are assumed to exist everywhere in the computational domain. For the spatial and temporal scales involved in the propagation of CSD, such a representation seems appropriate. Next, our simulations are based on a one-dimensional computational domain. As CSD waves spread through the whole depth of the cortical tissue, waves in two or three dimensions would be a more accurate representation. We would expect the computational CSD wave speeds to be reduced in center-initiated two-dimensional simulations. We also note that experimental studies often measure relative fluorescence intensity instead of glutamate concentrations directly, while substantial uncertainty is associated with those that do (Fabricius et al., 1993;Scheller et al., 2000). We have therefore not compared this model quantity with experimental data.
The astroglial membrane expresses a variety of voltagedependent and leak K 1 channels (Olsen, 2012), while our mathematical model only accounts for K ir 4.1. The lack of other K 1 model mechanisms limits the computational range of K ir 4.1 permeabilities: reductions further than 30% leads to numerical instabilities. We hypothesize that a further reduction of the K ir 4.1 expression could be obtained by including other glial K 1 membrane channels (e.g., K 1 leak channels). For further computational studies of K ir 4.1 dynamics, it would indeed be advantageous to reduce the K ir 4.1 expression to the point where the resting glial membrane potential is around -58 mV, as observed in vivo (Chever et al., 2010). Differences in CSD DC shift duration between different sites have been noted in experiments (Theis et al., 2003;Yao et al., 2015). In our computational model, because of the homogeneity of the material parameters, we expect (and observe) little spatial variation in the computed quantities of interest. Including heterogeneity and/or anisotropy in the model parameters would allow for studying how local tissue properties effect CSD dynamics.
The brain blood supply is continuously adjusted to match the tissue demand. CSD is, however, characterized by a transient supply-demand mismatch, sometimes referred to as vascular uncoupling (Ayata and Lauritzen, 2015). Mestre et al. (2020) demonstrated vasoconstriction associated with anoxic SDs in a stroke model facilitated CSF influx along the perivascular spaces via the so-called glymphatic pathway, and that this fluid was an important source of water accumulation in the ensuing brain edema. Moreover, the hemodynamic responses in CSD varies considerably across species, experimental approach, and more (Ayata and Lauritzen, 2015). An attempt to generalize this response, by for instance mimicking a variable Na 1 / K 1 -ATPase activity or K ir 4.1 conductance because of fluctuating ATP availability, is therefore outside the scope of our computational model.
Validating complex computational models with numerous state variables and model parameters against experimental findings with considerable variability is challenging. Although the methodology suggested here captures variability between the experimental studies, it fails in capturing the uncertainty (i.e., the SD) within each individual study. Conducting rigorous meta-analysis by combining data from the individual experiments, weighted by, e.g., the studies quality, size and/or other factors, to a pooled estimate for the experimental quantities of interest would be advantageous. However, the methodological differences (e.g., recording techniques, recording sites, in vivo, awake or anesthetized, vs in vitro) between the experimental studies makes meta-analysis difficult (Greco et al., 2013). The computational model considered within this work is complex with numerous model parameters, many of which are difficult to measure experimentally, giving rise to considerable uncertainty. A thorough uncertainty and sensitivity analysis of the computational CSD model would be interesting for future work. Importantly, the model constitutes a general yet detailed framework for capturing fundamental mechanisms of the interplay between ionic movement and volume control of the ECS. As such, applications of the framework are not limited to CSD.
We conclude that the Mori framework is a promising tool for predicting complex phenomena of ionic electrodiffusion and osmotic dynamics in brain tissue. Our findings indicate that the computational model predictions qualitatively matched well with experiments when applied to CSD. As a general framework for ionic dynamics and osmosis, applications are not limited to CSD. Mathematical modeling can be useful for isolating effects and, e.g., pointing at potential confounding effects in AQP4 knock-out mice and how they differ from WT in various aspects.