RESEARCH ARTICLE
A Kirchhoff-Nernst-Planck framework for
modeling large scale extracellular
electrodiffusion surrounding morphologically
detailed neurons
Andreas Solbrå
ID
1,2
, Aslak Wigdahl Bergersen
3
, Jonas van den Brink
3
,
Anders Malthe-Sørenssen
ID
1,2
, Gaute T. Einevoll
ID
1,2,4
, Geir Halnes
ID
4
*
1 Center for Integrative Neuroplasticity, University of Oslo, Oslo, Norway, 2 Department of Physics, University
of Oslo, Oslo, Norway, 3 Simula Research Laboratory, Fornebu, Norway, 4 Department of Mathematical
Sciences and Technology, Norwegian University of Life Sciences, Ås, Norway
* geir.halnes@nmbu.no
Abstract
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 cur-
rently unfeasible experimentally, it is in principle achievable with computational models
based on biophysical principles and constraints. Previous computational models of extracel-
lular 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 electroneu-
trality, 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 concentra-
tions 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.
PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1006510 October 4, 2018 1 / 26
a1111111111
a1111111111
a1111111111
a1111111111
a1111111111
OPEN ACCESS
Citation: Solbrå A, Bergersen AW, van den Brink J,
Malthe-Sørenssen A, Einevoll GT, Halnes G (2018)
A Kirchhoff-Nernst-Planck framework for modeling
large scale extracellular electrodiffusion
surrounding morphologically detailed neurons.
PLoS Comput Biol 14(10): e1006510. https://doi.
org/10.1371/journal.pcbi.1006510
Editor: Ernest Barreto, George Mason University,
UNITED STATES
Received: February 6, 2018
Accepted: September 12, 2018
Published: October 4, 2018
Copyright: © 2018 Solbra
˚
et al. This is an open
access article distributed under the terms of the
Creative Commons Attribution License, which
permits unrestricted use, distribution, and
reproduction in any medium, provided the original
author and source are credited.
Data Availability Statement: All model code will be
available at https://github.com/CINPLA/KNPsim.
Funding: This work was funded by the Research
Council of Norway (BIOTEK2021 Digital Life project
‘DigiBrain’, project 248828). The funders had no
role in study design, data collection and analysis,
decision to publish, or preparation of the
manuscript.
Competing interests: The authors have declared
that no competing interests exist.
Author 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 inves-
tigate the role of ion-concentration dynamics in the pathological conditions, one must
measure the spatial distribution of all ion concentrations over time. This remains chal-
lenging 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 mm
3
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 understand-
ing phenomena involving ion concentration shifts at the tissue level.
Introduction
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 activ-
ity or aberrant ion transport may lead to changes in ECS ion concentrations. Local concentra-
tion 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 depres-
sion [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 unfea-
sible experimentally, it is in principle achievable with computational models based on biophys-
ical 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 spatiotempo-
ral 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 mod-
els 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-com-
partmental 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
Modelling large scale electrodiffusion surrounding morphologically detailed neurons
PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1006510 October 4, 2018 2 / 26
incorrect, they still allow for efficient and fairly accurate predictions of the membrane-poten-
tial 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 varia-
tions in ECS ion concentrations on short timescales, the accumulative effects of endured neu-
ronal 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-con-
centration 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 reac-
tion-diffusion schemes (see e.g., [25–27]). In these schemes, concentration dynamics are simu-
lated 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 experi-
mental data [26]. However, the net transport of abundant charge carriers such as Na
+
, K
+
,
Ca
2+
, 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 tak-
ing place at the tissue scale [23].
A series of modeling schemes have been developed that circumvent the brief charge-relaxa-
tion 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 I
ˆ
¼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.
Modelling large scale electrodiffusion surrounding morphologically detailed neurons
PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1006510 October 4, 2018 3 / 26
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 com-
paring how close their predictions come to the KNP scheme.
We present the results of three distinct simulation setups, which we will refer to as Applica-
tion 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 (nanosec-
onds 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 real-
istic 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 similari-
ties 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 mm
3
was run on
a standard desktop computer within a day. We expect that the presented simulation frame-
work will be of great use for future studies, especially for modeling tissue dynamics in the con-
text of exploring pathological conditions associated with large shifts in ECS ion concentrations
[3, 5–7].
Materials and methods
This section is thematically split into three parts. We begin by explaining the necessary physi-
cal theory, stating and deriving the equations which we implemented. Then, we explain in
more detail how the models were implemented, including details such as numerical schemes
and boundary conditions. Finally, we give the specific details on each of the three applications
used in the study. The source code can be found online, at https://github.com/CINPLA/
KNPsim, and the results in this study can be reproduced by checking out the tag PLoS.
Modelling large scale electrodiffusion surrounding morphologically detailed neurons
PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1006510 October 4, 2018 4 / 26
Theory
The Nernst-Planck equation for electrodiffusion. The ion concentration dynamics of an
ion species in a solution is described by the continuity equation:
@c
k
@t
¼ r J
k
þ f
k
; in O;
ð1Þ
where c
k
is the concentration of ion species k, f
k
represent any source terms in the system, O is
the domain for which the concentrations are defined, and J
k
is the concentration flux of ion
species k. In the applications in this study, f
k
is implemented as a set of point sources at speci-
fied coordinates in the ECS. In the Nernst-Planck equation, J
k
consists of a diffusive and an
electric component:
J
k
¼ J
diff
k
þ J
field
k
: ð2Þ
The diffusive component is given by Fick’s first law,
J
diff
k
¼ D
k
rc
k
; ð3Þ
where D
k
is the diffusion coefficient of ion species k. The electric component is
J
field
k
¼
D
k
z
k
c
k
c
r;
ð4Þ
where ϕ is the electric potential, z
k
is the valency of ion species k, and ψ = RT/F is defined by
the gas constant (R), Faraday’s constant (F) and the temperature (T) which we assume to be
constant (cf. Table 1). Inserting Eqs 2–4 into Eq 1, yields the time evolution of the concentra-
tion of ion species k:
@c
k
@t
¼ r D
k
rc
k
þ
D
k
z
k
c
k
c
r
þ f
k
; in O: ð5Þ
We model the ECS as a continuous medium, while in reality, the ECS only takes up roughly
20% of the tissue volume [44] in the brain. To compensate for this, we use the porous medium
approximation [45]. This involves two changes to the model. The diffusion constants of the
ion species are modified as
~
D
k
¼
D
k
l
2
;
ð6Þ
where λ is the tortuosity, which accounts for various hindrances to free diffusion and electrical
migration through the ECS. We used the value λ = 1.6 [46]. We denote the fraction of tissue
volume belonging to the ECS by α, and set the value α = 0.2. The sources in the system are
Table 1. The physical parameters used in the simulations.
symbol explanation value
R gas constant 8.314 J/(K mol)
T temperature 300 K
F Faraday’s constant 9.648 × 10
4
C/mol
0
vacuum permittivity 8.854 × 10
−12
F/m
r
relative permittivity 80
https://doi.org/10.1371/journal.pcbi.1006510.t001
Modelling large scale electrodiffusion surrounding morphologically detailed neurons
PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1006510 October 4, 2018 5 / 26
![Table 2. Diffusion constants and baseline ECS concentrations for the ion species considered, with values as in [4]. All ion constants were modified as ~Dk ¼ Dk=l 2 , where λ = 1.6 is the tortuosity. The general anion X− was given the properties of Cl−.](/figures/table-2-diffusion-constants-and-baseline-ecs-concentrations-2je1r9ai.webp)
