Multiscale Computer Modeling of Spreading Depolarization in Brain Slices

Abstract Spreading depolarization (SD) is a slow-moving wave of neuronal depolarization accompanied by a breakdown of ion concentration homeostasis, followed by long periods of neuronal silence (spreading depression), and is associated with several neurologic conditions. We developed multiscale (ions to tissue slice) computer models of SD in brain slices using the NEURON simulator: 36,000 neurons (two voltage-gated ion channels; three leak channels; three ion exchangers/pumps) in the extracellular space (ECS) of a slice (1 mm sides, varying thicknesses) with ion (K+, Cl–, Na+) and O2 diffusion and equilibration with a surrounding bath. Glia and neurons cleared K+ from the ECS via Na+/K+ pumps. SD propagated through the slices at realistic speeds of 2–4 mm/min, which increased by as much as 50% in models incorporating the effects of hypoxia or propionate. In both cases, the speedup was mediated principally by ECS shrinkage. Our model allows us to make testable predictions, including the following: (1) SD can be inhibited by enlarging ECS volume; (2) SD velocity will be greater in areas with greater neuronal density, total neuronal volume, or larger/more dendrites; (3) SD is all-or-none: initiating K+ bolus properties have little impact on SD speed; (4) Slice thickness influences SD because of relative hypoxia in the slice core, exacerbated by SD in a pathologic cycle; and (5) SD and high neuronal spike rates will be observed in the core of the slice. Cells in the periphery of the slice near an oxygenated bath will resist SD.


Introduction
Spreading depolarization (SD) is a slow-moving (1.7-9.2 mm/min), long-lasting (minutes) wave of neuronal depolarization accompanied by a breakdown in homeostatic maintenance of intracellular and extracellular ion concentrations, and is associated with reduced neuronal activity (spreading depression; Dreier, 2011;Woitzik et al., 2013;Cozzolino et al., 2018;Newton et al., 2018). SD has been observed in a number of species, can be elicited experimentally both in vivo and in brain slices, and has been implicated in several neurologic conditions, including ischemia, migraine, traumatic brain injury (TBI), and epilepsy (Cozzolino et al., 2018). SD is difficult to detect in humans noninvasively (Drenckhahn et al., 2012;Hartings et al., 2014;Zandt et al., 2015;Hofmeijer et al., 2018), making it important to study SD in experimental preparations and computer simulation to better understand its role in human disease, and possible treatments.
SD is intimately related with hypoxia both in slice and in vivo. Ischemia plays a complex role in ischemic diseases: hypoxia can initiate SD, which then can contribute to the extent of the ischemic penumbra (Nedergaard, 1988); and SD can itself trigger ischemia. The scintillating scotomata of classical migraine is believed to be an SD wave and increases the risk of stroke (Øie et al., 2020); complicated migraine is caused by more severe ischemia, which produces more pronounced and long-lasting deficits (Santos et al., 2012). In in vivo experiments, SD also led to hypoxia (Takano et al., 2007;Piilgaard and Lauritzen, 2009). To begin examining this complexity in our simulations, we identified three types of hypoxia related to SD, comparable to different slice experiments: (1) we induced hypoxia to in turn induce depolarization, a phenomenon that has been called "hypoxic SD-like depolarization" (HSD; Balestrino et al., 1988;Aitken et al., 1991); (2) we used the "classical" SD initiation protocol (adding a K 1 bolus to a slice) under hypoxic conditions, but before HSD had initiated, to compare to SD with an oxygenated bath; and (3) we looked at how SD induced hypoxia as perfusion from the bath failed to keep up with the metabolic demand of overworked Na 1 /K 1 pumps.
In this article, we used multiscale computational modeling of SD to relate the microscopic levels of ion and O 2 diffusion, channels, and pumps to the neuronal level of cell spiking up to the macroscopic level of tissue activation patterns (Fig. 1). Our baseline model was composed of 36,000 biophysically detailed point-neurons in an extracellular space (ECS) of a square slice (1 mm sides, 400 mm thick) with O 2 perfusion and ion flux with a surrounding bath where relevant concentrations are held constant at their baseline values. We simulated SD in both perfused and hypoxic slices. Our model showed that SD speed was augmented by propionate and hypoxia and suggested that changing the ECS was the principle mechanism through which they influence SD. We predicted that SD speed changes with slice thickness because of core hypoxia and increases with the total neuronal surface area in the tissue. SD speeds in all conditions were enhanced by hypoxia. We also predicted that increasing the size of the ECS relative to the tissue will inhibit SD. Finally, we identified a depthdependent relationship with greater SD propagation through the core of the slice compared with the periphery.

Materials and Methods
We developed a tissue-scale model of SD in slices by extending the framework developed by Wei et al. (2014) from a single cell in its local microenvironment to 36,000 cells (baseline) embedded in an ECS. We used the NEURON simulator and its extracellular reaction-diffusion framework (RxD-ECS) to simulate the electrophysiology of those neurons; the exchange of ions between them and the ECS; the diffusion of ions and O 2 through the slice; and the exchange of ions and O 2 between the slice and the bath solution in which it was submerged (Newton et al., 2018). Our model is not specific to any particular brain area, as we aimed to reproduce general properties applicable to different brain regions and to different pathophysiologies.

Model
Model neurons were all point-neurons (one compartment), which included voltage-gated K 1 and Na 1 channels; Na 1 , K 1 , and Clleak channels; KCC2 (K 1 -Cl -; Payne et al., 2003) and NKCC1 (Na 1 -K 1 -Cl -) cotransporters; and Na 1 /K 1 pumps (Fig. 1), ported from the study by Wei et al. (2014). In order for cells to balance a steadystate resting membrane potential in a slice with dynamic ion and O 2 concentrations, we added a passive electrical conductance with reversal potential (E rev = -70 mV) and conductance (g = 0.1 mS/cm 2 ). Neurons were closedended cylinders (no flux from the ends). The Na 1 /K 1 pump activity was dependent on the local concentration of O 2 in the slice, rather than ATP, and Na 1 /K 1 pumps were the only consumers of O 2 in the model (Wei et al., 2014). ATP reserve in rat cells is estimated to be ;2.6 mM (Veech et al., 1979). A single-cell simulation using a Michaelis-Menten approximation for ATP dependence on O 2 (Noske et al., 2010) demonstrated that [K 1 ] ECS = 15 mM increased pump activity by ;1.8Â, which would reduce reserves to ,1% within 2 s under hypoxic conditions.
To explore the effects of surface area to volume ratio (S/V), we used RxD-ECS to independently define a neural surface entirely separated from its volume-hence not following the overall geometry of the structure. This is possible since we used the concept of fractional volumes, rather than providing actual volume-occupying structures (McDougal et al., 2013). Neuronal volume fraction (b nrn ) is defined as the ratio of total neuronal volume (Vol nrn ) to total tissue volume: b nrn ¼ Vol nrn Vol tissue (compare with a ECS , which is Vol ECS Vol tissue ; Rice and Russo-Menna, 1997). Given a chosen tissue volume, b nrn, and a total number of neurons, N nrn (based on cell density times Vol tissue ), the volume of a single neuron, vol nrn , is as follows: (note the case used: vol nrn for single cell; Vol nrn for cumulative neuronal volume). In the NEURON simulator, neural compartments are cylinders, omitting the ends. In the present case, our point-neurons are each a single cylinder, defined by length L and diameter d. Setting L = d for simplicity, the surface area is defined as S nrn ¼ p Á d 2 . The associated volume, calculated for this cylinder, was not used. We therefore used the FractionalVolume class of RxD-ECS to scale the volume of a cell with the desired S/V to vol nrn , while the cell retains its original S nrn (McDougal et al., 2013).
To establish a biologically realistic range for S/V, we analyzed morphologic reconstructions of neurons from neuromorpho.org. We used results for cells with intact soma and dendritic reconstructions in three dimensions from animals .2 weeks of age. In rat neocortex, the average S/ V was 3.4 6 1.2 mm -1 for both pyramidal cells (n = 96) and interneurons (n = 108; Vetter et al., 2001;Larkum et al., 2004;Radman et al., 2009;Meyer et al., 2010;Boudewijns et al., 2013;Kubota, 2014). Higher S/V was grossly associated with larger dendritic trees; S/V scales inversely with diameter in cylindrical structures (excluding ends), as follows: S/V = 4 d . For our baseline simulations, we used S/V = 3.0 mm -1 ; b nrn = 0.24; neuronal density = 90 Glia are not explicitly modeled, but instead were represented as a field of sinks in every ECS voxel. Cell scale: each neuron had ion channels, 2 coexchangers; Na 1 /K 1 pump (asterisk indicates ATP/O 2 dependence). Ions were well mixed within each neuron (no intracellular diffusion). Protein scale: table (right) indicates species that control the activity of the intrinsic mechanisms in neurons and in glial field. Ion scale: ions diffused between ECS voxels by Fick's law using diffusion coefficients in Table 1. s² 10 3 neurons/mm 3 , typical of neocortex (Rice and Russo-Menna, 1997;Keller et al., 2018). Simulated slices were 1 Â 1 mm and ranged in thickness from 100 to 800 mm with 45 s²10 3 -120 s²10 3 neurons/mm 3 . The baseline simulation was a 400 mm slice with 90 s²10 3 neurons/mm 3 (36 s²10 3 neurons in total). Neurons were situated randomly throughout ECS with diffusion of Na 1 ECS , K 1 ECS , Cl -ECS , and O 2 , with diffusion coefficients (D) given in Table 1. Extracellular volume fraction (a ECS : the ratio of extracellular space volume to total tissue volume) and tortuosity (l ECS : hindrance to diffusion in tissue compared with free medium) were the same for all ions. This results in a lower effective diffusion coefficient (D*) for ions but not O 2 , which diffused through the slice unhindered. Diffusion through the ECS was calculated with a voxel size of 25 Â 25 Â 25 mm. Simulated slices were submerged in simulated bath solution where ion and O 2 concentrations were equivalent to those estimated for ECS (Table 1; Cressman et al., 2009;Wei et al., 2014) with Dirichlet boundary conditions.
In the perfused slice, [O 2 ] = 0.1 mM [corresponding to a bath solution aerated with 95% O 2 :5% CO 2 (Wei et al., 2014), a ECS = 0.2, and l ECS = 1.6 (Syková and Nicholson, 2008)]. We modeled the effects of propionate, which "primed" the tissue for SD, by reducing a ECS to 0.12 and total Clcontent in the slice by 50%, but keeping [O 2 ] = 0.1 mM and l ECS = 1.6 (Hrabe and Hrabetova, 2019). In both cases, SD was initiated by elevating initial [K 1 ] ECS within a spherical bolus at the center of the slice at t = 0. Baseline simulations were run with K 1 bolus with radius = 100 mm; [K 1 ] ECS = 70 mM. Hypoxia can also induce a depolarization, termed by some HSD, without applying a K 1 bolus (Pérez-Pinzón et al., 1995;Aitken et al., 1998). Experimentally, switching the gas to 95% N 2 :5% CO 2 leads to HSD within minutes, but immediately preceding the depolarization, Pérez-Pinzón et al. (1995) identified a period where a ECS is reduced and l ECS is increased, which they termed the "preanoxic depolarization phase." To simulate these conditions, we reduced a ECS to 0.13, increased l ECS to 1.8 (Pérez-Pinzón et al., 1995), reduced [O 2 ] in the slice to 0.01 mM, and increased K 1 to 15 mM in a 100-mm-radius sphere in the center of the slice, providing a nidus for depolarization initiation. In studies by both Pérez-Pinzón et al. (1995) and Aitken et al. (1998), the slice was reperfused immediately after detecting the depolarization, so we set the boundary conditions such that [O 2 ] in the bath was 0.1 mM (full bath oxygenation). We also simulated a scenario in which a 70 mM bolus of K 1 was applied to the slice during the preanoxic depolarization phase but before the HSD occurs, and the slice was not reperfused to make direct comparisons between SD in a hypoxic slice and SD in a perfused slice. This could be accomplished experimentally by continuously monitoring ECS properties to determine when to apply the K 1 bolus.
Studies have shown that a ECS changes dynamically over the course of SD Hrabe and Hrabetova, 2019). Since the biophysics of ECS changes during SD have not been elucidated at the timescale of our simulations, we incorporated these changes phenomenologically for a subset of simulations, reducing a ECS to 0.05-0.1 during passage of the SD  in each ECS voxel, as follows: This model of dynamic a ECS only accounts for the drop in a ECS during SD, not during its recovery after SD, which occurs on the timescale of minutes Hrabe and Hrabetova, 2019).
Glia were modeled by a background voltage-gated K 1 current and Na 1 /K 1 pump in each ECS voxel (Cressman et al., 2009;Øyehaug et al., 2012;Wei et al., 2014), rather than as individual cells.

Simulations
SD was initiated with K 1 bolus with radius = 100 mm; [K 1 ] ECS = 70 mM, unless noted otherwise. To follow the position of SD over time, we tracked the position of the K 1 wavefront, defined as the farthest location where [K 1 ] ECS exceeded 15 mM, and the positions of spiking neurons. Propagation speed was indicated as 0 if the K 1 wave did not propagate. Most simulations ran for 10 s, which was sufficient for SD to propagate throughout the entire slice. Because of computational limitations, we restricted our simulations to 10 s, which meant that we could not continue to the termination phase of SD, which would require considerably greater temporal and spatial scales.

Code accessibility
All codes for simulation, data analysis, and visualization described in this article were written in Python/NEURON and are freely available on GitHub (Extended Data 1). In the course of this study, we ran .600 simulations covering a range of slice sizes, cell densities, and durations on a number of different architectures. Simulating a 1 Â 1 Â 400 mm slice with a cell density of 90,000 neurons/mm 3 (36 s²10 3 neurons) for 10 s of simulation time took ;12.5 h on a parallel Linux system using 48 nodes on a 2.40 GHz Intel Xeon E5-4610 CPU. Incorporating dynamic a ECS into the same model on the same machine increased simulation time to ;18 h. Simulations were run using Neuroscience Gateway (Sivagnanam et al., 2013), Google Cloud Platform, and the on-site high-performance computer at SUNY Downstate Health Sciences University. Smaller versions of the model can be run on a personal computer.

Code availability
The model is publicly available on GitHub and ModelDB. The access code for ModelDB is Airy1870.

Results
In our model of an O 2 -perfused slice, a small bolus of elevated K 1 (70 mM, 100 mm radius) initiated a propagating K 1 wave with associated SD producing neuronal spiking (Fig. 2). The K 1 wave traveled radially outward from the bolus in three dimensions toward the edges of the slice at 2.3 mm/min, comparable to optical and electrophysiological measurements of SD propagation velocity in brain slices (Aitken et al., 1998;Joshi and Andrew, 2001;Hrabe and Hrabetova, 2019). Within the K 1 bolus, most cells fired a single spike and went immediately into depolarization block. Outside the K 1 bolus, cells fired a 200-900 ms burst of action potentials as [K 1 ] ECS increased around them. During the course of the SD-associated burst, instantaneous firing rates increased to as high as 250 Hz with decreasing spike heights during the burst, comparable to experimental observations (Devin Brisson et al., 2013;Lemaire et al., 2021). Cells then remained in depolarization block for the remainder of the 10 s measured (see Materials and Methods for computational limitations; Devin Brisson et al., 2013;Andrew, 2016). Spreading depolarization, seen intracellularly, produced Na 1 channel inactivation and prevented further spiking. The absence of spiking would be seen extracellularly as spreading depression-a silent region in the slice. The K 1 wave and SD were coincident in time and space, with spreading depression following closely behind; we primarily followed the K 1 wave since this was easiest to localize across conditions. SD in perfused slices was an all-or-none process; it could only be initiated above a certain threshold measured either in concentration-[K 1 ] ECS ! 20 mM (bolus diameter, 200 mm) or [K 1 ] ECS = 70 mM (bolus diameter ! 100 mm; Fig.  3). Beyond these thresholds, different K 1 bolus concentrations and diameters had only a minimal effect on wave speed.
Underpinning SD was a wave of pronounced imbalance of transmembrane ion concentrations (Fig. 4, Movie 1). Excess K 1 is briefly eliminated from the ECS via neural and glial homeostatic mechanisms. Once the K 1 wave arrived, however, the cells dumped large quantities of K 1 into the ECS because of the burst and subsequent prolonged depolarization (note the rapid depletion of [K 1 ] i ;    4B). The Na 1 /K 1 pump activity that contributed to K 1 elimination from the ECS in the core of the slice created a high demand for O 2 , exceeding the rate at which it could diffuse in from the bath. This resulted in much of the tissue becoming hypoxic before the arrival of the K 1 wave (Fig. 4, compare E, rapid falloff of O 2 , F, much slower rise of extracellular K 1 ). The rapid spread of O 2 deficit explains the total pump failure at intermediate locations in the slice. There were small upward deflections in the first three traces in Figure 4B reflecting homeostatic inward pumping. There was no upward deflection in the other three traces-O 2 has disappeared before the K 1 wave arrives. The preservation of pumping in the final, most peripheral, trace of Figure 4B is because of this measurement being at the edge of the slice, neighboring the O 2 source of the bath. In reality, neuronal ATP reserves will maintain pumping for a limited duration in the absence of O 2 , but reserves will quickly run out because of the high metabolic demand imposed on the pumps by large shifts in K 1 and Na 1 concentrations. Once [K 1 ] ECS reached ;14 mM, cellular homeostatic mechanisms totally broke down. Changes in intracellular and extracellular Na 1 concentrations were actually larger than the shifts in K 1 accompanying SD, and a wave of extracellular Na 1 deficit traveled along with the K 1 wave (Fig. 4C,G). Shifts in Clconcentrations proceeded Movie 1. Extracellular ion and O 2 concentrations across the slice averaged over depth, as well as neuronal spiking (white dots) from 250 neurons during the course of SD in a perfused slice. The spread of spiking and the K 1 wave can be seen in real time. We recommend downloading the file and using slower playback to visualize the spread of hypoxia. [View online] through the slice slightly more slowly and were less pronounced (Fig. 4D,H).
SD was facilitated by incorporating the effects of hypoxia or propionate treatment on the slice. HSD (caused by hypoxia, rather than introducing a K 1 bolus) propagated similarly to SD in a perfused slice, as shown experimentally (Aitken et al., 1998;Fig. 5). Pérez-Pinzón et al. (1995) identified a period immediately preceding the HSD they termed a preanoxic depolarization phase, with reduced a ECS and increased l ECS . We modeled HSD by incorporating these hypoxia effects, using a small region of slightly elevated K 1 (100 mm radius, 15 mM) to provide a nidus for HSD initiation. The spatiotemporal distribution of neuronal spiking was similar to that of the control (perfused) slice, as seen experimentally (Fig. 5B, Movie 2). However, the K 1 wave during HSD was faster (3.7 mm/min) than during SD in the perfused slice (2.3 mm/min; Fig. 5A). We also simulated standard K 1 -initiated SD in hypoxic slices by applying a 70 mM K 1 bolus (as with SD in the perfused slice) to a slice during this preanoxic depolarization phase (Pérez-Pinzón et al., 1995), resulting in K 1 wave speed of 3.4 mm/ min. Simulating propionate application (decreased a ECS = 0.12; halving [Cl -] i and [Cl -] ECS ; Tao et al., 2002;Hrabe and Hrabetova, 2019), increased K 1 wave speed to 4.8 mm/ min (Fig. 5A). Comparable speedups were also observed in the depolarization waves (Fig. 5B). Since these manipulations included combined changes to [O 2 ], Cl -, a ECS , and/or l ECS , we investigated their individual contributions over Figure 5. Hypoxia, propionate, and dynamic ECS increased SD speed principally through a ECS reduction. A, Radial K 1 wave position over time during SD in perfused (Movie 1), hypoxic, HSD (Movie 2), propionate conditions, and with dynamic changes in a ECS . Hypoxia, propionate, and dynamic changes in a ECS facilitated propagation. B, Radial position of SD wave represented by time to first spike in 126 selected cells at different distances from center. C, K 1 wave speeds with individual parameter changes (Fig. 2). a ECS had the greatest impact on SD speed over a physiologically plausible range (x-axis ranges: Movie 2. Extracellular ion and O 2 concentrations across the slice averaged over depth, as well as neuronal spiking (white dots) from 250 neurons during the course of hypoxic SD-like depolarization. The spread of spiking and the K 1 wave can be seen in real time. We recommend downloading the file and using slower playback to visualize the delay between initiating the simulation and the spread of the K 1 wave. [View online] Research Article: New Research the relevant ranges. Since both hypoxia and propionate decrease a ECS , we also tested increasing a ECS to as high as 0.42, which has been observed when making the ECS artificially hypertonic (Kume-Kick et al., 2002). a ECS had the largest influence on propagation, changing K 1 wave speed by .2 mm/min over the range tested, while K 1 wave speed changed by ,0.5 mm/min for the ranges of [O 2 ], l ECS , and [Cl -] values tested (Fig. 5C).
Experimental studies have demonstrated dynamic changes in a ECS occurring during SD, with a ECS dropping to as low as 0.05 at the peak of the depolarization Hrabe and Hrabetova, 2019). Given the strong influence of constant a ECS on SD propagation (Fig. 5C, orange line), we also explored the influence of dynamic a ECS (Fig.  5A, purple line). Dynamically decreasing a ECS was modeled as a local function of increasing [K 1 ] ECS (see Materials and Methods), such that a ECS dropped to 0.06 in the wake of SD, within the experimentally observed range of 0.05-0.1 . Incorporating a dynamic ECS increased the speed of SD propagation in perfused slice from 2.3 to 3.3 mm/min. The stereotyped characteristics of neuronal firing patterns during the depolarization wave remained unchanged by dynamic changes in a ECS . Different brain areas have different cell densities, and their neurons have different morphologic characteristics. We manipulated our generic model so as to explore the following three properties of neural tissue organization and shape: neuronal S/V; the fraction of tissue volume occupied by neurons (b nrn ); and cell density (number of neurons per mm 3 ; Fig. 6). Neuronal S/V varies across cell types, brain regions, and species. Examination of representative morphologies showed that S/V values are generally in a range of 2-10 mm -1 (see Materials and Methods), with neocortical principal cell S/V of 3.4 6 1.2 mm -1 (n = 96), which is significantly greater than the brainstem principal cell S/V of 2.2 6 1.2 mm -1 (n = 74; p ,0.001, Mann-Whitney U test; Núñez-Abades et al., 1994;Vetter et al., 2001;Larkum et al., 2004;Radman et al., 2009;Meyer et al., 2010;Ros et al., 2010;Boudewijns et al., 2013;Raslan et al., 2014;Williams et al., 2019). Neuronal volume fraction b nrn may differ with different brain areas and will differ under the pathologic condition of cytotoxic edema. Cell density varies across different neural areas.
Realistic (.1 mm -1 ) S/V was necessary for initiating SD (Fig. 6A)-SD could not be initiated using the actual 3D geometry of single-cylinder point-neurons with diameter and height selected to produce a baseline b nrn of 0.24 (S/ V = 0.02 mm -1 with 90,000 neurons/mm 3 , perfused or hypoxic). Above this threshold, K 1 wave speed increased with S/V. K 1 wave speed also increased with increased b nrn while keeping S/V and cell density (number of neurons per mm 3 ) constant (Fig. 6B). Cell density effects were less marked, whether through keeping surface area and volume constant (Fig. 6C) or keeping b nrn and S/V constant (Fig.  6D). In all cases, change with altered parameters was more pronounced in the hypoxic slice. Pooled together, we found a near linear relationship of K 1 wave speed with total neuronal surface area.
Slice thickness (100-800 mm) influenced SD by altering the ability of O 2 to penetrate to the tissue core (Fig. 7). SD could Figure 6. K 1 wave propagation speed proportional to total neuronal surface area in slice. SD initiated in perfused and hypoxic slices (no reperfusion) by introducing a 100-mm-radius central 70 mM spherical K 1 bolus. A-E, Effects of varying S/V of each cell while maintaining a constant b nrn (A); varying b nrn while keeping S/V constant, thus allowing S nrn to vary (B); varying cell density while keeping constant S nrn and vol nrn , thus allowing b nrn to vary (C); and varying cell density while keeping b nrn and S/V constant, thus allowing S nrn and vol nrn to vary (D). E, Pooled results: K 1 wave speed increased linearly with total neuronal surface area in both perfused and hypoxic slices. Hypoxia increased wave speed across conditions (0 speed indicates no SD). not be initiated in the 100 mm perfused slice-SD was not sustainable with full O 2 availability, but it was initiated in an hypoxic slice of the same thickness. With increasing thickness, an increasingly hypoxic core (despite O 2 perfusion of the bath) allowed K 1 wave speed to increase from 1.2 to 2.1 mm/min over 200-400 mm thickness (Fig.  7A). Above 400 mm, there was no increased speed with increased thickness. Similar patterns were observed in hypoxic slices.
We also observed depth-dependent differences in propagation of the SD within a slice (Fig. 8). A wave of high K 1 propagated through the core of the slice (650 mm from the center; Fig. 8A), and neurons there exhibited the membrane dynamics of SD: a burst of spikes followed by depolarization block lasting the remainder of the simulation (Fig. 8C). The K 1 wave reached the periphery of the slice (within 50 mm of bath) but in lower concentration (Fig. 8B), and neurons there were resistant to SD because Figure 7. Slice thickness effects on SD propagation. K 1 wave speed during SD in perfused and hypoxic slices of various thicknesses. SD could not be initiated in very thin (100 mm) slices when perfused but could in hypoxic slices. For perfused slices, K 1 wave speed increased with slice thickness between 200 and 400 mm then saturated for slices of greater thickness. A similar pattern was observed in the hypoxic slices with consistently faster K 1 wave speeds compared with in perfused slices.
Figure 8. Depth dependence of SD propagation in 400-mm-thick perfused slice. A, Spread of K 1 wave through slice core (650 mm from center). B, Wave of mildly elevated K 1 reached the periphery (within 50 mm of bath) from the central bolus. C, Spread of SD through the slice core. Voltages in slice core (color map) with spike rasters for 120-cell subset overlaid in white. Neurons in the core showed typical SD voltage dynamics (bursting followed by depolarization block). D, Voltages and raster at periphery; cells show regular spiking patterns at 10-70 Hz (V memb averaged across all cells in 25 Â 25 Â 100 mm voxels at center or periphery). of the high availability of O 2 from the bath (Fig. 8D). Instead, neurons in the periphery of the slice only underwent a modest depolarization and fired regularly at 10-70 Hz, without the bursting followed by depolarization block characteristic of SD.

Discussion
Our model reproduced a number of the properties of SD observed in brain slices (Figs. 2, 4). Slice models that most resembled cortical gray matter (high cell density, high neuronal S/V) showed SD speeds and neuronal firing patterns comparable to in vitro optical and electrophysiological measurements (Aitken et al., 1998;Devin Brisson et al., 2013;Hrabe and Hrabetova, 2019;Lemaire et al., 2021). The all-or-none nature of SD initiation, as well as the observed bolus [K 1 ] threshold of ;20 mM, was also in agreement with experiments in brain slices (Andrew et al., 2017).
Our simulations identify a pathologic cycle whereby SD induces, and is also exacerbated by, hypoxia in slice. In vivo, SD-associated metabolic demand for O 2 can exceed supply, resulting in tissue hypoxia (Takano et al., 2007;Piilgaard and Lauritzen, 2009). However, in vivo experiments showed hypoxia following SD, rather than preceding it as predicted by our in vitro simulations. In vivo, O 2 is supplied by the vasculature, and increased neural activity during SD may lead to increased blood supply and tissue oxygenation (Balança et al., 2017), while in vitro, O 2 is only supplied by the bath, which is effectively unchanging. While K 1 was slowly spreading outward across the slice, O 2 spread inward rapidly, following the gradient caused by O 2 depletion from overworked pumps, but was quickly consumed (Fig. 4, compare E, rapid falloff of O 2 , B, much slower falloff of [K 1 ] i ). Depending on the distance to the bath and to the inciting bolus, intermediate locations in the slice suffered various degrees of pump demand and partial pump failure. The resistance to SD of neurons in the periphery of the slice compared with those in the core (Fig. 8B) was comparable to findings in vivo: tissue near capillaries resist SD, while tissue farther away is relatively hypoxic and prone to SD (Takano et al., 2007).
By comparing the effects of changing O 2 availability, total Clcontent, a ECS , and l ECS on SD in isolation, we determined that a ECS influenced SD most strongly (Fig.  5C), accounting for most of the speedup observed in hypoxic and propionate-treated slices. Our results support the hypothesis that the main mechanism in the priming of propionate for SD is through reducing a ECS (Hrabe and Hrabetova, 2019), and suggest that the main mechanism for hypoxia speedup is also reduced a ECS .
We modeled SD in brain slices as a reaction-diffusion system in an unconnected neuronal population. Although we initiated SD by elevating extracellular K 1 and tracked the position of K 1 waves, the model itself was agnostic as to the agent of SD propagation. In our simulations of HSD, we saw that the initial elevation of extracellular K 1 did not immediately initiate a spreading wave of K 1 , but led to delayed spread (Fig. 5A, Movie 2). In our simulations, as in the slice experiment, when K 1 was applied to the tissue or hypoxia was induced, the elevation of K 1 preceded the wavefront of SD-associated neuronal spiking. By contrast, when SD has been experimentally initiated in other ways, depolarization may precede, rather than follow, the K 1 wave (Somjen, 2001). Incorporating synaptic connections into the model will provide a saltatory forward activation that will tend to speed up depolarization relative to the K 1 wave; the balance of these factors will depend strongly on cell excitability, and the densities and lengths of connections.

Model limitations
An irony of this type of study is that we are computer modeling an in vitro model of an in vivo animal model of human pathophysiology (model of a model of a model). There are necessarily distortions at each step. At the computer-modeling level, the major limitations of this study are the limitations that are inherent in all computer modeling-we necessarily made choices as to what to include and what to leave out. In the current set of simulations, we left out the following: (1) all neural connectivity; (2) types of neurons, including distinction between inhibitory and excitatory neurons; (3) dendritic morphologies plus additional membrane mechanisms; (4) glia, except as a generalized field; (5) volume-occupying structures (instead using fractional volumes); and (6) intracellular handling of ions and second messengers with effects on pumps and other mechanisms. These are largely correctable limitations that we will gradually begin to address in future versions of the model. Some, however, represent limitations in experimental knowledge that need to be addressed.
Additionally, we purposely designed the model to be generic rather than to reproduce the properties of any particular brain region and species. We feel this allowed us to generalize more readily (e.g., comparing SD in brainstem vs in cortex). Several interdependent tissue properties were treated independently with the benefit of allowing additional investigation by isolating parameters.
A major limitation of the model was the simplification in which Na 1 /K 1 pumps consumed O 2 to drive their conductances rather than ATP (Wei et al., 2014). This was partly justified by Na 1 /K 1 pumps using 91% of available ATP in the brain under normal physiological conditions (Lennie, 2003); also, neurons synthesize ATP as needed (Davis, 2020). Because of ATP reserves in the tissue and the simplification using O 2 as a proxy for ATP, the spread of hypoxia during SD in perfused slices may be slower than predicted by our simulations, but we still expect it to precede the SD wavefront in vitro.
We also note that extracellular Na 1 concentration shifts following SD were smaller than observed experimentally, and intracellular concentration shifts were larger than expected (Somjen, 2001). This may reflect biophysical mechanisms absent from the model or a lack of appropriate volume changes in ECS and intraneuronal spaces. In particular, we did not model large, intracellular, negatively charged macromolecules that produce Gibbs-Donnan effects, which contribute to SD (Lemale et al., 2022).

Experimentally testable predictions
Several of our predictions relate to the effects of manipulations on SD speed. These effects could be most easily assessed electrophysiologically by using a series of extracellular electrodes in a slice to note the time of population burst passage and subsequent time of silence (the depression phase).

Slower SD in brainstem slice compared with cortical slice
Compared with cortex, brainstem has lower cell density, higher a ECS , and low expression of ECS matrix molecules/perineural nets, implying low l ECS (Ogawa et al., 1985;Hobohm et al., 1998;Syková and Nicholson, 2008). All of these factors contribute to slower propagation speeds in our model (Figs. 5C, 6) Our analysis of principal cell morphologies from brainstem also suggested an S/V lower than those of neocortical principal cells, another factor contributing to slowing.
Increased SD speed with cytotoxic edema in penumbra after stroke or TBI Cytotoxic edema will increase b nrn , producing speedup, which will be enhanced in the hypoxic condition (Fig.  6B). We note that some of the fastest SD speeds (;9 mm/ min) have been observed in patients after stroke (Woitzik et al., 2013).
Hypoxia will increase SD speed under multiple conditions SD spread faster in hypoxic than in perfused slices. Experimentally, one would assess as follows: (1) remove O 2 from bath; (2) monitor a ECS to determine the onset of preanoxic depolarization phase (a ECS , ;0.099-0.0179, depending on brain region; Pérez-Pinzón et al., 1995); and (3) add a K 1 bolus. The slice should not be reperfused, unlike the procedure for HSD experiments (Pérez-Pinzón et al., 1995;Aitken et al., 1998).

Increasing a ECS will attenuate SD propagation
This could be assessed using hypertonic saline to increase a ECS (Kume-Kick et al., 2002). Hypertonic saline is sometimes used in reducing intracranial pressure after TBI (Oddo et al., 2009;Kamel et al., 2011;Mangat et al., 2020;Shi et al., 2020), and might therefore also reduce SD in these patients.

SD speed will be reduced by antiepileptic drugs
We showed here that dynamic changes in a ECS due to SD itself speeds up the SD wave (Fig. 5A,B). Similar changes in a ECS have been seen with ictal phenomena in the study by Colbourn et al. (2021), allowing us to hypothesize that this may be the linkage between SD and a ECS , presumably mediated through excessive release of neurotransmitters whether classical, peptidergic, or nitric oxide.
SD velocity will correlate with neural density, dendritic complexity, and total neuronal volume across regions Measurements of SD in various brain regions and across species can be assessed. Density is determined with counts in Nissl stain. Dendritic complexity increases S/V and can be assessed on traced, biocytin-filled cells with measures such as we performed here. Total neuronal volume can be assessed by measuring ascorbic acid in tissue (Rice and Russo-Menna, 1997). These effects could also be evaluated in tissue culture or in organoids.
Increasing the diameter or concentration of the K 1 bolus used to initiate SD beyond their thresholds will not change SD speed (Fig. 3).
Ease of SD initiation and SD propagation speed will increase with increasing slice thickness due to relatively hypoxic core SD will be difficult or impossible to initiate in very thin slices (Fig. 7).
SD propagates through the core of the slice Neurons near the surface of a perfused slice will be resistant to SD due to the high availability of O 2 (Fig. 8). Extracellular recordings looking for bursting and subsequent depression at different slice depths could be performed to confirm this prediction. However, one may want to avoid measurement directly at the slice surface where neurons will have suffered damage due to slicing and therefore may exhibit additional pathology that could alter SD propagation.

Future directions
Our model incorporated quantitative data and simpler models from numerous sources and at multiple spatial scales to constitute a unified framework for simulating SD in brain slices. We propose the use of this framework as a community tool for researchers in the field to test hypotheses; guide the design of new experiments; and incorporate new physiological, transcriptomic, proteomic, or anatomic data into the framework. The open-source, branchable, versioned GitHub repository can provide a dynamic database for SD simulations or modeling brain slices in general.