Research ArticleResearch Article: New Research, Disorders of the Nervous System

Multiscale Computer Modeling of Spreading Depolarization in Brain Slices

Craig Kelley, Adam J. H. Newton, Sabina Hrabetova, Robert A. McDougal and William W Lytton
eNeuro 4 August 2022, 9 (4) ENEURO.0082-22.2022; https://doi.org/10.1523/ENEURO.0082-22.2022

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.

  • extracellular space
  • hypoxia
  • NEURON
  • reaction-diffusion
  • spreading depolarization
  • spreading depression

Significance Statement

Spreading depolarization (SD) is a slow-moving wave of electrical and ionic imbalances in brain tissue and is a hallmark of several neurologic disorders. We developed a multiscale computer model of brain slices with realistic neuronal densities, ions, and oxygenation. Our model shows that SD is exacerbated by and causes hypoxia, resulting in strong SD dependence on slice thickness. Our model also predicts that the velocity of SD propagation is not dependent on its initiation, but instead on tissue properties, including the amount of extracellular space and the total area of neuronal membrane, suggesting faster SD following ischemic stroke or traumatic brain injury.

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 has been studied in brain slices from a wide range of species and brain regions, including neocortex, hippocampus, brainstem, and retina (Balestrino et al., 1988; Aitken et al., 1998; Müller and Somjen, 1998; Martins-Ferreira et al., 2000; Devin Brisson et al., 2013; Andrew, 2016; Hrabe and Hrabetova, 2019). It can be triggered experimentally by various means, including electrical stimulation, mechanical insult, and K+ and ouabain application (Leao, 1944; Bures et al., 1974; Balestrino, 1995; Aitken et al., 1998; Joshi and Andrew, 2001). SD can be facilitated by applying propionate to the slice (Tao et al., 2002; Hrabe and Hrabetova, 2019).

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+ 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+/K+ pumps.

In this article, we used multiscale computational modeling of SD to relate the microscopic levels of ion and O2 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 μm thick) with O2 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 depth-dependent relationship with greater SD propagation through the core of the slice compared with the periphery.

Figure 1.
Figure 1.

Multiscale model expanded. Tissue scale: a few of the 36 ⋅ 103 neurons (pink circles) embedded in the ECS of a brain slice submerged in a bath solution where ion and O2 concentrations were held constant. 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+/K+ pump (asterisk indicates ATP/O2 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.

View this table:
Table 1

Diffusion coefficients and baseline concentrations for ions and O2 in perfused slice (Samson et al., 2003; Haynes, 2014; Wei et al., 2014)

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 O2 through the slice; and the exchange of ions and O2 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+ and Na+ channels; Na+, K+, and Cl– leak channels; KCC2 (K+-Cl–; Payne et al., 2003) and NKCC1 (Na+-K+-Cl–) cotransporters; and Na+/K+ pumps (Fig. 1), ported from the study by Wei et al. (2014). In order for cells to balance a steady-state resting membrane potential in a slice with dynamic ion and O2 concentrations, we added a passive electrical conductance with reversal potential (Erev = –70 mV) and conductance (g = 0.1 mS/cm2). Neurons were closed-ended cylinders (no flux from the ends). The Na+/K+ pump activity was dependent on the local concentration of O2 in the slice, rather than ATP, and Na+/K+ pumps were the only consumers of O2 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 O2 (Noske et al., 2010) demonstrated that [K+]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 (βnrn) is defined as the ratio of total neuronal volume (Volnrn) to total tissue volume: βnrn=VolnrnVoltissue (compare with αECS, which is VolECSVoltissue ; Rice and Russo-Menna, 1997). Given a chosen tissue volume, βnrn, and a total number of neurons, Nnrn (based on cell density times Voltissue), the volume of a single neuron, volnrn, is as follows: volnrn=βnrn·VoltissueNnrn, (note the case used: volnrn for single cell; Volnrn 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 Snrn=π·d2 . 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 volnrn, while the cell retains its original Snrn (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 ± 1.2 μm– 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 = 4d . For our baseline simulations, we used S/V = 3.0 μm–1; βnrn = 0.24; neuronal density = 90 ⋅ 103 neurons/mm3, 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 μm with 45 ⋅103–120 ⋅103 neurons/mm3. The baseline simulation was a 400 μm slice with 90 ⋅103 neurons/mm3 (36 ⋅103 neurons in total). Neurons were situated randomly throughout ECS with diffusion of NaECS+ , KECS+ , ClECS– , and O2, with diffusion coefficients (D) given in Table 1. Extracellular volume fraction (αECS: the ratio of extracellular space volume to total tissue volume) and tortuosity (λ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 O2, which diffused through the slice unhindered. Diffusion through the ECS was calculated with a voxel size of 25 × 25 × 25 μm. Simulated slices were submerged in simulated bath solution where ion and O2 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, [O2] = 0.1 mm [corresponding to a bath solution aerated with 95% O2:5% CO2 (Wei et al., 2014), αECS = 0.2, and λECS = 1.6 (Syková and Nicholson, 2008)]. We modeled the effects of propionate, which “primed” the tissue for SD, by reducing αECS to 0.12 and total Cl– content in the slice by 50%, but keeping [O2] = 0.1 mm and λECS = 1.6 (Hrabe and Hrabetova, 2019). In both cases, SD was initiated by elevating initial [K+]ECS within a spherical bolus at the center of the slice at t = 0. Baseline simulations were run with K+ bolus with radius = 100 μm; [K+]ECS = 70 mm. Hypoxia can also induce a depolarization, termed by some HSD, without applying a K+ bolus (Pérez-Pinzón et al., 1995; Aitken et al., 1998). Experimentally, switching the gas to 95% N2:5% CO2 leads to HSD within minutes, but immediately preceding the depolarization, Pérez-Pinzón et al. (1995) identified a period where αECS is reduced and λECS is increased, which they termed the “preanoxic depolarization phase.” To simulate these conditions, we reduced αECS to 0.13, increased λECS to 1.8 (Pérez-Pinzón et al., 1995), reduced [O2] in the slice to 0.01 mm, and increased K+ to 15 mm in a 100-μm-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 [O2] in the bath was 0.1 mm (full bath oxygenation). We also simulated a scenario in which a 70 mm bolus of K+ 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+ bolus.

Studies have shown that αECS changes dynamically over the course of SD (Mazel et al., 2002; 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 αECS to 0.05–0.1 during passage of the SD (Mazel et al., 2002) in each ECS voxel, as follows: dαECSdt=(5·10−7)·([K+]ECS−3.5)·([K+]ECS−70.0)·αECS.

This model of dynamic αECS only accounts for the drop in αECS during SD, not during its recovery after SD, which occurs on the timescale of minutes (Mazel et al., 2002; Hrabe and Hrabetova, 2019).

Glia were modeled by a background voltage-gated K+ current and Na+/K+ 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+ bolus with radius = 100 μm; [K+]ECS = 70 mm, unless noted otherwise. To follow the position of SD over time, we tracked the position of the K+ wavefront, defined as the farthest location where [K+]ECS exceeded 15 mm, and the positions of spiking neurons. Propagation speed was indicated as 0 if the K+ 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 μm slice with a cell density of 90,000 neurons/mm3 (36 ⋅103 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 α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.

Extended Data 1

Simulation code: our tissue-scale model of SD in brain slices is available on GitHub. We used the NEURON simulator reaction-diffusion framework to implement embed thousands of neurons (based on the model from Wei et al., 2014) in the extracellular space of a brain slice, which is itself embedded in an bath solution. We initiated SD in the slice by elevating extracellular K+ in a spherical region at the center of the slice. The effects of hypoxia and propionate on the slice were modeled by appropriate changes to the volume fraction and tortuosity of the extracellular space and oxygen/chloride concentrations. Users need to install NEURON, and we recommend using MPI to parallelize simulations. Download Extended Data 1, ZIP file.

Code availability

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

Results

In our model of an O2-perfused slice, a small bolus of elevated K+ (70 mm, 100 μm radius) initiated a propagating K+ wave with associated SD producing neuronal spiking (Fig. 2). The K+ 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+ bolus, most cells fired a single spike and went immediately into depolarization block. Outside the K+ bolus, cells fired a 200–900 ms burst of action potentials as [K+]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+ 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+ wave and SD were coincident in time and space, with spreading depression following closely behind; we primarily followed the K+ wave since this was easiest to localize across conditions.

Figure 2.
Figure 2.

A small bolus of applied K+ initiates spreading K+ and depolarization waves in perfused slice. A, [K+]ECS averaged across slice depth (400 μm) at 4 time points during SD. B, Spike raster plot of 250 randomly selected neurons (of 36 ⋅103) during SD. Cells are ordered on y-axis by their radial distance from the center of the K+ bolus. Blank area under spikes represents region of spreading depression. Baseline values: [O2] = 0.1 mm; αECS = 0.2; λECS = 1.6; [Cl–]ECS = 130.0 mm; [Cl–]i = 6.0 mm.

SD in perfused slices was an all-or-none process; it could only be initiated above a certain threshold measured either in concentration—[K+]ECS ≥ 20 mm (bolus diameter, 200 μm) or [K+]ECS = 70 mm (bolus diameter ≥ 100 μm; Fig. 3). Beyond these thresholds, different K+ bolus concentrations and diameters had only a minimal effect on wave speed.

Figure 3.
Figure 3.

Diameter (bottom, x-axis) and concentration (top, x-axis) of the K+ bolus had minor effects on K+ wave speed. Mean and SDs (n = 5) for K+ wave speed versus bolus diameter (solid line; bolus [K+]ECS = 70 mm; 5 random cell position initializations) and versus bolus concentration (dashed line; diameter, 200 μm).

Underpinning SD was a wave of pronounced imbalance of transmembrane ion concentrations (Fig. 4, Movie 1). Excess K+ is briefly eliminated from the ECS via neural and glial homeostatic mechanisms. Once the K+ wave arrived, however, the cells dumped large quantities of K+ into the ECS because of the burst and subsequent prolonged depolarization (note the rapid depletion of [K+]i; Fig. 4B). The Na+/K+ pump activity that contributed to K+ elimination from the ECS in the core of the slice created a high demand for O2, 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+ wave (Fig. 4, compare E, rapid falloff of O2, F, much slower rise of extracellular K+). The rapid spread of O2 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—O2 has disappeared before the K+ 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 O2 source of the bath. In reality, neuronal ATP reserves will maintain pumping for a limited duration in the absence of O2, but reserves will quickly run out because of the high metabolic demand imposed on the pumps by large shifts in K+ and Na+ concentrations. Once [K+]ECS reached ∼14 mm, cellular homeostatic mechanisms totally broke down. Changes in intracellular and extracellular Na+ concentrations were actually larger than the shifts in K+ accompanying SD, and a wave of extracellular Na+ deficit traveled along with the K+ wave (Fig. 4C,G). Shifts in Cl– concentrations proceeded through the slice slightly more slowly and were less pronounced (Fig. 4D,H).

Figure 4.
Figure 4.

Concentrations at 7 radial locations measured during SD in perfused 400 μm slice. A, Cell within the K+ bolus (37.5 μm) produced a single spike; cells farther out fired a burst. Cells remained in depolarization block for the remainder of the simulation (10 s). B–D, Intracellular ion concentrations. E–H, Extracellular O2 and ion concentrations in neighboring ECS voxels. Movie 1 shows extracellular ion and O2 concentrations across the slice, as well as neuronal spiking (white dots) from 250 neurons during the course of SD (most easily seen with slowed playback).

Movie 1.

Extracellular ion and O2 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+ wave can be seen in real time. We recommend downloading the file and using slower playback to visualize the spread of hypoxia.

SD was facilitated by incorporating the effects of hypoxia or propionate treatment on the slice. HSD (caused by hypoxia, rather than introducing a K+ 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 αECS and increased λECS. We modeled HSD by incorporating these hypoxia effects, using a small region of slightly elevated K+ (100 μm 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+ 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+-initiated SD in hypoxic slices by applying a 70 mm K+ 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+ wave speed of 3.4 mm/min. Simulating propionate application (decreased αECS = 0.12; halving [Cl–]i and [Cl–]ECS; Tao et al., 2002; Hrabe and Hrabetova, 2019), increased K+ 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 [O2], Cl–, αECS, and/or λECS, we investigated their individual contributions over the relevant ranges. Since both hypoxia and propionate decrease αECS, we also tested increasing αECS to as high as 0.42, which has been observed when making the ECS artificially hypertonic (Kume-Kick et al., 2002). αECS had the largest influence on propagation, changing K+ wave speed by >2 mm/min over the range tested, while K+ wave speed changed by <0.5 mm/min for the ranges of [O2], λECS, and [Cl–] values tested (Fig. 5C).

Figure 5.
Figure 5.

Hypoxia, propionate, and dynamic ECS increased SD speed principally through αECS reduction. A, Radial K+ wave position over time during SD in perfused (Movie 1), hypoxic, HSD (Movie 2), propionate conditions, and with dynamic changes in αECS . Hypoxia, propionate, and dynamic changes in α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+ wave speeds with individual parameter changes (Fig. 2). αECS had the greatest impact on SD speed over a physiologically plausible range (x-axis ranges: [O2] = 0.01–0.1 mm; αECS = 0.07–0.42; λECS = 1.4–2.0; [Cl–]ECS:[Cl–]I = 3.0:65.0–6.0:130.0 mm).

Movie 2.

Extracellular ion and O2 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+ 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+ wave.

Experimental studies have demonstrated dynamic changes in αECS occurring during SD, with αECS dropping to as low as 0.05 at the peak of the depolarization (Mazel et al., 2002; Hrabe and Hrabetova, 2019). Given the strong influence of constant αECS on SD propagation (Fig. 5C, orange line), we also explored the influence of dynamic αECS (Fig. 5A, purple line). Dynamically decreasing αECS was modeled as a local function of increasing [K+]ECS (see Materials and Methods), such that αECS dropped to 0.06 in the wake of SD, within the experimentally observed range of 0.05–0.1 (Mazel et al., 2002). 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 α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 (βnrn); and cell density (number of neurons per mm3; 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 μm– 1 (see Materials and Methods), with neocortical principal cell S/V of 3.4 ± 1.2 μm– 1 (n = 96), which is significantly greater than the brainstem principal cell S/V of 2.2 ± 1.2 μm– 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 βnrn may differ with different brain areas and will differ under the pathologic condition of cytotoxic edema. Cell density varies across different neural areas.

Figure 6.
Figure 6.

K+ wave propagation speed proportional to total neuronal surface area in slice. SD initiated in perfused and hypoxic slices (no reperfusion) by introducing a 100-μm-radius central 70 mm spherical K+ bolus. A–E, Effects of varying S/V of each cell while maintaining a constant βnrn (A); varying βnrn while keeping S/V constant, thus allowing Snrn to vary (B); varying cell density while keeping constant Snrn and volnrn, thus allowing βnrn to vary (C); and varying cell density while keeping βnrn and S/V constant, thus allowing Snrn and volnrn to vary (D). E, Pooled results: K+ 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).

Realistic (>1 μm– 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 βnrn of 0.24 (S/V = 0.02 μm– 1 with 90,000 neurons/mm3, perfused or hypoxic). Above this threshold, K+ wave speed increased with S/V. K+ wave speed also increased with increased βnrn while keeping S/V and cell density (number of neurons per mm3) constant (Fig. 6B). Cell density effects were less marked, whether through keeping surface area and volume constant (Fig. 6C) or keeping β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+ wave speed with total neuronal surface area.

Slice thickness (100–800 μm) influenced SD by altering the ability of O2 to penetrate to the tissue core (Fig. 7). SD could not be initiated in the 100 μm perfused slice—SD was not sustainable with full O2 availability, but it was initiated in an hypoxic slice of the same thickness. With increasing thickness, an increasingly hypoxic core (despite O2 perfusion of the bath) allowed K+ wave speed to increase from 1.2 to 2.1 mm/min over 200–400 μm thickness (Fig. 7A). Above 400 μm, there was no increased speed with increased thickness. Similar patterns were observed in hypoxic slices.

Figure 7.
Figure 7.

Slice thickness effects on SD propagation. K+ wave speed during SD in perfused and hypoxic slices of various thicknesses. SD could not be initiated in very thin (100 μm) slices when perfused but could in hypoxic slices. For perfused slices, K+ wave speed increased with slice thickness between 200 and 400 μm then saturated for slices of greater thickness. A similar pattern was observed in the hypoxic slices with consistently faster K+ wave speeds compared with in perfused slices.

We also observed depth-dependent differences in propagation of the SD within a slice (Fig. 8). A wave of high K+ propagated through the core of the slice (±50 μm 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+ wave reached the periphery of the slice (within 50 μm of bath) but in lower concentration (Fig. 8B), and neurons there were resistant to SD because of the high availability of O2 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.

Figure 8.
Figure 8.

Depth dependence of SD propagation in 400-μm-thick perfused slice. A, Spread of K+ wave through slice core (±50 μm from center). B, Wave of mildly elevated K+ reached the periphery (within 50 μm 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 (Vmemb averaged across all cells in 25 × 25 × 100 μm voxels at center or periphery).

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+] 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 O2 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, O2 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, O2 is only supplied by the bath, which is effectively unchanging. While K+ was slowly spreading outward across the slice, O2 spread inward rapidly, following the gradient caused by O2 depletion from overworked pumps, but was quickly consumed (Fig. 4, compare E, rapid falloff of O2, B, much slower falloff of [K+]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 O2 availability, total Cl– content, αECS, and λECS on SD in isolation, we determined that α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 αECS (Hrabe and Hrabetova, 2019), and suggest that the main mechanism for hypoxia speedup is also reduced α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+ and tracked the position of K+ 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+ did not immediately initiate a spreading wave of K+, but led to delayed spread (Fig. 5A, Movie 2). In our simulations, as in the slice experiment, when K+ was applied to the tissue or hypoxia was induced, the elevation of K+ 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+ 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+ 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+/K+ pumps consumed O2 to drive their conductances rather than ATP (Wei et al., 2014). This was partly justified by Na+/K+ 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 O2 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+ 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 αECS, and low expression of ECS matrix molecules/perineural nets, implying low λ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 β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 O2 from bath; (2) monitor αECS to determine the onset of preanoxic depolarization phase (αECS, ∼0.099–0.0179, depending on brain region; Pérez-Pinzón et al., 1995); and (3) add a K+ 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 αECS will attenuate SD propagation

This could be assessed using hypertonic saline to increase α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 αECS due to SD itself speeds up the SD wave (Fig. 5A,B). Similar changes in α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 α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+ 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 O2 (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.

Acknowledgments

Acknowledgment: We thank Drs. Rena Orman, Mark Stewart, Richard Kollmar, and John Kubie (SUNY Downstate Health Sciences University) for useful discussions on this topic. We also thank Dr. Michael Hines (Yale University) for support with the NEURON simulator.

Footnotes

  • The authors declare no competing financial interests.

  • This work was funded by National Institutes of Health/National Institute of Mental Health Grant R01-MH-086638 and National Science Foundation Grant Internet2 E-CAS 1904444. We thank Drs. Rena Orman, Mark Stewart, Richard Kollmar, and John Kubie (SUNY Downstate Health Sciences University) for useful discussions on this topic. We also thank Dr. Michael Hines (Yale University) for support with the NEURON simulator.

  • Received February 21, 2022.
  • Revision received June 28, 2022.
  • Accepted July 11, 2022.

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

  1. Aitken PG, Jing J, Young J, Somjen GG (1991) Ion channel involvement in hypoxia-induced spreading depression in hippocampal slices. Brain Res 541:7–11. pmid:1709386
  2. Aitken PG, Tombaugh GC, Turner DA, Somjen GG (1998) Similar propagation of SD and hypoxic SD-like depolarization in rat hippocampus recorded optically and electrically. J Neurophysiol 80:1514–1521. pmid:9744955
    1. Monti MM,
    2. Sannita WG
    Andrew RD (2016) The persistent vegetative state: evidence that the lower brain survives because its neurons intrinsically resist ischemia. In: Brain function and responsiveness in disorders of consciousness (Monti MM, Sannita WG, eds), pp 119–144. Cham, Switzerland: Springer.
  4. Andrew RD, Hsieh Y-T, Brisson CD (2017) Spreading depolarization triggered by elevated potassium is weak or absent in the rodent lower brain. J Cereb Blood Flow Metab 37:1735–1747. pmid:27381828
  5. Balança B, Meiller A, Bezin L, Dreier JP, Marinesco S, Lieutaud T (2017) Altered hypermetabolic response to cortical spreading depolarizations after traumatic brain injury in rats. J Cereb Blood Flow Metab 37:1670–1686. pmid:27356551
  6. Balestrino M (1995) Pathophysiology of anoxic depolarization: new findings and a working hypothesis. J Neurosci Methods 59:99–103. doi:10.1016/0165-0270(94)00199-Q
    1. Somjen G
    Balestrino M, Aitken PG, Jones LS, Somjen GG (1988) The role of spreading Depression-Like hypoxic depolarization in irreversible neuron damage, and its prevention. In: Mechanisms of cerebral hypoxia and stroke (Somjen G, ed), pp 291–301. Boston: Springer.
  8. Boudewijns ZSRM, Groen MR, Lodder B, McMaster MTB, Kalogreades L, de Haan R, Narayanan RT, Meredith RM, Mansvelder HD, de Kock CPJ (2013) Layer-specific high-frequency action potential spiking in the prefrontal cortex of awake rats. Front Cell Neurosci 7:99. pmid:23805075
  9. Bures J, Buresová O, Kr˘ivánek J (1974) The mechanism and applications of Leão’s spreading depression of electroencephalographic activity. Prague: Academia.
  10. Colbourn R, Hrabe J, Nicholson C, Perkins M, Goodman JH, Hrabetova S (2021) Rapid volume pulsation of the extracellular space coincides with epileptiform activity in mice and depends on the NBCe1 transporter. J Physiol 599:3195–3220. doi:10.1113/JP281544 pmid:33942325
  11. Cozzolino O, Marchese M, Trovato F, Pracucci E, Ratto GM, Buzzi MG, Sicca F, Santorelli FM (2018) Understanding spreading depression from headache to sudden unexpected death. Front Neurol 9:19. pmid:29449828
  12. Cressman JR Jr., Ullah G, Ziburkus J, Schiff SJ, Barreto E (2009) The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: I. Single neuron dynamics. J Comput Neurosci 26:159–170. doi:10.1007/s10827-008-0132-4 pmid:19169801
  13. Davis GW (2020) Not fade away: mechanisms of neuronal ATP homeostasis. Neuron 105:591–593. doi:10.1016/j.neuron.2020.01.024 pmid:32078791
  14. Devin Brisson C, Lukewich MK, David Andrew R (2013) A distinct boundary between the higher brain’s susceptibility to ischemia and the lower brain’s resistance. PLoS One 8:e79589. pmid:24223181
  15. Dreier JP (2011) The role of spreading depression, spreading depolarization and spreading ischemia in neurological disease. Nat Med 17:439–447. pmid:21475241
  16. Drenckhahn C, Winkler MKL, Major S, Scheel M, Kang E-J, Pinczolits A, Grozea C, Hartings JA, Woitzik J, Dreier JP (2012) Correlates of spreading depolarization in human scalp electroencephalography. Brain 135:853–868. doi:10.1093/brain/aws010 pmid:22366798
  17. Hartings JA, Wilson JA, Hinzman JM, Pollandt S, Dreier JP, DiNapoli V, Ficker DM, Shutter LA, Andaluz N (2014) Spreading depression in continuous electroencephalography of brain trauma. Ann Neurol 76:681–694. pmid:25154587
  18. Haynes WM (2014) CRC handbook of chemistry and physics. Boca Raton, FL: CRC.
  19. Hobohm C, Härtig W, Brauer K, Brückner G (1998) Low expression of extracellular matrix components in rat brain stem regions containing modulatory aminergic neurons. J Chem Neuroanat 15:135–142. pmid:9797071
  20. Hofmeijer J, van Kaam CR, van de Werff B, Vermeer SE, Tjepkema-Cloostermans MC, van Putten MJAM (2018) Detecting cortical spreading depolarization with full band scalp electroencephalography: an illusion? Front Neurol 9:17. doi:10.3389/fneur.2018.00017 pmid:29422883
  21. Hrabe J, Hrabetova S (2019) Time-resolved integrative optical imaging of diffusion during spreading depression. Biophys J 117:1783–1794. doi:10.1016/j.bpj.2019.08.031 pmid:31542225
  22. Joshi I, Andrew RD (2001) Imaging anoxic depolarization during ischemia-like conditions in the mouse hemi-brain slice. J Neurophysiol 85:414–424. pmid:11152742
  23. Kamel H, Navi BB, Nakagawa K, Hemphill, Iii JC, Ko NU (2011) Hypertonic saline versus mannitol for the treatment of elevated intracranial pressure: a meta-analysis of randomized clinical trials*. Crit Care Med 39:554–559. pmid:21242790
  24. Keller D, Erö C, Markram H (2018) Cell densities in the mouse brain: a systematic review. Front Neuroanat 12:83. pmid:30405363
  25. Kubota Y (2014) Untangling GABAergic wiring in the cortical microcircuit. Curr Opin Neurobiol 26:7–14. pmid:24650498
  26. Kume-Kick J, Mazel T, Vorisek I, Hrabětová S, Tao L, Nicholson C (2002) Independence of extracellular tortuosity and volume fraction during osmotic challenge in rat neocortex. J Physiol 542:515–527. pmid:12122149
  27. Larkum ME, Senn W, Lüscher H-R (2004) Top-down dendritic input increases the gain of layer 5 pyramidal neurons. Cereb Cortex 14:1059–1070. doi:10.1093/cercor/bhh065 pmid:15115747
  28. Leao AAP (1944) Spreading depression of activity in the cerebral cortex. J Neurophysiol 7:359–390. doi:10.1152/jn.1944.7.6.359
  29. Lemaire L, Desroches M, Krupa M, Pizzamiglio L, Scalmani P, Mantegazza M (2021) Modeling NaV1.1/SCN1A sodium channel mutations in a microcircuit with realistic ion concentration dynamics suggests differential GABAergic mechanisms leading to hyperexcitability in epilepsy and hemiplegic migraine. PLoS Comput Biol 17:e1009239. doi:10.1371/journal.pcbi.1009239 pmid:34314446
  30. Lemale CL, Lückl J, Horst V, Reiffurth C, Major S, Hecht N, Woitzik J, Dreier JP (2022) Migraine aura, transient ischemic attacks, stroke, and dying of the brain share the same key pathophysiological process in neurons driven by Gibbs-Donnan forces, namely spreading depolarization. Front Cell Neurosci 16:837650. doi:10.3389/fncel.2022.837650
  31. Lennie P (2003) The cost of cortical computation. Curr Biol 13:493–497. doi:10.1016/s0960-9822(03)00135-0 pmid:12646132
  32. Mangat HS, Wu X, Gerber LM, Schwarz JT, Fakhar M, Murthy SB, Stieg PE, Ghajar J, Härtl R (2020) Hypertonic saline is superior to mannitol for the combined effect on intracranial pressure and cerebral perfusion pressure burdens in patients with severe traumatic brain injury. Neurosurgery 86:221–230. doi:10.1093/neuros/nyz046
  33. Martins-Ferreira H, Nedergaard M, Nicholson C (2000) Perspectives on spreading depression. Brain Res Brain Res Rev 32:215–234. pmid:10751672
  34. Mazel T, Richter F, Vargová L, Syková E (2002) Changes in extracellular space volume and geometry induced by cortical spreading depression in immature and adult rats. Physiol Res 51 [Suppl 1]:S85–S93. pmid:12479789
  35. McDougal RA, Hines ML, Lytton WW (2013) Reaction-diffusion in the NEURON simulator. Front Neuroinform 7:28. pmid:24298253
  36. Meyer HS, Wimmer VC, Hemberger M, Bruno RM, de Kock CPJ, Frick A, Sakmann B, Helmstaedter M (2010) Cell type-specific thalamic innervation in a column of rat vibrissal cortex. Cereb Cortex 20:2287–2303. doi:10.1093/cercor/bhq069 pmid:20534783
  37. Müller M, Somjen GG (1998) Inhibition of major cationic inward currents prevents spreading depression-like hypoxic depolarization in rat hippocampal tissue slices. Brain Res 812:1–13. pmid:9813218
  38. Nedergaard M (1988) Mechanisms of brain damage in focal cerebral ischemia. Acta Neurol Scand 77:81–101. pmid:3284272
  39. Newton AJH, McDougal RA, Hines ML, Lytton WW (2018) Using NEURON for reaction-diffusion modeling of extracellular dynamics. Front Neuroinform 12:41. pmid:30042670
  40. Noske R, Cornelius F, Clarke RJ (2010) Investigation of the enzymatic activity of the Na+, K+-ATPase via isothermal titration microcalorimetry. Biochim Biophys Acta 1797:1540–1545. doi:10.1016/j.bbabio.2010.03.021 pmid:20362545
  41. Núñez-Abades PA, He F, Barrionuevo G, Cameron WE (1994) Morphology of developing rat genioglossal motoneurons studied in vitro: changes in length, branching pattern, and spatial distribution of dendrites. J Comp Neurol 339:401–420. pmid:8132869
  42. Oddo M, Levine JM, Frangos S, Carrera E, Maloney-Wilensky E, Pascual JL, Kofke WA, Mayer SA, LeRoux PD (2009) Effect of mannitol and hypertonic saline on cerebral oxygenation in patients with severe traumatic brain injury and refractory intracranial hypertension. J Neurol Neurosurg Psychiatry 80:916–920. doi:10.1136/jnnp.2008.156596
  43. Ogawa Y, Eins S, Wolff JR (1985) Oligodendrocytes in the pons and middle cerebellar peduncle of the cat. topographical relations to neurons and transverse axon bundles. Cell Tissue Res 240:541–552. pmid:2410128
  44. Øie LR, Kurth T, Gulati S, Dodick DW (2020) Migraine and risk of stroke. J Neurol Neurosurg Psychiatry 91:593–604. pmid:32217787
  45. Øyehaug L, Østby I, Lloyd CM, Omholt SW, Einevoll GT (2012) Dependence of spontaneous neuronal firing and depolarisation block on astroglial membrane transport mechanisms. J Comput Neurosci 32:147–165. pmid:21667153
  46. Payne JA, Rivera C, Voipio J, Kaila K (2003) Cation-chloride co-transporters in neuronal communication, development and trauma. Trends Neurosci 26:199–206. pmid:12689771
  47. Pérez-Pinzón MA, Tao L, Nicholson C (1995) Extracellular potassium, volume fraction, and tortuosity in rat hippocampal CA1, CA3, and cortical slices during ischemia. J Neurophysiol 74:565–573. pmid:7472364
  48. Piilgaard H, Lauritzen M (2009) Persistent increase in oxygen consumption and impaired neurovascular coupling after spreading depression in rat neocortex. J Cereb Blood Flow Metab 29:1517–1527. pmid:19513087
  49. Radman T, Ramos RL, Brumberg JC, Bikson M (2009) Role of cortical cell type and morphology in subthreshold and suprathreshold uniform electric field stimulation in vitro. Brain Stimul 2:215–228.e1–3. pmid:20161507
  50. Raslan A, Ernst P, Werle M, Thieme H, Szameit K, Finkensieper M, Guntinas-Lichius O, Irintchev A (2014) Reduced cholinergic and glutamatergic synaptic input to regenerated motoneurons after facial nerve repair in rats: potential implications for recovery of motor function. Brain Struct Funct 219:891–909. doi:10.1007/s00429-013-0542-6 pmid:23543131
  51. Rice ME, Russo-Menna I (1997) Differential compartmentalization of brain ascorbate and glutathione between neurons and glia. Neuroscience 82:1213–1223. doi:10.1016/S0306-4522(97)00347-3
  52. Ros H, Magill PJ, Moss J, Bolam JP, Mena-Segovia J (2010) Distinct types of non-cholinergic pedunculopontine neurons are differentially modulated during global brain states. Neuroscience 170:78–91. doi:10.1016/j.neuroscience.2010.06.068 pmid:20603194
  53. Samson E, Marchand J, Snyder KA (2003) Calculation of ionic diffusion coefficients on the basis of migration test results. Mat Struct 36:156–165. doi:10.1007/BF02479554
  54. Santos E, Sánchez-Porras R, Dohmen C, Hertle D, Unterberg AW, Sakowitz OW (2012) Spreading depolarizations in a case of migraine-related stroke. Cephalalgia 32:433–436. doi:10.1177/0333102412441414 pmid:22407661
  55. Shi J, Tan L, Ye J, Hu L (2020) Hypertonic saline and mannitol in patients with traumatic brain injury: a systematic and meta-analysis. Medicine 99:e21655. doi:10.1097/MD.0000000000021655
  56. Sivagnanam S, Majumdar A, Yoshimoto K, Astakhov V, Bandrowski A, Martone M, Carnevale N (2013) Introducing the neuroscience gateway. Proceedings of the 5th International Workshop on Science Gateways. CEUR Workshop Proceedings.
  57. Somjen GG (2001) Mechanisms of spreading depression and hypoxic spreading depression-like depolarization. Physiol Rev 81:1065–1096. pmid:11427692
  58. Syková E, Nicholson C (2008) Diffusion in brain extracellular space. Physiol Rev 88:1277–1340. pmid:18923183
  59. Takano T, Tian G-F, Peng W, Lou N, Lovatt D, Hansen AJ, Kasischke KA, Nedergaard M (2007) Cortical spreading depression causes and coincides with tissue hypoxia. Nat Neurosci 10:754–762. pmid:17468748
  60. Tao L, Masri D, Hrabetová S, Nicholson C (2002) Light scattering in rat neocortical slices differs during spreading depression and ischemia. Brain Res 952:290–300. pmid:12376191
  61. Veech RL, Lawson J, Cornell N, Krebs HA (1979) Cytosolic phosphorylation potential. J Biol Chem 254:6538–6547. doi:10.1016/S0021-9258(18)50401-4 pmid:36399
  62. Vetter P, Roth A, Häusser M (2001) Propagation of action potentials in dendrites depends on dendritic morphology. J Neurophysiol 85:926–937. pmid:11160523
  63. Wei Y, Ullah G, Schiff SJ (2014) Unification of neuronal spikes, seizures, and spreading depression. J Neurosci 34:11733–11743. pmid:25164668
  64. Williams PA, Dalton C, Wilson CG (2019) Modeling hypoglossal motoneurons in the developing rat. Respir Physiol Neurobiol 265:40–48. pmid:30056243
  65. Woitzik J, Hecht N, Pinczolits A, Sandow N, Major S, Winkler MKL, Weber-Carstens S, Dohmen C, Graf R, Strong AJ, Dreier JP, Vajkoczy P (2013) Propagation of cortical spreading depolarization in the human cortex after malignant stroke. Neurology 80:1095–1102. doi:10.1212/WNL.0b013e3182886932 pmid:23446683
  66. Zandt B-J, ten Haken B, van Putten MJAM, Dahlem MA (2015) How does spreading depression spread? Physiology and modeling. Rev Neurosci 26:183–198. pmid:25719306

Synthesis

Reviewing Editor: William Stacey, University of Michigan

Decisions are customarily a result of the Reviewing Editor and the peer reviewers coming together and discussing their recommendations until a consensus is reached. When revisions are invited, a fact-based synthesis statement explaining their decision and outlining what is needed to prepare a revision will be listed below. The following reviewer(s) agreed to reveal their identity: Jens Dreier, Ghanim Ullah.

This paper uses a computational-only approach to predict several physiological (potentially testable in the future) variables that are involved in spreading depression. The model makes five main predictions: 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 due to relative hypoxia in the slice core, exacerbated by SD in a pathological cycle; and 5. Higher neuronal spike rates and faster spread will be observed in the core than the periphery of perfused slice during SD.

While there was interest in this model, and the authors did acknowledge many of the limitations of the approach, what was most concerning was that many of the parameters chosen and conclusions seem to conflict with prior experimental evidence. The concern is that this work might cement beliefs that have evidence to the contrary, without reconciling these differences. Some of the examples below are from several years ago, and perhaps could be replicated. Would your results motivate a re-look at those older experiments? Were the conditions different? Or do those old results disprove your predictions, uncovering critical flaws in the simulation? That really is the main crux of the concerns here.

1. “We simulated SD in both perfused and hypoxic slices. Our model confirms that SD speed is augmented by propionate and hypoxia and suggests that changing the ECS is the principle mechanism through which they influence SD.” And “Hypoxia will increase SD speed under multiple conditions”.

This contrasts with the true findings in hypoxia. Aitken et al. 1998 wrote the following: “Both events (normoxic and hypoxic SD) arose focally, after which a clear advancing wave front of increased reflectance and DC shift spread along the CA1 stratum radiatum and s. oriens. The rate of spread was similar: conduction velocity of normoxic SD was 8.73 +/- 0.92 mm/min (mean +/- SE) measured electrically and 5.84 +/- 0.63 mm/min measured optically, whereas (hypoxic) HSD showed values of 7.22 +/- 1.60 mm/min (electrical) and 6.79 +/- 0.42 mm/min (optical). When initiated in CA1, normoxic SD consistently failed to enter the CA3 region (7/7 slices) and could not be initiated by direct KCl injection in the CA3 region (n = 3). Likewise, the hypoxic SD-like optical signal showed onset in the CA1 region and halted at the CA1/CA3 boundary (9/9 slices)...”

"Once the K+ wave arrived however, the cells dumped large quantities of K+ into the ECS due to the burst and subsequent prolonged depolarization (note rapid depletion of [K+]i - Fig. 4B).”

In other words, Aitken et al. 1998 found that the rate of spread was similar during normoxic and hypoxic SD. We agree that it appears highly plausible (without experiments) that SD speed could be higher during hypoxia-and may have believed this to be true as well. However, this was not the case in real measurements. Maybe, these measurements should be reproduced, but the Somjen group which did these measurements is trustworthy. Many of their findings could be reproduced by others.

In all this, it is important to realize that the extracellular space (ECS) shrinks by about 70% during SD. Normally, the ECS makes ∼20% of the total volume. During SD, the ECS makes ∼6% of the total volume. Even if not a single potassium ion left the intracellular space, the extracellular potassium concentration would increase from ∼3 to ∼10 mM. I.e., not many potassium ions have to leave the cells to increase the extracellular potassium concentration from ∼3 to ∼50 mM which is measured during SD. The truth is that only a small quantity of K+ is dumped into the ECS during SD in contrast to the large quantity of Na+ dumped into the intraneuronal space. Plausibility considerations based on volume ratios and volume changes suggest that more than 10 times as many sodium ions are taken up into the cells during SD than potassium ions leave the cells. Correspondingly, very large numbers of chloride ions are also taken up into the cell, and this is measured accordingly (Hansen and Zeuthen, 1981; Kraig and Nicholson, 1978; Windmuller et al., 2005). These relationships are important and should be better explained. The central ion to understand SD may not be potassium but sodium. At least the sodium shifts are much larger than the potassium shifts in absolute terms. Also, during recovery, while it is not wrong that the sodium pump also restores extracellular potassium to normal, it is primarily activated by intracellular sodium and must remove much more sodium from the intracellular space than potassium from the ECS. Sodium should not be forgotten in all this. The driving force leading to SD is related to the sodium gradient rather than the potassium gradient (Lemale et al., 2022).

2. “The Na+/K+ pump activity which contributed to K+ elimination from the ECS in the core of the slice created a high demand for O2, 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+ wave (compare rapid fall-off of O2 in Fig. 4E to much slower rise of extracellular K+ in Fig. 4F). The rapid spread of O2 deficit explains the total pump failure at intermediate locations in the slice.”

The authors seem to forget that the tissue has an ATP reserve. In vivo and also in vitro, complete ischemia/hypoxia usually requires a minute or even longer before SD occurs. Even severe hypoxia does not imply an immediate lack of ATP. Moreover, we are unaware of evidence that tissue oxygen falls before the SD front in normal tissue. It typically falls behind the wave front. For example, compare figure 2 in (Piilgaard and Lauritzen, 2009) or the figures in (Takano et al., 2007). The measured findings seem to directly contradict your findings of tissue hypoxia preceding the wave front. In addition, hyperbaric oxygenation does not seem to have profound effects on SD (e.g. figure 5 in (Lindauer et al., 2010)). In real measurements, tissue oxygen seems to drop behind the wave front, if it drops at all. Under certain conditions, tissue oxygen can also increase in response to SD (examples are given in figure 2 (Balanca et al., 2017) and figure 1 (Major et al., 2021).

3 “Our results support Grafstein’s hypothesis that SD is primarily a reaction-diffusion process mediated by a regenerative K+ wave [Grafstein, 1956].”

Somjen wrote in his review in 2001 with regard to a hypothesis based on a possible propagation mechanism based on potassium (Somjen, 2001): “Yet another problem is that, at a given point in the tissue, [K+]o does not start to increase ahead of the ΔVo, as it should, if K+ were the agent of the propagation of SD”. Thus, the original recordings contradict this hypothesis and they have been independently reproduced by several groups. The voltage shift is only preceded by an increase in [K+]o when K+ is either applied artificially or hypoxia is induced which activates ATP-sensitive and G protein-dependent Ca2+-sensitive K+ channels long before it induces SD (Erdemli et al., 1998; Muller and Somjen, 2000; Revah et al., 2016). There are good arguments that a reaction-diffusion process is the basis of SD, but there are also good arguments against K+ as the neuroactive factor responsible for SD propagation. The most important one is the one formulated by Somjen above. Original recordings should not be neglected. Maybe, this must formulated more cautiously. It is likely that the neuroactive factor responsible for the reaction diffusion process is of a different nature than potassium.

4. In Figure 4C, it appears too extreme that the intracellular sodium concentration can increase to values of up to 80 mM during SD. The intracellular space is four times larger than the ECS. Even if you remove astrocytes, it is impossible that intracellular sodium increases to 80 mM. Plausibility suggests that it increases to 36 mM (Lemale et al., 2022) and true measurements suggest increases from 10 to between 20 and 30 mM (Gerkau et al., 2017). So, this part of the figure is not convincing and suggests that you did not sufficiently consider the volume ratio and volume changes between intracellular space and ECS. In this respect, you may also have a look at Figure 4 in (Somjen, 2001).

Minor comments

The following two statements are synonyms: “SD can be inhibited by enlarging ECS volume” and “SD velocity will be greater in areas with greater neuronal density, total neuronal volume, or
larger/more dendrites”. Greater neuronal density result in low ECS volume, which would increase the propagation of SD as indicated by the first statement.

Figure 4 legend: Remove one “shows” in the text “Supporting Movie 4-1 shows shows extracellular ion”.

L 197: Should “S:V = 0.02 μm” be “S:V = 0.02/μm”?

L 198: “K+ wave speed also increased with increased βnrn while keeping S:V and cell density constant (Fig. 6B).” It is not clear how can the cell density remain constant when neuronal volume relative to the tissue volume is changing? This comment also applies to the statements that follow: “Cell density effects were less marked, whether through keeping surface area and volume constant (Fig. 6C), or keeping βnrn and S:V constant (Fig. 6D).”

Legend of figure 4: “Supporting Movie 4-1 shows shows extracellular ion and O 2 concentrations...”

"SD is difficult to detect in humans noninvasively [Zandt et al., 2015], ...”

The original papers should also be cited here: (Drenckhahn et al., 2012; Hartings et al., 2014; Hofmeijer et al., 2018).

"SD has been studied in brain slices from a wide range of species and brain regions, including neocortex, hippocampus, brainstem, and retina... It can be triggered experimentally by various means, including electrical stimulation; mechanical insult; K+, glutamate, and ouabain application [Cozzolino et al., 2018].”

There is an important error in this sentence, which is a common misunderstanding and should be eradicated. Bath exposure to glutamate (up to 1 mM, which is pathophysiological) or glutamate receptor agonists induces depolarization and a negative direct current (DC) shift wherever one chooses to record in a brain slice. Imaging light transmission (LT) changes in an entire slice shows that this depolarization arises in the majority of imaged gray matter almost simultaneously, i.e. the depolarization does not spread throughout neocortical and hippocampal gray matter. A range of concentrations of bath-applied glutamate (0.01 to 1.0 mM), NMDA, or kainate fails to induce a “spreading” depolarization (Andrew et al., 1996; Andrew et al., 1999; Obeidat and Andrew, 1998). So, the elevated LT signal evoked by glutamate or by NMDA (denoting cell swelling and depolarization (Rungta et al., 2015)) develops uniformly in the slice and does not propagate. This contrasts with bath exposure to high K+, ouabain, hypoxia or oxygen-glucose deprivation which produce a true “spreading” depolarization (Aitken et al., 1998; Joshi and Andrew, 2001; Major et al., 2017; Petzold et al., 2005). Thus, glutamate should be removed from the list above.

"This mirrors findings in vivo where SD-associated metabolic demand for O2 can exceed supply, resulting in tissue hypoxia [Takano et al., 2007].”

The following references could be added for rats and humans (Piilgaard and Lauritzen, 2009; Winkler et al., 2017), but again your traces do not mirror the traces in real experiments where the drop in tissue oxygen follows the SD wave front.

The figure legends should have better descriptions of what is happening in the legends-they are too short. Non-modelling people should also be able to easily understand what is shown in the figures. Help the reader with brief summaries of what is shown. If more people understand what the authors want to explain in the figure, more people will cite this paper.

References cited above:

Aitken PG, Tombaugh GC, Turner DA, Somjen GG. Similar propagation of SD and hypoxic SD-like depolarization in rat hippocampus recorded optically and electrically. J Neurophysiol 1998; 80: 1514-21.

Akerman S, Holland PR, Goadsby PJ. Mechanically-induced cortical spreading depression associated regional cerebral blood flow changes are blocked by Na+ ion channel blockade. Brain Res 2008; 1229: 27-36.

Andrew RD, Adams JR, Polischuk TM. Imaging NMDA- and kainate-induced intrinsic optical signals from the hippocampal slice. J Neurophysiol 1996; 76: 2707-17.

Andrew RD, Jarvis CR, Obeidat AS. Potential sources of intrinsic optical signals imaged in live brain slices. Methods 1999; 18: 185-96, 79.

Balanca B, Meiller A, Bezin L, Dreier JP, Marinesco S, Lieutaud T. Altered hypermetabolic response to cortical spreading depolarizations after traumatic brain injury in rats. J Cereb Blood Flow Metab 2017; 37: 1670-86.

Drenckhahn C, Winkler MK, Major S, Scheel M, Kang EJ, Pinczolits A, et al. Correlates of spreading depolarization in human scalp electroencephalography. Brain 2012; 135: 853-68.

Erdemli G, Xu YZ, Krnjevic K. Potassium conductance causing hyperpolarization of CA1 hippocampal neurons during hypoxia. J Neurophysiol 1998; 80: 2378-90.

Gerkau NJ, Rakers C, Petzold GC, Rose CR. Differential effects of energy deprivation on intracellular sodium homeostasis in neurons and astrocytes. J Neurosci Res 2017; 95: 2275-85.

Hansen AJ, Zeuthen T. Extracellular ion concentrations during spreading depression and ischemia in the rat brain cortex. Acta Physiol Scand 1981; 113: 437-45.

Hartings JA, Wilson JA, Hinzman JM, Pollandt S, Dreier JP, DiNapoli V, et al. Spreading depression in continuous electroencephalography of brain trauma. Ann Neurol 2014; 76: 681-94.

Hofmeijer J, van Kaam CR, van de Werff B, Vermeer SE, Tjepkema-Cloostermans MC, van Putten M. Detecting Cortical Spreading Depolarization with Full Band Scalp Electroencephalography: An Illusion? Front Neurol 2018; 9: 17.

Joshi I, Andrew RD. Imaging anoxic depolarization during ischemia-like conditions in the mouse hemi-brain slice. J Neurophysiol 2001; 85: 414-24.

Kraig RP, Nicholson C. Extracellular ionic variations during spreading depression. Neuroscience 1978; 3: 1045-59.

Lemale CL, Luckl J, Horst V, Reiffurth C, Major S, Hecht N, et al. Migraine Aura, Transient Ischemic Attacks, Stroke, and Dying of the Brain Share the Same Key Pathophysiological Process in Neurons Driven by Gibbs-Donnan Forces, Namely Spreading Depolarization. Frontiers in cellular neuroscience 2022; 16: 837650.

Luckl J, Lemale CL, Kola V, Horst V, Khojasteh U, Oliveira-Ferreira AI, et al. The negative ultraslow potential, electrophysiological correlate of infarction in the human cortex. Brain 2018; 141: 1734-52.

Major S, Petzold GC, Reiffurth C, Windmuller O, Foddis M, Lindauer U, et al. A role of the sodium pump in spreading ischemia in rats. J Cereb Blood Flow Metab 2017; 37: 1687-705.

Muller M, Somjen GG. Na(+) and K(+) concentrations, extra- and intracellular voltages, and the effect of TTX in hypoxic rat hippocampal slices. J Neurophysiol 2000; 83: 735-45.

Obeidat AS, Andrew RD. Spreading depression determines acute cellular damage in the hippocampal slice during oxygen/glucose deprivation. Eur J Neurosci 1998; 10: 3451-61.

Petzold GC, Windmuller O, Haack S, Major S, Buchheim K, Megow D, et al. Increased extracellular K+ concentration reduces the efficacy of N-methyl-D-aspartate receptor antagonists to block spreading depression-like depolarizations and spreading ischemia. Stroke 2005; 36: 1270-7.

Piilgaard H, Lauritzen M. Persistent increase in oxygen consumption and impaired neurovascular coupling after spreading depression in rat neocortex. J Cereb Blood Flow Metab 2009; 29: 1517-27.

Revah O, Lasser-Katz E, Fleidervish IA, Gutnick MJ. The earliest neuronal responses to hypoxia in the neocortical circuit are glutamate-dependent. Neurobiol Dis 2016; 95: 158-67.