An electrodiffusive, ion conserving Pinsky-Rinzel model with homeostatic mechanisms.

In most neuronal models, ion concentrations are assumed to be constant, and effects of concentration variations on ionic reversal potentials, or of ionic diffusion on electrical potentials are not accounted for. Here, we present the electrodiffusive Pinsky-Rinzel (edPR) model, which we believe is the first multicompartmental neuron model that accounts for electrodiffusive ion concentration dynamics in a way that ensures a biophysically consistent relationship between ion concentrations, electrical charge, and electrical potentials in both the intra- and extracellular space. The edPR model is an expanded version of the two-compartment Pinsky-Rinzel (PR) model of a hippocampal CA3 neuron. Unlike the PR model, the edPR model includes homeostatic mechanisms and ion-specific leakage currents, and keeps track of all ion concentrations (Na+, K+, Ca2+, and Cl-), electrical potentials, and electrical conductivities in the intra- and extracellular space. The edPR model reproduces the membrane potential dynamics of the PR model for moderate firing activity. For higher activity levels, or when homeostatic mechanisms are impaired, the homeostatic mechanisms fail in maintaining ion concentrations close to baseline, and the edPR model diverges from the PR model as it accounts for effects of concentration changes on neuronal firing. We envision that the edPR model will be useful for the field in three main ways. Firstly, as it relaxes commonly made modeling assumptions, the edPR model can be used to test the validity of these assumptions under various firing conditions, as we show here for a few selected cases. Secondly, the edPR model should supplement the PR model when simulating scenarios where ion concentrations are expected to vary over time. Thirdly, being applicable to conditions with failed homeostasis, the edPR model opens up for simulating a range of pathological conditions, such as spreading depression or epilepsy.


Introduction
The neuronal action potential (AP) is generated by a transmembrane influx of Na + , which depolarizes the neuron, followed by an efflux of K + , which repolarizes it. Likewise, all neurodynamics is fundamentally about the movement of ions, which are the charge carriers in the brain. Therefore, it might seem peculiar that most models of neuronal activity are based on the approximation that the concentrations of the main charge carriers (Na + , K + , and Cl − ) do not change over time. This approximation is, for example, incorporated in the celebrated Hodgkin-Huxley model [1], and a large number of later models based on a Hodgkin-Huxley type formalism (see, e.g., [2][3][4][5][6][7]).
Setting the ion concentrations to not change over time is often a fairly good approximation. The reason is that the number of ions that need to cross the membrane to charge up the neuron by, say, an AP worth of millivolts, is too small to have any notable impact on ion concentrations on either side of the membrane (see, e.g., Box 2.2 in [8]), meaning that concentration changes on a short time scale can be neglected. On a longer time-scale, the ionic exchange due to APs (or other neuronal events), is normally reversed by a set of homeostatic mechanisms such as ion pumps and cotransporters, which work to maintain constant baseline concentrations. In Hodgkin-Huxley type models, the large number of ion pumps, cotransporters and passive ionic leakages that strive towards maintaining baseline conditions are therefore not explicitly modeled. Instead, they are simply assumed to do their job and are grouped into a single passive and non-specific leakage current I leak = g leak (ϕ m −E leak ), which determines the cell's resting potential (for a critical study of this approximation, see [9]).
Another approximation commonly applied by modelers of neurons is that the extracellular potential is constant and grounded (ϕ e = 0) so that the only voltage variable that one needs to worry about when simulating neurodynamics is the transmembrane potential (ϕ m ). This assumption is implicit in the majority of morphologically explicit models of neurons, where the (spatial) signal propagation in dendrites and axons are computed using the cable equation (see, e.g., [10][11][12]). Cable-equation based, multicompartmental neuronal models are widely used within the field of neuroscience, both for understanding dendritic integration and neuronal response properties at the single neuron level (see, e.g., [3,4,6,7]) and for exploring the dynamics of large neuronal networks (see e.g., [13][14][15]). They are even used in the context of performing forward modeling of extracellular potentials, such as local field potentials (LFP), the electrocorticogram (ECoG), and electroencephalogram (EEG) (see, e.g., [16][17][18]), despite the evident inconsistency involved when first computing neurodynamics under the approximation that ϕ e = 0 (Fig 1A), and then in the next step using this dynamics to predict a nonzero ϕ e (Fig 1B). The approximation is nevertheless useful since ϕ e is typically so much smaller than Accurate modeling of such systems thus requires a unified, electrodiffusive framework that ensures a physically consistent relationship between ion concentrations, charge density, and electrical potentials.
Intra-or extracellular electrodiffusion is not an issue in single-compartment models, of which there are quite a few that incorporate ion concentration dynamics in a more or less consistent way [28,29,[33][34][35][36][37][38][39][40][41][42][43][44][45][46][47]. Single compartment models are useful in many aspects. However, in order to represent morphological features of neurons, such as e.g., differential expression of ion channels in the soma versus dendrites, or account for transport processes in the space inside or outside neurons, one needs models with more than a single compartment. Among the several morphologically explicit models that have included homeostatic machinery and explicitly simulated ion concentration dynamics (see e.g., [27,[48][49][50][51][52][53][54][55][56][57]), neither have accounted for the electrodiffusive coupling between the movement of ions and electrical potentials (see Results section titled Loss in accuracy when neglecting electrodiffusive effects on concentration dynamics). Hence, to our knowledge, no morphologically explicit neuron model has so far been developed that ensures biophysically consistent dynamics in ion concentrations and electrical potentials during long-time activity, although useful mathematical framework for constructing such models are available [58][59][60][61][62].
The goal of this work is to propose what we may refer to as "a minimal neuronal model that has it all". By "has it all", we mean that it (1) has a spatial extension, (2) considers both extracellular-and intracellular dynamics, (3) keeps track of all ion concentrations (Na + , K + , Ca 2+ , and Cl − ) in all compartments, (4) keeps track of all electrical potentials (ϕ m , ϕ e , and ϕ i -the latter being the intracellular potential) in all compartments, (5) has differential expression of ion channels in soma versus dendrites, and can fire somatic APs and dendritic calcium spikes, (6) contains the homeostatic machinery that ensures that it maintains a realistic dynamics in ϕ m and all ion concentrations during long-time activity, and (7) accounts for transmembrane, intracellular and extracellular ionic movements due to both diffusion and electrical migration, and thus ensures a consistent relationship between ion concentrations and electrical charge. Being based on a unified framework for intra-and extracellular dynamics (Fig 1C), the model thus accounts for possible ephaptic effects from extracellular dynamics, as neglected in standard feedforward models based on volume conductor theory (Fig 1A and 1B). By "minimal" we simply mean that we reduce the number of spatial compartments to the minimal, which in this case is four, i.e., two neuronal compartments (a soma and a dendrite), plus two extracellular compartments (outside soma and outside dendrite). Technically, the model was constructed by adding homeostatic mechanisms and ion concentration dynamics to an existing model, i.e., the two-compartment Pinsky-Rinzel (PR) model [3], and embedding in it a consistent electrodiffusive framework, i.e., the previously developed Kirchhoff-Nernst-Planck framework [31,32,60,62]. For the remainder of this paper, we refer to our model as the electrodiffusive Pinsky-Rinzel (edPR) model.
The remainder of this article is organized as follows. First, we present the edPR model and illustrate the numerous variables that it can simulate. Next, we show that the edPR model can reproduce the firing properties of the original PR model. By running long-time simulations (several minutes of biological time) on both models, we identify the firing conditions under which the two models maintained a similar firing pattern, and under which conditions concentration effects became important so that dynamics of the edPR model diverged from the original PR model over time. Finally, we use the edPR model to explore the validity of some important assumptions commonly made in the field of computational neuroscience, regarding the decoupling of electrical and diffusive signals. We believe that the edPR model will be of great value for the field of neuroscience, partly because it gives a deepened insight into the balance between neuronal firing and ion homeostasis, partly because it lends itself to explore under which conditions the common modeling assumption of constant ion concentrations is warranted, and most importantly because it opens for more detailed mechanistic studies of pathological conditions associated with large changes in ion concentrations, such as epilepsy and spreading depression [22][23][24][25].

An electrodiffusive Pinsky-Rinzel model
The here proposed electrodiffusive Pinsky-Rinzel (edPR) model is inspired by the original Pinsky-Rinzel (PR) model [3], which is a two-compartment (soma + dendrite) version of a CA3 hippocampal cell model, initially developed by Traub et al. [2]. In the original PR model, the somatic compartment contains Na + , and K + delayed rectifier currents (I Na and I K−DR ), while the dendritic compartment contains a voltage-dependent Ca 2+ current (I Ca ), a voltagedependent K + afterhyperpolarization current (I K−AHP ), and a Ca 2+ -dependent K + current (I K −C ). In addition, both compartments contain passive leakage currents. Despite its small number of compartments and conductances, the PR model can reproduce a variety of realistic firing patterns when responding to somatic or dendritic stimuli, including somatic APs and dendritic calcium spikes.
In the edPR model, we have adopted all mechanisms from the original PR model. In addition, we have (i) made all ion channels and passive leakage currents ion-specific, (ii) included 3Na + /2K + pumps (I pump ), K + /Cl − cotransporters (I KCC2 ), Na + /K + /2Cl − cotransporters (I NKCC1 ), and Ca 2+ /2Na + exchangers (I Ca−dec ), and (iii) included two extracellular compartments (outside soma + outside dendrite). To compute the dynamics of the edPR, we used an electrodiffusive KNP-framework for consistently computing the voltage-and ion concentration dynamics in the intra-and extracellular compartments [60]. The model is summarized in

Key dynamical variables in the electrodiffusive Pinsky-Rinzel model
While the key variable in the original PR model is the membrane potential ϕ m , the edPR model allows us to compute a multitude of variables relevant to neurodynamics. The functionality of the edPR model is illustrated in Fig 3, which shows a 60 s simulation where the model fires at 1 Hz for 10 s. We have plotted a selection of output variables, including the membrane potentials ( Fig 3A and 3B), extracellular potentials ( Fig 3C and 3D), the dynamics of all ion concentrations in all compartments (Fig 3E-3H), concentration effects on ionic reversal potentials (Fig 3I-3J), concentration effects on the electrical conductivity of the intra-and extracellular medium (Fig 3K), and ATP consumption (Fig 3L) of the 3Na + /2K + pumps and Ca 2+ /2Na + exchangers.
Unlike neuronal models based on cable theory, where ϕ e is assumed to be zero so that ϕ m = ϕ i , the edPR model computes ϕ m , ϕ i , and ϕ e from a consistent framework where ephaptic effects from ϕ e on ϕ m are accounted for (Fig 3C). Due to the electrical coupling between the soma and dendrite, the fluctuations in ϕ m were similar in these compartments, and a more detailed analysis of the AP shapes is found further below. While an action potential essentially gave a depolarization followed by a repolarization of ϕ m , its extracellular signature was essentially a voltage drop (to about -5 mV) followed by a voltage increase (to about +5 mV). This biphasic response of the extracellular AP signature has been seen in several studies (for an analysis, see [20,21]). In experimental recordings, amplitudes in ϕ e fluctuations are typically on the order of 100 μV, which is much smaller than that predicted by the edPR model. The discrepancy is an artifact that is mainly due to the 1D approximation in the edPR model (see Discussion). The dendritic extracellular potential (Fig 3D) was by definition zero at all times, as this compartment was used as the reference point for the potential.
The effect of neuronal firing on the ion concentration dynamics is illustrated in Fig 3E-3H. Before the stimulus onset, the cell was resting at approximately -68 mV, and ion concentrations remained at baseline values. During AP firing, the ion concentrations varied in a jigsawlike fashion in all compartments, except for Ca 2+ , which returned to baseline between each AP and showed notable variation only inside/outside the dendrite since the soma contained no Ca 2+ channels. As the extracellular volume was set to be half as big as the intracellular volume, changes in extracellular ion concentrations were about twice as big as the changes in intracellular ion concentrations. The jigsaw pattern was most pronounced for the K + and Na + concentrations, as these were the main mediators of the APs (Fig 3E-3H). The pattern reflects a cycle of (i) incremental steps away from baseline concentrations, which were mediated by the complex of mechanisms active during the APs, followed by (ii) slower decays back towards baseline, which were mediated by pumps and cotransporters working between the APs. In this simulation, the decay was incomplete, so that concentrations reached gradually larger peak Two plus two compartments (soma + dendrite), with intracellular space to the left and extracellular space to the right. Two kinds of fluxes of different ion species k are involved: transmembrane fluxes (j k,dm , j k,sm ) and intra-and extracellular fluxes (j k,i , j k,e ). The dynamics of the potential ϕ and ion concentration dynamics in all compartments were computed using an electrodiffusive framework, ensuring bulk electroneutrality and a consistent relationship between ion concentrations, electrical charge, and voltages. (B) Active currents were taken from the original PR model [3]. In the soma, these consisted of Na + and K + delayed rectifier currents (I Na and I K-DR ). In the dendrite, these consisted of a voltage-dependent Ca 2+ current (I Ca ), a Ca 2+ -dependent K + current (I K-C ), and a voltage-dependent K + afterhyperpolarization current (I K-AHP ). Ion specific passive (leakage-) currents and homeostatic mechanisms were taken from a previous model by Wei et al. [45], and were identical in the soma and dendrite. These included Na + , K + and Cl − leak currents, a 3Na + /2K + pump (I pump ), a K + /Cl − cotransporter (I KCC2 ), and a Na + /K + /2Cl − cotransporter (I NKCC1 ). In addition, the soma and dendrite included a Ca 2+ /2Na + exchanger (I Ca-dec ), providing an intracellular Ca 2+ decay similar to that in the PR model. values by each consecutive AP. However, as we show later (see Section titled The edPR model predicts homeostatic failure due to high firing rate), the concentrations did, in this case, approach a firing-frequency dependent steady state.
When the firing ceased in Fig 3, the pumps and cotransporters could work uninterruptedly to re-establish the baseline ion concentrations. The resting membrane potential of about -68 mV, was recovered quite rapidly (ms timescale). After this, the slower recovery process of the ion concentration was due to an electroneutral exchange of ions between the neuron and the extracellular space. A full recovery of the baseline concentrations took on the order of 80 s (confirmed by running a longer simulation than the one shown in Fig 3).
As ion concentrations varied during the simulation, so did the ionic reversal potentials, E k (Fig 3I-3J). The by far largest change was seen for the Ca 2+ reversal potential in the dendrite (E k,d ), which dropped by as much as -30 mV during an AP, (i.e., from a baseline value of 124 mV to 94 mV). The explanation is that the basal intracellular Ca 2+ -concentration is extremely low (100 nM) compared to the concentrations of other ion species (several mM), and therefore experienced a much larger relative change during the simulation. Among the main charge carriers (Na + , Cl − , K + ), the lowest concentration is found for K + in the extracellular space ( Table 5 in Methods). For that reason, the second largest change in reversal potential was found for E K , which increased by about 5 mV (i.e., from a basal value of -84 mV to -79 mV) in both the soma and dendrite. The changes in E Ca and E K had a relatively minor impact on the firing pattern in the shown simulations, as the relative change in the driving force ϕ m −E k was not that severe.
The conductivities (σ) of the intra-and extracellular bulk solutions depend on the availability of free charge carriers, and are in the edPR model functions of the ion concentrations and their mobility (cf. Eq 19). The changes in σ were minimal during the conditions simulated here (Fig 3K), i.e., σ varied by a few μS/m over the course of the simulation, while the basal levels were approximately 0.11 S/m and 0.59 S/m for the intra-and extracellular solutions, respectively.
Finally, the 3Na + /2K + pump and Ca 2+ /2Na + exchanger use energy in the form of ATP to move ions against their gradients. The 3Na + /2K + pump exchanges 3 Na + ions for 2 K + ions, and consumes one ATP per cycle [63], while we assumed that the Ca 2+ /2Na + exchanger consumed 1 ATP per cycle (i.e., per Ca 2+ exchanged, as in [64]). As the edPR model explicitly models these processes, we could compute the ATP (energy) consumption of the pumps during the simulation. Fig 3L shows the accumulative number of ATP consumed from the onset of the simulation. The 3Na + /2K + pump was constantly active, as it combated leakage currents and worked to maintain the baseline concentration even before stimulus onset. Before stimulus onset, it consumed ATP at a constant rate (linear curve), which increased only slightly at t = 10 s when the neuron started to fire (small dent in the curve). As the neuron did not contain any passive leakage of Ca 2+ , the Ca 2+ /2Na + exchangers were only active while the neuron was firing. During firing, the Ca 2+ /2Na + exchanger combated the Ca 2+ entering through the dendritic Ca 2+ channels, and then consumed approximately the same amount of energy as the 3Na + /2K + pump (parallel curves). A high metabolic cost of dendritic Ca 2+ spikes has previously been reported also in cortical layer 5 pyramidal neurons [64].
We note that the edPR model had a stable resting state before stimulus onset and that it returned to this resting state after the stimulus had been turned off. In this resting state, ion concentrations remained constant, and ϕ m was approximately -68 mV. This resting equilibrium was due to a balance between the ion-specific leakage channels, pumps, and cotransporters, which we adopted from previous studies (see Methods). However, the existence of a homeostatic equilibrium was not highly sensitive to the choice of model parameters. As we confirmed through a sensitivity analysis, varying membrane parameters (one by one) with ±15% of their default values did not abolish the existence of a stable resting state, but shifted the resting potential by maximally ±3 mV ( Fig 4A) and the resting concentrations by maximally 5% (Fig 4B-4E). Sensitivity of (A) the somatic membrane potential (ϕ sm ) and (B-E) ion concentrations outside the soma to variations of the leak conductances g Na;leak , g K;leak , and g Cl;leak , the ATPase pump strength ρ, and the co-transporter strengths U nkcc1 and U kcc2 . The model was run for 10 s with default parameters. At t = 10 s, selected parameters were changed, one per simulation, by ±15% of their default value. In all cases, the model approached a new steady state during the 3 min simulation, which was not dramatically different from the default steady state. The resting potential was most sensitive to g Na;leak . This was not surprising, as Na + has the reversal potential (57 mV) that is furthest away from the resting potential (� -68 mV), making the driving force (ϕ m −E k ) largest for Na + . All concentration variables were most sensitive either to g Na;leak or ρ. For [Ca 2+ ] se and [Cl − ] se the sensitivity to these parameters were indirect, i.e., through their effects on the resting potential and driving forces. (A-E) Results only shown for somatic compartments, as they were almost identical in the the dendritic compartments. Only extracellular concentrations were shown, but intracellular concentrations followed the same time coarse and intracellular deviances from default values were smaller (due to larger intracellular volume fraction). As we showed in Fig 3L, the Ca 2+ /2Na + exchanger is not active during rest, and it was therefore not included in the sensitivity analysis. https://doi.org/10.1371/journal.pcbi.1007661.g004

The edPR model reproduces the short term firing properties of the original PR model
A motivation behind basing the electrodiffusive (edPR) model on a previously developed (PR) model, was that we wanted to use the firing properties of the original PR model as a "ground truth" when constraining the edPR model. In particular, we wanted the edPR model to qualitatively reproduce the interplay between somatic action potentials and dendritic Ca 2+ spikes, as this was an essential feature of the original PR model [3]. In the PR model, this interplay depended strongly on the coupling strength (coupling conductance) between the soma and dendrite compartment. A weak coupling resulted in a wobbly ping-pong effect, where a somatic AP triggered a dendritic Ca 2+ spike, which in turn fed back to the soma, giving rise to secondary oscillations in ϕ m (Fig 5A). With a strong (about five times stronger) coupling, the somatic and dendritic responses became more similar in shape, as expected ( Fig 5B).
Since the edPR model contained membrane mechanisms and ephaptic effects not present in the PR model, selected parameters in the edPR model had to be re-tuned in order to obtain similar firing as the PR model (see Methods). With the selected parameterization of the edPR model (see the Parameterizations section), we were able to reproduce the characteristic features seen in the PR model for a weak ( Fig 5C) and strong (about five times stronger) coupling between the soma and dendrite ( Fig 5D).

The edPR model predicts homeostatic failure due to high firing rate
As previously discussed, the PR model was, as most existing neuronal models, constructed under the assumption that ion concentration effects are negligible, an assumption that is justified for short term neurodynamics, and for long term dynamics provided that the activity level is sufficiently low for the homeostatic mechanisms to maintain concentrations close to baseline over time. Hence, we expect there to be a scenario (S1) with a moderately low firing rate, where the PR and edPR models can fire continuously and regularly over a long time exhibiting similar firing properties, and another scenario (S2) with a higher firing rate, where the PR and edPR models exhibit similar firing properties initially in the simulation, but where the dynamics of the two models diverge over time due to homeostatic failure accounted for by the edPR model, but not the PR model (which ad hoc assume perfect homeostasis). Simulations of two such scenarios are shown in Figs 6 and 7.
To simulate scenario S1, the PR model ( Fig 6A and 6B) and edPR model ( Fig 6C-6J) were given a constant input (see figure caption) that gave them a firing rate of 1 Hz. Both models settled at a regular firing rate, and in neither of them the firing pattern changed over time, even in simulations of as much as an hour of biological time. For the edPR model, the S1 scenario is the same as that simulated for only a brief period in Fig 3. There, we observed that the ion concentrations gradually changed during the first seconds after the onset of stimulus ( Fig  3E-3H). However, for endured firing, the ion concentrations and reversal potentials settled on a (new) dynamic steady state (Fig 6E-6J), where they deviated by *1-5 mM from the baseline concentrations during rest (i.e., for edPR receiving no input). The apparent "thickness" of the curves (e.g., the red curve for K + in Fig 6H) is due to concentration fluctuations at the short time scale of AP firing. However, after each AP, the homeostatic mechanisms managed to reestablish ionic gradients before the next AP occurred, so that no slow concentration-dependent effect developed in the edPR model at a long time scale.
To simulate scenario S2, the PR model (Fig 7A and 7B) and edPR model (Fig 7C-7J) were given a constant input (see figure caption) that gave them a firing rate of about 3 Hz. The PR model, which included no concentration-dependent effects, settled on a regular firing rate that it could maintain for an arbitrarily long time. Unlike the PR model, the edPR model did not settle at a steady state, but had a firing rate of * 3 Hz only for a period of * 5 s after stimulus onset. During this period, the ion concentrations gradually diverged from the baseline levels ( Fig 7G-7J). The corresponding changes in ionic reversal potentials (Fig 7E and 7F) affected the neuron's firing properties and caused its firing rate to gradually increase before it eventually entered the depolarization block and got stuck at about ϕ m = −30mV. The main explanation behind the altered firing pattern was the change in the K + reversal potential, which, for example, at 9 s after stimulus onset (t = 19 s) had increased by as much as 20 mV from baseline. This shift led to a depolarization of the neuron, which explains both the (gradually) increased firing rate and the (final) depolarization block, i.e., the condition where ϕ m could no longer repolarize to levels below the firing threshold, and AP firing was abolished due to a permanent inactivation of active Na + channels. Neuronal depolarization block is a well-studied phenomenon, which is often caused by high extracellular K + concentrations [65].
The homeostatic failure in S2 was due to the edPR model having a too high firing rate for the ion pumps and cotransporters to maintain ion concentrations close to baseline. The firing rate of 3 Hz was the limiting case (found by trial and error), i.e., for lower firing rates than this, the model could maintain regular firing for an arbitrarily long time. As many neurons can fire at quite high frequencies, a tolerance level of 3 Hz might seem a bit low, and we here provide some comments to this. Firstly, we note that the edPR model could fire at 3 Hz (and gradually higher frequencies) for about 9 s, and could also maintain a higher firing rate than this for a limited time. Secondly, the PR model, and thus the edPR model, represented a hippocampal CA3 neuron, which has been found to have an average firing rate of less than 0.5 Hz [66], so that endured firing of � 3 Hz may be abnormal for these neurons. Thirdly, under biological conditions, glial cells, and in particular astrocytes, provide additional homeostatic functions [67] that were not accounted for in the edPR model, and the inclusion of such functions would probably increase the tolerance level of the neuron. Fourthly, the (3 Hz) tolerance level was a consequence of modeling choices and could be made higher, e.g., by increasing pump rates or compartment volumes. However, we did not do any model tuning in order to increase the tolerance level, as we, in light of the above arguments, considered a 3 Hz tolerance level to be acceptable.

The edPR model predicts homeostatic failure due to impaired homeostatic mechanisms
Above we simulated homeostatic failure occurring because the firing rate became too high for the homeostatic mechanisms to keep up (S2). Homeostatic failure may also occur due to impairment of the homeostatic mechanisms, either due to genetic mutations (see, e.g., [68]) or because the energy supply is reduced, such as after a stroke (see, e.g., [25]). Here, we have used the edPR model to simulate a version of this, i.e., a third scenario (S3) where the ATP-dependent mechanisms, that is, the 3Na + /2K + pumps and the Ca 2+ /2Na + exchangers, were turned off.
In S3, the neuron received no external input, so that the dynamics of the neuron was solely due to gradually dissipating transmembrane ion concentration gradients. After an initial transient, we observed a slow and gradual increase in the membrane potential for about 48 s ( Fig  8A). This coincided with a slow and gradual change in the ion concentrations (Fig 8D-8G) and ionic reversal potentials (Fig 8B and 8C) due to predominantly passive leakage over the membrane.
At about t = 48 s, the membrane potential reached the firing threshold, at which point the active channels started to use what was left of the concentration gradients to generate action potentials and Ca 2+ spikes. This resulted in a burst of activity. During this bursts of activity, the concentration gradients dissipated even faster, since both active and passive channels were then open. As a consequence, the "resting" membrane potential was further depolarized and the neuron went into depolarization block [65]. After this, the neuron continued to "leak" until it settled at a new steady state. The non-zero final equilibrium potential is known as the Donnan equilibrium or the Gibbs-Donnan equilibrium [69]. The reason why the cell did not approach an equilibrium with ϕ m = 0 and identical ion concentrations on both side of the membrane, is that the model contained static residual charges, representing negatively charged macromolecules typically residing in the intracellular environment (see Methods), the sum of which resulted in a final state with a negatively charged inside. In addition, since the Ca 2+ channel inactivated, and since the model had no passive Ca 2+ leakage, Ca 2+ could end up being trapped inside/outside the membrane and did not by necessity approach the Donnan equilibrium, although it was close to it.
As the Ca 2+ dynamics in Fig 8G may seem counterintuitive, we here give some additional explanation of it. During the burst and initial stages of the depolarization block, the dendritic Ca 2+ channels were open. Extracellular Ca 2+ then diffused from the soma towards the dendrite, where it flowed into the neuron. This resulted in a low Ca 2+ concentration in both extracellular compartments and a high Ca 2+ concentration in the intracellular dendritic compartment. The reason why the intracellular Ca 2+ equilibrated more slowly than the extracellular, was that, by assumption, only 1% of the intracellular Ca 2+ concentration was unbuffered and free to diffuse (see Methods), hence, the effective intracellular concentration gradient was a factor 100 lower than it "appears" in Fig 8G. A pattern resembling that in Fig 8A, i.e., a period of silence, followed by a burst of activity, and then silence again, has been seen in experimental EEG recordings of decapitated rats [70], where the activity burst was referred to as "the wave of death", and the phenomenon was ascribed to the lack of energy supply to homeostatic mechanisms. The simulation in Fig 8A represents the single-cell correspondence to this death wave. We note that this phenomenon has been simulated and analyzed thoroughly in a previous modeling study, using a simpler, single compartmental model with ion conservation [40]. We, therefore, do not analyze it further here.

Loss in accuracy when neglecting electrodiffusive effects on concentration dynamics
The concentration-dependent effects studied in the previous subsection were predominantly due to changes in ionic reversal potentials. Effects like this could, therefore, be accounted for by any model that in some way incorporates ion concentration dynamics [27][28][29], provided that the ion concentration dynamics is accurately modeled. As we argued in the Introduction, previous multicompartmental neuron models that do incorporate ion concentration dynamics have not done it in a complete, ion conserving way that ensures a biophysically consistent relationship between ion concentration, electrical charge, and electrical potentials (see, e.g., [27,[48][49][50][51][52][53][54][55][56][57]). To specify, the change in the ion concentration in a given compartment will, in reality, depend on (i) the transmembrane influx of ions into this compartment, (ii) the diffusion of ions between this compartment and its neighboring compartment(s), and (iii) the electrical drift of ions between this compartment and its neighboring compartment(s). Some of the cited models account for only (i) [27,49,51], others account for (i) and (ii) [48,50,[52][53][54][55][56][57], but neither account for (iii). When (iii) is not accounted for, electrical and diffusive processes are implicitly treated as independent processes, a simplifying assumption which is also incorporated in the reaction-diffusion module [71] in the NEURON simulation environment [72]. In models that apply this assumption, there will therefore be drift currents (along axons and dendrites) that affect ϕ m (through the cable equation), but not the ion concentration dynamics, although they should, since also the drift currents are mediated by ions.
Here, we use simulations on the edPR model to test the inaccuracy introduced when not accounting for the effect of drift currents on ion concentration dynamics. We do so by comparing how many ions that were transferred from the somatic to the dendritic compartment through the intracellular ( Fig 9A) and extracellular (Fig 9B) space, due to ionic diffusion (orange curves) versus electrical drift (blue curves), throughout the simulation in Fig 3. We note that Fig 9 shows the accumulatively moved number of ions (from time zero to t) due to axial fluxes exclusively. That is, the large number of, for example, Na + ions transported intracellularly from the dendrite to the soma (negative sign) in Fig 9A1, does not by necessity mean that Na + ions were piling up in the soma compartment, as the membrane efflux of Na + was not accounted for in the figure.
Although diffusion tended to dominate the intracellular transport of ions on the long time scale (Fig 9A1-9A4), the transport due to electrical drift was not vanishingly small. For example, the number of K + and Cl − ions transported by electrical drift was at the end of the stimulus period (t = 20 s) about 35% of the transport due to diffusion for both species. In the extracellular space, diffusion was the clearly dominant transporter of Na + and K + (Fig 9B1 and 9B2), while diffusion and electrical drift were of comparable magnitude for the other ion species (Fig  9B3-9A4). Of course, these estimates are all specific to the edPR model, as they will depend strongly on the included ion channels, ion pumps and cotransporters, and on how they are distributed between the soma and dendrite. In general, however, the simulations in Fig 9 suggest that electrical drift is likely to have a non-negligible effect on ion concentration dynamics, and that ignoring this effect will give rise to rather inaccurate estimates.
Finally, we also converted the sum of ionic fluxes in Fig 9 into an effective current, represented as the number of transported unit charges, e + (Fig 9A5-9B5). Interestingly, diffusion and drift contributed almost equally to the axial charge transport in the system. We note, however, that the movement of charges per time unit is indicated by the slope of the curves, which was much larger for the drift case (blue curve) than for diffusion (orange curve). The drift curve had a jigsaw shape, which shows that drift was moving charges back and forth in the system, while the diffusion always went in the same direction, explaining why it, despite being smaller than the drift current, had a comparably large accumulative effect on charge transport. The temporally averaged picture of charge transport that emerges from Fig 9A5 is that of a slow current loop where charge is transferred intracellularly from the soma to the dendrite (Fig 9A5), where it crosses the membrane (outward current), and then is transferred extracellularly back from the dendrite to the soma (Fig 9B5), before crossing the membrane again (inward current). This configuration is similar to the slow loop current seen during spatial buffering by astrocytes (see, e.g. Fig 1 in [67]).

Loss in accuracy when neglecting electrodiffusive effects on voltage dynamics
In the previous section, we investigated the consequences of neglecting (iii) the contribution of drift currents on ion concentration dynamics. Here, we investigate the consequences of neglecting the effect of ionic diffusion (along dendrites) on the electrical potential, focusing on the extracellular potential ϕ e . Forward modeling of extracellular potentials is typically based on volume conductor (VC) theory [16-18, 20, 21], which assumes that diffusive effects on electrical potentials are negligible. Being based on a unified electrodiffusive KNP framework (Fig 1), the edPR model accounts for the effects of ionic diffusion on the electrical potentials, and can thus be used to address the validity of this assumption.
To illustrate the effect of diffusion on ϕ e , we may split it into a component ϕ VC,e explained by standard VC-theory, and a component ϕ diff,e representing the additional contribution caused by diffusive currents (Eq 81). In the simulation in Fig 3, the diffusive contribution was found to be very small compared to the VC-component (Fig 10). However, while ϕ VC,e fluctuated rapidly from negative to positive values during neuronal activity, ϕ diff,e varied on a slower time scale and had the same directionality throughout the simulation. This is equivalent to what we saw in Fig 9B5, i.e., that diffusion always moved charge in the same direction. Moreover, if we take the temporal averages of the potentials over the time series in Fig 10A, we find that they are -0.0023 mV, 0.0037 mV, and -0.0060 mV for ϕ e , ϕ diff,e , and ϕ VC,e , respectively. This shows that the average diffusion-and VC-components of the total potential were of the same order of magnitude. As we also have demonstrated in previous studies, diffusion is thus likely to be important for the low-frequency components of extracellular potentials [31, 32, 73,74]. Albeit small, the slowly varying diffusion evoked shifts in ϕ e are putatively important for explaining the direct-current (DC) like shifts and long-time concentration dynamics reported during, e.g., spreading depression [25,26].

Discussion
The original Pinsky-Rinzel (PR) is a reduced model of a hippocampal neuron, which reproduces the essential somatodendritic firing properties of CA3 neurons despite having only two compartments [3]. Simplified neuron models like that are useful, partly because their reduced complexity makes them easier to analyze, and as such, can lead to insight in key neuronal mechanisms, and partly because they demand less computer power and can be used as modules in large scale network simulations. Whereas the PR model, as most available neuron models, assumes that ion concentrations remain constant during the simulated period, the electrodiffusive Pinsky-Rinzel (edPR) model proposed here models ion concentration dynamics explicitly. The edPR model may thus be seen as a supplement to the PR model, which should be applied to simulate conditions where ion concentrations are expected to vary with time.
In the results section, we showed that the edPR model closely reproduced the firing properties of the PR model for short term dynamics (Fig 5), and for long term dynamics provided that the firing rate was sufficiently low for the homeostatic mechanisms to maintain ion concentrations close to baseline (Fig 6). We also showed that if the firing rate became too high (Fig 7), or if the homeostatic mechanisms were impaired (Fig 8), unsuccessful homeostasis would cause ion concentrations to gradually shift over time, and lead to slowly developing changes in the firing properties of the edPR model, changes that were not accounted for by the original PR model. The edPR model was based on an electrodiffusive framework [60], which ensured a consistent relationship between ion concentrations, electrical charge, and electrical potential in four compartments. To our knowledge, the edPR model is the first multicompartmental neuronal model that ensures complete and consistent ion concentration and charge conservation.

Model assumptions
The construction of the edPR model naturally involved making a set of modeling choices, and the most important of these are discussed here. Firstly, in the construction of the model, we focused on morphological simplicity, biophysical rigor, and mechanistic understanding, rather than on replicating any specific biological scenario and incorporating biological details. Secondly, simultaneous data of variations in all intra-and extracellular concentrations during neuronal firing are not available, and it might not even be feasible to obtain such data. Consequently, computational modeling based on biophysical constraints may be the best means to estimate it. The concentration dynamics in the edPR model were thus not directly constrained to data but constrained so that there was, at all times, an internally consistent relationship between all ion concentrations and all electrical potentials, ensuring an electroneutral bulk solution. Thirdly, to include extracellular dynamics to models of neurons or networks of such is computationally challenging, since the extracellular space, in reality, is an un-confined three-dimensional continuum, locally affected by populations of nearby neurons and glial cells. As we wanted to keep things simple and conceptual, we chose to use closed boundary conditions, i.e., no ions and no charge were allowed to leave or enter the system consisting of the single (2-compartment) neuron and its local and confined (2-compartment) surrounding (Fig 2). Tecnically, it would be straightforward to increase the number of compartments (i.e., the spatial resolution) in the model.
A consequence of using closed boundary conditions was that the extracellular (like the intracellular) currents became one-dimensional (from soma to dendrite), while in reality, extracellular currents pass through a three-dimensional volume conductor. The edPR model could be made three dimensional if embedded in a bi-or tri-domain model (as discussed below). However, currently, it is 1D, and the effect of the 1D assumption was essentially an increase in the total resistance (fewer degrees of freedom) for extracellular currents, which gave rise to an artificially high amplitude in extracellular AP signatures (Fig 3). We note, however, that the closed boundary is actually equivalent to assuming periodic boundary conditions, so that the edPR model essentially simulates the hypothetical case of a population of perfectly synchronized neurons, i.e., one where all neurons are doing exactly the same as the simulated neuron, so that no spatial variation occurs. Likely, this may give accurate predictions for ion concentration shifts over time, as these reflect a temporal average of activity, but less accurate predictions for brief and unique electrical events, such as action potentials, which are not likely to be elicited in perfect synchrony by a large population [31].
Fourthly, to faithfully represent a morphologically complex neuron with a reduced number of compartments is a non-trivial task. Available analytical theory for collapsing branching dendrites into equivalent cylinders are generally based on certain assumptions about branching symmetries, and on preserving electrotonic distances [75]. However, it is unlikely that the length constants of electrodynamics and ion concentration dynamics scale in the same way. Hence, in the edPR model, the volumes and membrane areas of, and cross-section areas between, the two neuronal compartments were here introduced as rather arbitrary model choices, fixed at values that were verified to give agreement between the firing properties of the edPR model and the PR model.

Outlook
Being applicable to simulate conditions with failed homeostasis, the edPR model opens up for simulating a range of pathological conditions, such as spreading depression or epilepsy [22][23][24][25], which are associated with large scale shifts in extracellular ion concentrations. A particular context in which we anticipate the edPR model to be useful is that of simulating spreading depression. Previous spatial, electrodiffusive, and biophysically consistent models of spreading depression have targeted the problem at a large-scale tissue-level, using a mean-field approach [30, 76,77]. These models were inspired by the bi-domain model [78], which has been successfully applied in simulations of cardiac tissue [79,80]. The bi-domain model is a coarse-grained model, in which the tissue is considered as a bi-phasic continuum consisting of an intracellular and extracellular domain. That is, a set of intra-and extracellular variables (i.e., voltages and ion concentrations), and the ionic exchange between the intra-and extracellular domains, are defined at each point in space. This simplification allows for large scale simulations of signals that propagate through tissue but sacrifices morphological detail. In the context of spreading depression, a shortcoming with this simplification is that the leading edge of the spreading depression wave in both the hippocampus and cortex is in the layers containing the apical dendrites [22]. This suggests that the different expression of membrane mechanisms in deeper (somatic) and higher (dendritic) layers may be crucial for fully understanding the propagation and genesis of the wave. In this context, the edPR model could enter as a module in a, let us say, bi-times-two-domain model, where each point in (xy) space contains a set of (i) somatic intracellular variables, (ii) somatic extracellular variables, (iii) dendritic intracellular variables, and (iv) dendritic extracellular variables, and thus accounts for the differences between the higher and lower layers. We should note that the state of the art models of spreading depression are not bi-domain models but rather tri-domain models, as they also include a glial domain to account especially for the work done by astrocytes in K + buffering [30, 76,77]. Hence, to use the edPR model to expand the current spreading depression models, a natural first step would be to include a glial (astrocytic) compartment in it, so that it eventually could be implemented as a tri-times-two-domain model.

The Kirchoff-Nernst-Planck (KNP) framework
In the following section, we derive the KNP continuity equations for a one-dimensional system containing two plus two compartments (Fig 2A), with sealed boundary conditions (i.e., no ions can enter or leave the system). The geometrical parameters used in the edPR model were as defined in Table 1. Since typical neuronal/extracellular/glial volume fractions in neuronal tissue are 0.4/0.2/0.4 [82], we let the extracellular space be half as voluminous as the intracellular neuronal space.
Two kinds of fluxes are involved: transmembrane fluxes and intra-and extracellular fluxes. The framework is general to the choice of the transmembrane fluxes. A transmembrane flux of ion species k (j k,m ) represents the sum of all fluxes through all membrane mechanisms that allow ion k to cross the membrane. Intracellular flux densities are described by the Nernst-Planck equation: In Eq 1, 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 ER, λ i is the tortuosity, which represents the slowing down of diffusion due to obstacles, γ k ([k] di −[k] si )/Δx is the axial 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� i is the average concentration, that is, γ k ([k] di + [k] si )/2, and (ϕ di −ϕ si )/ Δx is the axial potential gradient. Similarly, the extracellular flux densities are described by In Eq 2, we do not include γ k , as all ions can move freely in the extracellular space. Diffusion constants, tortuosities, and intracellular fractions of freely moving ions used in the edPR model were as in Table 2.
Ion conservation. The KNP framework is based on the constraint of ion conservation. To keep track of ion concentrations we solve four differential equations for each ion species k: For each compartment, all flux densities are multiplied by the area they go through and divided by the volume they enter to calculate the change in ion concentration. If we insert the Nernst-Planck equation (Eq 1) for the intracellular flux density, the first of these equations becomes: where the voltage variables so far are undefined. Four constraints to derive ϕ. If we have four ion species (Na + , K + , Cl − , and Ca 2+ ) in four compartments, we have 20 unknown parameters (16 for [k] and four for ϕ), while Eqs 3-6 for four ion species give us only 16 equations. To solve this, we need to define additional constraints that allow us to express the potentials ϕ in terms of ion concentrations.
As we may define an arbitrary reference point for ϕ, we take as our first constraint, i.e., (i) the potential outside the dendrite is defined to be zero. The second constraint is that (ii) the membrane is a parallel plate capacitor that always separates a charge Q on one side from an opposite charge −Q on the other side, giving rise to a voltage difference Here, C m is the total capacitance of the membrane, i.e., C m = c m A m , where c m is the more commonly used capacitance per membrane area. As, by definition, ϕ m = ϕ i − ϕ e , we get: in the dendrite (since ϕ de = 0), and in the soma. The third constraint is that (iii) we assume bulk electroneutrality. This means that the net charge associated with the ion concentrations in a given compartment by constraint must be identical to the membrane charge in this compartment. The intracellular dendritic charge is By inserting this into Eq 10, we obtain the final expression for ϕ di : By inserting the equivalent expression for Q si into Eq 11, we get Here, the extracellular potential is not set to zero, so we need a fourth constraint to determine ϕ si and ϕ se separately. The fourth and final constraint is that (iv) we must ensure charge anti-symmetry. For the charge anti-symmetry between the two sides of the capacitive membrane (Q i = −Q e ) to be preserved in time, we must define our initial conditions so that this is the case at t = 0, and the system dynamics so that this stays the case. Hence, the system dynamics must ensure that dQ di /dt = −dQ de /dt and dQ si /dt = −dQ se /dt. The membrane currents (in isolation) will always fulfill this criterion, as any charge that crosses the membrane by definition disappears from one side of it and pops up at the other. Hence, we thus need to make sure that also the axial currents (in isolation) fulfill the criterion. The system must thus be constrained so that, if an extracellular current transports a charge δq into a given extracellular compartment, the intracellular current must transport the opposite charge −δq into the adjoint intracellular compartment. That is, we must have that: where i i and i e are the intra-and extracellular current densities, respectively. To find an expression for these, we multiply Eqs 1 and 2 by Fz k and sum over all ion species k. The expressions for the intra-and extracellular current densities then become: In Eq 15, the first term is the diffusion current density: which is defined by the ion concentrations. The second term is the field driven current density where we have identified the conductivity as Similarly, Eq 16 can be written in terms of i diff,e , i field,e , and σ e . By combining Eqs 14, 15 and 16, we obtain: In Eq 20, ϕ di and ϕ de are already known from Eqs 8 and 12, while i diff and σ are expressed in terms of ion concentrations. We may thus solve Eqs 13 and 20 for the last two voltage variables ϕ se and ϕ si : Membrane mechanics Leakage channels. In the original PR model, the membrane leak current represents the combined contribution from all ion species. When using the KNP framework, on the other hand, where we keep track of all ions separately, the leak current must be ion-specific. We modeled this as in [45], that is, for each ion species k, we implemented a passive flux density across the membrane where g k;leak 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. The reversal potential is a function of ion concentrations, and is calculated using the Nernst equation: Here, R is the gas constant, T is the absolute temperature, γ k is the intracellular fraction of free ions, and [k] e and [k] i are the concentrations of ion k outside and inside the cell, respectively. We included Na + , K + , and Cl − leak currents in both compartments.
Active ion channels. All active ion channel currents were adopted from the original PR model [3], as they were described in [8], and converted to ion channel fluxes. The soma compartment contained a Na + flux (j Na ) and a K + delayed rectifier flux (j K−DR ), while the dendrite contained a voltage-dependent Ca 2+ flux (j Ca ), a voltage-dependent K + AHP flux (j K−AHP ), and a Ca 2+ -dependent K + flux (j K−C ): The voltage-dependent conductances were modeled using the Hodkin-Huxley formalism with differential equations for the gating variables: The conductances and gating variables were given by: g Ca ¼ g Ca s 2 z; ð34Þ Thirdly, we modified the voltage-dependent Ca 2+ current to include an inactivation variable z (Eqs 31 and 34). We implemented this inactivation like they did in [84] (Eqs A2-A3), but set the time constant to 1 s, the half-activation voltage to -30 mV, and the slope of the steady-state Boltzmann fit to z 1 to 0.001. In the original PR model, inactivation was neglected due to the argument that it was too slow to have an impact on simulation outcomes [2]. However, in our simulations, we observed that it had a significant impact, and therefore we included it. Homeostatic mechanisms. To maintain baseline ion concentrations for low frequency activity we added 3Na + /2K + pumps, K + /Cl − cotransporters (KCC2), and Na + /K + /2Cl − cotransporters (NKCC1). Their functional forms were taken from [45].
where ρ, U kcc2 , and U nkcc1 are pump and cotransporter strengths. We assumed optimal pump functionality and set ρ to be the pump strength used in [45] for the fully oxygenated state with normal bath potassium (ρ max ). Intracellular Ca 2+ decay was modeled in a similar fashion as in [3], but to ensure ion conservation we modeled it as an electroneutral Ca 2+ /2Na + exchanger, exchanging one Ca 2+ (outward) for two Na + (inward). Putatively, this phenomenological model for the decay could represent the joint effect of several mechanisms in a real system, such as the Ca 2+ /3Na + exchanger, a Ca 2+ leakage current, SERCA pumps, etc. The decay flux density was defined as: where U Ca−dec is the decay rate, and [Ca 2+ ] i,b is the basal Ca 2+ concentration, set to 0.01 mM. Model summary. We summarize the model here for easy reference. In short, we solved four differential equations for all ion species k: At each time step, ϕ in all four compartments was derived algebraically: The total membrane flux densities were as follows: j Na;dm ¼ j Na;leak þ 3j pump þ j nkcc1 À 2j CaÀ dec ; ð71Þ

Original Pinsky-Rinzel model
We implemented the original Pinsky-Rinzel equations from Box 8.1 in [8]. The reversal potential of the leak current, not specified in [8], was set to -68 mV to ensure a resting potential close to that of the edPR model. We also used this as the initial potentials, that is, ϕ sm,0 =  Table 3, we get an absolute temperature of 309.14 K, corresponding to a body temperature of 36˚C.
https://doi.org/10.1371/journal.pcbi.1007661.t003 associated with the initial ion concentrations: ½X À � i;0 ¼ z Na ½Na þ � i;0 þ z K ½K þ � i;0 þ z Cl ½Cl À � i;0 þ z Ca ½Ca 2þ � i;0 À � m;0 c m A m V i F ; ð75Þ ½X À � e;0 ¼ z Na ½Na þ � e;0 þ z K ½K þ � e;0 þ z Cl ½Cl À � e;0 þ z Ca ½Ca 2þ � e;0 þ � m;0 In the next step, we calibrated the model by running it for 1800 s to obtain the post-calibrated values of all the state variables (Table 5). These post-calibrated values were written to file and used as initial conditions in all simulations shown throughout this paper. For technical reasons, we did not read the constant residual concentrations, [X − ] i,0 and [X − ] e,0 , to/from file, but re-calculated them from Eqs 75 and 76 in the beginning of each simulation to minimize rounding errors and ensure strict electroneutrality. While the pre-calibrated initial conditions were identical in the somatic and dendritic compartment, the post-calibration values were not strictly identical, but identical up to the decimal place included in Table 5. Hence, the indicated values apply for both the soma and dendrite compartments. The post-calibrated values of the ion concentrations gave us the following reversal potentials: E Na = 57mV, E K = −84mV, E Cl = −79mV, and E Ca = 124mV.
The pre-calibration values for the ion concentrations were taken from Table 2.1 in [85], which lists the ranges of typical intra-and extracellular concentrations in mammalian neurons. From the ranges given in this table, we selected values that made the edPR model (throughout the calibration period) reside close to the selected initial membrane potential of -68 mV, a typical value found for CA3 pyramidal cells in adult rats [86]. Stimulus current. We stimulated the cell by injecting a K + current i stim into the soma. Previous computational modeling of a cardiac cell has shown that stimulus with K + causes the least physiological disruption [33]. To ensure ion conservation, we removed the same amount of K + ions from the corresponding extracellular compartment: Analysis. Fig 9: To calculate the accumulative transport of ion species k in the intracellular solution (from time zero to t) due to diffusion, we integrated A i N A j k,diff, i from time zero to t, where N A is the Avogadro constant. Similarly, we integrated A e N A j k,diff,e to calculate the accumulative transport of ions in the extracellular solution due to diffusion. We did the same calculations with j k,drift to study the accumulative transport of ions due to drift. When knowing the accumulative transport of each ion species, k akkum , we calculated the total transport of e + from their weighted sum: where the last equality follows when we insert Eq 18 for the extracellular field-driven current density i field,e , and use that ϕ de = 0. As in [32], we may split ϕ se into two components: where ϕ VC,se is the potential as it would be predicted from standard volume conductor (VC) theory [20,21], and ϕ diff,se is the additional contribution from diffusion [32]. With this, Eq 80 can be written: We may split this into two equations if we recognize that is the standard formula used in VC theory, which is based on the assumption that the extracellular current is exclusively due to a drop in the extracellular VC-potential ϕ VC,se . The remainder of Eq 82 then leaves us with i diff;e ¼ À s e � diff;se Dx : ð84Þ Since we already knew i e and i diff,e from simulations with the KNP framework, we used Eqs 83 and 84 to calculate ϕ VC,se and ϕ diff,se separately.
Numerical implementation. We implemented the differential equations in Python 3.6 and solved them using the solve_ivp function from SciPy. We used its default Runge-Kutta method of order 5(4), and set the maximal allowed step size to 10 −4 . The code is made available at https://github.com/CINPLA/EDPRmodel and https://github.com/CINPLA/ EDPRmodel_analysis.