Impact Factor 1.821

The Frontiers in Neuroscience journal series is the 1st most cited in Neurosciences

Methods ARTICLE

Front. Comput. Neurosci., 17 July 2015 | https://doi.org/10.3389/fncom.2015.00094

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

  • Department of Simulation and Modeling, Faculty of Informatics, Goethe Center for Scientific Computing, Goethe University Frankfurt, Frankfurt am Main, Germany

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, 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 ℝ3 denoting the neuron's intra- and extracellular space, respectively, and Ω̄inΩ̄out=Γ the membrane, a two dimensional manifold embedded in ℝ3. Let Ω=ΩinΩout=3 be the whole space. Let, furthermore, Φin, Φout, and Vm 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 nin → 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 gNa+, the potassium conductivity gK+ and the leakage conductivity gL. The quantities ENa+, EK+, and EL 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 (≈ 3mScm) (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 (dBdt=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.

        ∇·E=ρϵϵ0=(ρ=0)0,        (1)
      ∇×E=-dBdt=!0,        (2)
          ⇒E=-Φ,        (3)
-ΔΦ=0.        (4)

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 jall(Vm)]. Hence,

-σinΦin·ninout=-σoutΦout·ninout=jall      onΓ.      (5)

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

-ΔΦout=0                                    inΩout,    (6)
 -ΔΦin=0                                     inΩin,    (7)
        Vm=Φin-Φout                     onΓ.    (8)

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

jall=cmdVmdt+m3hgNa+(VmENa+)+n4gK+(VmEK+)            +gL(VmEL)+jStm.    (9)

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

δdend(x)={1on the dendrite0else},    δactive(x)={1on the soma or nodes of Ranvier0else},    δsyn(x)={1on the postsynaptic density0else},    δstim(x)={1on the stimulation area0else}.

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

We define

jHH(n, m, h, Vm)=m3hgNa+(VmENa+)+n4gK+(VmEK+)                                     +gL(VmEL).    (10)

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 α(Vm|pre,t)

jsyn(Vm|pre,t)=H(Vm|pre,t)·α(Vm|pre,t),      (11)

where Vm|pre is the membrane potential at the presynaptic terminal.

Then the refined total transmembrane current has the form:

jall(x, Vm)=cmdVmdt+δactive(x)jHH(n, m, h, Vm)                           +δstim(x)jStm(t)+δsyn(x)jsyn(Vm|pre, t).    (12)
FIGURE 1
www.frontiersin.org

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 jall. Within jall 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 (Vm = Φ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 Vm = Φ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 ΦoutIN:

        Vm=Φin-Φout=Φin-ΦoutIN     onΓ,      (13)
Φout=ΦoutIN                                         onΓ.      (14)

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

                                            -ΔΦoutIN=0                       inΩin,    (15)
                                                   ΦoutIN=Φout                 onΓ,-Δ(Φin-ΦoutIN)=-ΔVm=0                         inΩin,    (16)
              -σinVm·ninout=jall(Vm)                                                          +σinΦoutIN·ninout                                                                              onΓ.

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

                             ΔΦout =0                                 inΩout,    (17)
σoutΦout · ninout =jall(Vm)                     onΓ, 
                             ΔΦoutIN =0                                 inΩin,    (18)
                                    ΦoutIN =Φout                          onΓ, 
                                ΔVm =0                                 inΩin,    (19)
       −σinVm · ninout=jall(Vm)                                                  +σinΦoutIN · ninout                                                                                       onΓ.
FIGURE 2
www.frontiersin.org

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

jeph:=σinΦoutIN·ninout.    (20)

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 Vmout, and ΦoutIN) with a linear combination of shape functions. Our shape functions bj(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

bj(xk)=δjkxjΩh      (21)
bj(x)is continuous and linear on each element      (22)

We represent our unknown functions with these

   Vm(x,t)=j=0Nvmjtbj(x),      (23)
Φout(x,t)=j=0Nϕoutjtbj(x),      (24)
ΦoutIN(x,t)=j=0NϕoutjIN,tbj(x).      (25)

Purpose of the discretization schema is to establish linear systems out of the differential Equations (17–19) which uniquely determine the unknowns coefficients vmjt,ϕoutjt,ϕoutjIN,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 Bk be the control volume of the k-th grid node. We integrate the differential equations over this control volume and apply Gauß' integral theorem:

-ΔΦout(x,t)=0      (26)
BkΔj=0Nϕoutjtbj(x)dx=Bkj=0NϕoutjtΔbj(x)dx=Bkj=0Nϕoutjtbj(x) · n(x)dS(x) =j=0NϕoutjtBkbj(x) · n(x)dS(x) =j=0Nϕoutjtakj.    (27)

Because ∂Bk is a polyhedron and bj(x) is analytically known, the integrals Bkbj(x)·n(x)dS(x)=akj 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 ddt. 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!

                  ΔΦin =0    (28)
ΩinΔΦindx=ΓσinΦin · ninoutdS(x)                               =Γjall(Vm)dS(x)=!0    (29)

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 = BR be a sphere with radius R and Γ = ∂BR its boundary. The spherical harmonics Ylm(θ,ϕ) satisfy the Laplace problem on this geometry:

         ΔYlm =0    (30)
Φout(r, θ, ϕ)=l0mll(blmr(l+1))Ylm(θ, ϕ)    (31)
ΔΦout  =0.    (32)

The solution Φout is concretized by the coefficients blm. These are determined by the transmembrane flux jall(Vm):

Φoutr|r=R=l0m-ll-(l+1)1Rl+2blmYlm(θ,ϕ)=jall(Vm)    (33)
          ⇒bkn=-Rl+2l+10π02πsin(θ)jall(Vm)Ykn(θ,ϕ)dθdϕ.    (34)

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

b00=Rl+2l+10π02πsin(θ)jall(Vm)14πdθdϕ =Rl+2(l+1)4πΓjall(Vm)dS(x)=0.    (35)

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 (dl2·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
www.frontiersin.org

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

||f||L2([a,b])=ab|f|2dx,    (36)

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 VmNEURON and the solution computed at refinement level x, denoted with VmLevel x, over the interval [0,T]

||VmNEURON-VmLevelx||L2([0,T])||VmLevelx||L2([0,T]).    (37)

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
www.frontiersin.org

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
www.frontiersin.org

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 5ms, 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
www.frontiersin.org

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 jall(x,Vm) (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 jall(x,Vm) 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=2mScm,σout=20mScm.

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
www.frontiersin.org

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

Keywords: models, theoretical, ephaptic coupling, dipole effect, detailed 3D-modeling, 3D-modeling, cable equation

Citation: Xylouris K and Wittum G (2015) A three-dimensional mathematical model for the signal propagation on a neuron's membrane. Front. Comput. Neurosci. 9:94. doi: 10.3389/fncom.2015.00094

Received: 17 March 2015; Accepted: 02 July 2015;
Published: 17 July 2015.

Edited by:

Tobias Alecio Mattei, Brain & Spine Center - InvisionHealth - Kenmore Mercy Hospital, USA

Reviewed by:

Ingo Bojak, University of Reading, UK
Le Wang, Boston University, USA
Xin Tian, Tianjin Medical University, China
Mikhail Katkov, Weizmann Institute of Science, Israel

Copyright © 2015 Xylouris and Wittum. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Konstantinos Xylouris, Department of Simulation and Modeling, Faculty of Informatics, Goethe Center for Scientific Computing, Goethe University Frankfurt, Kettenhofweg 139, Frankfurt am Main 60325, Germany, konstantinos.xylouris@gcsc.uni-frankfurt.de