article: 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.
  Cells swell during a wide variety of pathologies, including trauma, ischemia, hypoxia, seizures, and spreading depression [1–3]. Changes in osmolality can change the susceptibility to epileptiform activity [4–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–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 (O2) and glucose deprivation (OGD) that may or may not recover depending on the cell type after O2 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–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–31], characteristic ionic changes [13–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 thalamus-hypothalamus 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 O2 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.
  We investigate the role of cell size and relative (to intracellular volume) extracellular space on neuronal behavior by varying the radius of the cell, rin, keeping the total radius of cell and extracellular space, rtot, fixed. Thus changing rin 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, ∞ = 8mM—a value typically used for inducing SZ in in vitro [30]. The behavior of the cell changes dramatically as we increase rin. A bifurcation diagram showing the maximum and minimum of [K]o as a function of rin is shown in Fig 1A. Briefly, we fixed rin and ran the simulation generating a time-trace representing [K]o versus time. Depending on the rin 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 rin in Fig 1A represent the maxima and minima of [K]o oscillations for that particular value of rin. This process was repeated many times, each time incrementing rin by a small amount. For rin values, where [K]o does not oscillate (i. e. [K]o converges to a steady state), the maxima and minima have the same value and is represented by a line (stable, red; unstable, blue). To capture both stable and unstable behaviors, the steady states were simulated in XPPAUT. For all rin < 4. 66μm, [K]o remains unchanged and the cell remains in steady state with V ≈ −60mV (Fig 1B). As we increase rin 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 rin = 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 small-amplitude spikes on its way out of the depolarized state (Fig 1D). As we increase rin 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 rin > 4. 924μm where V remains fixed at a stable resting value. 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 rin. At rin > 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 rin 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 rin 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 rin if β is kept constant (S2 (A) Fig). Using β as a bifurcation parameter at fixed rin on the other hand causes the cell to make the transitions between steady state, SZ, and SD (S2 (B) Fig). Nevertheless, the cell size per se is an important parameter that together with β shapes the bifurcations and the parameter ranges where different behaviors are observed (compare Fig 1A and S2 Fig). 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 rin 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 rin can change enough during one SZ so that the cell would make the transition from SZ to SD regions. A cell with rin = 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 rin = 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. Fig 2B shows the variations in [Na]i during this transition. The thick arrow in Fig 2A indicates that the cell would go into a depolarized state similar to ischemia-induced AD if [K]o increases further either due to swelling or lack of O2. This point will be further elaborated later in this paper. 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 rin and glial K+ buffering strength, Bglia. In Fig 3A, we show the maximum of [K]o as a function of rin and Bglia. As is clear, for weaker glial buffering the cell’s behavior changes as a function of rin in the same manner as in Fig 1. However, for stronger glial buffering (Bglia > 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 rin and Bglia. For larger Bglia, 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 [O2] from simulations in Fig 3 shows that the resting state at large rin (right side in Fig 3B) is energetically favorable as compared to the one at small rin (left side in Fig 3B) (Fig 4). Interestingly, the local O2 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) (not shown) show that their behavior contrasts with [K]o. That is, the amplitudes of [Na]i, [K]i, and [Na]o oscillations decrease with increasing rin. Because the pumps are stimulated by [Na]i and [K]o (Eq 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]. To gain further insights into the mechanisms behind the three behaviors in Fig 1C-1E exhibited by the cell, we performed a bifurcation analysis of the model fixing the slower variables [K]o, [Na]i, [Cl]i, [O2], and volume (Vol). In Fig 5, we show the maxima and minima of V as a function of [K]o. For small [K]o values, the cell remains in stable steady state (SSS1; red line on the left) which collides with an unstable steady state (USS; black dotted line) giving rise to a periodic orbit (PO; green circles) through a SNIC bifurcation (Fig 5A). The unstable steady state gains stability (SSS2; red line on the right) through a Hopf bifurcation (HB) where the membrane potential of the cell is locked into a depolarized state. The SNIC moves to the right as we increase [Na]i (Fig 5B) and bifurcates into a saddle homoclinic (HC) bifurcation (Fig 5C, 5D). The locations of special points, HB, SNIC, the limit point (LP) (Fig 5A), and HC (Fig 5C) as [K]o and [Na]i vary simultaneously are shown through a two-parameter bifurcation diagram in Fig 6. The limit cycles and steady states from the full (but with fixed volume) model at different rin values are also shown in Fig 6. For clarity, we will restrict our discussion to the case of [Cl]i = 8mM but the argument applies to the other values of [Cl]i as well. During the SZ event like the one shown in Fig 1C (rin = 4. 82μm), the cell enters the periodic orbit from steady state on the left (SSS1 in Fig 5A) through a SNIC bifurcation and spikes for a few seconds before going back to SSS1. The [K]o values when the trajectory crosses the SNIC bifurcation curve indicate the beginning and end of the burst in the time trace. Since the burst terminates at a SNIC, the frequency of the burst is low at the beginning and end of the burst. This limit cycle with bursts never generates [K]o substantially above a physiological ceiling, about 12 mM experimentally [42] and about 13–14 mM in Fig 6, below which spikes and seizures are observed but not SD. Indeed, the model generates rather pure seizure dynamics with [K]o below the ceiling. At rin = 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 rin = 4. 84μm also applies to the case with rin = 4. 9μm (Fig 1E), except that the HB to HB trajectory loop is more dramatic than rin = 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 rin = 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). 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 O2 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, [O2]∞, 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 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 O2 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 space and bath solution (S3 (B) Fig). Because diffusion of ions within the extracellular space is a function of the size and geometry of extracellular space [47], we needed to make diffusion a function of O2 [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 O2 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 O2 and glucose. In the absence of O2, 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 (Kir) 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 Kir 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 O2, 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 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 O2 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 O2 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) to compare the cell dynamics during SZ, SD, and AD. The change in the membrane potential along with [K]o and percentage change in the cell volume during the simulation represented by the dashed line in Fig 9A is shown in Fig 9B. The gray planes in Fig 9A indicate the approximate regions for SZ, SD, and AD. The cell exhibits SZ below the physiological [K]o ceiling represented by the bottom plane, SD between bottom and top planes, and AD above top plane. The position of these planes will change with the size and density of neurons.
  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 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–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 mentioned above, there are four different time-scales in our model: fast (variables V, n, and h; Eq (1) ), intermediate ([K]o, [Na]i; Eqs (3,7), slow ([O2], [Cl]i; Eqs (5,10) ), and infra-slow (Vol; Eq (14) ). Strictly speaking, the ultra-fast m-gate (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 slow-fast 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 O2 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 O2 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 O2 [61], other pathways such as Ca2+-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–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 Ca2+ 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.
  The membrane potential V of the cell is modeled with the following set of modified Hodgkin-Huxley type equations [15,63] C d V d t=I N a + I K + I L - I p u m p + I r a n d, I N a=- g N a m 3 h (V - V N a), I K=- g K n 4 (V - V K), I L=I K L + I N a L + I C l L = - g K L (V - V K) - g N a L (V - V N a) - g C l L (V - V C l), d q / d t=α q (1 - q) - β q q, q = m, n, h. (1) Where n4 and m3 h represent the gating variables for delayed rectifier potassium (IK) and transient sodium (INa) currents. The leak current (IL) has three components: K+ (IKL), Na+ (INaL), and chloride (IClL) leak. Ipump is the net current due to the ATP-dependent pump that extrudes 3 Na+ for bringing 2 K+ in. A random current (Irand) 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] α m = 0. 1 (V + 30) 1 - e x p (- 0. 1 (V + 30) ), β m = 4 e x p (- V + 55 18), α n = 0. 01 (V + 34) 1 - e x p (- 0. 1 (V + 34) ), β n = 0. 125 e x p (- V + 44 80), α h = 0. 07 e x p (- V + 44 20), β h = 1 1 + e x p (- 0. 1 (V + 14) ). (2) We will use the instantaneous steady-state form of m, i. e. m = αm/ (αm + βm) [65]. 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+, O2, and Cl− concentrations and the rational for different terms within these equations are described in detail in [15,66–68], [18], and [33] respectively and are summarized below. Given IK, IKL, Ipump, diffusion of K+ to the microenvironment (Idiff), glial buffering (Iglia), and K+ − Cl− cotransporters (IKCC), the extracellular K+ dynamics, [K]o, can be represented in the model as d [ K ] o d t = 1 τ (γ β (I K + I K L - 2 I p u m p) - I d i f f - I g l i a + I K C C), I p u m p = 𝓕 ([ O 2 ]) ρ (1 1 + e x p ( (25 - [ N a ] i) / 3) ) (1 1 + e x p (8 - [ K ] o) ), I d i f f = 𝓕 ([ O 2 ]) ϵ K ([ K ] o - [ K ] o, ∞), I g l i a = 𝓕 ([ O 2 ]) B g l i a 1 + e x p ( (18 - [ K ] o) / 2. 5), (3) where τ = 1000 is used to convert seconds to milliseconds. γ converts ionic current to concentration rate of change and is calculated using γ = A F × V o l [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 Idiff 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]. Each sigmoidal term saturates for high values of [Na]i and [K]o respectively. The ATP required to keep the pump running depends on the local O2 availability. The ATP concentration and hence the pump strength decreases as the cell depletes its local O2 reservoir. We use a sigmoid function 𝓕 ([O2]) to model the O2 concentration ([O2]) -dependence of the pump activity [18] 𝓕 ([ O 2 ]) = 1 1 + e x p ( (16 - [ O 2 ]) / 4). (4) The justification for the [O2]-dependence of Idiff and Iglia terms is elaborated in [33] and further discussed in the Results section. For simplicity, we use the same functional form as in Ipump for the [O2]-dependence of Idiff and Iglia. Using separate functions for the [O2]-dependence of these three fluxes as in [33] does not qualitatively change our results. The local O2 concentration is controlled by the activity of the Na+ − K+ pump and diffusion of O2 from the bath solution (or vasculature) to the extracellular space. Thus the rate equation for [O2] as developed in [18] based on [O2] imaging experiments [17] is given as d [ O 2 ] d t = 1 τ (- α γ I p u m p + ϵ O ([ O 2 ] ∞ - O 2) ). (5) Where [O2]∞ is the oxygen concentration in the perfusion solution with a normal range of 30 − 32mg/L when aerated with 95%O2 and 5%CO2 at 32–34°C. α is a conversion factor that converts pump current in mM/s to O2 concentration change in mg/L/s and the diffusion constant of O2 (ϵO) is obtained from Fick’s law [18]. [K]o, ∞ 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, ∞ 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/Δx2, where D = 250 × 10−6cm2/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 Bglia 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–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 cation-chloride 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: [Cl−]i = [Cl−]o[K+]o/[K+]i (i. e. ECl = EK) [74], even a small increase in [K]o in the mature brain will reverse Cl− transport, from efflux to influx. [K]i, [Cl]i, and [Cl]o are intracellular K+, Cl−, and extracellular Cl− concentrations respectively. KCC2 in our model is formulated by a logarithmic function [26,33] I K C C = ρ K C C l o g ([ K ] i [ C l ] i [ K ] o [ C l ] o), (6) 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 [ N a ] i d t = 1 τ γ (I N a + I N a L - 3 I p u m p). (7) 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] [ K ] i = 140 mM + (18 mM - [ N a ] i). (8) We also assumed that the total amount of sodium is conserved, and hence only one differential equation for sodium is needed [15], so that [ N a ] o = 144 mM - β ([ N a ] i - 18 mM), (9) where 140mM, 18mM, and 144mM in the above equations are the normal resting values of [K]i, [Na]i, and [Na]o respectively. The intracellular Cl− concentration is modeled by Cl− leak and K+ − Cl− cotransporter [33] d [ C l ] i d t = 1 τ (γ I C l L - I K C C β), (10) and the extracellular Cl− concentration, [Cl]o, is set according to the conservation of charge due to Na+, K+, Cl−, and Ca+ ions in the extracellular space i. e. [ C l ] o = [ K ] o + [ N a ] o + 2. 0 [ C a ] o. (11) Where [Ca]o = 1mM is the extracellular calcium (Ca2+) concentration. The main conclusions in this paper qualitatively remain the same when [K]i, [Na]o, and [Cl]o 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 V k = 26. 64 ln ([ K ] o [ K ] i), V N a = 26. 64 ln ([ N a ] o [ N a ] i), V C l = 26. 64 ln ([ C l ] i [ C l ] o). (12) 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]. 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 V o l ^ = V o l i n i t i a l × (1. 1029 - 0. 1029 × e x p ( (π o - π i) / 20) ), (13) where V o l ^ and Volinitial 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. π i = [ N a ] i + [ C l ] i + [ K ] i + [ A ] i + [ C a ] i, π o = [ N a ] o + [ C l ] o + [ K ] o + [ A ] o + [ C a ] o. [A]i = 132. 1mM [79] and [A]o = 18mM are the intra- and extracellular concentrations of the impermeable anions. The value of [A]o is based on charge balance under resting conditions using [Cl]o = 132mM, and [Ca]i = 100nM is the intracellular Ca2+ concentration. Following [79], we implement the change in the cell volume as a first-order process. d V o l d t = 1 τ v (V o l ^ - V o l) (14) 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 (gNa 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 g X = g ¯ X (A s s / A i n s). (15) Where Ass and Ains are the steady state (initial) and instantaneous cell surface areas respectively. g ‾ X and gX 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. Thus the rate equation for [K]o (Eq 3) becomes d [ K ] o d t = 1 τ (γ β (I K + I K L - 2 I p u m p) - I d i f f - I g l i a + I K C C) - (d V o l o d t) ([ K ] o V o l o). Where Volo is the extracellular volume. Similarly, the rate equations for [Na]i (Eq 7) and [Cl]i (Eq 10) also have additional terms equal to (− d V o l d t) ([ N a ] i V o l i) and (− d V o l d t) ([ C l ] i V o l i) respectively. However, these fluxes do not change the results presented in this paper qualitatively and are not considered. Notice that there are four different time-scales in our model: fast (variables V, n, and h; Eq (1) ), intermediate ([K]o, [Na]i; Eqs (3,7) ), slow ([O2], [Cl]i; Eqs (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. 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.
summary: Massive rearrangement of ions across the plasma membrane and changes in cellular volume are common features of states such as seizures, spreading depression, and ischemia. In this paper, we focus on how volume itself influences neuronal activity. We build a unified computational framework for a wide range of neuronal behaviors by exploiting their previously unexplored common features. By combining the dynamic ion concentrations and volume, conservation of charge, and the energy requirements of the cell within a Hodgkin-Huxley type framework, we demonstrate the feasibility of a comprehensive framework encompassing a wide range of neuronal behaviors. We show the spontaneous transition of a neuron between seizure and spreading depression when the cell swells and contracts in response to varying osmotic pressure as a result of the rearrangement of different ions. Our model closely reproduces anoxic depolarization and relevant neuronal swelling during ischemia and reveals a dynamical definition for the experimentally described physiological ceilings that demarcate seizure from spreading depression and spreading depression from anoxic depolarization. This study opens up a new way of studying neuronal behavior where different states need not be treated separately but rather as a dynamical continuum of the neuronal membrane potential and its microenvironment.
title: The Role of Cell Volume in the Dynamics of Seizure, Spreading Depression, and Anoxic Depolarization