Ion dynamics at the energy-deprived tripartite synapse

The anatomical and functional organization of neurons and astrocytes at ‘tripartite synapses’ is essential for reliable neurotransmission, which critically depends on ATP. In low energy conditions, synaptic transmission fails, accompanied by a breakdown of ion gradients, changes in membrane potentials and cell swelling. The resulting cellular damage and cell death are causal to the often devastating consequences of an ischemic stroke. The severity of ischemic damage depends on the age and the brain region in which a stroke occurs, but the reasons for this differential vulnerability are far from understood. In the present study, we address this question by developing a comprehensive biophysical model of a glutamatergic synapse to identify key determinants of synaptic failure during energy deprivation. Our model is based on fundamental biophysical principles, includes dynamics of the most relevant ions, i.e., Na+, K+, Ca2+, Cl− and glutamate, and is calibrated with experimental data. It confirms the critical role of the Na+/K+-ATPase in maintaining ion gradients, membrane potentials and cell volumes. Our simulations demonstrate that the system exhibits two stable states, one physiological and one pathological. During energy deprivation, the physiological state may disappear, forcing a transit to the pathological state, which can be reverted when blocking voltage-gated Na+ and K+ channels. Our model predicts that the transition to the pathological state is favoured if the extracellular space fraction is small. A reduction in the extracellular space volume fraction, as, e.g. observed with ageing, will thus promote the brain’s susceptibility to ischemic damage. Our work provides new insights into the brain’s ability to recover from energy deprivation, with translational relevance for diagnosis and treatment of ischemic strokes.

Information transfer at synapses [1] critically depends on the cellular availability of unexplored, whereas experimental data shows that the duration of ED is a crucial A model for ion homeostasis at the tripartite synapse 96 We introduce a novel biophysical model, extending our previous work on a single 97 cell [11] by incorporating more biological detail. The model in [11] contained a single 98 neuronal compartment in an infinite extracellular space, whose volume is regulated by 99 changes in osmolarity. Here, we model a neuronal and astrocytic soma, a presynaptic 100 terminal, the perisynaptic astrocyte process, the synaptic cleft and a global extracellular 101 space (ECS), illustrated in Fig. 1. The volumes of the somatic compartments and ECS 102 can change; the cleft, the presynaptic terminal and the perisynaptic astrocyte process 103 have fixed volumes. In each of these compartments, we describe the energy-dependent 104 dynamics of the ions Na + , K + , Cl − , Ca 2+ and glutamate, including vesicle recycling, 105 by ordinary differential equations (ODEs). We further assume that Ca 2+ and glutamate 106 are confined to the presynaptic terminal, the perisynaptic astrocyte process and cleft. 107 The concentrations of Na + , K + and Cl − are assumed to be the same in the synaptic as 108 well as in the somatic compartments. This allows us to work with both small molar 109 changes of Ca 2+ and glutamate and large Na + , K + and Cl − gradients across the somatic 110 membranes, within the same compartmental framework. Compartments, ion channels and transporters included in the modelling of the glutamatergic synapse. Shown are the three main components representing a presynaptic neuron, an astrocyte and the extracellular space (ECS). Each of these compartments also contains a synaptic compartment as indicated by the different shading and the additional box (presynaptic terminal, perisynaptic process and synaptic cleft, respectively). The largest ATP consumption in the presynaptic neuron and the astrocyte is by the Na + /K + -ATPase (NKA). At the presynaptic terminal, ATP is also needed to energize glutamate uptake into vesicles. The key transporters at the cleft are the Na + /Ca 2+ -exchanger (NCX) and the Excitatory Acid Amino Transporter (EAAT). NKCC1: Na + -K + -Cl − -cotransporter. KCC: K + -Cl − -cotransporter. Kir4.1: K + inward rectifier channel 4.1.

111
Model equations 112 For an ion X in compartment i, we describe the dynamics of the number of moles N i concentration [X] i of ion X in compartment i by [X] i = N i X /W i . A summary of the 117 notation used is shown in Table 1. Thus we obtain the following system of differential 118 equations that describe the dynamics of N i X and W i where F is Faraday's constant, z X is the valence of ion X, λ i is the osmotic flux rate for 120 compartment i, and q denotes the Hodgkin-Huxley gating variables. To model 121 glutamate dynamics in the cleft, we combine existing models of vesicle recycling from 122 Tsodyks and Markram [42] and Walter et al. [43], as illustrated in Fig. 2. Non-releasable 123 glutamate pool (N ) is made readily-releasable through four Ca 2+ -dependent 124 intermediate steps (R, R 1 , R 2 , R 3 ). After release into the cleft from a fused vesicle, (F ), 125 glutamate is removed by presynaptic and astroglial EAATs and enters the 'inactive (I)' 126 state. As a final step, glutamate is again stored in a vesicle depot (D) and enters the 127 non-releasable pool (N ) again. The rate equations of these states are given by where Y and X span over {N, R, R 1 , R 2 , R 3 , I, D, F }. The rate constants ν X depend on 129 Ca 2+ concentrations in the neuronal synaptic compartment. Eqs. (1) and (2) describe 130 ion and volume dynamics in the neuronal and astrocyte compartments. To obtain 131 extracellular dynamics, we use three conservation laws, i.e., for the preservation of charge, volume, and mass, respectively. Our model is essentially 133 described by Equations (1), (2) and (3). Here, the constants C W and C X are the total 134 volume and total molar amount of ions present in the system, respectively. We present 135 further details in the section Materials and Methods. 136 Simulations and model calibration 137 In the following sections, we present and discuss in silico experiments simulating ED. 138 We simulate ED by temporary blocking and subsequently restoring of the NKA current 139 in the neuronal and astrocyte compartment. The expression for the NKA current is 140 given by, The expression f i NKA describes the inward NKA current. It depends on neuronal Na + concentration and extracellular K + concentration. Details can be found in Materials and Methods. ED is simulated by manipulating the function I max NKA (t), which is the amount of energy available (in %). The function I max NKA (t) is implemented as where we have the U-shaped function such that The parameter β controls the steepness, t start denotes the onset time, and t end denotes 142 the offset time. The parameter P min is the minimum available energy when ED is   ED begins at t = t start min. and ends at t = t end min. while being reduced to a minimum of P min . (B,C) Experimental traces from [8] (left) and the corresponding model simulations (right) with (B) showing intracellular sodium for neurons and astrocytes and (C) showing extracellular sodium and potassium. Empirically adjusted parameters: P min and α 0 e (initial extracellular volume fraction) were chosen by fitting model dynamics qualitatively to Na + and K + concentration time-traces. Here, P min = 50% and α 0 e = 80%. The difference in Na + increase between neurons and astrocytes is attributed to the presence of fast Na + influx through voltage-gated Na + channels in neurons, which are lacking in astrocytes. Please note that scaling axis in panels B and C may slightly differ between experiments and simulations for optimal display purposes.
All the parameters in this model have been sourced from previously published work. 145 However, we set the parameters α e , P min , P n NKA (= P a NKA ) to fit the experimental traces 146 in [8]  presented in this work. 160 We summarize the various simulations performed ahead in Table 2 The parameter α e is the initial extracellular volume ratio and is given by We fix the initial volumes W 0 n and W 0 a , and choose the extracellular volume to be

169
Astrocytes and ion homeostasis 170 We first assess the contribution of the astrocyte to the extra-and intracellular ion 171 homeostasis under physiological conditions. For this, we stimulate our neuron with a 172 10-second long pulse of square wave input with magnitude 25 pA, both with and 173 without a functional astrocyte (Fig. 4). Simulations are performed with a realistic 174 initial extracellular volume fraction for mature rodent cortex [36] by setting α e = 20% 175 and full ATP availability (P min =100%). When the astrocyte is functional, the current 176 injection (black trace) induces a burst of neuronal action potentials (in total 475 during 177 10 s) and transient depolarizations of the astrocyte. The burst firing of neurons is 178 accompanied by a transient decrease in extracellular Na + , while neuronal Na + 179 increases, consistent with experimental data [44]. On the other hand, in our simulation, 180 astrocytic Na + slightly decreases in response to neuronal burst firing. After stimulation 181 subsides, both membrane potentials and ion concentrations return to baseline, see following Na + entry via voltage-gated Na + channels. As astrocytic K + uptake via 187 Kir4.1, NKA and NKCC1, is blocked, K + accumulates in the extracellular space, 188 resulting in a sustained depolarization block of the neuron. As the continued Na + 189 influx is higher than NKA-mediated Na + efflux, the neuron accumulates even more Na +

202
Dynamics after varying periods of ED 203 We know from experimental observations that energy depletion results in accumulation 204 of Na + in neurons and astrocytes while extracellular K + increases see, e.g. [8,12]. At 205 the (thermodynamic) equilibrium, the Nernst potentials of Na + and K + reach equal 206 values. During this evolution towards equilibrium, neurons can generate oscillatory 207 depolarizations, known as anoxic oscillations [4,7]. These drive Ca 2+ into the 208 presynaptic terminal by activating of voltage-gated Ca 2+ channels (VGCCs) which 209 triggers the release of glutamate into the cleft. The efflux of K + does not fully 210 compensate the increase in intracellular Na + . Thus, intracellular Cl − increases to 211 preserve electroneutrality [45]. The resulting osmotic imbalance leads to water influx, 212 i.e., cell swelling [12]. 213 We first simulate such dynamics by setting the available energy P min to 50% of 214 baseline activity at an initial extracellular volume ratio α e = 80% for 5 minutes, shown 215 in Fig. 5-A. Similar to experimental observations, reduced NKA activity leads to an 216 increase in intracellular Na + and extracellular K + . Reactivated NKA restores the Na + Without a functional astrocyte, the neuron depolarizes after the burst, and remains in this state, even if the astrocyte function is restored. Here, we plot neuronal (blue), astrocyte (orange) and extracellular (green) traces against time for several quantities. The initial extracellular volume ratio is α e = 20%. The shaded red area corresponds to periods during which ion transport across the astrocytic plasma membrane is blocked. and K + gradients, so membrane potentials recover to baseline in agreement with 218 experiments [8,32]. Changes in cell volume are below 5%, i.e., they are negligible for 219 this simulation.   for 5 minutes (A) and 15 minutes (B) demonstrates the existence of two stable states: 1) before ED (baseline resting state) and 2) after prolonged ED of 15 minutes (stable depolarized state). Here, we plot neuronal (blue), astrocyte (orange) and extracellular (green) traces against time for several quantities. The initial extracellular volume ratio α e = 80% and minimal energy available P min = 50%. Shaded grey areas correspond to the period where Na + /K + -ATPase (NKA) activity is gradually reduced to P min and restored to baseline, identical to Fig. 3  Further, ion gradients do not return to physiological values, and permanent cell swelling 232 is observed (approximately 20% increase in neuronal, and 10% in astrocytic volume).

233
This occurs once again due to bistability. After a sufficiently long, but transient, period 234 of ED resulting in both neuron and astrocytic reduction of NKA activity, the 235 physiological state is not restored, as reflected by persistent membrane depolarization 236 and increased cell volume. This contrasts with the scenario presented in Fig. 5 Recovery from a depolarization block after energy restoration was addressed 246 experimentally in [32]. Neurons from two different brain regions were subjected to 15 247 minutes of zero energy conditions. Magnocellular neuroendocrine cells recovered, while 248 pyramidal cells remained depolarized post energy restoration (see also Fig. 8). As 249 argued in [46], this may result from different ECS volume fractions. We now explore 250 this hypothesis in our model by changing the baseline value of the parameter α e (the 251 extracellular volume fraction in %). In Fig. 6-A, we simulate ED for 5 minutes for 252 α e = 80% (large ECS) and α e = 20% (small, realistic ECS) and plot the relative 253 neuronal volume change after 30 minutes of energy restoration. The parameter P min is 254 set to 50%. For α e = 80%, we observe no volume increase, which corresponds to a 255 successful recovery after energy restoration, and vice-versa for α e = 80%. Further, we 256 simulate ED during 5 minutes for different values of α e and P min and plot the relative 257 neuronal volume change after 30 minutes of energy restoration. We find that 258 compartmental volumes (and thus, membrane potentials and ion concentrations) recover 259 for larger values of α e . This result implies that neurons and astrocytes surrounded by 260 smaller extracellular spaces exhibit a relatively higher vulnerability to transient ED. 261 We provide further insight into the selective sensitivity to ECS size using bifurcation 262 diagrams shown in Fig. 7-A. We show diagrams for both α e = 20% and α e = 80% as a 263 function of the available energy P min . We plot initial neuronal volume (in %) against 264 P min , after setting I max NKA (t) = P min as constant. We obtain two branches of resting 265 states. Solid (red/blue) lines correspond to stable resting states, while the dashed blue 266 branches correspond to unstable resting states.

267
The bifurcation diagram can be interpreted by following the stable resting states. 268 We start at the blue line at P min = 100%, which corresponds to the baseline 269 physiological state. As we lower NKA activity in the neuron and astrocyte, we move to 270 the left, along the solid blue curve. For values lower than P min = 63.93%, the blue curve 271 does not exist. Decreasing P min , the physiological branch disappears via a limit point curve. When we return to P min = 100%, the pathological state is still stable, implying 286 that recovery is impossible without further intervention. As the pump is already fully 287 activated, this resting state would become unstable only if its maximal strength P min is 288 increased above 145% of its initial value, at which a Hopf bifurcation occurs (inverted 289 triangle), which is biophysically unrealistic. The pathological and physiological resting 290 state curves are much farther apart for α e = 80% than for α e = 20%. Moreover, for 291 α e = 80%, the limit point bifurcation (star) is further to the left (i.e. more severe ED). 292 This shows that the transition to the pathological state occurs at larger values of P min 293 for synapses with a small ECS (small α e ).

NCX block Astrocyte EAAT block
Na v block Here, α e = 80%, P min = 0% and P scale = 2. With these parameter settings, the membrane potential recovers to baseline conditions after the period with ED. Note that during recovery, the experimentally observed oscillations are also faithfully simulated.
At baseline, P scale NKA = 1. When P scale NKA > 1, NKA becomes stronger by a factor P scale NKA . We 299 keep the baseline conditions (ion concentrations, volumes and membrane potentials) 300 fixed by changing the leak permeabilities of K + and Na + in neurons and astrocytes. 301 We first simulate ED for 15 minutes in large ECS (α e = 80%) and realistic ECS 302 (α e = 20%), shown in Fig. 6-B.1,2. Here, we set P scale NKA to 1 and P min to 0%. In both 303 cases, neuronal and astrocyte membranes depolarize and do not recover after energy is 304 restored. Note that it takes longer to reach the peak depolarization for neurons in large 305 ECS, in agreement with observations in the bifurcation diagram ( Fig. 7 NKA being twice as strong. In both cases, the membrane potentials transit to a stable 308 depolarization block during the period of ED. In both cases, the synapse now returns to 309 the physiological resting state after energy restoration. In the case of α e = 80%, the 310 neuronal membrane potential closely mimics experimental findings for magnocellular 311 neuroendocrine cells [32] (see Fig. 8 transits through the small parameter regime where the periodic orbit exists (Fig. 6-B.4). 319 When I max NKA reaches 100%, the periodic orbit disappears, and the membrane potential 320 returns to baseline.

321
Pharmacological blockers 322 We simulate various scenarios observed experimentally where neurons and astrocytes in 323 acute tissue slices of mouse brain were subjected to ED in the presence of an inhibitor 324 of EAATs (TFB-TBOA) or NCX [8]. We first set α e = 80% and P min = 50%, identical 325 to Fig. 5, and present four scenarios, illustrated in Fig. 7 reverses, and thus mediates Na + efflux in low energy conditions, a phenomenon also 331 suggested in experimental observations [8]. Blocking NCX transport promotes the 332 cellular depolarization of the neuron and the astrocyte, driving them to the irreversible 333 pathological state. Third, motivated by experimental evidence [12] and modelling work [11], we explore 335 the potential for recovery after blocking neuronal Na + influx pathways. We first 336 simulate a short duration of ED (P min = 50%) simultaneously blocking the voltage-gated 337 Na + channel, upon which the neuron barely depolarizes (Fig. 7-B.3). Next, we consider 338 the system after it has transited to a stable pathological state after longer ED, shown in 339 ( Fig. 7-B.4). Subsequently, we block the voltage-gated Na + channel of the neuron for 10 340 minutes. We observe that the neuron repolarizes and transits to the physiological 341 resting state, with the restoration of all other ion concentrations and cell volume.

342
Motivated by this, we block neuronal K + efflux pathways to explore the potential for 343 recovery. As before, we consider the system after it has transited to a stable 344 pathological state, after longer ED. Blocking the voltage-gated K + channel for 10 345 minutes ( Fig. 7-B.4), which is the major K + efflux pathway in the neuron, results in 346 neuron and astrocyte repolarisation, and the entire system transits to the stable 347 physiological resting state, as before ED.

348
In Fig. 9, we demonstrate synaptic recovery after blocking the voltage-gated Na + 349 channel and voltage-gated K + channel. We set α e = 80% and subject the neuron to ED 350 for 15 minutes (P min = 50%). Then, as the neuron has transited to the pathological  Finally, we mention that, apart from neuronal Na + and K + voltage-gated channels, 362 blocking any other ion channel did not assist in recovery from the pathological state, see 363 S3 Fig.   364

365
We present a detailed biophysical model of energy-dependent ion fluxes in different 366 compartments and of changes in cellular membrane potentials of the tripartite synapse 367 to further our understanding of their dynamics in low energy conditions. We calibrate 368 the model to Na + and K + concentration time-traces obtained from in-situ chemical 369 ischemia experiments [8]. We demonstrate that astrocyte function is instrumental in Neuronal stimulation upon recovery from ED shows different glutamate transients as compared to neuronal stimulation in a pathological state. Recovery from ED is achieved by blocking neuronal voltage-gated Na + channels (A) and blocking neuronal voltge-gated K + channels (B). Glutamate in the cleft (green trace) and neuronal membrane potential (blue trace) are shown, in response to neuronal excitation in physiological and pathological conditions. First, ED is simulated between t = 5 and t = 20 minutes (P min = 50%, α e = 80%). Then, neurons are subjected to 25 pA square wave input for 10 seconds, as indicated by the black trace. The system is then brought back to the physiological state by blocking voltage-gated Na + channels (shaded green area). After a little more than an hour, the neurons are subjected again, to 25 pA square wave input for 10 seconds.
maintaining physiological ion gradients for action potential generation and proper 371 synaptic transmission. Crucially, the model indicates that surrounding extracellular 372 volume size and baseline NKA pumping capability controls the ischemic vulnerability of 373 the neuron-astrocyte interaction. Further, bifurcation analysis shows how the bistability 374 depends on extracellular volume. Finally, we show that intervention through blocking 375 voltage-gated Na + channels can revive the system from a pathological state.

376
Loss of synaptic function depends on the depth and duration of 377 ischemia 378 In resting conditions, our model shows astrocyte and neuron membrane potentials close 379 to the K + Nernst potential, similar to experimental observations. Loss of 380 NKA-mediated Na + and K + transport results in intracellular Na + and extracellular K + 381 accumulation and membrane depolarization, in accordance with experimental 382 observations [8,32] and previous simulations [11]. As expected [47], astrocytic work performed in brain tissue slices [44] has shown that most astrocytes undergo an 387 increase in intracellular Na + during neuronal bursting activity. However, in about a 388 third of cells, the increase was rather brief and was followed by an undershoot in [Na + ] i . 389 This undershoot could be mimicked by increasing extracellular K + , indicating that it 390 was due to activation of astrocyte NKA. Moreover, the combined addition of glutamate 391 and high K + caused a biphasic increase-decrease in some cells, indicating that both 392 processes -activation of NKA resulting in an export of Na + and activation of glutamate 393 uptake resulting in an import of Na + -counteract each other. Notably in our model,

394
NKA-induced export of Na + is higher than EAAT-mediated Na + uptake in astrocytes. 395 If the partial and transient ED is sufficiently long, no recovery occurs. During low 396 energy conditions, glutamate in the cleft can rise to several hundred µM or even a few 397 mM [48] and excitotoxic cell death will eventually follow [16]. The persistent increase in 398 glutamate during the duration of ED is faithfully reproduced in our simulations 399 (Fig. 5-B). Our simulations in Fig. 6 show that a smaller ECS makes the neurons and astrocytes 409 more likely to depolarize as the extracellular ion gradients will change faster; further, a 410 smaller ECS makes it more likely that neurons and astrocytes remain in the 411 pathological state, which is further illustrated by the bifurcation diagrams in Fig. 7. We 412 conclude that synapses surrounded by smaller ECS are more vulnerable to ED. For a 413 larger ECS, the basins of attraction are farther apart in state space, making it less likely 414 for the trajectory to escape from the physiological equilibrium to the pathological one. 415 These simulations predict that synapses with smaller extracellular spaces are less likely 416 to recover from ED. As the ECS size significantly declines with ageing [33,36], this may 417 aggravate ischemic damage to the aged brain [49]. 418 We also simulated the recovery observed in magnocellular neurons by changing the 419 baseline pump strength with the parameter P scale NKA in Fig. 8-B. The neuron makes a few 420 transient oscillations before returning to baseline, which was indeed observed in [32]. In 421 this case, the Hopf bifurcation (inverted triangle) in Fig. 7 Andrew, showing that transient oxygen-glucose deprivation has differential effects on 426 the potential for recovery of neurons in the hypothalamus and thalamus [32]. In similar 427 experiments [50,51], the difference in NKA activity due to varying α-isoform expression 428 was highlighted as a possible governing factor for recovery. 429 We remark that our findings are at variance with an earlier modelling study, showing 430 that smaller ECS size improves the recovery of neurons from transient ATP 431 depletion [46]. In this model [46], that comprised Na + , K + and Cl − currents of a 432 presynaptic neuron and the Na + /K + -ATPase, the role of glia was limited to K + 433 buffering and passive extracellular K + diffusion. This modelling choice may underlie the 434 divergent findings in our simulations, where an extended role for astrocytes together 435 with ECS volume restrictions allows persistent K + accumulation in the ECS.

436
Neurons can be rescued from the pathological equilibrium state 437 In agreement with experimental observations [32], our simulations also show that a 438 pathological state may persist after energy restoration. Phenomenologically, this 439 persistence results from a persistent Na + current (the window current) at a membrane 440 voltage near -30 mV [11,52], that is too large to be counteracted by the Na + /K + pump. 441 The system can be rescued from this state, however, by a temporary blockage of 442 voltage-gated Na + and K + currents, see Fig. 7-B and Fig. 9. We show that blocking 443 Na + influx via voltage-gated Na + channels and K + efflux via voltage-gated K + 444 channels serves to recover neuronal and astrocytic membrane potentials and Na + /K + 445 homeostasis. Moreover, in this recovered state, healthy synaptic function is resumed as 446 demonstrated by exciting the neuron in Fig. 9, where glutamate transients also return 447 to the physiological range, identical to Fig. 4-A.

448
Our prediction regarding voltage-gated Na + channels agrees with previous modelling 449 work, see [11]. While a study in rats after middle cerebral artery occlusion showed 450 significantly improved neurological outcome after treatment with the Na + channel 451 blocker valproic acid [53], to our knowledge, no experimental data have been reported 452 that explicitly support this model prediction. Our prediction regarding voltage-gated 453 K + channels is novel with respect to current modelling literature. An experimental 454 study showed that using the K + channel blocker tetraethylammonium (TEA) was able 455 to attenuate ischemia-triggered apoptosis in neurons [54]. However, to our knowledge, 456 similar to the case with voltage-gated Na + channels, there is no experimental data that 457 directly supports our prediction regarding K + channels. 458 Decrease in NKA activity in neurons and astrocytes creates an imbalance of Na + 459 and K + gradients, causing net compartmental Na + influx and K + efflux. In that 460 context, it is not too surprising that our model predicts that blocking major Na + and 461 K + pathways in the neuron presents a potential pathway for recovery from the 462 pathological state post ED. However this is not the case for astrocytes. Blocking any of 463 the Na + or K + pathways post ED in the astrocyte (see supplementary figure S3 Fig)   464 does not perturb the pathological state. This can be attributed to the fact that in our 465 model, there is no astrocytic Na + influx or K + efflux process that contributes to large 466 gradients, compared to the neuronal voltage-gated channels.

467
Relation to other computational models 468 Our work extends the single-neuron formalism in [11] to a neuron-astrocyte interaction 469 that describes biophysical processes in synaptic and somatic compartments. There have 470 been several computational studies of neuronal dynamics in the context of energy 471 failure [7,26,28,29], while other studies have explicitly modelled astrocyte dynamics in 472 the context of spreading depression, astrocyte Ca 2+ signaling and physiological 473 function [22,40,55]. To the best of our knowledge, however, no other computational 474 study has explicitly modelled astrocyte dynamics nor did they include Na + , K + , Cl − , 475 Ca 2+ and glutamate dynamics into a single model in the context of ED. In our work, we 476 combine the dynamics of these five ions to provide a holistic description of ion and 477 volume dysregulation during low energy conditions. Further, the model considers 478 physical laws that arise in limiting cases, such as the Gibbs-Donnan equilibrium, which 479 is reached when pump activity is absent. This provides a platform to extend the current 480 formalism by introducing more ions, cellular compartments and transport mechanisms. concentrations to appear, which speeds up the pathological effects of ED. However, this 487 will not affect the bifurcation diagram structure, but only how tipping between the two 488 basins of attraction occurs during short term ED. We also do not include astrocyte gap 489 junctions which may change their permeability during low energy conditions and 490 operate in the same timescale as membrane depolarization [56]. However, their role in 491 mitigating extracellular K + uptake is debated.

492
Future directions 493 The current formulation of the model incorporates ion transients in synaptic and 494 somatic compartments and predicts bistable behaviour in response to ischemia. Adding 495 another neuron with a postsynaptic terminal will allow us to make predictions regarding 496 synaptic transmission in ischemic conditions and to compare presynaptic versus 497 postsynaptic vulnerability to such conditions. The first predictions from our model 498 about glutamate transients can benefit from models such as [57] to include the 499 glutamate-glutamine cycle, which is critical to the dynamics of neurotransmitter 500 replenishment. Moreover, by introducing pH regulation, also the sodium-bicarbonate 501 cotransporter (NBCe1) [58] and the sodium-proton exchanger (NHE) could be 502 incorporated. These mechanisms will increase Na + transport into astrocytes, and might 503 result in a net increase in astrocytic Na + with physiological stimulation, improving the 504 current model.  The model describes the dynamics of molar amounts N i X of the ions Na + , K + , Cl − , Ca 2+ and glutamate, and compartmental volumes W i , for i = {n, a}. The Table 3. Common model parameters along with sources. Units are presented in the same manner as they are implemented in the Python code. All adjusted parameters are in the same order of magnitude as their original counterparts. corresponding membrane potentials follow from the relation, where i = {n, a} and z X is the valence of the ion/species X. The ions X also include impermeable ions A − and B + in each of the somatic compartments. This is necessary to maintain a non-zero resting membrane potential across the semi-permeable neuronal and astrocyte membranes. The molar quantities of these ions are unknown, we estimate them from baseline conditions in the section Estimating parameters. Assuming that Na + , K + and Cl − concentrations are the same in the somatic and synaptic compartments, we get The volumes W ps , W c , W pap are constant. We further assume that all glutamate and Ca 2+ in the neurons and astrocytes are located in the synaptic compartments. Thus, In the following sections, we elaborate on how the dynamics of ion amounts N i X an 525 volumes W i are described in the model. The values of a few common parameters, such 526 as fixed volumes and physical constants, are presented in Table 3.

527
Neuronal dynamics 528 The following currents/fluxes are used to describe neuronal somatic dynamics: Parameters corresponding to all processes in the neuronal compartment can be found in 537 Table 4.
538 Table 4. Model parameters for the neuronal compartment, along with sources. Units are presented in the same manner as they are implemented in the Python code. All adjusted parameters are in the same order of magnitude as their original counterparts.

[dimensionless]
Ratio of extracellular to intracellular proton concentration L n Neuronal membrane water permeability [11] Voltage-gated currents and leak channels 539 The Goldman-Hodgkin-Katz (GHK) currents are solutions to the Nernst-Planck equations that describe the electrodiffusive flux of ions across a permeable membrane when we assume membrane homogeneity and instantaneous and independent movement of ions. We use the GHK currents to model voltage-gated and leak currents [62]. The gating variables are from [37] and are similar to those in the Hodgkin-Huxley model. The currents are as follows, The dynamical equations for q are given by Eq. (1) Note that this is different from the subscript n, which refers to the neuronal somatic compartment. They have the same expressions as in [37]. For the Cl − gated current, we adopt the choice from [11], All the leak currents are modelled as GHK currents Active transport across neuronal membrane: The NKA exchanges three Na + for two K + by consuming one molecule of adenosine triphosphate (ATP). It is modelled as a function of intracellular Na + and extracellular K + as in [61] by the following flux, where and where P n NKA is the NKA permeability or the pump strength. P scale NKA is the scaling factor which we vary in simulations to scale the strength of NKA and the function I max N KA (t) is used to simulate energy derivation for some time period. They are also explained in the section Simulations and model calibration. The corresponding Na + and K + currents are given by I Na + ,n NKA = 3I n NKA (t), March 19, 2021 23/41 Secondary active transport across neuronal membrane: The K + -Cl − -cotransporter (KCC) is a symporter that allows one Cl − ion and K + to leave the neuron, along its concentration gradient. It is the main extruder for Cl − ions in the neuron, thereby providing a counterforce to the gated Cl − channel, which mediates a massive influx of Cl − after Na + loading in neurons [63][64][65][66][67]. We model the flux as the difference of the K + and Cl − Nernst potentials as in [40], The corresponding K + and Cl − currents are given by Astrocyte soma 545 Astrocytes possess a wealth of membrane ion channels and transporters, which allow 546 them to detect, respond and modulate neuronal activity. Major tasks fulfilled by 547 astrocytes at glutamatergic synapses are the regulation of extracellular K + homeostasis 548 and the re-uptake of synaptically-released glutamate [68]. In this section, we describe 549 the incorporated relevant ion channels/cotransporters that govern astrocyte dynamics 550 during physiological conditions and in response to metabolic stress. We use the 551 following currents/fluxes to describe astrocyte somatic dynamics:  Table 5 lists all parameters corresponding to astrocyte fluxes/currents. The weakly inwardly rectifying K + channel Kir4.1 is highly expressed in astrocytes and 562 maintains the resting membrane potential [70,71], close to the K + reversal potential. 563 We choose the model from [69], as it allows the current to vanish at the Gibbs-Donnan 564 condition. This property is not present in recently published models on Kir4.1, such 565 as [72]. The current is given by, where E a K + is the K + reversal potential in the astrocyte,

[dimensionless]
Ratio of extracellular to intracellular proton concentration L a Neuronal membrane water permeability [11] Active transport: Na + /K + -ATPase

569
The NKA in the astrocyte follows the exact same model as that in the neuron, as in [40]. The NKA current in the astrocyte is given by, where and where P a NKA is the astrocyte NKA pump strength. Thus the corresponding Na + and K + currents are given by I Na + ,a NKA = 3I a NKA (t), March 19, 2021 25/41 The astrocyte K + removal mechanism is complemented by the inward 573 Na + -K + -2Cl − -cotransporter (NKCC1), which is highly expressed in astrocytes [68].

574
NKCC1 is a symporter and transports one Na + , one K + and two Cl − into the 575 astrocyte. It is activated by high extracellular K + and plays a major role in astrocyte 576 swelling [73,74]. Astrocyte Cl − regulation also crucially depends on NKCC1 [75,76].

577
The flux is proportional to the difference of Nernst potentials of the respective ions as is 578 done in [40] and is given by Thus the corresponding Na + , K + and Cl − currents are given by Leak currents 580 So far, we only have an inward NKCC1 flux to model movement of Cl − in the astrocyte. 581 We approximate the remaining fluxes as passive electrodiffusive currents via the 582 Goldman-Hodgkin-Katz formula for ion currents, where the formula for GHK (·, ·, ·) is given by Eq. (14), and X = {Na + , K + , Cl − }. EAAT. We assume that the volumes of these synaptic compartments remain small and 588 fixed during the first few hours of metabolic stress. Note that, as previously introduced 589 in Eq. (11) and (12), we assume that all of their ions are confined to the 'synaptic 590 compartments' only, i.e., the presynaptic terminal, synaptic cleft and astrocyte process. 591 We now describe the relevant channels/cotransporters acting in the synaptic 592 compartments and the mechanism of glutamate recycling in the cleft.

594
The re-uptake of synaptically released glutamate is mediated by high-affinity,

595
Na + -dependent glutamate transporters (EAATs). EAATs are expressed by both 596 presynaptic terminals and astrocytes, with astrocytes mediating about 90% of 597 glutamate uptake in the CNS [14]. The cotransporter protein carries one glutamate 598 molecule, three Na + and one H + into the cells in exchange for one K + . The transport 599 yields a net double positive charge influx. We model EAAT in the same way the KCC 600 and NKCC1 cotransporters are modelled. Thus, the EAAT current is written as, March 19, 2021 26/41 The corresponding ion currents are NCX is that it reverses when [Na + ] increases in the respective compartment [8]. 610 We follow the model from [61], which describes the NCX current by The corresponding Na + and Ca 2+ currents are given by Vesicular recycling

611
In response to action potentials that reach the presynaptic terminal and consequent 612 Ca 2+ elevations, glutamate is released into the synaptic cleft. This process involves 613 packing of the neurotransmitter into synaptic vesicles, which fuse with the presynaptic 614 membrane following formation of the SNARE complex [78]. The packing of glutamate 615 into vesicles by vesicular glutamate transporters (VGLUTs) depends on [Cl − ] and on a 616 proton gradient across the vesicular membrane, which is mediated by an 617 ATP-dependent proton pump. It was suggested that VGLUT expression declines with 618 age, although these ideas still remain inconclusive [79]. In this work, we assume that 619 vesicular packing and recycling is not directly energy-dependent. This allows us to see 620 what happens during partial ED when Na + /K + -ATPase is affected, but glutamate 621 continues to be efficiently packed into vesicles.

622
To model vesicular recycling, we combine models from [42] and [43], see Fig. 2. The 623 model proposed by Walter et al. [43] describes the sequential slow-fast mechanism of 624 packing glutamate into vesicles, depending on Ca 2+ elevations. The sequence models 625 the pathway of glutamate from a large storage pool (called depot) to the irreversible 626 fused state, which is when glutamate is released into the synaptic cleft by the 627 interaction of vesicles with SNARE proteins lined up on the presynaptic membrane.

628
The cycle is then completed by plugging in the model by Tsodyks and Markram [42] 629 which models glutamate recycling back into the depot. The equations for the inactive 630 (I) and fused (F) state are then adjusted to include the uptake of glutamate via EAAT 631 and electrodiffusive leak currents. This is done by removing the linear term for 632 recruitment of glutamate from the fused state back into the inactive state, and replacing 633 it with terms for EAAT and leak dynamics. The dynamical equations for the various 634 states of glutamate are given by, (Mod.1) The coefficient k 1 is taken from [43] 636 where K M is the half saturation Ca 2+ concentration in the presynaptic terminal to recruit vesicles from the depot into the non-releasable pool. The coefficients k 2 and k −2 are taken from [43] k 2 (Ca 2+ ) = k 20 + g(Ca 2+ )k 2cat , where the probability for a Ca 2+ -bound catalyst g(Ca 2+ ) is given by [43] 637 From the molar amounts of the various states of glutamate, we can define glutamate 638 concentrations in the neuron and cleft. Thus, and All parameters corresponding to vesicular recycling can be found in Table 6. The 641 expression for N F is derived from conservation laws, see section Conservation laws.

642
Volume dynamics

643
The exact channels for water movement between the extracellular space, neuron and 644 astrocytes is are still debated today [80,81]. In our work, we assume it to depend 645 where ∆π i is the osmotic pressure gradient given by, for X, Y ∈ {Na + , K + , Cl − } and i ∈ {n, a}. From Eq. (1), we then obtain the relation Model equations 650 The dynamics of individual ion amounts based on the described ion currents is given by I Na + ,n G + I Na + ,n NKA + I Na + ,n EAAT + I Na + ,n NCX + I Na + ,n  Table 7. However, we 654 have not described extracellular dynamics yet. These are obtained directly from 655 conservation laws, which we describe in the following section.

656
Conservation laws 657 The model has constant total volume W tot , i.e.
As a consequence we get a conservation law for ionic molar amounts giving where the sum is over all compartments, for each ion X. As the net charge in the 660 system must be zero, we have, at all times where Y contains the impermeable cations and anions. The equations (Cons.1-Cons. 3) give us the three conserved quantities. As a consequence, we can now describe extracellular dynamics from Obtained from experimental traces in [8].
Chosen to be the same as in the presynaptic terminal (Empirical).
Chosen to be the same as in the presynaptic terminal (Empirical where X accounts for all ions (including impermeable ones) in the system. At baseline 668 conditions, Eq. (10) and (Cons.3) provide three more rest conditions, dynamics, we assume constant volumes of the presynaptic terminal, synaptic cleft and 679 perisynaptic astrocytes processes.

680
The various parameters estimated in this section are laid out in Table 8.

681
Python implementation

682
The code is implemented in Python and will be available on GitHub upon publication 683 (for now, kindly use: mkalia94.github.io/TriSyn_11022020.zip). The simulations 684 were made with the CVode solver, implemented in the Python package assimulo. Supporting information S1 Fig. Extension of Fig. 4 . Here, we further plot Cl − and Ca 2+ concentration profiles, along with important ion gradients. (B) shows that astrocytes are critical for maintaining ion homeostasis when Na + /K + -ATPase (NKA) is fully functional. Here, we plot neuronal (blue), astrocyte (orange) and extracellular (green) traces against time for several quantities. The initial extracellular volume ratio α e = 20%. Shaded red area corresponds to periods during which ion transport across the astrocytic plasma membrane is blocked. Neurons are subjected to a 25 pA square wave input, as indicated by the black trace. The burst contains 475 action potentials. The green extracellular calcium trace has an offset of 1.8 mM.