An electrodiffusive neuron-extracellular-glia model with somatodendritic interactions

Computational modeling in neuroscience has largely focused on simulating the electrical activity of neurons, while ignoring other components of brain tissue, such as glial cells and the extracellular space. As such, most existing models can not be used to address pathological conditions, such as spreading depression, which involves dramatic changes in ion concentrations, large extracellular potential gradients, and glial buffering processes. We here present the electrodiffusive neuron-extracellular-glia (edNEG) model, which we believe is the first model to combine multicompartmental neuron modeling with an electrodiffusive framework for intra- and extracellular ion concentration dynamics in a local piece of neuro-glial brain tissue. The edNEG model (i) keeps track of all intraneuronal, intraglial, and extracellular ion concentrations and electrical potentials, (ii) accounts for neuronal somatic action potentials, and dendritic calcium spikes, (iii) contains a neuronal and glial homeostatic machinery that gives physiologically realistic ion concentration dynamics, (iv) accounts for electrodiffusive transmembrane, intracellular, and extracellular ionic movements, and (v) accounts for glial and neuronal swelling caused by osmotic transmembrane pressure gradients. We demonstrate that the edNEG model performs realistically as a local and closed system, i.e., that it maintains a steady state for moderate neural activity, but experiences concentration-dependent effects, such as altered firing patterns and homeostatic breakdown, when the activity level becomes too intense. Furthermore, we study the role of glia in making the neuron more tolerable to hyperactive firing and in limiting neuronal swelling. Finally, we discuss how the edNEG model can be integrated with previous spatial continuum models of spreading depression to account for effects of neuronal morphology, action potential generation, and dendritic Ca2+ spikes which are currently not included in these models. Author summary Neurons communicate by electrical signals mediated by the movement of ions across the cell membranes. The ionic flow changes the ion concentrations on both sides of the cell membranes, but most modelers of neurons assume ion concentrations to remain constant. Since the neuronal membrane contains structures called ion pumps and cotransporters that work to maintain close-to baseline ion concentrations, and the brain contains a cell type called astrocytes that contribute in keeping an appropriate ionic environment for neurons, the assumption is justifiable in many scenarios. However, for several pathological conditions, such as epilepsy and spreading depression, the ion concentrations may vary dramatically. To study these scenarios, we need models that account for changes in ion concentrations. In this paper, we present what we call the electrodiffusive neuron-extracellular-glia model (edNEG), which keeps track of all ions in a closed system containing a neuron, the extracellular space surrounding it, and an astrocytic “domain”. The edNEG model ensures a complete and consistent relationship between ion concentrations and charge conservation. We envision that the model can be used to study a range of pathological conditions such as spreading depression and, hence, be of great value for the field of neuroscience.

Computational modeling in neuroscience has largely focused on simulating the electrical 2 activity of neurons and networks of such, while ignoring other components of brain 3 tissue, such as glial cells and the extracellular space. Within that paradigm, 4 biophysically detailed neuron models are typically based on a combination of a 5 Hodgkin-Huxley type formalism for membrane mechanisms (see, e.g., [1,2]), and cable 6 theory for how signals propagate in dendrites and axons (see, e.g., [3,4]). Two 7 underlying assumptions in these standard models are that (i) the extracellular space 8 (ECS) is isopotential and grounded, and thus does not affect the neurodynamics, and 9 (ii) that the concentrations of the main charge carriers (Na + , K + , and Cl − ) remain 10 constant over the simulated period. 11 The assumptions (i-ii) are never strictly true. Neuronal activity does give rise to 12 electric fields in the ECS, and in principle, a field will affect the membrane potential 13 dynamics of both the neuron that gave rise to the field and of its neighbors. Such 14 so-called "ephaptic" effects have been the topic of many studies (see, e.g., [5][6][7][8][9][10][11][12]). 15 Furthermore, electrical signals in neurons are generated by transmembrane ion fluxes, 16 which will alter both intra-and extracellular ion concentrations. This may change ionic 17 reversal potentials, and the effect that this may have on neurodynamics has also been 18 the topic of many studies (see, e.g., [13][14][15][16][17]). 19 As homeostatic mechanisms, such as ion pumps and cotransporters, strive to 20 maintain ion concentrations close to constant baseline levels [18], and as ECS potentials 21 tend to be very small compared to the membrane potentials of neurons [9], the 22 assumptions (i-ii) are still warranted under many conditions. However, there are also 23 many conditions where these assumptions are not justified. For example, in tightly 24 packed bundles of axons, it is likely that the activity in one axon may affect its 25 neighbors both (ephaptically) through the electric field that it evokes [6,10], and 26 through the ion concentration changes it generates in the narrow ECS separating 27 them [16]. On the much larger spatial scale of brain tissue, spreading depression (SD) 28 and a number of related pathological conditions are associated with dramatic shifts in 29 the K + concentration and giant DC-like voltage gradients in the ECS, which may be as 30 large as several tens of millimolar and millivolts, respectively [19][20][21][22][23][24]. The 31 pathophysiology of SD is believed to largely depend on the dynamics of extracellular 32 neuron-extracellular-glia (edNEG) model, the neuron and glial domain interact through 72 a shared ECS. The edNEG model thus includes the main machinery responsible for ion 73 concentration dynamics in a "unit" piece of brain tissue, i.e., corresponding to a single 74 neuron, and the ECS and glial ion uptake that it has to its disposal. The edNEG model 75 has six compartments, two for each of the three domains. It has the functionality that it 76 (1) keeps track of all ion concentrations (Na + , K + , Ca 2+ , and Cl − ) in all compartments, 77 (2) keeps track of the electrical potential in all compartments, (3) has different ion 78 channels in neuronal soma and dendrites so that the neuron can fire somatic action 79 potentials (APs) and dendritic calcium spikes, (4) contains the neuronal and glial 80 homeostatic machinery that maintains a realistic dynamics of the membrane potential 81 and ion concentrations, (5) accounts for transmembrane, intracellular and extracellular 82 ionic movements due to both diffusion and electrical migration, and (6) accounts for 83 cellular swelling of neurons and glial cells due to osmotic pressure gradients. 84 In this first implementation, we study the edNEG model as a closed system, i.e., ions 85 are conserved and confined to stay within the six-compartment system. We focus on 86 The model for the glial domain was taken from a previous model for astrocytic spatial 115 buffering [38] and added to the edPR model so that both the neuron and glial domain 116 interacted with the ECS. Unlike the previous neuron [17] and glial [38] models that it 117 was based upon, the edNEG model was constructed so that it also accounted for cellular 118 swelling due to osmotic pressure gradients. We implemented the edNEG model using 119 the electrodiffusive KNP framework [17,38], which consistently outputs the voltage-and 120 ion concentration dynamics in all compartments. 121 The edNEG model is depicted in Fig 1. Both the neuron and glial domain contained 122 cell-specific and ion-specific passive leakage channels, cotransporters, and ion pumps 123 that ensured a homeostatic ion balance in the system. The neuron contained additional 124 active ion channels that were different in the somatic versus dendritic compartment, 125 making it susceptible to fire somatic action potentials and dendritic Ca 2+ spikes. Both 126 glial compartments contained inward rectifying K + channels. All included membrane Ions of species k were carried by two types of fluxes: transmembrane (index m) fluxes (j k,msn , j k,mdn , j k,msg , j k,mdg ) and intra-domain fluxes in the neuron (j k,in ), the ECS (j k,e ), and the glial domain (j k,ig ). An electrodiffusive framework was used to calculate ion concentrations and electrical potentials in all compartments. (B) The neuronal membrane contained the same mechanisms as in [17]. Active ion channels were taken from [33]. The soma contained Na + and K + delayed rectifier currents (I Na and I K−DR ), and the dendrite contained a voltage-dependent Ca 2+ current (I Ca ), a voltage-dependent K + afterhyperpolarization current (I K−AHP ), and a Ca 2+ -dependent K + current (I K−C ). Both compartments contained Na + , K + , and Cl − leak currents (I leak ), 3Na + /2K + pumps (I pump,n ), K + /Cl − cotransporters (I KCC2 ), and Na + /K + /2Cl − cotransporters (I NKCC1 ), modeled as in [44]. They also contained Ca 2+ /2Na + exchangers (I Ca−dec ), mimicking the Ca 2+ decay in [33] and modeled like in [17]. The glial membrane mechanisms were taken from [38], and they were the same in both compartments. They included Na + and Clleak currents (I leak ), inward rectifying K + currents (I K−IR ), and 3Na + /2K + pumps (I pump,g ).
To verify that we preserved the characteristic firing properties of the Pinsky-Rinzel 140 model when we made it ion conserving and embedded it within the edNEG model, we 141 implemented two versions of the edNEG model, one with a strong coupling between the 142 soma and dendrite layers, and one with a weak coupling (see Methods for definition of 143 weak and strong coupling). When we stimulated the two versions with constant current 144 injections to the neuronal soma, they elicited spikes that were similar to that of the default in all other simulations in this paper.

150
When we in the edNEG model combined a previous neuron model [17] and a previous 151 glial model [38], and made them share ECS, the original resting state of both the 152 previous models were disturbed. To obtain realistic resting potentials in the new system, 153 we had to re-tune selected parameters (see section titled Model tuning for details). The 154 existence of a realistic resting state for the tuned edNEG model is verified in Fig 2E. It 155 shows a simulation where the neuron was stimulated between t = 10 s and t = 20 s with 156 a constant current injection that made it fire at 1 Hz. Both the neuron and glial domain 157 stayed at their resting potentials of approximately −70 mV and −83 mV, respectively, 158 when unstimulated (t < 10 s), and returned to this resting state after the stimulus had 159 been turned off (t > 20 s). The dynamics of the membrane potentials and ion 160 concentrations during on-going activity is analyzed in further detail in the next section. 161 Steady-state firing in the edNEG model 162 In standard (Hodgkin-Huxley type) neuron models, which the original Pinsky-Rinzel 163 model [33] is an example of, the key dynamical variable is the membrane potential. In 164 addition to modeling the membrane potential, the edNEG model presented here keeps 165 track of all neuronal, glial, and extracellular ion concentrations, and accounts for 166 changes in cellular and extracellular volume fractions due to osmotic gradients. It also 167 accounts for the effect that changes in these variables may have on neuronal firing 168 properties.

169
When the neuron is active, the exchange of ions due to AP firing will be 170 counteracted by the homeostatic mechanisms striving to restore baseline concentration 171 gradients. Hence, we expect that for moderately low neuronal firing, the edNEG model 172 will enter a dynamic steady-state scenario (S1) where homeostasis is successful, and 173 firing can prevail for an arbitrarily long period of time without ion concentrations 174 diverging far off from baseline. We also expect that for a too-high neuronal activity 175 level, the edNEG model will enter a scenario (S2) where the homeostatic mechanisms 176 fail to keep up, and where gradual changes in ion concentrations will lead to gradual 177 changes in neuronal firing properties, and eventually to ceased AP firing.

178
The existence of a dynamic steady-state scenario (S1) is illustrated in Fig 3, which 179 shows how selected variables vary in the edNEG model during a 1400 s simulation. The 180 neuron received a stimulus from t = 1 s to t = 600 s that made it fire at 1 Hz ( Fig 3A). 181 To examine the steady-state scenario (S1), we have divided the simulation into four  Fig 3B-D), covering the last 10 s of the simulation, 187 when the system had returned to the original resting state. In all these phases, we 188 examined the temporal development of the neuronal ( Fig 3B) and glial ( Fig 3C) reversal 189 potentials, and the neuronal and glial swelling (Fig 3D).

190
In the initial phase, the concentrations of all ion species varied with time due to the 191 influxes and effluxes associated with cellular activity. In Fig 3B-   while the downstroke reflects the homeostatic mechanisms that were active between APs, 201 working to restore the baseline concentrations. The fact that E K showed the largest 202 variations was as expected, as the extracellular K + had the lowest baseline value of all 203 ion species, and therefore experienced the largest relative changes during AP firing.

204
The homeostatic recovery between APs was incomplete during the initial phase, and 205 the reversal potentials zig-zagged away from baseline for each consecutive AP (  activity, so that after a period of regular firing, the system entered a dynamic 208 steady-state phase where the zigs and the zags became equal in magnitude, and the 209 reversal potential did not deviate further from baseline (Fig 3B2,C2).

210
In this simulation, E K deviated by maximally ∼3 mV from the baseline reversal 211 potential, which was not enough to have a visible impact on the regular firing of the 212 neuron. Hence, the edNEG model supported a steady-state scenario (S1), where the 213 neuron could fire regularly and continuously without dissipating its concentration 214 gradients. For firing in scenario S1, the neuron in the edNEG model performs similarly 215 to the original Pinsky-Rinzel model [33], which does not model ion concentrations, but 216 assumes that they remain constant.

217
When the stimulus was turned off, the recovery phase started. The membrane 218 potentials returned rapidly to values very close to the resting potential (Fig 3A), while 219 the ion concentrations (and thus the reversal potentials) returned more slowly towards 220 baseline (  The fact that the membrane potentials were almost constant during the recovery of the 225 reversal potentials, indicates that the ion concentration recovery was due to a close-to 226 electroneutral exchange of ions over the neuronal and glial membranes. Hence, the 227 edNEG model predicts that "memories" of previous spiking history may linger in a 228 neuron for several minutes, in the form of altered concentrations, even if it appears to 229 have returned to baseline by judging from its membrane potential.

230
When the ion concentrations changed, so did the osmotic pressure gradients. This 231 caused the neuron and glial domain to swell over the simulated time course (Fig 3D). The existence of a scenario (S2) where the homeostatic mechanisms fail to keep up with 241 the neuronal exchange is illustrated in Fig 4. There, the neuron received a strong input 242 current (150 pA) for three seconds, which gave it a high firing rate ( Fig 4A1). While the 243 neuron fired, ion concentrations gradually changed, leading to changes in ionic reversal 244 potentials ( Fig 4B1,C1), which in turn caused a gradual depolarization of the neuron 245 and made it fire even faster. The neuron could tolerate this strong input for only a little 246 more than 2 s before it became unable to re-polarize to levels below the AP firing 247 threshold, and the firing ceased due to a permanent inactivation of the AP generating 248 Na + channels. This condition, when a neuron is depolarized to voltage levels making it 249 incapable of eliciting further APs, is known as depolarization block. It is a well-studied 250 phenomenon, often caused by high extracellular K + concentrations [46]. This kind of

265
In this simulation, the system never returned to its baseline resting state, and the 266 neuron never regained its ability to elicit APs. This has previously been referred to as a 267 wave-of-death-like dynamics [44,47]. It also resembles the neural dynamics seen under 268 the onset of SD [44], but during SD, neurons tend to recover baseline activity after 269 about one minute as the SD wave passes [23]. Putatively, this recovery depends on K + 270 being transported away from the local region by ECS electrodiffusion and spatial 271 buffering through the astrocytic network, and quite likely also vascular clearance. As the 272 edNEG model studied here represented a local and closed system, such spatial riddance 273 of K + did not occur, but we anticipate that recovery might be observed if the edNEG 274 model were expanded to a spatially continuous model (see Discussion for more on this). 275

Soma versus dendrite 276
It is known that the leading edge of the SD wave tends to occur in the layers containing 277 the apical dendrites [20,32]. Inspired from this, we wanted to explore if the edNEG block, after which it stayed highest in the soma layer ( Fig 5A2). Similarly, the neuronal 289 dendrite also swoll more than the soma during neuronal firing (Fig 5B1), whereas the 290 somatic swelling caught up after the neuron had entered depolarization block ( Fig 5B2). 291 To some degree, these observations are in agreement with the notion that SD initiates in 292 dendritic layers, although the differences between the layers were admittedly rather 293 small in the edNEG model (see Discussion for further comments).

294
The effect of the glial domain on neuronal tolerance levels 295 As the subsystem containing only the neuron and ECS was studied thoroughly in [17], 296 we here put an emphasis on exploring what difference the new glial domain made for 297 the system, especially in terms of how it (i) affected the neuronal tolerance level for AP 298 firing, and (ii) how it affected the dynamics of ECS K + concentrations and cellular 299 swelling. To do this, we compared two versions of the model, one being the full edNEG 300 model with glia included, and one without glia included. In the latter case, we removed 301 the glial influence by setting all glial membrane conductances and water permeabilities 302 to zero, i.e., we sealed the glial membrane.

303
When comparing, we wanted to make sure that the neuron fired with the same rate 304 in both versions, something that we could not control using a continuous step current  In the results shown, we stimulated the soma with an inward K + current, but the results were very similar when we stimulated with an inward Na + current, or when the stimulus was applied to the dendritic compartment instead of the soma. "sustained" in this context will mean "sustained for at least 90 s". We chose to stop the 316 simulations at t = 90 s, partly to reduce computation time, and partly because this is a 317 typical time window within which an SD wave passes [23].

318
Whereas the membrane potential dynamics at 4 Hz firing were very similar in the 319 versions with and without glia, the ion concentration dynamics were not (Fig 6E). In 320 the case without glia, the ECS K + concentration peaked at 11.4 mM, while it stayed 321 below 7.7 mM when glia was included. In both versions, however, the K + concentration 322 reached a ceiling level and eventually stabilized at a constant value, so that the system 323 entered a dynamic steady state (cf. scenario S1 in Fig 3).

324
To study homeostatic breakdown (cf. scenario S2 in Fig 4), we increased the ceiling level where the neuron entered depolarization block and AP firing ceased.

331
As the increase in the ECS K + concentration was much faster when glia was not 332 present, the version without glia entered depolarization block after less than 5 s of 333 activity (Fig 6D), while the version including glia maintained the 8 Hz AP firing for 334 almost 9 s (Fig 6B). The time points where depolarization block was reached can be 335 seen as a dent in the ECS K + concentration curves (Fig 6F), which occurred at a  We were surprised to observe that the afterhyperpolarization following APs 339 increased during the journey towards depolarization block in Fig 6B, despite the K + 340 reversal potential becoming more depolarized during the simulation. When exploring 341 this phenomenon, we found that the increased afterhyperpolarization was caused by the 342 electrogenic 3Na + /2K + pump, which increased its activity level as the intracellular Na + 343 concentration and ECS K + concentration increased during the simulation (cf. Eq 62), 344 leading to an increased outward current. Such hyperpolarization by the ATP-ase driven 345 3Na + /2K + pump has also been reported in other studies [48,49].  (Fig 6G,H,J). During steady-state firing, the neuron swoll by up to only 363 1.6 % when glia was present, and by 3.7 % when glia was not present (Fig 6J for   364 frequencies below the "dent"-frequency). We presented the edNEG model for local ion concentration dynamics in brain tissue 376 containing a neuronal, extracellular, and glial domain (Fig 1). The model contained 377 essential ion channels and homeostatic mechanisms, and accounted for somatodendritic 378 signaling by neurons, for electrodiffusive ion concentration dynamics within all domains, 379 as well as for neuronal and glial swelling due to concentration-dependent osmotic 380 pressure gradients. 381 We demonstrated that the edNEG model had realistic dynamical properties in the 382 sense that it supported a scenario (S1) when the homeostatic mechanisms could 383 maintain constant ion concentrations so that the neuron could maintain low firing 384 frequencies for an arbitrary long time (Fig 3), and a scenario (S2) when the neuron fired 385 too fast for the homeostatic mechanisms to keep up, so that ionic concentrations 386 gradually changed, leading eventually to the neuron entering depolarization block and 387 losing its ability to generate further action potentials (Fig 4). The first scenario (S1) 388 represents normal physiological conditions and could be modeled fairly well with simpler 389 and more conventional neuronal models assuming constant ion concentrations and 390 reversal potentials. The second scenario (S2) resembles the onset of pathological 391 conditions such as spreading depression (SD) and the wave of death [44] and requires 392 models that explicitly account for variations in ion concentrations. Of course, 393 concentration effects are not only relevant during homeostatic breakdown. As Fig 3   394 showed, concentrations vary at a much slower time course than the membrane potential, 395 and may give neurons "memories" of previous spiking history that may last for several 396 minutes. Furthermore, we also showed that the homeostatic machinery in itself can 397 affect the firing patterns of a neuron, as the afterhyperpolarizations seen in Fig 6A-D   398 were concentration-dependent effects evoked by the electrogenic 3Na + /2K + pump.

399
As the edNEG model was constructed by expanding a previous model [17] by (i) 400 adding a glial domain and (ii) accounting for cellular swelling due to osmotic gradients, 401 we put an extra emphasis on exploring how the presence of glia affected neuronal firing 402 and swelling. In Fig 6, we showed that the glial support increased the tolerance level for 403 neuronal firing and that the neuron could maintain steady-state firing for at least 90 s 404 at frequencies up to 7.6 Hz in the presence of glia, but only up to 4.7 Hz when the glia 405 domain was inactivated. Furthermore, the presence of glia reduced the swelling of the 1956, is that diffusion of K + through the ECS is the main propagator of the 432 wave [25,26]. However, buffering of K + through the glial syncitium [86,87], and 433 electrical drift of K + along the large DC-like voltage shifts [28,88] are also likely to 434 contribute to the wave propagation. Initiation of SD, and the leading edge of the SD 435 wave, are often seen to occur in superficial (dendritic) layers of cortex and hippocampus, 436 suggesting that dendritic membrane mechanisms play an important role for its 437 pathophysiology [19,20,24,[30][31][32].

438
To our knowledge, only one computational model exists that has combined spatial 439 propagation of SD with morphologically detailed neuron models [82]. However, this 440 model was not based on an electrodiffusive formalism, and did not account for effects of 441 extracellular potentials on neurodynamics and K + transport. Other spatial models of 442 SD [27-29, 40, 88] have been inspired by the coarse-grained bi-domain model [89], which 443 previously has been used to simulate cardiac tissue [90,91]. These models are 444 electrodiffusive, and treat brain tissue as a homogeneous, coarse-grained, continuum, 445 making them computationally efficient to allow for large scale simulations of SD 446 propagation. However, they are limited in terms of neuronal detail, as none of them 447 include fast neuronal mechanisms for action potential generation, or account for any 448 morphological aspects of neurons, i.e., they do not account for the differences between 449 dendritic and somatic layers.

450
The edNEG model was highly inspired by the previous tri-domain continuum model 451 by Tuttle et al. 2019 [27], which is the most advanced of the spatial SD models. It 452 includes neurons, glia, and extracellular space (Fig 7A), and it accounts for cellular 453 swelling, a number of (slow) membrane mechanisms, and electrodiffusive ion 454 concentration dynamics. Unlike other local models of ion concentration dynamics in 455 tissue (see, e.g., [13,44,60,71]), the edNEG model was based on the same kind of 456 electrodiffusive formalism as the model by Tuttle et al. 2019 [27], and should in that 457 regard be compatible with the tri-domain continuum framework used there (Fig 7A). 458 We envision that the edNEG model can be integrated with this framework to obtain a 459 tri-times-two-domain model (Fig 7B) that expands the functionality of the original neuron-ECS-glial-brain tissue at a large spatial scale. We anticipate that such a 468 framework will be of great value for the field of neuroscience, partly because it gives a 469 deepened insight into the balance between neuronal firing, ion homeostasis, and glial 470 buffering at a local scale, partly because it may lead to new insight in the physiology of 471 brain tissue in general, but most importantly because it invites detailed mechanistic 472 studies of a number of pathological conditions associated with shifts in extracellular 473 concentrations, such as SD, ischemic or hemorrhagic stroke, traumatic brain injury, 474 migraine, and epileptic seizures [18, 20-23, 85, 92].

475
As a preliminary result, using the increase in the ECS K + -concentration as an 476 indicator of SD initiation, the edNEG model was in agreement with the notion of earlier 477 SD initiation in the dendritic layer, although layer differences were admittedly quite 478 small in our simulations (Fig 5). We note, however, that the model was not in any way 479 tuned to reproduce explicit SD data, and that these layer specificities followed directly 480 from adopting a set of membrane mechanisms from a previous CA3 neuron model [33]. 481 Putatively, larger differences between the layers could be obtained by adjusting either 482 the somatic and dendritic ion channel conductances or the compartment sizes. In the 483 current version, the neuronal soma and dendrite compartment were implemented with 484 identical volumes and surface areas (cf. the version of the Pinsky-Rinzel model 485 presented in [45]). If, for example, the neuronal surface-to-volume area was assumed to 486 be larger in the dendritic layer, which would sound like a plausible assumption, we 487 would expect a faster dissipation of the concentration gradients, and thus a more rapid 488 increase in the ECS K + concentration in the dendritic layer. Hence, if embedded in a 489 tri-times-two continuum model for SD, there would be several approaches to retuning 490 the edNEG model to make it in agreement with specific experimental data. In addition, modules interact with neighboring modules via spatiotemporal electrodiffusion in the glial domain (spatial buffering through glial syncytium), and the ECS domain. Neurons are typically assumed not to interact laterally in such a spatially continuous fashion. The spatial dynamics can, in principle, occur in all directions (3D), but a 1D illustration was used in the figure. (B) Tri-times-two domain model, where the local module has been replaced with the edNEG model so that it has two layers: a somatic (deep) layer, and a dendritic (superficial) layer. The local module then contains six domains: a neuronal, ECS, and glial domain in each of the two layers. Within each layer (xy-plane), the dynamics are modeled in the same fashion as in (A). In addition, the tri-times-two domain model contains between-layer dynamics in terms of electrodiffusive transports intracellularly in neurons, intracellularly in the glial syncytium, and extracellularly. The key novelty is that the tri-times-two domain model can include different ion channels in the soma versus dendrites of neurons, and can capture somatodendritic signaling in neurons.

492
The Kirchoff-Nernst-Planck (KNP) framework for a 493 tri-times-two compartment model 494 In a previous study [17], we derived the Kirchhoff-Nernst-Planck (KNP) 495 framework [16,38,83,93] for a closed-boundary system containing 2 × 2 compartments, 496 representing a soma, a dendrite, and extracellular space (ECS) outside the soma and 497 dendrite. In the edNEG model, we expanded the KNP framework to also include a glial 498 domain, and to account for osmotically induced volume changes. The three domains 499 (neuron + ECS + glia) were all represented as two-compartment models (one 500 compartment in the soma layer and one in the dendrite layer). Within each layer, the 501 neuron and glial domain interacted with the ECS through transmembrane currents 502 (Fig 1). Volume changes were due to osmotic pressure gradients that we implemented as 503 functions of the ionic concentrations (see section titled Volume dynamics). Geometrical 504 parameters, including initial volumes, are listed in Table 1 1437 · 10 −18 m 3 † The parameter α determines the coupling strength between the soma and dendrite layers of the model. Its default value was 2, same as in [17].

506
Two kinds of fluxes transport ions in the system: transmembrane fluxes and axial fluxes. 507 The axial fluxes are driven by electrodiffusion, and we describe them using the 508 Nernst-Planck equation. It follows that the intracellular flux density of the neuron for 509 ion species k is expressed as: In Eq 1, we use the same notation as in [17] so that D k is the diffusion constant, γ k is the fraction of freely moving ions, that is, ions that are not buffered or taken up by the endoplasmatic reticulum, λ i is the tortuosity, which represents hindrances in free diffusion due to obstacles, γ k ([k] dn − [k] sn )/∆x is the longitudinal concentration gradient, z k is the charge number of ion species k, F is the Faraday constant, R is the gas constant, T is the absolute temperature, [k] n is the average concentration, that is, γ k ([k] dn + [k] sn )/2, and (φ dn − φ sn )/∆x is the longitudinal potential gradient. Likewise, the extracellular flux densities and the glial intracellular flux densities are described, respectively, by All ions can move freely in the extracellular and glial space, thus, γ k is not included in 511 Eqs 2 and 3. Diffusion constants, tortuosities, and intracellular fractions of freely 512 moving ions are listed in Table 2.  [36,38] γ Na , γ K , γ Cl (intraneuronal fractions of free ions) 1 [17] γ Ca (intraneuronal fraction of free ions) 0.01 [17] † The table is adopted from [17].

Ion conservation 514
To keep track of all ions in the system, we solve six differential equations for each ion species k. Conservation of ions gives: dN k,de dt = +j k,mdn A m + j k,e A e + j k,mdg A m , where N is the amount of substance, in units of mol. To find the change in N , all ion 515 flux densities are multiplied by the area they go through. The variable j k,m represents 516 the sum of all membrane flux densities of ion species k, and j k,in , j k,e , and j k,ig 517 represent the axial flux densities. To find the ion concentrations, which appear in many 518 of the following equations, we divide the amounts of substance in a compartment by the 519 compartment volume at the beginning of each time step. 520 We insert the Nernst-Planck equation for the axial flux density (Eq 1) into Eq 4 and 521 get: In Eq 10, [k] dn and [k] sn are the intraneuronal ion concentrations of the dendrite and 523 soma, defined as N k,dn /V dn and N k,sn /V sn , respectively. We define the voltage 524 variables φ dn and φ sn below.

525
Six constraints to derive φ 526 If we have four ion species (Na + , K + , Cl − , and Ca 2+ ) in six compartments, we get 24 527 equations to solve (Eqs 4-9 times four) and 30 unknowns (N and φ). We overcome this 528 by defining φ in terms of ion concentrations using a set of constraints similar to those 529 used in [17]. 530 1. Arbitrary reference point for φ. The first constraint is simple; we can choose an 531 arbitrary reference point for φ. We define it to be outside the dendrite, which 532 gives us: 2. Neuronal membrane is a parallel plate capacitor (dendrite). As the second 534 constraint, we use that the membrane is a parallel plate capacitor. This means 535 that it will always separate a charge Q on one side from an opposite charge −Q 536 on the other side. This gives rise to a voltage difference across the membrane where C m is the total capacitance of the membrane, i.e., C m = c m A m , where c m is 538 the capacitance per membrane area. We know, by definition, that 539 φ mdn = φ dn − φ de , and since φ de = 0, we get: We assume bulk electroneutrality, meaning that all net charge in the dendritic 541 compartment must be on the membrane. It follows that where F is the Faraday constant, z k is the charge number of ion species k, [k] dn is 543 the ion concentration, and V dn is the volume. By inserting this into Eq 13, we get 544 3. Neuronal membrane is a parallel plate capacitor (soma). The second constraint 545 also applies to the soma, and gives us the criterion: Here, the outside potential is not set to zero, so this constraint is not sufficient to 547 determine φ sn and φ se separately. 548 4. Glial membrane is a parallel plate capacitor (dendrite layer). The glial membrane 549 is no different than the neuronal membrane when it comes to acting as a parallel 550 plate capacitor, so we get: where we have used that φ de = 0.

5.
Glial membrane is a parallel plate capacitor (soma layer). Constraint number (iv) 553 also applies to the soma layer, and gives us: We can now calculate φ dn and φ dg from Eqs 14 and 16 but to determine φ sn , φ se , 555 and φ sg , we need a sixth constraint. 556 6. Current anti-symmetry. The sixth constraint is to ensure charge (anti-)symmetry. 557 We must define the initial conditions so that the membrane separates a charge Q 558 on one side from an opposite charge −Q on the other side, and the system 559 dynamics so that it stays this way. The membrane fluxes (alone) fulfill this 560 criterion, since a charge that leaves a compartment automatically pops up on the 561 other side of the membrane, making sure that dQ i /dt = −dQ e /dt. For the axial 562 fluxes to fullfill the criterion, we must have that: where i stands for current density. We find expressions for i in , i ig , and i e , by multiplying Eqs 1-3 by F z k and sum over all ion species k. Expressions for the current densities then become: The first term in Eq 19 is the diffusion current density and is defined by the ion 564 concentrations: The second term is the field driven current density where σ n is the conductivity: Likewise, Eq 20 can be written in terms of i diff,ig , i field,ig , and σ g , and Eq 21 in 568 terms of i diff,e , i field,e , and σ e . We combine Eqs 18-21 and obtain: (25) In Eq 25, we know φ dn , φ de , and φ dg from Eqs 14, 11, and 16, and i diff and σ from the ion concentrations. We solve Eqs 15, 17, and 25 to find φ sn , φ sg , and φ se : Neuronal membrane mechanisms 570 The neuronal membrane mechanisms were equal to those in [17]. We list them again 571 here for easy reference.

572
Leakage channels 573 Both neuronal compartments contained Na + , K + , and Cl − leak currents. The flux 574 densities were modeled as in [44]: where k denotes the ion species,ḡ k,leak is the ion conductance, φ m is the membrane 576 potential, E k is the reversal potential, F is the Faraday constant, and z k is the charge 577 number. Reversal potentials are given by the Nernst equation: where R is the gas constant, T is the absolute temperature, γ k is the intracellular

581
Active ion channels 582 We adopted all active ion channels from the Pinsky-Rinzel model [33]. This included Na + and K + delayed rectifier fluxes in the soma (j Na , j DR ), and a voltage-dependent Ca 2+ flux (j Ca ), a voltage-dependent K + afterhyperpolarization flux (j AHP ), and a Ca 2+ -dependent K + flux (j C ) in the dendrite: Here, g k is the ion conductance, φ m is the membrane potential, E k is the reversal potential, F is the Faraday constant, and z k is the charge number for ion species k. We used the Hodkin-Huxley formalism to model the voltage-dependent conductances, with differential equations for the gating variables: July 7, 2020 22/36 and g Na =ḡ Na m 2 ∞ h, g Ca =ḡ Ca s 2 z, , with φ 7 = φ m + 0.03 (54) with φ 8 = φ m + 0.05 and φ 9 = φ m + 0.0535 (57) Both neuronal compartments contained a 3Na + /2K + pump, a K + /Cl − cotransporter (KCC2), a Na + /K + /2Cl − cotransporter (NKCC1), and a 2Na + /Ca 2+ exchanger. The functional forms of the pumps and cotransporters were taken from [44], while the 2Na + /Ca 2+ exchanger was modeled as in [17]: Here, ρ n , U kcc2 , and U nkcc1 are pump and cotransporter strengths, U Ca−dec is the Ca 2+ 586 decay rate, and [Ca +2 ] n,b is the basal Ca 2+ concentration..

587
Glial membrane mechanisms 588 We modeled the glia as an astrocytic domain and adopted the membrane mechanisms 589 from [38]. They included Na + and Cl − leak channels, modeled as in Eq 31, an inward 590 rectifying K + channel, and a 3Na + /2K + pump: Here,ḡ K−IR is the K + ion conductance, φ m is the membrane potential, E K is the K + 592 reversal potential, F is the Faraday constant, z K is the K + charge number, [K + ] e,b is 593 the basal K + concentration in the extracellular space, ∆φ = φ m − E K , E K,b is the 594 reversal potential for K + at basal concentrations, ρ g is the pump strength, and 595 [Na + ] g,threshold and [K + ] e,threshold are the pump's threshold concentrations for Na + and 596 K + , respectively. We included the same set of membrane mechanisms in both astrocytic 597 compartments.

598
Volume dynamics

599
To calculate the osmotically induced volume changes dV /dt, we used the formalism outlined in [42]: Here, G n and G g are the neuronal and glial water permeabilities, respectively, given in 600 units of m 3 /Pa/s, and Ψ is the water potential, given in units of Pa. We assumed the 601 hydrostatic pressure differences to be zero, so that water flow was driven by osmotic 602 pressure differences only, and we calculated the solute potentials from: Here, i is the ionization factor (van't Hoff factor), which is 1 for ions, M is the osmotic 604 concentration of solutes measured in moles per cubic meter, R is the gas constant, and 605 T is the absolute temperature. Equations 74 and 75 follow from the assumption that 606 the total volume did not change, that is, the system was closed.

607
Like in [27], we only considered effects of transmembrane water flow, and 608 intra-domain water flow due to hydrostatic pressures were neglected. As predicted in a 609 previous study, bulk flow at physiological hydrostatic pressure is expected to be low [95]. 610
At each time step, we derived φ algebraically in all six compartments: Volume dynamics was given by:   The edNEG model combined two previous models, one consisting of a neuron and 615 ECS [17], and the other of a glial domain (astrocyte) and ECS [38]. When we combined 616 the models, we set the initial ionic concentrations in the neuron identical to those 617 in [17], the initial ionic concentrations in the glial domain identical to those in [38], and 618 made the two cells share the same ECS where we set the initial concentrations equal to 619 those in the previous glia model [38]. As these initial concentrations (Table 5) differed 620 from the initial ECS concentrations in the previous neuron model [17], the neuron was 621 not in equilibrium with the (new) ECS. This was because the altered ECS 622 concentrations gave rise to altered concentration-dependent activity of the ion pumps, 623 cotransporters, and ionic currents through ion channels. We found that the leakage 624 currents were most important, and that a re-tuning of the leak conductances (ḡ Na,leak,n , 625 g K,leak,n , andḡ Cl,leak,n ) in the neuron model was sufficient to obtain a system with a 626 plausible resting state. The tuning was done by requiring that the initial leakage 627 currents should be identical to those in [17], i.e., we set: July 7, 2020 27/36 with φ m being the resting potential in [17] (−67.7 mV), E k,old being the reversal 629 potential for ion species k at steady state in [17], and E k,new being the reversal 630 potential obtained by the new initial ion concentrations ( reason why the new and original values were not identical was that not only the leakage 638 currents, but also the ion pumps and cotransporters, and to a small extent the active 639 channels, were active at rest.

640
To obtain comparable spike shapes between the edNEG model and the original 641 Pinsky-Rinzel model (Fig 2), we manually tuned the Ca 2+ conductance of the neuron 642 (ḡ Ca ), as well as the coupling strength between the soma and dendrite layers. The 643 coupling strength was regulated by adjusting the intracellular cross-section areas αA i (cf. 644 Table 1) by adjusting the unit-less parameter α. The parameter α was 2 for simulations 645 with stong coupling between the soma and dendrite layers, and 0.51 for simulations with 646 weak coupling. Strong coupling was used in all simulations except in the simulation 647 shown in Fig 2C where it was 0.51. The Ca 2+ conductance (ḡ Ca ), along with the other 648 membrane mechanisms, were the same in all simulations (values as in Table 4).

649
Initial conditions

650
Before tuning the edNEG model, we defined its initial volumes ( Table 1), amounts of 651 ions, membrane potentials, and gating variables ( Table 5, Pre-calibrated column) using 652 values from the two previous models in [17] and [38]. After re-tuning selected 653 parameters (as described in the previous subsection), the system was close to, but not 654 strictly in equilibrium, and for this reason we calibrated the edNEG model for 5000 s.

655
The water permeabilities were set to zero during the calibration. 656 We wrote the final values from the calibration to file (see Table 5, Post-calibrated 657 column) and used them as initial conditions in all simulations shown throughout this 658 paper. Note that the edNEG model takes amounts of ions (in units of mol) as input, 659 while we have listed ion concentrations in Table 5. The post-calibrated values of the ion 660 concentrations correspond to the following reversal potentials: E Na,n = 55 mV, 661 E Na,g = 62 mV, E K,n = −96 mV, E K,g = −88 mV, E Cl,n = −90 mV, E Cl,g = −83 mV, 662 and E Ca,n = 124 mV.

663
To ensure charge symmetry and electoneutrality, we defined a set of static residual 664 charges, based on the initial amounts of ions. These represent negatively charged 665 macromolecules present in real cells. We defined them as constant amounts of ion 666 species X − with charge number z X = −1 and diffusion constant D X = 0. To ensure 667 strict electroneutrality, we did not read residual charges to/from file, but calculated 668 them at the beginning of each simulation.They were given by the following expressions: 669 N X,n = z Na N Na,n,0 + z K N K,n,0 + z Cl N Cl,n,0 + z Ca N Ca,n,0 − φ mn,0 c m A m F , N X,e = z Na N Na,e,0 + z K K K,e,0 + z Cl N Cl,e,0 + z Ca N Ca,e,0 + (φ mn,0 + φ mg,0 ) c m A m F , N X,g = z Na N Na,g,0 + z K N K,g,0 + z Cl N Cl,g,0 − φ mg,0 c m A m F .
Additionally, we introduced a set of static residual molecules to ensure zero osmotic pressure gradients across the membranes at the beginning of each simulation. These were defined as osmotic concentrations of a molecule M: [M] n = (N Na,n,0 + N K,n,0 + N Cl,n,0 + N Ca,n,0 )/V sn,0 , [M] e = (N Na,e,0 + N K,e,0 + N Cl,e,0 + N Ca,e,0 )/V se,0 , [M] g = (N Na,g,0 + N K,g,0 + N Cl,g,0 )/V sg,0 . 1 Values with more decimals included were read to/from file and used in the simulations.
(Available at https://github.com/CINPLA/edNEGmodel_analysis.) † φ m is not an independent state variable, but defined at each time point from the ion concentrations. * Only 1% of the total intracellular Ca 2+ , that is, a 100 nM, was assumed to be free (unbuffered).

Stimulus current 670
We stimulated the neuron like we did in [17], that is, by applying a K + injection current i stim into the soma, and removing the same amount of K + from the corresponding extracellular compartment to ensure ion conservation: Numerical implementation 671 We implemented the code in Python 3.6 and solved the differential equations using the 672 solve ivp function from SciPy with its Runge-Kutta method of order 3 (2). We set the 673