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

In order to be able to examine the extracellular potential's influence on network activity and to better understand dipole properties of the extracellular potential, we present and analyze a three-dimensional formulation of the cable equation which facilitates numeric simulations. When the neuron's intra- and extracellular space is assumed to be purely resistive (i.e., no free charges), the balance law of electric fluxes leads to the Laplace equation for the distribution of the intra- and extracellular potential. Moreover, the flux across the neuron's membrane is continuous. This observation already delivers the three dimensional cable equation. The coupling of the intra- and extracellular potential over the membrane is not trivial. Here, we present a continuous extension of the extracellular potential to the intracellular space and combine the resulting equation with the intracellular problem. This approach makes the system numerically accessible. On the basis of the assumed pure resistive intra- and extracellular spaces, we conclude that a cell's out-flux balances out completely. As a consequence neurons do not own any current monopoles. We present a rigorous analysis with spherical harmonics for the extracellular potential by approximating the neuron's geometry to a sphere. Furthermore, we show with first numeric simulations on idealized circumstances that the extracellular potential can have a decisive effect on network activity through ephaptic interactions.


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(Rall, , 1964Scott, 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 R 3 denoting the neuron's intra-and extracellular space, respectively, and¯ in ∩¯ out = Ŵ the membrane, a two dimensional manifold embedded in R 3 . Let = in ∪ out = R 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 Na + , the potassium conductivity g K + and the leakage conductivity g L . The quantities E Na + , E 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 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 nonmembrane 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 ( d 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.
The constants ǫ 0 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, −σ in ∇ in · n in→out = −σ out ∇ out · n in→out = j all on Ŵ. (5) With this boundary condition in mind, we arrive at the threedimensional cable equation (Figure 1): The flux j all contains all fluxes passing the membrane. Considering just the Hodgkin-Huxley model and some additional stimulus, it looks like: Since it is possible to have different dynamics on each region of the neuronal membrane, we furthermore introduce the following δ-functions δ dend (x) = 1 on the dendrite 0 else , δ active (x) = 1 on the soma or nodes of Ranvier 0 else , δ syn (x) = 1 on the postsynaptic density 0 else , δ stim (x) = 1 on the stimulation area 0 else . With the help of these δ-functions, we can define a more refined transmembrane flux considering where it precisely occurs. We define 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 where V m | pre is the membrane potential at the presynaptic terminal.
Then the refined total transmembrane current has the form:

Numeric Model
The three dimensional cable equation (Equations 6-8) is a nonsymmetric 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 IN out : 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 Thus, instead of solving the system (Equations 6-8) we solve (Figure 2): 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 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.

Numeric Discretization and Procedures
In space, we discretize this system  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 IN out ) 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 is continuous and linear on each element (22) We represent our unknown functions with these Purpose of the discretization schema is to establish linear systems out of the differential Equations (17-19) which uniquely determine the unknowns coefficients v t m j , φ t out j , φ IN,t out j 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: Because ∂B k is a polyhedron and b j (x) is analytically known, the integrals ∂B k ∇b j (x) · 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 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!
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 m l (θ, φ) satisfy the Laplace problem on this geometry: The solution out is concretized by the coefficients b lm . These are determined by the transmembrane flux j all (V m ): Especially, we obtain for the first coefficient b 00 which corresponds to the potential of a monopole: 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 ( d l ≈ 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).
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 crosssection 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 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 m NEURON and the solution computed at refinement level x, denoted with V m Level x , over the interval [0, T] 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.
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.
However as regards the emerging of the action potential (Table 1, Figure 4), the propagation speed of 5 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).
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 σ in = 2 mS cm , σ out = 20 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.
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 onedimensional 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 threedimensional 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 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.
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.