Finite Element Simulation of Ionic Electrodiffusion in Cellular Geometries

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.


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 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., 2013Halnes et al., , 2016Halnes et al., , 2017Niederer, 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 threedimensional 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(Halnes et al., , 2015Saetra et al., 2020), or components within a hybrid modeling scheme to compute extracellular dynamics (Halnes et al., 2016Solbrå 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.

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

Representation of the Computational Domain
We consider N domains i n ⊂ 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 = i 1 ∪ · · · ∪ i N ∪ 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 i n , by Ŵ n . We assume that Ŵ n ∩ Ŵ m = ∅ for all n = m and that Ŵ n ∩ ∂ = ∅ (It follows that ∂ i n ∩ ∂ e = ∅ for all n = 1, . . . , N.). For simplicity and clarity, we present the mathematical model for one intracellular region i 1 = i with membrane Ŵ below. The extension to multiple intracellular regions is immediate (but notationally cumbersome).

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] → R (mol/m 3 ) and the electrical potentials φ r : r × (0, T] → R (V). Conservation of ions for the bulk of each region r stipulates that for t ∈ (0, T]. Here, J k r : r × (0, T] → 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 J k r 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 k r 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 k r ∇[k] r and (ii) the ion concentrations that are transported via electrical potential gradients, i.e., the ion migration −D k r z k ψ −1 [k] r ∇φ r where D k r ψ −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., e φ e dx = 0.
(4) 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 k r .
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.

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) → R [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 I k ch : The channel currents I k ch are subject to modeling. Typical models for I k ch 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 J k r · 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 k cap,r in region r at the interface Ŵ be some fraction α k r of the total capacitive current I cap : Specifically, we assume that: and note that k∈K α k r = 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.

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 I k ch . 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).

Passive membrane dynamics
We model the passive membrane channel current for ion species k as (Sterratt et al., 2011): where g k L 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.

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.

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.

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.

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.

Weak Formulation of the Governing Equations
Multiplying (1) with test functions v k r (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 ∈ R along with a test function d e ∈ R and letting e φ e d e dx = 0.
Finally, multiplying (11) by a test function q, and integrating over Ŵ yields We remark that this is a weak formulation of a set of timedependent, non-linear equations. In particular, recall that I ch and I k ch depend on φ M and [k] r cf. (9) while α k r 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) • We evaluate α k r 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 φ n M = φ n i − φ n e and t = t n − t n−1 . • For the Hodgkin-Huxley model, we use the following twostep splitting procedure. Consider n ∈ [1, . . . , N] with t n − t n−1 = t, and assume that [k] n−1 r and φ n−1 M at time step t n−1 are known.
-In the first (ODE) step, we update the membrane potential φ n M 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 [k] n r , φ n r and I n M (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 φ n M 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 φ n M = φ n i − φ n e is the membrane potential solution at t n from the ODE step.
The steps are repeated until global end time t N is reached.

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 (lowerdimensional) mesh T Ŵ of this interface (cf. Figure 2).
Next, we introduce separate finite element spaces for approximating the unknown fields in the weak formulation (27) 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 Ŵ . We denote the finite element spaces for approximating [k] r by V k r , the spaces for approximating φ r by W r and the spaces for approximating I M by Q. Let u, v = uv dx. 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 [k] n r ∈ V k r and φ n M ∈ Q at time step t n , and the previously computed I k ch and α k [cf. (35) and (33)], find the ion concentrations [k] r ∈ V k r , 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 ∈ R) such that: for all v k i ∈ V k i , v k e ∈ V k e , w i ∈ W i , w e ∈ W e , d e ∈ R 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 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 Ŵ of the interface Ŵ (right). Note that the shared facets of the extracellular and intracellular meshes form the (codimension 1) mesh of the interface. denoted by · 0 and · 1 , respectively: for u : r → R, In addition, for I : Ŵ → R, we define the broken L 2 -norm by summing over the L 2 -norms over the mesh cells of the interface mesh T Ŵ :

Implementation
The numerical scheme was implemented using a mixed dimensional framework from the FEniCS finite element library (Alnaes 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.

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

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

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: [Na] e i = 0.7 + 0.3 sin(2πx) sin (2πy)   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 t = 1 64 · 10 −5 , and then quarter the timestep in each series. The errors are evaluated at t = 2 64 · 10 −5 .

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.

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.   Approximation errors are measured at t = 2 64 · 10 −5 , i.e., e.g., [Na]r − [Na] r,h 0 = [Na]r (t) − [Na] r,h (t) 0 and similarly for all reported values. The largest time step is t = 1/64 · 10 −5 ; the time step was quartered in each refinement level.

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.

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.

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.

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.

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

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

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 = 1 S/m and σ e = 0.1 S/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.1 S/m). In this case, simultaneous stimulation of the 8 peripheral axons (B-C) induced an action potential in axon A x: mesh resolution (for Model A), Size, linear system size (number of degrees of freedom); Memory, Maximal memory usage of simulation relative to baseline. T A , CPU time for finite element assembly; T S , CPU time for linear system solution; T, Total CPU time for simulation; N, number of timesteps.
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.
( Figure 7F). Stimulation of axon A alone did not induce an action potential in the 8 peripheral axons (Figure 7E).

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