An Evaluation of the Accuracy of Classical Models for Computing the Membrane Potential and Extracellular Potential for Neurons

1. Introduction

Computational modeling in neurophysiology is a rapidly developing field taking on problems of enormous complexity. This is illustrated in the recent paper by Markram et al. (2015) where the authors present results of amazingly detailed digital algorithmic reconstruction of a neocortical volume segment (about 0.29 mm³) of rat cortex, containing ~31,000 neurons with ~37 million synapses. The complexity of the project is astonishing, and it opens amazing perspectives for insight in the complexities of the brain. The paper also raises questions of more philosophical nature brilliantly examined in the accompanying perspective by Koch and Buice (2015).

The development of enormously complex computational models motivates closer examination of the basis of the mathematical models underpinning the computations. It is the purpose of this study to evaluate the accuracy of two basic models extensively used throughout the field of computational neurophysiology, and our main question is whether the popularity of these models is warranted by their accuracy.

The first model we consider is the celebrated Cable equation used to compute membrane potentials and transmembrane currents. This model is absolutely essential in computational neurophysiology, and is used in numerous papers every year. The derivation of the model is classical and can be found in any introduction to computational neurophysiology; see e.g., Sterratt et al. (2011), Ermentrout and Terman (2010), Scott (2002), Dayan and Abbott (2001), and Koch (1999). An important assumption in the most common derivation of the Cable equation is that the extracellular conductivity is very large, and that consequently the extracellular potential can be assumed to be constant. This assumption represents a major simplification of the model since the extracellular field does not have to be represented in the model, which means that a costly solution of a Poisson equation in the extracellular domain is avoided.

One way of interpreting the effect of ignoring the coupling to the extracellular potential is that (as we shall see below) we disregard the so-called ephaptic current; see e.g., Holt and Koch (1999). It is well known that neglecting this current represents a key assumption, and the validity of the assumption, and also the effect of ephaptic coupling, have previously been discussed by several authors; see e.g., Buzsáki et al. (2012), Bhalla (2012), Goldwyn and Rinzel (2016), Anastassiou et al. (2011), Anastassiou and Koch (2015), and Bokil et al. (2001). An analytical treatment of the effect of ephaptic currents on nerve pulses in parallel nerve fibers is given in Chapter 8 of Scott (2002). That exposition is motived by classical experiments performed by Katz and Schmitt (see e.g., Katz and Schmitt, 1940) and analyzed by an extension of the scalar Cable equation to a 2 × 2 system of partial differential equations governing the membrane potential of the neighboring nerve fibers. This work is followed up by Shneider and Pekker (2015) who suggest that the ephaptic current acts as a synchronization mechanism for action potentials along bundles of neurons. For axons, the coupling is particularly important in the unmyelinated case; see Bokil et al. (2001) for an analysis of bundles of olfactory nerve axons. Furthermore, Goldwyn and Rinzel (2016) recently studied ephaptic interactions in a bundle of neurons and found that the effects of the ephaptic currents were small but not negligible.

The question of ephaptic coupling between cells has been studied for a long time; 75 years ago (Arvanitaki, 1942) stated that there is no doubt that the activity of an element in the midst of a cell agglomeration can influence that of its neighbors, even when specialized contact surfaces for transmission, i.e., those loci traditionally known as synapses and which have been endowed with particular properties are lacking. In Holt and Koch (1999), Holt and Koch analyse the magnitude and possible consequences of ephaptic coupling. They observe that spikes from a neuron can cause an extracellular potential of a few mV near the cell body, and they analyse the effect of this on nearby cells. The impact of ephaptic coupling remains uncertain (Anastassiou et al., 2011; Anastassiou and Koch, 2015), but it seems to be acknowledged that ephaptic currents may be significant. However, it is usually not taken into account in most computational analyses of neurons, and the reason for this is clearly to improve computational efficiency. In this paper we will quantify the error introduced by this assumption. We will compare the results of the Cable equation to those of an accurate mathematical model which includes the ephaptic current. The more accurate model will be referred to as the EMI model since it builds on detailed representation of both the Extracellular space surrounding the neuron, the Membrane of the neuron and the Intracellular space of the neuron. EMI computations are typically much more CPU demanding than solving the Cable equation, but the model faithfully represents the physics of the neuron and its surroundings. Variants of the EMI model have been studied previously by e.g., Krassowska and Neu (1994), Ying and Henriquez (2007), Henríquez et al. (2013), Agudelo-Toro and Neef (2013), and Agudelo-Toro (2012). For linear membrane currents and specialized geometries, analytical solutions are available; see e.g., Rall (1962), Rall (1969), Klee and Rall (1977), Krassowska and Neu (1994), Ying and Henriquez (2007), and Agudelo-Toro and Neef (2013).

The second model we consider is the standard formalism for computing the extracellular potential based on solutions of the Cable equation. It is well known that if the current sources are given by Dirac delta functions, the solution of the Poisson equation, defined on an infinite domain, can be computed by an explicit formula, see e.g., Einevoll et al. (2013). Based on the solution of the Cable equation, the current sources can be defined for each compartment in the numerical solution, and the solution of the Poisson equation can (due to linearity) be given as the sum of contributions from all compartments. Note that in practice the so-called line-source approximation (Holt and Koch, 1999) where the current sources are assumed to be evenly distributed along cylindrical axes of dendritic compartments, is commonly used rather than the point-source approximation built on solutions of the Dirac delta functions. However, these two methods are directly related as a line source can be arbitrarily accurately approximated by a line of delta-function sources.

The combined use of the two models (a) and (b) in computing extracellular potentials is especially intriguing since (a) is solved based on the assumption that the extracellular field is constant, and then (b) is used to compute the non-constant extracellular potential.

We have evaluated the accuracy of these two basic models by comparing the results with the results obtained by solving the EMI model. Our findings can be summarized as follows:

1. We find that the membrane potential computed by the Cable equation qualitatively resembles the solution of the EMI model but may differ quantitatively (several millivolts) from the solution of the EMI model.

2. We find, using reasonable parameters, that the magnitude of the ephaptic current is comparable to the other currents in our example model, so that its omission is, in general, difficult to justify.

3. For our example application the error in neglecting the ephaptic effect when computing the extracellular potentials is found to be 10% or more, and stem from the inaccurate computation of the transmembrane currents when the extracellular potentials are assumed to be constant.

We have found the EMI model to be a useful framework for assessing the accuracy of the classical models. The EMI model is, however, much more computationally demanding, it is much more difficult to implement correctly, and therefore very challenging to apply to problems of greater complexity than the simple problems addressed in the present report.

The rest of this report is organized as follows: In the Methods Section, we derive the classical Cable equation and highlight what assumptions are needed to remove the extracellular potential from the model. Given the solution of the Cable equation, we show how to compute the extracellular potential by solving a boundary value problem, how to approximate the solution by solving a Poisson equation, and how to approximate the solution of the extracellular potential using a classical summation formula. Finally, we introduce the EMI model where the dynamics of the extracellular space, the cell membrane and the intracellular space are fully coupled, and we show how the EMI model can be solved numerically. In the Results Section, we study the error introduced in the model by ignoring the ephaptic currents and how the ephaptic current depends on the extracellular conductivity. Furthermore, we compare the extracellular potential around a single simplified neuron computed by various approximate models, and we also compare the extracellular potential between two simplified neurons. Finally, we show that the numerical solutions seem to converge under mesh refinement and that infinite domains can be reasonably well represented using large extracellular domains. Implications and relevance of the results are examined in the Discussion Section. In an Appendix in Supplementary Material we give a theoretical estimate of the error introduced by removing the ephaptic current.

2. Methods

The standard way of computing the extracellular potential surrounding a neuron is a two-step process: (a) solve the Cable equation, and (b) use the transmembrane currents defined by step (a) to compute the extracellular potential. Our aim is to assess the accuracy of the solution of these two steps. For comparison we will use an accurate model combining the Extracellular domain, the Membrane, and the Intracellular domain, referred to as the EMI model. Below, the EMI solution will be regarded as the reference solution, and therefore solutions computed by all other methods (derived below) will be compared to the EMI solution.

We will take care to try to explain exactly how the EMI model and the two-step models are defined and solved although the derivations presented here can, at least in part, be found elsewhere. The derivations will also help us clarify what assumptions underlie the various models.

2.1. The Classical Two-Step Method

We start by describing the two steps of the classical approach of computing the extracellular potential (Holt and Koch, 1999; Lindén et al., 2014). The first step is to compute the membrane potential and transmembrane currents. In the classical approach this problem is solved assuming a constant extracellular potential. We briefly review the derivation of the Cable equation in order to clarify exactly what assumption is made in order to remove the extracellular potential from the equation defining the membrane potential. By identifying what term is ignored in the equation, this term can be evaluated and used to illuminate the accuracy of the Cable equation.

The second step is to compute the extracellular potential by using the solution of the Cable equation to define the transmembrane current sources. This step can be done in numerous ways, and we will derive alternative methods starting with the approach considered to most faithfully represent the physics involved, and then derive simpler and more efficient methods in order to end up with the classical summation formula defining the extracellular potential.

2.1.1. The Cable Equation

Consider a simplified neuron geometry illustrated in Figure 1. The intracellular space of the neuron is denoted by Ω_i and the boundary of Ω_i is the membrane of the neuron, denoted by Γ. The size of Ω_i is given by l_x, l_y, and l_z. In the derivation of the Cable equation, the neuron is divided into compartments, see e.g., Sterratt et al. (2011), Scott (2002), and Ermentrout and Terman (2010), and it is assumed that the variations in the y- and z-directions are small and can be ignored. Our derivation is based on the version of the Cable equation used in Holt and Koch (1999). The compartments are denoted by Ω_i,k, and have length Δx. For the k-th compartment, the transmembrane current density (positive outward) is given by (see e.g., Sterratt et al., 2011)

$\begin{array}{ll}{I}_{\mathrm{\text{m}}}^k=C_{\mathrm{\text{m}}}\frac{d{v}^k}{dt}+I_{\mathrm{\text{ion}}}^k,& (1)\end{array}$

where v is the membrane potential, C_m is the cell membrane capacitance and I_ion is the ionic current density out of the cell. Furthermore, assuming ohmic resistance along the length of the neuron, we have

$\begin{array}{ll}\Delta x{I}_{k+1/2}={\sigma }_{\mathrm{\text{i}}}(u_{\mathrm{\text{i}}}^k-u_{\mathrm{\text{i}}}^{k+1}),& (2)\end{array}$

where $u_{\mathrm{\text{i}}}^k$ is the intracellular potential in compartment k, σ_i is the intracellular conductivity, and I_{k + 1/2} is the intracellular current density from compartment k to compartment k + 1. Applying Kirchhoff's current law, the sum of the currents flowing out of a compartment must equal the sum of currents flowing into a compartment, i.e.,

$\begin{array}{ll}\left|{\Gamma }^k\right|I_{\mathrm{\text{m}}}^k=l_y{l}_z(I_{k-1/2}-I_{k+1/2}),& (3)\end{array}$

where |Γ^k| is the membrane area associated with Ω_i,k. Therefore,

$\begin{array}{ll}\left|{\Gamma }^k\right|\left(C_{\mathrm{\text{m}}}\frac{d{v}^k}{dt}+I_{\mathrm{\text{ion}}}^k\right)=\frac{{\sigma }_{\mathrm{\text{i}}}{l}_y{l}_z}{\Delta x}(u_{\mathrm{\text{i}}}^{k-1}-2u_{\mathrm{\text{i}}}^k+u_{\mathrm{\text{i}}}^{k+1}).& (4)\end{array}$

To simplify notations, we assume that l_y = l_z : = h, and we have

$\begin{array}{ll}{C}_{\mathrm{\text{m}}}\frac{d{v}^k}{dt}+I_{\mathrm{\text{ion}}}^k=\frac{{\sigma }_{\mathrm{\text{i}}}}{4}\frac{h}{\Delta x^2}\left(u_{\mathrm{\text{i}}}^{k-1}-2u_{\mathrm{\text{i}}}^k+u_{\mathrm{\text{i}}}^{k+1}\right).& (5)\end{array}$

Certainly, in the limit of small compartments (Δx → 0), we have

$\begin{array}{ll}{C}_{\mathrm{\text{m}}}\frac{\partial v}{\partial t}+I_{\mathrm{\text{ion}}}=\eta \frac{{\partial }^2{u}_{\mathrm{\text{i}}}}{\partial x^2},& (6)\end{array}$

where we have introduced the conductance

$\begin{array}{ll}\eta =\frac{h{\sigma }_{\mathrm{\text{i}}}}{4}.& (7)\end{array}$

The membrane potential is defined as

$\begin{array}{ll}v=u_{\mathrm{\text{i}}}-u_{\mathrm{\text{e}}},& (8)\end{array}$

where u_e denotes the extracellular potential. Therefore, we can replace the intracellular potential u_i in Equation (6) by v + u_e to get

$\begin{array}{ll}{C}_{\mathrm{\text{m}}}\frac{\partial v}{\partial t}+I_{\mathrm{\text{ion}}}=\eta \left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^2{u}_{\mathrm{\text{e}}}}{\partial x^2}\right).& (9)\end{array}$

At this point a major assumption is introduced; it is assumed that the extracellular potential varies so little that it can be taken to be a constant (see e.g., Sterratt et al., 2011)¹;

$\begin{array}{ll}{u}_{\mathrm{\text{e}}}\approx \mathrm{\text{const}}.& (10)\end{array}$

Building on this assumption we arrive at the classical Cable equation

$\begin{array}{ll}{C}_{\mathrm{\text{m}}}\frac{\partial v}{\partial t}+I_{\mathrm{\text{ion}}}=\eta \frac{{\partial }^{2}v}{\partial x^2}.& (11)\end{array}$

Note that the term we ignored in the derivation of the Cable equation is

$\begin{array}{ll}{I}_{\mathrm{\text{eph}}}=\eta \frac{{\partial }^2{u}_{\mathrm{\text{e}}}}{\partial x^2},& (12)\end{array}$

which is referred to as the ephaptic current density (Holt and Koch, 1999). In the computations below we will compute this current using the EMI model and use it to quantify the effect of the assumption underlying the classical Cable equation.

FIGURE 1

Figure 1. Sketch of a simplified neuron of rectangular cuboid shape with dimensions l_x, l_y, and l_z. The intracellular domain is denoted Ω_i, the boundary is Γ, and the compartments of length Δx are denoted by Ω_{i, k}.

2.1.2. Computing the Transmembrane Current Based on the Solution of the Cable Equation

Next we address the problem of computing the transmembrane current based on the solution of the Cable equation. Suppose that the Cable equation is solved numerically using an implicit finite difference scheme of the form

$\begin{array}{ll}{C}_{\mathrm{\text{m}}}\frac{v_{n,k}-v_{n-1,k}}{\Delta t}+I_{\mathrm{\text{ion}},n,k}=\eta \frac{v_{n,k-1}-2v_{n,k}+v_{n,k+1}}{\Delta x^2},& (13)\end{array}$

where, as above, Δx denotes the spatial discretization in form of compartments, Δt denotes the time-step, and n is used to enumerate the time steps. Then, the associated transmembrane current density is given by

$\begin{array}{ll}{I}_{\mathrm{\text{m}}}^{k,n}=\eta \frac{v_{n,k-1}-2v_{n,k}+v_{n,k+1}}{\Delta x^2}.& (14)\end{array}$

All the methods discussed below for computing the extracellular potential rely on an approximation of this current density, but the methods differ in how the current is approximated and in the assumptions made on the extracellular domain.

2.1.3. Computing the Extracellular Potential in Terms of Solving a Boundary Value Problem; The CBV Method

Consider the simplified 2D geometry illustrated in Figure 2. Our aim is now to compute the extracellular potential in Ω_e for the given transmembrane currents computed as explained above. The problem we have to solve is given by

$\begin{array}{ll}{\nabla }^2{u}_{\mathrm{\text{e}}}=0,in{\Omega }_{\mathrm{\text{e}}},& (15)\end{array}$

$\begin{array}{ll}{\sigma }_{\mathrm{\text{e}}}\frac{\partial u_{\mathrm{\text{e}}}}{\partial n_{\mathrm{\text{e}}}}=I_{\mathrm{\text{m}}},at\Gamma ,& (16)\end{array}$

where I_m is computed by Equation (14) and n_e is the outward pointing normal vector of Ω_e. The boundary condition at the outer boundary of Ω_e will be described for the simulations presented below.

FIGURE 2

Figure 2. Sketch of a simplified neuron geometry and its surroundings; the extracellular domain Ω_e, the cell membrane Γ, and the intracellular domain Ω_i. The normal vector pointing out of Ω_i, is denoted by n_i and, similarly, n_e denotes the normal vector pointing out of Ω_e.

In our computations, the Laplace Equation (15) together with the boundary condition (16) is solved numerically using straightforward finite difference approximations leading to a linear system of algebraic equations. The finite difference stencil used for Equation (15) will be described below.

We will refer to this method for computing the extracellular potential as the CBV-method since it comprises the solution of the Cable equation (C) and the solution of a boundary value (BV) problem.

2.1.4. Computing the Extracellular Potential by Solving the Poisson Equation; The CP Method

In the CBV method the transmembrane currents setting up the extracellular potential are positioned at the interface between the intracellular and extracellular domains. In the standard method for computing extracellular potentials (referred to as the CS method below), the transmembrane currents are instead assumed to be positioned at the points (or lines) at the center of the intracellular domain (Holt and Koch, 1999). One step in this direction is to replace the boundary value problem (15, 16) of the CBV method with a Poisson equation of the form

$\begin{array}{ll}\nabla ·(\sigma \nabla u)=-C,in\Omega ,& (17)\end{array}$

where Ω = Ω_e ∪ Ω_i. Here σ equals σ_i and σ_e in Ω_i and Ω_e, respectively, and u equals u_i and u_e in Ω_i and Ω_e, respectively. The problem now is how to define the current source density C. In order to define C, we start by recalling that integration by parts gives

$\begin{array}{ll}{\displaystyle {\int }_{\Omega }}\nabla ·(\sigma \nabla u)\varphi dV={\displaystyle {\int }_{\Gamma }}\sigma \frac{\partial u}{\partial n}\varphi ds-{\displaystyle {\int }_{\Omega }}\nabla \varphi ·(\sigma \nabla u)dV& (18)\end{array}$

for any smooth functions u and ϕ, see e.g., Grossmann et al. (2007, p. 140). By choosing ϕ = 1, and using this identity, it follows from Equation (17) that the integral of C must be given by

$\begin{array}{ll}{\displaystyle {\int }_{{\Omega }_{\mathrm{\text{i}}}}}CdV=-{\displaystyle {\int }_{{\Omega }_{\mathrm{\text{i}}}}}\nabla ·({\sigma }_{\mathrm{\text{i}}}\nabla u_{\mathrm{\text{i}}})dV=-{\displaystyle {\int }_{\Gamma }}{\sigma }_{\mathrm{\text{i}}}\frac{\partial u_{\mathrm{\text{i}}}}{\partial n_{\mathrm{\text{i}}}}ds={\displaystyle {\int }_{\Gamma }}{I}_{\mathrm{\text{m}}}ds,& (19)\end{array}$

where n_i is the outward pointing normal vector of Ω_i. We now want to define the source term C such that the identity (Equation 19) holds. To this end, we define the constants²

$\begin{array}{ll}{C}_k=\frac{\left|{\Gamma }_k\right|}{\left|{\Omega }_{\mathrm{\text{i}},k}\right|}{I}_{\mathrm{\text{m}},k},& (20)\end{array}$

for every compartment Ω_i,k, where

$\begin{array}{ll}{I}_{\mathrm{\text{m}},k}=\frac{1}{\left|{\Gamma }_k\right|}{\displaystyle {\int }_{{\Gamma }_k}}{I}_{\mathrm{\text{m}}}ds& (21)\end{array}$

is the average transmembrane current density of the compartment. Based on these constants, we can define the source term

$\begin{array}{ll}C=C_k\text{\hspace{0.17em}for}x\in {\Omega }_{\text{i},k}.\hfill & (22)\hfill \end{array}$

With this definition of the source term, we have

${\displaystyle {\int }_{{\Omega }_{\mathrm{\text{i}}}}}Cdx={\displaystyle {\int }_{\Gamma }}{I}_{\mathrm{\text{m}}}ds$

and therefore Equation (19) holds provided that the current source density C is defined by (22).

We can now approximate the solution of the boundary value problem (15,16) defined on Ω_e with the Poisson problem (17) defined on the entire Ω = Ω_e ∪ Ω_i. It remains to be seen that the current flowing into the extracellular domain Ω_e defined by boundary condition (16) is the same as the amount of current flowing out of the intracellular domain Ω_i in the solution of the Poisson Equation (17). This holds, since by the definition of C we have

$\begin{array}{ll}{\displaystyle {\int }_{\Gamma }}{\sigma }_{\mathrm{\text{e}}}\frac{\partial u_{\mathrm{\text{e}}}}{\partial n_{\mathrm{\text{e}}}}ds={\displaystyle {\int }_{\Gamma }}{I}_{\mathrm{\text{m}}}ds=-{\displaystyle {\int }_{\Gamma }}{\sigma }_{\mathrm{\text{i}}}\frac{\partial u_{\mathrm{\text{i}}}}{\partial n_{\mathrm{\text{i}}}}ds.& (23)\end{array}$

Note that this method effectively assumes the transmembrane current to be homogeneously distributed in the intracellular domain in the computation of the extracellular potential. Again, the numerical solution of Equation (22) is obtained by the finite difference method where the right-hand side of the equation is evaluated in the mesh points. This leads to a linear system of algebraic equations.

The method of computing the extracellular potential by solving the Cable Equation (11), using the result to define the source term C by Equation (22), and then solving the Poisson Equation (17), will be referred to as the CP method (C is for Cable and P is for Poisson).

2.1.5. Computing the Extracellular Potential by the Point Source Method; The CS Method

The final method for computing the extracellular potential based on the solution of the Cable equation we will consider is the point source method. This method relies on two basic assumptions; first it is assumed that all the current can be gathered in the center of each compartment; and second, it is assumed that the extracellular space is infinite. Under these assumptions, the Poisson equation can be solved analytically, see e.g., Holt (1998), Holt and Koch (1999), Gold et al. (2006), and Einevoll et al. (2013). This dramatically increases computational efficiency and thus this approach is extremely popular and completely dominates computations of extracellular potentials around neurons. Again, our aim is to assess the accuracy of this method.

By using the notation introduced above, we define current sources for each compartment by

$\begin{array}{ll}{c}_k=\left|{\Omega }_{\mathrm{\text{i}},k}\right|C_k,& (24)\end{array}$

and define the associated Poisson problems

$\begin{array}{ll}{\sigma }_{\mathrm{\text{e}}}{\nabla }^2{u}_{\mathrm{\text{e}},k}=-c_k\delta (r-r_k),& (25)\end{array}$

where r = (x, y, z) and r_k is the center of the k−th compartment. The solution of this problem reads

$\begin{array}{ll}{u}_{\mathrm{\text{e}},k}=\frac{c_k}{4\pi {\sigma }_{\mathrm{\text{e}}}|r-r_k|},& (26)\end{array}$

and therefore, by linearity, the extracellular potential is given by

$\begin{array}{ll}{u}_{\mathrm{\text{e}}}={\displaystyle {\displaystyle \sum _k}}{u}_{\mathrm{\text{e}},k}=\frac{1}{4\pi {\sigma }_{\mathrm{\text{e}}}}{\displaystyle \sum _k}\frac{c_k}{|r-r_k|}.& (27)\end{array}$

Note that |r − r_k| denotes the Euclidean distance for r to the point r_k. In the computations below we will refer to this method of computing the extracellular potential as the CS-method (where C is for Cable and S is for sum).

2.2. The Extracellular-Membrane-Intracellular (EMI) Model

The dynamics of a neuron and its extracellular surroundings can be accurately modeled by explicitly considering the Extracellular space (Ω_e), the Membrane (Γ) and the Intracellular domain (Ω_i); as mentioned above we call this the EMI model. Analytical examples of solutions are given by Krassowska and Neu (1994), finite element formulations are provided by Ying and Henriquez (2007), Henríquez et al. (2013), and Agudelo-Toro and Neef (2013); see also Agudelo-Toro (2012) for a detailed derivation of the model.

The main elements of the model are sketched in Figure 2. Note that Ω = Ω_i ∪ Ω_e contains a single cell, where Ω_i is the intracellular domain of the cell and Ω_e is the extracellular space surrounding the cell. We let u_i and u_e denote the intra- and extracellular potentials, and at the interface between the intracellular and extracellular domains, given by Γ, we define the membrane potential by v = u_i − u_e. Then the electrical potential defined in Ω = Ω_i ∪ Ω_e is governed by the system

$\begin{array}{ll}\nabla ·{\sigma }_{\mathrm{\text{i}}}\nabla u_{\mathrm{\text{i}}}=0,in{\Omega }_{\mathrm{\text{i}}},& (28)\end{array}$

$\begin{array}{ll}\nabla ·{\sigma }_{\mathrm{\text{e}}}\nabla u_{\mathrm{\text{e}}}=0,in{\Omega }_{\mathrm{\text{e}}},& (29)\end{array}$

$\begin{array}{ll}{u}_{\mathrm{\text{e}}}=0,at\partial {\Omega }_{\mathrm{\text{e}}},& (30)\end{array}$

$\begin{array}{ll}{n}_{\mathrm{\text{e}}}·{\sigma }_{\mathrm{\text{e}}}\nabla u_{\mathrm{\text{e}}}=-n_{\mathrm{\text{i}}}·{\sigma }_{\mathrm{\text{i}}}\nabla u_{\mathrm{\text{i}}},at\Gamma ,& (31)\end{array}$

$\begin{array}{ll}{u}_{\mathrm{\text{i}}}-u_{\mathrm{\text{e}}}=v,at\Gamma ,& (32)\end{array}$

$\begin{array}{ll}{I}_{\mathrm{\text{m}}}=-n_{\mathrm{\text{i}}}·{\sigma }_{\mathrm{\text{i}}}\nabla u_{\mathrm{\text{i}}},at\Gamma ,& (33)\end{array}$

$\begin{array}{ll}\frac{\partial v}{\partial t}=\frac{1}{C_{\mathrm{\text{m}}}}(I_{\mathrm{\text{m}}}-I_{\mathrm{\text{ion}}}),at\Gamma ,& (34)\end{array}$

where σ_i and σ_e are intra- and extracellular conductivities, n_i and n_e are the normal vectors of Ω_i and Ω_e, C_m is the cell membrane capacitance, and the ion current density is given by I_ion.

2.2.1. Numerical Methods

We describe the finite difference scheme for solving the system (28)–(34) in the case of passive ion currents; i.e., for a case where I_ion is linear. In this case the problem (28)–(34) is linear and it is straightforward to define a fully implicit finite difference scheme. In the description of the solution method, we will consider the 2D case illustrated in Figure 3. The extension to 3D is notationally messy but conceptually straightforward.

FIGURE 3

Figure 3. Sketch of the computational mesh for Ω_e and Ω_i; the nodes of Ω_e are marked by “×,” the nodes of Ω_i are marked by “°,” and the membrane is defined as the intersection of Ω_e and Ω_i marked by “⊗.”

The system (28)–(34) can be triggered in several different ways. Since we want to compare results using the Cable equation and the EMI model, we will apply an initial condition that can be used in an identical manner for both methods. This is achieved by assuming that the membrane potential is given at time t = 0, and by adding a one dimensional synaptic input current.

We let $(u_{\mathrm{\text{e}}}^{n,j,k},v^{n,j,k},u_{\mathrm{\text{i}}}^{n,j,k})$ denote finite difference approximations of (u_e, v, u_i) at (t_n, x_j, y_k) = (nΔt, jΔx, kΔy) for given mesh parameters Δt, Δx and Δy. The computational nodes used for the discrete version of the system are shown in the right panel of Figure 3; nodes of the extracellular domain are marked by “×,” the intracellular nodes are marked by “◦,” and the nodes on the membrane are marked by “⊗.”

Suppose that v = vⁿ⁻¹ is known at time t = t_n−1. The update from t_n−1 to t_n is computed by solving a coupled linear system defined by a finite difference version of the system (28)–(34). In each node of the extracellular domain the elliptic Equation (29) is replaced by a finite difference scheme of the form

$\begin{array}{l}\frac{{\sigma }_{\mathrm{\text{e}}}^{j+1/2,k}(u_{\mathrm{\text{e}}}^{n,j+1,k}-u_{\mathrm{\text{e}}}^{n,j,k})-{\sigma }_{\mathrm{\text{e}}}^{j-1/2,k}(u_{\mathrm{\text{e}}}^{n,j,k}-u_{\mathrm{\text{e}}}^{n,j-1,k})}{\Delta x^2}\\ +\frac{{\sigma }_{\mathrm{\text{e}}}^{j,k+1/2}(u_{\mathrm{\text{e}}}^{n,j,k+1}-u_{\mathrm{\text{e}}}^{n,j,k})-{\sigma }_{\mathrm{\text{e}}}^{j,k-1/2}(u_{\mathrm{\text{e}}}^{n,j,k}-u_{\mathrm{\text{e}}}^{n,j,k-1})}{\Delta y^2}=0,& (35)\end{array}$

where ${\sigma }_{\mathrm{\text{e}}}^{j+1/2,k}={\sigma }_{\mathrm{\text{e}}}((j+1/2)\Delta x,k\Delta y).$ Likewise, the elliptic equation (28) is replaced by a finite difference scheme of similar form (u_e replaced by u_i and σ_e replaced by σ_i). The numerical scheme given by Equation (35) provides one equation for all nodes in the domain Ω_e\Γ and (as explained above) for all nodes in the domain Ω_i\Γ.

It remains to specify three equations for all nodes on the membrane Γ since there are three unknowns, (u_e, v, u_i), in each of the membrane nodes. One equation is clearly given by Equation (32); i.e., $u_{\mathrm{\text{i}}}^{n,j,k}-u_{\mathrm{\text{e}}}^{n,j,k}=v^{n,j,k}$ for all nodes (x_j, y_k) on the membrane Γ. The second equation is provided by replacing the flux-equality Equation (31) by a finite difference equation, and the third equation is the discrete version of Equation (34) in terms of an implicit scheme;

$v^{n,j,k}-\frac{\Delta t}{C_{\text{m}}}(I_{\text{m}}^{n,j,k}-I_{\text{ion}}^{n,j,k})=v^{n-1,j,k}.(36)$

Here, I_m is defined as a discrete version of Equation (33). Furthermore, in the passive case, the function I_ion is linear with respect to v and therefore the entire system is linear.

The four corners of the membrane mesh need special attention. In these nodes we define two extracellular and two intracellular flux terms; one term from the normal derivative in the x-direction and one from the normal derivative in the y-direction. Furthermore, we let the sum of the two intracellular fluxes equal the sum of the two extracellular fluxes in the flux-equality Equation (31) and let $I_{\mathrm{\text{m}}}^{n,j,k}$ in Equation (36) be the mean of the two intracellular fluxes. In the 3D extension we similarly define three extracellular and three intracellular flux terms for the corner nodes where three membrane planes intersect, and two extracellular and two intracellular flux terms for the edge nodes where two membrane planes intersect.

In the case of simple, rectangular geometries, this numerical strategy is straightforward. However, for more complex geometries, finite element or finite volume methods should be used.

3. Results

In this section we will report results using the methods described above. We will start the section by investigating the error in the membrane potential introduced by ignoring the ephaptic current Equation (12).

Secondly, we will compare the extracellular potential computed by the CBV, CP, and CS methods with the solution of the EMI model. Clearly, there are a set of different assumptions underlying these methods: The CS method is unique in assuming the extracellular domain to be infinite. In order to be able to compare the results of the CS method with the other methods, we have used large extracellular domains. In order to estimate how large the domain must be, we have systematically increased the size of the extracellular space until convergence of the EMI solutions and then used the largest domain for our comparisons.

For the CS method the transmembrane currents are gathered in the center of each compartment thus giving rise to the classical formula of the solution, whereas for the CP method the transmembrane currents are distributed over each compartment, and numerical methods are used to compute the solution of the associated Poisson equation. In contrast, in the CBV and EMI methods the transmembrane currents setting up the extracellular potentials are placed at the interface between the intracellular and extracellular domains. The CBV and EMI methods are thus defined on the same domain, and the only difference lies in the proper self-consistent modeling of ephaptic effects in the EMI method.

For convenience, the abbreviations (CBV, CP, CS, and EMI) and references to the methods are summarized in Table 1.

TABLE 1

Table 1. Definition of the methods used to compute the extracellular potential.

3.1. Model Parameters

We consider the Cable equation and the EMI model using the parameters given in Table 2 (unless otherwise stated). The domain Ω = Ω_i ∪ Ω_e is defined as

$\begin{array}{ll}\Omega =\left[0,L_x\right]\times \left[0,L_y\right]\times \left[0,L_z\right],& (37)\end{array}$

and the intracellular domain, Ω_i, is shaped as a rectangular cuboid of size l_x × l_y × l_z located in the center of Ω. The ionic current density I_ion is defined as

$\begin{array}{ll}{I}_{\mathrm{\text{ion}}}=I_{\mathrm{\text{leak}}}+I_{\mathrm{\text{syn}}},& (38)\end{array}$

where I_leak is the leak current density given by

$\begin{array}{ll}{I}_{\mathrm{\text{leak}}}=g_L(v-v_{\mathrm{\text{rest}}}),& (39)\end{array}$

and I_syn is the conductance-based synaptic current density with single-exponential dynamics (see Gerstner et al., 2014) given by

$\begin{array}{ll}{I}_{\mathrm{\text{syn}}}=g_s(x)e^{-\frac{t-t_0}{\alpha }}(v-v_{\mathrm{\text{eq}}}).& (40)\end{array}$

TABLE 2

Table 2. Parameters used in the computations of the Cable equation and the EMI model.

For the first 10% of the cell in the x-direction, g_s(x) is given by the value g_syn in Table 2. On the remaining part of the membrane g_s(x) is set to zero.

We use the initial condition v = v_rest = −90 mV for the membrane potential. In addition, we apply the boundary condition $\frac{\partial v}{\partial x}=0$ at the start and the end of the cell in the Cable equation and the boundary condition u_e = 0 on the outer boundary of Ω_e in the EMI, CBV and CP methods unless otherwise stated.

3.2. Numerical Assessment of the Error in Membrane Potential Introduced by Ignoring the Ephaptic Current

In Figure 4 we show the membrane potential computed by solving the Cable equation and the EMI model for different values of h, σ_i, σ_e, and g_L. The solutions are compared in the compartment 25 μm from the start of the cell (i.e., in the center of the cell in the x-direction). The difference is several millivolts, but it is reduced as the intracellular conductivity σ_i is reduced or the size h (recall that h = l_y = l_z) of the neuron is reduced. Furthermore, we observe that the difference is reduced as the extracellular conductivity, σ_e, or g_L is increased. These observations are consistent with our theoretical finding in an Appendix in Supplementary Material given below where we show that, under reasonable assumptions, the error introduced in the transmembrane potential by removing the ephaptic current goes like

$\begin{array}{ll}O\left(\frac{h{\sigma }_{\mathrm{\text{i}}}}{g_L{\sigma }_{\mathrm{\text{e}}}}\right).& (41)\end{array}$

FIGURE 4

Figure 4. Comparison of the membrane potential computed by solving the Cable equation (red) and the EMI model (blue) for some different values of h, σ_i, σ_e, and g_L, where we recall that h = l_y = l_z (the width of the neuron). The plots show how the membrane potential in the compartment 25 μm from the start of the cell changes with time from t = 0.1 to 0.5 ms. The parameters used in the computations are given in Table 2 except for the values given above each plot. We observe that the difference between the two solutions increases when the value of h or σ_i is increased, and the difference decreases when the value of σ_e or g_L is increased. Note that in order to observe any effect of changing the value of g_L, we increase the default value by a factor of order 100–1,000 in the lower panel of the figure.

To summarize, the error increases when h or σ_i are increased, and the error decreases if g_L or σ_e are increased.

3.3. The Magnitude of the Ephaptic Current Decreases as the Extracellular Conductivity is Increased

As mentioned above, the derivation of the Cable equation relies on the assumption that the extracellular potential is constant, and under that assumption, the ephaptic current defined by Equation (12) can be ignored. This can also be understood on biophysical grounds as a high extracellular conductivity implies a low extracellular resistance so that potential drops due to extracellular currents driven through the extracellular medium will be small. In the limit of very large extracellular conductivities these potential drops will become negligible, i.e., the assumption of constant extracellular potentials in the standard Cable equation will become fulfilled.

In Table 3 the maximum magnitude (absolute value) of the ephaptic current (computed by solving the EMI model) is given as a function of the extracellular conductivity σ_e, and we note that the magnitude decreases as σ_e is increased. In addition, we report the value of the maximum ephaptic current multiplied by the value of σ_e and observe that this value is close to a constant, so we have

$\begin{array}{ll}{I}_{\mathrm{\text{eph}}}~O(1/{\sigma }_{\mathrm{\text{e}}}).& (42)\end{array}$

TABLE 3

Table 3. Maximum absolute values of I_eph from time t = 0.02 to 1 ms as a function of σ_e as computed by the EMI method.

Therefore, for very large values of σ_e, the ephaptic current can be ignored, but the reported values of σ_e are in general not so large that this assumption can be generally trusted.

It is also interesting to compare the size of the ephaptic current with the size of the other currents involved in the dynamics of the model neuron. Figure 5 shows the time evolution of each of the terms in Equation (9) and we observe that the size of the ephaptic current is comparable to the size of the other terms in the equation. The peak of the ephaptic current is located at the jump in the synaptic input.

FIGURE 5

Figure 5. Values of each of the terms in Equation (9). In the (Upper panel), we show the time evolution of the terms in the point (8, 10, 7 μm) inside the synaptic input zone and the point (12, 10, 7 μm) outside the synaptic input zone. In the (Lower panel), we show the values of the terms for y = 10 μm, z = 7 μm, and x ∈ [5 μ m, 30μm] at time t = 0.02 ms (left) and t = 0.2 ms (right). The solution of the EMI model is used to compute each of the terms. In addition, we show $\eta \frac{{\partial }^{2}v}{\partial x^2}$ for the corresponding solution of the Cable equation, where I_eph is assumed to be zero. We observe that the size of I_eph is comparable to the size of the other terms in Equation (9) and that neglecting I_eph leads to a considerable difference in the value of the term $\eta \frac{{\partial }^{2}v}{\partial x^2}$. The parameters used in the computations are given in Table 2.

3.4. Comparing the Extracellular Potential Computed by the CBV, CP, CS, and EMI Methods

In this section we will compare the extracellular potentials (EPs) computed by the EMI, CBV, CP and CS methods described above (see Table 1 for definitions of the abbreviations). When comparing the predicted extracellular potentials for the various methods, observed differences will expectedly have different model origins. For the EMI and CBV methods the key physical difference is in the lack of inclusion of ephaptic effects in the CBV method. Compared to EMI and CBV where the transmembrane currents setting up the EP are at the true membrane interface between the intracellular and extracellular domains, the CP and CS methods assume that the EP-generating currents are defined as the right-hand side of the Poisson Equation (17). For the CS method the current source density is gathered in a single point in the center of the neuronal compartment, whereas for the CP method the current density is evenly distributed over the entire compartment (See Figure 1).

3.4.1. Convergence under Mesh Refinements

In Figure 6 we show the extracellular potential computed by the EMI method for four different values of the discretization parameter Δx = Δy = Δz. The solutions for the 0.5 μm resolution and the 0.25 μm resolution appear to be similar and we use a spatial discretization of Δx = Δy = Δz = 0.5 μm for the rest of our computations.

FIGURE 6

Figure 6. Extracellular potential computed by the stationary EMI model for four different values of Δx = Δy = Δz. We show the solution in a rectangle of size 60 × 30 μm on the plane in the center of the domain in the z-direction. The white area represents the cell. We use the parameters given in Table 2 except for an increased value of $g_L=3·10^{-5}$ μS/μm² and a domain of size 60 × 60 × 60 μm.

To reduce the computational cost in this case, we consider the stationary version of the model, i.e., we set the time derivative in Equation (34) to zero. We use the parameter values given in Table 2, except for an increased value of $g_L=3·10^{-5}$ μS/μm² and a domain of size 60 × 60 × 60 μm. We again let g_s(x) be g_syn for the first 10% of the cell in the x-direction and zero elsewhere and apply the boundary condition u_e = 0 on the outer boundary of the extracellular domain.

3.4.2. Convergence of the EMI Solution as the Domain Size is Increased

In the derivation of the CS method, the extracellular domain is assumed to be infinite (see Section 2.1.5). When comparing CS and EMI results, we therefore wish to compare the solution of the CS method to the solution of the EMI model as the size of the extracellular domain approaches infinity.

We again consider the stationary version of the model with the parameter values given in Table 2, except for an increased value of $g_L=3·10^{-5}$ μS/μm² and an increased domain size.

Figure 7 shows the stationary solution of the EMI model for four different sizes of the extracellular domain. We observe that as the size of the extracellular domain increases, the solution of the EMI model appears to converge, and we assume that the solution for a domain of size 120 × 120 × 120μm is sufficiently large to represent the EMI solution of an infinite domain.

FIGURE 7

Figure 7. Comparison of the extracellular potential around a neuron computed by the stationary EMI model for four different sizes of the extracellular domain. The plots to the left show the solution in a rectangle of size 60 × 30 μm on the plane in the center of the domain in the z-direction. The white area represents the neuron. The plot to the right shows the extracellular potential along a line 2 μm above the neuron in the y-direction and in the center of the domain in the z-direction. The parameters used in the computations are given in Table 2 except for L_x, L_y, and L_z, which are specified for each simulation, and g_L, which is set to 3 · 10⁻⁵ μS/μm².

3.4.3. One Single Simplified Neuron

Our first test case for comparing the methods for computing the extracellular potential is a single neuron of the form given above. The extracellular potential computed by the CBV, CP, CS, and EMI methods are presented in Figure 8 (see Table 1 for definitions of the abbreviations). In Table 4 we report the maximum difference between the extracellular potential computed by the EMI model and the extracellular potential computed by each of the other methods. The deviation of the CBV result from the EMI result is smaller than the difference to the CP and CS results. Thus the largest differences appear to come from the different assumptions of placement of the transmembrane currents in the EP-generating step (compare CBV vs. CP and CS). The effect of the ephaptic current (CBV vs. EMI) is smaller.

FIGURE 8

Figure 8. The extracellular potential around a neuron shaped as a rectangular cuboid computed by the stationary versions of the EMI, CBV, CP, and CS methods. The plots to the left show the solution in a rectangle of size 60 × 30 μm on the plane in the center of the domain in the z-direction. The white area represents the neuron. The plot to the right shows the extracellular potential along a line 2 μm above the neuron in the y-direction and in the center of the domain in the z-direction. We use the parameters specified in Table 2 except for L_x = L_y = L_z = 120 μm and $g_L=3·10^{-5}$ μS/μm². The abbreviations (EMI, CBV, CP, and CS) are summarized in Table 1.

TABLE 4

Table 4. Maximum difference between the solution for the extracellular potential in Ω_e\Γ computed by the EMI method and each of the other methods for the test case in Figure 8.

3.4.4. Two Simplified Neurons

In Figure 9 we show the extracellular potential around two neurons of the form given above computed by the CBV, CP, CS, and EMI methods. In the upper part of the figure the neurons are separated by a distance of 10 μ m in the y-direction and in the lower part the neurons are separated by a distance of 4 μm. In Table 5 we report the maximum difference between the extracellular potential computed by the EMI method and each of the other methods for the two test cases.

FIGURE 9

Figure 9. The extracellular potential around two neurons computed by the stationary versions of the EMI, CBV, CP, and CS methods. The plots to the left show the solution in a rectangle of size 60 × 40 μm on the plane in the center of the domain in the z-direction. The white areas represent the neurons. The plots to the right show the extracellular potential along the line in the center of the space between the two neurons. In the upper five plots, the neurons are separated by a distance of 10 μm in the y-direction, and in the lower five plots the neurons are separated by a distance of 4 μm. In all plots g_s(x) is given by g_syn for x ∈ [55, 60 μm] and is zero on the rest of the membrane for the lower neuron. For the upper neuron g_s(x) is given by g_syn for x ∈ [60 μm, 65 μm]. We use the parameters specified in Table 2 except for L_x = L_y = L_z = 120 μm and $g_L=3·10^{-5}$ μS/μm². The abbreviations (EMI, CBV, CP, and CS) are summarized in Table 1.

TABLE 5

Table 5. Maximum difference between the solution for the extracellular potential in Ω_e\Γ computed by the EMI method and each of the other methods for the test cases in Figure 9.

As for the case with a single simplified neuron above, the deviation of the CBV result from the EMI result is seen to be smaller than the difference to the CP and CS results. Interestingly, in the lower part of Figure 9 where the distance between the two neurons is very small (4 μm), the EMI and CBV results are essentially identical in the space between the cells.

3.4.5. Confined Extracellular Space

Figure 10 shows the extracellular potential around a neuron in a domain of size 60 × 20 × 20μm computed by the EMI, CBV, CP, and CS methods. The left panel shows the solution for a homogeneous Dirichlet boundary condition, and the right panel shows the solution for a homogeneous Neumann boundary condition on the outer boundary of the extracellular space.

FIGURE 10

Figure 10. Extracellular potential around a neuron computed by the EMI, CBV, CP, and CS methods. We consider the stationary version of the models and the parameter values given in Table 2 except for an increased value of $g_L=3·10^{-5}\mu$m. A Dirichlet boundary condition, u_e = 0, is applied in the simulation in the left panel and a Neumann boundary condition, $\frac{\partial u_{\mathrm{\text{e}}}}{\partial n_{\mathrm{\text{e}}}}=0$, is applied in the right panel. The (Upper panels) show the extracellular potential in the plane in the center of the domain in the z-direction for each of the methods. The (Lower panel) shows the solution along a line 2 μm above the cell in the y-direction and in the center of the domain in the z-direction. Note that in the case of Neumann boundary conditions, we include the additional constraint ${\int }_{{\Omega }_{\mathrm{\text{e}}}}{u}_{\mathrm{\text{e}}}dV=0$ for the EMI, CBV, and CP methods in order to obtain unique solutions. The abbreviations (EMI, CBV, CP, and CS) are summarized in Table 1.

As explained above, the CS method is founded on the assumption of an infinite extracellular space. We have therefore focused on a very large computational domain mimicking the properties of an infinite domain. Certainly, also limited domains are of interests and simulation results are given in Figure 10 using both Dirichlet and Neumann type boundary conditions. Although we present results for all four models, it is important to keep in mind that a confined domain breaks a basic assumption underlying the CS method and consequently we get very large errors, especially in the case of Neumann type boundary conditions.

Note that in the case of Neumann boundary conditions, the solution is not uniquely determined by the systems defining the EMI, CBV, and CP methods, and we expand the systems with the additional constraint

$\begin{array}{ll}{\displaystyle {\int }_{{\Omega }_{\mathrm{\text{e}}}}}{u}_{\mathrm{\text{e}}}dV=0& (43)\end{array}$

in order to obtain unique solutions of the methods.

3.4.6. Effects of the Size of the Synaptic Input Area

In Figure 11 we show the extracellular potential surrounding a neuron for four different sizes of the synaptic input area. The upper panel shows the extracellular potential computed by the EMI method, and the lower panel shows a comparison of the extracellular potentials computed by each of the methods along a line above the neuron. We note from the simulations that the results are qualitatively similar for all different sizes of the synaptic input region. Therefore, we choose to focus on the 10% synaptic input region as the base case for our simulations.

FIGURE 11

Figure 11. Extracellular potential around a neuron with a synaptic input area of length 5, 10, 20, and 30% of the total cell length. The (Upper panel) shows the extracellular potential computed by the EMI method in the plane in the center of the domain in the z-direction. The (Lower panel) shows the solution for each of the methods along a line 2 μm above the cell in the y-direction and in the center of the domain in the z-direction. The figure shows the solution of the stationary version of the models using the parameter values given in Table 2 except for an increased value of $g_L=3·10^{-5}\mu$m. We apply a homogeneous Dirichlet boundary condition on the outer boundary of the extracellular domain. The abbreviations (EMI, CBV, CP, and CS) are summarized in Table 1.

3.4.7. Simulation Time

Table 6 shows the CPU time for the simulations shown in the left panel of Figure 10 using a direct and an iterative solver.

TABLE 6

Table 6. CPU time in seconds for the EMI, CBV, CP, and CS methods.

The EMI model is clearly much more computationally expensive than the classical CS method. This is expected because the EMI model involves solving a large coupled system of equations, whereas the CS method only requires solving the Cable equation which involves a much smaller number of unknowns. After solving the Cable equation, the CS methods assumes that the extracellular potential may be found directly by the explicit formula (27), so no further equations has to be solved.

Moreover, in the computations reported in the table, the extracellular potential is computed for all nodes in the mesh. In the CS method this is not necessary, and the CPU time for the CS method could possibly be further reduced by only computing the values for the points of interest. This is not possible for the EMI method (or the CBV or CP methods) because the systems of equations has to be solved for all nodes in order to find the solution in a single point.

In contrast, the simulation time for the CBV and CP methods are more comparable to that of the EMI model, at least for the direct solver. This is because these methods also rely on solving a linear system of equations for all nodes in the extracellular domain or the entire domain for the CVB and CP methods, respectively.

The extra complexity introduced in the EMI model by solving for the membrane, intracellular and extracellular potentials simultaneously is apparent, however, when an iterative method is applied to solve the linear system. The fourth column of Table 6 shows the solution time for each of the methods using the bistable conjugate gradient stabilized method with an incomplete LU preconditioner. In this case, the CBV and CP methods are much faster than the EMI method.

4. Discussion

In the present paper we have compared four different methods for computing neural dynamics. In the numerically comprehensive EMI model the intracellular and extracellular dynamics are solved self-consistently and the membrane potentials and extracellular potential are computed simultaneously. For the other methods (CBV, CP, CS; see Table 1 for definitions of abbreviations) the membrane potential is first computed using the Cable equation, and the resulting transmembrane currents are used in a second step to compute the extracellular potential.

In the CBV method the transmembrane currents are placed on the interface between the intracellular and extracellular domains, and the only difference with the EMI model is the lack of self-consistency in the two-step computational scheme inherent in the CBV scheme, that is, the transmembrane currents are first computed using the Cable equation assuming a constant extracellular potential, while a non-constant potential (both in space and time) is computed in the second step.

For the CP and CS methods an additional assumption is made in the second step, namely that the effect of the transmembrane currents are assumed to be represented in terms of currents source densities. Specifically, for the CP method, the current source density is distributed evenly over a neuronal compartment and a numerical scheme is used to solve the resulting Poisson Equation (17), and for the CS method the source density is concentrated in a single point and thus the classical sum formula (27) of the solution can be applied.

4.1. Ignoring the Ephaptic Current

4.1.1. Error in Membrane Potential Introduced by Ignoring the Ephaptic Current

To study the error introduced by ignoring the ephaptic current in the Cable equation, we compared the membrane potential computed by solving the Cable equation to the corresponding solution of the EMI model. In our simple test case, we found that the membrane potential computed by the Cable equation could differ several millivolts from the solution of the EMI model and that the magnitude of the error seems to decrease with the value of the intracellular conductivity, σ_i, and the cell width, h. This suggests that the Cable equation is applicable for computing the membrane potential for sufficiently thin dendrites.

4.1.2. Ephaptic Current Decreases with Increasing Extracellular Conductivity

In the derivation of the Cable equation, it is assumed that the extracellular conductivity σ_e is so large that the extracellular potential varies very little in space and can be assumed to be a constant. As a result, the ephaptic current I_eph will be zero and may be removed from the model. In our numerical simulations of the EMI model, we confirmed that the size of I_eph decreases when the value of σ_e is increased (see Table 3). In fact, we found that the maximum absolute value of I_eph appeared to be inversely proportional to the value of σ_e. However, we also observed that the magnitude of I_eph was similar in size to the other currents involved in the model (see Figure 5). This suggests that I_eph is not negligible for the stylized neuron geometries and model parameters chosen here.

4.2. Error in Neglecting Ephaptic Currents

The CBV and EMI methods are defined on identical domains, and the key physical difference between the methods is the absence of ephaptic effects in the CBV method. Comparisons of computed extracellular potentials indeed show such ephaptic effects of varying magnitudes, both for the extracellular potential outside a single activated neuron (Figure 8) and between two activated neurons (Figure 9).

4.3. Effects of Position of Transmembrane Currents

To explore the effects of assumed positions of transmembrane currents, it is easiest to compare results from the three methods, i.e., CBV, CP, and CS, where the transmembrane currents in all cases are found from the Cable equation. Here effects from the ephaptic current are in all cases absent. For the present examples we observe that CS and CP results are typically quite similar, but both quite different from the CBV results (Figures 8, 9). From the point-source formula in (27) we see that the contribution to the extracellular potential from a point current source is inversely proportional to distance, and it is thus not surprising that this difference in assumed source positions has a sizeable effect on the predicted extracellular potentials.

4.4. Effects of Size of Extracellular Domain

Both in the CP and CBV methods (as well as in the EMI model) the extracellular domain is finite, while in the CS method the extracellular domain is infinite so that the solution of the Poisson equation can be given as an explicit sum. With a very small extracellular domain, corresponding to a small piece of brain tissue embedded in an insulator, large deviations from the infinite-domain results will be observed (Figure 7; see also Figure 10).

In order to compare results with the CS method, we here computed the EMI solution for gradually larger domains until the solutions appeared to converge. Further, we regarded the converged solution as the solution of the EMI problem for an infinite extracellular space, i.e., we estimated that the difference of the results for the largest considered domain and a (hypothetical) infinite domain was negligible for the present purposes. Roughly speaking, convergence was obtained for an extracellular space extending twice the length of the cable in every direction.

4.5. The Simplified Geometry

Today, simulations of neurons typically use much more complex and realistic geometries than what has been applied here. Already in Clark and Plonsey (1968) were able to analytically evaluate the extracellular potential of a cylindrical neuron, and it is certainly of interest to evaluate the models and methods discussed here in more realistic geometries. The generic limitation of the finite difference method used here is that it is hard to apply, correctly, to non-rectangular geometries. In Agudelo-Toro and Neef (2013), the finite element method is used and this gives much more freedom to represent realistic geometries. However, the code used in Agudelo-Toro and Neef (2013) required extremely fine times steps and we therefore focused on simplified geometries where the problem could be solved using a reasonable number of time steps. The solution of the EMI model using an implicit formulation will likely reduce the time step restrictions and this is subject for ongoing investigations.

Another limitation in the present study is the size of the extracellular space. This space is actually quite limited (see e.g., Syková and Nicholson, 2008), but the assumption of an infinite extracellular space is necessary for the application of the classical CS method (summation method), and thus we have used very large extracellular domains in order to provide fair comparisons with the classical model at the cost of simulating more realistic volumes.

4.6. Neural Tissue

The EMI model provides a useful framework for accurate computations of the electrophysiology of a small number of cells and their surroundings. It is, however, very hard to apply this methodology to neural tissue consisting of huge numbers of cells. In simulations of cardiac tissue, the Bidomain approach has successfully been applied to simulate the electrophysiology, see e.g., Keener and Sneyd (2009), Franzone et al. (2014), Sundnes et al. (2006); Roth (2013), and Trayanova (2011). Recently, a similar approach has been applied to neural tissue, see Meffin et al. (2014) and Tahayori et al. (2014). Most likely, some form of homogenization process is needed to derive tractable mathematical models for neural tissue.

4.7. Possible Additive Effects for Non-linear Membrane Dynamics

We have focused on a linear membrane model in order to highlight the effect of removing the ephaptic current in the simplest possible case. More generally, the question is whether ephaptic coupling would constitute a “feedback” mechanism with electric fields altering the activity of the same neural elements that gave rise to them in the first place, see Anastassiou and Koch (2015). For a linear model, this feedback mechanism was recently found to be the small but not negligible, see Goldwyn and Rinzel (2016), which clearly is consistent with our findings. However, the effect may very well be larger for non-linear models of the membrane dynamics; small electric fields can be amplified by non-linear effects, see Radman et al. (2007). At present, we have not conducted systematic simulations using a non-linear membrane model.

4.8. Other Assumptions

We note also that the current study was limited to standard simulation frameworks in neuroscience, where intra- and extracellular currents are assumed to be purely Ohmic, so that Equations (28) and (29) apply in the bulk solutions. That is, we did not include possible contributions from advective currents, displacement currents and ionic diffusion currents. These currents are typically neglected, as they are believed to play negligible roles for the system electrodynamics under most biophysically relevant conditions. However, computational studies have indicated that at least ionic diffusion could, in some scenarios, influence electrical potentials (see e.g., Qian and Sejnowski, 1989; Bédard and Destexhe, 2013; Halnes et al., 2013, 2016; Pods et al., 2013; Pods, 2017). These effects were not accounted for in the current study. We also confined our simulations to linear, passive membranes even if it known that active voltage-gated channels affect the extracellular potential; see e.g., Ness et al. (2016).

5. Conclusion

We have compared various methods for computing membrane potentials and extracellular potentials. For the simple test cases considered here, non-negligible errors were observed when neglecting ephaptic effects, i.e., when comparing results from the EMI model with the CBV model building on results from the Cable equation. Further, substantial differences in the predicted extracellular potentials were observed depending on whether transmembrane current sources were assumed to be placed in the center of the neural compartment or at the membrane interfaces. This study motivates further analysis of the errors for computations based on more realistic representations of the geometry and dynamics of the neurons using the EMI model.

Author Contributions

Writing paper: AT, KJ, GH, GE. Programming: KJ, LP, GL. Simulation: KJ, AT. Numerical methods: AT, KJ, JS, GL. Discussion: AT, KJ, TM, AE, GH, GE. Outline of research: all authors.

Conflict of Interest Statement

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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fncom.2017.00027/full#supplementary-material

Footnotes

1. ^Assuming that the extracellular potential is a linear function of position would result in the same simplified model.

2. ^Note that C_k is a constant defined on the compartment Ω_{i, k}; it is constant in space but varies in time.

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.

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Arvanitaki, A. (1942). Effects evoked in an axon by the activity of a contiguous one. J. Neurophysiol. 5, 89–108.

Google Scholar

Bédard, C., and Destexhe, A. (2013). Generalized cable theory for neurons in complex and heterogeneous media. Phys. Rev. E 88:022709. doi: 10.1103/PhysRevE.88.022709

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhalla, U. S. (2012). “Multi-compartmental models of neurons,” in Computational Systems Neurobiology, ed N. Le Novère (Dordrecht; Heidelberg; London; New York, NY: Springer), 193–225.

Google Scholar

Bokil, H., Laaris, N., Blinder, K., Ennis, M., and Keller, A. (2001). Ephaptic interactions in the mammalian olfactory system. J. Neurosci. 21, 1–5.

PubMed Abstract | Google Scholar

Buzsáki, G., Anastassiou, C. A., and Koch, C. (2012). The origin of extracellular fields and currents – EEG, ECoG, LFP and spikes. Nat. Rev. Neurosci. 13, 407–420. doi: 10.1038/nrn3241

PubMed Abstract | CrossRef Full Text | Google Scholar

Clark, J., and Plonsey, R. (1968). The extracellular potential field of the single active nerve fiber in a volume conductor. Biophys. J. 8, 842–864.

PubMed Abstract | Google Scholar

Dayan, P., and Abbott, L. F. (2001). Theoretical Neuroscience. MIT Press.

Google Scholar

Einevoll, G. T., Kayser, C., Logothetis, N. K., and Panzeri, S. (2013). Modelling and analysis of local field potentials for studying the function of cortical circuits. Nat. Rev. Neurosci. 14, 770–785. doi: 10.1038/nrn3599

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, G. B., and Terman, D. H. (2010). Mathematical Foundations of Neuroscience, Vol. 35. Berlin; Heidelberg; New York, NY: Springer Science & Business Media.

Google Scholar

Franzone, P. C., Pavarino, L. F., and Scacchi, S. (2014). Mathematical Cardiac Electrophysiology. Cham; Heidelberg; New York, NY; Dordrecht; London: Springer International Publishing.

Google Scholar

Gerstner, W., Kistler, W. M., Naud, R., and Paninski, L. (2014). Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition. Cambridge: Cambridge University Press.

Google Scholar

Gold, C., Henze, D. A., Koch, C., and Buzsáki, G. (2006). On the origin of the extracellular action potential waveform: a modeling study. J. Neurophysiol. 95, 3113–3128. doi: 10.1152/jn.00979.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Grossmann, C., Roos, H.-G., and Stynes, M. (2007). Numerical Treatment of Partial Differential Equations. Berlin; Heidelberg; New York, NY: Springer.

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Henríquez, F., Jerez-Hanckes, C., and Altermatt, M. F. R. (2013). “Dynamic finite-element model of axon extracellular stimulation,” in 6th International IEEE/EMBS Conference on Neural Engineering (NER), 2013, 589–592.

Holt, G. R. (1998). A Critical Reexamination of Some Assumptions and Implications of Cable Theory in Neurobiology. Ph.D. thesis, California Institute of Technology.

Holt, G. R., and Koch, C. (1999). Electrical interactions via the extracellular potential near cell bodies. J. Comput. Neurosci. 6, 169–184.

PubMed Abstract | Google Scholar

Katz, B., and Schmitt, O. H. (1940). Electric interaction between two adjacent nerve fibres. J. Physiol. 97, 471.

PubMed Abstract | Google Scholar

Keener, J. P., and Sneyd, J. (2009). Mathematical Physiology. New York, NY: Springer.

Google Scholar

Klee, M., and Rall, W. (1977). Computed potentials of cortically arranged populations of neurons. J. Neurophysiol. 40, 647–666.

PubMed Abstract | Google Scholar

Koch, C. (1999). Biophysics of Computation. New York, NY; Oxford: Oxford University Press.

Google Scholar

Koch, C., and Buice, M. A. (2015). A biological imitation game. Cell 163, 277–280.

PubMed Abstract | Google Scholar

Krassowska, W., and Neu, J. C. (1994). Response of a single cell to an external electric field. Biophys. J. 66, 1768–1776.

PubMed Abstract | Google Scholar

Lindén, H., Hagen, E., Leski, S., Norheim, E. S., Pettersen, K. H., and Einevoll, G. T. (2014). LFPy: A tool for biophysical simulation of extracellular potentials generated by detailed model neurons. Front. Neuroinform. 7:41. doi: 10.3389/fninf.2013.00041

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Meffin, H., Tahayori, B., Sergeev, E. N., Mareels, I. M., Grayden, D. B., and Burkitt, A. N. (2014). Modelling extracellular electrical stimulation: III. Derivation and interpretation of neural tissue equations. J. Neural Eng. 11:065004. doi: 10.1088/1741-2560/11/6/065004

CrossRef Full Text | Google Scholar

Ness, T. V., Remme, M. W., and Einevoll, G. T. (2016). Active subthreshold dendritic conductances shape the local field potential. J. Physiol. 594, 3809–3825. doi: 10.1113/JP272022

PubMed Abstract | CrossRef Full Text | Google Scholar

Pods, J. (2017). A comparison of computational models for the extracellular potential of neurons. J. Integr. Neurosci. 16, 19–32.

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Qian, N., and Sejnowski, T. (1989). An electro-diffusion model for computing membrane potentials and ionic concentrations in branching dendrites, spines and axons. Biol. Cybern. 62, 1–15.

Google Scholar

Radman, T., Su, Y., An, J. H., Parra, L. C., and Bikson, M. (2007). Spike timing amplifies the effect of electric fields on neurons: implications for endogenous field effects. J. Neurosci. 27, 3030–3036. doi: 10.1523/JNEUROSCI.0095-07.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Rall, W. (1962). Electrophysiology of a dendritic neuron model. Biophys. J. 2, 145–167.

PubMed Abstract | Google Scholar

Rall, W. (1969). Time constants and electrotonic length of membrane cylinders and neurons. Biophys. J. 9, 1483–1508.

PubMed Abstract | Google Scholar

Roth, B. J. (2013). Bidomain simulations of defibrillation: 20 years of progress. Heart Rhythm 10, 1218–1219. doi: 10.1016/j.hrthm.2013.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Scott, A. (2002). Neuroscience: A Mathematical Primer. New York, NY: Springer Science & Business Media.

Google Scholar

Shneider, M. N., and Pekker, M. (2015). Correlation of action potentials in adjacent neurons. Phys. Biol. 12:066009. doi: 10.1088/1478-3975/12/6/066009

PubMed Abstract | CrossRef Full Text | Google Scholar

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

Google Scholar

Sundnes, J., Lines, G., Cai, X., Nielsen, B., Mardal, K.-A., and Tveito, A. (2006). Computing the Electrical Activity of the Heart. Heidelberg: Springer.

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Tahayori, B., Meffin, H., Sergeev, E. N., Mareels, I. M., Burkitt, A. N., and Grayden, D. B. (2014). Modelling extracellular electrical stimulation: IV. Effect of the cellular composition of neural tissue on its spatio-temporal filtering properties. J. Neural Eng. 11:065005. doi: 10.1088/1741-2560/11/6/065005

CrossRef Full Text | Google Scholar

Trayanova, N. A. (2011). Whole-heart modeling: applications to cardiac electrophysiology and electromechanics. Circ. Res. 108, 113–128. doi: 10.1161/CIRCRESAHA.110.223610

PubMed Abstract | CrossRef Full Text | Google Scholar

Tveito, A., and Winther, R. (2009). Introduction to Partial Differential Equations: A Computational Approach, Vol. 29, 2nd Edn. New York, NY: Springer-Verlag.

Google Scholar

Ying, W., and Henriquez, C. S. (2007). Hybrid finite element method for describing the electrical response of biological cells to applied fields. IEEE Trans. Biomed. Eng. 54, 611–620. doi: 10.1109/TBME.2006.889172

PubMed Abstract | CrossRef Full Text | Google Scholar