# Finite Element Simulation of Ionic Electrodiffusion in Cellular Geometries

^{1}Department for Scientific Computing and Numerical Analysis, Simula Research Laboratory, Oslo, Norway^{2}Centre for Integrative Neuroplasticity, University of Oslo, Oslo, Norway^{3}Department of Physics, University of Oslo, Oslo, Norway^{4}Faculty of Science and Technology, Norwegian University of Life Sciences, Ås, Norway

Mathematical models for excitable cells are commonly based on cable theory, which considers a homogenized domain and spatially constant ionic concentrations. Although such models provide valuable insight, the effect of altered ion concentrations or detailed cell morphology on the electrical potentials cannot be captured. In this paper, we discuss an alternative approach to detailed modeling of electrodiffusion in neural tissue. The mathematical model describes the distribution and evolution of ion concentrations in a geometrically-explicit representation of the intra- and extracellular domains. As a combination of the electroneutral Kirchhoff-Nernst-Planck (KNP) model and the Extracellular-Membrane-Intracellular (EMI) framework, we refer to this model as the KNP-EMI model. Here, we introduce and numerically evaluate a new, finite element-based numerical scheme for the KNP-EMI model, capable of efficiently and flexibly handling geometries of arbitrary dimension and arbitrary polynomial degree. Moreover, we compare the electrical potentials predicted by the KNP-EMI and EMI models. Finally, we study ephaptic coupling induced in an unmyelinated axon bundle and demonstrate how the KNP-EMI framework can give new insights in this setting.

## 1. Introduction

The most common computational models for excitable cells are those based on cable theory (Rall, 1977; Koch, 1999). In its standard form, the cable model is based on several simplifying assumptions, most importantly that the extracellular potential and both intracellular and extracellular ion concentrations are constant in space and time. Multi-compartmental neuron models based on cable theory are widely used within the field of neuroscience to simulate large network of interacting neurons (see e.g., Markram et al., 2015). In such models, only synaptic interactions between neurons are considered, whereas changes in the extracellular field and extracellular ion concentrations associated with a neuron's activity are assumed to be too small to have any influence on its neighboring neurons (or itself). Although these assumptions are only approximations, the resulting models still give accurate predictions of neuronal electrodynamics in many scenarios. Indeed, concentration changes are often limited by neuronal and glial uptake mechanisms that strive toward maintaining concentrations close to basal levels.

However, there are also many scenarios that involve dramatic changes in extracellular ion concentrations. On a large spatial scale, ion concentration changes are a trademark of several pathological conditions, such as spreading depression or epilepsy (Dietzel et al., 1989; Somjen, 2001; Syková and Nicholson, 2008; Ayata and Lauritzen, 2015). Extracellular concentration shifts will lead to changes in neuronal reversal potentials, and can thus affect the dynamical properties of the neurons (Kager et al., 2000; Øyehaug et al., 2012; Wei et al., 2014). Under non-pathological conditions, concentration-dependent, electrodiffusive effects are hypothesized to be important in specific microdomains of the brain (Savtchenko et al., 2017). In general, the extracellular ion concentration changes resulting from a neuronal event can be expected to be largest in regions where the extracellular space is small and confined.

Similarly, there are several scenarios where the assumption of a constant extracellular potential may be questionable. For instance, ephaptic interactions have been reported to play a role for neural phenomena taking place at both small and large spatial scales (Holt and Koch, 1999; Bokil et al., 2001; Anastassiou et al., 2011; Anastassiou and Koch, 2015; Goldwyn and Rinzel, 2016; Tveito et al., 2017b; Han et al., 2018; Shifman and Lewis, 2019). Ephaptic interaction (or coupling) is a coupling between neurons via the extracellular potential, which is hard or impossible to represent under the aforementioned assumption.

The olfactory nerve is one example in which variations in ion concentrations and extracellular potentials may be important. Whereas most axons in the mammalian brain are coated in an insulating layer of *myelin*, the axons in the olfactory nerve are unmyelinated and organized in tight bundles (Doucette, 1984; Griff et al., 2000). In view of the tight packing, one might expect large ion concentration variations in the extracellular space between the olfactory nerve axons. Moreover, the olfactory nerve axon arrangement will maximize any ephaptic coupling, with a potential evolutionary purpose (Lowe, 2003). In addition, diffusion along extracellular ion concentration gradients can generate so-called *diffusion potentials* (Halnes et al., 2016; Savtchenko et al., 2017; Solbrå et al., 2018), which may constitute an additional ephaptic effect on membrane potentials.

There are several computational studies considering ephaptic interaction in the brain. Bokil et al. (2001) use a simplified model based on cable theory, and find that an action potential in a single axon can evoke action potentials in neighboring axons. A more detailed model for coupling intra- and extracellular currents is the *Extracellular-Membrane-Intracellular* (EMI) model (Krassowska and Neu, 1994; Ying and Henriquez, 2007; Agudelo-Toro, 2012; Agudelo-Toro and Neef, 2013; Tveito et al., 2017a,b). The EMI model incorporates explicit 3D shapes of the neuron, allowing for morphologically detailed descriptions of the neuropil. However, neither of the aforementioned frameworks explicitly model the ion concentrations and can therefore not capture ephaptic effects due to electrodiffusion, such as diffusive potentials.

The most physically detailed scheme for modeling electrodiffusion is the *Poisson-Nernst-Planck* (PNP) framework (Lopreore et al., 2008; Pods et al., 2013; Holcman and Yuste, 2015; Cartailler et al., 2017a,b; Sacco et al., 2017). The PNP framework is based on explicitly simulating charge relaxation processes taking place at small spatiotemporal scales (~nm and ~ns), and thus requires high resolutions in both time and space. Consequently, applications have been limited to studying dynamics at the ion channel and cell membrane level. An alternative approach is to assume that the bulk tissue is electroneutral, thus circumventing the need for explicit modeling of charge relaxation processes. Models based on the electroneutrality assumption are therefore numerically stable for coarser spatial and temporal resolutions, allowing for longer simulations on larger domains.

On this background, a series of electroneutral models for ionic electrodiffusion have been developed, both for homogenized domains (Mori et al., 2008; Halnes et al., 2013, 2016, 2017; Niederer, 2013; Pods, 2017; Solbrå et al., 2018), and for domains including an explicit geometrical representation of the cells and of the extracellular space (Mori and Peskin, 2009). In particular, Mori and Peskin (2009) presents a finite volume method for solving a system of equations describing cellular electrical activity accounting for both geometrical effects and ion concentration dynamics.

In this paper, we present a variation of the Mori and Peskin (2009) model and introduce a mortar-based finite element formulation of this model. Key advantages of the finite element formulation are (i) the independence of dimension: the same scheme is applicable for one-, two-, or three-dimensional domains (with zero-, one-, or two-dimensional cell membranes/interfaces); (ii) the handling of complicated interface geometries; and (iii) the straightforward use of more accurate, i.e., higher order polynomial schemes. The framework can be viewed as a combination of the EMI framework and the electroneutral *Kirchhoff-Nernst-Planck* (KNP) framework (Solbrå et al., 2018), and will henceforth be referred to as the KNP-EMI framework. Previous numerical schemes for the KNP framework are restricted to simplified 1D geometries (Halnes et al., 2013, 2015; Sætra et al., 2020), or components within a hybrid modeling scheme to compute extracellular dynamics (Halnes et al., 2016, 2017; Solbrå et al., 2018).

The KNP-EMI framework can be viewed as an extension of the EMI framework by the explicit modeling of ion concentrations and the effects of ionic electrodiffusion. We here evaluate the effect of these extensions by comparing the KNP-EMI and EMI solutions in idealized axon domains, and find that the solutions are qualitatively similar but differ locally. However, the KNP-EMI simulations give further insights into the importance of extracellular bulk conductivities for ephaptic couplings in neural tissue: KNP-EMI simulations of idealized, unmyelinated axon bundles reveal increased extracellular bulk conductivities and, as a result, a reduced tendency toward induction of action potentials in neighboring axons.

## 2. Methods

We present the governing equations for ionic electrodiffusion in neural tissue with a geometrically explicit representation of the cellular membranes in section 2.1 below. To take full advantage of this framework, a numerical solution scheme capable of efficiently handling three-dimensional, complicated geometries is required. We here propose a novel numerical solution scheme using a mortar finite element method (Bernardi et al., 1993; Agudelo-Toro and Neef, 2013) and a two-step splitting scheme, described in section 2.2. This solution algorithm flexibly allows for arbitrary geometries and efficient solution of the separate subproblems. Our implementation of this algorithm is openly available at: https://zenodo.org/record/3492075#.XahQOhh9g5k.

### 2.1. A Mathematical Framework for Electrodiffusion With Explicit Membrane Representation

#### 2.1.1. Representation of the Computational Domain

We consider *N* domains ${\Omega}_{{i}^{n}}\subset {\mathbb{R}}^{d}$ (*d* = 1, 2, 3) for *n* = 1, …, *N* representing disjoint intracellular regions (physiological cells, e.g., neurons) and an extracellular region Ω_{e}, and let the complete domain $\Omega ={\Omega}_{{i}^{1}}\cup \cdots \cup {\Omega}_{{i}^{N}}\cup {\Omega}_{e}$ with boundary ∂Ω. See Figure 1 (Right) for an illustration of a sample domain configuration. We denote the cell membrane associated with cell *i*^{n}, i.e., the boundary of the physiological cell ${\Omega}_{{i}^{n}}$, by Γ_{n}. We assume that Γ_{n} ∩ Γ_{m} = ∅ for all *n*≠*m* and that Γ_{n} ∩ ∂Ω = ∅ (It follows that $\partial {\Omega}_{{i}^{n}}\cap \partial {\Omega}_{e}=\varnothing $ for all *n* = 1, …, *N*.). For simplicity and clarity, we present the mathematical model for one intracellular region ${\Omega}_{{i}^{1}}={\Omega}_{i}$ with membrane Γ below. The extension to multiple intracellular regions is immediate (but notationally cumbersome).

**Figure 1**. Overview of the computational domains. **(Left)** Idealized axon bundle consisting of 9 cuboid-shaped axons. **(Middle)** Cross-section of the axon bundle, where the axons are labeled with repeated labels for symmetric positions. **(Right)** Idealized 2D computational domain with one intracellular region Ω_{i} and extracellular region Ω_{e}.

#### 2.1.2. Intracellular and Extracellular Governing Equations

We will here derive a system of coupled, time-dependent, non-linear partial differential equations to describe ionic electrodiffusion in this domain. We consider a set of ion species *K*. Typically *K* will include sodium *Na*^{+}, potassium *K*^{+}, and chloride *Cl*^{-}. For each ion species *k* ∈ *K* and each region *r* ∈ {*i, e*}, we model the *ion concentrations* [*k*]_{r} : Ω_{r} × (0, *T*] → ℝ (mol/m^{3}) and the *electrical potentials* ϕ_{r} : Ω_{r} × (0, *T*] → ℝ (V). Conservation of ions for the bulk of each region Ω_{r} stipulates that

for *t* ∈ (0, *T*]. Here, ${\text{J}}_{r}^{k}:{\Omega}_{r}\times \left(0,T\right]\to {\mathbb{R}}^{d}$ is the regional ion flux density [mol/(m^{2}s)] of ion *k*. To proceed, we invoke the KNP assumption of bulk electroneutrality. In this case, the ion flux densities ${\text{J}}_{r}^{k}$ satisfy:

where *z*^{k} is the valence of ion species *k* and *F* is Faraday's constant. The assumption (2) states that the total net flow of ions (weighted by the respective valences) out of any infinitesimal representative bulk volume is zero. Furthermore, we assume that the each regional ion flux density can be expressed by a Nernst-Planck equation as follows:

Here, ${D}_{r}^{k}$ denotes the effective diffusion coefficient (m^{2}/s) of ion species *k* in the region *r*. The constant ψ = *RTF*^{−1} combines Faraday's constant *F*, the absolute temperature *T*, and the gas constant *R*. The ion flux density, i.e., the flow rate of ions per unit area, is thus modeled as the sum of two terms: (i) the diffusive movement of ions due to ionic gradients $-{D}_{r}^{k}\nabla {\left[k\right]}_{r}$ and (ii) the ion concentrations that are transported via electrical potential gradients, i.e., the ion migration $-{D}_{r}^{k}{z}^{k}{\psi}^{-1}{\left[k\right]}_{r}\nabla {\varphi}_{r}$ where ${D}_{r}^{k}{\psi}^{-1}$ is the electrochemical mobility. This model ignores convective effects, and thus assumes that the underlying material (typically fluid) is at rest. As the potential ϕ_{e} is only determined up to a constant in Equations (1–3)„ an additional constraint is required, e.g.,

By inserting (3) into (2) we recognize (from volume conductor theory) the following expression for the bulk conductivity σ_{r}:

Notably, the bulk conductivity σ_{r} depends on the ion concentrations [*k*]_{r} and the diffusion coefficients ${D}_{r}^{k}$. Electrodiffusive models without explicit modeling of the ion concentrations typically set the bulk conductivity as a independent parameter (e.g., Krassowska and Neu, 1994; Bokil et al., 2001; Tveito et al., 2017a).

Inserting (3) into (1) and into (2), we thus obtain a system of *N*|*K*|+*N* equations (*N*|*K*| parabolic, *N* elliptic) for the *N*|*K*|+*N* unknown scalar fields. The system remains to be closed by appropriate initial conditions, boundary conditions, and importantly interface conditions.

#### 2.1.3. Interface Conditions

We next turn to modeling the cell membrane currents and membrane potential across the interface Γ. We denote the membrane potential ϕ_{M} as the jump in the electrical potential over the membrane:

We introduce the *total ionic current density* *I*_{M} : Γ × (0, *T*) → ℝ [C/(m^{2}s)] across the interface Γ. By definition and by conservation of total charge, we have that

where **n**_{r} denotes the boundary normal pointing out of Ω_{r} for *r* ∈ {*i, e*}. Next, we assume that *I*_{M} consists of two components: (i) a total channel current *I*_{ch} and (ii) a capacitive current *I*cap:

We further assume that the total channel current *I*_{ch} is the sum of the ion specific channel currents ${\text{I}}_{\text{ch}}^{k}$:

The channel currents ${\text{I}}_{\text{ch}}^{k}$ are subject to modeling. Typical models for ${\text{I}}_{\text{ch}}^{k}$ notably includes an synaptic input current *I*syn, leaky passive neuron, Hodgkin-Huxley etc., and will be detailed further below in section 2.1.4. On the other hand, the capacitive current *I*cap is defined over to be the capacitance *C*_{M} times the rate of change of the voltage (Sterratt et al., 2011), hence:

Inserting (10) into (8) and rearranging gives the following relation for the membrane potential ϕ_{M}:

It remains to specify a set of interface conditions for the specific ion fluxes ${\text{J}}_{r}^{k}\xb7{\text{n}}_{r}$ for *r* ∈ {*i, e*}. Here, we propose a heuristic approach via ion specific capacitive current modeling. An alternative approach is presented in Mori and Peskin (2009). As for the total current, we assume that the capacitive current can be represented as a sum of ion-associated currents:

Without loss of generality, we let the ion specific capacitive current *I*cap, *r*^{k} in region Ω_{r} at the interface Γ be some fraction ${\alpha}_{r}^{k}$ of the total capacitive current *I*cap:

Specifically, we assume that:

and note that $\sum _{k\in K}{\alpha}_{r}^{k}=1$ for *r* ∈ {*i, e*}. By definition of the ion currents and the expression for the capacitive current given by (8), we let the intracellular and extracellular ion fluxes across the membrane be given by:

for *k* ∈ *K*.

#### 2.1.4. Modeling Specific Ion Channels

The framework presented thus far allows for general representations of the ion channel current dynamics. In particular, the framework admits different choices of ion specific channel current models ${\text{I}}_{\text{ch}}^{k}$. An advantage of the geometrically explicit framework is that it allows for different channel currents models for individual cells and e.g., geometrically heterogeneous material properties. We here summarize two examples of ion specific channel currents: a passive membrane model (Sterratt et al., 2011) and the Hodgkin-Huxley model (Hodgkin and Huxley, 1952).

##### 2.1.4.1. Passive membrane dynamics

We model the passive membrane channel current for ion species *k* as (Sterratt et al., 2011):

where ${g}_{L}^{k}$ is a constant leak conductivity, and *E*^{k} is the ion specific reversal potential, given by

with valence *z*^{k}, Faraday's constant *F*, absolute temperature *T*, and gas constant *R*.

##### 2.1.4.2. Hodgkin-Huxley membrane dynamics

In order to model active membrane dynamics, we use the standard Hodgkin-Huxley membrane model (Hodgkin and Huxley, 1952). The ion species under consideration are sodium Na^{+}, potassium K^{+}, and chloride Cl^{-}, and the model additionally introduces three gating variables *m, h, n* associated with sodium channel activation, potassium channel activation, and potassium channel inactivation, respectively. The membrane potential ϕ_{M} is then modeled by the following specialization of (11):

with ion specific membrane channel currents:

Here, ḡ^{k} is the maximal conductivity for ion species *k*. The gating variables are governed by the following ODE:

for *p* ∈ {*m, h, n*}. The rate constants α_{p} and β_{p} take the form

where *p*_{∞} is the steady state value for activation and τ_{p} is the time constant.

#### 2.1.5. Initial and Boundary Conditions

We assume that initial conditions are given for all ion concentrations, both intracellularly and extracellularly:

Furthermore, we assume that these conditions are compatible with the assumption of bulk electroneutrality, i.e., that the initial state of the system satisfies:

In addition, we assume that an initial condition is given for the membrane potential:

Finally, a set of boundary conditions will close the system. We describe specific boundary conditions in the numerical experiments in section 2.3.

#### 2.1.6. Summary of Governing Equations

In summary, the mathematical framework for electrodiffusion with explicit geometrical representation of the cell membranes is comprised of the bulk equations (1), (2) with (3), the interface conditions (7), (11) with (6) and (9), and (15) with (14), the initial conditions (24) and (26), and additional boundary conditions. We will refer to this set of equations as the KNP-EMI framework.

### 2.2. Numerical Methods

To solve the KNP-EMI framework numerically, we consider a finite difference time integration scheme, a splitting scheme, and a mortar finite element method in space. We derive the new finite element scheme and describe the splitting algorithm in the sections below.

#### 2.2.1. Weak Formulation of the Governing Equations

Multiplying (1) with test functions ${v}_{r}^{k}$ (for *r* ∈ {*i, e*}), integration over the intracellular and extracellular domains Ω_{i} and Ω_{e} separately, integration by parts, and inserting (15) for the ion fluxes across the membrane, yields

Similarly, multiplying (2) by test functions *w*_{r} for *r* ∈ {*i, e*}, integration by parts and inserting (7) for the total membrane current, yields

The zero average constraint (4) is enforced by introducing an additional unknown (a Lagrange multiplier) *c*_{e} ∈ ℝ along with a test function *d*_{e} ∈ *R* and letting

Finally, multiplying (11) by a test function *q*, and integrating over Γ yields

We remark that this is a weak formulation of a set of time-dependent, non-linear equations. In particular, recall that *I*_{ch} and ${\text{I}}_{\text{ch}}^{k}$ depend on ϕ_{M} and [*k*]_{r} cf. (9) while ${\alpha}_{r}^{k}$ depends on [*k*]_{r} cf. (14).

To solve this system numerically, we consider the following approximations.

• We discretize the time derivatives in (27)–(28) and (32) using a finite difference method.

• We approximate ${\text{J}}_{r}^{k}$ at time *t*^{n} by the linearized ion flux density [cf. (3)]:

• We evaluate ${\alpha}_{r}^{k}$ at time *t*^{n} by the previous value [cf. (14)]:

Moreover, we evaluate *I*_{ch} and the discretization of (32) depending on the choice of ion channel model (cf. section 2.1.4) as follows.

• For the passive model, we insert the linear relation given by (16) directly in (27)–(28) and (32). Moreover, the implicit discretization of (32) reads as:

at time *t*^{n} with ${\varphi}_{M}^{n}={\varphi}_{i}^{n}-{\varphi}_{e}^{n}$ and Δ*t* = *t*^{n} − *t*^{n−1}.

• For the Hodgkin-Huxley model, we use the following two-step splitting procedure. Consider *n* ∈ [1, …, *N*] with *t*^{n} − *t*^{n−1} = Δ*t*, and assume that ${\left[k\right]}_{r}^{n-1}$ and ${\varphi}_{M}^{n-1}$ at time step *t*^{n−1} are known.

- In the first (ODE) step, we update the membrane potential ${\varphi}_{M}^{n}$ at time step *t*^{n} by solving the ODE system (17)–(23), with *I*_{M} set to zero, using 25 explicit (forward) Euler steps of size Δ*t*^{*} = Δ*t*/25.

- In the second (PDE) step, we solve for ${\left[k\right]}_{r}^{n}$, ${\varphi}_{r}^{n}$ and ${I}_{M}^{n}$ (for *r* ∈ {*i, e*}) in the linear system arising from spatial discretization of (27)–(32), with *I*_{ch} set to zero in (32), and *I*_{ch} approximated by

in (27)–(28), where ${\varphi}_{M}^{n}$ is the membrane potential solution at *t*^{n} from the ODE step (see section 2.2.2 for details). The implicit discretization of (32) reads as:

where ${\varphi}_{M}^{n}={\varphi}_{i}^{n}-{\varphi}_{e}^{n}$ is the membrane potential solution at *t*^{n} from the ODE step.

The steps are repeated until global end time *t*^{N} is reached.

#### 2.2.2. Spatial Discretization

To numerically solve the PDE part of the governing equations defined on the domain Ω = Ω_{i}∪Ω_{e}, we use a mortar finite element method. We discretize each subdomain Ω_{r} by a conforming mesh ${{T}}_{r}$ for *r* ∈ {*i, e*}. We assume that the meshes ${{T}}_{i}$ and ${{T}}_{e}$ match at the common interface Γ, and define a (lower-dimensional) mesh ${{T}}_{\Gamma}$ of this interface (cf. Figure 2).

**Figure 2**. Schematic representation of meshes for the discretization of the PDE part of the governing electrodiffusive equations using a mortar finite element method. Mesh ${{T}}_{e}$ of the extracellular subdomain Ω_{e} **(left)**, mesh ${{T}}_{i}$ of the intracellular subdomain Ω_{i} **(middle)** and mesh ${{T}}_{\Gamma}$ of the interface Γ **(right)**. Note that the shared facets of the extracellular and intracellular meshes form the (codimension 1) mesh of the interface.

Next, we introduce separate finite element spaces for approximating the unknown fields in the weak formulation (27)–(32), [*k*]_{r} : Ω_{r} → ℝ, ϕ_{r} : Ω_{r} → ℝ for *r* ∈ {*i, e*}, *I*_{m} : Γ → ℝ. We approximate the ion concentrations [*k*]_{r} and potentials ϕ_{r} using continuous piecewise linear polynomials (linear Lagrange finite elements) over the meshes ${{T}}_{r}$. These fields thus have degrees of freedom defined on the vertices of the extracellular and intracellular meshes. The Lagrange multiplier *c*_{e} is approximated using a single real number. Furthermore, the transmembrane current *I*_{M} is approximated using continuous piecewise linear polynomials over the facet mesh ${{T}}_{\Gamma}$. We denote the finite element spaces for approximating [*k*]_{r} by ${V}_{r}^{k}$, the spaces for approximating ϕ_{r} by *W*_{r} and the spaces for approximating *I*_{M} by *Q*. Let ${\langle u,v\rangle}_{\Omega}={\int}_{\Omega}uv\text{d}x$. For notational simplicity, we denote the approximation of [*k*]_{r} by [*k*]_{r}, the approximation of ϕ_{r} by ϕ_{r}, and the approximation of *I*_{M} by *I*_{M} below. We here use linear polynomials for concreteness, but the formulation also applies directly for higher order polynomials.

We then solve the PDE step in the two-step splitting scheme described in section 2.2.1 as follows: given ${\left[k\right]}_{r}^{n}\in {V}_{r}^{k}$ and ${\varphi}_{M}^{n}\in Q$ at time step *t*^{n}, and the previously computed ${\text{I}}_{\text{ch}}^{k}$ and α^{k} [cf. (35) and (33)], find the ion concentrations ${\left[k\right]}_{r}\in {V}_{r}^{k}$, the potentials ϕ_{r} ∈ *W*_{r}, the total transmembrane current density *I*_{M} ∈ *Q* at time step *t*^{n+1} (and the Lagrange multiplier *c*_{e} ∈ ℝ) such that:

for all ${v}_{i}^{k}\in {V}_{i}^{k}$, ${v}_{e}^{k}\in {V}_{e}^{k}$, *w*_{i} ∈ *W*_{i}, *w*_{e} ∈ *W*_{e}, *d*_{e} ∈ ℝ and *q* ∈ *Q*. The ion flux terms on the right-hand side are replaced by appropriate boundary conditions in the subsequent sections.

To evaluate the accuracy of the numerical solutions defined over Ω_{r} for *r* ∈ {*i, e*}, we use the standard *L*^{2} and *H*^{1} norms denoted by ||·||_{0} and ||·||_{1}, respectively: for *u* : Ω_{r} → ℝ,

In addition, for *I* : Γ → ℝ, we define the broken *L*^{2}-norm by summing over the *L*^{2}-norms over the mesh cells of the interface mesh ${{T}}_{\Gamma}$:

#### 2.2.3. Implementation

The numerical scheme was implemented using a mixed dimensional framework from the FEniCS finite element library (Alnæs et al., 2015). The linear systems arising in the numerical experiments were solved using a direct (MUMPS) solver. The code is publicly available at: https://zenodo.org/record/3492075#.XahQOhh9g5k.

#### 2.2.4. Comparison With EMI Framework

In the numerical experiments comparing the KNP-EMI and the EMI models, the EMI model is discretized using the mortar finite element formulation as presented in Tveito et al. (2017b).

### 2.3. Computational Models and Parameters

We consider two model set-ups for testing the presented methodology (Model A and B), a model (Model C) for comparing simulation results between the KNP-EMI and EMI frameworks, and a model for studying ephaptic coupling (Model D). The model set-ups are described in detail here. The model parameters are given in Table 1, unless otherwise stated in the text. We assume that all axons in each simulation have the same membrane channel current *I*_{ch}. We denote the spatial coordinates in this and subsequent sections by (*x, y, z*).

#### 2.3.1. Model A: One Axon With a Passive Membrane Model

For Model A, we consider a two-dimensional domain $\Omega ={\Omega}_{i}\cup {\Omega}_{e}=\left[0,6.0\xb71{0}^{-5}\right]\times \left[0,6.0\xb71{0}^{-5}\right]$ m, with one intracellular domain (cell) ${\Omega}_{i}=\left[6.0\xb71{0}^{-6},5.6\xb71{0}^{-5}\right]\times \left[2.8\xb71{0}^{-5},3.4\xb71{0}^{-5}\right]$ m. We mesh this domain by dividing the domain into *n* × *m* rectangles, with Δ*x* = 6.0·10^{−5}/*n* and Δ*y* = 3.0·10^{−5}/*m*, and dividing each rectangle into two triangles by a diagonal, for a series of Δ*x* = Δ*y* = 2.0·10^{−6}, 1.0·10^{−6}, 5.0·10^{−7}, 2.5·10^{−7} m. We model *I*_{ch} using the passive model, as described in section 2.1.4, and prescribe a synaptic input *I*_{syn} of the form

where α is the synaptic time constant, *H*(*x*) = {1 for *x* ∈ *Z* and 0 elsewhere} for an interval *Z*. We let *Z* = [5.0 · 1.0^{−5}, 1.0 · 10^{−5}] m, and set *t*_{0} = 0, ${g}_{\text{syn}}=1.25\xb71{0}^{3}$ S/m^{2}. At the exterior boundary ∂Ω, we apply the boundary condition

describing that no ions can leave or enter the system.

#### 2.3.2. Model B: One Axon With a Passive Membrane Model and Non-physical Parameters

To evaluate the numerical accuracy of the mortar finite element scheme presented in section 2.2, we construct an analytical solution using the method of manufactured solutions (Roache, 1998). In particular, we let the analytical solution to (1)–(15) be given by:

with the passive model *I*_{ch} = ϕ_{M} and with *I*syn = 0. We assume that the parameter values all equal one: ${C}_{m}={D}_{i}^{k}={D}_{e}^{k}=F=G=R=1$, and that *K* = {Na^{+}, K^{+}, Cl^{−}}. We consider a two-dimensional domain Ω = Ω_{i}∪Ω_{e} = [0, 1] × [0, 1], with one intracellular domain Ω_{i} = [0.25, 0.75] × [0.25, 0.75]. The domain is meshed as for Model A (cf. section 2.3.1) for a series of *n* = *m* = 8, 16, 32, 64, 128, 256. In the numerical experiments for this test case, we initially let $\Delta t=\frac{1}{64}\xb71{0}^{-5}$, and then quarter the timestep in each series. The errors are evaluated at $t=\frac{2}{64}\xb71{0}^{-5}$.

#### 2.3.3. Model C: Multiple Axons With a Passive Membrane Model

For Model C, we define three different two-dimensional domains: (C1) a domain with one intracellular region (cell), (C2) a domain with two intracellular regions with a distance of 4.0·10^{−6}m in the *y*-direction between the cells, and (C3) a domain with two intracellular regions with a distance of 1.0·10^{−5}m in the *y*-direction between the cells. More precisely, we let

**Model C1:** $\Omega ={\Omega}_{i}\cup {\Omega}_{e}=\left[0,1.2\xb71{0}^{-4}\right]\times \left[0,1.2\xb71{0}^{-4}\right]$ m, ${\Omega}_{i}=\left[3.5\xb71{0}^{-5},8.5\xb71{0}^{-5}\right]\times \left[5.7\xb71{0}^{-5},6.3\xb71{0}^{-5}\right]$ m.

**Model C2:** $\Omega ={\Omega}_{i}^{1}\cup {\Omega}_{i}^{2}\cup {\Omega}_{e}=\left[0,1.2\xb71{0}^{-4}\right]\times \left[0,1.2\xb71{0}^{-4}\right]$ m, with two cells ${\Omega}_{i}^{1}=\left[3.5\xb71{0}^{-5},8.5\xb71{0}^{-5}\right]\times \left[5.2\xb71{0}^{-5},5.8\xb71{0}^{-5}\right]$ m and ${\Omega}_{i}^{2}=\left[3.5\xb71{0}^{-5},8.5\xb71{0}^{-5}\right]\times \left[6.2\xb71{0}^{-5},6.8\xb71{0}^{-5}\right]$ m.

**Model C3:** Ω as in Model C2 but with ${\Omega}_{i}^{1}=\left[3.5\xb71{0}^{-5},8.5\xb71{0}^{-5}\right]\times \left[4.9\xb71{0}^{-5},5.5\xb71{0}^{-5}\right]$ m and ${\Omega}_{i}^{2}=\left[3.5\xb71{0}^{-5},8.5\xb71{0}^{-5}\right]\times \left[6.5\xb71{0}^{-5},7.1\xb71{0}^{-5}\right]$ m.

The ion channel currents ${\text{I}}_{\text{ch}}^{k}$ are modeled using the passive model described in section 2.1.4. The synaptic input current model (37) is applied with ${g}_{\text{syn}}=1.25\xb71{0}^{3}$ S/m^{2}, *t*_{0} = 0, and with *Z* = [3.5·10^{−5}, 4.0·10^{−5}] m for Model C1, ${Z}_{1}=\left[6.0\xb71{0}^{-5},6.5\xb71{0}^{-5}\right]$ m for Model C2, and ${Z}_{2}=\left[5.5\xb71{0}^{-5},6.0\xb71{0}^{-5}\right]$ m for Model C3. At the exterior boundary ∂Ω, we apply the boundary condition (38). In order to compare the KNP-EMI and the EMI framework, we set the bulk conductivity σ_{r} in the EMI model by (5) with initial values ${\left[k\right]}_{r}^{0}$ for the ion concentration [*k*]_{r}. Note that σ_{r} will generally change over time.

#### 2.3.4. Model D: Axon Bundle With Active Hodgkin-Huxley Membrane Model

For Model D, we consider a domain $\Omega ={\Omega}_{i}^{1}\cup \cdots \cup {\Omega}_{i}^{9}\cup {\Omega}_{e}=\left[0,4.0\xb71{0}^{-4}\right]\xb7\left[0,1.4\xb71{0}^{-6}\right]\times \left[0,1.4\xb71{0}^{-6}\right]$ m, where 9 cuboidal cells of size 3.9·10^{−4} × 2.0·10^{−7} × 2.0·10^{−7} m are placed uniformly throughout Ω (cf. Figure 1). The distance between the cells is 1.0·10^{−7} m. The domain is meshed as in section 2.3.1 with Δ*y* = Δ*z* = 1.0·10^{−7} m and Δ*x* = 1.25·10^{−6} m. The ion channel currents are modeled using the Hodgkin-Huxley model as described in section 2.1.4. An action potential is induced every 20 ms throughout the simulations by applying the synaptic input current model (37) with *g*_{syn} = 40 S/m^{2}, α = 0.002 s, and *t*_{0} = 0, 0.02, 0.04 s. We ran two sets of simulations: (1) stimulating, i.e., applying the synaptic current to the membrane of, the middle axon only (axon A in Figure 1), and (2) stimulating the 8 axons around axon A (axons B,C in Figure 1). At the exterior boundary, we apply the boundary condition (38).

## 3. Results

We here present results from numerical experiments using the KNP-EMI framework and the numerical method presented above. We start by assessing the accuracy (Model A and B) and performance (Model A and D) of the numerical method. Next, we compare the KNP-EMI and EMI frameworks in idealized 2D axons (Model C), before we finally investigate ephaptic coupling in unmyelinated axons bundles (Model D). Note that the values in this section are given in physiologically reasonable units for the sake of readability; i.e., the results have been converted from the SI base units, e.g., seconds to milliseconds.

### 3.1. Numerical Verification and Accuracy

To evaluate the numerical accuracy and convergence of the proposed numerical approach, we consider two sets of experiments. First, we examine the convergence of the model under mesh refinement by visual inspection. Second, we perform a formal convergence analysis for a smooth test case with manufactured solution.

#### 3.1.1. Inspection of Convergence Under Mesh Refinement

The extracellular potential and sodium (Na^{+}) concentration of Model A for four different mesh resolutions are shown at *t* = 10 ms in Figure 3. The system quickly (after 3 ms) reaches a semi-steady state where the membrane potential does not change notably over time, but there is a slow exchange of sodium (Na^{+}) and potassium (K^{+}) ions due to the leak currents. The extracellular sodium concentration does not appear to change visibly under mesh refinement, and the extracellular potential seems to reach a converged state for the finest mesh resolution. More precisely, the mean relative difference (over time) between the solutions for the finest mesh (Δ*x* = 0.25) and the coarser mesh resolutions Δ*x* = 2.0, Δ*x* = 1.0, and Δ*x* = 0.5 are 4.8·10^{−5}, 4.6·10^{−6}, and 1.6·10^{−6} for [Na]_{e}, respectively, and 2.0·10^{−2}, 2.3·10^{−3}, and 4.3·10^{−4} for ϕ_{e}, respectively, at the point (4, 31). We conclude the differences between the solutions are small and decreasing, indicating convergence of the method.

**Figure 3**. Model A: Comparison (under mesh refinement) of extracellular potential **(A–D)** and extracellular sodium concentration **(E–H)** in the surroundings of a single simplified axon at *t* = 10 ms.

#### 3.1.2. Convergence Rates of Numerical Solutions

Using Model B, we analyzed the rates of convergence for the approximations of all solution variables under refinement in space and time. Based on properties of the approximation spaces and the time discretization, the optimal theoretical rate of convergence is 1 in the *H*^{1}-norm and 2 in the *L*^{2}-norm and the broken *L*^{2}-norm. Our numerical findings (Table 2) are in agreement with the theoretically optimal rates. We observe second order convergence in the *L*^{2}-norm for the approximation of the extracellular and intracellular concentrations and potentials, and first order convergence in the *H*^{1}-norm. For the transmembrane ionic current *I*_{M}, we observe a convergence rate of 1.5 in the broken *L*^{2}-norm. The loss of convergence of ~0.5 for *I*_{M} is likely due to a lack of smoothness of the interface in the test domain.

**Table 2**. Approximation errors (with convergence rates in parenthesis) for the extracellular and intracellular concentrations and potentials, and transmembrane current, under simultaneous refinement in time and space.

#### 3.1.3. Effect of Boundary Conditions

To examine whether the size of the extracellular space (ECS) affects the solution near the axon, we consider Model C1 with both the default size of the ECS (120 × 120 μm) and with an extended ECS (240 × 240 μm). The axon is placed in the same position in both cases. The two cases do not differ substantially in the extracellular Na^{+} concentrations or the extracellular potentials near the axon (Figure 4). The maximal difference between the default model and the model with extended extracellular space, measured 2μm above the axon at *t* = 10 ms, for ϕ_{e} and [Na]_{e} are 6.33·10^{−3} mV and 5.72·10^{−6} mM, respectively.

**Figure 4**. Comparison of results in Model C1 with the default size of ECS (120 × 120 μm), and with an extended ECS (240 × 240 μm) at *t* = 10 ms. The upper panels display the ECS potential, both at 2μm above the axon **(A)** and in the surrounding ECS for the default model **(B)** and for the model with extended ECS **(C)**. The lower panels display the ECS sodium concentration, both at 2μm above the axon **(D)** and in the surrounding ECS for the default model **(E)** and for the model with extended ECS **(F)**.

### 3.2. Numerical Performance

To evaluate the performance and scalability of the implementation of the presented framework, we consider an additional set of experiments measuring the memory usage and CPU timings for simulations of Model A (2D) and D (3D). We observe that the memory usage increases sublinearly with the system size: increasing the system size by a factor of four for Model A leads to an increase in memory of a factor 1−3 (Table 3). We observe that the CPU time for the simulations grows superlinearly with the system size: increasing the system size by a factor of four for Model A leads to an increase in total simulation time of a factor 3−5 (Table 3). This behavior is expected as the linear systems are solved using a direct solver. The total simulation time is dominated by the cost of finite element assembly and linear solves (70–94%). For small system sizes, the time required for finite element assembly is comparable to that of the linear solves. However, for larger system sizes, and especially in 3D, the linear solution time dominates the total simulation cost.

### 3.3. Comparison of the KNP-EMI and EMI Frameworks in Idealized Axons

The KNP-EMI framework extends the EMI framework by explicitly modeling the ionic concentrations and incorporating ionic electrodiffusion. A key question is when and to what extent the solutions from the two (KNP-EMI and EMI) frameworks differ. To compare the two frameworks, we consider three models (Model C1, C2, and C3) and compare the corresponding solution of the KNP-EMI equations (the KNP-EMI solution) with the solution of the EMI equations (the EMI solution).

We first consider the extracellular potential resulting from stimulating a single axon (Model C1) using the KNP-EMI and EMI frameworks (Figure 5). We observe that the KNP-EMI and EMI solutions are qualitatively very similar: an extracellular potential difference of ~0.12 mV along the length of the axon develops in both (Figures 5A–C).

**Figure 5**. A comparison of the KNP-EMI and the EMI frameworks using Model C1 at *t* = 10 ms. Extracellular potentials measured 2μm above the cell **(A)**. Extracellular potentials from the KNP-EMI **(B)** and the EMI framework **(C)** surrounding the cell.

Next, we compare the extracellular potentials resulting from stimulating two neighboring axons (Model C2 and C3) using the KNP-EMI and the EMI frameworks (Figure 6). The two models differ by the distance between the axons. For Model C2, we again observe that the KNP-EMI and EMI solutions match closely, but differ locally. The maximal difference between the extracellular potential solutions is 0.016 mV (Figures 6A–C). For Model C3, we observe the analogous behavior, but note that the extracellular field is weaker than for Model C2 (Figures 6D–F).

**Figure 6**. A comparison of the KNP-EMI and the EMI frameworks using Model C2 and C3 at *t* = 10 ms. The extracellular potentials from Model C2 **(A)** and C3 **(D)** on the midline between the neurons (*y* = 60μm). The extracellular potentials surrounding the cells as calculated by KNP-EMI **(B)** and EMI **(C)** using Model C2, and by KNP-EMI **(E)** and EMI **(F)** using Model C3.

### 3.4. Ephaptic Coupling in Unmyelinated Axon Bundles

We now turn to explore the effect of ephaptic coupling in an idealized axon bundle with 9 axons using the KNP-EMI framework. We consider two sets of simulations using Model D: (1) stimulating, i.e., applying the synaptic current to the membrane of, the middle axon only (axon A, Figure 1), and (2) stimulating the 8 axons around axon A (axons B,C, Figure 1).

#### 3.4.1. Electrodiffusion Effects in Unmyelinated Axon Bundles

To investigate ephaptic coupling, we first apply a synaptic current to stimulate the cell membrane of the middle axon of the axon bundle (axon A, Figure 1, Model D). The synaptic current induces a series of action potentials in axon A and also induces substantial changes in the surrounding extracellular potential (Figures 7A,B). The extracellular potential fluctuations further spread to axon B. However, the ephaptic effect on the membrane potential in axon B is relatively small (1–2 mV), and is not sufficient to reach the threshold for inducing an action potential (Figure 7C).

**Figure 7**. Effects of ephaptic coupling in a bundle of axons at *x* = 200 μm. The membrane potential (ϕ_{M}) of axon A **(A)**, and the extracellular potential (ϕ_{e}) measured 0.05 μm away from the membrane of axon A **(B)** during stimuli of axon A only. Ephaptic coupling measured in axon B when only axon A is stimulated **(C)**, and measured in axon A when all peripheral axons (B–C) are stimulated **(D)**. Setting σ_{i} = 1.0 S/m and σ_{e} = 0.1 S/m in the EMI framework increases the ephaptic coupling to the point where simultaneous action potentials in all eight surrounding axons will induce an action potential in the central axon **(F)**. However, only stimulating the middle axon (A) will not induce action potentials in the peripheral axons **(E)**.

The ephaptic effect is stronger if we simultaneously stimulate the cell membranes of all eight peripheral axons (axons B–C). Again, we observe a series of action potentials in the eight stimulated axons. Moreover, the combined ephaptic currents have a pronounced excitatory effect on axon A, but again fail to induce an action potential there (Figure 7D).

The difference between the EMI and KNP-EMI simulations are due to the time evolution of the intracellular and extracellular ion concentrations, accounted for by the KNP-EMI model but not by the EMI model. For each action potential fired, the Nernst potential will change due to alterations in the ionic concentrations using the KNP-EMI framework (Figure 8), whereas in the EMI framework the Nernst potential is constant.

**Figure 8**. Ion concentration dynamics in an axon bundle measured at the middle axon (A) at *x* = 200 μm using the KNP-EMI framework, both when middle axon (A) is stimulated, and when all peripheral axons (B–C) are stimulated. Extracellular sodium **(A)** and extracellular potassium **(B)** concentrations evaluated 0.05 μm away from axon A. Intracellular sodium **(C)** and intracellular potassium **(D)** concentrations evaluated at the center of axon A. Reversal potentials for sodium **(E)** and potassium **(F)** at the membrane of axon A.

Our predictions differ from those made in a similar study by Bokil et al. (2001), who found that a single active neighbor can induce action potentials in all nearby axons. We hypothesize that the main explanation for these differences is that the bulk conductivities differ between the two studies. Here, in the KNP-EMI framework, the bulk conductivities are functions of the ion concentrations [cf. (5)]. Using realistic values for the intra- and extracellular ion concentrations, we obtained bulk conductivities values of σ_{i}≈2.01μS/μm and σ_{e}≈1.31 S/m. In contrast, Bokil et al. set the bulk conductivities as free parameters, with σ_{i} = 1S/m and σ_{e} = 0.1S/m as the corresponding effective bulk conductivities in the EMI model. Tveito et al. (2017b) found that the ephaptic current was inversely proportional to σ_{e}, which suggests that the ephaptic current was more than seven times stronger in Bokil et al. (2001) than here.

In light of this, we repeated the simulations of the EMI model using the lower effective bulk conductivity values (σ_{i} = 1 S/m and σ_{e} = 0.1S/m). In this case, simultaneous stimulation of the 8 peripheral axons (B–C) induced an action potential in axon A (Figure 7F). Stimulation of axon A alone did not induce an action potential in the 8 peripheral axons (Figure 7E).

## 4. Discussion

We have presented a finite element-based numerical method for a revised mathematical model of ionic electrodiffusion with explicit geometrical representation of the extracellular space, the intracellular space and the cell membrane. Our numerical scheme is based on the mortar finite element method and is capable of efficiently handling complex geometries in one, two, or three spatial dimensions. Our numerical investigations demonstrate that the scheme is accurate and yields optimal convergence rates in the relevant norms.

Further, we compared the KNP-EMI framework and the EMI framework by computationally studying (i) extracellular fields surrounding passive idealized axons, and (ii) membrane potentials in a bundle of unmyelinated axons under Hodgkin-Huxley membrane mechanisms. The potentials predicted by the two frameworks are essentially identical during the first period (~5 ms) of the simulations, but the predictions later differ due to changes in ion concentrations (only accounted for by the KNP-EMI framework). We note that the strongest ephaptic coupling is due to changes in the Nernst potentials (ionic ephaptic coupling), and not via extracellular potentials (electric ephaptic coupling).

The predictions of ephaptic coupling made in this study differs from those made by Bokil et al. (2001) using cable theory. This discrepancy is likely due to differences in the extracellular bulk conductivities. Indeed, an important difference between geometrically explicit frameworks (e.g., PNP, EMI, and KNP-EMI) and homogenized frameworks (e.g., cable theory) is the interpretation of the bulk conductivities σ_{i} and σ_{e}. In homogenized frameworks based on volume-conductor theory, the bulk conductivity σ is interpreted as the tissue average, i.e., the effective bulk conductivity for currents propagating over distances in brain tissue (Holt and Koch, 1999; Pettersen and Einevoll, 2008; Reimann et al., 2013). Importantly, this tissue-averaged bulk conductivity is smaller than the actual conductivity of the extracellular solution, largely due to the fact that the extracellular space only constitutes about 20% of the total tissue volume. On the other hand, in the KNP-EMI framework, the bulk conductivities are defined in terms of the local ion concentrations and will thus vary consistently across the domain.

The KNP-EMI framework is comparable to the PNP framework (Lopreore et al., 2008; Pods et al., 2013; Holcman and Yuste, 2015; Cartailler et al., 2017a,b; Sacco et al., 2017) as both frameworks can account for the explicit morphology of neural tissue (Noguchi et al., 2005; Biess et al., 2007). The difference between the frameworks is in the way that the electrical potential ϕ is computed. In the more physically detailed PNP framework, ϕ is computed from the Poisson equation ∇^{2}ϕ = −ρ/ϵ, where ρ is the charge density and ϵ is the permittivity of the medium. In neural tissue, the charge relaxation time, i.e., the typical time-scale that ρ varies on, is in the order of 1 ns, and most of the local net charge density is resolved in nanometer thick layers surrounding neuronal membranes (Grodzinsky, 2011; Gratiy et al., 2017). Hence, simulations on the PNP framework requires a spatiotemporal resolution smaller than nanoseconds and nanometers, and become unstable otherwise. In contrast, the KNP-EMI framework circumvents the need for explicit modeling of charge accumulation near the membrane by assuming bulk electroneutrality, so that all net charge is associated as a membrane charge and not resolved spatially; that is, the membrane interface conditions (6)–(15) ensure that the mesh elements bordering the membrane contain the charges consistent with the membrane potential, regardless of mesh size. The electroneutrality condition has been shown to be a good approximation on spatiotemporal scales larger than micrometers and microseconds (Pods, 2017; Solbrå et al., 2018), and the application of it results in a more numerically stable framework for coarser time and space resolutions, allowing for longer simulations on larger domains. The differences between the PNP framework and (bulk-) electroneutral frameworks, such as KNP have been discussed extensively in previous works (Mori and Peskin, 2009; Pods, 2017; Solbrå et al., 2018).

An example of a phenomenon where large ion concentration changes in brain tissue build up over time, is (cortical) spreading depression. During spreading depression, the extracellular K^{+} concentration can change from a basal level of 3–5 mM to peak values at tens of mM over a period of several minutes (Somjen, 2001). As such, we advocate that the KNP-EMI model would be suitable for studying cellular level aspects of spreading depression computationally. However, for simulations of longer duration (>50 ms), the membrane mechanism model should be chosen carefully. The Hodgkin-Huxley formalism used in this paper to describe the membrane mechanisms does not account for the effect of ion pumps and co-transporters, which generally will strive to restore concentrations to baseline. As a consequence, the concentration changes taking place in our simulations are likely overestimates of what could be expected in a real biological system. Adding ion pumps and co-transporters to the membrane model would be relatively straightforward (Hübel and Dahlem, 2014; Hübel et al., 2014).

In conclusion, the KNP-EMI framework presented in this paper allows for detailed computational studies of the interplay between ion movement, membrane mechanisms and electrical potential in healthy neural tissue and under pathological conditions. The computational expense of KNP-EMI simulations compared to, e.g., homogenized models calls for further research into efficient and scalable solution methods.

## Data Availability Statement

The datasets analyzed for this study are publicly available and can be found at the following source: https://zenodo.org/record/3492075#.XahQOhh9g5k.

## Author Contributions

AE, AS, GE, GH, and MR designed the study, derived the mathematical model and numerical method, and wrote the manuscript. AE and AS developed the software implementation and performed the numerical experiments.

## Funding

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement 714892 (Waterscales), and from the Research Council of Norway (BIOTEK2021 Digital Life project DigiBrain, project 248828).

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

We thank Klas Pettersen (NORA) and Cécile Daversin-Catty (Simula Research Laboratory) for useful discussions.

## References

Agudelo-Toro, A. (2012). *Numerical simulations on the biophysical foundations of the neuronal extracellular space* (Ph.D. thesis), Niedersächsische Staats-und Universitätsbibliothek Göttingen, Göttingen, Germany.

Agudelo-Toro, A., and Neef, A. (2013). Computationally efficient simulation of electrical activity at cell membranes interacting with self-generated and externally imposed electric fields. *J. Neural Eng*. 10:026019. doi: 10.1088/1741-2560/10/2/026019

Alnæs, M. S., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., et al. (2015). The FEniCS project version 1.5. *Archiv. Numer. Softw*. 3, 9–23. doi: 10.11588/ans.2015.100.20553

Anastassiou, C. A., and Koch, C. (2015). Ephaptic coupling to endogenous electric field activity: why bother? *Curr. Opin. Neurobiol*. 31, 95–103. doi: 10.1016/j.conb.2014.09.002

Anastassiou, C. A., Perin, R., Markram, H., and Koch, C. (2011). Ephaptic coupling of cortical neurons. *Nat. Neurosci*. 14:217. doi: 10.1038/nn.2727

Ayata, C., and Lauritzen, M. (2015). Spreading depression, spreading depolarizations, and the cerebral vasculature. *Physiol. Rev*. 95, 953–993. doi: 10.1152/physrev.00027.2014

Bernardi, C., Maday, Y., and Patera, A. T. (1993). “Domain decomposition by the mortar element method,” in *Asymptotic and Numerical Methods for Partial Differential Equations With Critical Parameters. NATO ASI Series (Series C: Mathematical and Physical Sciences)*, Vol. 384, eds H. G. Kaper, M. Garbey, and G. W. Pieper (Dordrecht: Springer), 384.

Biess, A., Korkotian, E., and Holcman, D. (2007). Diffusion in a dendritic spine: the role of geometry. *Phys. Rev. E* 76:021922. doi: 10.1103/PhysRevE.76.021922

Bokil, H., Laaris, N., Blinder, K., Ennis, M., and Keller, A. (2001). Ephaptic interactions in the mammalian olfactory system. *J. Neurosci*. 21:RC173. doi: 10.1523/JNEUROSCI.21-20-j0004.2001

Cartailler, J., Schuss, Z., and Holcman, D. (2017a). Analysis of the Poisson-Nernst-Planck equation in a ball for modeling the voltage-current relation in neurobiological microdomains. *Phys. D Nonlin Phenom*. 339, 39–48. doi: 10.1016/j.physd.2016.09.001

Cartailler, J., Schuss, Z., and Holcman, D. (2017b). Electrostatics of non-neutral biological microdomains. *Sci. Rep*. 7:11269. doi: 10.1038/s41598-017-11590-6

Dietzel, I., Heinemann, U., and Lux, H. (1989). Relations between slow extracellular potential changes, glial potassium buffering, and electrolyte and cellular volume changes during neuronal hyperactivity in cat. *Glia* 2, 25–44. doi: 10.1002/glia.440020104

Doucette, J. R. (1984). The glial cells in the nerve fiber layer of the rat olfactory bulb. *Anat. Rec*. 210, 385–391. doi: 10.1002/ar.1092100214

Goldwyn, J. H., and Rinzel, J. (2016). Neuronal coupling by endogenous electric fields: cable theory and applications to coincidence detector neurons in the auditory brain stem. *J. Neurophysiol*. 115, 2033–2051. doi: 10.1152/jn.00780.2015

Gratiy, S. L., Halnes, G., Denman, D., Hawrylycz, M. J., Koch, C., Einevoll, G. T., et al. (2017). From Maxwell's equations to the theory of current-source density analysis. *Eur. J. Neurosci*. 45, 1013–1023. doi: 10.1111/ejn.13534

Griff, E. R., Greer, C. A., Margolis, F., Ennis, M., and Shipley, M. T. (2000). Ultrastructural characteristics and conduction velocity of olfactory receptor neuron axons in the olfactory marker protein-null mouse. *Brain Res*. 866, 227–236. doi: 10.1016/S0006-8993(00)02291-5

Grodzinsky, F. (2011). *Fields, Forces, and Flows in Biological Systems*. London; New York, NY: Garland Science; Taylor & Francis Group.

Halnes, G., Mäki-Marttunen, T., Keller, D., Pettersen, K. H., Andreassen, O. A., and Einevoll, G. T. (2016). Effect of ionic diffusion on extracellular potentials in neural tissue. *PLoS Comput. Biol*. 12:e1005193. doi: 10.1371/journal.pcbi.1005193

Halnes, G., Mäki-Marttunen, T., Pettersen, K. H., Andreassen, O. A., and Einevoll, G. T. (2017). Ion diffusion may introduce spurious current sources in current-source density (CSD) analysis. *J. Neurophysiol*. 118, 114–120. doi: 10.1152/jn.00976.2016

Halnes, G., Østby, I., Pettersen, K. H., Omholt, S. W., and Einevoll, G. T. (2013). Electrodiffusive model for astrocytic and neuronal ion concentration dynamics. *PLoS Comput. Biol*. 9:e1003386. doi: 10.1371/journal.pcbi.1003386

Halnes, G., Østby, I., Pettersen, K. H., Omholt, S. W., and Einevoll, G. T. (2015). “An electrodiffusive formalism for ion concentration dynamics in excitable cells and the extracellular space surrounding them,” in *Advances in Cognitive Neurodynamics (IV). Advances in Cognitive Neurodynamics*, ed H. Liljenström (Dordrecht: Springer), 353–360. doi: 10.1007/978-94-017-9548-7_50

Han, K.-S., Guo, C., Chen, C. H., Witter, L., Osorno, T., and Regehr, W. G. (2018). Ephaptic coupling promotes synchronous firing of cerebellar Purkinje cells. *Neuron* 100, 564–578. doi: 10.1016/j.neuron.2018.09.018

Hodgkin, A. L., and Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. *J. Physiol*. 117, 500–544. doi: 10.1113/jphysiol.1952.sp004764

Holcman, D., and Yuste, R. (2015). The new nanophysiology: regulation of ionic flow in neuronal subcompartments. *Nat. Rev. Neurosci*. 16:685. doi: 10.1038/nrn4022

Holt, G., and Koch, C. (1999). Electrical interactions via the extracellular potential near cell bodies. *J. Comput. Neurosci*. 6, 169–184. doi: 10.1023/A:1008832702585

Hübel, N., and Dahlem, M. A. (2014). Dynamics from seconds to hours in Hodgkin-Huxley model with time-dependent ion concentrations and buffer reservoirs. *PLoS Comput. Biol*. 10:e1003941. doi: 10.1371/journal.pcbi.1003941

Hübel, N., Schöll, E., and Dahlem, M. A. (2014). Bistable dynamics underlying excitability of ion homeostasis in neuron models. *PLoS Comput. Biol*. 10:e1003551. doi: 10.1371/journal.pcbi.1003551

Kager, H., Wadman, W. J., and Somjen, G. G. (2000). Simulated seizures and spreading depression in a neuron model incorporating interstitial space and ion concentrations. *J. Neurophysiol*. 84, 495–512. doi: 10.1152/jn.2000.84.1.495

Koch, C. (1999). *Biophysics of Computation: Information Processing in Single Neurons. 1st Edn*. New York, NY: Oxford University Press.

Krassowska, W., and Neu, J. C. (1994). Response of a single cell to an external electric field. *Biophys. J*. 66, 1768–1776. doi: 10.1016/S0006-3495(94)80971-3

Lopreore, C. L., Bartol, T. M., Coggan, J. S., Keller, D. X., Sosinsky, G. E., Ellisman, M. H., et al. (2008). Computational modeling of three-dimensional electrodiffusion in biological systems: application to the node of Ranvier. *Biophys. J*. 95, 2624–2635. doi: 10.1529/biophysj.108.132167

Lowe, G. (2003). Electrical signaling in the olfactory bulb. *Curr. Opin. Neurobiol*. 13, 476–481. doi: 10.1016/S0959-4388(03)00092-8

Markram, H., Muller, E., Ramaswamy, S., Reimann, M. W., Abdellah, M., Sanchez, C. A., et al. (2015). Reconstruction and simulation of neocortical microcircuitry. *Cell* 163, 456–492. doi: 10.1016/j.cell.2015.09.029

Mori, Y., Fishman, G. I., and Peskin, C. S. (2008). Ephaptic conduction in a cardiac strand model with 3D electrodiffusion. *Proc. Natl. Acad. Sci. U.S.A*. 105, 6463–6468. doi: 10.1073/pnas.0801089105

Mori, Y., and Peskin, C. (2009). A numerical method for cellular electrophysiology based on the electrodiffusion equations with internal boundary conditions at membranes. *Commun. Appl. Math. Comput. Sci*. 4, 85–134. doi: 10.2140/camcos.2009.4.85

Niederer, S. (2013). Regulation of ion gradients across myocardial ischemic border zones: a biophysical modelling analysis. *PLoS ONE* 8:e60323. doi: 10.1371/journal.pone.0060323

Noguchi, J., Matsuzaki, M., Ellis-Davies, G. C., and Kasai, H. (2005). Spine-neck geometry determines NMDA receptor-dependent Ca^{2+} signaling in dendrites. *Neuron* 46, 609–622. doi: 10.1016/j.neuron.2005.03.015

Øyehaug, L., Østby, I., Lloyd, C. M., Omholt, S. W., and Einevoll, G. T. (2012). Dependence of spontaneous neuronal firing and depolarisation block on astroglial membrane transport mechanisms. *J. Comput. Neurosci*. 32, 147–165. doi: 10.1007/s10827-011-0345-9

Pettersen, K. H., and Einevoll, G. T. (2008). Amplitude variability and extracellular low-pass filtering of neuronal spikes. *Biophys. J*. 94, 784–802. doi: 10.1529/biophysj.107.111179

Pods, J. (2017). A comparison of computational models for the extracellular potential of neurons. *J. Integr. Neurosci*. 16, 19–32. doi: 10.3233/JIN-170009

Pods, J., Schönke, J., and Bastian, P. (2013). Electrodiffusion models of neurons and extracellular space using the Poisson-Nernst-Planck equations-numerical simulation of the intra- and extracellular potential for an axon model. *Biophys. J*. 105, 242–254. doi: 10.1016/j.bpj.2013.05.041

Rall, W. (1977). “Core conductor theory and cable properties of neurons,” in *Handbook of Physiology*, Chapter 3, eds E. Kandel, J. Brookhardt, and V. M. Mountcastle (Bethesda, MD: American Physiological Society), 39–97.

Reimann, M. W., Anastassiou, C. A., Perin, R., Hill, S. L., Markram, H., and Koch, C. (2013). A biophysically detailed model of neocortical local field potentials predicts the critical role of active membrane currents. *Neuron* 79, 375–390. doi: 10.1016/j.neuron.2013.05.023

Roache, P. J. (1998). *Verification and Validation in Computational Science and Engineering*, Vol. 895. Albuquerque, NM: Hermosa.

Sætra, M. J., Einevoll, G. T., and Halnes, G. (2020). An electrodiffusive, ion conserving Pinsky-Rinzel model with homeostatic mechanisms. *bioRxiv*. doi: 10.1101/2020.01.20.912378

Sacco, R., Airoldi, P., Mauri, A. G., and Jerome, J. W. (2017). Three-dimensional simulation of biological ion channels under mechanical, thermal and fluid forces. *Appl. Math. Model*. 43, 221–251. doi: 10.1016/j.apm.2016.10.053

Savtchenko, L. P., Poo, M. M., and Rusakov, D. A. (2017). Electrodiffusion phenomena in neuroscience: a neglected companion. *Nat. Rev. Neurosci*. 18:598. doi: 10.1038/nrn.2017.101

Shifman, A. R., and Lewis, J. E. (2019). Elfenn: a generalized platform for modeling ephaptic coupling in spiking neuron models. *Front. Neuroinform*. 13:35. doi: 10.3389/fninf.2019.00035

Solbrå, A., Wigdahl, B. A., van den Brink, Jonas Anders, M.-S. T E. G, and Geir, H. (2018). A Kirchhoff-Nernst-Planck framework for modeling large scale extracellular electrodiffusion surrounding morphologically detailed neurons. *PLoS Comput. Biol*. 14:e261107. doi: 10.1101/261107

Somjen, G. G. (2001). Mechanisms of spreading depression and hypoxic spreading depression-like depolarization. *Physiol. Rev*. 81, 1065–1096. doi: 10.1152/physrev.2001.81.3.1065

Sterratt, D., Graham, B., Gillies, A., and Willshaw, D. (2011). *Principles of Computational Modelling in Neuroscience*. Cambridge: Cambridge University Press. doi: 10.1017/CBO9780511975899

Syková, E., and Nicholson, C. (2008). Diffusion in brain extracellular space. *Physiol. Rev*. 88, 1277–1340. doi: 10.1152/physrev.00027.2007

Tveito, A., Jæger, K. H., Kuchta, M., Mardal, K.-A., and Rognes, M. E. (2017a). A cell-based framework for numerical modeling of electrical conduction in cardiac tissue. *Front. Phys*. 5:48. doi: 10.3389/fphy.2017.00048

Tveito, A., Jæger, K. H., Lines, G. T., Paszkowski, Ł., Sundnes, J., Edwards, A. G., et al. (2017b). An evaluation of the accuracy of classical models for computing the membrane potential and extracellular potential for neurons. *Front. Comput. Neurosci*. 11:27. doi: 10.3389/fncom.2017.00027

Wei, Y., Ullah, G., and Schiff, S. J. (2014). Unification of neuronal spikes, seizures, and spreading depression. *J. Neurosci*. 34, 11733–11743. doi: 10.1523/JNEUROSCI.0516-14.2014

Keywords: finite element, electrodiffusion, ion concentrations, cell membrane, ephaptic coupling, KNP-EMI

Citation: Ellingsrud AJ, Solbrå A, Einevoll GT, Halnes G and Rognes ME (2020) Finite Element Simulation of Ionic Electrodiffusion in Cellular Geometries. *Front. Neuroinform.* 14:11. doi: 10.3389/fninf.2020.00011

Received: 22 October 2019; Accepted: 05 March 2020;

Published: 25 March 2020.

Edited by:

Arnd Roth, University College London, United KingdomReviewed by:

Stavros I. Dimitriadis, Cardiff University, United KingdomGillian Queisser, Temple University, United States

Adam John Hunter Newton, Yale University, United States

Copyright © 2020 Ellingsrud, Solbrå, Einevoll, Halnes and Rognes. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Marie E. Rognes, meg@simula.no

^{†}These authors have contributed equally to this work