PLoS Computational Biology
Public Library of Science
An electrodiffusive, ion conserving Pinsky-Rinzel model with homeostatic mechanisms
DOI 10.1371/journal.pcbi.1007661, Volume: 16, Issue: 4,

Table of Contents




Neurons generate their electrical signals by letting ions pass through their membranes. Despite this fact, most models of neurons apply the simplifying assumption that ion concentrations remain effectively constant during neural activity. This assumption is often quite good, as neurons contain a set of homeostatic mechanisms that make sure that ion concentrations vary quite little under normal circumstances. However, under some conditions, these mechanisms can fail, and ion concentrations can vary quite dramatically. Standard models are thus not able to simulate such conditions. Here, we present what to our knowledge is the first multicompartmental neuron model that accounts for ion concentration variations in a way that ensures complete and consistent ion concentration and charge conservation. In this work, we use the model to explore under which activity conditions the ion concentration variations become important for predicting the neurodynamics. We expect the model to be of great value for the field of neuroscience, as it can be used to simulate a range of pathological conditions, such as spreading depression or epilepsy, which are associated with large changes in extracellular ion concentrations.

Sætra, Einevoll, Halnes, and Lytton: An electrodiffusive, ion conserving Pinsky-Rinzel model with homeostatic mechanisms


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., [27]).

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 Ileak = gleak(ϕmEleak ), 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., [1012]). 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., [1315]). 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., [1618]), 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 ϕm that the (ephaptic) effect of ϕe on neurodynamics can be neglected without severe loss in accuracy [19].

Modeling intra- and extracellular dynamics: Standard theory vs. unified framework.
Fig 1
(A) The dynamics of the membrane potential (ϕm) and transmembrane currents of neurons are typically modeled using cable theory. It is then assumed that the extracellular environment is grounded (ϕe = 0). Typically, it is also assumed that ion concentrations both in the intra- and extracellular space are constant, so that also ionic reversal potentials remain constant. (B) When knowing the transmembrane neuronal currents (as computed in (A) ), standard volume conductor theory [20, 21] allows us to estimate the extracellular potential, which is computed as the sum of neuronal point-current sources weighted by their distance to the recording location. An underlying assumption is that fluctuations in ϕe (as computed in (B)) are so small that they have no effect on the neurodynamics (as computed in (A)), i.e., there is no ephaptic coupling. Another underlying assumption (cf. constant ion concentrations) is that extracellular diffusive currents do not affect electrical potentials. (C) We propose a unified, electrodiffusive framework for intra- and extracellular ion concentration and voltage dynamics, assuring a consistent relationship between ion concentrations, electrical charge, and electrical potential in all compartments.Modeling intra- and extracellular dynamics: Standard theory vs. unified framework.

There are, however, scenarios where the assumptions of constant ion concentrations and a grounded extracellular space are not justifiable. Notably, large-scale extracellular ion concentration changes are a trademark of several pathological conditions, including epilepsy and spreading depression [2225]. In these cases, neurons are unable to maintain their baseline conditions because they for various reasons are too active and/or their homeostatic mechanisms are too slow. During spreading depression, the extracellular K+ concentration can change from a baseline value of about 3-5 mM to pathological levels of several tens of mM, and the increased K+ concentration tends to coincide with a slow, direct-current (DC) like drop in the extracellular potential, which may be several tens of millivolts in amplitude [25, 26], and can give rise to large spatial gradients. For example, one experiment saw the extracellular K+-concentration and ϕe vary by as much as 30 mM and 20 mV, respectively, over the hippocampal depth [26]. Such dramatic gradients in the extracellular environment are likely to have a strong impact on the dynamical properties of neurons, both through the concentration-dependent changes in ion-channel reversal potentials [2729] and putatively through a direct ephaptic effect from ϕe on the membrane potential.

The construction of accurate neuron models that include ion concentration dynamics (and conservation) poses two key challenges. Firstly, ion conserving models need a finely adjusted balance between the homeostatic machinery and all passive and active ion-specific currents so that all ion concentrations, as well as voltages, vary in a biophysically realistic way over time when the neuron is active. Secondly, in spatially extended models, ions will not move only across membranes, but also within the extracellular and intracellular space. Such ionic movement may be propelled both by diffusion and electrical drift. Ionic diffusion can, in principle, affect the electrical potential (since ions carry charge), and the electrical potential can, in principle, affect ion concentration dynamics (since ions drift along potential gradients) [3032]. 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, 3347]. 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, 4857]), 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 [5862].

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+, Ca2+, 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 [2225].


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 (INa and IK−DR), while the dendritic compartment contains a voltage-dependent Ca2+ current (ICa), a voltage-dependent K+ afterhyperpolarization current (IK−AHP), and a Ca2+-dependent K+ current (IK−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 (Ipump), K+/Cl cotransporters (IKCC2), Na+/K+/2Cl cotransporters (INKCC1), and Ca2+/2Na+ exchangers (ICa−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 Fig 2 and described in details in the Methods section.

edPR model architecture.
Fig 2
(A) 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 (jk,dm, jk,sm) and intra- and extracellular fluxes (jk,i, jk,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 (INa and IK-DR). In the dendrite, these consisted of a voltage-dependent Ca2+ current (ICa), a Ca2+-dependent K+ current (IK-C), and a voltage-dependent K+ afterhyperpolarization current (IK-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 (Ipump), a K+/Cl cotransporter (IKCC2), and a Na+/K+/2Cl cotransporter (INKCC1). In addition, the soma and dendrite included a Ca2+/2Na+ exchanger (ICa-dec), providing an intracellular Ca2+ decay similar to that in the PR model.edPR model architecture.

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 Ca2+/2Na+ exchangers.

Output of the edPR model.
Fig 3
A 27 pA step-current injection was applied to the somatic compartment between t = 10 s and t = 20 s, and the model responded with a firing rate of 1 Hz. (A-B) The membrane potential ϕm of the soma and the dendrite, respectively. (C-D) The extracellular (index e) potential ϕe of the soma (index s) and the dendrite (index d), respectively. The dendritic extracellular compartment was chosen as the reference point when calculating potentials, so ϕde was zero by definition. Since amplitudes in ϕm were so much larger than for ϕe, intracellular (index i) potentials (ϕi = ϕe+ ϕm) were similar to ϕm, and therefore not shown. (E-H) Ion concentrations dynamics of all ion species k (Na+, Cl, K+, Ca2+) in all four compartments shown in terms of their deviance from baseline concentrations. (I-J) Changes in reversal potentials for all ion species in the soma and the dendrite, respectively. (K) Change in conductivity of the intra- and extracellular media (σi and σe, respectively). (L) Accumulative number of ATP molecules consumed by the 3Na+/2K+ pumps and Ca2+/2Na+ exchangers.Output of the edPR model.

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 jigsaw-like fashion in all compartments, except for Ca2+, which returned to baseline between each AP and showed notable variation only inside/outside the dendrite since the soma contained no Ca2+ 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 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, Ek (Fig 3I–3J). The by far largest change was seen for the Ca2+ reversal potential in the dendrite (Ek,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 Ca2+-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 EK, 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 ECa and EK had a relatively minor impact on the firing pattern in the shown simulations, as the relative change in the driving force ϕmEk 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 Ca2+/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 Ca2+/2Na+ exchanger consumed 1 ATP per cycle (i.e., per Ca2+ 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 Ca2+, the Ca2+/2Na+ exchangers were only active while the neuron was firing. During firing, the Ca2+/2Na+ exchanger combated the Ca2+ entering through the dendritic Ca2+ channels, and then consumed approximately the same amount of energy as the 3Na+/2K+ pump (parallel curves). A high metabolic cost of dendritic Ca2+ 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 analysis.
Fig 4
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 Unkcc1 and Ukcc2. 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 (ϕmEk) largest for Na+. All concentration variables were most sensitive either to g¯Na,leak or ρ. For [Ca2+]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 Ca2+/2Na+ exchanger is not active during rest, and it was therefore not included in the sensitivity analysis.Sensitivity analysis.

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 Ca2+ 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 Ca2+ 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).

Short term dynamics of the PR and edPR models.
Fig 5
The original PR model (top row) and the edPR model (bottom row) exhibit the same spike shape characteristics. (A) Spike shape in PR model for weak coupling (low coupling conductance) between the soma and the dendrite. (B) Spike shape in PR model for strong coupling (high intracellular conductivity) between the soma and the dendrite. (C) Spike shape in edPR model for weak coupling between the soma and the dendrite. (D) Spike shape in edPR model for strong coupling between the soma and the dendrite. (A-D) A step-stimulus current was turned on at t = 10 s, with stimulus strength being 1.35 μA/cm2 in (A), 0.78 μA/cm2 in (B), 31 pA in (C), and 27 pA in (D). The panels show snapshots of a selected spike. See the Parameterizations section in Methods for a full description of the parameters used.Short term dynamics of the PR and edPR models.

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.

Model comparison for scenario with low frequency firing.
Fig 6
Simulations on the PR model and edPR model when both models are driven by a constant input, giving them a firing rate of about 1 Hz. Simulations covered one hour (3600 s) of biological time. (A-D) A 10 s sample of the dynamics of the somatic membrane potential ϕsm and dendritic (free) Ca2+ concentration in the PR model (A-B) and edPR model (C-D). This regular firing pattern was sustained over the full 3600 s simulation in both models (inset panels). (D) Of the total amount of intracellular Ca2+, only 1% (as plotted) was assumed to be free (unbuffered). (E-F) Ionic reversal potentials and (G-J) ion concentrations in the edPR model did not vary on a long time scale. Indices i, e, s, and d indicate intracellular, extracellular, soma, and dendrite, respectively. (A-J) Stimulus onset was t = 10 s in both models, and stimulus strength was istim = 0.78μA/cm2 in the PR model (A-B) and istim = 27pA in the edPR model (C-J) . See the Parameterizations section in Methods for a full description of the parameters used.Model comparison for scenario with low frequency firing.
Model comparison for scenario with high frequency firing.
Fig 7
Simulations on the PR model and edPR model when both models are driven by a constant input, giving them a firing rate of about 3 Hz. Simulations covered 200 s of biological time. (A-D) A 12 s sample of the dynamics of the somatic membrane potential ϕsm and dendritic (free) Ca2+ concentration in the PR model (A-B) and edPR model (C-D). The regular firing pattern in the PR model (A-B) was sustained over the full 200 s simulation (inset panels), while the edPR model stopped firing and entered depolarization block around t = 20 s. (D) Of the total amount of intracellular Ca2+, only 1% (as plotted) was assumed to be free (unbuffered). (E-F) Ionic reversal potentials and (G-J) ion concentrations in the edPR model varied throughout the simulation, and gradually diverged from baseline conditions. Indices i, e, s, and d indicate intracellular, extracellular, soma, and dendrite, respectively. Main panels show 12 s samples of the ion concentration dynamics, while insets show the dynamics over the full 200 s simulations. (A-J) Stimulus onset was t = 10 s in both models, and stimulus strength was istim = 1.55μA/cm2 in the PR model (A-B) and istim = 48pA in the edPR model (C-J) . See the Parameterizations section in Methods for a full description of the parameters used.Model comparison for scenario with high frequency firing.

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 re-establish 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 Ca2+/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.

The wave of death.
Fig 8
Simulation on the edPR model when the 3Na+/2K+ pumps and the Ca2+/2Na+ exchangers were turned off. The model received no external stimulus. The simulation covered 10 minutes of biological time. (A) A 60 s sample of the dynamics of the somatic membrane potential ϕsm. Inset shows a close-up of the burst of activity occurring at about t = 48 s. (B-C) Reversal potentials in the soma (B) and dendrite (C). (D-G) Ion concentrations in all four compartments. Somatic and dendritic concentrations were almost identical for all ion species except for Ca2+. Indices i, e, s, and d indicate intracellular, extracellular, soma, and dendrite , respectively. See the Parameterizations section in Methods for a full description of the parameters used.The wave of death.

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 Ca2+ 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 Ca2+ channel inactivated, and since the model had no passive Ca2+ leakage, Ca2+ 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 Ca2+ 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 Ca2+ channels were open. Extracellular Ca2+ then diffused from the soma towards the dendrite, where it flowed into the neuron. This resulted in a low Ca2+ concentration in both extracellular compartments and a high Ca2+ concentration in the intracellular dendritic compartment. The reason why the intracellular Ca2+ equilibrated more slowly than the extracellular, was that, by assumption, only 1% of the intracellular Ca2+ 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 [2729, 3357], 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, 4857]). 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, 5257], 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.

Axial transport of ions and charge due to drift versus diffusion.
Fig 9
(A1-A4) The number of ions transported intracellularly from soma to dendrite from time zero to t by electrical drift versus ionic diffusion. (B1-B4) The number of ions transported extracellularly from (outside) soma to (outside) dendrite from time zero to t. (A5) Net charge transported intracellularly from soma to dendrite, represented as the number of unit charges e+. (B5) Net charge transported extracellularly from soma to dendrite, represented as the number of unit charges e+. (A-B) The simulation was the same as in Fig 3. See the Analysis section in Methods for a description of how we did the calculations.Axial transport of ions and charge due to drift versus diffusion.

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 [1618, 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].

Effect of diffusion on extracellular potential.
Fig 10
The extracellular potential ϕe in the edPR model, split (cf. Eq 81) into a component explained by standard VC-theory (ϕVC,e) and a “correction” (ϕdiff,e) when diffusive contributions are accounted for. (A-B) The simulation was the same as in Fig 3. (B) Close-up of selected AP in (A) . See the Analysis section in Methods for a description of how we calculated ϕVC,e and ϕdiff,e.Effect of diffusion on extracellular potential.


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.


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 [2225], 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.

Table 1
Geometrical parameters.
Δx (distance between the two compartments)667⋅10−6 m
As (somatic membrane area)616⋅10−12 m2 *
Ad (dendritic membrane area)616⋅10−12 m2 *
Ai (intracellular cross-section area)α · As
Ae (extracellular cross-section area)Ai/2
Vsi (somatic intracellular volume)1437⋅10−18 m3 *
Vdi (dendritic intracellular volume)1437⋅10−18 m3 *
Vse (somatic extracellular volume)718.5⋅10−18 m3 *
Vde (dendritic extracellular volume)718.5⋅10−18 m3 *
* The intracellular volumes (Vsi, Vdi) and membrane areas (As, Ad) correspond to spheres with radius 7 μm.
The parameter α describes the coupling strength of the model and is defined in the Parameterizations section. Its default value was 2.

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 (jk,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, Dk 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, zk 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.

Table 2
Diffusion constants, tortuosities, and intracellular fractions of freely moving ions.
DNa (Na+ diffusion constant)1.33⋅10−9 m2/s[31, 81]
Dk (K+ diffusion constant)1.96⋅10−9 m2/s[31, 81]
DCl (Cl diffusion constant)2.03⋅10−9 m2/s[31, 81]
DCa (Ca2+ diffusion constant)0.71⋅10−9 m2/s[31, 81]
λi (intracellular tortuosity)3.2[60, 82]
λe (extracellular tortuosity)1.6[60, 82]
γNa, γK, γCl (intracellular fractions of free ions)1
γCa (intracellular fraction of free ions)0.01

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 Ca2+) in four compartments, we have 20 unknown parameters (16 for [k] and four for ϕ ), while Eqs 36 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, Cm is the total capacitance of the membrane, i.e., Cm = cmAm, where cm 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 thus Qdi=Fkzk[k]diVdi. By inserting this into Eq 10, we obtain the final expression for ϕdi:

By inserting the equivalent expression for Qsi 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 (Qi = −Qe) 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 dQdi/dt = −dQde/dt and dQsi/dt = −dQse/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 ii and ie are the intra- and extracellular current densities, respectively. To find an expression for these, we multiply Eqs 1 and 2 by Fzk 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 idiff,e, ifield,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 idiff 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, Ek is the reversal potential, F is the Faraday constant, and zk 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 (jNa) and a K+ delayed rectifier flux (jK−DR), while the dendrite contained a voltage-dependent Ca2+ flux (jCa), a voltage-dependent K+ AHP flux (jK−AHP), and a Ca2+-dependent K+ flux (jK−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:

All these equations were taken from [8] (with errata [83]) and converted so that values are given in SI units: units for rates (α’s, β’s) are 1/s, unit for τz is s, and units for voltages ϕ are V. The equations were used in their original form, except those related to Ca2+ dynamics, where we made the following changes: Firstly, as a large fraction of intracellular Ca2+ is buffered or taken up by the ER, we multiplied [Ca2+ ] in Eqs 51 and 52 by a factor γCa, which refers to the fraction of free Ca2+ within the cell, and set this to be 0.01. As [Ca2+ ] in Eqs 51 and 52 were multiplied with 0.01, only the free Ca2+ could affect the Ca2+ activated ion channels. We further assumed that only the free Ca2+ could move between the intracellular compartments (Eq 1) and affect the Ca2+ reversal potential (Eq 24). Secondly, the original PR model had an abstract and unitless variable for the intracellular Ca2+ concentration, with a basal concentration of 0.2, while we defined a (biophysically realistic) baseline concentration of 0.01 mM, which corresponds to a concentration of free Ca2+ of 100 nM. In Eqs 51 and 52 we therefore subtracted 99.8⋅10−6(mol/m3) from the Ca2+ concentration to correct for the shift in baseline. Thirdly, we modified the voltage-dependent Ca2+ 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 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 ρ, Ukcc2, and Unkcc1 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 Ca2+ decay was modeled in a similar fashion as in [3], but to ensure ion conservation we modeled it as an electroneutral Ca2+/2Na+ exchanger, exchanging one Ca2+ (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 Ca2+/3Na+ exchanger, a Ca2+ leakage current, SERCA pumps, etc. The decay flux density was defined as:

where UCa−dec is the decay rate, and [Ca2+]i,b is the basal Ca2+ 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:

Fig 2 summarizes the model. The parameters involved in this model and their values used in this study are listed in Tables 14.

Table 3
Temperature and physical constants.
T (absolute temperature)309.14 K[45]*
F (Faraday constant)9.648⋅104 C/mol
R (gas constant)8.314 J/(mol.K)
* The temperature is not explicitly given in [45], but from Eq 3 in [45] we know that RTF=26.64·10-3V. By using the values of R and F listed in Table 3, we get an absolute temperature of 309.14 K, corresponding to a body temperature of 36°C.
Table 4
Membrane parameters.
cm3⋅10−2 F/m2[3, 8]
g¯Na,leak0.247 S/m2[45]
g¯K,leak0.5 S/m2[45]
g¯Cl,leak1.0 S/m2[45]
g¯Na300 S/m2[3, 8]
g¯DR150 S/m2[3, 8]
g¯Ca118 S/m2
g¯AHP8 S/m2[3, 8]
g¯C150 S/m2[3, 8]
ρ1.87⋅10−6 mol/(m2s)[45]*
Ukcc27.0⋅10−7 mol/(m2s)[45]*
Unkcc12.33⋅10−7 mol/(m2s)[45]*
UCa−dec75 s−1[3, 8]
* We multiplied the original values from [45] by a conversion factor 73·10-6 m to convert the units from mM/s to mol/m2 s. The conversion factor equals the initial inverse surface area to volume ratio from [45].

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 = −68mV and ϕdm,0 = −68mV. The other initial conditions were n0 = 0.001, h0 = 0.999, s0 = 0.009, c0 = 0.007, q0 = 0.01, and [Ca2+]0 = 0.2, same as in [3].



The parameters listed in Tables 14 were used in all the simulations of the electrodiffusive Pinsky-Rinzel (edPR) model. We tuned the Ca2+ conductance g¯Ca manually to obtain comparable spike shapes between the edPR model and the original PR model, as well as the fraction of free Ca2+ inside the cell, and the coupling strength between the soma and the dendrite.

In the edPR model, the coupling strength between the soma and dendrite was proportional to the ratio Ai/Δx, and all model outputs depended on this ratio, and not on Ai or Δx in isolation. By choice, we adjusted the coupling strength by varying Ai = αAm through adjusting the parameter α. We could have obtained the equivalent effect by varying Δx instead. The default value of α was set to 2. All simulations were run using this value, except in Fig 5C where α was set to 0.43.

In the original PR model, the coupling strength between the soma and dendrite was represented by a coupling conductance gc, which had a default value of 10.5mS/cm2 . In Fig 5A, gc was set to 2.26mS/cm2.

Initial conditions

The initial conditions for the edPR model were obtained through a two-step procedure. In the first step, we specified a set of pre-calibration initial values: We set (i) the initial membrane potential, ϕm,0 , to -68 mV, (ii) the concentrations to the pre-calibrated values in Table 5, and (iii) the gating variables (Table 5) to the same initial values as in [3]. Based on the initial concentration values, we also defined (iv) a set of static intracellular and extracellular residual charges, representing negatively charged macromolecules present in real neurons. We represented these as constant concentrations ([X]i,0 and [X]e,0) of anions with charge number zX = −1 (assuming this to be the mean charge number of the real macromolecules) and diffusion constants DX = 0 (assuming immobility). The residual charges were introduced to ensure consistency between the initial membrane potential and the charge density associated with the initial ion concentrations:


Table 5
Initial conditions.
ϕm,0 -68 mV-67.7 mV
[Na+]i,015 mM16.9 mM
[Na+]e,0145 mM141.2 mM
[K+]i,0140 mM139.5 mM
[K+]e,05 mM5.9 mM
[Cl]i,04 mM5.4 mM
[Cl]e,0110 mM107.1 mM
[Ca2+]i,00.01 mM*0.01 mM*
[Ca2+]e,01.1 mM1.1 mM
[X]i,0 151.0 mM151.0 mM
[X]e,0 42.2 mM42.2 mM
1 Preciser values (with more decimals included) were read to/from file and used in the simulations. (Available at
ϕm is not an independent state variable, but at each time point an algebraic function of ion concentrations, as computed through the KNP formalism.
* Only 1% of the total intracellular Ca2+, that is, a 100 nM, was assumed to be free (unbuffered).
Not state variables, but constants, derived (initially) from Eqs 75 and 76.

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: ENa = 57mV, EK = −84mV, ECl = −79mV, and ECa = 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 istim 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:



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 AiNAjk,diff,i from time zero to t, where NA is the Avogadro constant. Similarly, we integrated AeNAjk,diff,e to calculate the accumulative transport of ions in the extracellular solution due to diffusion. We did the same calculations with jk,drift to study the accumulative transport of ions due to drift. When knowing the accumulative transport of each ion species, kakkum, we calculated the total transport of e+ from their weighted sum:


Fig 10: To calculate ϕVC,e and ϕdiff,e, we looked at the extracellular axial current as it is given in the KNP formalism:

where the last equality follows when we insert Eq 18 for the extracellular field-driven current density ifield,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

Since we already knew ie and idiff,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 and



AL Hodgkin, AF Huxley. . A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol. 1952;117:, pp.500–544. , doi: 10.1113/jphysiol.1952.sp004764


RD Traub, RK Wong, R Miles, H Michelson. . A model of a CA3 hippocampal pyramidal neuron incorporating voltage-clamp data on intrinsic conductances. Journal of Neurophysiology. 1991;66(2):, pp.635–650. , doi: 10.1152/jn.1991.66.2.635


PF Pinsky, J Rinzel. . Intrinsic and network rhythmogenesis in a reduced Traub model for CA3 neurons. Journal of computational neuroscience. 1994;1(1-2):, pp.39–60. , doi: 10.1007/bf00962717


ZF Mainen, J Joerges, JR Huguenard, TJ Sejnowski. . A model of spike initiation in neocortical pyramidal neurons. Neuron. 1995;15(6):, pp.1427–1439. , doi: 10.1016/0896-6273(95)90020-9


M Pospischil, M Toledo-Rodriguez, C Monier, Z Piwkowska, T Bal, Y Fr?gnac, et al. Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern. 2008;99(4-5):, pp.427–441. , doi: 10.1007/s00422-008-0263-8


G Halnes, S Augustinaite, P Heggelund, GT Einevoll, M Migliore. . A multi-compartment model for interneurons in the dorsal lateral geniculate nucleus. PLoS Comput Biol. 2011;7(9):, pp.e1002160, doi: 10.1371/journal.pcbi.1002160


E Hay, S Hill, F Schürmann, H Markram, I Segev. . Models of neocortical layer 5b pyramidal cells capturing a wide range of dendritic and perisomatic active properties. PLoS computational biology. 2011;7(7):, pp.e1002107, doi: 10.1371/journal.pcbi.1002107


D Sterratt, B Graham, A Gillies, D Willshaw. Principles of computational modelling in neuroscience. Cambridge University Press; 2011.


FF Offner. . Ion flow through membranes and the resting potential of cells. The Journal of membrane biology. 1991;123(2):, pp.171–182. , doi: 10.1007/bf01998087


C Koch. Biophysics of computation: information processing in single neurons. Oxford university press; 2004.


W Rall. . Core conductor theory and cable properties of neurons. Comprehensive physiology. 2011; p., pp.39–97.


A Tveito, KH Jæger, GT Lines, Ł Paszkowski, J Sundnes, AG Edwards, et al. An Evaluation of the Accuracy of Classical Models for Computing the Membrane Potential and Extracellular Potential for Neurons. Frontiers in Computational Neuroscience. 2017;11(April):, pp.1–18.


RD Traub, D Contreras, MO Cunningham, H Murray, FE LeBeau, A Roopun, et al. Single-column thalamocortical network model exhibiting gamma oscillations, sleep spindles, and epileptogenic bursts. Journal of neurophysiology. 2005;93(4):, pp.2194–2232. , doi: 10.1152/jn.00983.2004


H Markram, E Muller, S Ramaswamy, MW Reimann, M Abdellah, CA Sanchez, et al. Reconstruction and simulation of neocortical microcircuitry. Cell. 2015;163(2):, pp.456–492. , doi: 10.1016/j.cell.2015.09.029


A Arkhipov, NW Gouwens, YN Billeh, S Gratiy, R Iyer, Z Wei, et al. Visual physiology of the layer 4 cortical circuit in silico. PLoS computational biology. 2018;14(11):, pp.e1006535, doi: 10.1371/journal.pcbi.1006535


G Buzsáki, CA Anastassiou, C Koch. . The origin of extracellular fields and currents?EEG, ECoG, LFP and spikes. Nature reviews neuroscience. 2012;13(6):, pp.407, doi: 10.1038/nrn3241


GT Einevoll, C Kayser, NK Logothetis, S Panzeri. . Modelling and analysis of local field potentials for studying the function of cortical circuits. Nature Reviews Neuroscience. 2013;14(11):, pp.770–785. , doi: 10.1038/nrn3599


E Hagen, S Næss, TV Ness, GT Einevoll. . Multimodal modeling of neural network activity: computing LFP, ECoG, EEG and MEG signals with LFPy 2.0. Frontiers in neuroinformatics. 2018;12:, pp.92, doi: 10.3389/fninf.2018.00092


CA Anastassiou, C Koch. . Ephaptic coupling to endogenous electric field activity: why bother?Current opinion in neurobiology. 2015;31:, pp.95–103. , doi: 10.1016/j.conb.2014.09.002


GR Holt, C Koch. . Electrical interactions via the extracellular potential near cell bodies. Journal of computational neuroscience. 1999;6(2):, pp.169–184. , doi: 10.1023/a:1008832702585


KH Pettersen, GT Einevoll. . Amplitude variability and extracellular low-pass filtering of neuronal spikes. Biophysical journal. 2008;94(3):, pp.784–802. , doi: 10.1529/biophysj.107.111179


GG Somjen. . Mechanisms of Spreading Depression and Hypoxic Spreading Depression-Like Depolarization. Physiol Rev. 2001;81(3):, pp.1065–1096. , doi: 10.1152/physrev.2001.81.3.1065


F Fröhlich, M Bazhenov. . Potassium dynamics in the epileptic cortex: new insights on an old topic. Neuroscientist. 2008;14(5):, pp.422–433. , doi: 10.1177/1073858408317955


BJ Zandt, B ten Haken, MJ van Putten, MA Dahlem. . How does spreading depression spread? Physiology and modeling. Reviews in the Neurosciences. 2015;26(2):, pp.183–198. , doi: 10.1515/revneuro-2014-0069


C Ayata, M Lauritzen. . Spreading depression, spreading depolarizations, and the cerebral vasculature. Physiological Reviews. 2015;95(3):, pp.953–993. , doi: 10.1152/physrev.00027.2014


O Herreras, G Somjen. . Analysis of potential shifts associated with recurrent spreading depression and prolonged unstable spreading depression induced by microdialysis of elevated K+ in. Brain research. 1993;610:, pp.283–294. , doi: 10.1016/0006-8993(93)91412-l


H Kager, WJ Wadman, GG Somjen. . Simulated seizures and spreading depression in a neuron model incorporating interstitial space and ion concentrations. Journal of neurophysiology. 2000;84(1):, pp.495–512. , doi: 10.1152/jn.2000.84.1.495


E Barreto, JR Cressman. . Ion concentration dynamics as a mechanism for neuronal bursting. Journal of biological physics. 2011;37(3):, pp.361–373. , doi: 10.1007/s10867-010-9212-6


L Øyehaug, I Østby, CM Lloyd, SW Omholt, GT Einevoll. . Dependence of spontaneous neuronal firing and depolarisation block on astroglial membrane transport mechanisms. Journal of computational neuroscience. 2012;32(1):, pp.147–165. , doi: 10.1007/s10827-011-0345-9


Y Mori. . A multidomain model for ionic electrodiffusion and osmosis with an application to cortical spreading depression. Physica D: Nonlinear Phenomena. 2015;308:, pp.94–108. , doi: 10.1016/j.physd.2015.06.008


G Halnes, T Mäki-Marttunen, D Keller, KH Pettersen, OA Andreassen, GT Einevoll. . Effect of ionic diffusion on extracellular potentials in neural tissue. PLoS computational biology. 2016;12(11):, pp.e1005193, doi: 10.1371/journal.pcbi.1005193


A Solbrå, AW Bergersen, J Van Den Brink, A Malthe-Sørenssen, GT Einevoll, G Halnes. . A Kirchhoff-Nernst-Planck framework for modeling large scale extracellular electrodiffusion surrounding morphologically detailed neurons. PLoS computational biology. 2018;14(10):, pp.e1006510, doi: 10.1371/journal.pcbi.1006510


J Kneller, RJ Ramirez, D Chartier, M Courtemanche, S Nattel. . Time-dependent transients in an ionically based mathematical model of the canine atrial action potential. American Journal of Physiology-Heart and Circulatory Physiology. 2002;282(4):, pp.H1437–H1451. , doi: 10.1152/ajpheart.00489.2001


G Somjen, H Kager, W Wadman. . Computer simulations of neuron-glia interactions mediated by ion flux. Journal of computational neuroscience. 2008;25(2):, pp.349–365. , doi: 10.1007/s10827-008-0083-9


G Florence, Ma Dahlem, ACG Almeida, JWM Bassani, J Kurths. . The role of extracellular potassium dynamics in the different stages of ictal bursting and spreading depression: a computational study. Journal of theoretical biology. 2009;258(2):, pp.219–28. , doi: 10.1016/j.jtbi.2009.01.032


JR Cressman, G Ullah, J Ziburkus, SJ Schiff, E Barreto. . The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: I. Single neuron dynamics. Journal of computational neuroscience. 2009;26(2):, pp.159–170. , doi: 10.1007/s10827-008-0132-4


G Ullah, JR Cressman Jr, E Barreto, SJ Schiff. . The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: II. Network and glial dynamics. Journal of computational neuroscience. 2009;26(2):, pp.171–183. , doi: 10.1007/s10827-008-0130-6


J Lee, SJ Kim. . Spectrum measurement of fast optical signal of neural activity in brain tissue and its theoretical origin. Neuroimage. 2010;51(2):, pp.713–722. , doi: 10.1016/j.neuroimage.2010.02.076


J Lee, DA Boas, SJ Kim. . Multiphysics neuron model for cellular volume dynamics. IEEE Transactions on Biomedical Engineering. 2011;58(10):, pp.3000–3003. , doi: 10.1109/TBME.2011.2159217


BJ Zandt, B ten Haken, JG van Dijk, MJ van Putten. . Neural dynamics during anoxia and the “wave of death”. PLoS One. 2011;6(7):, pp.e22127, doi: 10.1371/journal.pone.0022127


N Hübel, E Schöll, MA Dahlem. . Bistable dynamics underlying excitability of ion homeostasis in neuron models. PLoS computational biology. 2014;10(5):, pp.e1003551, doi: 10.1371/journal.pcbi.1003551


MA Dahlem, J Schumacher, N Hübel. . Linking a genetic defect in migraine to spreading depression in a computational model. PeerJ. 2014;2:, pp.e379, doi: 10.7717/peerj.379


N Hübel, MA Dahlem. . Dynamics from seconds to hours in Hodgkin-Huxley model with time-dependent ion concentrations and buffer reservoirs. PLoS computational biology. 2014;10(12):, pp.e1003941, doi: 10.1371/journal.pcbi.1003941


Y Wei, G Ullah, J Ingram, SJ Schiff. . Oxygen and seizure dynamics: II. Computational modeling. Journal of neurophysiology. 2014;112(2):, pp.213–223. , doi: 10.1152/jn.00541.2013


Y Wei, G Ullah, SJ Schiff. . Unification of Neuronal Spikes, Seizures, and Spreading Depression. Journal of Neuroscience. 2014;34(35):, pp.11733–11743. , doi: 10.1523/JNEUROSCI.0516-14.2014


N Hübel, G Ullah. . Anions govern cell volume: a case study of relative astrocytic and neuronal swelling in spreading depolarization. PloS one. 2016;11(3):, pp.e0147060, doi: 10.1371/journal.pone.0147060


N Hübel, MS Hosseini-Zare, J Žiburkus, G Ullah. . The role of glutamate in neuronal ion homeostasis: A case study of spreading depolarization. PLoS computational biology. 2017;13(10):, pp.e1005804, doi: 10.1371/journal.pcbi.1005804


H Kager, W Wadman, G Somjen. . Conditions for the triggering of spreading depression studied with computer simulations. Journal of neurophysiology. 2002;88(5):, pp.2700–2712. , doi: 10.1152/jn.00237.2002


E Cataldo, M Brunelli, JH Byrne, E Av-Ron, Y Cai, DA Baxter. . Computational model of touch sensory cells (T Cells) of the leech: role of the afterhyperpolarization (AHP) in activity-dependent conduction failure. Journal of computational neuroscience. 2005;18(1):, pp.5–24. , doi: 10.1007/s10827-005-5477-3


H Kager, W Wadman, G Somjen. . Seizure-like afterdischarges simulated in a model neuron. Journal of computational neuroscience. 2007;22(2):, pp.105–128. , doi: 10.1007/s10827-006-0001-y


MD Forrest, MJ Wall, DA Press, J Feng. . The sodium-potassium pump controls the intrinsic firing of the cerebellar Purkinje neuron. PloS one. 2012;7(12):, pp.e51169, doi: 10.1371/journal.pone.0051169


JC Chang, KC Brennan, D He, H Huang, RM Miura, PL Wilson, et al. A mathematical model of the metabolic and perfusion effects on cortical spreading depression. PLoS One. 2013;8(8):, pp.e70469, doi: 10.1371/journal.pone.0070469


G Le Masson, S Przedborski, L Abbott. . A computational model of motor neuron degeneration. Neuron. 2014;83(4):, pp.975–988. , doi: 10.1016/j.neuron.2014.07.001


MD Forrest. . Simulation of alcohol action upon a detailed Purkinje neuron model and a simpler surrogate model that runs> 400 times faster. BMC neuroscience. 2015;16(1):, pp.27, doi: 10.1186/s12868-015-0162-6


GP Krishnan, G Filatov, A Shilnikov, M Bazhenov. . Electrogenic properties of the Na+/K+ ATPase control transitions between normal and pathological brain states. Journal of neurophysiology. 2015;113(9):, pp.3356–3374. , doi: 10.1152/jn.00460.2014


A Zylbertal, A Kahan, Y Ben-Shaul, Y Yarom, S Wagner. . Prolonged intracellular Na+ dynamics govern electrical activity in accessory olfactory bulb mitral cells. PLoS biology. 2015;13(12):, pp.e1002319, doi: 10.1371/journal.pbio.1002319


A Zylbertal, Y Yarom, S Wagner. . The Slow Dynamics of Intracellular Sodium Concentration Increase the Time Window of Neuronal Integration: A Simulation Study. Frontiers in computational neuroscience. 2017;11:, pp.85, doi: 10.3389/fncom.2017.00085


N Qian, TJ Sejnowski, S Diego. Biological Cybernetics. 1989;15:, pp.1–15. , doi: 10.1007/BF00217656


Y Mori, C Peskin. . A numerical method for cellular electrophysiology based on the electrodiffusion equations with internal boundary conditions at membranes. Communications in Applied Mathematics and Computational Science. 2009;4.1:, pp.85–134. , doi: 10.2140/camcos.2009.4.85


G Halnes, I Østby, KH Pettersen, SW Omholt, GT Einevoll. . Electrodiffusive model for astrocytic and neuronal ion concentration dynamics. PLoS computational biology. 2013;9(12):, pp.e1003386, doi: 10.1371/journal.pcbi.1003386


S Niederer. . Regulation of ion gradients across myocardial ischemic border zones: a biophysical modelling analysis. PloS one. 2013;8(4):, pp.e60323, doi: 10.1371/journal.pone.0060323


A Ellingsrud, A Solbrå, G Einevoll, G Halnes, M Rognes. . Finite element simulation of ionic electrodiffusion in cellular geometries. Frontiers in Neuroinformatics. 2020; 14:11, doi: 10.3389/fninf.2020.00011


D Attwell, SB Laughlin. . An energy budget for signaling in the grey matter of the brain. Journal of Cerebral Blood Flow & Metabolism. 2001;21(10):, pp.1133–1145. , doi: 10.1097/00004647-200110000-00001


G Yi, Y Fan, J Wang. . Metabolic cost of dendritic Ca2+ action potentials in layer 5 pyramidal neurons. Frontiers in neuroscience. 2019;13, doi: 10.3389/fnins.2019.01221


In: MD Binder, N Hirokawa, U Windhorst, editors. Depolarization Block. Berlin, Heidelberg: Springer Berlin Heidelberg; 2009 p. , pp.943–944. Available from: , doi: 10.1007/978-3-540-29678-2_1453.


K Mizuseki, S Royer, K Diba, G Buzsáki. . Activity dynamics and behavioral correlates of CA3 and CA1 hippocampal pyramidal neurons. Hippocampus. 2012;22(8):, pp.1659–1680. , doi: 10.1002/hipo.22002


A Gardner-Medwin. . Analysis of potassium dynamics in mammalian brain tissue. The Journal of physiology. 1983;335:, pp.393–426. , doi: 10.1113/jphysiol.1983.sp014541


C Reiffurth, M Alam, M Zahedi-Khorasani, S Major, JP Dreier. . Na+/K+-ATPase α isoform deficiency results in distinct spreading depolarization phenotypes. Journal of Cerebral Blood Flow & Metabolism. 2019; p. 0271678X19833757.


P Nelson. Biological physics. WH FreemanNew York; 2004.


CM Van Rijn, H Krijnen, S Menting-Hermeling, AM Coenen. . Decapitation in rats: latency to unconsciousness and the?wave of death?PloS one. 2011;6(1):, pp.e16514.


RA McDougal, ML Hines, WW Lytton. . Reaction-diffusion in the NEURON simulator. Frontiers in neuroinformatics. 2013;7:, pp.28, doi: 10.3389/fninf.2013.00028


M Hines, AP Davison, E Muller. . NEURON and Python. Frontiers in neuroinformatics. 2009;3:, pp.1, doi: 10.3389/neuro.11.001.2009


G Halnes, T Mäki-Marttunen, KH Pettersen, OA Andreassen, GT Einevoll. . Ion diffusion may introduce spurious current sources in current-source density (CSD) analysis. Journal of Neurophysiology. 2017;118(1):, pp.114–120. , doi: 10.1152/jn.00976.2016


SL Gratiy, G Halnes, D Denman, MJ Hawrylycz, C Koch, GT Einevoll, et al. From Maxwell’s equations to the theory of current-source density analysis. European Journal of Neuroscience. 2017;45(8):, pp.1013–1023. , doi: 10.1111/ejn.13534


W Rall. Core conductor theory and cable properties of neurons In: ER Kandel, JM Brookhardt, V M Mountcastle, editors. Handbook of Physiology. Bethesda: American Physiological Society; 1977 p. , pp.39–97. Available from:


R O’Connell, Y Mori. . Effects of Glia in a Triphasic Continuum Model of Cortical Spreading Depression. Bulletin of Mathematical Biology. 2016;78(10):, pp.1943–1967. , doi: 10.1007/s11538-016-0206-9


A Tuttle, J Riera Diaz, Y Mori. . A computational study on the role of glutamate and NMDA receptors on cortical spreading depression using a multidomain electrodiffusion model. PLoS computational biology. 2019;15(12). , doi: 10.1371/journal.pcbi.1007455


RS Eisenberg, EA Johnson. . Three-dimensional electrical field problems in physiology. Progress in biophysics and molecular biology. 1970;20:, pp.1–65. , doi: 10.1016/0079-6107(70)90013-1


CS Henriquez. . Simulating the electrical behavior of cardiac tissue using the bidomain model. Critical reviews in biomedical engineering. 1993;21(1):, pp.1–77.


J Sundnes, BF Nielsen, KA Mardal, X Cai, GT Lines, A Tveito. . On the computational complexity of the bidomain and the monodomain models of electrophysiology. Annals of biomedical engineering. 2006;34(7):, pp.1088–1097. , doi: 10.1007/s10439-006-9082-z


SE Lyshevski. Nano and molecular electronics handbook. CRC Press; 2016.


KC Chen, C Nicholson. . Spatial buffering of potassium ions in brain extracellular space. Biophysical journal. 2000;78(6):, pp.2776–2797. , doi: 10.1016/S0006-3495(00)76822-6


Errata, Principles of computational modelling in neuroscience;.


JA Wolf, JT Moyer, MT Lazarewicz, D Contreras, M Benoit-Marand, P O’Donnell, et al. NMDA/AMPA ratio impacts state transitions and entrainment to oscillations in a computational model of the nucleus accumbens medium spiny projection neuron. Journal of Neuroscience. 2005;25(40):, pp.9080–9095. , doi: 10.1523/JNEUROSCI.2220-05.2005


D Purves, G Augustine, D Fitzpatrick, W Hall, A LaMantia, L White. Neuroscience. Sinauer Associates, Inc.: Sunderland, MA; 2012.


R Tyzio, A Ivanov, C Bernard, GL Holmes, Y Ben-Ari, R Khazipov. . Membrane potential of CA3 hippocampal pyramidal cells during postnatal development. Journal of neurophysiology. 2003;90(5):, pp.2964–2972. , doi: 10.1152/jn.00172.2003

25 Feb 2020

Dear Dr. Halnes,

Thank you very much for submitting your manuscript "An electrodiffusive, ion conserving Pinsky-Rinzel model with homeostatic mechanisms" for consideration at PLOS Computational Biology. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. The reviewers appreciated the attention to an important topic. Based on the reviews, we are likely to accept this manuscript for publication, providing that you modify the manuscript according to the review recommendations.

Please prepare and submit your revised manuscript within 30 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. 

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to all review comments, and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Thank you again for your submission to our journal. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.


William W Lytton, MD

Guest Editor

PLOS Computational Biology

Kim Blackwell

Deputy Editor

PLOS Computational Biology


A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact immediately:


Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The paper applies the Kirchoff-Nernst-Plank framework (a method for modeling electrodiffusion) to a reduced CA3 cell model. The original two compartment model is augmented with two additional extracellular compartments, modeling of both the intracellular and extracellular potentials, electrodiffusion of K+, Na+, Cl- and Ca2+ and several homoeostatic mechanisms. The model qualitatively matches the original, but unlike the simple model they started from the electrodiffusive model can be driven into depolarization block by excessive stimulation or depolarize when the homoeostatic mechanisms are impaired. This novel and innovative approach will be useful for the studying the mechanisms underling pathological conditions like ischemia and epilepsy, that can lead to spreading depolarization.

The sodium calcium exchanges takes in 3 sodium for every calcium ion at the cost of 1 ATP. Why do you have it exchanging 2 sodium ions instead? Eq 67 & 71.

The tortuosity is mentioned in the methods section, the value results from an increase in path length as ions diffusion around obstacles. Does this make sense for your model? If so, why haven't you also included the free volume fraction?

This seems like it could be an issue for interpreting the models boundary conditions as equivalent to period boundary conditions for a synchronized population, as the population would have a relatively large extracellular space.

Figure 4: Are the ion concentrations as stable as the somatic membrane potential?

Do all the simulations remain stable over the whole +/-15% parameter change or are some slowly depolarizing, like a less extreme 8A? (or hyperpolarizing).

Figure 6 and 7: Plotting changes freely diffusion calcium on the same scale as sodium is not useful. Why not plot the states as you did for figure 8?

Figure 8:

Why do you turn off all of the homoeostatic mechanisms, rather than just those two that depend on ATP?

G] Why does the intracellular dendritic calcium differ so much from the intracellular somatic calcium?

Initial conditions: How did you adjust the ion concentrations to arrive at stable initial condition? Would the model be unstable if you used initial with concentration from the literature? e.g. Ions in the Brain, Somjen. GG, 2002.

Also, why didn't you use the post calibration states as your initial conditions? What do the first 15s of the simulations look like?

Does this mean the initial (post-calibration) conditions for the sensitivity analysis were different for each simulations?

Minor issues/typos;

Line 392: ICRP model -- only mentioned here.

Figure 10: figure legend is not the same as the caption.

Table 1 caption: Two As, both intracellular volumes don't correspond to a sphere with radius 7um.

Reviewer #2: This paper proposed an electrodiffusive Pinsky-Rinzel (edPR) model, which is “a minimal neuronal model” with two neuronal compartments (a soma and a dendrite), plus two extracellular compartments (outside soma and outside dendrite), by adding homeostatic mechanisms and ion concentration dynamics. The edPR model doesnot only reproduce the membrane potential dynamics of the PR model, but also accounts for changes in neuronal firing properties due to deviations from baseline ion concentrations. They expect the model to be important because it opens for more detailed mechanistic studies of the pathological conditions associated with large changes in ion concentrations such as epilepsy and spreading depression.

Overall, the article is well organized and presented. I have some minor comments here.

1. The statement of “the first multicompartmental neuron model that in a biophysically consistent way does account for the effects of ion concentration variations ” is too strong. For example, Kager et al 2000 is a morphologically explicit multicompartmental model with ion concentration variations.

2. What is the advantage of using the two-compartmental model compared with one-compartment model or morphologically explicit model?

3. The model has two neuronal compartment compartments (a soma and a dendrite), it is a little confused how did the authors investigate the ionic diffusion along the axons in the statement of “investigate the consequences of neglecting the effect of ionic diffusion (along dendrites and axons) on the electrical potential”.

4. Table 3, T is body temperature at 36C (309.15K).


Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biologydata availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: None


PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here:


To enhance the reproducibility of your results, PLOS recommends that you deposit laboratory protocols in, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see

27 Mar 2020

Submitted filename: Response_letter.pdf

7 Apr 2020

Dear Dr. Halnes,

We are pleased to inform you that your manuscript 'An electrodiffusive, ion conserving Pinsky-Rinzel model with homeostatic mechanisms' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology. 

Best regards,

William W Lytton, MD

Guest Editor

PLOS Computational Biology

Kim Blackwell

Deputy Editor

PLOS Computational Biology


Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Thank you for addressing and answering all my issues and questions.


Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biologydata availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes


PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

20 Apr 2020


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

Dear Dr Halnes,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Sarah Hammond

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom | Phone +44 (0) 1223-442824 | | @PLOSCompBiol electrodiffusive, ion conserving Pinsky-Rinzel model with homeostatic mechanisms&author=Marte J. Sætra,Gaute T. Einevoll,Geir Halnes,William W Lytton,William W Lytton,Kim T. Blackwell,William W Lytton,Kim T. Blackwell,William W Lytton,Kim T. Blackwell,&keyword=&subject=Research Article,Biology and Life Sciences,Cell Biology,Cellular Types,Animal Cells,Neurons,Neuronal Dendrites,Biology and Life Sciences,Neuroscience,Cellular Neuroscience,Neurons,Neuronal Dendrites,Biology and Life Sciences,Cell Biology,Cellular Types,Animal Cells,Neurons,Biology and Life Sciences,Neuroscience,Cellular Neuroscience,Neurons,Biology and Life Sciences,Physiology,Electrophysiology,Membrane Potential,Medicine and Health Sciences,Physiology,Electrophysiology,Membrane Potential,Research and Analysis Methods,Simulation and Modeling,Biology and Life Sciences,Physiology,Physiological Processes,Homeostasis,Homeostatic Mechanisms,Medicine and Health Sciences,Physiology,Physiological Processes,Homeostasis,Homeostatic Mechanisms,Biology and Life Sciences,Cell Biology,Cellular Structures and Organelles,Cell Membranes,Intracellular Membranes,Biology and Life Sciences,Physiology,Electrophysiology,Membrane Potential,Action Potentials,Medicine and Health Sciences,Physiology,Electrophysiology,Membrane Potential,Action Potentials,Biology and Life Sciences,Physiology,Electrophysiology,Neurophysiology,Action Potentials,Medicine and Health Sciences,Physiology,Electrophysiology,Neurophysiology,Action Potentials,Biology and Life Sciences,Neuroscience,Neurophysiology,Action Potentials,Biology and Life Sciences,Cell Biology,Extracellular Space,