A three-dimensional mathematical model for the signal propagation on a neuron's membrane

Introduction

The membrane potential belongs to the most important quantities of a neuron. Its function of time and space describes neuronal activity. It is a voltage across the membrane defined by the difference between the intra- and extracellular potential.

Since the neuron is embedded in ionic milieus, potential gradients in the off-membrane spaces result in electric fluxes, which are conserved according to the first principles. This conservation law is the basis of the standard cable equation which describes the unfolding and propagation of an action potential (Rall, 1962, 1964; Scott, 1975) very efficiently. The standard cable equation maps a neuron to a tree of lines, each of which corresponds to a cylindric compartment with mean diameter. On these structures, it computes the evolution of the membrane potential according to its diffusion equation.

The resulting extracellular potentials can be theoretically computed with the line source method (Holt and Koch, 1999; Gold et al., 2006), once the transmembrane currents have been determined with the aid of the cable equation's solution.

These extracellular potentials in turn can be exploited to examine ephaptic feedbacks on other neurons (Holt and Koch, 1999). Indeed, the distribution of the extracellular potential can elicit transmembrane currents which may have decisive effects on the membrane potential of neighboring cells (Anastassiou et al., 2011; Buzsáki et al., 2012).

The goal of the current paper is to develop and implement an integrated three-dimensional model which synchronously captures both quantities, the membrane potential and the extracellular potential, during activity and which uses the neuron's geometry as it is instead of reducing it to cylindric compartments. The aim of such a model is to deepen the knowledge in signal processing and to carry out simulations on small networks of realistic neurons while having all these influences in action.

The work of Voßen et al. (2007) did a first step in the development of a generalized cable equation. It was built on the principle of the continuity of electric fluxes. Although the core model with the intra-, extracellular and membrane potential was correctly derived, the subsequent approach used to couple these unknowns and to solve them numerically resulted in major difficulties. The limit case to the standard cable equation evoked greater challenges and the simulations themselves were restricted to a very small time period of hundreds of micro seconds on a small part of a passive membrane.

The study of Xylouris et al. (2010) used a more direct approach for the coupling and generalized the existing model to active membranes. Nonetheless, although it was capable to reproduce action potentials, it still lacked in many characteristics of the signal processing, like the width of the propagating signal, the waveform of extracellular potential at activity and the computation on more complicated geometries. Indeed, computations on more complicated geometries diverged numerically. Furthermore, the membrane potential's defining equation in Xylouris et al. (2010) was the transmembrane current, which contains the time derivative of the membrane's capacitive property as only differential operator. The membrane potential's propagation was provided indirectly through the difference between the intra- and extracellular potential- thus making it actually hard to expect correct results for the spacial distribution. Moreover, as consequence, it produced vanishing transmembrane currents causing zero extracellular potentials and zero ephaptic interactions. This is why, the solving procedure with this direct coupling was of little use.

This paper introduces a completely new coupling of the unknowns. Therein, the defining equation for the membrane potential contains its own spacial differential operator. For the first time, we could carry out simulations on three-dimensionally resolved ideal neurons and on a small network of cells. This description, furthermore, allows for a proof that the extracellular potential distributes in the extracellular space like a current multipole. It will show that the only current monopole for a neuron exists at rest.

Model

Three-Dimensional Cable Equation

Let Ω_in and Ω_out be domains in ℝ³ denoting the neuron's intra- and extracellular space, respectively, and ${\bar{\Omega }}_{\text{in}}\cap {\bar{\Omega }}_{\text{out}}=\Gamma$ the membrane, a two dimensional manifold embedded in ℝ³. Let $\Omega ={\Omega }_{\text{in}}\cup {\Omega }_{\text{out}}={ℝ}^3$ be the whole space. Let, furthermore, Φ_in, Φ_out, and V_m be the intra-, extracellular, and membrane potential, respectively. Φ will represent either Φ_in or Φ_out.

The quantities σ_in and σ_out denote the intra- and extracellular conductivities, respectively. The normal n_{in → out} is the normal on the membrane Γ pointing from the intracellular space to the extracellular. We will need this quantities in order to define the fluxes. For the active transmembrane flux, we will just consider the Hodgkin–Huxley model for the sake of a simpler writing. There we have the sodium conductivity $g_{{\text{Na}}^{+}}$, the potassium conductivity $g_{{\text{K}}^{+}}$ and the leakage conductivity g_L. The quantities $E_{{\text{Na}}^{+}}$, $E_{{\text{K}}^{+}}$, and E_L denote the reversal potentials of the indexed ions. The gating parameters n, m, h obey ordinary differential equations (Hodgkin and Huxley, 1952) and calibrate how much of the maximal possible ionic flux passes through the channel.

Considering the non-membrane conductivity (≈ 3$\frac{mS}{cm}$) (López-Aguado et al., 2001) and the dielectricity of water (≈1), Gary Holt demonstrated in his Ph.D. Thesis (Holt, 1997) that a possible non-membrane capacitor would discharge with a time constant of approximately 3 ns. Because this time scale is much faster than the one of the phenomena considered—the fast channel dynamics react on a μs-time scale—, it appears as good approximation to assume no capacitive properties for the non-membrane spaces (ρ = 0 in Ω_in and Ω_out). Indeed, this is the basis of the derivation for the three dimensional cable equation. In addition, we will assume to have time invariant magnetic fields ($\frac{d\overrightarrow{B}}{dt}=0$). Then, Gauß's and Faraday's law satisfy root equations in the intra- and extracellular space, so that the conservative electric field can be expressed with the aid of a potential gradient. Combining this gradient with Gauß's law immediately leads to the Laplace equation for the potentials in the non-membrane spaces.

$\begin{array}{llll}\hfill \nabla ·\overrightarrow{E}& =\frac{\rho }{ϵ{ϵ}_0}\stackrel{(\rho =0)}{=}0,& & \left(1\right)\end{array}$

$\begin{array}{llll}\hfill \nabla \times \overrightarrow{E}& =-\frac{d\overrightarrow{B}}{dt}\stackrel{!}{=}0,& & \left(2\right)\end{array}$

$\begin{array}{llll}\hfill \Rightarrow \overrightarrow{E}& =-\nabla \Phi ,& & \left(3\right)\end{array}$

$\begin{array}{llll}\hfill \Rightarrow -\Delta \Phi & =0.& & \left(4\right)\end{array}$

The constants ϵ₀ and ϵ are the dielectricities in vacuum and material, respectively.

Because of flux continuity, the flux across the membrane is continuous and must correspond to the flux emerging from the membrane dynamics [denoted with j_all(V_m)]. Hence,

$\begin{array}{lll}\hfill -{\sigma }_{\text{in}}\nabla {\Phi }_{\text{in}}·n_{\text{in}\to \text{out}}=-{\sigma }_{\text{out}}\nabla {\Phi }_{\text{out}}·n_{\text{in}\to \text{out}}=j_{\text{all}}\text{ on}\Gamma .& & \left(5\right)\end{array}$

With this boundary condition in mind, we arrive at the three-dimensional cable equation (Figure 1):

$\begin{array}{llll}\hfill -\Delta {\Phi }_{\text{out}}& =0& \hfill \text{in}{\Omega }_{\text{out}},& \left(6\right)\end{array}$

$\begin{array}{llll}\hfill -\Delta {\Phi }_{\text{in}}& =0& \hfill \text{in}{\Omega }_{\text{in}},& \left(7\right)\end{array}$

$\begin{array}{llll}\hfill V_m& ={\Phi }_{\text{in}}-{\Phi }_{\text{out}}& \hfill \text{on}\Gamma .& \left(8\right)\end{array}$

The flux j_all contains all fluxes passing the membrane. Considering just the Hodgkin–Huxley model and some additional stimulus, it looks like:

$\begin{array}{l}{j}_{\text{all}}=c_m\frac{d{V}_m}{dt}+m^{3}h{g}_{{\text{Na}}^{+}}(V_m-E_{{\text{Na}}^{+}})+n^4{g}_{{\text{K}}^{+}}(V_m-E_{{\text{K}}^{+}})\\ +g_{\text{L}}(V_m-E_{\text{L}})+j_{\text{Stm}}.& \left(9\right)\end{array}$

Since it is possible to have different dynamics on each region of the neuronal membrane, we furthermore introduce the following δ-functions

$\begin{array}{lll}\hfill {\delta }_{\text{dend}}(x)& =\left\{\begin{array}{ll}1& \text{on the dendrite}\\ 0& \text{else}\end{array}\right\},& \\ \hfill {\delta }_{\text{active}}(x)& =\left\{\begin{array}{ll}1& \text{on the soma or nodes of Ranvier}\\ 0& \text{else}\end{array}\right\},& \\ \hfill {\delta }_{\text{syn}}(x)& =\left\{\begin{array}{ll}1& \text{on the postsynaptic density}\\ 0& \text{else}\end{array}\right\},& \\ \hfill {\delta }_{\text{stim}}(x)& =\left\{\begin{array}{ll}1& \text{on the stimulation area}\\ 0& \text{else}\end{array}\right\}.\end{array}$

With the help of these δ-functions, we can define a more refined transmembrane flux considering where it precisely occurs.

We define

$\begin{array}{cc}\begin{array}{l}{j}_{\text{HH}}(n,m,h,V_m)=m^{3}h{g}_{{\text{Na}}^{+}}(V_m-E_{{\text{Na}}^{+}})+n^4{g}_{{\text{K}}^{+}}(V_m-E_{{\text{K}}^{+}})\\ +g_{\text{L}}(V_m-E_{\text{L}}).\end{array}& (10)\end{array}$

The synaptic activity is simply modeled with the aid of a modified Heaviside function H(x,t). This function should be one as soon as the membrane potential at the pre-synapse exceeds a certain value, say 2 mV, and it remains one for the time the synapse is active regardless of the presynaptic membrane potential. Additional activation at the pre-synapse should integrated by the synaptic function α(V_m|_pre,t)

$\begin{array}{lll}\hfill j_{\text{syn}}(V_m{|}_{\text{pre}},t)=H(V_m{|}_{\text{pre}},t)·\alpha (V_m{|}_{\text{pre}},t),& & \left(11\right)\end{array}$

where V_m|_pre is the membrane potential at the presynaptic terminal.

Then the refined total transmembrane current has the form:

$\begin{array}{cc}\begin{array}{l}{j}_{\text{all}}(x,V_m)=c_m\frac{d{V}_m}{dt}+{\delta }_{\text{active}}(x)j_{\text{HH}}(n,m,h,V_m)\\ +{\delta }_{\text{stim}}(x)j_{\text{Stm}}(t)+{\delta }_{\text{syn}}(x)j_{\text{syn}}(V_m{|}_{\text{pre}},t).\end{array}& (12)\end{array}$

FIGURE 1

Figure 1. Compact scheme of the three-dimensional cable equation. Assuming pure resistivity for the non-membrane spaces, the respective potentials distribute according to the Laplace equation therein. These Laplace problems fulfill at the membrane an interface-condition which complies with the conservation of fluxes. The emerging flux from the potential equals the total transmembrane flux, denoted with j_all. Within j_all all transmembrane currents are accumulated: capacitive, channel, any stimulation or synaptic currents.

Numeric Model

The three dimensional cable equation (Equations 6–8) is a non-symmetric system (Φ_in does not couple with Φ_out the same way as Φ_out with Φ_in) of PDEs which couples two Laplace equations in the intra- and extracellular space with the transmembrane flux. This flux depends on the membrane potential. One difficulty in solving this system is the coupling of the membrane potential, which lives on a lower dimensional manifold, with the quantities, which live in full space. Since the discretization of this system is carried out with the help of integrals, the lower dimensional quantity cannot be measured the same way as the quantities in space (because the space integrals do not see it at all). In order to get rid of this particularity, we will extend the membrane potential, which is defined by the difference between the intra- and extracellular potential (V_m = Φ_in − Φ_out) on the membrane, to the intracellular space. To that end, we extend the extracellular potential to the intracellular space and combine its extension with the intracellular potential equation. So, we arrive at a problem for the membrane potential in the intracellular space.

Because V_m = Φ_in − Φ_out on the membrane Γ, we will extend Φ_out to the intracellular space continuously so that the following identity holds. Let this extension be denoted with ${\Phi }_{\text{out}}^{IN}$:

$\begin{array}{llll}\hfill V_m={\Phi }_{\text{in}}-{\Phi }_{\text{out}}={\Phi }_{\text{in}}-{\Phi }_{\text{out}}^{IN}& \hfill \text{on}\Gamma ,& & \left(13\right)\end{array}$

$\begin{array}{llll}\hfill \Rightarrow {\Phi }_{\text{out}}={\Phi }_{\text{out}}^{IN}& \hfill \text{on}\Gamma .& & \left(14\right)\end{array}$

At this point we have some freedom to choose the right hand side of the extracellular potential extension equation. We choose it to be zero. Then it can be easily combined with the intracellular problem (Equation 7), which is a Lapalcian, too. We have

$\begin{array}{rc}-\Delta {\Phi }_{\text{out}}^{IN}=0\text{ in}{\Omega }_{\text{in}},& \left(15\right)\end{array}$

$\begin{array}{ll}{\Phi }_{\text{out}}^{IN}={\Phi }_{\text{out}}\text{ on}\Gamma ,& \\ \Rightarrow -\Delta ({\Phi }_{\text{in}}-{\Phi }_{\text{out}}^{IN})=-\Delta V_m=0\text{ in}{\Omega }_{\text{in}},& \left(16\right)\end{array}$

$\begin{array}{l}-{\sigma }_{\text{in}}\nabla V_m·n_{\text{in}\to \text{out}}=j_{\text{all}}(V_m)\\ +{\sigma }_{\text{in}}\nabla {\Phi }_{\text{out}}^{IN}·n_{\text{in}\to \text{out}}& \\ \text{ on}\Gamma .& \end{array}$

Thus, instead of solving the system (Equations 6–8) we solve (Figure 2):

$\begin{array}{cc}-\Delta {\Phi }_{\text{out}}=0\text{ in}{\Omega }_{\text{out}},& (17)\end{array}$

$-{\sigma }_{\text{out}}\nabla {\Phi }_{\text{out}}·n_{\text{in}\to \text{out}}=j_{\text{all}}(V_m)\text{ on}\Gamma ,$

$\begin{array}{cc}-\Delta {\Phi }_{\text{out}}^{IN}=0\text{ in}{\Omega }_{\text{in}},& (18)\end{array}$

${\mathrm{ \Phi }}_{\text{out}}^{IN}={\Phi }_{\text{out}}\text{ on}\Gamma ,$

$\begin{array}{cc}-\Delta V_m=0\text{ in}{\Omega }_{\text{in}},& (19)\end{array}$

$\begin{array}{l} -{\sigma }_{\text{in}}\nabla V_m·n_{\text{in}\to \text{out}}=j_{\text{all}}(V_m)\\ +{\sigma }_{\text{in}}\nabla {\Phi }_{\text{out}}^{IN}·n_{\text{in}\to \text{out}}\\ \text{ on}\Gamma .\end{array}$

FIGURE 2

Figure 2. By extending the extracellular potential to the intracellular space (green) continuously, an extension of the membrane potential into the intracellular space is established. By means of this trick we obtain a coupling between the extracellular and the membrane potential which can be directly used for numerics and simulations.

For referencing reasons, we will call the additional current, which is considered in the boundary condition of the membrane potential equation (Equation 19), as ephaptic current

$\begin{array}{ll}\hfill j_{\text{eph}}:={\sigma }_{\text{in}}\nabla {\Phi }_{\text{out}}^{IN}·n_{\text{in}\to \text{out}}.& \left(20\right)\end{array}$

Numeric Discretization and Procedures

In space, we discretize this system (Equations 17–19) with the finite volume method (Versteeg and Malalasekera, 2007). This method guarantees the local conservation of fluxes. This is necessary, because the model has been derived on this principle. Furthermore, important characteristics of the solution, as we will see in the following section depend on this conservation. In time, an implicit method is used while the non-linearity is resolved with the Newton method.

Similarly to the finite element method, we discretize the domain Ω with volume elements, for example tetrahedrals, whose edge points and edges form the grid Ω_h, and we approximate the unknown functions (in our case V_m,Φ_out, and ${\Phi }_{\text{out}}^{IN}$) with a linear combination of shape functions. Our shape functions b_j(x) have the property to be continuous and linear on each elements (j = 0,…,#Ω_h = N). They are as many as our grid points (#Ω_h = N) and are uniquely determined by the following defining conditions

$\begin{array}{lll}\hfill b_j(x_k)& ={\delta }_{jk}{x}_j\in {\Omega }_h& \left(21\right)\end{array}$

$\begin{array}{lll}\hfill b_j(x)& \text{is continuous and linear on each element }& \left(22\right)\end{array}$

We represent our unknown functions with these

$\begin{array}{lll}\hfill V_m(x,t)& ={\displaystyle \sum _{j=0}^N}{v}_{m_j}^t{b}_j(x),& \left(23\right)\end{array}$

$\begin{array}{lll}\hfill {\Phi }_{\text{out}}(x,t)& ={\displaystyle \sum _{j=0}^N}{\varphi }_{{\text{out}}_j}^t{b}_j(x),& \left(24\right)\end{array}$

$\begin{array}{lll}\hfill {\Phi }_{\text{out}}^{IN}(x,t)& ={\displaystyle \sum _{j=0}^N}{\varphi }_{{\text{out}}_j}^{IN,t}{b}_j(x).& \left(25\right)\end{array}$

Purpose of the discretization schema is to establish linear systems out of the differential Equations (17–19) which uniquely determine the unknowns coefficients $v_{m_j}^t,{\varphi }_{{\text{out}}_j}^t,{\varphi }_{{\text{out}}_j}^{IN,t}$ of these linear combinations. The upper index t should indicate that these coefficients are time dependent.

For the finite volume method, we need to construct a so called dual grid, which arises from the domain discretization and which is used in order to discretize the differential space operators. We call the elements of the dual grid control volumes. The volume elements of the dual grid are defined by the edge points which correspond to the barycenters of the initial tetrahedrals and the barycenters of its sides and edges. By this construction, we create as many control volumes as we have nodes in the grid Ω_h. Let B_k be the control volume of the k-th grid node. We integrate the differential equations over this control volume and apply Gauß' integral theorem:

$\begin{array}{ll}\hfill -\Delta {\Phi }_{\text{out}}(x,t)=0& \left(26\right)\end{array}$

$\begin{array}{cc}\begin{array}{ll}{\displaystyle {\int }_{B_k}-}\Delta {\displaystyle \sum _{j=0}^N{\varphi }_{{\text{out}}_j}^t}{b}_j(x)dx\hfill & =-{\displaystyle {\int }_{B_k}{\displaystyle \sum _{j=0}^N{\varphi }_{{\text{out}}_j}^t}}\Delta b_j(x)dx\hfill \\ \hfill & ={\displaystyle {\int }_{\partial B_k}{\displaystyle \sum _{j=0}^N{\varphi }_{{\text{out}}_j}^t}}\nabla b_j(x)·\overrightarrow{n}(x)dS(x)\hfill \\ \hfill & ={\displaystyle \sum _{j=0}^N{\varphi }_{{\text{out}}_j}^t}{\displaystyle {\int }_{\partial B_k}\nabla }{b}_j(x)·\overrightarrow{n}(x)dS(x)\hfill \\ \hfill & ={\displaystyle \sum _{j=0}^N{\varphi }_{{\text{out}}_j}^t}{a}_{kj}.\hfill \end{array}& (27)\end{array}$

Because ∂B_k is a polyhedron and b_j(x) is analytically known, the integrals ${\int }_{\partial B_k}\nabla b_j(x)·\overrightarrow{n}(x)dS(x)=a_{kj}$ can be analytically computed. Furthermore, on the membrane these integrals equal to the transmembrane flux (Equation 12) which in general can also depend on other unknowns, like the gating variables or the membrane potential. Furthermore, this is the term which includes the time operator $\frac{d}{dt}$. We discretize our equation fully implicit and because this flux is not linear, we apply Newton's method to solve the emerging equations for each time step. Therein, the Jacobian of the system needs to be inverted, which we accomplish with high efficient iterative solvers. More precisely, we use a parallel ILU-preconditioned BiCGstab method (Barrett et al., 1987). All of this has been implemented with the use of the C++-library ug4 (Vogel et al., 2012), providing flexible numerical tools for these purposes.

Results

The intracellular problem (Equation 7) is a Laplace problem with a Neumann boundary. We referred this to the approximation of purely resistive non-membrane spaces (i.e., the intra- and extracellular space do not contain any free charges). Thus, the driving force of the intracellular potential is given by its Neumann-flux on the boundary (i.e., the membrane). Now, integrating the Laplace equation over the whole neuron and applying Gauß's theorem yields an important constrain for the transmembrane currents: The fluxes are balanced out over the whole membrane at each point of time!

$\begin{array}{cc}-\Delta {\Phi }_{\text{in}}=0& (28)\end{array}$

$\begin{array}{cc}\begin{array}{l}\Rightarrow {\displaystyle {\int }_{{\Omega }_{\text{in}}}-}\Delta {\Phi }_{\text{in}}dx={\displaystyle {\int }_{\Gamma }-}{\sigma }_{\text{in}}\nabla {\Phi }_{\text{in}}·n_{\text{in}\to \text{out}}dS(x)\\ ={\displaystyle {\int }_{\Gamma }{j}_{\text{all}}}(V_m)dS(x)\stackrel{!}{=}0\end{array}& (29)\end{array}$

There are at least two important implications of this situation. First, an influx at some point of the membrane, necessarily leads to an out-flux at some other point of the membrane with the same total amount of current. Moreover, this must happen simultaneously, since otherwise the condition is violated.

Second, the extracellular potential distributes like a multipole in the extracellular space.

Dipole-like Distribution of the Extracellular Potential for a Idealized Sphere Neuron

Regardless of the neuron's shape, the extracellular potential equation (Equation 17) demonstrates that its only source is the transmembrane flux as expressed through its boundary condition. A current monopole of the extracellular potential would be defined by the overall transmembrane flux. Yet, this flux is always zero as shown before (Equation 29). Thus, there is no monopole component and the extracellular potential distributes in space like a current multipole. To get some quantitative idea of its distribution, we approximate the neuron's geometry to a sphere. Then, we are able to express the extracellular potential with a generalized Fourier series of spherical harmonics.

Let Ω_in = B_R be a sphere with radius R and Γ = ∂B_R its boundary. The spherical harmonics $Y_l^m(\theta ,\varphi )$ satisfy the Laplace problem on this geometry:

$\begin{array}{cc}-\Delta Y_l^m=0& (30)\end{array}$

$\begin{array}{cc}{\Phi }_{\text{out}}(r,\theta ,\varphi )={\displaystyle \sum _{l\ge 0}{\displaystyle \sum _{m\ge -l}^l(b_{lm}{r}^{-(l+1)})}}{Y}_l^m(\theta ,\varphi )& (31)\end{array}$

$\begin{array}{cc}\Rightarrow -\Delta {\Phi }_{\text{out}}=0.& (32)\end{array}$

The solution Φ_out is concretized by the coefficients b_lm. These are determined by the transmembrane flux j_all(V_m):

$\begin{array}{ll}\hfill \frac{\partial {\Phi }_{\text{out}}}{\partial r}{|}_{r=R}={\displaystyle \sum _{l\ge 0}}{\displaystyle \sum _{m\ge -l}^l}-(l+1)\frac{1}{R^{l+2}}{b}_{lm}{Y}_l^m(\theta ,\varphi )=j_{\text{all}}(V_m)& \left(33\right)\end{array}$

$\begin{array}{ll}\hfill \Rightarrow b_{kn}=-\frac{R^{l+2}}{l+1}{\displaystyle {\int }_0^{\pi }}{\displaystyle {\int }_0^{2\pi }}sin(\theta )j_{\text{all}}(V_m)Y_k^n(\theta ,\varphi )d\theta d\varphi .& \left(34\right)\end{array}$

Especially, we obtain for the first coefficient b₀₀ which corresponds to the potential of a monopole:

$\begin{array}{cc}\begin{array}{ll}{b}_{00}\hfill & =-\frac{R^{l+2}}{l+1}{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_0^{2\pi }s}}in(\theta )j_{\text{all}}(V_m)\frac{1}{\sqrt{4\pi }}d\theta d\varphi \hfill \\ \hfill & =-\frac{R^{l+2}}{(l+1)\sqrt{4\pi }}{\displaystyle {\int }_{\Gamma }{j}_{\text{all}}}(V_m)dS(x)=0.\hfill \end{array}& (35)\end{array}$

Thus, the solution of the extracellular potential does not contain any monopole-part and behaves like a multipole falling in space with higher powers of the distance.

Numerical Error Analysis and Verification by a Comperison with NEURON

NEURON (Hines and Carnevale, 1997) is a highly sophisticated simulation environment for modeling a wide range of neuronal networks with the aid of the standard cable equation. Since the current three-dimensional model generalizes the one dimensional cable equation and since there are no non-trivial analytic solutions of an active neuron for our equations, we want to use this software environment in order verify both our model and our implementation. Our results should be very similar with these of NEURON for comparable computational domains. In order to keep the three-dimensional computation fast and in order to be able to create suitable three-dimensional computational domains, we carry out this comparison on a very long cylinder l = 9.9 mm with small diameter d = 200 μm in relation to its length ($\frac{d}{l}\approx 2·10^{-4}$). Such cases approximately comply with the assumption of the one-dimensional model (of infinite cylinders). No significant differences in the rise and propagation of an arising action should be visible.

We use proMesh (Reiter, 2014) to construct the three dimensional cylindric soma with a length 9.8 mm and a diameter 200 μm (Figure 3).

FIGURE 3

Figure 3. Computational domain constructed with proMesh (Reiter, 2014) for the comparison of the 3D-model's results with NEURON. The cylinder represents a soma having a length of 9.8 mm and a diameter of 200 μm. The purple is the intracellular and the blue the extracellular space. The domain is discretized with tetrahedrals.

This test domain we now use in order to first verify the the correct implementation of our discretization schema and second in order to see that we indeed obtain almost identical solutions in comparison with those produced by NEURON.

First is obtained, if the computed solution converges as the computational grid fineness is increased. In order to assess the second point, we have to compare the one dimensional solution of NEURON with the three-dimensional solution of our model. By construction of the one dimensional cable equation, each quantity, although computed on every point of a line, actually represents a volumetric quantity. Thus, the one-dimensional model assumes for all quantities to be radial symmetric and iso-potential on cross-sections of a three-dimensional cylinder. Considering this particularity, we can blow up the solution of NEURON to a three-dimensional solution and compare it with the solution of our model or we compare NEURON's solution with our solution recorded on the cylinder axis. For the sake of simplicity, we use the second way considering that its difference with the volumetric comparison is just the factor of the cross-section area.

Because for three dimensional numeric computations, domains have to be discretized, even simple cylinders never correspond to ideal cylinders, which, however, are the basis of the one-dimensional model. Thus, we will always expect small quantitative differences in such a comparison and, therefore, we are already satisfied to evaluate the differences with NEURON with the aid of an Euclidean integral norm

$\begin{array}{ll}\hfill \left|\right|f|{|}_{L^2(\left[a,b\right])}=\sqrt{{\displaystyle {\int }_a^b}|f{|}^{2}dx},& \left(36\right)\end{array}$

where the interval [a,b] corresponds to the time interval of the simulation. Furthermore, in order to get this measure dimensionless, we will consider the relative error between the solution of neuron V_mNEURON and the solution computed at refinement level x, denoted with V_{mLevel x}, over the interval [0,T]

$\begin{array}{ll}\hfill \frac{||V_{m_{\text{NEURON}}}-V_{m_{\text{Level}x}}|{|}_{L^2([0,T])}}{||V_{m_{\text{Level}x}}|{|}_{L^2([0,T])}}.& \left(37\right)\end{array}$

Yet, qualitative measures like propagation speed and signal width should be identical.

Concerning the numeric convergence at grid refinement, we computed the solution on our cylinder, composed by a tetrahedral grid, at two levels of refinement and observed the desired convergence (Figure 4). This behavior should serve as benchmark for the right implementation of the finite volume discretization schema.

FIGURE 4

Figure 4. Comparison of the three-dimensional model with NEURON. (A) Computational domain with the marked areas (B–D) where the membrane potential is recorded. (B–D) Time courses of the membrane potential at the corresponding areas. The solution of the three dimensional model at refinement level 0 is the blue line. After two refinements the solution converges—the red line representing the solution at refinement level 1 coincides with the solution represented at refinement level 2 (dotted green line). This implies the correct implementation of the applied finite volume discretization schema. We see that the solution produced by the three-dimensional model (dotted green line) is almost the same as the solution produced by NEURON (purple line). The small differences are due to the nature of the three-dimensional modeling procedure (see text).

The solution between the standard cable equation and the three dimensional model are qualitatively undistinguishable (Figure 3). The small numerical differences (Table 1) are due to the aforementioned reasons: the cylinder in the computation is a disretization of an ideal one, the cylinder's length is finite (the standard cable equation assumes infinite cylinders). Moreover, since the three-dimensional model additionally considers the coupling of the extracellular potential on the membrane, so that there are always to be expected some subtile differences in the solutions, which are reflected in Table 1.

TABLE 1

Table 1. Relative error of the computed solution in comparison with NEURON.

However as regards the emerging of the action potential (Table 1, Figure 4), the propagation speed of $5\frac{m}{s}$, and the signal width (Table 1, Figure 4) we receive identical results.

Simulation on a Small Network of Four Idealized Neurons

With a computationally quite demanding simulation, we also solve the Equations (17–19) on a more complicated geometry representing four idealized neurons with chemical synapses (Figure 5).

FIGURE 5

Figure 5. Computational domain of four idealized neurons consisting of a soma, an axon, dendrites and chemical synapses. Computationally, the synapses are modeled by a postsynaptic current at the postsynaptic site as soon the membrane potential exceeds some threshold indicating an action potential at the pre-synapse. Due to the complexity of the computation, a small time period of around 14 ms is covered.

The simulation is demanding, because we have a non-linear time-dependent domain problem in three dimensions. It means we solve several a huge linear systems in each time step within Newton's method. Thereby, the time step to be chosen is constrained by the fast dynamics of the active membrane's gating variables, which in our case is chosen with 10 μs, while we aim to simulate the time period of 14 ms. This means we need to compute the solution for 1400 time steps, which is time demanding despite parallel procedures due to the geometry's complexity.

We constructed the computational domain given by a small network of four neurons with the help of an algorithm developed in Niklas Antes' master thesis (Antes, 2009). Each cell consists of a myelinated axon (diameter d ≈ 5μm), a soma (d ≈ 20μm) and dendrites (d ≈ 10μm). The cells are several hundred micrometer separated among each other.

As regards the transmembrane current j_all(x,V_m) (Equation 12) for the different cell parts, we just considered passive properties on the dendrites while an active membrane reflecting Hodgkin–Huxley dynamics for the soma as well as for the nodes of Ranvier. On the myelinated sheaths, the transmembrane current j_all(x,V_m) is composed of the first term in Equation (12) only, the capacitive current. Furthermore, two of the cells (cell 1 and cell 4, see Figure 5) own external input areas by which the network can be stimulated.

Because we simulate the relatively small time period of 14 ms, we let the synapses work as pre-defined strong post-synaptic current pulses of some nA, which are triggered as soon as the membrane potential at the pre-synapse indicates that an action potential has arrived. This is assumed to happen when the membrane potential at the pre-synapse exceeds the value of 5 mV.

For the sake of simplicity, we choose a constant intra- and extracellular conductivity ${\sigma }_{\text{in}}=2\frac{mS}{cm},{\sigma }_{\text{out}}=20\frac{mS}{cm}$.

We activate the network by stimulating cell number one (see Figure 5) with approximately 30 pA at each of its input areas over the whole simulation period of 14 ms. At the moment of 8ms, we then stimulate cell number four with a current pulse of approximately 0.5 nA over 20μs. Although this stimulation of the fourth cell is not enough to generate an action potential alone, within the regime of this network and with the ephaptic current activated (Equation 20), an action potential arises (see Figure 6). This demonstrates that ephaptic interactions can have a decisive effect as to whether a neuron fires.

FIGURE 6

Figure 6. Effect of ephaptic interactions. Two simulations on an idealized small network of four cells and idealized paradigm. One simulation (left column), in which the ephaptic current (Equation 20) is neglected and one simulation (right column) in which it is included. Both simulations are carried out with the same stimulation paradigm. An initial signal spreads through the network. Additionally around the moment of 8 ms, the forth cell (the upper cell of this network) is activated slightly with a current pulse so that it depolarizes just below the threshold for and action potential. Although the effects of ephaptic interactions are very small, we see that they can determine whether a neuron activates in particular circumstances.

The model integrates the impact of the extracellular potential into the signal processing. Though its impact is rather small, it still can have a significant effect when combined with the right stimulation at the right time. Action potentials can arise, which otherwise would not show up (Figure 6).

Discussion

The three-dimensional passive model of Voßen et al. (2007) has been extended to a model with active membrane dynamics and has been reformulated mathematically with the aid of an extension of the membrane potential into the intracellular space. This reformulation, for the first time, facilitated numeric simulations of neuronal activity on three-dimensionally resolved idealized neurons generalizing the one dimensional cable equation by fully incorporating the three-dimensional extension of the neurons' geometry and by automatically considering the extracellular potential's influence on the membrane. As shown, the latter influence -though it is quite small- in combination with additional stimulation at the right timing can lead to an action potential which otherwise would not have arisen.

For the sake of verifying the correct implementation of this model and because it should deliver similar results as the one-dimensional cable equation for the limit case of long and thin cylinders, we carried out a comparison with NEURON and obtained very good agreement between the two models.

Based on the assumption of charge-free non-membrane spaces -an assumption also used for the derivation of the standard cable equation-, we could provide strong theoretical evidence (to our knowledge for the first time) with the aid of the three-dimensional model that there aren't any current monopoles as the overall out-flux across the membrane balances out. A significant consequence of this behavior is that the leading term of the extracellular potential's multipole expansion vanishes so that it falls in space with higher powers of its distance to the transmembrane current source. In the work of Lindén et al. (2011), this very assumption has been applied for the extracellular potential in order to arrive at converging LFPs. The authors in Lindén et al. (2011) showed that a monopole behavior would lead to a diverging LFP.

We consider the ability to carry out realistic simulations with the cable equation on three-dimensionally resolved ideal neurons as important step and milestone on the way of refining and generalizing existing models for neuronal activity. This three dimensional model facilitates gaining a better understanding of all the processes involved in the signal processing, especially the influence of the extracellular potential activity on the membrane and the impact of the precise three-dimensional shape of the neuron's geometry. Concerning the ephaptic communication, it would be interesting to further investigate its influence on synchronous firing within networks. The latter point also seems to be very promising since lots of precise experimental geometric data are produced. Questions connecting function with geometry can be directly tackled with this model.

However, there is still a long way to go on this path, as the biggest challenge at the moment for our model is its computational demand. Further algorithmic and computational analysis needs to be invested in order to make applicable cutting edge solvers of linear systems arising from partial differential equations -like algebraic multi grid methods- on highly parallel machines, even on graphic card clusters. As next steps, we want to focus on these improvements.

On the other hand, the computational efficiency is a big advantage for standard one dimensional cable equation. Once we accomplished this efficiency for the three-dimensional model, there are still lots of interesting applications which we wish to address- especially concerning backward modeling with questions like which are the underlying network properties in order to reproduce a given a extracellular potential activity wave.

Furthermore, we see the need of a deeper theoretical analysis of this model with the purpose to provide a mathematical proof that it converges to the standard cable equation for the limit case of infinite cylinders and vanishing extracellular resistivity.

Our long-range purpose is to generalize this model with homogenization and multi-scale techniques so that to be able to simulate the activity of bigger clusters of neuronal networks while also considering the detail in processing on the small scale.

Realized steps on this path will be hopefully items of future publications.

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.

Acknowledgments

The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007–2013) under grant agreement number 650003 (Human Brain Project).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fncom.2015.00094

References

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

Antes, N. (2009). Ein Werkzeug zur Erzeugung von Oberächengeometrien von Neuronen. Master's Thesis, University of Heidelberg.

Barrett, R., Berry, M., Chan, T. F., Demmel, J., Donato, J., Dongarra, J., et al. (1987). Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, Vol. 43. Philadelphia, PA: Society for Industrial and Applied Mathematics.

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

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

Hines, M. L., and Carnevale, N. T. (1997). The neuron simulation environment. Neural Comput. 9, 1179–1209. doi: 10.1162/neco.1997.9.6.1179

PubMed Abstract | CrossRef Full Text | Google Scholar

Hodgkin, A. L., and Huxley, A. F. (1952). Propagation of electrical signals along giant nerve fibres. Proc. R. Soc. Lond. B 140, 177–183. doi: 10.1098/rspb.1952.0054

PubMed Abstract | CrossRef Full Text | Google Scholar

Holt, G. R. (1997). 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. doi: 10.1023/A:1008832702585

PubMed Abstract | CrossRef Full Text | Google Scholar

Lindén, H., Tetzlaff, T., Potjans, T. C., Pettersen, K. H., Grün, S., Diesmann, M., et al. (2011). Modeling the spatial reach of the LFP. Neuron 72, 859–872. doi: 10.1016/j.neuron.2011.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

López-Aguado, L., Ibarz, J., and Herreras, O. (2001). Activity-dependent changes of tissue resistivity in the ca1 region in vivo are layer-specific: modulation of evoked potentials. Neuroscience 108, 249. doi: 10.1016/S0306-4522(01)00417-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Rall, W. (1962). Theory of physiological properties of dendrites. Ann. N. Y. Acad. Sci. 96, 1071–1092. doi: 10.1111/j.1749-6632.1962.tb54120.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Rall, W. (1964). “Theoretical significance of dendritic trees for neuronal input-output relations,” in Neural Theory and Modeling, ed R. F. Reiss (Stanford, CA: Stanford University Press), 73–97.

Reiter, S. (2014). Effiziente Algorithmen und Datenstrukturen für die Realisierung von Adaptiven, Hierarchischen Gittern auf Massiv Parallelen Systemen. Ph.D. dissertation, Universität Frankfurt.

Scott, A. C. (1975). The electrophysics of a nerve fiber. Rev. Mod. Phys. 47, 487–533. doi: 10.1103/RevModPhys.47.487

CrossRef Full Text | Google Scholar

Versteeg, H. K., and Malalasekera, W. (2007). An Introduction to Computational Fluid Dynamics: The Finite Volume Method. Harlow: Prentice Hall.

Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. (2012). Ug 4: A novel flexible software system for simulating pde based models on high performance computers. Comput. Vis. Sci. 16, 165–179. doi: 10.1007/s00791-014-0232-9

CrossRef Full Text | Google Scholar

Voßen, C., Eberhard, J. P., and Wittum, G. (2007). Modeling and simulation for three-dimensional signal propagation in passive dendrites. Comput. Vis. Sci. 10, 107–121. doi: 10.1007/s00791-006-0036-7

CrossRef Full Text | Google Scholar

Xylouris, K., Queisser, G., and Wittum, G. (2010). A three-dimensional mathematical model of active signal processing in axons. Comput. Vis. Sci. 13, 409–418. doi: 10.1007/s00791-011-0155-7

CrossRef Full Text | Google Scholar