PLoS Computational Biology

Public Library of Science

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

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

Volume:
16,
Issue: 4

DOI 10.1371/journal.pcbi.1007661

Abstract

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.

The neuronal action potential (AP) is generated by a transmembrane influx of Na^{+}, which depolarizes the neuron, followed by an efflux of K^{+}, which repolarizes it. Likewise, all neurodynamics is fundamentally about the movement of ions, which are the charge carriers in the brain. Therefore, it might seem peculiar that most models of neuronal activity are based on the approximation that the concentrations of the main charge carriers (Na^{+}, K^{+}, and Cl^{−} ) do not change over time. This approximation is, for example, incorporated in the celebrated Hodgkin-Huxley model [1], and a large number of later models based on a Hodgkin-Huxley type formalism (see, e.g., [2–7]).

Setting the ion concentrations to not change over time is often a fairly good approximation. The reason is that the number of ions that need to cross the membrane to charge up the neuron by, say, an AP worth of millivolts, is too small to have any notable impact on ion concentrations on either side of the membrane (see, e.g., Box 2.2 in [8]), meaning that concentration changes on a short time scale can be neglected. On a longer time-scale, the ionic exchange due to APs (or other neuronal events), is normally reversed by a set of homeostatic mechanisms such as ion pumps and cotransporters, which work to maintain constant baseline concentrations. In Hodgkin-Huxley type models, the large number of ion pumps, cotransporters and passive ionic leakages that strive towards maintaining baseline conditions are therefore not explicitly modeled. Instead, they are simply assumed to do their job and are grouped into a single *passive* and non-specific leakage current *I*_{leak} = *g*_{leak}(*ϕ*_{m}−*E*_{leak} ), which determines the cell’s resting potential (for a critical study of this approximation, see [9]).

Another approximation commonly applied by modelers of neurons is that the extracellular potential is constant and grounded (*ϕ*_{e} = 0) so that the only voltage variable that one needs to worry about when simulating neurodynamics is the transmembrane potential (*ϕ*_{m} ). This assumption is implicit in the majority of morphologically explicit models of neurons, where the (spatial) signal propagation in dendrites and axons are computed using the cable equation (see, e.g., [10–12]). Cable-equation based, multicompartmental neuronal models are widely used within the field of neuroscience, both for understanding dendritic integration and neuronal response properties at the single neuron level (see, e.g., [3, 4, 6, 7]) and for exploring the dynamics of large neuronal networks (see e.g., [13–15]). They are even used in the context of performing forward modeling of extracellular potentials, such as local field potentials (LFP), the electrocorticogram (ECoG), and electroencephalogram (EEG) (see, e.g., [16–18]), despite the evident inconsistency involved when first computing neurodynamics under the approximation that *ϕ*_{e} = 0 (Fig 1A), and then in the next step using this dynamics to predict a nonzero *ϕ*_{e} (Fig 1B). The approximation is nevertheless useful since *ϕ*_{e} is typically so much smaller than *ϕ*_{m} that the (ephaptic) effect of *ϕ*_{e} on neurodynamics can be neglected without severe loss in accuracy [19].

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 [22–25]. 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 [27–29] 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) [30–32]. Accurate modeling of such systems thus requires a unified, electrodiffusive framework that ensures a physically consistent relationship between ion concentrations, charge density, and electrical potentials.

Intra- or extracellular electrodiffusion is not an issue in single-compartment models, of which there are quite a few that incorporate ion concentration dynamics in a more or less consistent way [28, 29, 33–47]. Single compartment models are useful in many aspects. However, in order to represent morphological features of neurons, such as e.g., differential expression of ion channels in the soma versus dendrites, or account for transport processes in the space inside or outside neurons, one needs models with more than a single compartment. Among the several morphologically explicit models that have included homeostatic machinery and explicitly simulated ion concentration dynamics (see e.g., [27, 48–57]), neither have accounted for the electrodiffusive coupling between the movement of ions and electrical potentials (see Results section titled Loss in accuracy when neglecting electrodiffusive effects on concentration dynamics). Hence, to our knowledge, no morphologically explicit neuron model has so far been developed that ensures biophysically consistent dynamics in ion concentrations and electrical potentials during long-time activity, although useful mathematical framework for constructing such models are available [58–62].

The goal of this work is to propose what we may refer to as “a minimal neuronal model that has it all”. By “has it all”, we mean that it (1) has a spatial extension, (2) considers both extracellular- and intracellular dynamics, (3) keeps track of all ion concentrations (Na^{+}, K^{+}, Ca^{2+}, and Cl^{−}) in all compartments, (4) keeps track of all electrical potentials (*ϕ*_{m}, *ϕ*_{e}, and *ϕ*_{i}—the latter being the intracellular potential) in all compartments, (5) has differential expression of ion channels in soma versus dendrites, and can fire somatic APs and dendritic calcium spikes, (6) contains the homeostatic machinery that ensures that it maintains a realistic dynamics in *ϕ*_{m} and all ion concentrations during long-time activity, and (7) accounts for transmembrane, intracellular and extracellular ionic movements due to both diffusion and electrical migration, and thus ensures a consistent relationship between ion concentrations and electrical charge. Being based on a unified framework for intra- and extracellular dynamics (Fig 1C), the model thus accounts for possible ephaptic effects from extracellular dynamics, as neglected in standard feedforward models based on volume conductor theory (Fig 1A and 1B). By “minimal” we simply mean that we reduce the number of spatial compartments to the minimal, which in this case is four, i.e., two neuronal compartments (a soma and a dendrite), plus two extracellular compartments (outside soma and outside dendrite). Technically, the model was constructed by adding homeostatic mechanisms and ion concentration dynamics to an existing model, i.e., the two-compartment Pinsky-Rinzel (PR) model [3], and embedding in it a consistent electrodiffusive framework, i.e., the previously developed Kirchhoff-Nernst-Planck framework [31, 32, 60, 62]. For the remainder of this paper, we refer to our model as the electrodiffusive Pinsky-Rinzel (edPR) model.

The remainder of this article is organized as follows. First, we present the edPR model and illustrate the numerous variables that it can simulate. Next, we show that the edPR model can reproduce the firing properties of the original PR model. By running long-time simulations (several minutes of biological time) on both models, we identify the firing conditions under which the two models maintained a similar firing pattern, and under which conditions concentration effects became important so that dynamics of the edPR model diverged from the original PR model over time. Finally, we use the edPR model to explore the validity of some important assumptions commonly made in the field of computational neuroscience, regarding the decoupling of electrical and diffusive signals. We believe that the edPR model will be of great value for the field of neuroscience, partly because it gives a deepened insight into the balance between neuronal firing and ion homeostasis, partly because it lends itself to explore under which conditions the common modeling assumption of constant ion concentrations is warranted, and most importantly because it opens for more detailed mechanistic studies of pathological conditions associated with large changes in ion concentrations, such as epilepsy and spreading depression [22–25].

The here proposed electrodiffusive Pinsky-Rinzel (edPR) model is inspired by the original Pinsky-Rinzel (PR) model [3], which is a two-compartment (soma + dendrite) version of a CA3 hippocampal cell model, initially developed by Traub et al. [2]. In the original PR model, the somatic compartment contains Na^{+}, and K^{+} delayed rectifier currents (*I*_{Na} and *I*_{K−DR}), while the dendritic compartment contains a voltage-dependent Ca^{2+} current (*I*_{Ca}), a voltage-dependent K^{+} afterhyperpolarization current (*I*_{K−AHP}), and a Ca^{2+}-dependent K^{+} current (*I*_{K−C}). In addition, both compartments contain passive leakage currents. Despite its small number of compartments and conductances, the PR model can reproduce a variety of realistic firing patterns when responding to somatic or dendritic stimuli, including somatic APs and dendritic calcium spikes.

In the edPR model, we have adopted all mechanisms from the original PR model. In addition, we have (i) made all ion channels and passive leakage currents ion-specific, (ii) included 3Na^{+}/2K^{+} pumps (*I*_{pump}), K^{+}/Cl^{−} cotransporters (*I*_{KCC2}), Na^{+}/K^{+}/2Cl^{−} cotransporters (*I*_{NKCC1}), and Ca^{2+}/2Na^{+} exchangers (*I*_{Ca−dec} ), and (iii) included two extracellular compartments (outside soma + outside dendrite). To compute the dynamics of the edPR, we used an electrodiffusive KNP-framework for consistently computing the voltage- and ion concentration dynamics in the intra- and extracellular compartments [60]. The model is summarized in Fig 2 and described in details in the Methods section.

While the key variable in the original PR model is the membrane potential *ϕ*_{m} , the edPR model allows us to compute a multitude of variables relevant to neurodynamics. The functionality of the edPR model is illustrated in Fig 3, which shows a 60 s simulation where the model fires at 1 Hz for 10 s. We have plotted a selection of output variables, including the membrane potentials (Fig 3A and 3B), extracellular potentials (Fig 3C and 3D), the dynamics of all ion concentrations in all compartments (Fig 3E–3H), concentration effects on ionic reversal potentials (Fig 3I–3J), concentration effects on the electrical conductivity of the intra- and extracellular medium (Fig 3K), and ATP consumption (Fig 3L) of the 3Na^{+}/2K^{+} pumps and Ca^{2+}/2Na^{+} exchangers.

Unlike neuronal models based on cable theory, where *ϕ*_{e} is assumed to be zero so that *ϕ*_{m} = *ϕ*_{i}, the edPR model computes *ϕ*_{m}, *ϕ*_{i}, and *ϕ*_{e} from a consistent framework where ephaptic effects from *ϕ*_{e} on *ϕ*_{m} are accounted for (Fig 3C). Due to the electrical coupling between the soma and dendrite, the fluctuations in *ϕ*_{m} were similar in these compartments, and a more detailed analysis of the AP shapes is found further below. While an action potential essentially gave a depolarization followed by a repolarization of *ϕ*_{m} , its extracellular signature was essentially a voltage drop (to about -5 mV) followed by a voltage increase (to about +5 mV). This biphasic response of the extracellular AP signature has been seen in several studies (for an analysis, see [20, 21]). In experimental recordings, amplitudes in *ϕ*_{e} fluctuations are typically on the order of 100 *μ* V, which is much smaller than that predicted by the edPR model. The discrepancy is an artifact that is mainly due to the 1D approximation in the edPR model (see Discussion). The dendritic extracellular potential (Fig 3D) was by definition zero at all times, as this compartment was used as the reference point for the potential.

The effect of neuronal firing on the ion concentration dynamics is illustrated in Fig 3E–3H. Before the stimulus onset, the cell was resting at approximately -68 mV, and ion concentrations remained at baseline values. During AP firing, the ion concentrations varied in a jigsaw-like fashion in all compartments, except for Ca^{2+}, which returned to baseline between each AP and showed notable variation only inside/outside the dendrite since the soma contained no Ca^{2+} channels. As the extracellular volume was set to be half as big as the intracellular volume, changes in extracellular ion concentrations were about twice as big as the changes in intracellular ion concentrations. The jigsaw pattern was most pronounced for the K^{+} and Na^{+} concentrations, as these were the main mediators of the APs (Fig 3E–3H). The pattern reflects a cycle of (i) incremental steps away from baseline concentrations, which were mediated by the complex of mechanisms active during the APs, followed by (ii) slower decays back towards baseline, which were mediated by pumps and cotransporters working between the APs. In this simulation, the decay was incomplete, so that concentrations reached gradually larger peak values by each consecutive AP. However, as we show later (see Section titled The edPR model predicts homeostatic failure due to high firing rate), the concentrations did, in this case, approach a firing-frequency dependent steady state.

When the firing ceased in Fig 3, the pumps and cotransporters could work uninterruptedly to re-establish the baseline ion concentrations. The resting membrane potential of about -68 mV, was recovered quite rapidly (ms timescale). After this, the slower recovery process of the ion concentration was due to an electroneutral exchange of ions between the neuron and the extracellular space. A full recovery of the baseline concentrations took on the order of 80 s (confirmed by running a longer simulation than the one shown in Fig 3).

As ion concentrations varied during the simulation, so did the ionic reversal potentials, *E*_{k} (Fig 3I–3J). The by far largest change was seen for the Ca^{2+} reversal potential in the dendrite (*E*_{k,d}), which dropped by as much as -30 mV during an AP, (i.e., from a baseline value of 124 mV to 94 mV). The explanation is that the basal intracellular Ca^{2+}-concentration is extremely low (100 nM) compared to the concentrations of other ion species (several mM), and therefore experienced a much larger relative change during the simulation. Among the main charge carriers (Na^{+}, Cl^{−}, K^{+}), the lowest concentration is found for K^{+} in the extracellular space (Table 5 in Methods). For that reason, the second largest change in reversal potential was found for *E*_{K}, which increased by about 5 mV (i.e., from a basal value of -84 mV to -79 mV) in both the soma and dendrite. The changes in *E*_{Ca} and *E*_{K} had a relatively minor impact on the firing pattern in the shown simulations, as the relative change in the driving force *ϕ*_{m}−*E*_{k} was not that severe.

The conductivities (*σ* ) of the intra- and extracellular bulk solutions depend on the availability of free charge carriers, and are in the edPR model functions of the ion concentrations and their mobility (cf. Eq 19). The changes in *σ* were minimal during the conditions simulated here (Fig 3K), i.e., *σ* varied by a few *μ*S/m over the course of the simulation, while the basal levels were approximately 0.11 S/m and 0.59 S/m for the intra- and extracellular solutions, respectively.

Finally, the 3Na^{+}/2K^{+} pump and Ca^{2+}/2Na^{+} exchanger use energy in the form of ATP to move ions against their gradients. The 3Na^{+}/2K^{+} pump exchanges 3 Na^{+} ions for 2 K^{+} ions, and consumes one ATP per cycle [63], while we assumed that the Ca^{2+}/2Na^{+} exchanger consumed 1 ATP per cycle (i.e., per Ca^{2+} exchanged, as in [64]). As the edPR model explicitly models these processes, we could compute the ATP (energy) consumption of the pumps during the simulation. Fig 3L shows the accumulative number of ATP consumed from the onset of the simulation. The 3Na^{+}/2K^{+} pump was constantly active, as it combated leakage currents and worked to maintain the baseline concentration even before stimulus onset. Before stimulus onset, it consumed ATP at a constant rate (linear curve), which increased only slightly at *t* = 10 s when the neuron started to fire (small dent in the curve). As the neuron did not contain any passive leakage of Ca^{2+}, the Ca^{2+}/2Na^{+} exchangers were only active while the neuron was firing. During firing, the Ca^{2+}/2Na^{+} exchanger combated the Ca^{2+} entering through the dendritic Ca^{2+} channels, and then consumed approximately the same amount of energy as the 3Na^{+}/2K^{+} pump (parallel curves). A high metabolic cost of dendritic Ca^{2+} spikes has previously been reported also in cortical layer 5 pyramidal neurons [64].

We note that the edPR model had a stable resting state before stimulus onset and that it returned to this resting state after the stimulus had been turned off. In this resting state, ion concentrations remained constant, and *ϕ*_{m} was approximately -68 mV. This resting equilibrium was due to a balance between the ion-specific leakage channels, pumps, and cotransporters, which we adopted from previous studies (see Methods). However, the existence of a homeostatic equilibrium was not highly sensitive to the choice of model parameters. As we confirmed through a sensitivity analysis, varying membrane parameters (one by one) with ±15% of their default values did not abolish the existence of a stable resting state, but shifted the resting potential by maximally ±3 mV (Fig 4A) and the resting concentrations by maximally 5% (Fig 4B–4E).

A motivation behind basing the electrodiffusive (edPR) model on a previously developed (PR) model, was that we wanted to use the firing properties of the original PR model as a “ground truth” when constraining the edPR model. In particular, we wanted the edPR model to qualitatively reproduce the interplay between somatic action potentials and dendritic Ca^{2+} spikes, as this was an essential feature of the original PR model [3]. In the PR model, this interplay depended strongly on the coupling strength (coupling conductance) between the soma and dendrite compartment. A weak coupling resulted in a wobbly ping-pong effect, where a somatic AP triggered a dendritic Ca^{2+} spike, which in turn fed back to the soma, giving rise to secondary oscillations in *ϕ*_{m} (Fig 5A). With a strong (about five times stronger) coupling, the somatic and dendritic responses became more similar in shape, as expected (Fig 5B).

Since the edPR model contained membrane mechanisms and ephaptic effects not present in the PR model, selected parameters in the edPR model had to be re-tuned in order to obtain similar firing as the PR model (see Methods). With the selected parameterization of the edPR model (see the Parameterizations section), we were able to reproduce the characteristic features seen in the PR model for a weak (Fig 5C) and strong (about five times stronger) coupling between the soma and dendrite (Fig 5D).

As previously discussed, the PR model was, as most existing neuronal models, constructed under the assumption that ion concentration effects are negligible, an assumption that is justified for short term neurodynamics, and for long term dynamics provided that the activity level is sufficiently low for the homeostatic mechanisms to maintain concentrations close to baseline over time. Hence, we expect there to be a scenario (S1) with a moderately low firing rate, where the PR and edPR models can fire continuously and regularly over a long time exhibiting similar firing properties, and another scenario (S2) with a higher firing rate, where the PR and edPR models exhibit similar firing properties initially in the simulation, but where the dynamics of the two models diverge over time due to homeostatic failure accounted for by the edPR model, but not the PR model (which ad hoc assume perfect homeostasis). Simulations of two such scenarios are shown in Figs 6 and 7.

To simulate scenario S1, the PR model (Fig 6A and 6B) and edPR model (Fig 6C–6J) were given a constant input (see figure caption) that gave them a firing rate of 1 Hz. Both models settled at a regular firing rate, and in neither of them the firing pattern changed over time, even in simulations of as much as an hour of biological time. For the edPR model, the S1 scenario is the same as that simulated for only a brief period in Fig 3. There, we observed that the ion concentrations gradually changed during the first seconds after the onset of stimulus (Fig 3E–3H). However, for endured firing, the ion concentrations and reversal potentials settled on a (new) dynamic steady state (Fig 6E–6J), where they deviated by ∼1-5 mM from the baseline concentrations during rest (i.e., for edPR receiving no input). The apparent “thickness” of the curves (e.g., the red curve for K^{+} in Fig 6H) is due to concentration fluctuations at the short time scale of AP firing. However, after each AP, the homeostatic mechanisms managed to 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.

Above we simulated homeostatic failure occurring because the firing rate became too high for the homeostatic mechanisms to keep up (S2). Homeostatic failure may also occur due to impairment of the homeostatic mechanisms, either due to genetic mutations (see, e.g., [68]) or because the energy supply is reduced, such as after a stroke (see, e.g., [25]). Here, we have used the edPR model to simulate a version of this, i.e., a third scenario (S3) where the ATP-dependent mechanisms, that is, the 3Na^{+}/2K^{+} pumps and the Ca^{2+}/2Na^{+} exchangers, were turned off.

In S3, the neuron received no external input, so that the dynamics of the neuron was solely due to gradually dissipating transmembrane ion concentration gradients. After an initial transient, we observed a slow and gradual increase in the membrane potential for about 48 s (Fig 8A). This coincided with a slow and gradual change in the ion concentrations (Fig 8D–8G) and ionic reversal potentials (Fig 8B and 8C) due to predominantly passive leakage over the membrane.

At about *t* = 48 s, the membrane potential reached the firing threshold, at which point the active channels started to use what was left of the concentration gradients to generate action potentials and Ca^{2+} spikes. This resulted in a burst of activity. During this bursts of activity, the concentration gradients dissipated even faster, since both active and passive channels were then open. As a consequence, the “resting” membrane potential was further depolarized and the neuron went into depolarization block [65]. After this, the neuron continued to “leak” until it settled at a new steady state. The non-zero final equilibrium potential is known as the Donnan equilibrium or the Gibbs-Donnan equilibrium [69]. The reason why the cell did not approach an equilibrium with *ϕ*_{m} = 0 and identical ion concentrations on both side of the membrane, is that the model contained static residual charges, representing negatively charged macromolecules typically residing in the intracellular environment (see Methods), the sum of which resulted in a final state with a negatively charged inside. In addition, since the Ca^{2+} channel inactivated, and since the model had no passive Ca^{2+} leakage, Ca^{2+} could end up being trapped inside/outside the membrane and did not by necessity approach the Donnan equilibrium, although it was close to it.

As the Ca^{2+} dynamics in Fig 8G may seem counterintuitive, we here give some additional explanation of it. During the burst and initial stages of the depolarization block, the dendritic Ca^{2+} channels were open. Extracellular Ca^{2+} then diffused from the soma towards the dendrite, where it flowed into the neuron. This resulted in a low Ca^{2+} concentration in both extracellular compartments and a high Ca^{2+} concentration in the intracellular dendritic compartment. The reason why the intracellular Ca^{2+} equilibrated more slowly than the extracellular, was that, by assumption, only 1% of the intracellular Ca^{2+} concentration was unbuffered and free to diffuse (see Methods), hence, the effective intracellular concentration gradient was a factor 100 lower than it “appears” in Fig 8G.

A pattern resembling that in Fig 8A, i.e., a period of silence, followed by a burst of activity, and then silence again, has been seen in experimental EEG recordings of decapitated rats [70], where the activity burst was referred to as “the wave of death”, and the phenomenon was ascribed to the lack of energy supply to homeostatic mechanisms. The simulation in Fig 8A represents the single-cell correspondence to this death wave. We note that this phenomenon has been simulated and analyzed thoroughly in a previous modeling study, using a simpler, single compartmental model with ion conservation [40]. We, therefore, do not analyze it further here.

The concentration-dependent effects studied in the previous subsection were predominantly due to changes in ionic reversal potentials. Effects like this could, therefore, be accounted for by any model that in some way incorporates ion concentration dynamics [27–29, 33–57], provided that the ion concentration dynamics is accurately modeled. As we argued in the Introduction, previous multicompartmental neuron models that do incorporate ion concentration dynamics have not done it in a complete, ion conserving way that ensures a biophysically consistent relationship between ion concentration, electrical charge, and electrical potentials (see, e.g., [27, 48–57]). To specify, the change in the ion concentration in a given compartment will, in reality, depend on (i) the transmembrane influx of ions into this compartment, (ii) the diffusion of ions between this compartment and its neighboring compartment(s), and (iii) the electrical drift of ions between this compartment and its neighboring compartment(s). Some of the cited models account for only (i) [27, 49, 51], others account for (i) and (ii) [48, 50, 52–57], but neither account for (iii). When (iii) is not accounted for, electrical and diffusive processes are implicitly treated as independent processes, a simplifying assumption which is also incorporated in the reaction-diffusion module [71] in the NEURON simulation environment [72]. In models that apply this assumption, there will therefore be drift currents (along axons and dendrites) that affect *ϕ*_{m} (through the cable equation), but not the ion concentration dynamics, although they should, since also the drift currents are mediated by ions.

Here, we use simulations on the edPR model to test the inaccuracy introduced when not accounting for the effect of drift currents on ion concentration dynamics. We do so by comparing how many ions that were transferred from the somatic to the dendritic compartment through the intracellular (Fig 9A) and extracellular (Fig 9B) space, due to ionic diffusion (orange curves) versus electrical drift (blue curves), throughout the simulation in Fig 3. We note that Fig 9 shows the accumulatively moved number of ions (from time zero to *t*) due to axial fluxes exclusively. That is, the large number of, for example, Na^{+} ions transported intracellularly from the dendrite to the soma (negative sign) in Fig 9A1, does not by necessity mean that Na^{+} ions were piling up in the soma compartment, as the membrane efflux of Na^{+} was not accounted for in the figure.

Although diffusion tended to dominate the intracellular transport of ions on the long time scale (Fig 9A1–9A4), the transport due to electrical drift was not vanishingly small. For example, the number of K^{+} and Cl^{−} ions transported by electrical drift was at the end of the stimulus period (*t* = 20 s) about 35% of the transport due to diffusion for both species. In the extracellular space, diffusion was the clearly dominant transporter of Na^{+} and K^{+} (Fig 9B1 and 9B2), while diffusion and electrical drift were of comparable magnitude for the other ion species (Fig 9B3–9A4). Of course, these estimates are all specific to the edPR model, as they will depend strongly on the included ion channels, ion pumps and cotransporters, and on how they are distributed between the soma and dendrite. In general, however, the simulations in Fig 9 suggest that electrical drift is likely to have a non-negligible effect on ion concentration dynamics, and that ignoring this effect will give rise to rather inaccurate estimates.

Finally, we also converted the sum of ionic fluxes in Fig 9 into an effective current, represented as the number of transported unit charges, e^{+} (Fig 9A5–9B5). Interestingly, diffusion and drift contributed almost equally to the axial charge transport in the system. We note, however, that the movement of charges per time unit is indicated by the slope of the curves, which was much larger for the drift case (blue curve) than for diffusion (orange curve). The drift curve had a jigsaw shape, which shows that drift was moving charges back and forth in the system, while the diffusion always went in the same direction, explaining why it, despite being smaller than the drift current, had a comparably large accumulative effect on charge transport. The temporally averaged picture of charge transport that emerges from Fig 9A5 is that of a slow current loop where charge is transferred intracellularly from the soma to the dendrite (Fig 9A5), where it crosses the membrane (outward current), and then is transferred extracellularly back from the dendrite to the soma (Fig 9B5), before crossing the membrane again (inward current). This configuration is similar to the slow loop current seen during spatial buffering by astrocytes (see, e.g. Fig 1 in [67]).

In the previous section, we investigated the consequences of neglecting (iii) the contribution of drift currents on ion concentration dynamics. Here, we investigate the consequences of neglecting the effect of ionic diffusion (along dendrites) on the electrical potential, focusing on the extracellular potential *ϕ*_{e} . Forward modeling of extracellular potentials is typically based on volume conductor (VC) theory [16–18, 20, 21], which assumes that diffusive effects on electrical potentials are negligible. Being based on a unified electrodiffusive KNP framework (Fig 1), the edPR model accounts for the effects of ionic diffusion on the electrical potentials, and can thus be used to address the validity of this assumption.

To illustrate the effect of diffusion on *ϕ*_{e}, we may split it into a component *ϕ*_{VC,e} explained by standard VC-theory, and a component *ϕ*_{diff,e} representing the additional contribution caused by diffusive currents (Eq 81). In the simulation in Fig 3, the diffusive contribution was found to be very small compared to the VC-component (Fig 10). However, while *ϕ*_{VC,e} fluctuated rapidly from negative to positive values during neuronal activity, *ϕ*_{diff,e} varied on a slower time scale and had the same directionality throughout the simulation. This is equivalent to what we saw in Fig 9B5, i.e., that diffusion always moved charge in the same direction. Moreover, if we take the temporal averages of the potentials over the time series in Fig 10A, we find that they are -0.0023 mV, 0.0037 mV, and -0.0060 mV for *ϕ*_{e}, *ϕ*_{diff,e}, and *ϕ*_{VC,e} , respectively. This shows that the average diffusion- and VC-components of the total potential were of the same order of magnitude. As we also have demonstrated in previous studies, diffusion is thus likely to be important for the low-frequency components of extracellular potentials [31, 32, 73, 74]. Albeit small, the slowly varying diffusion evoked shifts in *ϕ*_{e} are putatively important for explaining the direct-current (DC) like shifts and long-time concentration dynamics reported during, e.g., spreading depression [25, 26].

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.

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

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.

Parameter | Value |
---|---|

Δx (distance between the two compartments) | 667⋅10^{−6} m |

A_{s} (somatic membrane area) | 616⋅10^{−12} m^{2} * |

A_{d} (dendritic membrane area) | 616⋅10^{−12} m^{2} * |

A_{i} (intracellular cross-section area) | α · A_{s} ^{†} |

A_{e} (extracellular cross-section area) | A_{i}/2 |

V_{si} (somatic intracellular volume) | 1437⋅10^{−18} m^{3} * |

V_{di} (dendritic intracellular volume) | 1437⋅10^{−18} m^{3} * |

V_{se} (somatic extracellular volume) | 718.5⋅10^{−18} m^{3} * |

V_{de} (dendritic extracellular volume) | 718.5⋅10^{−18} m^{3} * |

* The intracellular volumes (*V*_{si}, *V*_{di}) and membrane areas (*A*_{s}, *A*_{d}) correspond to spheres with radius 7 μm.

Two kinds of fluxes are involved: transmembrane fluxes and intra- and extracellular fluxes. The framework is general to the choice of the transmembrane fluxes. A transmembrane flux of ion species k (*j*_{k,m}) represents the sum of all fluxes through all membrane mechanisms that allow ion k to cross the membrane.

Intracellular flux densities are described by the Nernst-Planck equation:

$\begin{array}{c}\hfill {j}_{\mathrm{k},\mathrm{i}}=-\frac{{D}_{\mathrm{k}}}{{\lambda}_{\mathrm{i}}^{2}}\frac{{\gamma}_{\mathrm{k}}({\left[\mathrm{k}\right]}_{\text{di}}-{\left[\mathrm{k}\right]}_{\text{si}})}{\Delta x}-\frac{{D}_{\mathrm{k}}{z}_{\mathrm{k}}F}{{\lambda}_{\mathrm{i}}^{2}RT}{\overline{\left[\mathrm{k}\right]}}_{\mathrm{i}}\frac{{\varphi}_{\text{di}}-{\varphi}_{\text{si}}}{\Delta x}.\end{array}$

In Eq 1, *D*_{k} is the diffusion constant, *γ*_{k} is the fraction of freely moving ions, that is, ions that are not buffered or taken up by the ER, λ_{i} is the tortuosity, which represents the slowing down of diffusion due to obstacles, *γ*_{k}([k]_{di}−[k]_{si})/Δx is the axial concentration gradient, *z*_{k} is the charge number of ion species k, *F* is the Faraday constant, *R* is the gas constant, *T* is the absolute temperature, ${\overline{\left[\mathrm{k}\right]}}_{\mathrm{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

$\begin{array}{c}\hfill {j}_{\mathrm{k},\mathrm{e}}=-\frac{{D}_{\mathrm{k}}}{{\lambda}_{\mathrm{e}}^{2}}\frac{{\left[\mathrm{k}\right]}_{\text{de}}-{\left[\mathrm{k}\right]}_{\text{se}}}{\Delta x}-\frac{{D}_{\mathrm{k}}{z}_{\mathrm{k}}F}{{\lambda}_{\mathrm{e}}^{2}RT}{\overline{\left[\mathrm{k}\right]}}_{\mathrm{e}}\frac{{\varphi}_{\text{de}}-{\varphi}_{\text{se}}}{\Delta x}.\end{array}$

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.

Parameter | Value | Reference |
---|---|---|

D_{Na} (Na^{+} diffusion constant) | 1.33⋅10^{−9} m^{2}/s | [31, 81] |

D_{k} (K^{+} diffusion constant) | 1.96⋅10^{−9} m^{2}/s | [31, 81] |

D_{Cl} (Cl^{−} diffusion constant) | 2.03⋅10^{−9} m^{2}/s | [31, 81] |

D_{Ca} (Ca^{2+} diffusion constant) | 0.71⋅10^{−9} m^{2}/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 |

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:

$\begin{array}{c}\hfill \frac{d{\left[\mathrm{k}\right]}_{\text{si}}}{dt}=-{j}_{\mathrm{k},\text{sm}}\xb7\frac{{A}_{\mathrm{s}}}{{V}_{\text{si}}}-{j}_{\mathrm{k},\mathrm{i}}\xb7\frac{{A}_{\mathrm{i}}}{{V}_{\text{si}}},\end{array}$

$\begin{array}{c}\hfill \frac{d{\left[\mathrm{k}\right]}_{\text{di}}}{dt}=-{j}_{\mathrm{k},\text{dm}}\xb7\frac{{A}_{\mathrm{d}}}{{V}_{\text{di}}}+{j}_{\mathrm{k},\mathrm{i}}\xb7\frac{{A}_{\mathrm{i}}}{{V}_{\text{di}}},\end{array}$

$\begin{array}{c}\hfill \frac{d{\left[\mathrm{k}\right]}_{\text{se}}}{dt}=+{j}_{\mathrm{k},\text{sm}}\xb7\frac{{A}_{\mathrm{s}}}{{V}_{\text{se}}}-{j}_{\mathrm{k},\mathrm{e}}\xb7\frac{{A}_{\mathrm{e}}}{{V}_{\text{se}}},\end{array}$

$\begin{array}{c}\hfill \frac{d{\left[\mathrm{k}\right]}_{\text{de}}}{dt}=+{j}_{\mathrm{k},\text{dm}}\xb7\frac{{A}_{\mathrm{d}}}{{V}_{\text{de}}}+{j}_{\mathrm{k},\mathrm{e}}\xb7\frac{{A}_{\mathrm{e}}}{{V}_{\text{de}}}.\end{array}$

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:
$\begin{array}{c}\hfill \frac{d{\left[\mathrm{k}\right]}_{\text{si}}}{dt}=-{j}_{\mathrm{k},\text{sm}}\xb7\frac{{A}_{\mathrm{s}}}{{V}_{\text{si}}}+\frac{{A}_{\mathrm{i}}{D}_{\mathrm{k}}}{{V}_{\text{si}}{\lambda}_{\mathrm{i}}^{2}\Delta x}[{\gamma}_{\mathrm{k}}({\left[\mathrm{k}\right]}_{\text{di}}-{\left[\mathrm{k}\right]}_{\text{si}})+\frac{{z}_{\mathrm{k}}F}{RT}{\overline{\left[\mathrm{k}\right]}}_{\mathrm{i}}({\varphi}_{\text{di}}-{\varphi}_{\text{si}})],\end{array}$

where the voltage variables so far are undefined.If we have four ion species (Na^{+}, K^{+}, Cl^{−}, and Ca^{2+}) in four compartments, we have 20 unknown parameters (16 for [k] and four for *ϕ* ), while Eqs 3–6 for four ion species give us only 16 equations. To solve this, we need to define additional constraints that allow us to express the potentials *ϕ* in terms of ion concentrations.

As we may define an arbitrary reference point for *ϕ*, we take

$\begin{array}{c}\hfill {\varphi}_{\text{de}}=0,\end{array}$

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

$\begin{array}{c}\hfill {\varphi}_{\mathrm{m}}=Q/{C}_{\mathrm{m}}.\end{array}$

Here, $\begin{array}{c}\hfill {\varphi}_{\text{dm}}={\varphi}_{\text{di}}={Q}_{\text{di}}/{C}_{\mathrm{m}},\end{array}$

in the dendrite (since $\begin{array}{c}\hfill {\varphi}_{\text{sm}}={\varphi}_{\text{si}}-{\varphi}_{\text{se}}={Q}_{\text{si}}/{C}_{\mathrm{m}},\end{array}$

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 ${Q}_{\text{di}}=F{\displaystyle \sum _{\mathrm{k}}}{z}_{\mathrm{k}}{\left[\mathrm{k}\right]}_{\text{di}}{V}_{\text{di}}$. By inserting this into Eq 10, we obtain the final expression for *ϕ*_{di}:

$\begin{array}{c}\hfill {\varphi}_{\text{di}}=\left(F{\displaystyle \sum _{\mathrm{k}}}{z}_{\mathrm{k}}{\left[\mathrm{k}\right]}_{\text{di}}{V}_{\text{di}}\right)/\left({c}_{\mathrm{m}}{A}_{\mathrm{d}}\right).\end{array}$

By inserting the equivalent expression for $\begin{array}{c}\hfill {\varphi}_{\text{si}}-{\varphi}_{\text{se}}={Q}_{\text{si}}/{C}_{\mathrm{m}}=\left(F{\displaystyle \sum _{\mathrm{k}}}{z}_{\mathrm{k}}{\left[\mathrm{k}\right]}_{\text{si}}{V}_{\text{si}}\right)/\left({c}_{\mathrm{m}}{A}_{\mathrm{s}}\right).\end{array}$

Here, the extracellular potential is not set to zero, so we need a fourth constraint to determine The fourth and final constraint is that (iv) we must ensure charge anti-symmetry. For the charge anti-symmetry between the two sides of the capacitive membrane (*Q*_{i} = −*Q*_{e}) to be preserved in time, we must define our initial conditions so that this is the case at *t* = 0, and the system dynamics so that this stays the case. Hence, the system dynamics must ensure that *dQ*_{di}/*dt* = −*dQ*_{de}/*dt* and *dQ*_{si}/*dt* = −*dQ*_{se}/*dt*. The membrane currents (in isolation) will always fulfill this criterion, as any charge that crosses the membrane by definition disappears from one side of it and pops up at the other. Hence, we thus need to make sure that also the axial currents (in isolation) fulfill the criterion. The system must thus be constrained so that, if an extracellular current transports a charge *δq* into a given extracellular compartment, the intracellular current must transport the opposite charge −*δq* into the adjoint intracellular compartment. That is, we must have that:

$\begin{array}{c}\hfill {A}_{\mathrm{i}}{i}_{\mathrm{i}}=-{A}_{\mathrm{e}}{i}_{\mathrm{e}},\end{array}$

where $\begin{array}{c}\hfill {i}_{\mathrm{i}}=-\frac{F}{{\lambda}_{\mathrm{i}}^{2}\Delta x}\sum _{\mathrm{k}}{D}_{\mathrm{k}}{z}_{\mathrm{k}}{\gamma}_{\mathrm{k}}({\left[\mathrm{k}\right]}_{\text{di}}-{\left[\mathrm{k}\right]}_{\text{si}})-\frac{{F}^{2}}{RT{\lambda}_{\mathrm{i}}^{2}\Delta x}\sum _{\mathrm{k}}{D}_{\mathrm{k}}{z}_{\mathrm{k}}^{2}{\overline{\left[\mathrm{k}\right]}}_{\mathrm{i}}({\varphi}_{\text{di}}-{\varphi}_{\text{si}}),\end{array}$

$\begin{array}{c}\hfill {i}_{\mathrm{e}}=-\frac{F}{{\lambda}_{\mathrm{e}}^{2}\Delta x}\sum _{\mathrm{k}}{D}_{\mathrm{k}}{z}_{\mathrm{k}}({\left[\mathrm{k}\right]}_{\text{de}}-{\left[\mathrm{k}\right]}_{\text{se}})-\frac{{F}^{2}}{RT{\lambda}_{\mathrm{e}}^{2}\Delta x}\sum _{\mathrm{k}}{D}_{\mathrm{k}}{z}_{\mathrm{k}}^{2}{\overline{\left[\mathrm{k}\right]}}_{\mathrm{e}}({\varphi}_{\text{de}}-{\varphi}_{\text{se}}).\end{array}$

In Eq 15, the first term is the diffusion current density:

$\begin{array}{c}\hfill {i}_{\text{diff},\mathrm{i}}=-\frac{F}{{\lambda}_{\mathrm{i}}^{2}\Delta x}\sum _{\mathrm{k}}{D}_{\mathrm{k}}{z}_{\mathrm{k}}{\gamma}_{\mathrm{k}}({\left[\mathrm{k}\right]}_{\text{di}}-{\left[\mathrm{k}\right]}_{\text{si}}),\end{array}$

which is defined by the ion concentrations. The second term is the field driven current density
$\begin{array}{c}\hfill {i}_{\text{field},\mathrm{i}}=-{\sigma}_{i}\frac{({\varphi}_{\text{di}}-{\varphi}_{\text{si}})}{\Delta x},\end{array}$

where we have identified the conductivity as
$\begin{array}{c}\hfill {\sigma}_{\mathrm{i}}=\frac{{F}^{2}}{RT{\lambda}_{\mathrm{i}}^{2}}\sum _{\mathrm{k}}{D}_{\mathrm{k}}{z}_{\mathrm{k}}^{2}{\overline{\left[\mathrm{k}\right]}}_{\mathrm{i}}.\end{array}$

Similarly, Eq 16 can be written in terms of *i*_{diff,e}, *i*_{field,e}, and *σ*_{e} . By combining Eqs 14, 15 and 16, we obtain:

$\begin{array}{c}\hfill -{A}_{\mathrm{i}}{i}_{\text{diff},\mathrm{i}}+{A}_{\mathrm{i}}{\sigma}_{\mathrm{i}}\xb7\frac{({\varphi}_{\text{di}}-{\varphi}_{\text{si}})}{\Delta x}={A}_{\mathrm{e}}{i}_{\text{diff},\mathrm{e}}-{A}_{\mathrm{e}}{\sigma}_{\mathrm{e}}\xb7\frac{({\varphi}_{\text{de}}-{\varphi}_{\text{se}})}{\Delta x}.\end{array}$

In Eq 20, *ϕ*_{di} and *ϕ*_{de} are already known from Eqs 8 and 12, while *i*_{diff} and *σ* are expressed in terms of ion concentrations. We may thus solve Eqs 13 and 20 for the last two voltage variables *ϕ*_{se} and *ϕ*_{si}:

$\begin{array}{c}\hfill {\varphi}_{\text{se}}=({\varphi}_{\text{di}}-\frac{\Delta x}{{\sigma}_{\mathrm{i}}}\xb7{i}_{\text{diff},\mathrm{i}}-\frac{{A}_{\mathrm{e}}\Delta x}{{A}_{\mathrm{i}}{\sigma}_{\mathrm{i}}}\xb7{i}_{\text{diff},\mathrm{e}}-\frac{{Q}_{\text{si}}}{{c}_{\mathrm{m}}{A}_{\mathrm{s}}})/(1+\frac{{A}_{\mathrm{e}}{\sigma}_{\mathrm{e}}}{{A}_{\mathrm{i}}{\sigma}_{\mathrm{i}}}),\end{array}$

$\begin{array}{c}\hfill {\varphi}_{\text{si}}=\frac{{Q}_{\text{si}}}{{c}_{\mathrm{m}}{A}_{\mathrm{s}}}+{\varphi}_{\text{se}}.\end{array}$

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

$\begin{array}{c}\hfill {j}_{\mathrm{k},\text{leak}}={\overline{g}}_{\mathrm{k},\text{leak}}({\varphi}_{\mathrm{m}}-{E}_{\mathrm{k}})/\left(F{z}_{\mathrm{k}}\right),\end{array}$

where ${\overline{g}}_{\mathrm{k},\text{leak}}$ is the ion conductance, $\begin{array}{c}\hfill {E}_{\mathrm{k}}=\frac{RT}{{z}_{\mathrm{k}}F}ln\frac{{\left[\mathrm{k}\right]}_{\mathrm{e}}}{{\gamma}_{\mathrm{k}}{\left[\mathrm{k}\right]}_{\mathrm{i}}}.\end{array}$

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.

All active ion channel currents were adopted from the original PR model [3], as they were described in [8], and converted to ion channel fluxes. The soma compartment contained a Na^{+} flux (*j*_{Na}) and a K^{+} delayed rectifier flux (*j*_{K−DR}), while the dendrite contained a voltage-dependent Ca^{2+} flux (*j*_{Ca}), a voltage-dependent K^{+} AHP flux (*j*_{K−AHP}), and a Ca^{2+}-dependent K^{+} flux (*j*_{K−C}):

$\begin{array}{c}\hfill {j}_{\text{Na}}={g}_{\text{Na}}({\varphi}_{\text{sm}}-{E}_{\text{Na},\mathrm{s}})/\left(F{z}_{\text{Na}}\right),\end{array}$

$\begin{array}{c}\hfill {j}_{\mathrm{K}-\text{DR}}={g}_{\text{DR}}({\varphi}_{\text{sm}}-{E}_{\mathrm{K},\mathrm{s}})/\left(F{z}_{\mathrm{K}}\right),\end{array}$

$\begin{array}{c}\hfill {j}_{\text{Ca}}={g}_{\text{Ca}}({\varphi}_{\text{dm}}-{E}_{\text{Ca},\mathrm{d}})/\left(F{z}_{\text{Ca}}\right),\end{array}$

$\begin{array}{c}\hfill {j}_{\mathrm{K}-\text{AHP}}={g}_{\text{AHP}}({\varphi}_{\text{dm}}-{E}_{\mathrm{K},\mathrm{d}})/\left(F{z}_{\mathrm{K}}\right),\end{array}$

$\begin{array}{c}\hfill {j}_{\mathrm{K}-\mathrm{C}}={g}_{\mathrm{C}}({\varphi}_{\text{dm}}-{E}_{\mathrm{K},\mathrm{d}})/\left(F{z}_{\mathrm{K}}\right).\end{array}$

The voltage-dependent conductances were modeled using the Hodkin-Huxley formalism with differential equations for the gating variables:

$\begin{array}{c}\hfill \frac{dx}{dt}={\alpha}_{\mathrm{x}}(1-x)-{\beta}_{\mathrm{x}}x,\phantom{\rule{10pt}{0ex}}\text{with}\phantom{\rule{4pt}{0ex}}x=m,h,n,s,c,q,\end{array}$

$\begin{array}{c}\hfill \frac{dz}{dt}=\frac{{z}_{\infty}-z}{{\tau}_{\mathrm{z}}}.\end{array}$

The conductances and gating variables were given by:

$\begin{array}{c}\hfill {g}_{\text{Na}}={\overline{g}}_{\text{Na}}{m}_{\infty}^{2}h,\end{array}$

$\begin{array}{c}\hfill {g}_{\text{DR}}={\overline{g}}_{\text{DR}}n,\end{array}$

$\begin{array}{c}\hfill {g}_{\text{Ca}}={\overline{g}}_{\text{Ca}}{s}^{2}z,\end{array}$

$\begin{array}{c}\hfill {g}_{\mathrm{C}}={\overline{g}}_{\mathrm{C}}c\chi \left(\left[{\text{Ca}}^{2+}\right]\right),\end{array}$

$\begin{array}{c}\hfill {g}_{\text{AHP}}={\overline{g}}_{\text{AHP}}q,\end{array}$

$\begin{array}{c}\hfill {\alpha}_{\mathrm{m}}=-\frac{3.2\xb7{10}^{5}\xb7{\varphi}_{1}}{exp(-{\varphi}_{1}/0.004)-1},\phantom{\rule{0.277778em}{0ex}}\text{with}\phantom{\rule{4pt}{0ex}}{\varphi}_{1}={\varphi}_{\mathrm{m}}+0.0469\end{array}$

$\begin{array}{c}\hfill {\beta}_{\mathrm{m}}=\frac{2.8\xb7{10}^{5}\xb7{\varphi}_{2}}{exp({\varphi}_{2}/0.005)-1},\phantom{\rule{0.277778em}{0ex}}\text{with}\phantom{\rule{4pt}{0ex}}{\varphi}_{2}={\varphi}_{\mathrm{m}}+0.0199\end{array}$

$\begin{array}{c}\hfill {m}_{\infty}=\frac{{\alpha}_{\mathrm{m}}}{{\alpha}_{\mathrm{m}}+{\beta}_{\mathrm{m}}}\end{array}$

$\begin{array}{c}\hfill {\alpha}_{\mathrm{h}}=128exp\frac{-0.043-{\varphi}_{\mathrm{m}}}{0.018},\end{array}$

$\begin{array}{c}\hfill {\beta}_{\mathrm{h}}=\frac{4000}{1+exp(-{\varphi}_{3}/0.005)},\phantom{\rule{0.277778em}{0ex}}\text{with}\phantom{\rule{4pt}{0ex}}{\varphi}_{3}={\varphi}_{\mathrm{m}}+0.02\end{array}$

$\begin{array}{c}\hfill {\alpha}_{\mathrm{n}}=-\frac{1.6\xb7{10}^{4}\xb7{\varphi}_{4}}{exp(-{\varphi}_{4}/0.005)-1},\phantom{\rule{0.277778em}{0ex}}\text{with}\phantom{\rule{4pt}{0ex}}{\varphi}_{4}={\varphi}_{\mathrm{m}}+0.0249\end{array}$

$\begin{array}{c}\hfill {\beta}_{\mathrm{n}}=250exp(-{\varphi}_{5}/0.04),\phantom{\rule{0.277778em}{0ex}}\text{with}\phantom{\rule{4pt}{0ex}}{\varphi}_{5}={\varphi}_{\mathrm{m}}+0.04\end{array}$

$\begin{array}{c}\hfill {\alpha}_{\mathrm{s}}=\frac{1600}{1+exp(-72({\varphi}_{\mathrm{m}}-0.005))},\end{array}$

$\begin{array}{c}\hfill {\beta}_{\mathrm{s}}=\frac{2\xb7{10}^{4}\xb7{\varphi}_{6}}{exp({\varphi}_{6}/0.005)-1},\phantom{\rule{0.277778em}{0ex}}\text{with}\phantom{\rule{4pt}{0ex}}{\varphi}_{6}={\varphi}_{\mathrm{m}}+0.0089\end{array}$

$\begin{array}{c}\hfill {z}_{\infty}=\frac{1}{1+exp({\varphi}_{7}/0.001)},\phantom{\rule{0.277778em}{0ex}}\text{with}\phantom{\rule{4pt}{0ex}}{\varphi}_{7}={\varphi}_{\mathrm{m}}+0.03\end{array}$

$\begin{array}{c}\hfill {\tau}_{\mathrm{z}}=1,\end{array}$

$\begin{array}{c}\hfill {\alpha}_{\mathrm{c}}=\{\begin{array}{cc}52.7exp(\frac{{\varphi}_{8}}{0.011}-\frac{{\varphi}_{9}}{0.027}),\hfill & \text{if}\phantom{\rule{4pt}{0ex}}{\varphi}_{\mathrm{m}}\le -0.01\phantom{\rule{4pt}{0ex}}\text{V}\hfill \\ 2000exp(-{\varphi}_{9}/0.027),\hfill & \text{otherwise}\hfill \end{array}\end{array}$

$\begin{array}{cc}& \text{with}\phantom{\rule{4pt}{0ex}}{\varphi}_{8}={\varphi}_{\mathrm{m}}+0.05\phantom{\rule{0.277778em}{0ex}}\text{and}\phantom{\rule{4pt}{0ex}}{\varphi}_{9}={\varphi}_{\mathrm{m}}+0.0535\end{array}$

$\begin{array}{c}\hfill {\beta}_{\mathrm{c}}=\{\begin{array}{cc}2000exp(-{\varphi}_{9}/0.027)-{\alpha}_{\mathrm{c}},\hfill & \text{if}\phantom{\rule{4pt}{0ex}}{\varphi}_{\mathrm{m}}\le -0.01\phantom{\rule{4pt}{0ex}}\text{V}\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}\end{array}$

$\begin{array}{c}\hfill \chi =min(\frac{{\gamma}_{\text{Ca}}\left[{\text{Ca}}^{2+}\right]-99.8\xb7{10}^{-6}}{2.5\xb7{10}^{-4}},1),\end{array}$

$\begin{array}{c}\hfill {\alpha}_{\mathrm{q}}=min(2\xb7{10}^{4}({\gamma}_{\text{Ca}}\left[{\text{Ca}}^{2+}\right]-99.8\xb7{10}^{-6}),10),\end{array}$

$\begin{array}{c}\hfill {\beta}_{\mathrm{q}}=1.\end{array}$

All these equations were taken from [8] (with errata [83]) and converted so that values are given in SI units: units for rates (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].

$\begin{array}{c}\hfill {j}_{\text{pump}}=\frac{\rho}{1.0+exp((25-{\left[{\text{Na}}^{+}\right]}_{\mathrm{i}})/3)}\xb7\frac{1.0}{1.0+exp(3.5-{\left[{\mathrm{K}}^{+}\right]}_{\mathrm{e}})},\end{array}$

$\begin{array}{c}\hfill {j}_{\text{kcc}2}={U}_{\text{kcc}2}ln\left(\frac{{\left[{\mathrm{K}}^{+}\right]}_{\mathrm{i}}{\left[{\text{Cl}}^{-}\right]}_{\mathrm{i}}}{{\left[{\mathrm{K}}^{+}\right]}_{\mathrm{e}}{\left[{\text{Cl}}^{-}\right]}_{\mathrm{e}}}\right),\end{array}$

$\begin{array}{c}\hfill {j}_{\text{nkcc}1}={U}_{\text{nkcc}1}f\left({\left[{\mathrm{K}}^{+}\right]}_{\mathrm{e}}\right)(ln(\frac{{\left[{\mathrm{K}}^{+}\right]}_{\mathrm{i}}{\left[{\text{Cl}}^{-}\right]}_{\mathrm{i}}}{{\left[{\mathrm{K}}^{+}\right]}_{\mathrm{e}}{\left[{\text{Cl}}^{-}\right]}_{\mathrm{e}}})+ln(\frac{{\left[{\text{Na}}^{+}\right]}_{\mathrm{i}}{\left[{\text{Cl}}^{-}\right]}_{\mathrm{i}}}{{\left[{\text{Na}}^{+}\right]}_{\mathrm{e}}{\left[{\text{Cl}}^{-}\right]}_{\mathrm{e}}}\left)\right),\end{array}$

$\begin{array}{c}\hfill f\left({\left[{\mathrm{K}}^{+}\right]}_{\mathrm{e}}\right)=\frac{1}{1+exp(16-{\left[{\mathrm{K}}^{+}\right]}_{\mathrm{e}})},\end{array}$

where Intracellular Ca^{2+} decay was modeled in a similar fashion as in [3], but to ensure ion conservation we modeled it as an electroneutral Ca^{2+}/2Na^{+} exchanger, exchanging one Ca^{2+} (outward) for two Na^{+} (inward). Putatively, this phenomenological model for the decay could represent the joint effect of several mechanisms in a real system, such as the Ca^{2+}/3Na^{+} exchanger, a Ca^{2+} leakage current, SERCA pumps, etc. The decay flux density was defined as:

$\begin{array}{c}\hfill {j}_{\text{Ca}-\text{dec}}={U}_{\text{Ca}-\text{dec}}({\left[{\text{Ca}}^{2+}\right]}_{\mathrm{i}}-{\left[{\text{Ca}}^{2+}\right]}_{\mathrm{i},\mathrm{b}})\xb7\frac{{V}_{\mathrm{i}}}{{A}_{\mathrm{m}}}\end{array}$

where We summarize the model here for easy reference. In short, we solved four differential equations for all ion species k:

$\begin{array}{c}\hfill \frac{d{\left[\mathrm{k}\right]}_{\text{si}}}{dt}=-{j}_{\mathrm{k},\text{sm}}\xb7\frac{{A}_{\mathrm{s}}}{{V}_{\text{si}}}-{j}_{\mathrm{k},\mathrm{i}}\xb7\frac{{A}_{\mathrm{i}}}{{V}_{\text{si}}},\end{array}$

$\begin{array}{c}\hfill \frac{d{\left[\mathrm{k}\right]}_{\text{di}}}{dt}=-{j}_{\mathrm{k},\text{dm}}\xb7\frac{{A}_{\mathrm{d}}}{{V}_{\text{di}}}+{j}_{\mathrm{k},\mathrm{i}}\xb7\frac{{A}_{\mathrm{i}}}{{V}_{\text{di}}},\end{array}$

$\begin{array}{c}\hfill \frac{d{\left[\mathrm{k}\right]}_{\text{se}}}{dt}=+{j}_{\mathrm{k},\text{sm}}\xb7\frac{{A}_{\mathrm{s}}}{{V}_{\text{se}}}-{j}_{\mathrm{k},\mathrm{e}}\xb7\frac{{A}_{\mathrm{e}}}{{V}_{\text{se}}},\end{array}$

$\begin{array}{c}\hfill \frac{d{\left[\mathrm{k}\right]}_{\text{de}}}{dt}=+{j}_{\mathrm{k},\text{dm}}\xb7\frac{{A}_{\mathrm{d}}}{{V}_{\text{se}}}+{j}_{\mathrm{k},\mathrm{e}}\xb7\frac{{A}_{\mathrm{e}}}{{V}_{\text{de}}}.\end{array}$

At each time step, $\begin{array}{c}\hfill {\varphi}_{\text{de}}=0,\end{array}$

$\begin{array}{c}\hfill {\varphi}_{\text{di}}={Q}_{\text{di}}/\left({c}_{\mathrm{m}}{A}_{\mathrm{d}}\right)\end{array}$

$\begin{array}{c}\hfill {\varphi}_{\text{se}}=({\varphi}_{\text{di}}-\frac{\Delta x}{{\sigma}_{\mathrm{i}}}\xb7{i}_{\text{diff},\mathrm{i}}-\frac{{A}_{\mathrm{e}}\Delta x}{{A}_{\mathrm{i}}{\sigma}_{\mathrm{i}}}\xb7{i}_{\text{diff},\mathrm{e}}-\frac{{Q}_{\text{si}}}{{c}_{\mathrm{m}}{A}_{\mathrm{s}}})/(1+\frac{{A}_{\mathrm{e}}{\sigma}_{\mathrm{e}}}{{A}_{\mathrm{i}}{\sigma}_{\mathrm{i}}}),\end{array}$

$\begin{array}{c}\hfill {\varphi}_{\text{si}}=\frac{{Q}_{\text{si}}}{{c}_{\mathrm{m}}{A}_{\mathrm{s}}}+{\varphi}_{\text{se}}.\end{array}$

The total membrane flux densities were as follows:

$\begin{array}{c}\hfill {j}_{\text{Na},\text{sm}}={j}_{\text{Na}}+{j}_{\text{Na},\text{leak}}+3{j}_{\text{pump}}+{j}_{\text{nkcc}1}-2{j}_{\text{Ca}-\text{dec}},\end{array}$

$\begin{array}{c}\hfill {j}_{\mathrm{K},\text{sm}}={j}_{\mathrm{K}-\text{DR}}+{j}_{\mathrm{K},\text{leak}}-2{j}_{\text{pump}}+{j}_{\text{nkcc}1}+{j}_{\text{kcc}2},\end{array}$

$\begin{array}{c}\hfill {j}_{\text{Cl},\text{sm}}={j}_{\text{Cl},\text{leak}}+2{j}_{\text{nkcc}1}+{j}_{\text{kcc}2},\end{array}$

$\begin{array}{c}\hfill {j}_{\text{Ca},\text{sm}}={j}_{\text{Ca}-\text{dec}},\end{array}$

$\begin{array}{c}\hfill {j}_{\text{Na},\text{dm}}={j}_{\text{Na},\text{leak}}+3{j}_{\text{pump}}+{j}_{\text{nkcc}1}-2{j}_{\text{Ca}-\text{dec}},\end{array}$

$\begin{array}{c}\hfill {j}_{\mathrm{K},\text{dm}}={j}_{\mathrm{K}-\text{AHP}}+{j}_{\mathrm{K}-\mathrm{C}}+{j}_{\mathrm{K},\text{leak}}-2{j}_{\text{pump}}+{j}_{\text{nkcc}1}+{j}_{\text{kcc}2},\end{array}$

$\begin{array}{c}\hfill {j}_{\text{Cl},\text{dm}}={j}_{\text{Cl},\text{leak}}+2{j}_{\text{nkcc}1}+{j}_{\text{kcc}2},\end{array}$

$\begin{array}{c}\hfill {j}_{\text{Ca},\text{dm}}={j}_{\text{Ca}}+{j}_{\text{Ca}-\text{dec}}.\end{array}$

Fig 2 summarizes the model. The parameters involved in this model and their values used in this study are listed in Tables 1–4.Parameter | Value | Reference |
---|---|---|

c_{m} | 3⋅10^{−2} F/m^{2} | [3, 8] |

$g\xafNa,leak$ | 0.247 S/m^{2} | [45] |

$g\xafK,leak$ | 0.5 S/m^{2} | [45] |

$g\xafCl,leak$ | 1.0 S/m^{2} | [45] |

$g\xafNa$ | 300 S/m^{2} | [3, 8] |

$g\xafDR$ | 150 S/m^{2} | [3, 8] |

$g\xafCa$ | 118 S/m^{2} | |

$g\xafAHP$ | 8 S/m^{2} | [3, 8] |

$g\xafC$ | 150 S/m^{2} | [3, 8] |

ρ | 1.87⋅10^{−6} mol/(m^{2}s) | [45]* |

U_{kcc2} | 7.0⋅10^{−7} mol/(m^{2}s) | [45]* |

U_{nkcc1} | 2.33⋅10^{−7} mol/(m^{2}s) | [45]* |

U_{Ca−dec} | 75 s^{−1} | [3, 8] |

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 *n*_{0} = 0.001, *h*_{0} = 0.999, s_{0} = 0.009, *c*_{0} = 0.007, *q*_{0} = 0.01, and [Ca^{2+}]_{0} = 0.2, same as in [3].

The parameters listed in Tables 1–4 were used in all the simulations of the electrodiffusive Pinsky-Rinzel (edPR) model. We tuned the Ca^{2+} conductance ${\overline{g}}_{\text{Ca}}$ manually to obtain comparable spike shapes between the edPR model and the original PR model, as well as the fraction of free Ca^{2+} 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 *A*_{i}/Δx, and all model outputs depended on this ratio, and not on *A*_{i} or Δ*x* in isolation. By choice, we adjusted the coupling strength by varying *A*_{i} = *αA*_{m} 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 *g*_{c}, which had a default value of 10.5mS/cm^{2} . In Fig 5A, *g*_{c} was set to 2.26mS/cm^{2}.

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 *z*_{X} = −1 (assuming this to be the mean charge number of the real macromolecules) and diffusion constants *D*_{X} = 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:

$\begin{array}{c}\hfill {\left[{\mathrm{X}}^{-}\right]}_{\mathrm{i},0}={z}_{\text{Na}}{\left[{\text{Na}}^{+}\right]}_{\mathrm{i},0}+{z}_{\mathrm{K}}{\left[{\mathrm{K}}^{+}\right]}_{\mathrm{i},0}+{z}_{\text{Cl}}{\left[{\text{Cl}}^{-}\right]}_{\mathrm{i},0}+{z}_{\text{Ca}}{\left[{\text{Ca}}^{2+}\right]}_{\mathrm{i},0}-{\varphi}_{\mathrm{m},0}\frac{{c}_{\mathrm{m}}{A}_{\mathrm{m}}}{{V}_{\mathrm{i}}F},\end{array}$

$\begin{array}{c}\hfill {\left[{\mathrm{X}}^{-}\right]}_{\mathrm{e},0}={z}_{\text{Na}}{\left[{\text{Na}}^{+}\right]}_{\mathrm{e},0}+{z}_{\mathrm{K}}{\left[{\mathrm{K}}^{+}\right]}_{\mathrm{e},0}+{z}_{\text{Cl}}{\left[{\text{Cl}}^{-}\right]}_{\mathrm{e},0}+{z}_{\text{Ca}}{\left[{\text{Ca}}^{2+}\right]}_{\mathrm{e},0}+{\varphi}_{\mathrm{m},0}\frac{{c}_{\mathrm{m}}{A}_{\mathrm{m}}}{{V}_{\mathrm{e}}F}.\end{array}$

Variables | Pre-calibrated | Post-calibrated^{1} |
---|---|---|

ϕ_{m,0} ^{†} | -68 mV | -67.7 mV |

[Na^{+}]_{i,0} | 15 mM | 16.9 mM |

[Na^{+}]_{e,0} | 145 mM | 141.2 mM |

[K^{+}]_{i,0} | 140 mM | 139.5 mM |

[K^{+}]_{e,0} | 5 mM | 5.9 mM |

[Cl^{−}]_{i,0} | 4 mM | 5.4 mM |

[Cl^{−}]_{e,0} | 110 mM | 107.1 mM |

[Ca^{2+}]_{i,0} | 0.01 mM* | 0.01 mM* |

[Ca^{2+}]_{e,0} | 1.1 mM | 1.1 mM |

[X^{−}]_{i,0} ^{‡} | 151.0 mM | 151.0 mM |

[X^{−}]_{e,0} ^{‡} | 42.2 mM | 42.2 mM |

n_{0} | 0.001 | 0.0003 |

h_{0} | 0.999 | 0.999 |

s_{0} | 0.009 | 0.007 |

c_{0} | 0.007 | 0.005 |

q_{0} | 0.010 | 0.011 |

z_{0} | 1.0 | 1.0 |

* Only 1% of the total intracellular Ca^{2+}, that is, a 100 nM, was assumed to be free (unbuffered).

In the next step, we calibrated the model by running it for 1800 s to obtain the post-calibrated values of all the state variables (Table 5). These post-calibrated values were written to file and used as initial conditions in all simulations shown throughout this paper. For technical reasons, we did not read the constant residual concentrations, [*X*^{−}]_{i,0} and [*X*^{−}]_{e,0} , to/from file, but re-calculated them from Eqs 75 and 76 in the beginning of each simulation to minimize rounding errors and ensure strict electroneutrality. While the pre-calibrated initial conditions were identical in the somatic and dendritic compartment, the post-calibration values were not strictly identical, but identical up to the decimal place included in Table 5. Hence, the indicated values apply for both the soma and dendrite compartments. The post-calibrated values of the ion concentrations gave us the following reversal potentials: *E*_{Na} = 57mV, *E*_{K} = −84mV, *E*_{Cl} = −79mV, and *E*_{Ca} = 124mV.

The pre-calibration values for the ion concentrations were taken from Table 2.1 in [85], which lists the ranges of typical intra- and extracellular concentrations in mammalian neurons. From the ranges given in this table, we selected values that made the edPR model (throughout the calibration period) reside close to the selected initial membrane potential of -68 mV, a typical value found for CA3 pyramidal cells in adult rats [86].

We stimulated the cell by injecting a K^{+} current *i*_{stim} into the soma. Previous computational modeling of a cardiac cell has shown that stimulus with K^{+} causes the least physiological disruption [33]. To ensure ion conservation, we removed the same amount of K^{+} ions from the corresponding extracellular compartment:

$\begin{array}{c}\hfill \frac{d{\left[{\mathrm{K}}^{+}\right]}_{\text{si}}}{dt}+=\frac{{i}_{\text{stim}}}{F{z}_{\mathrm{K}}{V}_{\text{si}}},\end{array}$

$\begin{array}{c}\hfill \frac{d{\left[{\mathrm{K}}^{+}\right]}_{\text{se}}}{dt}-=\frac{{i}_{\text{stim}}}{F{z}_{\mathrm{K}}{V}_{\text{se}}}.\end{array}$

Fig 9: To calculate the accumulative transport of ion species k in the intracellular solution (from time zero to *t*) due to diffusion, we integrated *A*_{i}*N*_{A}*j*_{k,diff,i} from time zero to t, where *N*_{A} is the Avogadro constant. Similarly, we integrated *A*_{e}*N*_{A}*j*_{k,diff,e} to calculate the accumulative transport of ions in the extracellular solution due to diffusion. We did the same calculations with *j*_{k,drift} to study the accumulative transport of ions due to drift. When knowing the accumulative transport of each ion species, k_{akkum}, we calculated the total transport of e^{+} from their weighted sum:

$\begin{array}{c}\hfill {\mathrm{e}}_{\text{akkum}}^{+}={z}_{\text{Na}}{\text{Na}}_{\text{akkum}}^{+}+{z}_{\mathrm{K}}{\mathrm{K}}_{\text{akkum}}^{+}+{z}_{\text{Cl}}{\text{Cl}}_{\text{akkum}}^{-}+{z}_{\text{Ca}}{\text{Ca}}_{\text{akkum}}^{2+}.\end{array}$

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

$\begin{array}{c}\hfill {i}_{\mathrm{e}}={i}_{\text{diff},\mathrm{e}}+{i}_{\text{field},\mathrm{e}}={i}_{\text{diff},\mathrm{e}}+{\sigma}_{\mathrm{e}}\frac{{\varphi}_{\text{se}}}{\Delta x},\end{array}$

where the last equality follows when we insert Eq 18 for the extracellular field-driven current density $\begin{array}{c}\hfill {\varphi}_{\text{se}}={\varphi}_{\text{VC},\text{se}}+{\varphi}_{\text{diff},\text{se}},\end{array}$

where $\begin{array}{c}\hfill {i}_{\mathrm{e}}={i}_{\text{diff},\mathrm{e}}+{\sigma}_{\mathrm{e}}\frac{{\varphi}_{\text{VC},\text{se}}}{\Delta x}+{\sigma}_{\mathrm{e}}\frac{{\varphi}_{\text{diff},\text{se}}}{\Delta x}.\end{array}$

We may split this into two equations if we recognize that
$\begin{array}{c}\hfill {i}_{\mathrm{e}}={\sigma}_{\mathrm{e}}\frac{{\varphi}_{\text{VC},\text{se}}}{\Delta x},\end{array}$

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 $\begin{array}{c}\hfill {i}_{\text{diff},\mathrm{e}}=-{\sigma}_{\mathrm{e}}\frac{{\varphi}_{\text{diff},\text{se}}}{\Delta x}.\end{array}$

Since we already knew *i*_{e} and *i*_{diff,e} from simulations with the KNP framework, we used Eqs 83 and 84 to calculate *ϕ*_{VC,se} and *ϕ*_{diff,se} separately.

We implemented the differential equations in Python 3.6 and solved them using the `solve_ivp` function from SciPy. We used its default Runge-Kutta method of order 5(4), and set the maximal allowed step size to 10^{−4}. The code is made available at https://github.com/CINPLA/EDPRmodel and https://github.com/CINPLA/EDPRmodel_analysis.

AL Hodgkin, AF Huxley. . A quantitative description of membrane current and its application to conduction and excitation in nerve.

RD Traub, RK Wong, R Miles, H Michelson. . A model of a CA3 hippocampal pyramidal neuron incorporating voltage-clamp data on intrinsic conductances.

PF Pinsky, J Rinzel. . Intrinsic and network rhythmogenesis in a reduced Traub model for CA3 neurons.

ZF Mainen, J Joerges, JR Huguenard, TJ Sejnowski. . A model of spike initiation in neocortical pyramidal neurons.

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.

G Halnes, S Augustinaite, P Heggelund, GT Einevoll, M Migliore. . A multi-compartment model for interneurons in the dorsal lateral geniculate nucleus.

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.

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.

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.

H Markram, E Muller, S Ramaswamy, MW Reimann, M Abdellah, CA Sanchez, et al . Reconstruction and simulation of neocortical microcircuitry.

A Arkhipov, NW Gouwens, YN Billeh, S Gratiy, R Iyer, Z Wei, et al . Visual physiology of the layer 4 cortical circuit in silico.

G Buzsáki, CA Anastassiou, C Koch. . The origin of extracellular fields and currents?EEG, ECoG, LFP and spikes.

GT Einevoll, C Kayser, NK Logothetis, S Panzeri. . Modelling and analysis of local field potentials for studying the function of cortical circuits.

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.

KH Pettersen, GT Einevoll. . Amplitude variability and extracellular low-pass filtering of neuronal spikes.

GG Somjen. . Mechanisms of Spreading Depression and Hypoxic Spreading Depression-Like Depolarization.

F Fröhlich, M Bazhenov. . Potassium dynamics in the epileptic cortex: new insights on an old topic.

BJ Zandt, B ten Haken, MJ van Putten, MA Dahlem. . How does spreading depression spread? Physiology and modeling.

C Ayata, M Lauritzen. . Spreading depression, spreading depolarizations, and the cerebral vasculature.

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.

H Kager, WJ Wadman, GG Somjen. . Simulated seizures and spreading depression in a neuron model incorporating interstitial space and ion concentrations.

L Øyehaug, I Østby, CM Lloyd, SW Omholt, GT Einevoll. . Dependence of spontaneous neuronal firing and depolarisation block on astroglial membrane transport mechanisms.

Y Mori. . A multidomain model for ionic electrodiffusion and osmosis with an application to cortical spreading depression.

G Halnes, T Mäki-Marttunen, D Keller, KH Pettersen, OA Andreassen, GT Einevoll. . Effect of ionic diffusion on extracellular potentials in neural tissue.

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.

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.

G Somjen, H Kager, W Wadman. . Computer simulations of neuron-glia interactions mediated by ion flux.

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.

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.

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.

J Lee, SJ Kim. . Spectrum measurement of fast optical signal of neural activity in brain tissue and its theoretical origin.

BJ Zandt, B ten Haken, JG van Dijk, MJ van Putten. . Neural dynamics during anoxia and the “wave of death”.

N Hübel, E Schöll, MA Dahlem. . Bistable dynamics underlying excitability of ion homeostasis in neuron models.

MA Dahlem, J Schumacher, N Hübel. . Linking a genetic defect in migraine to spreading depression in a computational model.

N Hübel, MA Dahlem. . Dynamics from seconds to hours in Hodgkin-Huxley model with time-dependent ion concentrations and buffer reservoirs.

N Hübel, G Ullah. . Anions govern cell volume: a case study of relative astrocytic and neuronal swelling in spreading depolarization.

N Hübel, MS Hosseini-Zare, J Žiburkus, G Ullah. . The role of glutamate in neuronal ion homeostasis: A case study of spreading depolarization.

H Kager, W Wadman, G Somjen. . Conditions for the triggering of spreading depression studied with computer simulations.

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.

MD Forrest, MJ Wall, DA Press, J Feng. . The sodium-potassium pump controls the intrinsic firing of the cerebellar Purkinje neuron.

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.

MD Forrest. . Simulation of alcohol action upon a detailed Purkinje neuron model and a simpler surrogate model that runs> 400 times faster.

GP Krishnan, G Filatov, A Shilnikov, M Bazhenov. . Electrogenic properties of the Na+/K+ ATPase control transitions between normal and pathological brain states.

A Zylbertal, A Kahan, Y Ben-Shaul, Y Yarom, S Wagner. . Prolonged intracellular Na+ dynamics govern electrical activity in accessory olfactory bulb mitral cells.

A Zylbertal, Y Yarom, S Wagner. . The Slow Dynamics of Intracellular Sodium Concentration Increase the Time Window of Neuronal Integration: A Simulation Study.

Y Mori, C Peskin. . A numerical method for cellular electrophysiology based on the electrodiffusion equations with internal boundary conditions at membranes.

G Halnes, I Østby, KH Pettersen, SW Omholt, GT Einevoll. . Electrodiffusive model for astrocytic and neuronal ion concentration dynamics.

S Niederer. . Regulation of ion gradients across myocardial ischemic border zones: a biophysical modelling analysis.

A Ellingsrud, A Solbrå, G Einevoll, G Halnes, M Rognes. . Finite element simulation of ionic electrodiffusion in cellular geometries.

G Yi, Y Fan, J Wang. . Metabolic cost of dendritic Ca2+ action potentials in layer 5 pyramidal neurons.

K Mizuseki, S Royer, K Diba, G Buzsáki. . Activity dynamics and behavioral correlates of CA3 and CA1 hippocampal pyramidal neurons.

C Reiffurth, M Alam, M Zahedi-Khorasani, S Major, JP Dreier. . Na+/K+-ATPase *α* isoform deficiency results in distinct spreading depolarization phenotypes.

CM Van Rijn, H Krijnen, S Menting-Hermeling, AM Coenen. . Decapitation in rats: latency to unconsciousness and the?wave of death?

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.

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.

W Rall. Core conductor theory and cable properties of neurons In: ER Kandel, JM Brookhardt, V M Mountcastle, editors.

R O’Connell, Y Mori. . Effects of Glia in a Triphasic Continuum Model of Cortical Spreading Depression.

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.

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.

Errata, Principles of computational modelling in neuroscience;. http://www.compneuroprinciples.org/errata.

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.

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.

Sincerely,

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 ploscompbiol@plos.org immediately:

[LINK]

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 Biology*data 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, https://pacev2.apexcovantage.com. 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 figures@plos.org.

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: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, PLOS recommends that you deposit laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see http://journals.plos.org/ploscompbiol/s/submission-guidelines#loc-materials-and-methods

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 Biology*data 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

PCOMPBIOL-D-20-00066R1

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 ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol

Citing articles via

Tweets

https://www.researchpad.co/tools/openurl?pubtype=article&doi=10.1371/journal.pcbi.1007661&title=An 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,

© 2020 Newgen KnowledgeWorks | Privacy & Cookie Policy | Powered by: Nova