article: Many pathological conditions, such as seizures, stroke, and spreading depression, are associated with substantial changes in ion concentrations in the extracellular space (ECS) of the brain. An understanding of the mechanisms that govern ECS concentration dynamics may be a prerequisite for understanding such pathologies. To estimate the transport of ions due to electrodiffusive effects, one must keep track of both the ion concentrations and the electric potential simultaneously in the relevant regions of the brain. Although this is currently unfeasible experimentally, it is in principle achievable with computational models based on biophysical principles and constraints. Previous computational models of extracellular ion-concentration dynamics have required extensive computing power, and therefore have been limited to either phenomena on very small spatiotemporal scales (micrometers and milliseconds), or simplified and idealized 1-dimensional (1-D) transport processes on a larger scale. Here, we present the 3-D Kirchhoff-Nernst-Planck (KNP) framework, tailored to explore electrodiffusive effects on large spatiotemporal scales. By assuming electroneutrality, the KNP-framework circumvents charge-relaxation processes on the spatiotemporal scales of nanometers and nanoseconds, and makes it feasible to run simulations on the spatiotemporal scales of millimeters and seconds on a standard desktop computer. In the present work, we use the 3-D KNP framework to simulate the dynamics of ion concentrations and the electrical potential surrounding a morphologically detailed pyramidal cell. In addition to elucidating the single neuron contribution to electrodiffusive effects in the ECS, the simulation demonstrates the efficiency of the 3-D KNP framework. We envision that future applications of the framework to more complex and biologically realistic systems will be useful in exploring pathological conditions associated with large concentration variations in the ECS.
  The brain mainly consists of a dense packing of neurons and neuroglia, submerged in the cerebrospinal fluid which fills the extracellular space (ECS). Neurons generate their electrical signals by exchanging ions with the ECS through ion-selective channels in their plasma membranes. During normal signaling, this does not lead to significant changes in local ion concentrations, as neuronal and glial transport mechanisms work towards maintaining ion concentrations close to baseline levels. However, endured periods of enhanced neuronal activity or aberrant ion transport may lead to changes in ECS ion concentrations. Local concentration changes often coincide with slow shifts in the ECS potential [1–3], which may be partly evoked by diffusive electrical currents, i. e. , currents carried by charged ions moving along ECS concentration gradients [2,4]. While concentration gradients can influence electrical fields, the reverse is also true, since ions move not only by diffusion but also by electric drift. A better understanding of the electrodiffusive interplay between ECS ion dynamics and ECS potentials may be a prerequisite for understanding the mechanisms behind many pathological conditions linked to substantial concentration shifts in the ECS, such as epilepsy and spreading depression [3,5–7]. A simultaneous and accurate knowledge of the concentration of all ion species is needed to make reliable estimates of electrodiffusive effects in the ECS. Although this is currently unfeasible experimentally, it is in principle achievable with computational models based on biophysical principles and constraints. However, in most computational models in neuroscience ion-concentration dynamics are only partially modeled, or are ignored altogether. One reason for this is the challenge involved in keeping track of all ion concentrations and their spatiotemporal dynamics. Another reason may be the strong focus within the community on modeling the neuronal membrane dynamics at short timescales, during which both intra- and extracellular concentration changes are relatively small and putatively negligible. Although there exist models that account for ion concentration shifts and their effects on neuronal and glial reversal potentials [8–11], the most common computational models for excitable cells, the multi-compartmental models and the cable equation, are based on the assumptions that (i) the ECS potential is constant (ground), and (ii) the ion concentrations are constant [12,13]. The NEURON simulator [14,15] is based on these assumptions, and although they are physically incorrect, they still allow for efficient and fairly accurate predictions of the membrane-potential dynamics. Because of assumption (i), multi-compartmental models are unsuited for modeling ECS dynamics, and several approaches have been taken to construct models which include ECS effects. A majority of computational studies of ECS potentials are based on volume conductor (VC) theory [16–21]. VC-schemes link neuronal membrane dynamics to its signatures in the ECS potential. In contrast to the multi-compartmental models, VC-schemes are derived by allowing the ECS potential to vary, but still assuming that the ion concentrations are constant. VC-schemes are attractive, because they offer closed-form solutions, and allow the calculation of the electric field for arbitrarily large systems. Although it may be reasonable to neglect variations in ECS ion concentrations on short timescales, the accumulative effects of endured neuronal activity may lead to significant concentration changes in the ECS, which are related to the aforementioned pathological conditions. Naturally, models that do not include ion-concentration dynamics are not applicable for exploring such pathologies. Furthermore, VC-schemes neglect the effects from diffusive currents on the ECS potentials [4,22,23], and in previous computational studies we have found the low-frequency components of the ECS potential to be dominated by diffusion effects [4,24]. A simplified approach to modeling concentration dynamics in brain tissue, is to use reaction-diffusion schemes (see e. g. , [25–27]). In these schemes, concentration dynamics are simulated under the simplifying assumption that ions move due to diffusion only. This approach has been used for many specific applications, giving results in close agreement with experimental data [26]. However, the net transport of abundant charge carriers such as Na+, K+, Ca2+, and Cl−, is also influenced by electric forces, which is not incorporated in diffusion only (DO) -schemes. Furthermore, DO-schemes do not include the influence that diffusing ions can have on the electrical potential. To account for the electric interactions between the different ion species, as well as the effect of such electric forces on the ECS potential, an electrodiffusive modeling framework is needed. The most detailed modeling scheme for electrodiffusion is the Poisson-Nernst-Planck (PNP) scheme [28–34]. The PNP-scheme explicitly models charge-relaxation processes, that is, tiny deviations from electroneutrality involving only about 10—9 of the total ionic concentration [35]. This requires a prohibitively high spatiotemporal resolution, which makes the PNP-scheme too computationally expensive for modeling the ECS on the tissue scale. Even the state-of-the-art simulations in the literature are on the order of milliseconds on computational domains of micrometers. The PNP-scheme is therefore not suited for simulating processes taking place at the tissue scale [23]. A series of modeling schemes have been developed that circumvent the brief charge-relaxation processes, and solve directly for the ECS potential when the system is in a quasi-steady state [4,23,36–42]. Circumventing charge-relaxation allows for simulations on spatiotemporal scales which are larger, compared to what is possible with the PNP-scheme, by several orders of magnitude. The charge-relaxation can be bypassed by replacing Poisson’s equation with the constraint that the bulk solution is electroneutral. These schemes have been shown to deviate from the PNP-scheme very close to the cell membrane (less than 5 Î¼m), but to give a good agreement in the bulk solution [23]. The simplest electroneutral modeling scheme is the Kirchhoff-Nernst-Planck (KNP) scheme, previously developed in our group [41,42]. A similar scheme was developed in parallel in the heart cell community [40]. The KNP-scheme has previously been used to study electrodiffusive phenomena such as spatial K+ buffering by astrocytes [41], effects of ECS diffusion on the local field potential [4], and the implication for current-source density analysis [24]. For simplicity, these previous applications were limited to idealized 1-D setups with a relatively coarse spatial resolution. Furthermore, a comparison between the KNP framework and other simulation frameworks was not included in previous studies. In the present study, we introduce a 3-D version of the KNP framework which can be used to simulate the electrodiffusive dynamics of ion-concentrations and the electrical potential in the ECS on large spatiotemporal scales. We establish in which situations the assumptions used in the KNP scheme are warranted by comparing it to the more physically detailed PNP scheme. Furthermore, we identify the conditions under which an electrodiffusive formalism is needed by comparing the KNP scheme to the VC and DO schemes. The simplified schemes can be derived from the KNP scheme by assuming, respectively, that (for VC) diffusive effects on the membrane potential and (for DO) migratory effects on the concentration dynamics are negligible. Accordingly, the accuracy of the simplifying assumptions can be assessed by comparing how close their predictions come to the KNP scheme. We present the results of three distinct simulation setups, which we will refer to as Application 1, Application 2, and Application 3 for the remainder of this study: In Application 1, we consider an idealized 1-D domain filled with a salt solution, starting with a nonzero ion concentration gradient. We solve the system using the PNP-scheme, the KNP-scheme, and a DO-scheme. We compare results on short and long timescales (nanoseconds and seconds), to highlight the similarities and differences between the schemes. In Application 2, we consider a 3-D domain with an ion concentration point source and a point sink, of equal magnitude, embedded in a standard ECS ion solution. We compare results obtained with the VC- and KNP-schemes to highlight their similarities and differences. In Application 3, we consider a morphologically realistic pyramidal neuron model [43]) embedded in a 3-D ECS solution. The neuronal morphology is inserted as a 1-D branching tree, which means that it does not occupy any volume, but gives rise to a morphologically realistic spatial distribution of neuronal membrane current sources or sinks. The ECS dynamics is computed using the KNP-scheme, and show how concentration gradients gradually build up in the ECS due to the neural activity, and how this influences the local potential in the ECS. We compare results obtained with the VC-, DO-, and KNP-schemes to highlight their similarities and differences. The first two applications are simplified simulation setups, used to better understand the differences between the schemes introduced above, while the third application is the main result of this study, as it illustrates the scales at which the KNP-scheme can be used. To our knowledge, the KNP-scheme is the first simulation framework which can handle 3-D electrodiffusion in neuronal tissue at relatively large spatiotemporal scales without demanding an insurmountable amount of computer power. For Application 3, the long-term ECS ion-concentration dynamics (about 100 s) in a spatial region of about 1 mm3 was run on a standard desktop computer within a day. We expect that the presented simulation framework will be of great use for future studies, especially for modeling tissue dynamics in the context of exploring pathological conditions associated with large shifts in ECS ion concentrations [3,5–7].
  The solver for the above modeling schemes was implemented utilizing FEniCS, an open-source platform for solving partial differential equations using the finite element method [51]. The reader is referred to [52] for an extensive introduction to the finite element method. The time derivative was approximated using an implicit Euler time-stepping scheme. We chose to use this method as the PNP equations are highly unstable, and the implicit Euler step offers superior stability to other methods [33]. Employing the implicit Euler step makes the PNP-scheme fairly stable, at the cost of numerical efficiency. Piecewise linear Lagrangian elements were used for all unknowns in all simulations [52]. The system was solved monolithically, using Newton’s method to solve for each time step. Due to limitations of FEniCS’ built-in nonlinear solver, we implemented our own Newton solver, which is found in the source code. Details on the FEniCS implementation are found in S1 Appendix. For readability, we include the complete set of equations used in the KNP-scheme. ∂ c k ∂ t - ∇ · [ D ˜ k ∇ c k + D ˜ k z k c k ψ ∇ ϕ ] = f ˜ k, in Ω, (27) ∇ · (σ ∇ ϕ VC) + F ∑ k ∈ ions f ˜ k z k + i ˜ cap = 0, in Ω, (28) ∇ · (σ ∇ ϕ diff + ∇ b) = 0, in Ω, (29) c k = c k, out, at ∂ Ω, (30) (σ ∇ ϕ diff + ∇ b) · n = 0, at ∂ Ω, (31) ∇ ϕ VC · n = 0, at ∂ Ω. (32) where we use the following definitions ϕ = ϕ VC + ϕ diff, (33) σ = F ∑ k ∈ ions D ˜ k z k 2 ψ c k, (34) b = F ∑ k ∈ ions D ˜ k z k c k. (35) Furthermore, f ˜ k and i ˜ cap are assumed to be known quantities. In the following, we go through the details of the three applications presented.
  The differences and similarities between KNP, PNP and DO were illustrated using the simplified 1-D domain representing a box filled with a solution containing two ion species, Na+ and X−, with an initial concentration gradient (see Methods). When we simulated Application 1 using the PNP-scheme, two distinct modes of behavior were revealed. Firstly, an initial charge-relaxation mode occurred. During this mode, the X− concentration spread towards the left side of the box faster than the Na+ concentration, since X− had a larger diffusion coefficient than Na+. Due to the charge separation caused by this process, an electrical potential difference rapidly built up in the system (see blue line in Fig 2A). The charge separation process only went on for a few nanoseconds, before the system reached a quasi-steady state. Secondly, a quasi-steady state occurred, in which the potential difference was such that diffusion and electrical drift were in equilibrium, and further charge separation was prevented. The quasi-steady state potential changed on a slow timescale of seconds (Fig 2B). This was the timescale at which the concentration differences in the system evened out. A comparison between the two lines in Fig 2A and 2B highlights the key difference between the PNP (blue line) and KNP (yellow dashed line) schemes. The KNP-scheme does not model the relaxation of the charge concentration, but derives the quasi-steady state potential directly. The results using KNP therefore deviated from those using PNP only on very short timescales during which the quasi-steady state potential was built up. After the initial charge-relaxation (t> 5 ns), the two schemes were in excellent agreement. In theoretical studies of binary systems of ions, it has been shown that the ion concentrations will develop approximately as a single particle species moving by simple diffusion with a modified diffusion coefficient [54], D ˜ NaX = 2 D ˜ Na D ˜ X D ˜ Na + D ˜ X. (42) In Fig 2C the ion concentration profiles obtained with this theoretical approximation have been compared to a KNP-based simulation, as well as a DO-simulation where the electric charge of the ions was ignored, so that they were allowed to diffuse independently. The KNP solution was in agreement with Eq 42, and shows that the anion and cation concentrations coincided closely due to the strong electrical forces that would be associated with charge separation. This comparison highlights the shortcomings of the DO-model, which (i) gave inaccurate predictions of the two ion-concentration profiles, and (ii) predicted that the twin ion species were not locally balanced, so that the system contained a charge density ρ = F (cNa−cX) which varied over the x-axis. According to the Poisson equation (Eq 8), the charge density associated with the concentration profiles in Fig 2C would amount voltage differences of 247 kV across the system, which is clearly not physically realistic. We next compared the KNP-scheme with the simpler VC-scheme for computing ECS potentials. We first did this for an idealized setup, using the computational domain of a 3-D box consisting of the extracellular domain. The initial ion concentrations were uniform, with values from Table 2. The box contained one K+ point source, and one K+ point sink, with equal magnitudes (see Methods). The electrical potential ϕ at the end of the stimulus period (t = 1 s) is shown in Fig 3A. The resulting electrical potential could be split into two contributions, i. e. , the contribution from the source and sink currents, ϕVC (Fig 3B), and the contribution from the diffusive currents in the ECS, ϕdiff (Fig 3C). In this simulation, ϕVC had a constant value which was reached immediately when the source and sink were turned on, and remained constant throughout the stimulus period since the source and sink were constant. This was not the case for ϕdiff, which was close to zero at t = 0, but increased gradually throughout the simulation, as concentration gradients built up in the system. The diffusion potential ϕdiff was quite local, and gave strongest contribution close to the source and sink (Fig 3C). The differences between the KNP scheme and VC scheme are illustrated in Fig 3D, which shows the potential difference between two points close to the source and sink. VC and KNP gave similar estimates of the potential ϕ only for a brief period after the source and sink had been turned on, when the contribution from diffusion was still small. As the simulation progressed, the local build-up of ion concentrations evoked a shift in the KNP-simulated potential ϕ, but left ϕVC unaffected. The difference between the schemes ϕ − ϕVC is identical to the diffusion potential ϕdiff, which is not accounted for by the VC scheme. The diffusion component had a sort of screening effect, reducing the potential difference between the source and sink compared to what we would predict in the absence of diffusion (i. e. |ϕ| < |ϕVC|). After roughly 0. 1 seconds, ϕdiff had lowered ϕ by about 5%, compared to what the VC model predicted. This effect lasted until the source and sink were turned off (Fig 3D). After this, ϕVC immediately dropped to zero, while in the KNP-scheme, a diffusion potential ϕdiff remained. In the absence of sources, ϕdiff decayed gradually until system reached a steady state. Finally, we used the KNP-scheme to model electrodiffusion in a 3-D cylinder containing a morphologically detailed neuron (Fig 1A). The neuron was represented as a set of point sources and sinks that had a spatial distribution representing the neuronal morphology, and which exchanged ions with the system as illustrated in Fig 1B (see Methods for details). Since ECS ion concentrations generally vary on a slow timescale (order of seconds), the system was simulated for 80 s. Due to the efficiency of the KNP-scheme, we could run this simulation in less than 15 hours using 4 CPU cores in parallel on a standard desktop computer.
  The presented implementation of the 3-D KNP-scheme was limited to a relatively small piece of neural tissue, which included only a single pyramidal neuron modeled with the NEURON simulator. Although this approach provided clarity regarding the possible single-neuron contribution to the ECS electro- and concentration dynamics, there are many reasons why this simulation setup was far from representative for any biologically realistic scenario. Firstly, the brain is packed with neurons and glial cells, and even the small volume considered here would in reality be populated by thousands of neurons that would contribute to the ECS dynamics. In a realistic scenario, ECS concentration shifts may become much larger than for the single-neuron case, where the maximal changes in ion concentration were roughly 1% (for K+ near the soma). Accordingly, the effects of electrodiffusion in the single neuron-setup used here were small, and the simulations in Application 3 should predominantly be regarded as a proof of principle, and as a precursor for what could be expected in more realistic scenarios with additive contributions from networks of neurons to extracellular ion dynamics. Secondly, the simulations were performed under the simplifying assumption that there was no feedback from the ECS dynamics to the neurodynamics. This assumption gave us the advantage that the neuronal output did not depend on the modelling scheme used for the ECS dynamics, meaning that VC-, DO- and KNP-schemes could be compared under identical input/output conditions. In reality, however, shifts in extracellular (or intracellular) concentrations would influence the reversal potential of transmembrane currents, which in turn would have an impact on the neurodynamics (see [8–11]), and these effects are not accounted for in the current model. In addition, simulations based on the NEURON simulator do not account for ephaptic feedback from extracellular fields onto neurons, which could also have moderate impact on neurodynamics (see [55]). Thirdly, the neuron model used in our study did not include the Na+/K+-exchanger pump [43]. This is the case for most neuron models currently available in NEURON (for a model with ion pumps, see [56]). Together with astrocytic uptake mechanisms, the exchanger pump would strive towards maintaining extra- and intracellular ion concentrations close to baseline levels. For a biophysical modeling scheme that derives the transmembrane transport through channels and pumps from first principles, see e. g. , [57,58]. In the presence of such mechanisms, the single-neuron contribution to ECS concentration shifts would likely be smaller than in the simulations presented here, or would require a higher neural activity level in order to occur. The model limitations mentioned above were also present in the previous 1-D implementation of the KNP-scheme, and we refer to this previous work for a more thorough discussion [4]. Several previous studies have explored ECS electrodiffusion on small spatiotemporal scales [23,28–34,36–39]. Electrodiffusive models tailored to explore larger spatiotemporal scales have to our knowledge so far only been implemented in 1-D [4,24,40–42]. We have previously used a 1-D implementation of the KNP-scheme to explore the effect of diffusive currents on ECS potentials [4]. Qualitatively, the results from the 1-D simulations were similar to those found in the current 3-D model, and in both cases it was concluded that diffusive currents affected the low-frequency components of ECS potentials (Fig 7). In the 1-D implementation, diffusion and electrical drift were confined to the spatial direction along the neuronal extension, i. e. along the basal dendrite-soma-apical dendrite axis, and did not occur in the lateral direction. As discussed in [4], the 1-D assumption is equivalent to assuming lateral homogeneity, which means that any neuron is surrounded by neighboring neurons with identical activity patterns as itself. Clearly, this is would not apply to most biological scenarios, which means that the results obtained with the 3-D implementation are generally more reliable. In contrast to the 1-D model, the 3-D model predicts how ion concentrations vary in all spatial directions surrounding the neuronal sources, and can thus be applied to a broader range of biological problems. The diffusion potentials addressed by the KNP-scheme are those that arise due to ECS concentration gradients on relatively large spatial scales. Diffusion potentials of this kind are often referred to as liquid junction potentials, since they are most pronounced at the boundary between two solutions of different ion compositions [35,59]. These diffusion potentials are unrelated to filtering effects hypothesized to arise due to diffusion in the vicinity of the membrane when electric charge is transferred from the intracellular to the extracellular space [60,61]. The presented version of the KNP-scheme was developed for use in a hybrid simulation setup where the dynamics of ion concentrations and the electrical potential in the ECS were computed with KNP, while the neurodynamics was computed with the NEURON simulator tool. By necessity, this scheme shares the limitations of the NEURON simulator in terms of handling intracellular ion dynamics, which by default is not electrodiffusive in the NEURON environment [14]. A natural future endeavor will therefore be to derive a computational scheme that in a consistent way couples both the intra- and extracellular ion-concentration and voltage dynamics based on the KNP-scheme. Such a scheme will represent a generalization of the previously developed extracellular-membrane-intracellular (EMI) model [55,62], which in a consistent way couples the intra- and extracellular electrodynamics, but so far does not include ion-concentration dynamics and thus not diffusive currents. This being said, the hybrid KNP/NEURON version presented here is valuable in its own right, since it allows the KNP framework to be combined with the many models that are already available in the NEURON software. Although we here only studied a single neuron’s contribution to the ECS electro- and concentration dynamics, we envision future applications of the 3-D KNP to more complex systems. The 3-D version of the KNP-scheme is generally applicable, and could be used to simulate the dynamics of ion concentrations and the electrical potential in the ECS surrounding arbitrary complex models of neurons or populations of neurons. It could, for example, be used to account for extracellular transport processes in the Blue Brain Simulator [63]. If combined with a cellular system of neurons and glial cells tailored to represent a specific experimental condition, the KNP-scheme could be used to explore the mechanisms behind spreading depression [3,7, 64], epilepsy [5,65], or other pathological conditions associated with large shifts in extracellular ion concentrations. In such applications, however, the hybrid KNP framework would need to be expanded to also account for the effect of ECS concentration dynamics on neuronal reversal potentials. This would require the development of an application programming interface that allows two-way communication between the simulator tools NEURON (for neurons) and FEniCS (for ECS). Since neurodynamics typically occurs at a much faster time scale (millisceonds) than changes in ECS ion concentrations (seconds), the interface could be based on using multiple time scales in order to keep the computational cost manageable. For example, a millisecond resolution could be used in NEURON simulations, while ECS concentrations could be computed (and updated in NEURON) at a slower temporal resolution of e. g. , seconds.
summary: Many pathological conditions, such as epilepsy and cortical spreading depression, are linked to abnormal extracellular ion concentrations in the brain. Understanding the underlying principles of such conditions may prove important in developing treatments for these illnesses, which incur societal costs of tens of billions annually. In order to investigate the role of ion-concentration dynamics in the pathological conditions, one must measure the spatial distribution of all ion concentrations over time. This remains challenging experimentally, which makes computational modeling an attractive tool. We have previously introduced the Kirchhoff-Nernst-Planck framework, an efficient framework for modeling electrodiffusion. In this study, we introduce a 3-dimensional version of this framework and use it to model the electrodiffusion of ions surrounding a morphologically detailed neuron. The simulation covered a 1 mm3 cylinder of tissue for over a minute and was performed in less than a day on a standard desktop computer, demonstrating the framework’s efficiency. We believe this to be an important step on the way to understanding phenomena involving ion concentration shifts at the tissue level.
keywords: medicine and health sciences
  membrane potential
  ions
  electrophysiology
  computerized simulations
  neuroscience
  computational neuroscience
  neuronal dendrites
  physical chemistry
  computer and information sciences
  extracellular space
  animal cells
  chemistry
  cellular neuroscience
  neuronal morphology
  cell biology
  physiology
  neurons
  single neuron function
  biology and life sciences
  cellular types
  physical sciences
  computational biology
title: A Kirchhoff-Nernst-Planck framework for modeling large scale extracellular electrodiffusion surrounding morphologically detailed neurons