Abstract
The pathological phenomena of seizures and spreading depression have long been considered separate physiological events in the brain. By incorporating conservation of particles and charge, and accounting for the energy required to restore ionic gradients, we extend the classic Hodgkin–Huxley formalism to uncover a unification of neuronal membrane dynamics. By examining the dynamics as a function of potassium and oxygen, we now account for a wide range of neuronal activities, from spikes to seizures, spreading depression (whether high potassium or hypoxia induced), mixed seizure and spreading depression states, and the terminal anoxic “wave of death.” Such a unified framework demonstrates that all of these dynamics lie along a continuum of the repertoire of the neuron membrane. Our results demonstrate that unified frameworks for neuronal dynamics are feasible, can be achieved using existing biological structures and universal physical conservation principles, and may be of substantial importance in enabling our understanding of brain activity and in the control of pathological states.
Introduction
The electrical hyperactivity of seizures was first described by Berger (1933), whereas the propagating silencing of electrical activity in spreading depression (SD) was discovered by Leão (1944). These two pathological phenomena have long been considered separate physiological events (Somjen, 2004), characterized by their different patterns of neuronal activities (Leão, 1944; Traynelis and Dingledine, 1988; Žiburkus et al., 2006), characteristic ionic changes (Kager et al., 2000, 2002, 2007; Cressman et al., 2009; Barreto and Cressman, 2011; Krishnan and Bazhenov, 2011), and known interactions with oxygen (Czéh et al., 1993; Ingram et al., 2013). Seizures propagate rapidly, whereas the depolarization in SD propagates slowly (Somjen, 2004). In humans, SD triggers migraine, and is increasingly recognized as adversely affecting the outcome of head trauma, cerebral hypoxia, and stroke (Lauritzen et al., 2011). Here, we show that these pathological dynamics, as well as normal spiking behavior of neurons, share a unified dynamical underpinning.
Both seizures and SD are induced by similar mechanisms: hypoxia, hypoglycemia, neural injury, high concentrations of K+, or Na+/K+ pump inhibition (Dreier, 2011). In in vitro experiments, increasing potassium in the bath perfusate to 8.5 mm (normal 3.5 mm) induces spontaneous periodic seizures (Traynelis and Dingledine, 1988). Elevating potassium further to 26 mm (Anderson and Andrew, 2002) or 40 mm (Zhou et al., 2010) evokes SD. SD can also be initiated by injecting highly concentrated KCl into the tissue, or by placing a KCl crystal on the brain (Somjen, 2001). Hypoxic SD (HSD) is observed after reducing blood flow or oxygen (Czéh et al., 1993; Jing et al., 1994). If blood or oxygen is restored, HSD terminates and normal function returns; if not restored, an electrical “wave of death” ensues (Zandt et al., 2011). In the intact brain or slice in fully oxygenated perfusate, local extracellular tissue potassium concentration ([K+]o) stays below a “ceiling” of 8–15 mm during seizures (Somjen, 2004). Seizures can be associated with the leading wave-front of SD, and such mixed states are typically seen in HSD (Czéh et al., 1993).
Seizures and SD are observed as network phenomena, but as our biophysical models of these phenomena have improved, the concept that single cells can generate the underlying patterns of these pathological activities has become increasingly accepted (Kager et al., 2000; Cressman et al., 2009). In intracellular recordings during seizures, excitatory and inhibitory cell types can simultaneously exhibit either seizure or SD dynamics, respectively (Žiburkus et al., 2006).
By incorporating conservation of particles and charge, and accounting for the energy required to restore ionic gradients, we extend the classic Hodgkin–Huxley formalism to uncover a unification of neuronal membrane dynamics. We now account for a wide range of neuronal activities, from spikes to seizures, spreading depression (whether high potassium or hypoxia induced), mixed seizure and spreading depression states, and the terminal anoxic wave of death.
Materials and Methods
Our mathematical model builds on our previous work. In Cressman et al. (2009), we extended the Hodgkin–Huxley formalism to incorporate the dynamic concentrations of Na+ and K+ ions. Intracellular Na+ concentration ([Na+]i) was controlled by transient Na+ current, Na+ leak, and Na+/K+ exchange pump. Extracellular K+ concentration ([K+]o) was controlled by delayed rectified K+ current, K+ leak, Na+/K+ pump, glial buffering, and diffusion to the blood vessels. The leak current also included a component due to Cl−. Intracellular K+ ([K+]i) and extracellular Na+ ([Na+]o) were modeled by using the electroneutrality and conservation of Na+ ions respectively. Although Cl− was included in the model, the concentration of Cl− inside and outside the cell was assumed to be constant. In Cressman et al. (2009), we characterized the bifurcation structure of this extension to the Hodgkin–Huxley formalism based on potassium.
In recent work (Wei et al., 2014), we demonstrated that oxygen could also serve as a bifurcation parameter in the Hodgkin–Huxley formalism, by introducing an oxygen dependency on the Na+/K+ exchange pumps.
In this present work, we sought to examine the combination of potassium and oxygen dependency studied in the previous works by Cressman et al. (2009) and Wei et al. (2014). We add oxygen dynamics and incorporate the energy (oxygen) dependence of the Na+/K+ exchange pumps that contributes to the membrane current balance equation and glial K+ buffering. We also add dynamics to the concentrations of intracellular and extracellular Cl− concentrations ([Cl−]i and [Cl−]o, respectively) that involves cotransporters in addition to Cl− leak current. By keeping track of all of the ion fluxes, we calculate the osmolality changes and the resulting cellular volume changes. The rate equations for membrane potential and ion concentrations are modified from Cressman et al. (2009) accordingly to incorporate these changes. We consider three versions of the model: (1) the full model where the [Na+]o, [K+]i, and [Cl+]o are modeled as differential equations (unlike Cressman et al., 2009), (2) a simplified model where these three concentrations are given by electroneutrality and conservation laws (similar to Cressman et al., 2009), and (3) a model version where various current (fluxes) are given by Goldman–Hodgkin–Katz formalism. In the following, we discuss these three models in detail.
The full model
Membrane potential dynamics.
We used a single compartment model containing transient sodium currents, delayed rectifier potassium currents, specific leak currents for sodium, potassium and chloride ions, and sodium/potassium pump current. The membrane potential, V, of the neuron is modeled by modified Hodgkin-Huxley equations:
where GNa, GK, GNaL, GKL, and GClL represent Na+ voltage-gated maximal conductance, K+ voltage-gated maximal conductance, Na+ leak conductance, K+ leak conductance, and Cl− leak conductance, respectively. INa and IK are the sodium and potassium currents passing through the voltage-gated ion channels including the sodium and potassium leak current, respectively. ICl is the chloride leak current. Ipump is the net current due to the Na+/K+ ATP-dependent pump that stoichiometrically intrudes 2 K+ and extrudes 3 Na+, which is a missing part in the classic Hodgkin–Huxley formalism. γ = S/(Fυi) is a conversion factor from the current units (μA/cm2) into the concentration units (mm/s), where S, υi, and F are the surface area of the cell, intracellular volume, and Faraday constant, respectively.
The activation and inactivation variables m, h, and n vary between 0 and 1, and represent the fraction of ion selective channels in the closed and open states. The parameters αm, βm, αh, βh, αn, and βn are opening and closing rates of the ion channel state transitions that are dependent on V. The equations of these rates are from a pyramidal cell model (Gloveli et al., 2005), originally derived from a model of hippocampal neurons (Traub et al., 1991), shown as follows:
The reversal potentials of Na+ (ENa), K+ (EK), and Cl− (ECl) are given by Nernst equations:
where [·]i and [·]o represent concentrations inside and outside the cell, respectively. The ECl quotient was reversed to account for the different valencies. The units and description of all parameters are summarized in Table 1. Unlike the Hodgkin–Huxley equations where various ion concentrations are fixed, we have incorporated potassium, sodium, and chloride ion dynamics that govern the reversal potentials.
Units and description of the parameters used in the model
Ion concentration dynamics.
The concentration of each ion type is continuously updated by integrating the relevant ion currents and fluxes in the manner of (Lee and Kim, 2010). The rate of change of the number of extracellular potassium ions, 
where τ = 1000 is used to convert second to millisecond. υi and υo are the intracellular and extracellular volumes, respectively. β = υi/υo is the ratio of intra-/extracellular volume. The ion concentrations are calculated by the ion number over the volume within the compartment, that is [·]i = Ni/υi, [·]o = No/υo, where i indicates inside and o indicates outside of the cell.
The functional forms of neuronal Na+/K+ pump current (Ipump), glial buffering current (Iglia), and potassium diffusion current (Idiff) are taken from our previous work (Cressman et al., 2009) and modified below to incorporate the O2 dependence of these fluxes. For the Na+/K+ pump on glia (Igliapump), we use the same functional form as used for neuronal pump but with 1/3 the strength of neuronal pump. We use this ratio due to the fact that the relative resting energy consumption in neurons versus glia is ∼3:1 (Attwell and Laughlin, 2001). Iglia represents the combined effect of glial K+ uptake through inward rectified K+ channels and Na+/K+/2Cl− cotransporters.
where ρ, Gglia, ϵk, and [K+]bath represent pump strength, glial uptake strength, potassium diffusion coefficient, and bath potassium concentration, respectively. ρ = ρmax, Gglia = Gglia,max, ϵk = ϵk,max for the fully oxygenated state with normal bath potassium. We assume a fixed intracellular sodium concentration [Na+]gi = 18 mm in the glial compartment.
Chloride is the primary permeant anion and its homeostasis is important for a range of neurophysiological processes. Neurons regulate the intracellular chloride ([Cl−]i) by cation-chloride cotransporters, especially the Na+/K+/2Cl− cotransporter (NKCC1) and K+/Cl− cotransporter (KCC2; Payne et al., 2003). In the embryonic and early postnatal brain, neurons show robust expression of NKCC1 but minimal expression of KCC2 (Payne et al., 2003). In the mature brain, the expression of KCC2 increases, accompanied by a concurrent downregulation of NKCC1 expression (Payne et al., 2003). KCC2 is important in maintaining low [Cl−]i, resulting in hyperpolarizing GABA responses. Because KCC2 operates close to its thermodynamic equilibrium: [Cl−]i = [Cl−]o[K+]o/[K+]i (i.e., ECl = EK; Payne et al., 2003), even a small increase in extracellular K+ in the mature brain will reverse Cl− transport, from efflux to influx.
The NKCC1 cotransporter (Inkcc1) and KCC2 cotransporter currents (Ikcc2) are modeled in a Nernst-like fashion (Østby et al., 2009) as follows:
where Ukcc2 = 0.3 mm/s and Unkcc1 = 0.1 mm/s are cotransporter strength (Payne et al., 2003; Østby et al., 2009, 2012) estimated using the peak conductances given by Lauf and Adragna (2001). The NKCC1 and KCC2 transporters mediate Cl− transport without any net charge movement across the membrane. In the model, KCC2 is expressed at high relative levels as in an adult brain, maintaining low [Cl−]i in the resting steady-state. However, with only KCC2 in the model, [Cl−]o increased during spreading depression as the extracellular space decreased in volume, in contradiction with experimental findings (Hansen and Zeuthen, 1981). We therefore needed to have KCC2 and NKCC1 balance each other realistically as they do to regulate [Cl−]i in real neurons (Nardou et al., 2013). It is known that high levels of [K+]o causes activation of NKCC1 in glial cells (Su et al., 2002). We assumed that neuronal NKCC1 has similar dynamics in our neurons, and incorporated this behavior by making NKCC1 activated by high [K+]o as given by function f([K+]o). This helped to maintain a physiologically plausible and experimentally consistent [Cl−]o during spreading depression in our model. Incorporating NKCC1 on the glial membrane and permitting it to add to the clearance of high [K+]o into the glial compartment would further enhance the physiological fidelity of future iterations of this model.
Oxygen concentration dynamics.
Support of neural spiking consumes the most oxygen in the brain (Lennie, 2003). In our single compartment model, the oxygen consumption can be estimated by the activity of Na+/K+ pumps that transport 3 Na+ outward with 2 K+ inward against their concentration gradients for each ATP hydrolyzed (Erecińska and Dagani, 1990; Attwell and Laughlin, 2001; Lennie, 2003). Typically, a spike may consume 2.4 × 109 molecules of ATP, 91% of which is consumed on Na+/K+ pumps (Lennie, 2003). As neural activity increases so does oxygen and ATP utilization.
The extracellular oxygen concentration, [O2]o, around a single neuron is supplied diffusively by the bath solution in in vitro experiments. Therefore, [O2]o can be modeled as follows:
where [O2]bath is the bath oxygen concentration in the perfusion solution with a normal value ∼32 mg/L, when aerated with 95% O2: 5% CO2 at 32−34°C. The diffusion constant, ϵo, is obtained from Fick's law, ϵo = D/Δx2. We used a diffusion coefficient D = 1.7 × 10− 5 cm2/s for oxygen in brain tissue (Homer et al., 1983) and Δx = 100 μm for the average distance from electrode tip to the surface of the slice. α converts charge utilization in the pump current (mm/s) to the rate of oxygen concentration change (mgL−1s−1).
The detailed calculation of α is shown as follows. When the Na+/K+-ATPase pump current is 1 mm/s, it transports 3 mm/s Na+ outward and 2 mm/s K+ inward. The amount of ATP required to be hydrolyzed for this process is 1 mm/s. The pump is fueled primarily by oxidative phosphorylation, which yields up to 36 molecules of ATP from the complete oxidation of 1 glucose with 6 oxygen molecules: C6H12O6 + 6 O2 → 6 CO2 + 6 H2O + 36 ATP.
The amount of oxygen needed to generate 1 mm/s ATP is 1/6 mm/s. Because the molar mass of O2 is 32 g/mol, the concentration of oxygen expended on 1 mm/s pump current is 5.3 mgL−1s−1 O2. Therefore, the conversion factor α between pump current (mm/s) and oxygen concentration change (mgL−1s−1) is 5.3 g/mol.
With oxygen dynamics in the model, we have modified the Na+/K+-ATP pump rate (ρ) in Equation 5 using a sigmoid function of [O2]o (Petrushanko et al., 2007) as follows:
Both the neuronal and glial Na+/K+-ATP pumps are modified so that they depend on the available oxygen concentration in the extracellular space. The pump strength decreases as the cell depletes its local oxygen reservoir. This dynamic microenvironment is essential to account for experimental findings.
In the intact brain, glial cells are not only in close proximity to neurons, but also contact blood microvessels and form a substantial component of the blood brain barrier. Ion transport between blood vessels and glial cells is dependent upon active transport through the Na+/K+ pump (Gloor, 1997). To simulate the anoxic wave of death, the pump strength (ρ), the uptake of K+ ions by the glial cells (Gglia), and diffusion of K+ to the blood (ϵk) were previously set to be zero (Zandt et al., 2011). To better reflect the oxygen dependency of Gglia and ϵk, we modified them in Equation 5 according to a sigmoid function of oxygen concentration in the blood vessels, or in our case in vitro the bath concentration, as follows:
Thus, ρ, Gglia, and ϵk will approach zero as seen physiologically when [O2]bath is extremely low.
Several observations support the way we have modeled the dependence of K+ glial buffering and extracellular space diffusion in hypoxia. In the absence of O2, astrocytes attempt to buffer the increased extracellular K+ by switching to anaerobic glycolysis and may swell substantially (Arumugam et al., 2010), restricting diffusion of K+ and limiting glial energy reserves. Astrocytic inward rectifying K+ channels (Kir) also contribute to K+ siphoning. Kir3.x, a variant of Kir, gates K+ through interaction with G-protein coupled receptors (GPCRs). The activation of GPCRs by its ligand results in the release of intracellular effector molecules Gα and Gβγ from inactive G-protein complexes, Gαβγ. Gβγ complexes in turn interact with Kir3. x to open them, allowing K+ to flow into astrocytes (and neurons; for review of Kir channels, see Hibino et al., 2010). Na+/K+/Cl− cotransporters 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 (Witthoft et al., 2013) 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 two-way K+ trafficking at the blood–brain barrier occurs at the junctions between astrocytic end-feet and blood vessels (Gloor, 1997). Astrocytes release K+ next to tight-junction sealed endothelial cells that are tightly packed around blood vessels. Na+/K+ pumps transfer that K+ into the endothelial cells and it is then extruded into the blood vessels through K+ channels. A reverse process transfers K+ from blood vessels to astrocytes and finally to neurons. The lack of ATP would disrupt this pathway, consequently impairing the K+ diffusion between blood vessels and extracellular space. In an in vitro brain slice, there is no measurement data to support the precise nature of the glial buffer extrusion to the bath external to the slice, although the K+ transfer to the vascular system must be shut down. To a degree, diffusion and glial buffering can substitute for each other in clearing K+ from the extracellular space (Cressman et al., 2009; Wei et al., 2014). Our strategy is to keep to in vivo physiological fidelity in this modeling, and we anticipate that there will be some residual leakage of K+ to bath unaccounted for in our model at extremely low oxygen pressures for slice experiments. We note that none of the unification phenomena that we focus on in this paper is dependent upon how we model this extremely low oxygen regime.
Volume dynamics.
The dynamics of cell volume is modified from the work of Kager et al. (2002, 2007). The extracellular volume is 14.29% of the initial cellular volume at the resting state in the absence of osmotic pressure gradients and maximally shrinks to 4% of the initial cellular volume at spreading depression (Jing et al., 1994). The cellular volume was therefore not allowed to expand beyond the limits (110.29% of the initial cell volume) imposed by the extracellular volume in the model, so that the total sum of extracellular volume and intracellular neuronal volume is constant at 114.29% of the initial neuronal volume after (Kager et al., 2002, 2007):
where υ̂i represents the expected intracellular neuronal volume based on the difference between extracellular osmotic pressure (πo) and intracellular osmotic pressure (πi). The extracellular fluid and intracellular fluid are composed of water and several ionic species: sodium ions (Na+), potassium ions (K+), chloride ions (Cl−), and miscellaneous impermeant ions (proteins, amino acids, phosphates, etc.), with a net negative charge, A−.
where [A−]i = 132 mm and [A−]o = 18 mm are the intra- and extracellular concentration of the impermeable anions, calculated by assuming that the initial osmotic pressure gradient is zero. We implemented the change of cell volume as a first-order process with a time-constant of 250 ms (Kager et al., 2002, 2007):
The extracellular volume υo can be calculated by assuming constant total volume.
where β0 = 7 is the initial ratio of initial intra-/extracellular volume, and υi0 is the initial volume of the cell with a radius of 7 μm. The instantaneous ratio of intracellular/extracellular volume β = υi/υo is updated depending on the instantaneous volume values.
With volume dynamics in the model, diffusion of potassium in the extracellular space of the brain is constrained by the volume fraction β (Syková and Nicholson, 2008). To reflect this observation, we modify the functional form of ϵk in Equation 9 by multiplying it with a sigmoidal function of β:
The simplified model
In performing bifurcation analyses using numerical solvers using so called continuation methods (Allgower and Georg, 1980), it is very helpful to simplify the equations, as well as eliminate some of the fast dynamics (Rubin and Terman, 2004; Cressman et al., 2009) so that the solver can more effectively follow the continuous variation in a parameter (Ermentrout, 2002). In neuronal models, this is especially important whenever fast bursting or seizing dynamics may be encountered. We also need to hold the parameters not being followed constant, so that the equations are near a steady-state fixed point or limit cycle. Therefore to perform the analytical bifurcation analysis of our model using XPPAUT (Ermentrout, 2002) in Figure 2E,F, we used a simplified model where three ion concentrations, [K+]i, [Na+]o, and [Cl−]o, are constrained by the electroneutrality of the medium inside the cell and conservation of Na+ and Cl− inside and outside the cell.
In addition to replacing three of our dynamical equations (Eq. 4) by these simpler conservation equations (Eq. 15), we held the volume fixed in the simplified model. Similar simplifying conservation equations for Na+ and K+ were used by Cressman et al., (2009) and Wei et al. (2014). All other equations remain the same as in the full model.
The Goldman–Hodgkin–Katz model
To show the robustness of our results, we reproduce some of the figures using Goldman–Hodgkin–Katz (GHK) formalism. The Hodgkin–Huxley model uses linear relationships for current and voltage (Ohm's law). Nevertheless, it is well known that the actual current flowing through ion-selective permeability channels in membranes is nonlinear, rectifying (Hille, 2001), and can be more accurately accounted for using the GHK equations (Hodgkin and Katz, 1949). The equations modeling the current due to the flow of a given ion across the membrane given by the Hodgkin–Huxley and GHK formalism, respectively, are as follows:
where E, R, T, z, and F, represent the Nernst potential, gas constant, absolute temperature, ion valence, and Faraday's constant, respectively. Ci and Co are the concentrations of a specific ion outside and inside the cell. P is the permeability of the membrane to the flow of a specific ion, which depends only on the types and numbers of ion channels present in the membrane. G is the ion conductance that measures the ability of the membrane to carry electrical current. The detailed calculations for the conversion between P and G are as follows:
G (in units of S/m2), is equal to σ/L, where σ (in units of S/m) is electrical conductivity and L (in units of m) is thickness of the cell membrane.
P (in units of m/s), is equal to D/L, where D (in units of m2/s) is diffusion coefficient.
The relation between ionic mobility μ (in units of m2/Vs) and electrical conductivity σ (in units of S/m) is σ = neμ, where n (in units of m−3) is the number density of monovalent ions and e is the electronic charge. For any substance, n can be expressed in terms of molar concentration C (in units of mol/m3) as follows: n = CNA, where NA is Avogadro's number. σ can also be expressed as follows: σ = CNAeμ = CFμ, where F is Faraday's constant.
Einstein's relation between D and μ is D = μRT/F.
Therefore, P/G = D/σ =
When we used GHK equations for voltage gated channels, we used the following substitutions in the Hodgkin–Huxley model:
retaining the Ohmic relationship for the leak current ICl = GClL(V − ECl), and where PNa and PK are the maximal permeability of the membrane for Na+ and K+ respectively. Thus, the critical difference between this version of the model and the full model described above is that in this version the GHK formalism is used for ion channel currents instead of the Hodgkin–Huxley formalism. The equations for ion concentrations and volume dynamics remain unchanged.
In the results reported, the bifurcation results in Figure 2E,F are produced using the simplified model; Figure 2H,I,J are produced by the GHK model; whereas all other figures are produced using the full model. In all three models, we examined the double-bifurcation properties of varying the parameters for potassium and oxygen to demonstrate unification.
Numerical methods
The bifurcation analysis of the model was performed using widely available software XPPAUT (Ermentrout, 2002) and MATLAB (MathWorks). We used the fourth-order Runge–Kutta method for integrating the differential equations. To facilitate the dissemination of these results, the MATLAB code required to reproduce the full model, shown in Figures 2A and 4A of this paper, and the XPPAUT code required to reproduce the bifurcation diagram of the simplified model in Figure 2E of this paper will be archived at the Penn State University Scholarsphere website.
Results
In our extension of the Hodgkin–Huxley equations, we apply principles of conservation of particles and charge, and account for the energy required to transport charges. First, we include the fact that oxygen diffuses from a distant reservoir (perfusion bath in vitro or microvasculature in vivo), and that pump rates depend on local tissue O2. Second, we fully account for ion concentration changes across the neuronal membrane and extracellular space, permitting the neuron's volume to change in response to osmotic pressure. In Figure 1, we contrast the traditional Hodgkin–Huxley model (Hodgkin and Huxley, 1952; Fig. 1, left) with our extended biophysical model (Fig. 1, right).
Model schematics. Left, Classic Hodgkin–Huxley model. Right, Extension of Hodgkin–Huxley model permitting volume of the neuron to change as transmembrane osmolarity changes (heavy red arrows), incorporating energy dependent pumps on both neuron and glia (red circles), and letting O2 and K+ diffuse with finite time to the cell (heavy gray arrows). Transporters KCC2 (green) and NKCC1 (yellow) regulate intracellular Cl−.
We consider two compartments for concentration of K+ and O2 (Cressman et al., 2009; Barreto and Cressman, 2011): the distant reservoir and the extracellular space (averaging 25 μm distant from capillaries in vivo, and typically 100 μm distant from bath when recording in vitro). At the cellular level, seizures will be characterized by paroxysmal high-frequency spiking where extracellular [K+]o does not exceed the physiological ceiling of 8–12 mm, whereas cellular SD will be characterized by a rapid depolarization sufficient to block spiking activity with [K+]o exceeding the physiological ceiling (Somjen, 2001). We discover a double-bifurcation (seizures between 8 and 12 mm K+, and SD >18 mm K+) at normal oxygen pressure (32 mg/L), accounting for many experiments at comparable ion concentration exposures (Traynelis and Dingledine, 1988; Anderson and Andrew, 2002; Zhou et al., 2010; Fig. 2A). By moderately lowering O2 pressure, we observe fusion of seizure and spreading depression dynamics (Fig. 2B). The two-parameter bifurcation diagram for [K+]o as a function of [O2]bath and [K+]bath (Fig. 2C,D), shows a parameter space representing separate regions for steady-state excitability as in the original Hodgkin–Huxley model, and extended regimes corresponding to seizures, tonic firing, spreading depression, and wave of death.
Bifurcations of the model. Minimum and maximum [K+]o as a function of [K+]bath, when [O2]bath = 32 mg/L (A), 24 mg/L (B, solid black), and 18 mg/L (B, dashed blue). A, Insets, Details around bifurcation points, with steady-state (left) and tonic spiking (double lines, right). C, Two-parameter bifurcations for full model shows minimum and maximum [K+]o as function of [O2]bath and [K+]bath. D, Projection of two-parameter bifurcation: SS, steady-state; SZ, seizure; SD, spreading depression; TF, tonic firing; WoD, wave of death. Analytic bifurcation analysis (XPPAUT) of membrane potential V as a function of [K+]bath when [O2]bath = 32 mg/L (E), and 18 mg/L (F) for the simplified model. Stable (red) and unstable (black) fixed points, and stable (green) and unstable (blue) limit cycles depicted. SN, Saddle node; HB, Hopf bifurcation. Numerical bifurcation solvers such as XPPAUT frequently do not resolve small amplitude stable oscillations accompanying the unstable (blue) lines in E and F, reflecting the slow rise and fall of the membrane potential baseline value as the K+ reversal potential oscillates during periodic spontaneous seizures and spreading depression. G, The minimum and maximum of extracellular space [K+]o changes as a function of [K+]bath, when [O2]bath = 32 mg/L in the full model. Increasing the [K+]bath to higher levels demonstrates a depolarization block steady-state >80 mm. Current–voltage curve of GHK (red: P/G = 3 × 10−6; blue: P/G = 2 × 10−6; where G = 1 mS/cm2) and HH (black) equations, with [Na+]i = 18 mm and [Na+]o = 144 mm (H), [K+]i = 140 mm, and [K+]o = 4 mm (I). J, Bifurcation diagram using GHK (red, blue) and HH (black) equations when [O2]bath = 32 mg/L.
Bifurcation analysis with a reduced model (fixed volume and enforced electroneutrality; see Materials and Methods) demonstrates the underlying mathematical structure of such dynamics (Fig. 2E,F). As [K+]bath increases, the cell transitions from steady-state (Fig. 2E, red) to unstable cycle representing seizure (Fig. 2E, blue) through a saddle-node bifurcation near [K+]bath = 8 mm. The seizure converts to tonic firing at higher [K+]bath (Fig. 2E, green), and transitions to SD (Fig. 2E, blue) near [K+]bath = 20 mm. Increasing [K+]bath further, the cell transitions to steady-state depolarization block through a Hopf bifurcation (the full model demonstrates a similar bifurcation >80 mm; Fig. 2G). Decreasing oxygen leads to a fusion of seizure and SD (Fig. 2F), as in the full model (Fig. 2B, blue).
In Figure 2H,I,J, we demonstrate the difference between the current–voltage relationships between Ohm's law (black lines) as in the original Hodgkin–Huxley equations, and the use of the nonlinear GHK current relationships. Note that for a substantial range of permeability values, the double-bifurcation in [K+]o is retained whether using the more complex GHK or the simpler Hodgkin–Huxley equations (Fig. 2A).
The full spectrum of unified dynamics is shown in Figure 3. In normal [O2]bath and [K+]bath, the neuron is at stable rest (Fig. 3G, black). A short pulse of current induces a single spike (Fig. 3G, red). Increasing Na+ leak current leads to periodic spikes (Fig. 3G, blue). This excitable steady-state, singlet spiking, and periodic spiking with increased leak are classic Hodgkin–Huxley dynamics. As seen experimentally (Traynelis and Dingledine, 1988), between 8 and 12 mm [K+]bath, seizure activity appears (Fig. 3B,E), accompanied by slow K+ oscillations (Fig. 4A).
Full range of neuronal behaviors produced by the model. A, Two-parameter bifurcation as function of [O2]bath and [K+]bath, where: SS, steady-state; SZ, seizure; SD, spreading depression; TF, tonic firing; WoD, wave of death. B, periodic seizures. C, Tonic firing in parameter region between seizures and SD at sufficiently high [O2]bath. D, SD at [K+]bath = 26 mm. E, Single seizure. F, SD at [K+]bath = 40 mm. G, At normal [O2]bath and [K+]bath, the neuron has a stable resting state (black), generates a single spike (red) from step current (15 ms, 5 μA/cm2). Black bar represents 5 ms. Increasing Na+ leak (INaL = 0.0557 mS/cm2), periodic singlet firing develops (blue). H, When [O2]bath is near zero, the wave of death is seen as the membrane potential approaches the Donnan equilibrium. I, Mixed dynamics of seizure and spreading depression when [O2]bath is reduced.
Volume regulation during seizures and SD. Membrane potentials reflect seizures (A) and SD (B, top traces). Extracellular potassium concentration was elevated (middle traces), and cells swell (bottom traces) during both pathological conditions. β is defined as the intra- to extracellular volume ratio. Cell volume is maximal during SD. C, The dynamics of membrane potential (top trace), [K+]o (middle trace), and cell volume (bottom trace) during hypoxia induced spreading depression. Red bar represents the period of hypoxia. D, High-bath potassium- (26 mm) induced periodic spreading depression. Red, blue, and green colors represent the reversal potential of sodium, potassium, and chloride, respectively. All intracellular and extracellular ion concentration changes are shown, along with local tissue oxygen concentration, [O2]o, volume fraction, β, and pump current, Ipump.
Near the [K+]o physiological ceiling (12–15 mm), seizures end (Barreto and Cressman, 2011). Above this ceiling tonic firing is predicted, the only parameter region of this model for which we are unaware of existing experimental exploration. Elevating [K+]bath further, the second bifurcation begins with slow, large amplitude [K+]o oscillations characteristic of periodic SD (Fig. 3D), leading to massive ion redistribution, oxygen depletion and cell swelling (Fig. 4B; full details in Fig. 4D). SD is experimentally provoked with bath potassium in this same range: 26 mm (Anderson and Andrew, 2002; Fig. 3D) and 40 mm (Zhou et al., 2010; Fig. 3F). Notice that the excitability at the start of SD in Figure 3D,F is consistent with SD single-cell recordings (Canals et al., 2005).
Lowering oxygen in the perfusate reduces Na+/K+ pump activity, converting seizure into periodic SD or mixed seizure-SD states (Fig. 3I; Czéh et al., 1993). As [O2]bath is further reduced, the oscillation range of [K+]o decreases because the Na+/K+ pump cannot recover ion gradients fully before subsequent SD episodes. When [O2]bath is very low, there is insufficient pump strength to support firing, and steady-state resting voltage returns.
Sufficient reduction in [O2]bath in severe hypoxia in vitro leads to single-episode HSD (Czéh et al., 1993). In the model, setting [O2]bath = 0 mg/L leads to HSD, where gradients can only be restored with oxygen restoration (Fig. 4C). During HSD, osmotic imbalance leads to cell swelling until reaching the minimal extracellular space constraint (Fig. 4C, bottom). If oxygen supply is not restored, the Na+/K+ ATPase, diffusion and glial buffering become inactive, leading to the wave of death (Fig. 3H; Zandt et al., 2011). Increasing oxygen availability by increasing [O2]bath decreases the parameter region where both pathological seizure and SD dynamics are observed (Fig. 5).
Bifurcation diagram with different maximal Na/K pump rate ρ. The minimum and maximum of extracellular space [K+]o changes as a function of [K+]bath when [O2]bath = 32 mg/L and ρ = 0.8 mm/s (black, identical to results in Fig. 2A), compared with the one when [O2]bath = 40 mg/L and ρ = 1.5 mm/s (red), simulating an increase in available oxygenation. Because the original sigmoid function is saturated at 32 mg/L, we raise the value of ρ to simulate the increased energy supply to the pumps by an increase in oxygen availability. Note that increasing oxygen supply decreases the range of pathological activity of seizures and spreading depression.
Discussion
We have found within a single model of the biophysics of neuronal membranes that we can account for a broad range of experimental observations, from spikes to seizures and spreading depression. We are particularly struck by the apparent unification possible between the dynamics of seizures and spreading depression. We define unification as the finding that seizures and spreading depression lay along a dynamical continuum of the neuronal membrane. That such unity of dynamics might occur was hinted at with the increasing discovery of monogenic mutations in humans that lead to both seizures and migraines (Rogawski, 2010), and the experimental description of mixture states of seizure and SD dynamics in the same cells in hypoxic (Czéh et al., 1993) or immature physiological conditions (Haglund and Schwartzkroin, 1990).
Migraines occasionally trigger seizures, seizures often initiate postictal headache similar to migraines, and certain antiepileptic drugs demonstrate efficacy in migraine prophylaxis (Rogawski, 2010). Recent data reveals that monogenic mutations (CACNA1A, ATP1A2, and SCN1A) can cause epilepsy, migraines, or both (Rogawski, 2010).
We did not intend to model unification but rather more accurately reflect known biophysical properties of neurons embedded within mammalian brain. The original Hodgkin–Huxley framework was that of a perpetual motion machine within an infinite bath containing constant ionic levels. We extended the Hodgkin-Huxley equations using conservation principles to account for the ionic fluxes and the energetics required to restore ionic gradients. Our results require an explicit extracellular space, within which the neuron dynamically adjusts volume, and is surrounded by a glial syncytium that both buffers extracellular K+, and has a share of the energy dependent electrogenic pumps. Our observation of unification turned out to be an emergent property of the fundamental conservation principles applied to such a biophysical model.
Emergence is the appearance of features in nonlinear dynamical systems (Lewes, 1874; Albowitz, 1939; Scott et al., 1999) that cannot be predicted from the properties of the elemental components. In our case, emergence of unification comes from adding in essential physiology, such as oxygen and potassium dynamics, as well as conservation principles in terms of mass, charge, and energy.
We would anticipate that alternative models might lead to similar complex dynamics (Prinz et al., 2004). Our conjecture is that as in other fields seeking unification of physical phenomena, there may be multiple routes to unified theories of neuronal dynamics. Nevertheless, ours is a minimalist model, and we observed that we could not remove oxygen or volume homeostasis within this biophysical model and maintain unification. Furthermore, the underlying mathematical bifurcation structure exhibits the properties of unification in the full model, and this further demonstrates that the properties of unification are not dependent on all of the specific details of the full model.
Our model opens up an entirely new way of modeling stimulation of the nervous system, by accounting for all of the charge species that compose the simulation currents applied. The consideration for charge conservation in excitable tissue simulation is a topic that has been studied in the field of heart tissue modeling. In Endresen et al. (2000), the principles of charge and mass conservation, and energy balance in the process of moving ions in and out of the cell, in heart tissue, were laid out clearly. They created an algebraic formulation for membrane voltage, leaving out anions. Hund et al. (2001) found that the algebraic formulation did not have a substantial advantage over the differential (Hodgkin–Huxley type) formulation, but they did determine for long-term stimulation of such tissue that the use of conservation of charge with explicit ionic charge carriers, which they termed “conservative stimulation,” was essential to prevent model drift (they modeled K+ as the charge carrier of stimulation). Last, Kneller et al. (2002) explored the use of Na+, K+, and Cl− as explicit charge carriers, demonstrating that there are substantial differences in the amount of physiological disruption during long-term stimulation. In their model, stimulation using K+ as the charge carrier was least disruptive.
Our model, like the original Hodgkin–Huxley equations, is that of a small membrane patch where propagation of a spike requires modification to account for spatial extent and traveling wave dynamics (Hodgkin and Huxley, 1952). Similarly, our model would require a network with spatial extent to account for the propagation of seizures, and modeling a spatial excitable medium with reaction–diffusion mechanisms to account for propagation of SD. Our representation is that of a neuron within native cellular architecture such as a brain slice and is highly consistent with slice physiology experiments, but for in vivo, brain representation would require extension to include an autoregulating vascular system and an explicit blood–brain barrier interposed between extracellular space and capillary blood.
In addition to conceptualizing spikes, seizures, and SD as states along a continuum, our unified framework offers a novel strategy for model-based stimulation control of neuronal activity (Schiff, 2012). In addition to electrical stimulation, the common clinical maneuvers of increasing oxygen supply and modifying osmotic pressure are clinical interventions that our model can help orchestrate through closed-loop feedback. Our model suggests optimal routes to direct brain states out of seizures or SD, minimizing the trajectories that would produce increased pathological activity injurious to the brain, especially in reduced oxygen availability. Our model further opens up the possibility of model-based prediction when approaching a bifurcation boundary, a topic gaining increasing attention in geophysics and ecology (Scheffer et al., 2012), as well as with human seizures (Kramer et al., 2012; Jirsa et al., 2014). Our bifurcation results on the simplified model suggests that there may be underlying normal forms for the unification dynamics we have here observed (Ros et al., 2014). We are now in a position to model the biophysics of the monogenic origins of familial hemiplegic migraine, where mutations in sodium and calcium channels and ATPase (Rogawski, 2010) can be incorporated into our model framework. Finally, our framework permits us to model the underlying stimulation charge carriers, and permits the energy cost of stimulating neurons to be optimized and damage from overstimulation avoided.
Our results demonstrate that unified frameworks for neuronal dynamics are feasible, can be achieved using existing biological structures and universal physical conservation principles, and may be of substantial importance in enabling our understanding of brain activity and in the control of pathological states.
Notes
Supplemental material for this article is available at https://scholarsphere.psu.edu/files/6395wc365. To facilitate the dissemination of these results, the MATLAB code required to reproduce the full model, shown in Figs. 2A and 4A of this paper, and the XPPAUT code required to reproduce the bifurcation diagram of the simplified model in Fig. 2E of this paper will be archived at the Penn State University Scholarsphere website. This material has not been peer reviewed.
Footnotes
This work was supported by NIH US-German Collaborative Research in Computational Neuroscience 1R01EB014641-01 (S.J.S.), and an Early Career Award (G.U.) and long-term visitor support (S.J.S.), from the Mathematical Biosciences Institute and the National Science Foundation Grant DMS 0931642. We thank M. Dahlem, G. B. Ermentrout, J. R. Cressman, E. Barreto, S. Van Wert, and B. Gluckman for helpful discussions.
The authors declare no competing financial interests.
- Correspondence should be addressed to Dr Steven J. Schiff, The Pennsylvania State University, W311 Millennium Science Complex, Pollock Road, University Park, PA 16802-2131. sschiff{at}psu.edu
