The Role of Cell Volume in the Dynamics of Seizure, Spreading Depression, and Anoxic Depolarization

Cell volume changes are ubiquitous in normal and pathological activity of the brain. Nevertheless, we know little of how cell volume affects neuronal dynamics. We here performed the first detailed study of the effects of cell volume on neuronal dynamics. By incorporating cell swelling together with dynamic ion concentrations and oxygen supply into Hodgkin-Huxley type spiking dynamics, we demonstrate the spontaneous transition between epileptic seizure and spreading depression states as the cell swells and contracts in response to changes in osmotic pressure. Our use of volume as an order parameter further revealed a dynamical definition for the experimentally described physiological ceiling that separates seizure from spreading depression, as well as predicted a second ceiling that demarcates spreading depression from anoxic depolarization. Our model highlights the neuroprotective role of glial K buffering against seizures and spreading depression, and provides novel insights into anoxic depolarization and the relevant cell swelling during ischemia. We argue that the dynamics of seizures, spreading depression, and anoxic depolarization lie along a continuum of the repertoire of the neuron membrane that can be understood only when the dynamic ion concentrations, oxygen homeostasis,and cell swelling in response to osmotic pressure are taken into consideration. Our results demonstrate the feasibility of a unified framework for a wide range of neuronal behaviors that may be of substantial importance in the understanding of and potentially developing universal intervention strategies for these pathological states.


Introduction
Cells swell during a wide variety of pathologies, including trauma, ischemia, hypoxia, seizures, and spreading depression [1][2][3]. Changes in osmolality can change the susceptibility to epileptiform activity [4][5][6], and affect the amplitude of intra-and extracellularly recorded electrical signals [7]. Cells also change their volume during normal activity, and the change in cell size during individual action potentials has been estimated [8,9]. Despite this ubiquity of observed phenomena, the effect of cell swelling on single cell behavior is incompletely understood.
It is now accepted that the dynamic microenvironment within the extracellular space (ECS), modified by ionic fluxes from neurons, glia, and blood vessels, plays a critical role in neuronal behavior [1]. In particular, pathological states involving excessive neuronal depolarization such as epileptic seizure (SZ), spreading depression (SD), and anoxic depolarization (AD) during ischemia are characterized by major rearrangements of various ions across the cell membrane and neuronal microenvironment [1,[10][11][12][13][14][15][16]. In each of these three conditions, collapse of transmembrane ionic gradients requires enhanced oxygen and glucose consumption required by active transport systems to reestablish the gradients [17,18]. For the purpose of this paper, we define SZ, SD, and AD respectively as the ion concentrations-induced high-frequency bursts not usually seen in the normal condition of the same cell [1,19], the nearly complete depolarization of the cell's membrane potential that recovers spontaneously on the scale of seconds [13,19], and the nearly complete depolarization of the cell's membrane potential triggered by oxygen (O 2 ) and glucose deprivation (OGD) that may or may not recover depending on the cell type after O 2 and glucose is restored [1,20,21].
During the pathological states mentioned above, the massive rearrangement of ions across plasma membrane drives water molecules from the extra-to intracellular space leading cell swelling. For example, pyramidal neurons in slices from cortical layer V swell by as much as 60% during AD caused by 20 minutes OGD [20]. Although lacking functional aquaporins, neurons swell significantly in response to OGD and extracellular K + elevations [22]. Although still debated, the K + /Cl − and Na + /K + /2Cl − cotransporters are suspected to mediate the entry of water molecules into neurons [23,24]. Astrocytes on the other hand, express aquaporins [25]. The clearance of excessive K + due to high neuronal activity by astrocytes leads to osmotic gradients resulting in water influx through aquaporins and astrocytic dilation [26][27][28]. This paper focusses on the role of neuronal swelling in response to osmolality changes in cell behavior without considering the specifics of pathways involved in the water influx.
Despite the fact that SD was first observed by Leão as the silencing of spontaneous electrical activity during experiments on epileptic SZs, the two phenomena have long been considered separate physiological events [1]. They are characterized by different patterns of neuronal activities [29][30][31], characteristic ionic changes [13][14][15][16]19], and known interactions with oxygen [17,32]. By expanding the Hodgkin-Huxley type framework to incorporate conservation of particles and charge, and accounting for the energy required to restore ionic gradients, we recently uncovered a unified mechanism for SZ and SD [33]. Specifically, we showed that a wide-range of neuronal behaviors can be accounted for as a function of the cell's extracellular potassium concentration and oxygen supply. More recently, Hübel and Dahlem performed a detailed bifurcation analysis of different time-scales arising from the consideration of dynamic ion concentrations in conjunction with Hodgkin-Huxley type framework [34].
Extensive work has been done on the role of ion pumps, channels, and transporters in stroke [35]. Recently, Andrew and colleagues showed that neuronal populations in lower brain regions such as the hypothalamus are resistant while those in higher brain regions such as neocortex are more susceptible to ischemic injury [20]. They further showed that the thalamushypothalamus interface represents a discrete boundary where neurons in thalamus are more vulnerable than hypothalamus to ischemia, generating stronger AD in response to OGD, and do not recover as readily after restoring normal O 2 and glucose supply [21]. The authors postulated that the variability of ATP-dependent Na + − K + pumps in these regions could lead to the contrasting neuronal response in OGD conditions.
In this paper, we explore the effect of cell swelling on neuronal behavior by demonstrating the ability of cell volume to act as a bifurcation parameter. We seek a better understanding of how human brain cells respond to osmotic pressure-induced swelling in states such as SZ, SD, and AD, and universal intervention strategies for controlling these conditions. Here, we show that spontaneous transition between SZ and SD can be seen in a model neuron if volumetric changes in response to intense neuronal activity and ionic fluxes are taken into account. Without any adjustments, our model behaves in similar fashion as in vitro experiments under OGD. We further show that the variability in the geometry and microenvironment of neurons could play a significant part in their differential response in OGD conditions observed in in vitro experiments in different brain regions. Based on our results, we conclude that combining ion concentration dynamics during spiking with the sizes of intra-and extracellular spaces supports a unified framework for epileptic SZ, SD, and AD.

Transition between SZ and SD states
We investigate the role of cell size and relative (to intracellular volume) extracellular space on neuronal behavior by varying the radius of the cell, r in , keeping the total radius of cell and extracellular space, r tot , fixed. Thus changing r in is largely equivalent to changing the ratio of intra-to extracellular volume. Since we are interested in the pathological states of the cell, we use K + concentration in the distant reservoir [K] o,1 = 8mM-a value typically used for inducing SZ in in vitro [30]. The behavior of the cell changes dramatically as we increase r in . A bifurcation diagram showing the maximum and minimum of [K] o as a function of r in is shown in Fig 1A. Briefly, we fixed r in and ran the simulation generating a time-trace representing [K] o versus time. Depending on the r in value, [K] o either oscillates or converge to a steady state value. The initial few hundred seconds of the time-trace were discarded as a transient period and the maximum and minimum of [K] o values in the remaining trace were recorded. The lower and upper markers (green) respectively at a given r in in Fig 1A represent  For all r in < 4.66μm, [K] o remains unchanged and the cell remains in steady state with V % − 60mV ( Fig 1B). As we increase r in above 4.66μm, [K] o enters a periodic orbit via a Hopf bifurcation and starts oscillating with a small amplitude and the cell exhibits spontaneous periodic SZs ( Fig 1C) similar to those observed in experiments [36]. There is a sudden increase in the amplitude of [K] o oscillations at r in = 4.826μm where the peak [K] o goes well over 26mM-a concentration that is often used for inducing SD [37]. The periodic SZs transform to a behavior where the cell is locked into a depolarized state after burst-spiking and exhibits a few smallamplitude spikes on its way out of the depolarized state ( Fig 1D). As we increase r in further, this state disappears making way for mixed SZ and SD behavior where the high-frequency spiking is followed by the locking of V into a depolarized state and the cell comes out of the depolarized state without spiking (Fig 1E). Such mixed states are typically seen in the cells in hypoxic SD [32] or immature physiological conditions [38]. It is worth mentioning that this locking of neuronal membrane into depolarized state is the condition for SD at the single cell level [39]. At the network or tissue level the depolarization may also propagate [40]. This unification of SZ and SD dynamics is supported by the increasing discovery of monogenic mutations in humans that lead to both SZs and migraines [41]. The cell exhibits SZ-SD mixed behavior until it makes a transition to a silent state via another Hopf bifurcation at r in > 4.924μm where V remains fixed at a stable resting value. Cell shows a variety of behaviors as we vary its size. Here we consider Vol (Eq 14) as a bifurcation parameter and simulate Eqs (1,3,5,7,10). (A) Bifurcation diagram of [K] o as a function of r in where green circles, red solid line, and blue dashed line represent periodic orbit, stable, and unstable steady states respectively. The four regions in the bifurcation diagram marked by SZ, SD, and RS represent the parameter-regions where seizure, spreading depression, and resting states are observed respectively. The black solid and dashed limit cycles represent the change in r in as different ion concentrations vary during a single SZ at [K] o,1 = 8mM and 9mM respectively. For the limit cycles, the instantaneous r in values in the limit cycles are obtained from Eq (13). The four panels on the right show membrane potential (black), reversal potential of K + (blue), Cl − (green), and Na + (red) of the cell with different r in values given in the right corner of each panel.
The bifurcation diagrams for [Na] i , [K] i , and [Cl] i (similar to Fig 1A) (not shown) show that their behavior contrasts with [K] o . That is, the amplitude (defined as the difference between the maximum and minimum) of [Na] i , [K] i , and [Cl] i oscillations decreases as we increase r in . At r in > 4.924μm, the relatively larger intracellular volume dominates the extracellular space and these 3 concentrations drop to resting values leading to a return to resting membrane potential. That is, the extremely large intracellular volume leads to low intracellular concentrations that overshadow the effect of the changes in extracellular ion concentrations. The bifurcations in the cell's behavior described above are qualitatively preserved when [K] i , [Na] o , and [Cl] o are modeled by rate equations [33] instead of conservation equations (S1 Fig). It is important to point out that the kind of changes in r in shown in Fig 1A and the rest of the paper are physiologically relevant. For example, the surface area of cortical layer V pyramidal neurons increases by more than 50% in response to OGD [20]. Thus a spherical cell with initial radius of 4.75μm would swell to a final radius of 5.81μm, shrinking the extracellular space significantly.
The change in the ratio of intra-to extracellular volume as a result of changing r in plays a major role in the transition between SZ and SD states. Depending on the value of β, the cell exhibits steady state, SZ, or SD without any transition between these behaviors as we change The results in Fig 1 indicate that the size of the cells and how tightly they are packed in the tissue can play a significant role in their dynamics. An important followup question would be: could a cell swell enough so that it would spontaneously go through the transitions shown in Fig 1? To answer this question, we add the volume dynamics given by Eq (14) to our model, where the cell volume depends on the instantaneous ion concentrations. We compute the spontaneous change in r in as the ion concentrations inside and outside the cell vary during a single SZ (solid black line in Fig 1A). The limit cycle shows that r in can change enough during one SZ so that the cell would make the transition from SZ to SD regions. A cell with r in = 4.82μm (which is in the SZ region and will exhibit spontaneous SZs similar to Fig 1C) before a SZ starts would swell to a final r in = 4.84μm (solid black line in Fig 1A), well within the SD region ( Fig 1D). The crossover to the SD region is more prominent for higher K + in the bath solution (dashed black line in Fig 1A). Fig 2 shows the pathway to the spontaneous transition from SZ to SD caused by cell swelling obtained by simulating the full model (Eqs 1-15). The arrows in Figs 2, 6, and 9 indicate the direction of the trajectory. Initially [Na] i and [K] o slowly build up leading to an increase in intracellular volume (hence a decrease in the relative extracellular space) and higher excitability of the cell. The microenvironment reaches a point where V makes a transition from steady state to limit cycle via a saddle node on invariant cycle (SNIC) bifurcation and the cell exhibits a SZ (Fig 2A). After exiting the SZ state, the cell does not have sufficient time to reverse the swelling (normally pumps would restore ionic gradients reversing the swelling) and the cell exhibits a SD event by entering a second periodic orbit via a SNIC, with progressively decreasing amplitude followed by entrance into a depolarized state via a Hopf bifurcation. What physiological mechanisms help to regulate the brain so that most of the time, even in people with chronic recurring seizures and migraines, their brains are operating normally? To begin to address this question we show a two-parameters bifurcation diagram for [K] o in Fig 3 where the two parameters are r in and glial K + buffering strength, B glia . In Fig 3A, we show the maximum of [K] o as a function of r in and B glia . As is clear, for weaker glial buffering the cell's behavior changes as a function of r in in the same manner as in Fig 1. However, for stronger glial buffering (B glia > 21mM/sec) the cell becomes biased towards SZ behavior and is less likely to go to SD. The stronger glial buffering siphons away [K] o fast enough so that the cell recovers from swelling during SZ before entering SZ again. For even stronger glial buffering, the cell neither shows SZ nor SD behavior. This point is emphasized further in Fig 3B where we show the three regions described in Fig 1 as a function of r in and B glia . For larger B glia , both SZ and SD regions disappear. This demonstrates that glial K + buffering can play a neuroprotective role against SZ and SD behaviors and, if strong enough, will constrain the neuron to physiological dynamics. From our previous work, we know that extracellular diffusion shares similar dynamical properties with the glial buffer (see, e.g., Fig 7 in [18]).
A two parameter bifurcation diagram for [O 2 ] from simulations in Fig 3 shows that the resting state at large r in (right side in Fig 3B) is energetically favorable as compared to the one at small r in (left side in Fig 3B) (Fig 4). Interestingly, the local O 2 consumption during SD in a cell with smaller radius is much larger than the consumption in a cell that has swollen. This indicates that although the [K] o changes during SD in smaller cells are smaller, Na + − K + exchange pumps consume more energy to restore physiological [K] o and [Na] i in these cells due to their relatively larger extracellular reservoirs. A cell with an infinitesimally small extracellular space requires an infinitesimally small ionic flux to substantially change its extracellular potential, and an infinitesimally small work load on the membrane pumps to restore those gradients. As mentioned above, the bifurcation diagrams for [Na] i , [K] i , and [Na] o (similar to Fig 1A) (3)), higher concentrations of [Na] i in a smaller cell overwhelms the relatively smaller concentrations of [K] o and the Na + − K + exchange pumps increase their rates accordingly to rearrange those ions (Fig 4). It is important to point out that the amount of oxygen available to the cell also acts as a bifurcation parameter causing the cell to transition between different states [18,33].

Bifurcation analysis of different neuronal behaviors
To gain further insights into the mechanisms behind the three behaviors in Fig 1C-  At r in = 4.84μm (Fig 1D), there is a substantial change in [Na] i and [K] o and the one-parameter bifurcation diagram goes through the transitions shown in Fig 5A-5D. The beginning of the first burst starts when trajectory crosses the SNIC curve in Fig 6 and terminates when it crosses the HB curve. Again, the frequency within the burst is low at the beginning. The amplitude shrinks to zero at the end as the trajectory passes through the HB curve. After the first burst terminates, the trajectory follows the stable branch of SSS2 ( Fig 5) for a while before it crosses the HB curve again. The solution does not start to burst at this moment, because of the delayed loss of stability phenomenon (delayed HB) [43] that occurs when a parameter passes slowly through a HB point and the system's response changes from a slowly varying steady state to a slowly varying oscillation. So, the trajectory traces the unstable branch of SSS2 for a while before it starts to burst between HB and HC curves. Also referred to as "delayed" or "memory effect", this transition happens when the parameter is considerably beyond the value predicted from a straightforward bifurcation analysis which neglects the dynamic aspect of the parameter variation. This memory effect has been studied for different problems including the FitzHugh-Nagumo model [44]. The second burst terminates when the trajectory crosses the HC curve (notice the slight delay in the termination of the second burst for higher [Cl] i values). At this moment, the solution drops to the lower stable branch SSS1. Since [Na] i is approximately 30mM for the second burst, the amplitude of the second burst is smaller according to Fig 5. The delay also explains why the burst onset has a nonzero amplitude since the solution is away from the HB curve. The amount of delay depends on the time spent near the attracting branch of SSS2 (the loop from HB to HB).
The explanation for r in = 4.84μm also applies to the case with r in = 4.9μm (Fig 1E), except that the HB to HB trajectory loop is more dramatic than r in = 4.84μm. The trajectory spends more time on the attracting branch of SSS2 that increases the delay. Here, the solution moves along the unstable branch through the entire HB-HC regime without oscillating. It finally loses stability, but since there is no stable limit cycle anymore, it can only drop to the lower stable branch of SSS1. SSS2 (the depolarized state of V) ends to the left of HC curve. The drop to SSS1 (end of depolarized state) is faster (as in Fig 1D) for a slightly smaller cell (for example when r in = 4.88μm) because the trajectory is close to the SNIC-HC bifurcation point and hence the upper and lower branches are closer (not shown).

Anoxic depolarization
In addition to SZ, SD, and SZ-SD transition, our model closely reproduces AD. AD, the initial electrophysiological event during ischemia, is a front of depolarization that drains residual stored energy in compromised gray matter. Recently, Brisson and Andrew studied AD at the single cell level in neocortex and hypothalamus slices deprived of O 2 and glucose [20]. They used two-photon microscopy to image cell swelling simultaneously with patch-clamp membrane potential measurements in the OGD condition. Although, our model does not have the glucose component, it behaves in the same manner as that observed experimentally during energy substrate deprivation. We induced energy deprivation (ED) in the model by putting oxygen in the perfusion solution, [O 2 ] 1 , equal to zero, which is model-equivalent to OGD or ischemia-induced AD [45]. In Fig 7A we show the cell behavior in response to 5 min ED. During AD, the cell swells qualitatively in the same manner as observed in experiments (Fig 7B and  7E). The extra-and intracellular ion concentrations go through a massive redistribution during The Role of Volume in Neuronal Dynamics AD ( Fig 7C). Consistent with the experimental observations (shown in panels D-F), the cell depolarizes and swells significantly during AD. The behavior of the model-cell is reminiscent of the magnocellular neuroendocrine cell (MNC) in the hypothalamus with weak AD that depolarizes only transiently close to 0mV (Fig 7D). However, as we change the initial radius of the cell from 4.0μm to 3.0μm, it exhibits strong AD, approaching a steady Donnan equilibrium of − 2mV (red line in Fig 7A). This behavior is reminiscent of strong AD in pyramidal neurons from neucortical layer V where the cell's membrane depolarizes to near 0mV during OGD and does not recover to the pre-OGD state when O 2 and glucose supply is restored (Fig 7F). Our result sheds light on the importance of cell size and packing density in the strength of AD and the potential damage to the cell where the membrane potential is permanently brought into a depolarized state with no recovery to the pre-ED state.
During the modeling of AD, we made some additional observations. The main conclusion about the transition between SZ and SD states qualitatively remain the same for fixed Cl − concentrations. However, the cell cannot sustain a prolonged depolarization similar to AD for fixed Cl − concentrations (S3 (A) Fig), confirming an important role of Cl − in AD [46]. A similar behavior is observed when we maintain the normal K + diffusion between extracellular  [47], we needed to make diffusion a function of O 2 [33] so that there is decreased exchange of K + between bath or blood vessels and extracellular space when the cell is exposed to reduced oxygen. In addition, we had to make the K + glial buffering a function of O 2 to reproduce AD (as in [48]). Several observations support these assumptions. A significant portion of [K] o is buffered by astrocytes through ATP-dependent Na + − K + pumps that do not function in the absence O 2 and glucose. In the absence of O 2 , astrocytes attempt to buffer the increased extracellular K + by switching to anaerobic glycolysis and swell substantially [49], further restricting K + diffusion and limiting glial energy reserves. Astrocytic inward rectifying K + channels (K ir ) also contribute to K + siphoning, gating through interaction with G-protein coupled receptors pathway that is dependent on ATP (see [50] for a review on K ir channels). Similarly, Na + − K + − Cl − cotransporters (NKCC) that are found in astrocytes play a significant role in transferring K + (together with Na + and Cl − ) from extracellular space to astrocytes and are dependent on ion gradients [51] and thus indirectly on ATP. Hence ATP plays a crucial role in these pathways that would be disrupted in the absence of O 2 , leading to reduced K + buffering.
The other important effect is the reduced transport of ions at glial end-feet adjoining the pericapillary space. The two-way K + trafficking at the blood-brain barrier occurs at the junctions between astrocytic end-feet and blood vessels (see for example [52]). Astrocytes release K + next to tight junction sealed endothelial cells in blood vessels. Na + − K + pumps transfer that K + to the endothelial cells and it is then delivered into the blood vessels through K + channels. A reverse process transfers K + from blood vessels to astrocytes and finally to neurons ( Fig   Fig 8. K + exchange at the blood-brain barrier. Extracellular K + absorbed by astrocytes is transferred to the blood vessels through a sequential functioning of K + channels and Na + − K + exchange pumps at the junction between astrocytes and endothelial cells surrounding the blood vessel lumen. A reverse process transfers K + from blood vessels to astrocytes. Omitted from this picture is the K + exchange between astrocytes through gap junctions. doi:10.1371/journal.pcbi.1004414.g008 The Role of Volume in Neuronal Dynamics 8). The lack of ATP in AD would disrupt this pathway, consequently impairing the K + diffusion between blood vessels and extracellular space.
In Fig 9, we show the pathways of microenvironment changes leading to AD. Immediately after setting O 2 to zero, there is a rapid increase in [K] o and [Na] i followed by a slow increase in intracellular volume (Fig 9A).
[K] o drops slightly from its peak value once it enters the AD due to the slight delay in the [Cl] i rise (see Fig 7C). After the initial drop, [K] o stabilizes at fixed value when [Cl] i plateaus at its peak value. Once O 2 is restored, [K] o is restored to normal values followed by slow restoration of [Na] i and intracellular volume. The blue curve is replotted from Fig (2B)

Discussion
This is, to our knowledge, the first detailed study of the effects of cell volume on neuronal dynamics. We employed a recently discovered unification framework, where extending the Hodgkin-Huxley equations for mammalian neurons with energy balance and conservative principles demonstrated that a broad variety of neuronal states lie along a continuum of the repertoire of the neuronal membrane [33]. Using a variety of model simplifications, we were able to perform detailed bifurcation analyses that explained the full model effects as a function  (2) and (7) respectively. %ΔVol represents the percent change in volume and is defined as ððr 3 in;ins À r 3 in;ss Þ=r 3 in;ss Þ Â 100, where r in,ss and r in,ins represent initial steady state (ss) and instantaneous (ins) radius of the cell respectively. of volume itself as an order parameter. We can now better understand and unify a range of effects and factors that are critical in the transitions to the pathological states of SZ, SD, and AD.
Our study of the role of volume as an order parameter revealed a dynamical definition for the experimentally described physiological ceiling that separates seizure from SD activity [42]. Furthermore, we have delineated and predict a second ceiling, one that demarcates SD from AD. Our observations unveil a new way of investigating neuronal behavior where different states need not be treated separately but rather as a dynamical continuum of the neuronal membrane potential, ion concentrations, metabolic energy, and volume.
Previous experiments support our observations about the volume as an order parameter regulating the neuronal excitability. Exposure of brain cells to hyposmolality and the resultant shrinkage of extracellular space at clinically relevant levels promotes epilepiform activity in hippocampus and neocortex [53][54][55] and clinically [5]. Hyperosmolality (dehydration), on the other hand inhibits epilepiform activity. Mannitol, which reverses the hyposmolar state, abolishes synchronous neuronal bursting in the dentate gyrus and CA1 area of the hippocampus [4]. Furthermore, hyposmolality increases the amplitude of evoked field potentials and excitatory postsynaptic potentials recorded intracellularly in rat neocortical slices [7]. Conversely, mannitol-induced hyperosmolality reverses these features [7]. Similarly, decreasing the extracellular osmotic pressure converted non-bursting neurons to bursting neurons and decreased the stimulus requirements for evoking burst firing in native bursters, while increasing extracellular osmotic pressure suppressed burst firing [56].
Glial reactivity and scarring is a prominent feature of a broad variety of brain pathologies [57], and is prominent in epilepsy [58]. We found that when glial K + buffering is impaired, the cells can swell enough to cause transition from SZ to SD. Previous observations by Foley et al. [59] support this mechanism. The excitability of neurohypophysial neurons due to accumulating K + in the extracellular space decreased significantly by increasing the size of interstitial space [59]. The decrease in action potential amplitude (showing the cell's transition towards a depolarized state (SSS2)) in control cells in response to a train of 40 stimuli at 25Hz, diminished when a hypertonic solution containing 100mM sucrose was added to the normal ringer solution perfusate [59]. The membrane-impermeable sucrose increases extracellular osmotic pressure causing neurons to shrink. A similar protocol could be used to test the predictions of our model. Our hypothesis is that adding hypertonic solution to the tissue going though SD would transition the cell to SZ and possibly to the resting state. Taking glial swelling into account would make the shrinkage of interstitial space more dramatic, further supporting our claim that the reduction in extracellular space can be strong enough to cause the cell to transition from SZ to SD. Our results also highlight the importance of glial K + buffering strength in pathological states and shows that the cell is less likely to go into SZ and SD and the transition between the two when glia are functioning efficiently.
The different time-scales in our model provide a parsimonious explanation for the transition between SZ and SD states in Fig 1A. As (14)). Strictly speaking, the ultra-fast mgate (m) if modeled with the rate equation instead of the steady state approximation can be considered as a fifth time-scale. We believe that this could underlie the dramatic increase in the amplitude of [K] o oscillations-reminiscent of canard explosion-which is a feature of a slowfast system such as ours [60]. Testing this claim using geometric singular perturbation theory is beyond the scope of this study and will be the subject of future research.
In addition to the spontaneous transition between SZ and SD, our model closely reproduces AD, the hyperpolarization after O 2 is restored, and the cell swelling in AD. The close resemblance of our results to the experimental data makes our model a good candidate for future studies to understand hypoxia and ischemia and ways to protect against metabolic insults.
As pointed out above, Brisson et al. [20,21] emphasized the importance of variable pump rates in the stronger AD observed in the upper brain regions as compared to the lower brain regions in response to OGD. Our study shows that the cell volume and microenvironment play an important role in the strength of AD (Fig 7A). Our preliminary phase space analysis for the cell's failure in the higher brain regions to recover from AD reveals that the cell volume affects [K] o and [Na] i in such a way that they lead to lower pump activity and consequently the cell's failure to recover from AD after O 2 and glucose is restored. A complete phase space analysis of the recovery failure from AD is beyond the scope of this paper and the subject of our future research.
This current model is a simplified picture of a very complex reality. Although Na + − K + -ATPase consumes 91% of total available O 2 [61], other pathways such as Ca 2+ -ATPase and synaptic communication expend significant amounts of metabolic energy which is not taken into account in the current model. Similarly, the clearance of excessive K + accompanied by Cl − uptake and Na + expulsion by astrocytes in face of high neuronal activity would lead to significant dilation of astrocytes [26][27][28] during the pathological states discussed in this paper, which is missing from our model. Also missing from the model is the dynamic intra-and extracellular Ca 2+ concentrations, which are suspected to play a crucial role in neuronal excitotoxicity (see for example, chapter 4 of [1]). Osmotic pressure also affects the re-depolarizing component of spike after-depolarization and apparent membrane time-constant that are attributed to changes in the persistent Na + current [56]. These factors are important for developing a more comprehensive understanding of the pathological states discussed in this paper, and is the subject of future research.
To conclude, our finding further explored a unified mechanism that accounts for SZ, SD, and AD. The detailed biophysical models that take neuronal microenvironment such as ion concentrations, glia, blood vessels, metabolic energy requirements, and volume homeostasis into account will provide a better understanding of these conditions and may lead to unified therapies for SZ, SD, and possibly ischemic stroke-like injury prevention.

Methods
Our model builds on our previous work [15,18,33,62]. We consider a spherical cell of radius r in placed inside a spherical shell of fixed radius r tot = 5μm. So by changing r in , we are in effect changing the ratio of intra-to extracellular volume. The schematic of the model is shown in Fig 1 of [33].

Hodgkin-Huxley type equations
The membrane potential V of the cell is modeled with the following set of modified Hodgkin-Huxley type equations [15,63] C dV dt ¼ I Na þ I K þ I L À I pump þ I rand ; I Na ¼ Àg Na m 3 hðV À V Na Þ; I K ¼ Àg K n 4 ðV À V K Þ; I L ¼ I KL þ I NaL þ I ClL ¼ Àg KL ðV À V K Þ À g NaL ðV À V Na Þ À g ClL ðV À V Cl Þ; dq=dt ¼ a q ð1 À qÞ À b q q; q ¼ m; n; h: Where n 4 and m 3 h represent the gating variables for delayed rectifier potassium (I K ) and transient sodium (I Na ) currents. The leak current (I L ) has three components: K + (I KL ), Na + (I NaL ), and chloride (I ClL ) leak. I pump is the net current due to the ATP-dependent pump that extrudes 3 Na + for bringing 2 K + in. A random current (I rand ) representing the background input from the other neurons is included in some simulations and is modeled by a zero mean Gaussian processes. The meaning and values of parameters used in the model are given in Table 1. The rate equations for the gating variables are [64] a m ¼ 0:1ðV þ 30Þ 1 À expðÀ0:1ðV þ 30ÞÞ ; We will use the instantaneous steady-state form of m, i.e. m = α m /(α m + β m ) [65].

Ion concentration dynamics
The ion current equations are augmented with dynamic intra-and extracellular concentrations of K + , Na + , and Cl − . These concentrations are modulated by the neuron's intrinsic ionic currents, Na + − K + pump current, and K + − Cl − cotransporters. The glial buffering and diffusion into the microenvironment of the cell (from the bath solution in slice preparation and vasculature in vivo) also modulate the K + concentration in the extracellular space. The ion concentrations inside and outside the cell are coupled to the membrane voltage equations via the Nernst equation. The rate equations for K + and Na + , O 2 , and Cl − concentrations and the rational for different terms within these equations are described in detail in [15,[66][67][68], [18], and [33] respectively and are summarized below. Given I K , I KL , I pump , diffusion of K + to the microenvironment (I diff ), glial buffering (I glia ), and K + − Cl − cotransporters (I KCC ), the extracellular K + dynamics, [K] o , can be represented in the model as where τ = 1000 is used to convert seconds to milliseconds. γ converts ionic current to concentration rate of change and is calculated using g ¼ A FÂVol [15], where A, Vol, and F represent cell area, volume, and Faraday constant respectively. β is the intra-to extracellular volume ratio. γ and β change when the volume and surface area of the cell change. In simulations where volume is treated as a dynamical variable, both γ and β change dynamically. Modeling glial K + buffering by a rate equation as done in [16] qualitatively did not change our results. Similarly, making I diff a function of β explicitly as done previously [33] did not change our conclusions (but see the discussion in "Anoxic depolarization" section).
The Na + − K + pump is modeled as a product of sigmoidal functions, ρ is the maximum pump strength, and [Na] i is the intracellular Na + concentration [15] [17] is given as Where [O 2 ] 1 is the oxygen concentration in the perfusion solution with a normal range of 30 − 32mg/L when aerated with 95%O 2 and 5%CO 2 at 32-34°C. α is a conversion factor that converts pump current in mM/s to O 2 concentration change in mg/L/s and the diffusion constant of O 2 ( O ) is obtained from Fick's law [18].
[K] o,1 in the diffusion equation (Eq 3) is the K + concentration in the nearby reservoir. Physiologically, this would correspond to either the bath solution in a slice preparation, or the vasculature in the intact brain (noting that [K] o is kept below the plasma level by trans-endothelial transport). In physiological conditions, [K] o,1 is set equal to 4mM. We use a previously published value for the diffusion constant of [K] o to the nearby reservoir K [66] which was obtained from Fick's law. That is, K = 2D/Δx 2 , where D = 250 × 10 − 6 cm 2 /sec is the K + diffusion constant in neocortex [69], and Δx % 20μm for brain reflects the average distance between neurons and capillaries [70].
Both active and passive K + uptake into glia is incorporated into a simplified single sigmoidal response function that depends on extracellular K + concentration with B glia representing the buffering strength [15]. Two factors allow the glia to provide a nearly insatiable buffer for the extracellular space. The first is the large size of the glial network. Second, the glial endfeet surround the pericapillary space, which through interaction with arteriole walls affecting blood flow, and transport of ions to the vascular space, amplifies the effective buffering capability of the glia [71][72][73].
Chloride is the primary permeant anion and its homeostasis is important for a range of neurophysiological processes. Neurons regulate intracellular chloride ([Cl − ] i ) through cationchloride cotransporters, especially the Na + /K + /2Cl − cotransporter (NKCC1) and K + /Cl − cotransporter (KCC2) [74]. In the embryonic and early postnatal brain, neurons show robust expression of NKCC1 but minimal expression of KCC2 [74]. In the mature brain, the expression of KCC2 increases, accompanied by a concurrent down regulation of NKCC1 expression [74]. KCC2 is important in maintaining low [Cl − ] i , resulting in hyperpolarizing GABA responses. Since KCC2 operates close to its thermodynamic equilibrium: [74], even a small increase in [K] o in the mature brain will reverse Cl − transport, from efflux to influx.
where ρ KCC is the maximum strength of KCC2 and estimated using the peak conductance given in [75]. In [33], we also included the Na + /K + /2Cl − (NKCC1) cotransporter, however, the results in this paper qualitatively remain the same without NKCC1 and hence it is excluded from the model. Nevertheless, NKCC1 should be included while modeling neurons from the embryonic or early postnatal brain [33,74] and perhaps from pathological conditions where there may be aberrant transporter expression [76]. Intracellular Na + concentration is controlled by transient Na + , Na + leak, and Na + − K + pump currents [15] d½Na i dt Previously, we assumed that the flow of Na + into the cell is compensated by the flow of K + out of the cell so that [K] i can be approximated by [15] We also assumed that the total amount of sodium is conserved, and hence only one differential equation for sodium is needed [15], so that are modeled by rate equations governed by various fluxes (see "Results" section) instead of the simplified version due to the conservations shown above. This is consistent with our previous findings where we demonstrated that the unification of SZ and SD is qualitatively preserved when using simplified conservation equations [33]. However, the simplification due to conservations is helpful in performing bifurcation analysis using numerical solvers based on continuation methods such as XPPAUT [77]. Such simplifications are common and especially important in neuronal models where multiple time-scales are involved as in our model [15,33,78].
The reversal potentials for K + , Na + and Cl − are updated based on the instantaneous ion concentrations using the Nernst equations The dynamic ion concentrations feedback into the Hodgkin-Huxley type Eqs (1, 2) through Eq (12). A more complex formulation would employ Goldman-Hodgkin-Katz voltage and current equations, which were found qualitatively similar in terms of unification in recent work [33].

Dynamic cell volume
Intense neuronal firing during SZ, SD, AD and the relevant changes in ion concentrations cause cell swelling. We use the formalism of [14] to model the change in cell volume where c Vol and Vol initial are the instantaneous and initial volumes of the cell. π o and π i are the sums of all ion concentrations outside and inside the cell respectively, i.e.  [79], we implement the change in the cell volume as a first-order process.
The time-scale of volume change, τ v , is important for the spontaneous SZ-SD transition. For the transition to happen, the recovery of the cell from the swollen state would have to be slower than the electrical response changes during a SZ. The cell would also have to be in a SZ state for long enough to produce this effect, and that lower levels of [K] o , with shorter duration seizures, would keep the cell out of SD. If the cell were to swell and shrink very fast so that it would recover its pre-seizure volume before going into next seizure, the cell would only exhibit spontaneous periodic seizures without going through the SZ-SD transition. To ensure the robustness of our results and that the spontaneous SZ-SD transition seen in our model occurs even with faster volume changes, we considered τ v = 50ms, significantly smaller than the previously used phenomenological value of 250ms [33,79]. Nevertheless, we tested a wide range of τ v values and found that the model exhibits spontaneous SZ-SD transition for 40ms τ v 360ms. For τ v < 40ms, the cell recovers from swelling before going into seizure for the second time. For τ v > 360ms, cell still exhibits spontaneous SZ-SD transition, however, it seizes more than once before going into SD state due to the slower change in volume. The number of SZs before the cell transitions into SD increases as we increase τ v . For example, for τ v = 400ms, 500ms, and 1000ms, the cell seizes twice, thrice, and seven times respectively before making a transition to SD. While, the spontaneous SZ-SD transition itself is a robust phenomenon, τ v > 360ms causes the cell to seize multiple times prior to the transition. The rest of the results including all bifurcation diagrams in the paper qualitatively remain the same as we change τ v . In simulations where intracellular volume is treated as variable, the conductance densities for Na + , K + , and Cl − currents (g Na etc) are modified so that the total conductance for a channel type over the whole cell remains fixed even though the conductance per unit area is changing Where A ss and A ins are the steady state (initial) and instantaneous cell surface areas respectively. g X and g X are the conductance densities in case of fixed and changing volumes respectively. The effect of dynamic volume changes on membrane capacitance is negligible and not included. For example, the total capacitance of the cell during the spontaneous transitions between SZ and SD shown in Fig 2 increased from initial 2.9074pF to 3.0295pF as the cell swelled, a variation that did not change the results qualitatively.
The dynamic change in volume causes the concentration of a given ion specie to change for a fixed number of ions. Thus, the dynamic volume leads to an additional flux term in the ion concentration equations.  (5,10), and infra-slow (Vol; Eq 14). Although we focus on the implications of dynamic intracellular volume for cell behavior in pathological states in this study, we will also analyze the evolution of other slow variables during these states and how they effect cell pathology.

Numerical methods
The coupled differential equations are solved in Fortran 90 using 4th order Runge-Kutta method. The steady states in Fig 1A, and the bifurcation diagrams in Figs 5 and 6 are generated using XPPAUT software [77]. In the bifurcation diagrams of [K] o (and other slow variables), XPPAUT could only capture the small amplitude oscillations in [K] o due to individual membrane potential spikes. It could not capture the high amplitude low frequency oscillations due to periodic SZ and SD events. Therefore, we simulated the periodic orbits in Figs 1A, 3, S1, and S2 Figs using Fortran 90. The codes reproducing key results in the paper are given in Supporting Information Text S1.