Finite Element Simulation of Ionic Electrodiffusion in Cellular Geometries

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 {ℝ}^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ⁿ, 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).

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³) 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 {ℝ}^d$ is the regional ion flux density [mol/(m²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²/s) of ion species k in the region r. The constant ψ = RTF⁻¹ 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²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 Icap:

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 Isyn, leaky passive neuron, Hodgkin-Huxley etc., and will be detailed further below in section 2.1.4. On the other hand, the capacitive current Icap 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·{\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 Icap, r^k in region Ω_r at the interface Γ be some fraction ${\alpha }_r^k$ of the total capacitive current Icap:

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ⁿ by the linearized ion flux density [cf. (3)]:

• We evaluate ${\alpha }_r^k$ at time tⁿ 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ⁿ with ${\varphi }_M^n={\varphi }_i^n-{\varphi }_e^n$ and Δt = tⁿ − tⁿ⁻¹.

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

- In the first (ODE) step, we update the membrane potential ${\varphi }_M^n$ at time step tⁿ 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ⁿ 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ⁿ 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).

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ⁿ, 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ⁿ⁺¹ (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² and H¹ norms denoted by ||·||₀ and ||·||₁, respectively: for u : Ω_r → ℝ,

In addition, for I : Γ → ℝ, we define the broken L²-norm by summing over the L²-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·10^{-5}\right]\times \left[0,6.0·10^{-5}\right]$ m, with one intracellular domain (cell) ${\Omega }_i=\left[6.0·10^{-6},5.6·10^{-5}\right]\times \left[2.8·10^{-5},3.4·10^{-5}\right]$ m. We mesh this domain by dividing the domain into n × m rectangles, with Δx = 6.0·10⁻⁵/n and Δy = 3.0·10⁻⁵/m, and dividing each rectangle into two triangles by a diagonal, for a series of Δx = Δy = 2.0·10⁻⁶, 1.0·10⁻⁶, 5.0·10⁻⁷, 2.5·10⁻⁷ 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⁻⁵, 1.0 · 10⁻⁵] m, and set t₀ = 0, $g_{\text{syn}}=1.25·10^3$ S/m². 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 Isyn = 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}·10^{-5}$, and then quarter the timestep in each series. The errors are evaluated at $t=\frac{2}{64}·10^{-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⁻⁶m in the y-direction between the cells, and (C3) a domain with two intracellular regions with a distance of 1.0·10⁻⁵m in the y-direction between the cells. More precisely, we let

Model C1: $\Omega ={\Omega }_i\cup {\Omega }_e=\left[0,1.2·10^{-4}\right]\times \left[0,1.2·10^{-4}\right]$ m, ${\Omega }_i=\left[3.5·10^{-5},8.5·10^{-5}\right]\times \left[5.7·10^{-5},6.3·10^{-5}\right]$ m.

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

Model C3: Ω as in Model C2 but with ${\Omega }_i^1=\left[3.5·10^{-5},8.5·10^{-5}\right]\times \left[4.9·10^{-5},5.5·10^{-5}\right]$ m and ${\Omega }_i^2=\left[3.5·10^{-5},8.5·10^{-5}\right]\times \left[6.5·10^{-5},7.1·10^{-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·10^3$ S/m², t₀ = 0, and with Z = [3.5·10⁻⁵, 4.0·10⁻⁵] m for Model C1, $Z_1=\left[6.0·10^{-5},6.5·10^{-5}\right]$ m for Model C2, and $Z_2=\left[5.5·10^{-5},6.0·10^{-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·10^{-4}\right]·\left[0,1.4·10^{-6}\right]\times \left[0,1.4·10^{-6}\right]$ m, where 9 cuboidal cells of size 3.9·10⁻⁴ × 2.0·10⁻⁷ × 2.0·10⁻⁷ m are placed uniformly throughout Ω (cf. Figure 1). The distance between the cells is 1.0·10⁻⁷ m. The domain is meshed as in section 2.3.1 with Δy = Δz = 1.0·10⁻⁷ m and Δx = 1.25·10⁻⁶ 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², α = 0.002 s, and t₀ = 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⁻⁵, 4.6·10⁻⁶, and 1.6·10⁻⁶ for [Na]_e, respectively, and 2.0·10⁻², 2.3·10⁻³, and 4.3·10⁻⁴ for ϕ_e, respectively, at the point (4, 31). We conclude the differences between the solutions are small and decreasing, indicating convergence of the method.

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¹-norm and 2 in the L²-norm and the broken L²-norm. Our numerical findings (Table 2) are in agreement with the theoretically optimal rates. We observe second order convergence in the L²-norm for the approximation of the extracellular and intracellular concentrations and potentials, and first order convergence in the H¹-norm. For the transmembrane ionic current I_M, we observe a convergence rate of 1.5 in the broken L²-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.

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⁻³ mV and 5.72·10⁻⁶ mM, respectively.

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).

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).

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).

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.

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 ∇²ϕ = −ρ/ϵ, 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