Simultaneous changes in ion concentrations, glutamate, and cell volume together with exchange of matter between cell network and vasculature are ubiquitous in numerous brain pathologies. A complete understanding of pathological conditions as well as normal brain function, therefore, hinges on elucidating the molecular and cellular pathways involved in these mostly interdependent variations. In this paper, we develop the first computational framework that combines the Hodgkin–Huxley type spiking dynamics, dynamic ion concentrations and glutamate homeostasis, neuronal and astroglial volume changes, and ion exchange with vasculature into a comprehensive model to elucidate the role of glutamate uptake in the dynamics of spreading depolarization (SD) —the electrophysiological event underlying numerous pathologies including migraine, ischemic stroke, aneurysmal subarachnoid hemorrhage, intracerebral hematoma, and trauma. We are particularly interested in investigating the role of glutamate in the duration and termination of SD caused by K+ perfusion and oxygen-glucose deprivation. Our results demonstrate that glutamate signaling plays a key role in the dynamics of SD, and that impaired glutamate uptake leads to recovery failure of neurons from SD. We confirm predictions from our model experimentally by showing that inhibiting astrocytic glutamate uptake using TFB-TBOA nearly quadruples the duration of SD in layers 2-3 of visual cortical slices from juvenile rats. The model equations are either derived purely from first physical principles of electroneutrality, osmosis, and conservation of particles or a combination of these principles and known physiological facts. Accordingly, we claim that our approach can be used as a future guide to investigate the role of glutamate, ion concentrations, and dynamics cell volume in other brain pathologies and normal brain function.
Spreading depolarization (SD) is a self-propagating wave characterized by a near-complete breakdown of transmembrane ion gradients in cells, sustained depolarization in individual neurons, and swelling of neuronal and glia cells [1,2, 3,4, 5]. It is now well accepted that SD is relevant to many neurological disorders. Several studies have shown that SD is the pathophysiological correlate of the symptoms of migraine aura [6,7, 8,9, 10,11], and occurs frequently in acutely injured brain caused, for example, by ischemic stroke, aneurysmal subarachnoid hemorrhage, and trauma [12,3, 13,14,4, 5]. Several clinical studies by COSBID group [15] and others suggest that SD mediates cortical lesion development and secondary brain damage in patients with acute brain injury, impairs clinical recovery, and triggers new deficits [12,13,3]. Furthermore, significant evidence indicates that SD and epileptic seizures might have some shared mechanisms [16,17,18,19]. The local processes during SD are understood as the interplay of neurons, astrocytes, and the vascular system. The neuron releases large amounts of K+ and glutamate into the extracellular space (ECS) together with significant drop in extracellular Ca2+, Na+, Cl−, and pH when it depolarizes. Consequently, SD is accompanied by significant extracellular K+ and glutamate accumulation, activation of N-methyl-D-aspartate (NMDA) receptors, a general loss of ion homeostasis, and cytotoxic edema [20,21,22,23,4, 5,24,25,26,27]. Excitotoxicity and SD are largely overlapping phenomena. Glutamate is of particular interest because of its role in excitotoxicity and its synchronous extracellular rise with the onset of SD [13,28,27]. Activation of NMDA receptors by glutamate triggers the release of further glutamate and K+ that will diffuse to neighboring cells thus causing the propagation and sustainment of SD. This hypothesis is backed by significant evidence of glutamate receptors antagonists inhibiting SD. Slices experiments showed that ischemic cells with NMDA and non-NMDA receptors blocked, did not exhibit the fatal form of SD [26,29]. An NMDA receptor antagonist, Ketamine was shown to inhibit SD in swine cortex [30] and reduced the number of SD incidences in patients [29,31]. Astrocytes and vasculature are other key players regulating many aspects of SD [4,14]. In addition to coordinating matter transport between vasculature and neurons and playing a major role in the observed metabolic and hemodynamics effects that are key to our understanding of numerous neurovascular diseases, astrocytes protect against SD initiation due to their high capacity for K+ and glutamate uptake [14]. Increasing the expression of astrocytic glutamate transporters reduces the infarct volumes following ischemia [32] and protects against the onset of ischemia-induced SD [33]. Astrocytic swelling together with changes in neuronal volume can exacerbate SD and may lead to severe brain damage [34,35]. In astrocytes, volume–activated anion channels may release large amounts of glutamate leading to excitotoxic damage [36]. The knockouts of aquaporin 4 channels that mediate astrocytic swelling [37], protect against ischemia [38]. To summarize, SD is accompanied by an array of immense changes from molecular to network level. A better understanding of SD and a spectrum of related pathologies, therefore, hinges on elucidating the pathways involved in these changes. However, existing techniques are too limited to investigate all these pathways. To overcome this void, we develop a comprehensive model that takes into account these key variables to quantify the role of glutamate dynamics in SD. We are particularly interested in SD caused by K+ perfusion and oxygen glucose–deprivation (OGD). The model equations are either derived purely from first physical principles of electroneutrality, osmosis, and conservation of particles, or by a phenomenological combination of these principles and known physiological facts. Our model is successful in explaining experimental results about the role of glutamate in SD. We confirm the predictions of our model by showing that astrocytic glutamate transporters blocker (2S, 3S) -3-[3-[4- (trifluoromethyl) benzoylamino]benzyloxy]aspartate (TFB-TBOA) significantly elongates the duration of SD in cortical slices from 15-24 days old rats. While our discussion is focussed on glutamate, the model can be used to explore the role of other key pathways and swelling in the dynamics of SD. Furthermore, the framework can be applied to investigate the role of ion concentrations, glutamate, and cellular volume dynamics in other pathological conditions and normal brain function. Numerous single neuron models for investigating SD have been developed. The phenomenon is rather generic and is found in models with great physiological details [39,40,41,42] as well as in simplified HH based descriptions of the neuron [43,44,45,46,47,39,48,49]. With the help of these models, thresholds for SD ignition and recovery can be assessed. In particular, it can be analyzed how energy and oxygen supply, morphological parameters, and blood pressure affect the course of SD, how SD can be prevented, and when it is non–recoverable [50,44,51,52,43,53]. Only few of these models deal with swelling. Some incorporate neuronal swelling alone [54,44,43,55], while one model [49] deals only with the astrocytic volume. Only two models Ref. [40,56] include neuronal and astrocytic swelling simultaneously. The models in Refs. [57] and [58] for regular neuronal spiking and epileptic seizures respectively deal with astrocytic glutamate uptake with no ion concentration dynamics or swelling. The model in Ref. [58] does not include glutamate release from neurons during spiking. To our knowledge, no neuronal model (SD or otherwise) deals simultaneously with ion concentrations, neuronal and glial volume changes, and glutamate dynamics. As discussed above, the extreme changes and interdependence of these pathways during SD warrants a comprehensive computational framework encompassing all these key pathways—the subject of this paper.
Rate equations for the membrane potential of the neuron, gating variables for K+ and Na+ channels, ion concentrations inside the neuron, glia, and ECS, and volumes of the neuron, glia cell, and ECS are based on our previous work [45,46,53,56,59,43,60,61,44,51,63,64]. These equations together with the modifications due to the inclusion of glutamate dynamics, and the morphology used in this model are described in S1 Text. Here we outline the details about modeling the glutamate homeostasis. Glutamate is a neurotransmitter that is released into the cleft of a synaptic connection when the presynaptic, i. e. signal–sending, neuron depolarizes. Glutamate binds to the NMDA and AMPA receptors of the postsynaptic neuron and can thereby initiate an action potential (AP). After binding to a receptor the transmitter is free again and can bind another time or diffuse into the ECS. Neurons and glia cells clear glutamate by taking it up from the cleft or from the ECS. For an overview of glutamate–related processes we refer the reader to reviews by Benarroch [66], and Kandel et al. [67] (see part III). Several components of the computational model presented in this section are adapted from various computational studies [68,69,70,58,71] and have been modified or extended for the application to SD.
We look at SD caused by perfusion of brain slices with high K+ and SD caused by OGD. The first case is modeled by increasing Kbath (K+ in the bath) from 4 mM to 15 mM at the start of simulation and stays elevated throughout the experiment. For the second situation all pump and glia functions slowly cease within 15 sec, remain interrupted for 40 sec, and are then slowly reactivated within 15 sec. Fig 2 shows the evolution of the membrane potential, Nernst potentials, ion concentrations, and the volumes. The OGD protocol is indicated by the orange bar in Fig 2b and 2d. The light parts at the beginning and end of the bar mark the smooth cessation and reactivation, respectively. In OGD, the uptake parameter v c → n m a x is also slowly set to zero and glutamate clearance is interrupted as well. The normal value of v c → n m a x is 0. 03 mM/msec. This is in the range of values suggested by Rusakov, and Slichenko and Tass in Refs. [70,58]. The main result of Fig 2 is that our addition of glutamate–related processes does not change the familiar course of events of SD. In both cases SD begins with a short burst of spikes driving the neuron into depolarization block. The burst causes huge changes in ion concentrations, most prominently a drop in extracellular Na+ and a huge rise in extracellular K+ (see curves for Nae and Ke in Fig 2c and 2d). The extracellular space is also rapidly flooded with high concentration of glutamate that prohibits the ion channels from closing. As a consequence of these ion changes, the membrane potential differences become very small (see Fig 2c and 2d) and the neuron enters a phase of sustained depolarization very quickly. The time at which depolarization begins is marked by a red triangle pointing upwards. Thereafter the neuron remains depolarized for about 60 sec before it suddenly repolarizes and ion concentration begin to recover. Note that this abrupt transition leads to a brief overshoot into hyperpolarization. The time of this potential drop is marked by a green triangle pointing downwards. The main difference between SD caused by K+ perfusion and OGD in the model is that in the perfusion experiments extracellular K+ (and other concentrations because of it) builds slowly as K+ diffuses from bath to ECS till it reaches a point where the cell starts spiking and enters SD. The cell comes out of SD spontaneously as the K+ clearance mechanism overcomes the release processes. If the simulation is allowed to run for longer time, the cell repeats this process till Kbath is reduced back to physiological values. Another difference between SD caused by K+ perfusion and OGD is that in the latter case there is no spiking in the beginning of SD. The insets of Fig 2c and 2d show how the volumes of the neuron, the glia cell, and the ECS change during these ion fluxes. The ECS shrinks dramatically while the neuron is depolarized. During OGD, ECS shrinkage gets much faster when ion pumps and glial functions are slowly reactivated. The reason is that glia swelling is blocked during OGD since there is no particle uptake. After reactivation of regulatory functions, the glia cell starts swelling because of K+ uptake and ECS shrinkage is accelerated immediately. The potential and ion dynamics for high K+ perfusion and OGD are very similar. However, as Fig 3 shows the glutamate dynamics clearly differ. Fig 3a and 3b show extracellular glutamate concentrations in the synaptic cleft and in the ECS. De– and repolarization times are indicated by triangles again. We note that depolarization increases the glutamate concentration in the cleft to more than 10 mM in both cases. The insets of Fig 3a and 3b show more details on finer scales. The sharp peak of the concentration in the cleft decays within 2 to 3 sec (upper insets). For OGD, this goes along with a rise in glutamate in the ECS due to diffusion. The two concentrations are equal after 3 sec. For K+ perfusion, Ge remains practically zero at all times. In the lower inset of Fig 3a, we see that after the peak Gc goes to values between 10 and 20 μM, while in the OGD simulation, the level remains 5 mM in the cleft and the ECS. Concentrations go back to zero with repolarization. The lower inset of Fig 3b shows that the neurotransmitter is cleared quickly when neuronal and glial uptake functions are reactivated. Recall that both are impaired during OGD. Concentrations are reduced to values in the range of 10 to 20 μM within a few seconds before repolarization takes them to zero. Fig 3c and 3d show the different uptake pathways. Main plots present uptake from all of the synaptic clefts of the neuron. The total height of the colored region is total uptake, the yellow and the orange portions are the specific contributions of the neuron and the glia cell. By buffered glutamate we mean molecules that have been taken up by the cells, but have not been recycled. Vesicle reproduction reduces N u p G, which is why the amounts of buffered glutamate do not strictly grow, but can also shrink (see Eq (30) ). In the OGD run, uptake dynamics only start with uptake reactivation by the end of OGD, so the different time signature is simply dictated by the OGD protocol. What these uptake plots show us is that glial uptake is dominant for clearance from both the clefts and the ECS (see insets). Moreover, more glutamate is cleared from the ECS than from the cleft in OGD–caused SD. The reason is that the amount of glutamate that diffuses into the ECS is larger than the amount that remains in the clefts. Also in perfusion–caused SD, large amounts of the neurotransmitter are taken up from the ECS as well. This type of uptake is fast enough to maintain a concentration Ge that is nearly zero at all times. We like to stress that Fig 3 should only be seen as an overview of the glutamate dynamics that our model provides. Details like the concentration plateau that follows the depolarization peak depend on the fine balance of release, uptake, and glutamate recycling. A different release function and a more detailed incorporation of the recycling process could lead to a different plateau. The most reliable aspects of our model are the amount of glutamate release during the depolarization burst, the diffusion process, and cellular glutamate uptake, since these parts of the model have been developed based on experimental data [58,70,67]. As pointed out above, increasing the expression of astrocytic glutamate transporters has been shown to reduce the infarct volumes following ischemia [32] and protects against the onset of ischemia-induced SD [33]. A recent study shows reduced rates of glutamate and K+ clearance by cortical astrocytes during neural activity and reduced density of excitatory amino acid transporters 1a (EAAT-1a) in cortical perisynaptic astrocytes in heterozygous FHM type 2-knockin mice [91]. By partial inhibition of glutamate transporters in wild-type mice, this study provides clear evidence that defective glutamate clearance can account for most of the facilitation of SD initiation in FHM type 2-knockin mice [91]. In vitro studies showed impaired glutamate uptake in hippocampal mixed astrocyte-neuron cultures from mice expressing FHM type 2-causing α2 Na+/K+ ATPase. Induction of SD in these animals resulted in reduced recovery [92]. In the following, we provide a complete understanding of the glutamate uptake processes and its role in SD by reducing v c → n m a x (rates for other uptake pathways are implied). The simulations are shown in Fig 4 and they are consistent with our experimental observations (see below). Fig 4a and 4b show the membrane potential dynamics for maximal uptake rates of 20% and 18% of the v c → n m a x–value used in Figs 2 and 3. The repolarization time for normal uptake is about 143 sec as indicated by the vertical dashed line. With the lower uptake rates, recovery is delayed by about 65 and 105 sec, respectively. Before repolarizing, the membrane potential oscillates with a low amplitude. We note that the potentials for the first 143 sec into SD are the same in Figs 2a and 4a and 4b. This is also true for the ion concentrations, and we conclude that for the first 143 sec the ion fluxes are almost not affected by the smaller vmax–values. This implies that the contribution of the cotransport currents I i o n c o from Eq. (32S) – (34S) is negligible. The only difference for the first 143 sec between normal and impaired uptake lies in the much higher glutamate concentrations Gc (see Fig 4c vs the lower inset of Fig 3a). In other words, too much glutamate in the cleft prevents recovery. The neuron only repolarizes when the concentration is small enough. An uptake rate of 20% achieves this sooner than the one at 18%. In Fig 4d, we look at this delaying effect systematically and compare SD durations for a range of uptake rates between 100% (normal value) and 16%. The effect becomes noticeable at 35% and lower with a maximal duration of almost 500 sec for 16% uptake. Below 16% recovery fails. For uptake rates between 100% and 50% the duration of SD is nearly constant. In this range of uptake rates the evolution of Gc is comparable to Fig 3a. The initial jump in Gc is quickly reversed and recovery is unaffected by glutamate. Only for uptake rates of 35% and less do we see a slow down–regulation of Gc as in Fig 4c. We expect critical uptake rates to depend very sensitively on all processes of the glutamate cycle. The numerical values in Fig 4d are likely to differ between models and should not be assumed for a real system. However, we believe that the basic correlation between the duration of SD and glutamate regulation is a real effect that explains the delayed recovery from perfusion-induced SD in tissues that are exposed to TBOA in our experiments (see below). A recent study showed that 0. 5 and 1 mM TBOA prolonged SD by 148% and 426% respectively [28]. A respective increase of 167% and 374% in glutamate concentration was observed in the same experiments. We can understand delayed recovery and recovery failure in SD from a phase space perspective. Let us briefly review the general method for the case of intact glutamate uptake (see Figs 2 and 3). The role of glutamate clearance is addressed in a second step. One can determine the timescales of the different processes in the model by a dimensionality analysis, and it turns out that Δ N b a t h K and Δ N g l i a K are the slowest variables (see Ref. [46] for a complete timescale analysis of a very similar model). Please note that in our model presentation it is not possible to read off the timescales from the coupling parameters. They have different units and consequently different orders of magnitude. A dimensionless presentation can be obtained by appropriate rescaling of the dynamical variables, but this is not our focus here. The timescales of vascular coupling and glial buffering have been derived rigorously in the above mentioned study [46]. These processes are the slowest ones in the system and they occur on similar timescales. This can be used to apply a so–called slow–fast analysis which allows us to derive the threshold condition for repolarization. Let us combine the two quantities to define a single slow variable: Δ N K ≔ Δ N g l i a K + Δ N b a t h K (31) In a slow–fast analysis, we ask which states of the system are possible when a certain value of the slowest variable is given. To find these states we treat this variable as a parameter, which formally defines a subsystem of only the fast dynamics. The processes of the fast subsystem are membrane dynamics and ion fluxes across the neural membrane. We refer to them as transmembrane dynamics. For the fast subsystem we also treat glutamate as a parameter. This is not an approximation, because glutamate is fast. However, it is the only way to systematically study how elevated glutamate levels affect the other processes. For now we set the level extremely low (0. 0001 mM), which practically means we ignore glutamate. This assumption is good enough at the de– and repolarization points. Depending on our choice of ΔNK, the fast subsystem has a certain number of stable and unstable fixed points. With the help of the software tool AUTO [83], these fixed points can be found and followed, while ΔNK is varied within a certain range. Every stable fixed point of the ‘fast subsystem’ is stable in the full system except for dynamics of ΔNK. When the timescale separation is sufficiently large the fast variables equilibrate to one of these stable fixed points, and this way the dynamics of the whole system will be guided by the fixed point structure of the subsystem. Formally, this concept is known as a ‘quasi–steady–state reduction’ [93]. In Fig 5, we see how this works. The left panel (Fig 5a) shows the evolution of ΔNK over 200 sec. It is indeed slower than the other variables (see Fig 2) with the exception of Cli/e—a detail that has been addressed previously and can be neglected here [56]. Note that Δ N g l i a K and Δ N b a t h K are two independent dynamical variables, they affect volumes differently. So knowing ΔNK alone is not enough information and we need to specify Δ N g l i a K too. That is why, Δ N g l i a K is also included in Fig 5a. We have marked the crucial values at the de– and repolarization points by a pink and a turquoise X respectively. We see that the neuron depolarizes at the minimal value and repolarizes at the maximal value of ΔNK. This means that SD starts when the K+ content of the neuron and its ECS is elevated, and that recovery is on the other hand only possible when enough K+ is taken away (recall the sign of the difference terms in Eq. (22S) ). The relation between ΔNK and these events becomes more clear in Fig 5b, which shows the location of the fixed points of the fast subsystem in the (ΔNK, V) –plane and how the trajectory of the full system is guided by them. Near the depolarization point, we use the fixed point curves for Δ N g l i a K ≈ 129 fmol and near the repolarization point we use Δ N g l i a K ≈ 226 fmol. Note that we do not show the complete fixed point curves. The entire curves are both z–shaped and overlap strongly. Hence for clarity, we only show two disconnected portions of these curves. The trajectory of the full system is a closed loop in the phase space and the transition points can now be understood through the fixed point structure. Let us focus on the repolarization process. Depolarization can be explained analogously, and we refer the reader to a previous study for more details on both transitions [46]. As the trajectory approaches the repolarization point, it is closely guided by the stable depolarized fixed point curve. The potassium content decreases until the curve becomes unstable. Stability changes of fixed points generically occur in Hopf bifurcations or limit point bifurcations. The AUTO software provides this information and shows that the above change of stability is due to a Hopf bifurcation. For completeness, we mention that around a Hopf bifurcation there are always stable or unstable limit cycles. These limit cycles imply oscillation, which we often see shortly before repolarization. Typically, stable limit cycles in SD models only exist in a narrow range of ΔNK values around the Hopf bifurcation before they change stability in a limit point bifurcation of limit cycles (see Figs 2 and 3 in Ref. [46] where K ˜ e = - Δ N K / ω e). When there is no stable upper fixed point or limit cycle the trajectory drops back onto the polarized stable fixed point branch. Since the transition from depolarization to repolarization happens at or very close to the Hopf bifurcation, the two can be related to each other. We remark that the situation will be different when the upper and lower fixed point branches do not overlap (see below). The tracking of limit cycles in our model is numerically very involved, because of the many different timescales. The continuation of limit cycles is hence beyond the scope of this study, but a complete analysis for a similar model has been performed before [46]. Note that our bifurcation diagrams only contain the bifurcations that change the stability of the fully stable fixed points. Subsequent bifurcations that change the degree of instability are not relevant to our analysis here and are omitted. From the phase space perspective, we understand that recovery from SD relies on the existence of a stable repolarized state in the fast subsystem. In Fig 5b, such a state is available and also seems to exist with higher glutamate concentrations as in Fig 4a and 4b. We now increase Gc and see how the fixed point curve changes. Unlike the treatment of ΔNK as a parameter, this is not an approximation in the sense of a slow–fast analysis. Glutamate–related processes happen on fast timescales (see Fig 3a and 3b). We rather treat Gc as a parameter to study the effect of glutamate in some extreme scenarios. For example, how does high glutamate concentration in the cleft—that may occur with impaired clearance—affect the dynamics of the system? By fixing Gc, we can obtain qualitative insights and answers to such questions. To discuss recovery, we look at the fixed point curves for high Δ N g l i a K. Fig 6a shows these curves for three different values of Gc. The two lower values 0. 02 mM and 0. 05 mM are chosen to show that the fixed point structure is very sensitive to Gc. The highest value 0. 3 mM is included because it is near the glutamate levels that prevent early recovery in Fig 4c (see values at the dashed vertical line). The fixed point curves have a stable depolarized (upper) branch that ends at a maximal value of ΔNK and a stable polarized (lower) branch that begins at a minimal value of ΔNK. These points are defined by Hopf bifurcations. The red line contains all the lower branch Hopf bifurcations for Gc–values between 0 mM and 0. 35 mM, the green line contains the upper branch Hopf bifurcations. The polarized branch shifts towards higher ΔNK–values as Gc increases. For the depolarized branch, this shift is much smaller and the green bifurcation line is consequently shorter in the (ΔNK, V) –plane. For low Gc, the upper and lower fixed point branch overlap, and the sharp transition from a depolarized state to a polarized state is possible. For higher values there is a gap between the branches and instead of repolarization the system goes into persistent low amplitude oscillations after the upper branch becomes unstable. We have not included the corresponding limit cycles in Fig 6a, however, the time series in Fig 4 shows such oscillations. In Fig 4, the oscillations are unstable, causing the cell to transition to the polarized state. In case of recovery-failure, these small amplitude oscillations persist for the duration of simulations. Since we understand repolarization as the transition between overlapping fixed point branches, the interesting question is, at what glutamate level this overlap disappears. Above this level, recovery is no longer possible. In Fig 6b, the end of the upper stable fixed point branch and the beginning of the lower branch are shown in the (Gc, ΔNK) –plane. As long as the depolarized branch ends after the polarized branch begins, recovery is possible. The critical Gc–value is hence at the intersection of the two lines, which occurs near 0. 035 mM. If Gc was not a dynamical variable, but a system parameter, this value would be the threshold for recovery failure in SD. In our simulation of the full system, however, we cannot separate the dynamics of V and Gc. As V decreases, glutamate release slows down and consequently Gc decreases too. On the other hand, an increasing glutamate level depolarizes the neuron, and accordingly the two effects amplify each other. This leads to the glutamate drop from about 0. 1 mM to nearly zero at the repolarization point (see Fig 4c). Because of this fast interplay of V and Gc, the critical value derived in Fig 6b is not an obvious threshold in the full system, but only gives us a rough idea about glutamate levels near the repolarization point. The fixed point curves in Fig 6a do not approximate the dynamics of the whole system. Moreover, even the fast scale glutamate dynamics of the whole system are only a rough approximation of a real system. Nevertheless, we have learned that critical glutamate levels exist, beyond which the neuron will not repolarize. This observation is based on a model with a parametrical glutamate concentration in the cleft and glutamate coupling through NMDA and AMPA receptors. These parts of our model are based on a more accurate biophysical description than the release and uptake mechanisms of glutamate. Accordingly, we are confident that the effect we have found in Fig 6 is relevant: too much glutamate prevents recovery from SD. In the next section we will provide more insights into this effect. To understand how glutamate interferes with recovery, it is helpful to look at the membrane model, because the first step towards recovery is a change in the membrane state. A given set of ion concentrations determines the reversal potentials Eion and the pump current Ip. These quantities define the membrane model belonging to this ion configuration. We are now interested in the membrane models of the ion configurations around the repolarization point. That means that ion concentrations are now model parameters and we vary them such that we obtain the ion configurations on the depolarized fixed point branch near the repolarization point in Fig 5b. This parameter variation is naturally parametrized by ΔNK. The result of this continuation is shown in Fig 7a. There are two fixed points in the membrane model. One is depolarized and coincides with the fixed point of the whole transmembrane model. The other fixed point is in fact hyperpolarized. That is, it is more strongly polarized than EK. At the repolarization point, the depolarized state becomes unstable, while the hyperpolarized state continues to exist. The membrane potential drops very close to this point before ion concentrations and Nernst potentials re–adjust and bring the system to a slightly higher potential. The existence of a stable hyperpolarized membrane state is what initially drags the membrane potential down and is crucial to the neuron’s recovery. Let us have a closer look at the membrane states and compare the two shortly before the repolarization point. In Fig 7, we have marked the two stable membrane states that exist for ΔNK = 285 fmol. Some quantities that characterize the membrane state are listed in Table D in S1 Text. The conductances in the depolarized state are dominated by the gated channels which can be seen from g K ≫ g K l (and similar for Na+). For the hyperpolarized state, the opposite is true and we have g K ≈ g K l instead. Accordingly, the hyperpolarized fixed point condition can be approximated as g K l (V - E K) + g N a l (V - E N a) + g C l l (V - E C l) = - I p. (32) We can draw the following conclusion from this relation. The pump current Ip is nearly maximal, because we assume a high concentration of extracellular K+. Since the leak conductances are rather small, we conclude that the potential difference terms (V − Eion) must be sufficiently negative in the polarized fixed point. In particular, the large pump current forces the membrane potential below the K+ Nernst potential, which is what we call hyperpolarization. In summary, hyperpolarization is the result of a large pump current and very small conductances. If the leak conductances were not as small, the potential difference terms in Eq (32) would be less negative and depolarization would become weaker. In turn less depolarization will violate the approximation g K ≈ g K l and instead gK would be larger than g K l, which leads to even less depolarization. So we have an understanding of repolarization that demonstrates how important it is that the conductances of the neuron collapse strongly enough. The above consideration only took into account leak conductances and the normal gated ion channels. In addition to that, an elevated glutamate concentration in the cleft implies increased conductances of the NMDA and AMPA receptors as well and the same argument as above holds—the increased conductance implies a less polarized lower fixed point. Again, we increase Gc as a parameter and study its effect. In Fig 7b, we track the hyperpolarized state as we follow the end of the depolarized fixed point branch. For values from 0 mM to 0. 0336 mM, there is a stable polarized membrane state when the upper fixed point branch ends. It becomes less polarized with increasing Gc which is consistent with the above reasoning on the importance of collapsing conductances. For very high Gc values, a hyperpolarized membrane state no longer exists and repolarization becomes impossible. The critical value is consistent with the value we derived from the transmembrane model in Fig 6b. In summary, the membrane model teaches us two things. First, a breakdown of neural conductances is needed for the neuron to repolarize sufficiently. Second, over-stimulated synapses can prevent this process and lead to recovery failure. To confirm model predictions, we recorded SD episodes in layers 2-3 of visual cortex slices from 15-24 days old, male wild type Sprague Dawley rats, both at individual neuron and network levels. To initiate SD, we replaced control ACSF with high KCl (26mM) ACSF (see Methods section). EC and IC electrodes were placed about 500μm apart. In high K+ ACSF, SD typically occurred 30-45s after high KCl application. SD in the single cells in control groups started with a rapid depolarization followed by several spikes and a slow return to the resting membrane potential. The majority of cell spiking occurred before maximum depolarization was reached. SD in the EC recording was typically noticeable a few seconds later. SD in single cells typically lasted for 30-180 seconds. To evaluate the effect of impaired astrocytic and higher glutamate concentration in the ECS on SD, slices were incubated with 50nM astrocytic glutamate transporter blocker TFB-TBOA for 20 minutes. Fig 8 shows a summary of these experiments. Example traces representing the membrane potential of individual pyramidal cell (bottom traces) and EC recording at the network level (top traces) during SD from control (n = 9) and TFB-TBOA-treated (n = 9) slices are shown in Fig 8a–8c. Unlike control slices where several spikes were observed before the cell entered a depolarization block, TFB-TBOA completely blocked action potentials at the single cell level. The resting membrane potential was not affected by the application of TFB-TBOA. The mean duration of SD averaged over many events at the network and single cell level is shown in Fig 8d and 8e respectively. SD duration was defined and measured as the time from the initiation of rapid depolarization in individual neurons and network to the time when the membrane potential repolarized to its pre-SD potential value (dashed lines in Fig 8a–8c). In line with model predictions, blocking astrocytic glutamate transporters almost quadruples the duration of SD. It is also worth noticing that single neurons exhibit a stronger depolarization as compared to neurons in control slices. Furthermore, some neurons in control slices exhibit a behavior that is termed as mixed seizure and SD state [94,44,43] (see for example lower trace for control cell in Fig 8c), while strong SD was observed in almost all single neurons in the slices treated with TFB-TBOA. We would like to remark that, although not exactly the same, the average SD duration given by the model is comparable to experimentally observed values. A 100% glutamate uptake in the model gives over a minute long SD as compared to ∼ 2. 5 minutes long SD in control slices. This difference could be due to the higher K+ in the perfusion solution used in the experiment as compared to the model. It could also be due to our overestimation of the glutamate uptake (or underestimation of glutamate release) as lower uptake leads to longer SDs. Application of TFB-TBOA increases the duration of SD to almost 10 minutes. Decreasing the glutamate uptake by transporters to about 17. 5% of the control value in the model leads to a comparable (about 4 times) increase in the duration of SD as compared to control simulation. Reducing the uptake to 16% of the control value will lead to SD duration comparable to the observed values. We would like to point out that TFB-TBOA targets only astrocytic glutamate transporters, while the reduction in glutamate uptake in the model applies to both neuron and astrocyte. So, the application of TFB-TBOA does not necessarily mean the complete inhibition of glutamate uptake. Given that glia cell has about eight times more binding sites for glutamate and has about half the surface area available for uptake (see Diffusion and Glutamate Uptake section) as compared to neuron, there would still be roughly one-fourth glutamate transporters intact (ignoring other complications due to morphology, ion concentration dynamics etc.) even in the presence of TFB-TBOA. We also observed time-dependent changes in the neuronal properties due to TFB-TBOA. Action potentials were evoked by current pulses and studied under the current-clamp conditions (n = 9). TFB-TBOA (50nM) increased the cell membrane resistance from 150MΩ to 260MΩ (p<0. 05) after 10 minutes. Current pulses evoked fewer APs (7 versus 11 on average) in slices pretreated with 50nM TFB-TBOA (p<0. 05) (0-100pA current injection). After 20min treatment with 50nM TFB-TBOA, the AP threshold increased from -45mV to -32mV (p<0. 05), while the amplitude of AP decreased from 80mV to 30mV (p<0. 001). AP amplitude was measured as the voltage difference between the threshold and the peak value of the membrane potential during the AP. We believe that these changes in the neuronal properties could be due to the reduction in Na+ influx due to TFB-TBOA. In line with this argument, Bozzo et al. [95] demonstrated the inhibitory effects of TFB-TBOA on astrocytic Na+ responses to glutamate. They also claimed that TFB-TBOA has no effect on the membrane properties of cultured cortical neurons recorded in the whole-cell patch clam. Recently, Hosseini-Zare et al [90] on the other hand, claimed that fast voltage-gated Na+ currents are reduced by TFB-TBOA. Whether the reduction in Na+ currents is caused by a direct interaction of TFB-TBOA with fast voltage-gated Na+ channels, through inhibition of Na+ cotransport (notice that three Na+ are cotransported with one glutamate molecule), or through some other mechanism is not entirely clear. While our model includes the inhibition of Na+ cotransport through glutamate transporters, incorporating the effect of TFB-TBOA on fast voltage-gated Na+ currents, if proven unequivocally, in the model is beyond the scope of this study and will be investigated in the future.
Significant experimental and clinical data suggest that SD is involved in numerous brain pathologies including migraine, stroke, subarachnoid hemorrhage, and traumatic brain injury [13,12,3, 4]. The initiation, propagation, sustainment, and termination of SD involve immense changes in many molecular and cellular pathways shaping the interplay between neurons, extracellular space, glial cells, and vasculature. What complicates things further is that most of these modifications are dependent on each other, and some may have biphasic role [96]. For example, over-activation of NMDA receptors in the early phase of stroke is detrimental, but in delayed phase, they might mediate neuroprotection through neuroplasticity [97]. That is probably why all trials testing NMDA antagonists for stroke treatment have failed [96,98]. In clinical trials, SD episodes were discontinued in two patients treated with ketamine on one hand [29], while on the other hand, a cluster of SD occurred in another patient despite the presence of ketamine [99]. Thus a complete understanding of SD and finding clinically useful therapeutic interventions for the related pathologies hinge on elucidating this wide array of changes. However, the current experimental and clinical tools are too limited to simultaneous investigate all these changes, which necessitates physiologically relevant detailed computational models. In this paper, we developed a comprehensive model that incorporates many key elements involved in the dynamics of SD including: neuronal membrane potential dynamics, ion concentration dynamics in neuron, extracellular space, and glial cell, ion exchange with vasculature, swelling of neuron and glia, and detailed formalism of glutamate release and uptake processes. Although, we explore the effect of glutamate uptake and extracellular levels on the initiation, sustainment, and termination of SD, our model allows us to investigate the role of all these factors in the dynamics of SD simultaneously or one by one. Our results show that glutamate signaling plays a key role in the dynamics of SD since impaired glutamate uptake prolongs the duration of SD and leads to significant neuronal and glial swelling. Reducing glutamate uptake by transporters below 16% of the control value leads to the failure of cell’s recovery from SD. We verified this prediction experimentally by showing that SD in layers 2-3 of visual cortex from 15-24 days old rats are significantly prolonged by inhibiting glial glutamate uptake using TFB-TBOA. Our computational results are also consistent with a recent study, which showed that 0. 5 and 1 mM TBOA prolonged SD by 148% and 426% respectively [28]. A respective increase of 167% and 374% in glutamate concentration was observed in the same experiments. Our result is also in line with conclusions from in vivo and in vitro studies in Ref. [80], where the elimination of glial glutamate transporters were shown to lead to tonic increase in extracellular glutamate, resulting in widespread swelling and neuronal degeneration. Furthermore, increasing the expression of glial glutamate transporter EAAT2 through application of β-lactam antibodies significantly reduced extracellular glutamate in animal studies and protected against ischemic injury and neurodegeneration [33,32]. We would like to remark that the pattern of neuronal and glial swelling, and the dynamics of various ion concentrations observed in our model explains several experimental observations. We skip such details here and refer the interested reader to our recent work for a detailed discussion about these observations and the role of swelling in ischemic injury [56,53]. Furthermore, the dysfunction of glial glutamate transporters is also implicated in other acute and chronic neurological disorders including amyotrophic lateral sclerosis [100], brain tumors [101], epilepsy [102,103], Alzheimer’s disease [104], and motor discoordination [105]. Our approach can be adopted to quantify the role of different pathways involved in glutamate dynamics in these conditions. As mentioned above, our model skips several factors that could be key for glutamate homeostasis. For example, glutamate flux through astrocytic glutamate transporters reverses direction in the presence of high EC K+ or high IC glutamate concentration [106,107]. Thus, during SD where K+ and glutamate are both high, reversal of glutamate transport in astrocytes could become an additional source of EC glutamate build-up. The blockade of GLT transporters during SD in our experiments may have prevented glutamate transport reversal and caused glutamate to be trapped inside the glia, decreasing EC glutamate levels. Under different experimental conditions, however, it was observed that a different glutamate transporter blocker (D, L-threo-beta-hydroxyaspartate (THA) ) increased the amount of depolarization and duration of SD in the presence of high K+ [108]. The changes in the single cell properties observed in our experiments are also ignored in our model. Furthermore, our study is concerned mainly with the behavior of a single neuron during SD. Investigating the effect of changes in glutamate homeostasis on the spatial spread of SD will require a network model. Incorporating these key factors in the model is beyond the scope of the current manuscript and is the subject of our future studies. We would also like to point out that there are different types of SDs with different features and probably different mechanisms for induction and propagation [109,3]. Our study focuses on SD caused by K+ perfusion and OGD. SD due to OGD is particularly relevant for stroke. Whether glutamate is necessary for propagation of all kinds of SDs is still debated. Interestingly, very little glutamate diffuses out of the cleft when uptake is not impaired in our model, consistent with a glutamate-independent propagation of SD in our modeling conditions. To conclude, by combining an established semi–phenomenological neuron–glia description and first physical principles in a consistent way, we have developed a physiologically relevant, comprehensive model that incorporates many key components involved in the dynamics of SD. This new mathematical framework describes many aspects of neuronal membrane, ion concentration dynamics, cell swelling, and glutamate dynamics during SD accurately and provides deep insights into the mechanisms through which glutamate interferes with neuronal recovery from SD. We present strong experimental evidence in support of our study, and emphasize that most of our explanations come from general physical principles and biophysical reasoning. The theory is general and the components included are key to both normal and pathological brain function. Accordingly, we claim that our approach can be used as a future guide to investigate the role of ion concentrations, ion exchange with glia and blood vessels, cell swelling, and glutamate dynamics in other brain pathologies and normal brain function. For example, to investigate glutamate homeostasis in regular neuronal firing, one would assume that only a few synapses are involved in the release of glutamate unlike SD where we assume all 10,000 synapses releasing glutamate. Similarly, to investigate these variables in seizures induced by high K+, higher K+ concentration in the bath, Kbath (typically 8mM) should be used. Furthermore, the additional currents involved in the specific neuronal behavior in question should be included in the membrane potential ion concentrations dynamics.
Pathological conditions such as seizure, migraine, traumatic brain injury, and stroke are associated with extreme changes in ion concentrations and glutamate, cell swelling, and heavy exchange of matter between neurons, glia, and vasculature. However, current experimental tools are capable of measuring only a few of these variables, which necessitates the development of biophysically relevant models. This study provides a comprehensive computational framework derived from first physical principles and physiological facts that enables us to investigate a wide range of key variables in SD caused by K+ perfusion and oxygen-glucose deprivation. While we use SD as an example, our approach can be extended to other pathological conditions and normal brain function.
medicine and health sciences
neurochemistry
action potentials
nervous system
depolarization
membrane potential
electrophysiology
neuroscience
physiological processes
homeostasis
neurotransmitters
neuronal dendrites
animal cells
glutamate
biochemistry
cellular neuroscience
cell biology
anatomy
synapses
physiology
neurons
biology and life sciences
cellular types
neurophysiology
The role of glutamate in neuronal ion homeostasis: A case study of spreading depolarization