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.
Significance Statement
Movement of ions and molecules in and between cellular compartments is fundamental for brain function. Cortical spreading depression (CSD) is associated with dramatic failure of brain ion homeostasis. Better understanding the sequence of events in CSD could thus provide new insight into physiological processes in the brain. Despite extensive experimental research over the last decades, even basic questions related to mechanisms underlying CSD remain unanswered. Computational modeling can play an important role going forward, since simulation studies can address hypotheses that are difficult to target experimentally. Here, we assess the physiological validity of a novel mathematical framework for detailed modeling of brain electrodiffusion and osmosis, and provide a platform for in silico studies of CSD and other cerebral electromechanical phenomena.
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+ 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+ 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 Kir 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+, 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 Kir 4.1 expression results in reduced glial swelling and depolarization of the glial membrane during CSD.
Materials and Methods
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 α, the electrical potential ϕ, and the concentrations
Model overview. The tissue is represented as a 1D domain of length 10 mm including neurons, ECS, and glial cells. Within each compartment, the model describes the dynamics of the volume fraction (α), the Na+, K+, Cl– and glutamate concentrations (
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, voltage-gated K+ and Na+ channels, NMDA receptors, and the Na+/K+-ATPase. Across the glia-extracellular interface, we model leak channels, the inward rectifying K+ channel Kir 4.1, the Na/K/Cl co-transporter, and the Na+/K+-ATPase. The precise mathematical formulation is detailed in Tuttle et al. (2019), and for the sake of completeness included in Extended Data 1.
Extended Data 1
Supplementary material. Download Extended Data 1, TEX file.
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 (Richter et al., 2002; Rosic et al., 2019), or by topical application of K+ (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+/K+-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+/K+-ATPase activity (Pietrobon and Moskowitz, 2014). All these scenarios result in a local extracellular K+ 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+, K+, and Cl– is introduced over the neuronal membrane at the first (left-most) 0.02 mm of the computational domain for the first 2 s of simulation time.
Topical application of K+: the initial values for the extracellular K+ and Cl– concentrations are increased in the first 1 mm of the computational domain.
Disabled Na+/K+-ATPase: Na+/K+-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+ 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 xi at which the neuron potential peaks at time ti, we define the wave speed vi as
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+, Cl–, and glutamate: these are defined similarly as for the DC shift, with lower thresholds of 8 mm for K+ 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). 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+ 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 Kir 4.1 channels (Djukic et al., 2007; Enger et al., 2015; Ohno, 2018).
Overview of the computational models with parameter values: χg: glial gap junction factor, γne: neuronal membrane area-to-volume, γge: glial membrane area-to-volume, ηge: glial membrane water permeability,
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+ through spatial buffering in the hippocampus (Wallraff et al., 2006). In our computational models, the glial gap junction factor χ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 (χ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 χ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 (
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 η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 ηge is set to zero, we let ηge be reduced by 25%, 50%, and 75%.
Reduced Kir 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 Kir 4.1 resting conductance of the glial membrane
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 Δx = 1.25 μm and Δx = 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).
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., 2015, 2017; Yao 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+ 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 (
Code accessibility
The simulation software described in the paper is freely available online at https://bitbucket.org/adaje/supplementary-code-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.
Extended Data 2
Code. Download Extended Data 2, ZIP file.
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 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+ of 76.4 mm (Fig. 2B,F), and glutamate of 1.38 mm (Fig. 2A,E), and decreases in extracellular Na+ and Cl– concentrations (Fig. 2B,F). The increased levels of extracellular K+ and glutamate persist for 26 and 20 s, respectively, whereas the drop in extracellular Cl– lasts 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 vi used to calculate the mean wave propagation speed, varies between 5.78 and 5.85 mm/min.
Summary of computational quantities of interest for different models (A, B, C, D, E)
Simulated CSD wave in the WT model (A) triggered by excitatory fluxes. The upper panels display a snapshot in time (plot of the respective field vs spatial coordinate x) of ECS glutamate (A), ECS ion concentrations (B), potentials (C), and change in volume fractions (D) at 60 s. The lower panels display time evolution of ECS glutamate (E), ECS ion concentrations (F), potentials (G), and change in volume fractions (H) at x = 2.0 mm.
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+ and disabled Na+/K+-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 χ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).
Comparison of Model A (solid) and Model B (stippled) CSD wave. The upper panels display snapshots (plot of field vs spatial coordinate x) of ECS glutamate (A), ECS ion concentrations (B), potentials (C), and change in volume fractions (D) at 60 s. The lower panels display time evolution of ECS glutamate (E), ECS ion concentrations (F), potentials (G) and change in volume fractions (H) evaluated at x = 2.0 mm.
Comparison of electrical potentials during simulated CSD where the glial gap junction factor χg is reduced by 0% (Model A), 25%, 50%, 75%, and 100% (Model B). The panels display neuronal potentials (A), glial potentials (B), and ECS potentials (C) evaluated at x = 2.0 mm.
Reduced membrane area-to-volume reduces CSD wave speed and amplitudes
Different values for the ratios of cell membrane area to unit tissue volume γ 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 γne and glial cells γge by 75% (Model C) substantially alters the CSD wave characteristics (Fig. 5). In particular, the amplitudes of the ECS K+ 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 γne and the glial cells γge (by 25%, 50%, and 75%) results in a gradual reduction in the amplitudes of the ECS K+ 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).
Comparison of Model A (solid) and Model C (stippled) CSD wave. The upper panels display snapshots (plot of field vs spatial coordinate x) of ECS glutamate (A), ECS ion concentrations (B), potentials (C), and change in volume fractions (D) at 60 s. The lower panels display time evolution of ECS glutamate (E), ECS ion concentrations (F), potentials (G), and change in volume fractions (H) evaluated at x = 2.0 mm.
Comparison of simulated CSD where the neuronal and glia membrane area-to-volume factors γne and γge are reduced by 0% (Model A), 25%, 50%, and 75% (Model C). The panels display ECS glutamate concentrations (A), ECS potassium concentrations (B), ECS potentials (C), change in ECS volume fractions (D), neuronal potentials (E), glial potentials (F), change in neuronal volume fractions (G), and change in glial volume fractions (H) evaluated at x = 2.0 mm.
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.
Summary and overview of experimentally reported propagation speeds, DC shift amplitudes (DC) and their duration (DC dur.), peak in extracellular K+ levels (Peak
Comparison of intervals (boxes) and min/max-ranges (whiskers) for experimentally and computationally measured values for relevant quantities in WT mice. Red lines indicate median values, and open circles denote outliers. The panels display the mean wave propagation speed (A), amplitude (B), and duration (C) of the DC shift, amplitude of elevated extracellular potassium (D), duration of elevated extracellular glutamate concentration (E), amplitude (F), and duration (G) of neuronal swelling, and amplitude of extracellular shrinkage (H).
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] s, excluding measurements at half maximum amplitude) overlap with the computational interval ([41, 68] s; Fig. 7C), and the results are thus in good agreement.
The computational values for the peak in extracellular K+ concentration overlap with the experimental reports (min/max-range of [54.4, 76.46] vs [34.3, 75.54] mm; Fig. 7D). Experimentally, elevated glutamate levels are typically indicated via the relative change in mean fluorescence (ΔF/F) over time (Enger et al., 2015, 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 ΔF/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.
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).
Comparison of changes in volume fractions during simulated CSD where the glial water permeability ηge is reduced by 0% (Model A), 25%, 50%, 75%, and 100% (Model D). The panel displays the largest (over time) change in neuronal, glial, and ECS volume fractions relative to baseline at x = 2 mm for each of the reductions.
Experimentally, Thrane et al. (2013), Yao et al. (2015), and Enger et al. (2017) have studied CSD in AQP4 knock-out mice. Enger et al. (2017) observe no difference in the CSD wave propagation speeds between the WT and AQP4 knock-outs (4.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 ± 1.7 s in WT vs 15.7 ± 1.2 s in AQP4 knock-outs). Moreover, Thrane et al. (2013) report that the amplitude of elevated levels of extracellular K+ 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 Kir 4.1 expression changes glial and DC shift dynamics
Reducing the Kir 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 Kir 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 Kir 4.1 also results in a slightly higher ECS K+ 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
Comparison of Model A (solid) and Model E (stippled) CSD wave. The upper panels display snapshots in time of ECS glutamate (A), ECS ion concentrations (B), potentials (C), and change in volume fractions (D) at 60 s. The lower panels display time evolution of ECS glutamate (E), ECS ion concentrations (F), potentials (G), and change in volume fractions (H) evaluated at x = 2.0 mm.
Comparison of glial dynamics without stimuli (A) and during induced CSD (B, C) where the Kir 4.1 expression is reduced by 0% (Model A), 10%, 20%, and 30% (Model E), 40%, 50%, 60%, 70%, and 80%. The panels display the resting glial potential (A), the glial potential during simulated CSD (B), and the change in glial volume fractions during simulated CSD (C) evaluated at x = 2.0 mm. Asterisks (*) indicate that CSD could not be induced successfully in the computational models.
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+ 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 Kir 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, 2002; Shapiro, 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 Kir 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+, depending on the values of the gap junction and Kir 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. (2008, 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+ 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
Spatial K+ buffering by astrocytes is hypothesized to be a key mechanism in controlling extracellular K+ levels. Astrocytes buffer and redistribute extracellular K+ via the astrocytic networks by transferring K+ ions from regions with high K+ concentration to regions with lower K+ levels (Orkand et al., 1966; Kofuji and Newman, 2004). There are open questions related to the role of spatial K+ buffering in CSD. Spatial buffering may prevent buildup of extracellular K+, or it could potentially increase the CSD wave speed through K+ distribution in the direction of the wave propagation. The influx of K+ during spatial buffering is mainly mediated by the inward rectifying K+ channel Kir 4.1 (Ohno, 2018; Rasmussen et al., 2020). In general, studying the role of Kir 4.1 channels in vivo is challenging: the lack of Kir 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 Kir 4.1 knock-outs computationally could potentially aid in investigating the role of Kir 4.1. Our findings show that a reduced Kir 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; Sætra 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 voltage-dependent and leak K+ channels (Olsen, 2012), while our mathematical model only accounts for Kir 4.1. The lack of other K+ model mechanisms limits the computational range of Kir 4.1 permeabilities: reductions further than 30% leads to numerical instabilities. We hypothesize that a further reduction of the Kir 4.1 expression could be obtained by including other glial K+ membrane channels (e.g., K+ leak channels). For further computational studies of Kir 4.1 dynamics, it would indeed be advantageous to reduce the Kir 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+/K+-ATPase activity or Kir 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.
Acknowledgments
Acknowledgements: We thank the vibrant discussions and insights offered by Prof. Erlend A. Nagelhus, his contributions were essential in the early phases of this work.
Footnotes
The authors declare no competing financial interests.
This work was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programme Grant 714892, the EU Project (Horizon 2020) HBP SGA3 Grant 945539, Research Council of Norway Grants 249988 and 202326, the Olav Thon Foundation, and the Letten Foundation.
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.
References
In this issue
