242

Biophysical Journal

Volume 105

July 2013

242–254

Electrodiffusion Models of Neurons and Extracellular Space Using the
Poisson-Nernst-Planck Equations—Numerical Simulation of the Intraand Extracellular Potential for an Axon Model
Jurgis Pods,†‡* Johannes Schönke,†‡ and Peter Bastian†‡†Interdisciplinary Center for Scientific Computing, University of Heidelberg, Heidelberg, Germany; and ‡Bernstein Center for Computational
Neuroscience Heidelberg-Mannheim, Mannheim, Germany

ABSTRACT In neurophysiology, extracellular signals—as measured by local field potentials (LFP) or electroencephalography—are of great significance. Their exact biophysical basis is, however, still not fully understood. We present a threedimensional model exploiting the cylinder symmetry of a single axon in extracellular fluid based on the Poisson-Nernst-Planck
equations of electrodiffusion. The propagation of an action potential along the axonal membrane is investigated by means of
numerical simulations. Special attention is paid to the Debye layer, the region with strong concentration gradients close to
the membrane, which is explicitly resolved by the computational mesh. We focus on the evolution of the extracellular electric
potential. A characteristic up-down-up LFP waveform in the far-field is found. Close to the membrane, the potential shows amore intricate shape. A comparison with the widely used line source approximation reveals similarities and demonstrates thestrong influence of membrane currents. However, the electrodiffusion model shows another signal component stemming directlyfrom the intracellular electric field, called the action potential echo. Depending on the neuronal configuration, this might have asignificant effect on the LFP. In these situations, electrodiffusion models should be used for quantitative comparisons withexperimental data.

INTRODUCTION
Simulations of neuronal signal propagation have a long
history and are well established in neuroscience as a useful
tool to study brain function. Most models today are based on
the seminal work of Hodgkin and Huxley (1) for the membrane currents and the application of cable theory (2) toaccount for the neuron morphology. Among the mostwell-known simulators for these kinds of models are
NEURON (3) and GENESIS (4).
Although experimental setups with intracellular recordings can readily be replicated by compartment models,
extracellular measurements—an important tool in neurophysiology to study multiunit and network activity—cannot
be included directly.
There are models for obtaining the extracellular field of aneuron that are based on the membrane current source
density (5,6), based on the line source approximation(LSA) (7), which account for effects like frequency filteringof a complex extracellular space (8,9). In Pettersen et al.
(10), an inverse method was applied to estimate current
source densities from extracellular potentials. In thesemodels, the relevant parameters for the extracellular
medium are conductivity and permittivity. The problem isreduced to the solution of the electrostatic part of Maxwell’sequations, where the membrane is the only current source.
The local changes in ion concentrations caused by drift
and diffusion and their contribution to the electric field are
not considered.
Submitted December 10, 2012, and accepted for publication May 17, 2013.
*Correspondence: jurgis.pods@iwr.uni-heidelberg.de
Editor: Richard Bertram.
 2013 by the Biophysical Society
0006-3495/13/07/0242/13 $2.00

In Lopreore et al. (11), a detailed three-dimensional
numerical simulation of electrodiffusion has been carriedout for the node of Ranvier, showing the accumulationand depletion of ions close to the membrane, and therefore
the invalidity of the electroneutrality approximation close to
the membrane, as used in cable equation models. The study
focused on deviations from the cable equation, not on theextracellular signal. However, the membrane thickness
was greatly overestimated in this study as a consequenceof the coarseness of the spatial discretization, presumably
dictated by the available computational resources.
In Mori (12), the reason for the high computational
demand of electrodiffusion models based on the PoissonNernst-Planck (PNP) equations is discussed: The presenceof a Debye layer close to the membrane—over which concentrations change significantly—necessitates a fine spatial
resolution. A clever approximation is suggested that represents the Debye layer as a charge density boundary condition. In Mori et al. (13), this methodology was applied to
study the effect of gap junction conductances on cardiac
action potential (AP) propagation.
A recent study (14) investigated the AP propagation in areconstructed three-dimensional axonal structure intracellularly. Numerical methods for electrodiffusion-reaction
equations were analyzed in a comprehensive way in Lu
et al. (15), with special regard to surface potentials of
biomolecules.
To our knowledge, a model does not yet exist that explicitly resolves the Debye layer to study membrane
dynamics following neuronal excitation on the detailed levelof electrodiffusion.
http://dx.doi.org/10.1016/j.bpj.2013.05.041

Electrodiffusion Models of Axon and Extracellular Space243

In this work, we model the detailed evolution of the concentrations of the most relevant ion species and the resulting
electric field inside and (particularly) outside the cell during
the spread of an AP along an axonal membrane. To this end,
the Poisson-Nernst-Planck equations are solved numericallyby application of the finite element method. We propose an
efficient numerical scheme to tackle the computational
demand by using a suitable spatial grid that resolves themultiple spatial scales accurately while still using only aminimal number of unknowns.
We strive to elucidate how the complex interaction of ion
movements explains the evolution of the extracellularsignal. Finally, we compare our results with the LSA modeland explain the observed differences in signal shape.
THEORY AND METHODSThe Poisson-Nernst-Planck equations ofelectrodiffusion
The PNP system describes the movement of ions in a static solvent due todiffusion and electrostatic drift. It consists of the Nernst-Planck equationvni
þ V$Fi ¼ 0;
vt

(1a)

Fi ¼ Di ðVni þ zi ni VfÞ ;(1b)with the ion flux

and the Poisson equatione 2 n X
V$ðeVfÞ ¼ 
z i ni :e0 kT i(2)

The ni, i ¼ 1,., N, are relative concentrations (with respect to a scaling
concentration n* ¼ NA, equal to the Avogadro constant) with units inmillimolar for the N different ion species; zi is the valence; Di is the

position-dependent diffusion coefficient of the ion species i; f is thedimensionless relative electric potential energy with respect to the thermalenergy (f ¼ eU/kT) with U in volts; for room temperature, f ¼ 1 corresponds to ~25 mV; ε is the relative permittivity (which again may beposition-dependent); and T is the temperature of the solvent.
Equation 1 describes the time-dependent change in concentrations due todiffusion and drift through an electrical field. Equation 2 gives the electric
potential f at any point in space. We now use the PNP system to set up a
detailed model of an axon in an extracellular medium.

A three-dimensional model of axonal membraneand extracellular space
Our simulation model focuses on the axonal part of the neuron. Here, we areespecially interested in the membrane dynamics and its effects on the extracellular signal. As the solution of a full three-dimensional system is computationally expensive, we use a cylindrical coordinate system. Exploiting therotational symmetry of an idealized axon, the computational domain can be
reduced to two dimensions, as there is no change in angular direction. This
enables us to calculate valid three-dimensional results with a drasticallyreduced computational complexity. For the numerical solution, therectangular elements will be treated as (hollow) cylinders and the volumes
calculated accordingly. In the following, the two-dimensional geometry
will be used for illustration purposes, while implicitly meaning the cylindrical geometry.
In the upper part of Fig. 1, the cylindrical geometry is shown with thetwo-dimensional subset highlighted, constituting the effective computational domain. The x axis represents the domain’s symmetry axis, eliminating the angular coordinate q from the equations. The domain consistsof three partitions: Cytosol, membrane, and extracellular space. Cytosoland extracellular space are electrolytes, yielding the electrolyte domain
Uelec. It may contain an arbitrary number N of concentrations of different
ion species. Here, however, we will restrict ourselves to the minimal setof N ¼ 3 monovalent species sodium (Naþ), potassium (Kþ), and chloride(Cl). Sodium and potassium are needed for the ion channel dynamics
triggering an action potential; chloride is a representative of the anions
needed for electroneutrality in the bulk solution and does not cross themembrane in this model.
The membrane domain Umemb separates the two parts of the electrolyte
domain; therefore, Uelec is not connected. This setup necessitates the introduction of additional boundary conditions on the membrane-electrolyteinterface Gint next to obligatory boundary conditions on the domain

FIGURE 1 The two-dimensional computational
domain for the cylinder symmetric axon model.
The upper part shows the cylinder geometry into
which the two-dimensional computational grid isembedded, assuming symmetry in the angular
direction. (Lower part, solid lines) Domain boundary Gex. (Lower part, dashed lines) Interior (electrolyte-membrane) boundary Gint. This divides
the domain into the (unconnected) electrolyte
ðCYÞ
ðESÞ
domain Uelec ¼ Uelec WUelec and the separating
membrane Umemb domains. The Debye layer ofðESÞ
Uelec close to the membrane is highlighted(shading), followed by the near-field and far-fieldparts. The lower boundary represents the innercell symmetry axis. Note that this scheme is not
to scale, as the actual mesh sizes in y direction
differ by several orders of magnitude between the
Debye layer and far-field making the grid veryanisotropic.

Biophysical Journal 105(1) 242–254

244

Pods et al.

boundary Gext, such that the set of all boundary points is given by G ¼
Gint X Gext. The measures Uelec and Umemb form a partition of the computational domain: U ¼ Uelec WUmemb .
With this multidomain setup, we now have the Nernst-Planck expression
from Eq. 1 defined on Uelec and boundary conditions as

Fi , n ¼ jNP

on Gext;N WGint ;(3a)on Gext;D ;(3b)ni ¼ gNP

Membrane flux
The membrane flux boundary conditions are taken from an effective modelof ion channel currents. Given the time- and potential-dependent channel
conductances gi(f,t) from the Hodgkin-Huxley scheme (1)—here we use
the version from Koch (16), see the Supporting Material for details—thetotal conductance for each ion species is given bygK ¼ gKv þ gKL ;(5a)

gNa ¼ gNav þ gNaL ;

(5b)

(4a)

gCl ¼ 0;

(5c)

(4b)

with voltage-gated conductances gKv , gNav and fixed leak conductances gKL
and gNaL .
For the definition of the membrane flux, each point x ˛ GCYint on thecytosol-membrane interface is mapped to a point m(x) ˛ GESint on the
opposite membrane-extracellular space interface by a map m(x) ¼ x þdmemb , n, where dmemb is the membrane thickness and n is in this casethe unit outer normal at x pointing in the direction from cytosol tomembrane. The values of potential and concentrations evaluated at thesepoints are called fCY ¼ f(x), niCY ¼ ni(x), fES ¼ f(m(x)), and niES ¼
ni(m(x)). With these values on both sides of the membrane, the concentration-, potential-, and time-dependent membrane flux fimemb of speciesi ˛ {Na, K, Cl}is defined as

and the Poisson expression from Eq. 2 defined on the whole domain U withboundary conditions as

eVf , n ¼ jP
f ¼ gP

on Gext;N ;on Gext;D ;where n denotes the unit outer normal and G,,D and G,,N denote
Dirichlet and Neumann boundaries, respectively. For consistency, of course
Gext,N W Gext,D ¼ Gext holds.
Equations 1 and 2 are obviously defined on different domains. The
Poisson equation is defined on the whole domain, whereas the NernstPlanck equation is only defined on the electrolyte subdomain. This means
we assume the membrane to be free of mobile charge carriers for simplicity.
Membrane surface charges could easily be added as an additional source
term in the Poisson equation, but are not included here.
As can be seen in Fig. 1, the exterior boundary Gext ¼
right
top
WGleftGbottomext
ext WGext WGext consists of four parts, and the interior boundary
Gint of two nonconnected parts.
For the boundary conditions of Eqs. 3a, 3b, 4a, and 4b, we use membðni ; f; tÞ on ¼ Gint
fi
jNP ¼
;
0
on Gext;N ¼ Gext yGtop
ext
gNP ¼ n0i
jP ¼ 0on Gext;D ¼ Gtop
ext ;on Gext;N ¼ Gext yGtop
ext ;gP ¼ 0 on Gext;D ¼ Gtop
ext ;where we call the constant ni0 the bulk concentration of species i. The
membrane flux fimemb is defined in the next section. The Dirichletboundary conditions for the concentrations on the upper exterior boundary
serve to model an infinite reservoir for each ion species. The potential isclamped to zero at the upper extracellular boundary, which introduces anb max Þ calculated forerror equal to the value of the real potential value fðyan unrestricted domain (potential is 0 at infinity). For a point charge on
the membrane, this error would correspond to an absolute shift at eachpoint in the domain. Because the potential of a point charge falls off as
1/r in three dimensions, increasing the domain size by a factor 100
will reduce the error at a fixed point by a factor 1/100. For a finite line
charge (as in our case of an active membrane), this is only true if theradial distance is large compared to the length of the line charge (compare
to LSA, e.g., in Gold et al. (5)). In any case, increasing the domain size inthe y direction will reduce the error introduced by the upper Dirichletboundary. We chose a sufficiently large domain size of ymax ¼ 10 mm
to account for this.
At the lower boundary representing the intracellular symmetry axis, thepotential gradient and ion fluxes are vanishing such that no boundary
artifacts are introduced.
Biophysical Journal 105(1) 242–254



kT
nES
i
;
fimemb ðxÞ ¼ fimemb ðmðxÞÞ ¼ gi ðf; tÞ 2 2  z½f þ ln CY
ezn
ni
(6)

where [f] denotes the potential jump across the membrane, i.e., [f] ¼
fCY  fES. Note that two opposite points on the membrane interface are
identified with each other here, i.e., the membrane thickness is essentiallyneglected.
The flux was derived from the Hodgkin-Huxley membrane current
Iimemb ¼ gi(f,t)([U] – E) by replacing the constant battery E by a variable
concentration-dependent reversal potential calculated from the Nernst
equation (see below) and adding necessary scaling factors to bring fimemb ¼
Iimemb/ezn* to SI units mol/m2 s. This interior boundary condition fits nicely
into our framework, as it unifies the potential-dependent Hodgkin-Huxley
system with the concentration-dependent Nernst equation to arrive at an
expression that represents all the features of the potential- and concentration-dependent ion flux Fi of the PNP system.

Equilibrium states
The system’s steady state is equivalent to the neuronal resting state. The
equilibrium membrane potential can be easily validated: For one ion
species, the Nernst equation

E ¼

RT nESln
zF nCY(7)

gives the potential jump across the membrane.
The Goldman equation usually employed to calculate a neuron’s resting
potential is based on the assumption that the two electrolytes separated by
the membrane are bulk solutions, i.e., the ion concentrations are constantand only change within the membrane. As our model is based on the opposite assumptions—there are no ions within the membrane, but nonzero ion
fluxes in the neighboring electrolytes—the Goldman equation cannot be
used to validate the equilibrium state of our model.

Electrodiffusion Models of Axon and Extracellular Space

245

The model is initialized by setting the ion concentrations within
one electrolyte domain uniformly to their intra- and extracellular bulk
values ni0 (see Table S1 in the Supporting Material). The bulk
concentrations are chosen such that each electrolyte initially is electroneutral, i.e., the net charge is zero, which is a reasonable assumption
both physically and biologically with respect to energy minimization
principles.
To generate the resting state, the membrane leak channels are opened.
The active channels stay closed during the equilibration phase. Aftersome time (in the range of tens of milliseconds), an equilibrium state will
establish with a net zero membrane flux and the concentrations following
a Boltzmann distribution.
The Poisson-Boltzmann equation

V$ðeVfÞ ¼ 

e 2 n X 0
zi ni expðzi fÞ
e0 kT i(8)

predicts the concentration profile at equilibrium state. This equation will be
used for validation of the simulation later on, it is not explicitly used in thealgorithm.

Numerical methods
Finite element discretization
We use a finite element approach using the continuous Galerkin method tosolve the equations stated in The Poisson-Nernst-Planck Equations ofElectrodiffusion. Multiplying by test functions v, integrating over thedomain and applying integration by parts, we obtain the weak formulationof the equations in residual formulation.
For the Nernst-Planck equation, the temporal partZ

RNP;T ¼

ni v dt

i ¼ 1; .; N

(9a)

Uelec

Adaptive time-stepping

and the spatial partZ

RNP;S ¼

Z
Di ðVni þ zi ni VfÞ$Vv dx þUelec

jNP v dx i
GN

¼ 1; .; N(9b)
are combined, yielding

RNP ¼ RNP;T þ RNP;S :(9c)

For the Poisson equation, we get

Z
RP ¼

eVf$v þUimplicit Euler time-stepping scheme. As suggested in Fig. 1, the grid isrefined toward the membrane in the y direction. This is essential to resolve
the Debye length, the characteristic length scale over which the electrolyte
ion concentrations deviate significantly from their bulk values close to themembrane. In the x direction, the grid is equidistant (dx ¼ 100 mm); in the y
direction, it ranges from dymin ¼ 0.5 nm over the Debye layer to dymax ¼
100 mm, resulting in a very anisotropic grid.
It can easily be seen that each equation alone is linear in its unknowns.
One could therefore use an operator-split approach and solve the equations
alternately until convergence in each time step. We observed that a verysmall time step (on the order of nanoseconds) is necessary to solve thesystem this way.
Because the nonlinearity of the whole system results from the coupling ofboth equations, it seems reasonable to represent this crucial feature in the
numerical method. The system is therefore solved fully coupled using
Newton’s method, requiring the solution of one single linear system ineach iteration. For small systems, SuperLU (17) was used; for larger
systems with >50,000 unknowns, a stabilized BiConjugate Gradient
iterative solver, preconditioned by an inexact LU decomposition, turned
out to be faster while maintaining the same accuracy. We observe a significantly higher stability leading to possible time steps of approximately tensof microseconds for the fully coupled approach.
The cylinder symmetry introduces another subtle difficulty for thenumerical treatment of the system. Because the cell volumes increase in
the positive y direction (roughly with y dy), the entries of the full residual
R differ by several orders of magnitude (109 for the chosen domain sizeof 10 mm) solely by the presence of volume integrals in the weak formof the equations. We apply a threshold volume scaling to account for
this: At a certain distance from the membrane, a reference volume Vref iscalculated. All residuals belonging to an unknown at node i having
a volume Vi > Vref, where Vi is the minimum volume of all adjacent cellsof node i, are scaled by a factor Vref/Vi. This mathematically corresponds
to multiplying a diagonal matrix from left to the linear system matrix,
meaning that the same linear system is solved in each Newton iteration.
This scaling greatly increases the convergence properties of the Newton
algorithm for this cylinder geometry setup.

!
Z
e 2 n X
zi ni v dx þ
jP v dx; (10)e0 kT iGN

As the system has a great variability during the course of an action
potential, a small time step is needed to capture the large fluxes andpotential differences during this period. On the other hand, potential
differences during interspike intervals are small and so are themagnitudes of ion fluxes, allowing for the use of a larger time step.
Therefore, an adaptive time-stepping strategy is used to speed up thesimulation by controlling the time-step dt depending on the dynamics
of the system.
The time step is bounded by a minimum time step of tmin ¼ 0.05 ms and amaximum time step of tmax ¼ 50 ms. During an action potential (membranepotential EM > 50 mV) or when an external stimulation is present, themaximum time step is limited to tmax,AP ¼ 10 ms.
The change of the time step depends on the number of Newton iterations
itk needed to complete the previous time-step k:8< dtk  1:1
dtkþ1 ¼
dt =1:2: kdtk

itk <10^itk %itk1
itk >30
else:

(12)

and, therefore, the residual for the full system reads


R ¼


RNP
:
RP

(11)

By applying the method of lines, each equation is discretized in space firstby representing the unknown functions ni and f as well as the test functions
v by Q1 nodal basis functions on a tensor grid, and then in time using the

Coupling with the Hodgkin-Huxley equations
As described in the section Membrane Flux, the dynamic channel conductances from the Hodgkin-Huxley scheme are needed to calculate the total
membrane flux fimemb. For each channel type i, one additional ordinary differential equation per gating particle has to be solved in each time step. Thisis done using a simple Implicit Euler step. The membrane flux is calculated
Biophysical Journal 105(1) 242–254

246

Pods et al.

once at the beginning of each time step and kept fixed for the whole Newton
iteration. This way, the unknowns from the HH scheme do not enter thematrix for the PNP system, avoiding convergence issues arising fromchanging boundary conditions. However, this introduces a small splitting
error of O(dt).

Implementation
The implementation was done in Cþþ using the DUNE framework (18,19)
and the DUNE discretization module PDELab (20). Additionally, the
DUNE modules dune-multidomaingrid (21) and dune-multidomain (22)were used. The dune-multidomaingrid provides a metagrid that allows forthe partition of the grid into subdomains (Uelec and Umemb, in this case);dune-multidomain is an add-on to PDELab that lets the user define different
operators on different subdomains and specify coupling between subdomains. In other words, different equations can be solved on different partsof the grid, which is a highly suitable tool for our case of a nonlinearly
coupled system of PDEs defined on different domains.
Essentially, there has to be one operator for the assembly of the residual
expressions from Eqs. 9a, 9b, and 10. To couple Eqs. 1 and 2, however, theresiduals for the spatial part of the subdomain Uelec are treated together by asingle operator assembling the combined residual

RUelec ¼


RNP;S:
RP;Uelec

Here, RP;Uelec denotes those entries of the full Poisson residual Rp that comesfrom elements belonging to the Uelec subdomain.
The remaining entries RP;Umemb of the full Poisson residual are assembled
by a separate operator taking only contributions from the subdomain Umemb
into account.
A third operator handles the temporal part of the residual, RNP,T. A combined operator is obtained automatically by the DUNE module dunemultidomain, which can be used by PDELab’s Newton implementation toassemble the full residual to be minimized in each iteration.

RESULTS
Validation by analytical solutions
To test the numerical algorithms and to validate the implementation, a number of simulations have been run andcompared to analytical solutions. We used recently foundunsteady analytical solutions to the PNP equations (23)for the case of a single electrolyte domain. The comparisonsin one and two spatial dimensions showed the expectedorder of convergence in space and time (regarding ournumerical implementation) to the exact solutions.
Simulation parameters
In the following, we consider a square computational
domain of dimension 10  10 mm, where the membrane
ranges from ymemb ¼ 500 ms to ymemb þ dmemb, and themembrane thickness was chosen to be dmemb ¼ 5 nm.
For the mesh width, a uniform spacing of hx ¼ 100 mm
was used in the x direction. The Debye length for this electrolyte setup is ~0.8 nm, so we chose a minimum grid
spacing of hymin ¼ 0.5 nm at the membrane in the y direction
to account for this. A geometric spacing was used, smoothly
Biophysical Journal 105(1) 242–254

increasing the grid spacing in the y direction with increasing
distance from the membrane boundary to a maximum of
hymax ¼ 100 mm. The large difference between these lengths
underlines the multiscale character of this model.
The diffusion coefficients Di were chosen to be the diffusivity in water for each ion species. The relative permittivity
ε was 80 in the electrolytes and 2 on the membrane, in
accordance with Lu et al. (15). The temperature was fixedat T ¼ 6.3 C. Table S1, Table S2, and Table S3 containthe values chosen for initial concentrations and channel
parameters.
For these simulation parameters with 73,124 unknowns
and a simulated time of 20 ms, and using an average timestep size of 13.65 ms, a total computation time of ~25 h
was needed.
Equilibrium
When selectively opening only the leak channels for one ion
species, the system’s equilibrium membrane potential is
expected to be equal to the corresponding ionic reversal
potential, as predicted by Eq. 7. Table S3 shows the calculated equilibrium potentials, which indeed match the value
calculated by Nernst’s equation.
When opening both Na and K leak channels, the equilibrium membrane potential will reach a value lying between
the two channels’ reversal potentials. The relative leak conductances that result in a resting potential of ~65 mV
can be found in Table S3. The total leak conductance(0.5 mS/cm2) was always kept constant.
Fig. S1 in the Supporting Material shows the intra- andextracellular charge density profile at equilibrium. As
predicted by the Poisson-Boltzmann expression in Eq. 8,
both electrolytes adjust their concentrations to follow a
Boltzmann distribution toward the membrane. Fig. S2shows the evolution of membrane fluxes during the equilibration phase. The sum of inward- and outward-directed
fluxes tends to zero, marking the neuron’s resting state.
Action potential
Membrane potential

To evoke an action potential (AP), a sodium rectangle pulseof 0.965 nA is injected into the cell by adding a fixedamount of sodium at the stimulation site xstim ¼ (150 mm,
0 mm) for 2 ms. Due to the ion channel kinetics from
Coupling with the Hodgkin-Huxley Equations section, apotential wave travels along the axon, opening more channels along the way and keeping the AP alive. The conductance velocity is depending on the time constants of theion channel kinetics, but also on the intra- and extracellular
ion diffusion coefficients, and had a value of ~0.93 m/s forthis setup. Fig. 2 shows the potential time courses atdifferent positions along the axon. The first AP has a higher

Electrodiffusion Models of Axon and Extracellular Space247
x=50µmx=2050µm
x=4050µm
x=6050µm
x=8050µm

4020Potential [mV]

0

−20

−40

−60

−80
0

2

4

6

8

10
Time [ms]

12

14

16

18

20

FIGURE 2 Action potentials (APs) evoked at equidistant positions along the axon. The leftmost AP has a higher amplitude than the others due to thevicinity to the stimulus site. The following curves at equidistant positions along the axon are identical and show constant onset delay, indicating a wave
traveling with constant speed.

amplitude than the following APs, caused by theproximity to the stimulation site. Also, switching off thestimulus is reflected by an artifact in the repolarization
phase of the first AP. In Movie S1 in the Supporting
Material, the AP propagation along the axon in space andtime can be seen.
Extracellular potential

We now focus on the time evolution of local field potential(LFP) signals at any point in the extracellular domain. In
Fig. 3, the potential time courses are plotted for the samex coordinate, at increasing distances from the membrane.
Some major features can be identified from these curves:
A first positive peak (P1) followed by a larger negative
peak (N1), then a (very) small second positive peak (P2)with a subsequent longer phase of slowly varying potential
with negative curvature (S), and a last peak (P3). This
characteristic up-down-up shape is maintained at various
distances from the membrane. The potential time course
after P1 seems to largely follow the total membrane flux
at the same x coordinate (Fig. 4), in accordance with the
LSA model. However, the small positive deflection at the
beginning of the membrane flux curve alone cannot explain

the large amplitude (comparable in magnitude to N1) of P1.
We will now try to explain this difference by looking at
snapshots of the extracellular potential and ion concentrations at a fixed time point. Because during an AP thisprofile simply moves through space with a constant (known)
velocity (~0.93 m/s in this simulation), the complete information about the potential and the LFP dynamics at any
point in space can be gained from such plots. As the signal
moves in positive x direction for all following snapshots, weread from right to left as opposed to Fig. 3.
Fig. 5 a shows the spatial structure of the extracellularpotential profile within the Debye layer, i.e., at most a few
nanometers from the membrane. It is almost exclusively
dominated by the intracellular potential (compare to
Fig. 5 b), meaning that the intracellular potential spreads
across the membrane into the extracellular space, even if
only with a greatly reduced amplitude. We call this theecho of the AP.
The near-field potential profile (a few nanometers to ~10mm from the membrane) in Fig. 6 (see also Movie S2) showsthat this echo of the AP has a significant effect on the LFP in
the extracellular domain. The profile begins to the right witha rise in the potential (corresponding to P1) followed by a
Biophysical Journal 105(1) 242–254

248

Pods et al.
x 10

−3

P1

1

0.5

Potential [mV]

0
P3
P2

S

−0.5

−1

0.62597 [µm]11.0227 [µm]26.1498 [µm]34.8784 [µm]46.5215 [µm]62.0522 [µm]82.7686 [µm]110.4022 [µm]

−1.5

N1

3

4

5

6

7
Time [ms]

8

9

10

FIGURE 3 Extracellular potential at various distances from the membrane. The potential time courses for a fixed x coordinate x ¼ 5.05 mm is shown forseveral different y positions. Although the characteristic up-down-up shape is maintained, the amplitude spans several orders of magnitude as the signal isstrongly attenuated with distance.

sharp drop (N1) and another rise (P2). After this first phase,the potential has a longer phase of low variation (S) until
another, less pronounced peak is observed at the rear endof the traveling action potential (P3).
x 10

−3

5

0

Ion flux [mol/m2 s]

−5

−10

−15

−20
Total membrane flux, x=5100 [µm]
3

4

5

6

7
Time [ms]

8

9

10

FIGURE 4 Total membrane flux at fixed point at the membrane. The
membrane flux at the same x coordinate as the potential curves from
Fig. 3 is plotted over time. The influence of membrane flux on the potential
can be clearly seen, but it cannot explain the first peak in the potential.
Biophysical Journal 105(1) 242–254

When comparing this snapshot profile with the profile
immediately at the membrane (Fig. 5 a) by shape andposition, it can be observed that the single peak from
Fig. 5 a we termed the AP echo has been divided into twoparts P1 and P2 in Fig. 6 by an interrupting negative peakN1. The value N1 is a result of the opening of voltage-gatedsodium channels and the following massive depletion of
sodium, as observable in Fig. 7 a: The total charge density(compare to Fig. 7 d) follows the decrease in sodium andgenerates the potential drop N1. Roughly speaking, theactivation of sodium channels and the resulting negative
peak N1 splits up the single positive peak from the APecho into two peaks P1 and P2.
The difference between the Debye layer (Fig. 5 a) andnear-field (Fig. 6) potential profiles is striking. The former
is dominated by the AP echo; the latter shows approximately equal contributions of membrane currents and intracellular AP. This can be explained by the faster attenuationof the AP echo over the Debye layer (note the strongly
reduced amplitude of P1 from Fig. 5 a to Fig. 6) in relation
to the potential evoked by membrane fluxes.
It is also notable that the second peak P2 is attenuatedmuch more quickly than P1 and can barely even be seenElectrodiffusion Models of Axon and Extracellular Space

249

a

b

ðESÞ

FIGURE 5 Snapshot of the Debye layer extracellular potential profile. (a) The potential values for a narrow stripe of the extracellular domain Uelec(corresponding to the shaded area from Fig. 1) just above the membrane are plotted at a fixed time. This turns out to be the intracellular AP (b) propagatingover the membrane, with a significantly reduced amplitude.

at a distance of only a few micrometers away from themembrane. Instead, the potential largely follows the membrane flux from Fig. 4, which becomes the prevailingdriving force after its activation, as could already be seenin Fig. 3. This does not necessarily mean that the potential
echo only has an effect on the extracellular field until thechannels are opened, but rather that the rest of the APecho (i.e., the falling edge) is hidden in the LFP signal, as
it is superimposed by the negative potential generated
by activating membrane currents. In this setup, the APadvances the activation of channels by a few hundred
microseconds. A different channel timing might change
the resulting LFP waveform to a great extent.
The distance between the end of P1 and the beginning of
P3 close to the membrane (as in Fig. 6) is a good measure for

2

x 10

−3

the timescale of the extracellular field that we term the ‘‘LFP
valley length’’ (the region with negative potential values). It
also gives a characteristic length scale for the range ofsimultaneous ion channel activity along the axon. In our
simulation, this length is ~2000 mm at the membrane.
Fig. 8 and Movie S3 show the potential for a large part ofthe extracellular space, demonstrating that the near-fieldpotential profile essentially continues into distant space,albeit attenuating quickly with distance. The LFP valley
length (the diameter of the green area) increases notably
with distance, which we attribute to the diffusive characterof the system. The general pattern of the LFP stays thesame: a positive upwind domain (P1) just in front of the
opening channels, followed by a negative middle region(N1, S), and then again a positive rear domain (P3).

time = 5.7167 [ms]1
Potential [mV]

0
−1
−2

FIGURE 6 Snapshot of the near-field extracellular potential profile. Plotted are the potential
values for a stripe of the extracellular domain
ðESÞ
Uelec just above the extracellular Debye layer(i.e., above the shaded area from Fig. 1), at a fixed
time. It shows a more complex structure comparedto the Debye layer potential from Fig. 5.

−3
−4
−5

30
20
10
y position [µm]

0

0

2000

4000
x position [µm]

6000

8000

10000

Biophysical Journal 105(1) 242–254

250

Pods et al.

a

b

c

d

FIGURE 7 Snapshot of the Debye layer and near-field concentration profiles. The concentration profiles of all three ion species (a–c) and the charge
density (d) are plotted up to a distance of ~10 mm from the membrane. The strong concentration gradients close to the membrane demonstrate the influenceof the Debye layer. The charge density can be seen to be dominated by the sodium concentration.

The amplitudes of the potential are differing by several
orders of magnitude among the intracellular values(~100 mV), the extracellular values close to the membrane(up to a few 100 mV), and some millimeters away from themembrane (fractions of 1 mV).
The influence of potential and concentration changes on
the local charge redistribution can be observed separately
by studying the ion flux from Fig. S3 and Fig. S4, foreach ion species. It gives a good impression about the complex interplay of drift and diffusion causing ion movement
close to the membrane.
Comparison with line source approximation
model
The line source approximation introduced in Holt and Koch(7) is widely used as an effective model to compute theextracellular potential at any point in space, using only
the values of the membrane currents at a finite number ofline segments,
Biophysical Journal 105(1) 242–254

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
 h2 þ r 2  h rIlog pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
Fðr; hÞ ¼
;
4pDs  l2 þ r 2  l (13)where r is the radial distance from the line of length Ds, his the longitudinal distance from the end of the line, andl ¼ Ds þ h is the distance from the start of the line. The
parameter r describes the resistivity of the extracellular
medium. This model is very convenient, especially whenusing cable equation models based on a line segment
approximation of the original neuron geometry.
For the electrodiffusion model, the computed extracellular field can easily be compared with the LSA results,
as the membrane currents are explicitly known from theinternal membrane flux boundary conditions, so they canbe inserted into the LSA expression from Eq. 13 directly.
The resistivity r was chosen manually in such a way thatthe negative peaks N1 of electrodiffusion and LSA match
for a fixed position close to the membrane (z1 mm in thiscase). This resulted in a value of r ¼ 37.4 Ucm,

Electrodiffusion Models of Axon and Extracellular Space

251

Snapshot time = 5.7167 [ms]
+1e−2
+1e−3

7000

+1e−4
+1e−5
+1e−6

6000

+1e−7
+1e−8

y position [µm]

5000

+1e−9
+1e−10+1e−11

4000

−1e−1
−1e−2

3000

−1e−3
−1e−4
−1e−5

2000

−1e−6
−1e−7

1000

−1e−8
−1e−9

1000

2000

3000

4000

5000
x position [µm]

6000

7000

8000

9000

FIGURE 8 Snapshot of the whole extracellular potential. A contour plot of the logarithmically scaled extracellular potential at a fixed time point. The
domain was cut at the left, right, and upper boundary to exclude artifacts introduced by the boundary where the potential is almost zero and switches signs
due to small numerical errors.

corresponding to a conductivity of s ¼ 2.67 S/m, which is~30–50% higher than experimentally obtained values for
cerebrospinal fluid (24,25). Fig. 9 shows the comparisonof potential time courses for LSA and electrodiffusion(ED) models at various distances from the membrane. Additionally, the difference between the two curves is plotted.
One can see that both signals show a set of common
features: The negative peak N1 following the dominating
sodium current in the beginning of an AP, followed by atransition phase of slower variation (S) and the less pronounced positive peak P3 during the repolarization phase.
But the electrodiffusion signal shows another sharp positive
peak (P1) right in the beginning, which is not present in the
LSA model.
The difference between the LSA and ED curves resembles a distorted AP waveform that gets smeared out overdistance. Our explanation for the difference is that the ED
potential results from two contributors: The potential generated by membrane currents (as in LSA) and, in addition, an
AP echo part directly coming from the intracellular potential wave traveling along the axon. The difference plottedin Fig. 9 would then largely represent the latter part of thesignal stemming from the AP echo alone.

Looking closely, there is also a slight delay in the ED
signal—i.e., the N1 peak is time-shifted by ~90 ms to the
LSA signal, with the same holding for the P3 peak. These
differences cannot be a result of the membrane current sources, as these are identical for both models, so it has to be aneffect of the extracellular medium. We attribute the delay tothe local effects of redistributing ion concentrations,especially in the Debye layer with its large concentration
gradients, leading to a delayed transduction of the electricfield, wheras the LSA model assumes an instantaneous
response to the membrane current sources.
CONCLUSIONSIn this article we presented results from the numericalsimulation of an AP traveling along an active axonal membrane and its spread into the extracellular space. To ourknowledge, this is the first time that such detailed information about the dynamics of ion concentrations and the electric potential in time and space were given on this scale.
Recent applications of PNP theory to neurons have been
done in Lopreore et al. (11) and Mori et al. (13). In contrast
to these studies, we paid special attention in resolving the
Biophysical Journal 105(1) 242–254

252

Pods et al.

a

b

c

d

e

f

FIGURE 9 Comparison of electrodiffusion and line source approximation. The time courses of extracellular potentials calculated by electrodiffusion(ED, dashed line) and line source approximation (LSA, solid line) models are compared at different distances from the membrane (a–f) for a fixed x coordinate.
The similar shape for the part between N1 and P3 is apparent, but the first peak P1 is missing in LSA. (Dash-dotted line) Difference between the two curves.

Debye layer explicitly, thereby accounting for the effects oflarge concentration gradients. This also adds the benefit thatboth steady-state concentrations and equilibriumchannel reversal potentials can be validated directly, using
Poisson-Boltzmann and Nernst equations, respectively.
The evolution of the LFP signal and its various features
have been analyzed particularly. We showed that the large
concentration gradients have a significant impact on thenear-field extracellular field. A detailed evaluation of the
Biophysical Journal 105(1) 242–254

results shows the complex interplay of potential and concentration changes on the local ion redistribution. The
intricate structure of the results (compare to Fig. 6) reveals
that the neural membrane dynamics—even in this very
simple case of a minimal set of only two channel
types—is highly nontrivial and that a detailed study is
indeed justified.
The main finding is that the electrodiffusion modelshows significant deviations from the well-establishedElectrodiffusion Models of Axon and Extracellular Space

line source approximation model, in particular a positive
peak at the beginning of the LFP signal, directly causedby the intracellular potential propagating into extracellularspace. This implies that there is a second component next
to membrane currents contributing to the generation of the
LFP signal, which we termed the ‘‘AP echo’’. This component had a strong effect in our simplified model, but
this might look differently in other situations, e.g., in
myelinated fibers. In Gold et al. (5), the LSA showed agood agreement with experimental data, although theauthors reported difficulties in fitting the intra- andextracellular data simultaneously. However, there is also
experimental evidence that the LFP waveforms emerging
from the electrodiffusion model occur in unmyelinatedafferent nerve fibers (26) as well as in hippocampal CA1
neurons (27).
It remains to be shown in which situations the LSA model
gives a valid approximation, and under which conditions theeffect of the AP echo on the extracellular potential cannot be
neglected. An experimental study using juxta- and extracellular recordings of axon fibers should be able to validate thereported LFP attenuation and change in signal shape fromthe Debye layer to near-field and the more distant extracellular regime.
We emphasize that our model could numerically be simulated with standard methods, i.e., with conforming finite
elements for the spatial and an implicit Euler method forthe time discretization. However, some essential precautions
had to be taken to ensure a stable and efficient method,
namely:The grid had to resolve the Debye layer close to the membrane in the direction normal to the membrane. For thex direction, a much coarser mesh size was sufficient.
The PNP system had to be solved in a fully coupledfashion, as an operator-splitting approach dramatically reduced the time-step size needed for computing
a nonoscillatory solution.
The system had to be carefully equilibrated before setting
a stimulus, as the concentration profile toward themembrane had to reach steady state to get meaningful
reversal potentials for each ion channel.
The choice of an implicit time-stepping scheme took
advantage of higher stable time-step values andaccounted for the diffusion-dominance in the NernstPlanck equation that we observed for parameters in
the physiological range.
The application of a threshold volume scaling was
crucial to compensate the large differences of residual
magnitudes introduced by the drastically varying cellvolumes in this cylinder geometry. With this, theNewton iteration was able to converge even for large
domain sizes—using the same error tolerances as for aplain two-dimensional simulation—with an only
slightly lower average time-step size.

253

The model used here does not claim to be fully realistic
nor complete from the biological point of view. Instead it
should be considered as a first attempt to model the
dynamics of neural systems on this spatial scale. Particularly, the idealized setup of a single axon fiber in an extracellular space consisting exclusively of fluid will rarely occur
in reality. Additionally, no synaptic currents were regarded,
which are thought to be the main contributors of electroencephalography and LFP signals (28,29). However, the characteristic features of the LFP signal could be explained by
the interplay of ion channel currents and electrolyte
dynamics alone in this model.
To arrive at more realistic models, the greatly simplified
geometry from this model will have to be replaced in a later
study by either explicitly incorporating a complex reconstructed extracellular geometry, or by finding effective
material parameters (diffusivity of ions, electric permittivity) for representative intra- and extracellular media
without explicitly modeling nearby cells.
For the axon geometry, assessing the impact of branching
structures and ultimately, a complete attached single cellmorphology with hundreds of compartments and synapses,
is of great interest. Because these complicated geometries
cannot be handled by a cylinder symmetry, a full threedimensional simulation has to be carried out for this, necessitating the parallelization of our simulator.
Further model refinement concerning the detailed
structure and properties of the active membrane (channel
types and conductivities, surface charges, myelination)
will have to be done when comparing the results with experimental data.
SUPPORTING MATERIAL
Supporting analysis, three tables, four figures, reference (30), and threemovies are available at http://www.biophysj.org/biophysj/supplemental/
S0006-3495(13)00628-0.
We thank Prof. Andreas Draguhn and Dr. Martin Both for fruitful discussions and information about extracellular signals in slice recordings.
Many thanks also to Steffen Müthing for helping with the dune-multidomain setup and the anonymous reviewers for their valuable comments.
This work was funded through a grant from the German Ministry of
Education and Research (BMBF No. 01GQ1003A).

REFERENCES
1. Hodgkin, A. L., and A. F. Huxley. 1952. A quantitative description of
membrane current and its application to conduction and excitation in
nerve. J. Physiol. 117:500–544.
2. Rall, W. 1989. Cable theory for dendritic neurons. In Methods in
Neuronal Modeling. MIT Press, Cambridge, MA. pp. 9–92.
3. Hines, M. L., and N. T. Carnevale. 1997. The NEURON simulation
environment. Neural Comput. 9:1179–1209.
4. Bower, J., D. Beeman, and A. Wylde. 1995. The Book of GENESIS:
Exploring Realistic Neural Models with the General Neural Simulation
System. Telos Publishing, Candor, New York.
Biophysical Journal 105(1) 242–254

254

Pods et al.

5. Gold, C., D. A. Henze, ., G. Buzsáki. 2006. On the origin of the extracellular action potential waveform: a modeling study. J. Neurophysiol.
95:3113–3128.

18. Bastian, P., M. Blatt, ., O. Sander. 2008a. A generic grid interface forparallel and adaptive scientific computing. I. Abstract framework.
Computing. 82:103–119.

6. Bédard, C., H. Kröger, and A. Destexhe. 2004. Modeling extracellular
field potentials and the frequency-filtering properties of extracellularspace. Biophys. J. 86:1829–1842.
7. Holt, G. R., and C. Koch. 1999. Electrical interactions via the extracellular potential near cell bodies. J. Comput. Neurosci. 6:169–184.
8. Bédard, C., H. Kröger, and A. Destexhe. 2006. Model of low-passfiltering of local field potentials in brain tissue. Phys. Rev. E Stat.
Nonlin. Soft Matter Phys. 73:051911.

19. Bastian, P., M. Blatt, ., O. Sander. 2008b. A generic grid interface forparallel and adaptive scientific computing. II. Implementation and tests
in DUNE. Computing. 82:121–138.
20. Bastian, P., F. Heimann, and S. Marnach. 2010. Generic implementation of finite element methods in the distributed and unified numerics
environment (DUNE). Kybernetika. 46:294–315.

9. Bédard, C., and A. Destexhe. 2009. Macroscopic models of local fieldpotentials and the apparent 1/f noise in brain activity. Biophys. J.
96:2589–2603.
10. Pettersen, K. H., E. Hagen, and G. T. Einevoll. 2008. Estimation of
population firing rates and current source densities from laminar electrode recordings. J. Comput. Neurosci. 24:291–313.

21. Müthing, S., and P. Bastian. 2012. Dune-multidomaingrid: a metagrid
approach to subdomain modeling. In Advances in DUNE. A. Dedner,
B. Flemisch, and R. Klöfkorn, editors. Springer, Berlin, Germany.
pp. 59–73.
22. Müthing, S. 2012. Dune-multidomain. Accessed December 4, 2012.
http://users.dune-project.org/projects/dune-multidomain.
23. Schönke, J. 2012. Unsteady analytical solutions to the Poisson-NernstPlanck equations. J. Phys. A. 45:455204.

11. Lopreore, C. L., T. M. Bartol, ., T. J. Sejnowski. 2008. Computationalmodeling of three-dimensional electrodiffusion in biological systems:
application to the node of Ranvier. Biophys. J. 95:2624–2635.

24. Baumann, S. B., D. R. Wozny, ., F. M. Meno. 1997. The electricalconductivity of human cerebrospinal fluid at body temperature. IEEE
Trans. Biomed. Eng. 44:220–223.

12. Mori, Y. 2009. From three-dimensional electrophysiology to the cablemodel: an asymptotic study. Arxiv preprint. arXiv: 0901.3914.

25. Gabriel, S., R. W. Lau, and C. Gabriel. 1996. The dielectric propertiesof biological tissues: III. Parametric models for the dielectric spectrumof tissues. Phys. Med. Biol. 41:2271–2293.

13. Mori, Y., G. I. Fishman, and C. S. Peskin. 2008. Ephaptic conduction ina cardiac strand model with 3D electrodiffusion. Proc. Natl. Acad. Sci.
USA. 105:6463–6468.
14. Xylouris, K., G. Queisser, and G. Wittum. 2010. A three-dimensional
mathematical model of active signal processing in axons. Comput.
Visualization Sci. 13:409–418.
15. Lu, B., M. J. Holst, ., Y. C. Zhou. 2010. Poisson-Nernst-Planck
equations for simulating biomolecular diffusion-reaction processes.
I. Finite element solutions. J. Comput. Phys. 229:6979–6994.
16. Koch, C. 2004. Biophysics of Computation: Information Processing inSingle Neurons. Oxford University Press, Oxford, UK.
17. Demmel, J. W., S. C. Eisenstat, ., J. W. H. Liu. 1999. A supernodal
approach to sparse partial pivoting. SIAM J. Matrix Anal. Appl.
20:720–755.

Biophysical Journal 105(1) 242–254

26. Weidner, C., R. Schmidt, ., H. E. Torebjörk. 2000. Time course of
post-excitatory effects separates afferent human C fiber classes.
J. Physiol. 527:185–191.
27. Csicsvari, J., H. Hirase, ., G. Buzsáki. 1999. Oscillatory coupling of
hippocampal pyramidal cells and interneurons in the behaving Rat.
J. Neurosci. 19:274–287.
28. Buzsáki, G., C. A. Anastassiou, and C. Koch. 2012. The origin of
extracellular fields and currents—EEG, ECoG, LFP and spikes. Nat.
Rev. Neurosci. 13:407–420.
29. Nunez, P., and R. Srinivasan. 2005. Electric Fields of the Brain: theNeurophysics of EEG. Oxford University Press, Oxford, UK.
30. Carnevale, N. 2012. Vtrap derived. Accessed December 4, 2012. http://
www.neuron.yale.edu/phpBB/viewtopic.php?f¼15&t¼1075.