Computationally efficient simulation of electrical activity at cell membranes interacting with self-generated and externally imposed electric fields

The model represents the electrical activity of a biological cell inside a conductive medium (figure 1(A)). Space is segregated into an extracellular domain Ω_e and possibly multiple intracellular domains Ω_i (Ω = Ω_e∪Ω_i). The domains do not overlap (Ω_e∩Ω_i = ∅) and are each characterized by space-dependent conductivity tensors σ_e and σ_i. Unit normal vectors at domain boundaries point outward and are denoted n_{e, i}. Intra- and extracellular domains are separated by membranes Γ, assumed of zero thickness. Membranes accumulate charge due to capacitance and can harbor voltage-dependent ion channels. The voltage channels can follow nonlinear differential equations, for instance, Hodgkin–Huxley-like kinetics.

Zoom In Zoom Out Reset image size

For neural electric phenomena, a few simplifications in Maxwell's equations are possible. Typical spatial and temporal dimensions of neural current fluxes, together with the dielectric and magnetic parameters of biological media, render the feedback from the induced magnetic field onto the electric field negligible. Electromagnetic waves do not play an important role in the scales observed. The high polarizability of water, combined with the typical conductance of biological tissue, ensures that any free, unbalanced charge is balanced within fractions of a nanosecond. This is faster than the scale at which the electric processes of interest unfold.

These two relations justify the use of the quasi-static approximation of Maxwell's equations. In this approximation, the charge density in aqueous media is assumed to be zero [51, 52]. This however does not exclude the possibility of current volume density sources, and they are allowed in the extra- and intracellular spaces (ρ_{e, i}). In the quasi-static approximation, Poisson's equation governs the extra- and intracellular potentials Φ_{e, i}:

$$\begin{equation} -\nabla \cdot \sigma _{e}(x)\nabla \Phi _{e}(x,t) = \rho _{e}(x,t)\mbox{ in }\Omega _{e} \end{equation} \tag{ 1 }$$

$$\begin{equation} -\nabla \cdot \sigma _{i}(x)\nabla \Phi _{i}(x,t) = \rho _{i}(x,t)\mbox{ in }\Omega _{i}. \end{equation} \tag{ 2 }$$

In the absence of current sources, this simplifies to Laplace's equation

$$\begin{eqnarray*} &&-\nabla \cdot \sigma _{e}\nabla \Phi _{e} = -\nabla \cdot \sigma _{i}\nabla \Phi _{i} = 0. \end{eqnarray*}$$

On the exterior boundary ∂Ω_e, Dirichlet, Neumann or a mixture of both boundary conditions (BCs) can be applied (figures 1(B) and (C)):

$$\begin{equation} \Phi _{e} = \Phi _{D}(x,t)\mbox{ on }\partial \Omega _{D} \end{equation} \tag{ 3 }$$

$$\begin{equation} \sigma _{e}\nabla \Phi _{e}\cdot \mathbf {n}_{e} = I_{N}(x,t)\mbox{ on }\partial \Omega _{N}. \end{equation} \tag{ 4 }$$

If a grounding Dirichlet boundary is specified, then the solution is unique. When only Neumann BCs are specified, the problem is solvable up to an arbitrary constant if the total current through the boundary is conserved.

Membrane voltage is defined for any given point on the membrane as the potential difference for the same point in intra- and extracellular spaces:

$$\begin{equation} V_{m} = \Phi _{i}(x)-\Phi _{e}(x)\mbox{ on }\Gamma . \end{equation} \tag{ 5 }$$

Extracellular and intracellular current toward Γ is continuous for any small volume intersecting it and is denoted by the membrane current I_m:

$$\begin{equation} I_{m} = \mathbf {n}_{e}\cdot \sigma _{e}\nabla \Phi _{e} = -\mathbf {n}_{i}\cdot \sigma _{i}\nabla \Phi _{i}\mbox{ on }\Gamma . \end{equation} \tag{ 6 }$$

As current reaches one patch of membrane, it can either trigger the same amount of current being released in the opposite side or cross it through ionic flux. Total membrane current corresponds to the exchange of capacitive current plus the sum of transmembrane ionic currents (I_ion):

$$\begin{equation} I_{m} = C_{m}\frac{\partial V_{m}}{\partial t}+I_{\rm ion}, \end{equation} \tag{ 7 }$$

with C_m being the membrane capacitance per unit area. Ionic currents I_ion(V_m, u) can be determined by nonlinear ordinary differential equations depending on the membrane voltage and a vector of gating variables u_Γ at every node of Γ:

$$\begin{eqnarray*} &&\frac{\partial \mathbf {u}}{\partial t} = \mathbf {f}\left(V_{m},\mathbf {u_{\Gamma }}\right). \end{eqnarray*}$$

A simple passive membrane current can also be used:

$$\begin{eqnarray*} &&I_{\rm ion} = \frac{V_{m}}{R_{m}}, \end{eqnarray*}$$

with R_m being a constant membrane resistance. The main variables and parameters of the model are summarized in table A1.

To solve the complete space and membrane problem, the equations governing the potentials in the intra- and extracellular domains and the equations governing the time development of the membrane currents are solved alternatingly. Equations (1), (2), and (6) with BCs (3) and (4) are solved for a fixed moment in time with a fixed V_m using the FEM. A complete discretization and proof of existence of the solution for this spatial problem was presented by Lamichhane and Wohlmuth [53]. The FEM discretization is shown in appendix A.1.

Equation (7) is solved for a discrete time step after the solution of equations (1), (2) and (6) with the procedure specified in section A.1. At each time step, the discrete space equations (A.4) and (A.3) are solved for a fixed voltage vector V_m. The resulting I_m is used together with the vector of ionic currents per node $\mathbf {I}_{\rm ion}(\mathbf {V}_{m}^{\mathrm{n}})$ to find the next membrane voltage. The time iteration can be either explicit with the Euler method or implicit in the case of the CN method. The explicit scheme is easily obtainable with the previous values:

$$\begin{equation} \mathbf {V}_{m}^{\mathrm{n}+1} = \mathbf {V}_{m}^{\mathrm{n}}+\frac{\Delta t}{C_{m}}\left(\mathbf {I}_{m}^{\mathrm{n}}-\mathbf {I}_{\rm ion}^{\mathrm{n}}\right). \end{equation} \tag{ 8 }$$

The new $\mathbf {V}_{m}^{n+1}$ is then ready to be used again to solve the spatial equations.

The Euler scheme, however, is limited (see section 3) and an implicit CN was implemented. The CN scheme requires the next membrane current value $\mathbf {I}_{m}^{\mathrm{n}+1}$:

$$\begin{equation} \mathbf {V}_{m}^{\mathrm{n}+1} = \mathbf {V}_{m}^{\mathrm{n}}+\frac{\Delta t}{C_{m}}\left(\frac{1}{2}\left(\mathbf {I}_{m}^{\mathrm{n}+1}+\mathbf {I}_{m}^{\mathrm{n}}\right)-\mathbf {I}_{\rm ion}^{\mathrm{n}}\right). \end{equation} \tag{ 9 }$$

$\mathbf {I}_{m}^{\mathrm{n}+1}$ is an unknown future value of I_m and might prompt for an additional numerical strategy (e.g. a Newton solver). Still, through a careful re-arrangement of the terms, the same linear solver used to find Φ can be employed. Defining $c = \frac{\Delta t}{2C_{m}}$ and $\mathbf {g}^{\mathrm{n}} = \mathbf {V}_{m}^{\mathrm{n}}+\frac{\Delta t}{C_{m}}\left(\frac{1}{2}\mathbf {I}_{m}^{\mathrm{n}}-\mathbf {I}_{\rm ion}^{\mathrm{n}}\right)$, one obtains

$$\begin{eqnarray*} &&\mathbf {V}_{m}^{\mathrm{n}+1} = c\mathbf {I}_{m}^{\mathrm{n}+1}+\mathbf {g}^{\mathrm{n}}. \end{eqnarray*}$$

Combining the time iteration with the linear system as a vector in the advanced time step Mu^{k + 1} = b^{k + 1} results in

$$\begin{eqnarray*} \left[\begin{array}{@{}cc@{}}\mathbf {A} & \mathbf {B}\\ \mathbf {B}^{T} & \mathbf {0} \end{array}\right]\left[\begin{array}{@{}c@{}}\mathbf {\Phi }^{\mathrm{n}+1}\\ \mathbf {I}_{m}^{\mathrm{n}+1} \end{array}\right] = \left[\begin{array}{@{}c@{}}\mathbf {f}\\ \mathbf {G}\left(c\mathbf {I}_{m}^{\mathrm{n}+1}+\mathbf {g}^{\mathrm{n}}\right) \end{array}\right]. \end{eqnarray*}$$

After matrix manipulation and C = −cG, the new linear system is

$$\begin{equation} \left[\begin{array}{@{}cc@{}}\mathbf {A} & \mathbf {B}\\ \mathbf {B}^{T} & \mathbf {C} \end{array}\right]\left[\begin{array}{@{}c@{}}\mathbf {\Phi }^{\mathrm{n}+1}\\ \mathbf {I}_{m}^{\mathrm{n}+1} \end{array}\right] = \left[\begin{array}{@{}c@{}}\mathbf {f}\\ \mathbf {G}\mathbf {g}^{\mathrm{n}} \end{array}\right]. \end{equation} \tag{ 10 }$$

For the solution of (10), the system was factorized to eliminate I_m as done with system (A.11). However after the same factorization procedure [53], the left-hand side is non-symmetric (due to matrix C). The generalized minimal residual method was used instead of CG in this case.

Contrary to the Euler method, any solution to (10) depends on the previous I_m (implicit in gⁿ). This makes the method particularly weak at initialization and fast changing exterior stimuli. A better solution can be obtained by pre-calculating $\mathbf {I}_{m}^{\mathrm{n+\frac{1}{2}}}$ in a predictor–corrector scheme. This will be referred to as the Euler–Crank–Nicholson (ECN) method (see also figures 3(E) and (F)).

The software solver was written in C++ code, linked to the GNU-LGPL library CHASTE. CHASTE is being actively developed by the Computational Biology Group at the Oxford Computing Laboratory [50]. This framework provides functionalities for the solution of PDEs and ODEs for the heart and other tissues but does not provide cell membrane representations. To include the concept of insulating membranes, various additions to the CHASTE package were added using standard object-oriented programming techniques. These extended CHASTE's bidomain and ion channel functionalities, implementing the model, are presented in section 2. In its current version, the tool has been in development for 21 months. We named our tool CHASTE-Membrane and can be obtained for free on http://www.cs.ox.ac.uk/chaste/download.html under the 'Bolt-on Projects' section.

Two- and three-dimensional meshes were generated with the GNU-GPL tool Gmsh [54]. The neuron reconstruction presented in section 3 used NEURON's [55] Python's interface to parse the 1D geometry and Gmsh to generate the mesh. A custom-written Python code was used for the demarcation of the cell membranes and post-processing. Extra- and intracellular spaces are effectively independent meshes that only share node positions at the interfaces. Example meshes are provided in the software toolbox online. Paraview (Kitware Inc.) was used for visualization and sampled extraction of potentials over space. We used PETSc (Argonne National Lab.), integrated in CHASTE, to solve the linear systems produced by the FEM method. For the convergence simulations, an absolute tolerance value of 10⁻¹⁴ was used. For the simulations of figures 8 and 9, an absolute tolerance value of 10⁻⁶ was used. For the rest of the simulations, a relative tolerance value of 10⁻⁶ was used. All simulations, except the example in figures 8 and 9, were executed in a single processing core of an Intel® Core 2 Quad 2.83 GHz CPU, in a desktop computer with 4 GB RAM.

The solution scheme, introduced in the previous section, consists of two separate procedures: one solves equation (A.11) finding the electric potential Φ and the current I_m, under BCs on ∂Ω; the other adds up the membrane currents that result from the electric potential and ion fluxes, updating the membrane voltage V_m. The new V_m represents an updated condition on Γ. Given stationary BCs, if V_m over Γ corresponds to a steady-state value, the electric potential will correspond to the steady-state distribution of potential in Ω. In this stationary case, the implementation of the spatial solver (equation (A.11)) can be tested independent of time evolution (equations (8) and (9)).

Making use of this principle, it was assessed how the discretization of space influences the precision of the solution. For sufficiently fine spacing, the numerical solution should eventually converge toward the analytical solution. One of the cases where the analytical solution is known is for a 2D circular cell stimulated by a homogeneous field started at time t = 0 (see appendix section A.3. for this solution). Figures 2(A) and (B) show simulations for such a cell with diameter d = 15 µm and compare it to the analytical solution. The numerical precision was judged by comparison of the numerical solution ϕ_n and the analytic solution ϕ_a, measured by the normalized root mean square deviation (NRMSD):

$$\begin{eqnarray*} &&{\rm NRMSD} = {\sqrt{{\sum _{i = 1}^{n}(\phi _{n,i}-\phi _{a,i})^{2}}/{n}}} \bigg/{(\phi _{a,{\rm max}}-\phi _{a,{\rm min}}).} \end{eqnarray*}$$

As the typical grid spacing h went from 8 to 0.5 µm, the magnitude of the computed electric potential approached the analytic solution (figures 2(A) and (B)).

Zoom In Zoom Out Reset image size

The analytic solution assumes an infinite extracellular space, so the analytical and numerical results disagreed in the proximity of the domain boundaries ∂Ω. In a similar simulation, a distance from the membrane to the boundary of two times the cell diameter produces an error of 2.96% (figure 3(C)). At a distance of nine times the diameter, the error drops to negligible 0.69% (figure 3(C)).

Zoom In Zoom Out Reset image size

After confirming the accuracy of the spatial part of the numerical scheme, convergence was tested for the time development of the potential and membrane voltage as the time step was reduced. A configuration similar to that of figure 2 was used in figures 3(A) and (B). The extracellular field was imposed with Dirichlet BCs and was turned on at t = 0. For a mesh with a typical spacing of h = 1 µm, the numerical solution comes very close to the analytical solution after the transition phase. This is equivalent to the observation of figure 2. However, the initial phase of the V_m buildup is not captured properly (figure 3(A), left).

To achieve an accurate representation of the membrane voltage development, it is necessary to reduce the time step. The numerical solution converges to the analytic solution as the time step is reduced first to 5 ns and then to 0.5 ns (figures 3(A) and (B)). Because a smaller time step also allows for finer grid spacing, the grid was improved moderately. This had no influence on the convergence (data not shown). It was found empirically that the machine time per mesh node and per simulation time step (MTNT) was largely independent of the specific problem or the dimension. For the Euler forward solver used for figures 3(A) and (B), the MTNT amounts to 80–90 µs.

A similar MTNT (120 µs) was found for the next example: a 1 µm thin and 80 µm long cable under an axial field. The grid was generated with an intended grid spacing of 0.5 µm. A time step as small as 2 ns was necessary to achieve a stable solution. Due to the different morphology, the typical time of the V_m buildup is much larger than in the previous case. To compute a 300 ms time course, the simulation took 26 h (explicit Euler method). With many hours of computation time for a single 80 µm long cable, a timely treatment of realistic problems would not be possible with this approach. Using our insight into the different stability conditions for the Euler and CN solvers (see appendix), the latter was implemented. The CN solver tolerates orders of magnitude longer time steps keeping numerical stability (figure 3(D) CN 100 ns, CN 1 µs and appendix). CN and Euler solvers reach comparable precision when compared to a similar solution with the CEq in perfectly conducting extracellular space [15] (figure 3(D), dotted line). A combination of the computational steps that are used in the Euler and CN methods allows us to compensate the shortcomings of the CN alone (section 2). The resulting ECN solution is precise, even at times when the stimuli change, and allows for large time steps that would lead to numerical instability with the explicit Euler method alone (figure 3(E), ECN 100 ns). The instability of the explicit Euler method can be seen in the diverging zigzag lines in figures 3(E) and (F) (Euler 100 ns).

Stability of the explicit Euler scheme requires time steps in the order of hC_m/σ, where C_m/σ is typically ∼10 ns µm⁻¹. This means that a morphology with details in the 100 nm range demands 1 ns time steps (see appendix A.2.). The implicit CN scheme, on the other hand, is unconditionally stable for real-world parameters and can, theoretically, tolerate time steps close to the membrane time constant. Convergence results and simulation times for the different methods can be seen in table A2.

With the application of the CN solver to FEM problems with individual domains and insulating membranes, it seemed possible to treat realistic geometries and timescales within computation times of hours. This was put to a more rigorous test by the simulation of a 600 µm long and 0.5 µm thick axon-like fiber under stimulation by a strong field (1000 V m⁻¹). Cable endings are especially interesting for the effects of neuron stimulation because under spatially homogeneous conditions, endings are the places where the largest depolarization occurs [37, 56]. The exact shape of the endings and the influence of this shape on the magnitude of the induced depolarization have not been studied.

To show the potential of our approach to represent fine details in the morphology and analyze their effect, three different meshes representing differently shaped presynaptic boutons were created (figure 4(A)). The response of a finite neurite to an external field can be divided into two phases: a very rapid, initial depolarization (few microseconds) and a later equilibration of membrane currents (milliseconds). To represent both phases, a simulation duration of 8 ms was chosen. With a typical mesh spacing h = 0.25 µm (imposed by the cable radius), the explicit Euler method required a time step of 4 ns to avoid divergence due to numerical instability. The total computation time with the Euler forward method was eight weeks for cable I (figure 4(B)). The same precision was obtained more than 230 times faster, when the ECN solver was used with a time step of 1 µs, demonstrating the importance of this approach (figure 4(B)). Note that the numerical solution displays a small difference to the analytical solution in the steady-state membrane potential. This is likely an effect of the finite diameter of the 3D cable. Overall, the NRMSD was 1%.

Zoom In Zoom Out Reset image size

Although the ECN method was fast, to resolve the influence of the bouton shape on the depolarization in the initial phase, the time step of 1 µs was too large. For this purpose, the use of adaptive time steps is more appropriate. In this way, the fast response at the beginning of the simulation and the slow response at the steady state can be captured. With this approach, computation times can be reduced further (figures 4(C) and (D)). For this example, the points at which the time steps were switched had been chosen manually, guided by the magnitude of the membrane currents I_m. This approach is motivated by the fact that the change in the potential depends on I_m in equation (9). For large membrane currents, the time steps have to be chosen short in order to capture the rapid changes in the potential. This scheme can be extended into an automatically adapting time step: the values of I_m are available during the computation, as equation (A.12) is solved. Depending on the recent I_m values, the next time step can be chosen.

We found that the different bouton shapes influence the membrane potential V_m and the EP. The presence of a synaptic bouton increases the membrane depolarization during the initial phase, reflecting the rapid redistribution of charge in the bouton compartment. In the steady state, after several milliseconds, the larger surface of the bouton leads to a larger ionic current. This diminishes the depolarization, but increases the EPs (figure 4(E)).

The vast majority of studies of neuronal morphology is based on calculations that approximate the cell by cylindrical sections. The membrane potential in each section then develops according to the CEq. This approximation is very powerful and its application dominates the picture of electrical processes in excitable cells that we have today. Even the values for a basic parameter such as the cytosolic conductivity (σ_i) are obtained by matching measured voltage differences with voltage differences predicted from cable models [57]. The approximation of neurons as concatenations of cylindrical cables breaks down on the micrometer scale, for instance, when fine structures such as synaptic boutons or dendritic spines deviate from axial symmetry. While attempts can be made to account for the presence of extra membrane in the spines (by adjusting the surface capacitance and conductance obtaining valid expressions for the average voltage across a longer segment) on the microscopic scale, only a right representation of the morphological details can give correct simulations of the membrane voltages.

Even larger problems arise, when the relevant physical processes do not evolve within the axial/radial coordinates of the cable segments. To describe the interaction between a cylindrical segment and the extracellular field, only the field's component perpendicular to the cylinders faces is used to compute the CEq [38]. This works well as long as the approximation of a very thin cable is appropriate but it completely fails for spherical structures like cell bodies (figure 5). The key problem is the fact that only the area of the cylinder's faces determines the magnitude of the stimulus, while the area of the cylinder side determines the capacitance and conductance of the membrane. Face area and side area cannot at the same time match the effective areas of the sphere. This causes an erroneous result for the effect of extracellular stimulation of a cell body, even in the case that the field is oriented along the cylinders' axes (approximate 30% error in figure 5(B), 0°). When the orientation is changed, the approximation breaks down completely (figures 5(A) and (B)). This is a fundamental problem of the approximation by cylindrical segments and can only be solved when the membrane surface itself determines the orientation of the electric fluxes and not the axis of a segment.

Zoom In Zoom Out Reset image size

In the brain, neurons are not surrounded by empty, infinitely conductive space. Neurons and supporting cells are densely packed so that only about 20% of the brain volume is extracellular space [60]. Determining the influence of other cell bodies in the potential signal of the neuron is the key to understand the shape of extracellularly recorded action potentials (APs) and in general of local potentials. Collections of cells surrounding a neuron represent resistance and capacitance that have to be accounted for. The effect of this distribution in the extracellular space is still debated [61].

The tool presented in this work can be used to model collections of passive and active cells at close distances, resolving the EP during their activity. To illustrate the importance of representing other cell bodies, a model of a tight packing of 2D cells under stimulation was implemented (figure 6).

Zoom In Zoom Out Reset image size

Under a homogeneous electric field, a collection of cells in this configuration behaves essentially as a series of capacitors and resistors fed by a constant voltage source. When the voltage is turned on, current flows freely from one end of the circuit to another and the potential decays in the resistors (extra- and intracellular spaces). As charge builds up in the capacitors (the cell membranes), the current flow decreases and the membrane voltage increases. At the steady state, the potential difference across the entire equivalent 'circuit' is determined by the first and last membranes along the field direction. If the various capacitors and resistors have similar properties, the potential drops are distributed more or less equally across them (figure 6(D)). The main effect is that the membrane voltage of the cell targeted for stimulation is lower, compared to that of an isolated cell (figure 6(C)). For neuro-stimulation techniques such as transcranial magnetic stimulation that aims at supra-threshold activation of the cell and utilizes rapidly rising, strong stimuli as in this example, the shadowing effect and its consequences (e.g. secondary fields caused by the inhomogeneities of tissue [32]) can determine success or failure of stimulation. The formation of local secondary fields as visible in the potential gradients along the gaps between cells is not represented by activating function and the 'far-field' approach of the line-source method.

The case presented in figure 6 is equivalent to the transversal stimulation of a bundle of 'infinitely' long cells in 3D. In CEq solvers, such as NEURON, transversal stimulation exerts no effect on the membrane voltage, just as shown in figure 5. As the FEM does not impose any symmetry assumptions, it can represent the distortions of the EP caused by even more local obstructions, like a single small spherical cell body next to an active cell. The example in figure 7 demonstrates the influence of neighboring cells on the EPs during an AP. The slightly larger sphere in the center represents a neuronal soma that would receive a depolarizing current during the initial phase of an AP (in reality that would be the lateral current originating at the axon initial segment). This typically depolarizes the soma at a rate of up to 400 V s⁻¹. In the simulation, the current is supplied by an injection inside the cell with a current that varies over time, mimicking a real AP current shape. When the waveform of the resulting EPs between the active cell and the passive cell is compared to the EP at the opposite site, the amplifying influence of the adjacent membrane can be seen (figures 7(A) and (C)). At the peak amplitude, the potential is amplified by 8 µV (figure 7(D)). This difference is determined by the nearby cell whose own response and membrane voltage (figure 7(C), inset) distorts the potential in this region. Amplitudes of tens of microvolts can be typical of extracellular action potential (eAP) recordings at the given conductivity (figure 7(B)). With realistic complementary shapes of neighboring cells, the magnitude of the amplification would be much larger. The 8 µV effect on the potential is well resolved by the method although its magnitude is much smaller than the typical membrane voltage (tens of millivolts).

Zoom In Zoom Out Reset image size

Note that in this simulation, a current injection was given to the central cell. This illustrates the possibility to stimulate by current injection at arbitrary locations, which is a feature that would give the user the possibility to more precisely model stimulation via a pipette.

Comprehending effects of EPs and in general of the electrical functions of the cell requires realistic geometries and a realistic distribution of voltage-dependent channels. The construction of new stimulation and recording devices [25–28] such as those required to advance nerve–computer and brain–computer interfaces demands modeling environments that can model this heterogeneity. Simulation methods able to model biological tissue and artificial materials more realistically will be fundamental to lower development and experimentation costs.

To demonstrate the capabilities of the method and the tool to represent realistic setups, in figure 8, simulation of a cultured neuron equipped with voltage-gated ion channels on top of a glass coverslip is presented. This preparation is interesting for the development of multi-unit micro electrode arrays [62, 20].

Zoom In Zoom Out Reset image size

The morphology presented in figure 8 was generated from parts of the reconstruction of cell D151 from [30] (see section 2). The bottom of the cell was flattened to represent the interface with the glass; the cell and the glass were separated by a distance of 0.8 µm. The complete mesh consisted of 21 732 nodes and 86 716 tetrahedra in a domain of dimensions 340×200×60 µm³. An adaptive mesh was used, with refinements near and inside the cell down to 300 nm. A heterogeneous distribution of sodium and potassium voltage-dependent channels was set along the cell. The sodium conductances for the dendrite, soma/hillock, initial segment and the rest of the axon were, respectively, 10, 20, 120 and 20 mS cm⁻². The potassium conductances, for the same sections, were 10, 20, 30 and 20 mS cm⁻². Mammal kinetics at 23° were used [63].

The simulation of 8 ms of activity in this setup with the reference software required only 45 min on a single processor core. The solver was configured to use the ECN method with a fixed time step of 40 µs. A longer 100 ms simulation was executed with the same parameters producing a spike train (figure 8(E)). This required 6 h 27 min with a time step Δt = 100 µs. Computation time could be reduced to 4 h 16 min when the computation was distributed to two processors, and to 3 h 21 min in four processors. However, the mesh was not optimally partitioned (see section 4). A video of this simulation is provided online as the supplementary material, available from stacks.iop.org/JNE/10/026019/mmedia; there, only every second time step is displayed.

When the extracellular space is altered by the presence of other cells, the magnitude of EPs changes drastically. To demonstrate this, the layout of the realistic configuration was extended to include a passive cell that covers the soma, leaving a 1 µm cleft. This does not change the waveform of the intracellular potential during the AP (not shown) but the EP excursion increases 18-fold (figure 8(F)). This solves a longstanding problem of unrealistically small eAPs that occur when the line-source approximation is used to compute AP-induced EPs (see section 4).

A number of experimental studies investigated the influence of weak external fields on the neuronal population activity [5, 4]. Strikingly, fields strengths as low as 1 V m⁻¹ or even 0.2 V m⁻¹ [1, 3] were able to alter the temporal structure of ongoing spike activity. A full investigation of possible mechanisms of the observed effects with FEM simulations is outside the scope of our study and some aspects, such as polarization of thin processes in cultured neurons, are well captured in cable simulators (see, for instance, [64]). There is however one aspect our tool is particularly well suited to investigate: the membrane potential changes in the somatic region. To this end, we used a more complete reconstruction of the neuron from above (figure 8, cell D151 from [30]) and exposed it to a stationary field of 1 V m⁻¹. In figure 9, the resulting change of the resting membrane voltage of the entire neuron is shown in the left panels; the region around the soma is magnified on the right. Given the small magnitude of the induced field, the voltage changes about 400 µV and one can assume linearity. Fields of opposite orientation would cause effects of opposite sign. As the reconstructed neuron has a rather symmetric layout of the dendritic tree and the axon emerges from the base of the soma, only a week polarization of the axon hillock occurs. This would change for asymmetric configurations, for instance, for axons emerging from a primary basal dendrite.

Zoom In Zoom Out Reset image size

At any given moment, the membrane voltage V_m is assumed fixed and the system of partial differential equations is solved for variables I_m and Φ_{e, i} (denoted Φ when the domain is indifferent). The classic weak FEM formulation is obtained by multiplying equations (1) and (2) by test functions v_e, v_i, integrating and applying Green's first identity:

$$\begin{eqnarray} &&\fl\int \nolimits _{\Omega _{e}}\!\sigma _{e}\nabla \Phi _{e}\cdot \nabla v_{e}\,{\rm d}\mathbf {x}-\!\int \nolimits _{\Gamma }I_{m}v_{e}{\,\rm d}s =\! \int \nolimits _{\partial \Omega _{N}}I_{N}v_{e}{\,\rm d}s+\!\!\int \nolimits _{\Omega _{e}}\!\rho _{e}v_{e}{\,\rm d}\mathbf {x},\nonumber\\ \end{eqnarray} \tag{ A.1 }$$

$$\begin{eqnarray} &&\fl\int \nolimits _{\Omega _{i}}\sigma _{i}\nabla \Phi _{i}\cdot \nabla v_{i}{\,\rm d}\mathbf {x}+\int \nolimits _{\Gamma }I_{m}v_{i}{\,\rm d}s = \int \nolimits _{\Omega _{i}}\rho _{i}v_{i}{\,\rm d}\mathbf {x}. \end{eqnarray} \tag{ A.2 }$$

In equation (A.2), the assumption that no cell is open to the exterior boundary is made. The interface condition on Γ is comparable to a Neumann condition, with the difference that in this model, I_m is an unknown, sub-dimensional function on membrane space. Condition (5) is expressed in the weak form as

$$\begin{equation} \int \nolimits _{\Gamma }(\Phi _{i}-\Phi _{e})u{\,\rm d}s = \int \nolimits _{\Gamma }V_{m}u{\,\rm d}s, \end{equation} \tag{ A.3 }$$

with u being a test function also restricted to membrane space. Equation (A.3) can be understood as a special Dirichlet BC where the values must comply with a fixed potential difference V_m. The complete system of spatial equations is formed by (A.3) and the sum of equations (A.1) and (A.2):

$$\begin{eqnarray} &&\fl\int \nolimits _{\Omega }\sigma \nabla \Phi \cdot \nabla v{\,\rm d}\mathbf {x}+\int \nolimits _{\Gamma }I_{m}[v]_{e,i}{\,\rm d}s = \int \nolimits _{\partial \Omega _{N}}I_{N}v{\,\rm d}s+\int \nolimits _{\Omega }\rho v{\,\rm d}\mathbf {x}. \nonumber\\ \end{eqnarray} \tag{ A.4 }$$

Here, σ, Φ, v represent the corresponding symbol in any of the non-overlapping domains Ω_{e, i}. In (A.4), the operator[ ]_{i, e} changes the sign of the operand according to the face Γ that is being evaluated (integration is performed in the intracellular and extracellular face of Γ).

The matrix form of the problem is obtained by choosing a set of piecewise linear basis functions $\lbrace v^{1},v^{2},\ldots ,v^{N_{\Omega }}\rbrace$ for each of the N_Ω nodes in a triangulation Ω^h, and $\lbrace u^{1},u^{2},\ldots ,u^{N_{\Gamma }}\rbrace$ for each of the N_Γ nodes in a triangulation Γ^h. The result is a system of M = N_Ω + N_Γ equations:

$$\begin{eqnarray} &&\fl\int \nolimits _{\Omega }\sigma \nabla \Phi \cdot \nabla v^{\mathrm{i}}{\,\rm d}\mathbf {x}\,{+}\,\int \nolimits _{\Gamma }I_{m}[v^{\mathrm{i}}]_{e,i}{\,\rm d}s \,{=}\, \int \nolimits _{\Omega }\rho v^{\mathrm{i}}{\,\rm d}\mathbf {x}\quad\! \mathrm{i} \,{=}\, 1,..,N_{\Omega },\hspace{-6pt} \nonumber\\ \end{eqnarray} \tag{ A.5 }$$

$$\begin{eqnarray} \fl\int \nolimits _{\Gamma }u^{\mathrm{j}}[\Phi ]_{e,i}{\,\rm d}s & = \int \nolimits _{\Gamma }u^{\mathrm{j}}V_{m}{\,\rm d}s\quad \mathrm{j} = \mathrm{1},..,N_{\Gamma }. \end{eqnarray} \tag{ A.6 }$$

The Neumann term $\int \nolimits _{\partial \Omega _{N}}I_{N}v{\,\rm d}s$ was dropped but can be added in any future step.

Let Φ^k, $I_{m}^{\mathrm{j}}$ and $V_{m}^{l}$ be the approximations of values Φ(x^k), I_m(x^j) and V_m(x^l) at triangulation nodes x^k ∈ Ω^h, x^j ∈ Γ^h and x^l ∈ Γ^h. The continuous functions can then be replaced by $\Phi \approx \sum _{\mathrm{k} = 1}^{N_{\Omega }}\Phi ^{\mathrm{k}}v^{\mathrm{k}}$, $I_{m}\approx \sum _{\mathrm{j} = 1}^{N_{\Gamma }}I_{m}^{\mathrm{j}}u^{\mathrm{j}}$ and $V_{m}\approx \sum _{\mathrm{l} = 1}^{N_{\Gamma }}V_{m}^{\mathrm{l}}u^{\mathrm{l}}$, and the gradient by $\nabla \Phi \approx \sum _{\mathrm{k} = 1}^{N_{\Omega }}\Phi ^{\mathrm{k}}\nabla v^{\mathrm{k}}.$

Equations (A.5) and (A.6) for a given $\mathit {i}$ and a given $\mathit {j}$ are now

$$\begin{eqnarray} &&\fl\sigma \sum _{\mathrm{k}}\Phi ^{\mathrm{k}}\!\int \nolimits _{\Omega }\nabla v^{\mathrm{k}}\cdot \nabla v^{\mathrm{i}}{\,\rm d}\mathbf {x}+\sum _{\mathrm{j}}I_{m}^{\mathrm{j}}\int \nolimits _{\Gamma }u^{\mathrm{j}}[v^{\mathrm{i}}]_{i,e}{\,\rm d}s = \!\int \nolimits _{\Omega }\rho v^{\mathrm{i}}{\,\rm d}\mathbf {x},\nonumber\\ \end{eqnarray} \tag{ A.7 }$$

$$\begin{eqnarray} &&\fl\sum _{\mathrm{k}}\Phi ^{\mathrm{k}}\int \nolimits _{\Gamma }u^{\mathrm{j}}[v^{\mathrm{k}}]_{i,e}{\,\rm d}s = \sum _{\mathrm{l}}V_{m}^{\mathrm{l}}\int \nolimits _{\Gamma }u^{\mathrm{l}}u^{\mathrm{j}}{\,\rm d}s. \end{eqnarray} \tag{ A.8 }$$

Left-hand side entries of the matrix linear system form can then be defined as A^ik = σ∫_Ω∇v^k · ∇vⁱ dx and $B^{\mathrm{ij}} = \int \nolimits _{\Gamma }u^{\mathrm{j}}[v^{\mathrm{i}}]_{i,e}{\,\rm d}s$. Right-hand side vector entries can be defined as $f^{\mathrm{i}} = \int \nolimits _{\Omega }\rho v^{\mathrm{i}}{\,\rm d}\mathbf {x}$ and $G^{\mathrm{jl}} = \int \nolimits _{\Gamma }u^{\mathrm{l}}u^{\mathrm{j}}{\,\rm d}s$. The system of equations,

$$\begin{eqnarray} &&\sum _{\mathrm{k}}\Phi ^{\mathrm{k}}A^{\mathrm{ik}}+\sum _{\mathrm{j}}I_{m}^{\mathrm{l}}B^{\mathrm{ij}} = f^{\mathrm{i}} \end{eqnarray} \tag{ A.9 }$$

$$\begin{eqnarray} &&\sum _{\mathrm{k}}\Phi ^{\mathrm{k}}B^{\mathrm{jk}} = \sum _{\mathrm{l}}V_{m}^{\mathrm{l}}G^{\mathrm{jl}}, \end{eqnarray} \tag{ A.10 }$$

can be written in the matrix form as

$$\begin{equation} \mathbf {Mu} = \left[\begin{array}{@{}cc@{}}\mathbf {A} & \mathbf {B}\\ \mathbf {B}^{T} & \mathbf {0} \end{array}\right]\left[\begin{array}{@{}c@{}}\mathbf {\Phi }\\ \mathbf {I}_{m} \end{array}\right] = \left[\begin{array}{@{}c@{}}\mathbf {f}\\ \mathbf {G}\mathbf {V}_{m} \end{array}\right] = \mathbf {b}. \end{equation} \tag{ A.11 }$$

Matrix equation (A.11) has to be solved to obtain the potential Φ and membrane current I_m at any time step.

Although the Laplacian A matrix is positive definite, in general the matrix M is symmetric but indefinite. In [53], a factorization procedure is presented which can produce an equivalent positive-definite system after the elimination of vector I_m through static condensation. In that form, the conjugate gradient method can be used to solve the linear system. This form was used for the explicit Euler solver presented in the next sections. Although the details of the factorization in [53] will not be presented in this paper, we will still motivate the expression for I_m, as it is required for the stability analysis. This expression can be readily extracted from equation (A.11):

$$\begin{eqnarray*} &&\mathbf {BI}_{m} = \mathbf {f}-\mathbf {A\Phi }. \end{eqnarray*}$$

Due to symmetry between the membrane current expression and the voltage potential expressions (note that the I_m part of (A.9) and the Φ part of (A.10) produce the same matrix B), and a special node re-organization (nodes that touch the interface are moved to the bottom of the matrix), a left pseudo-inverse matrix W^T can be constructed for B such that

$$\begin{equation} \mathbf {I}_{m} = \mathbf {W}^{T}(\mathbf {f}-\mathbf {A\Phi }). \end{equation} \tag{ A.12 }$$

This expression depends only on Φ. Replaced back into equation (A.11), a linear system exclusively in terms of Φ can be obtained.

The explicit Euler form of the problem, although easily implementable, is only conditionally stable. As presented in section 2.3, a Crank–Nicolson (CN) scheme can be adapted to the solving algorithm providing excellent results (figures 3 and 4). In this section, the strict stability requirements for the Euler scheme and the more flexible CN requirements are shown.

The stability analysis for the general discretized version of the problem (equations (A.7), (A.8), (8) and (9)) is complex for two reasons: first, the piecewise nature of the problem resulting from the membrane, and second, the separation of the space and time equations in two steps. This results in a situation, where actually the difference between potentials, the membrane voltage, is the relevant parameter on the membrane, while in the classic stability analysis of the diffusion equation, it is the potential itself. The membrane is the place where local source terms, the membrane currents, are much more likely to cause instabilities than the solution of the Laplace or Poisson equation. Therefore, we restrict our study to this area and neglect possible variations in potentials of the non-membrane nodes.

The problem was analyzed for a regular 2D triangulation of typical element length h (figure A1(A)). At the membrane interface, an expression for I_m can be obtained and replaced in (8) and (9). Beginning with (A.7), for a membrane node x^j, assuming no external current sources, and u^j = 1, one obtains

$$\begin{equation} \sigma \sum _{\mathrm{k}}\Phi ^{\mathrm{k}}\int \nolimits _{\Omega }\nabla v^{\mathrm{k}}\cdot \nabla v^{\mathrm{j}}{\,\rm d}\mathbf {x}+I_{m}^{\mathrm{j}}\int \nolimits _{\Gamma }[v^{\mathrm{j}}]_{i,e}{\,\rm d}s = 0. \end{equation} \tag{ A.13 }$$

Note that u^j = 1 is a common assumption in FEM implementations of the Neumann BC. This assumption was also used in the membrane test functions of the computer solver. The membrane current at node x^j can be calculated from either the extracellular or intracellular potentials. With the adjacent extracellular nodes, this value is

$$\begin{equation} I_{m}^{\mathrm{j}} = -\sigma \frac{\sum _{\mathrm{k}}\Phi _{e}^{\mathrm{k}}\int \nolimits _{\Omega }\nabla v^{\mathrm{k}}\cdot \nabla v^{\mathrm{i}}{\,\rm d}\mathbf {x}}{-\int \nolimits _{\Gamma }v^{\mathrm{i}}{\,\rm d}s}. \end{equation} \tag{ A.14 }$$

In the two-dimensional triangulation (figure A1(A)), expansion and integration with the classic 'hat' test functions in (A.14) produce for the extra- and intracellular parts:

$$\begin{equation} I_{m}^{j} = \frac{\sigma }{h}\left(2\Phi _{e}^{\mathrm{j}}-\frac{1}{2}\Phi _{e}^{\mathrm{j}-1}-\frac{1}{2}\Phi _{e}^{\mathrm{j}+1}-\Phi _{e}^{\mathrm{p}}\right), \end{equation} \tag{ A.15 }$$

$$\begin{equation} -I_{m}^{j} = \frac{\sigma }{h}\left(2\Phi _{i}^{\mathrm{j}}-\frac{1}{2}\Phi _{i}^{\mathrm{j}-1}-\frac{1}{2}\Phi _{i}^{\mathrm{j}+1}-\Phi _{i}^{\mathrm{q}}\right). \end{equation} \tag{ A.16 }$$

Subtracting (A.16) from (A.15) and employing the definition of membrane voltage $V_{m}^{\mathrm{j}} = \Phi _{i}^{\mathrm{j}}-\Phi _{e}^{\mathrm{j}}$, one obtains

$$\begin{eqnarray} &&\fl I_{m}^{\mathrm{j}} = -\frac{\sigma }{h}\left(V_{m}^{\mathrm{j}}-\frac{1}{4}V_{m}^{\mathrm{j}-1}-\frac{1}{4}V_{m}^{\mathrm{j}+1}-\frac{1}{2}\Phi _{i}^{\mathrm{q}}+\frac{1}{2}\Phi _{e}^{\mathrm{p}}\right). \end{eqnarray} \tag{ A.17 }$$

As mentioned above at this point, the fluctuations of the exterior values $\Phi _{i}^{\mathrm{q}},\Phi _{e}^{\mathrm{p}}$ around the 'true' values are ignored. In this case, the standard ansatz of the von Neumann stability analysis can be used: a linear combination of periodic solutions of the form

$$\begin{equation} V_{m}^{\mathrm{j},\mathrm{n}} = \xi ^{\mathrm{n}}{\rm e}^{\mathtt {i}\gamma \mathrm{j}h}, \end{equation} \tag{ A.18 }$$

with j, n being the space and time indexes, $\mathtt {i}$ the imaginary unit and γ a variable wave number. The starting condition for V_m is on average zero. This implies that potentials in the extra- and intracellular regions of a region of membrane are on average the same. For simplicity, we choose this value to be zero. Under these conditions, potentials are allowed to fluctuate but only at larger scales than the fluctuations of V_m. Inserting (A.17) in the Euler scheme (8) and using a passive membrane resistance produces

$$\begin{eqnarray} &&V_{m}^{\mathrm{j},\mathrm{n}+1} = (1-\alpha -\beta )V_{m}^{\mathrm{j},\mathrm{n}}+\frac{\beta }{4}\left(V_{m}^{\mathrm{j}-1,\mathrm{n}}+V_{m}^{\mathrm{j}+1,\mathrm{n}}\right) \nonumber \\ &&\alpha = \frac{\Delta t}{C_{m}}\frac{1}{R_{m}},\quad \beta = \frac{\sigma \Delta t}{C_{m}h}. \end{eqnarray} \tag{ A.19 }$$

Replacing (A.18) into (A.19) leads to the amplification factor

$$\begin{eqnarray*} &&\xi _{E} = 1-\alpha -\frac{\beta }{2}-\beta \sin ^{2}\frac{\gamma h}{2}. \end{eqnarray*}$$

Stability is obtained under the condition |ξ_E| ⩽ 1, which is true for

$$\begin{eqnarray*} &&\Delta t\le \frac{2}{\frac{1}{R_{m}C_{m}}+\frac{3}{2}\frac{\sigma }{C_{m}h}}. \end{eqnarray*}$$

In cell biology, it is usual that C_mh/σ ≪ R_mC_m so the condition can be approximated by

$$\begin{equation} \Delta t<\frac{4C_{m}h}{3\sigma }. \end{equation} \tag{ A.20 }$$

This was tested in numerical experiments in 2D and 3D (figures A1(B)–(D)). The factor 4/3 is not exact as we dropped the contribution of the exterior values $\Phi _{i}^{\mathrm{q}},\Phi _{e}^{\mathrm{p}}$. Their influence should be on the order of the other potential-difference terms in (A.17) so that the result presented here is of the correct order.

Zoom In Zoom Out Reset image size

The restriction on the time step is laxer in the CN scheme. Inserting (A.17) into the CN scheme (9) expands to

$$\begin{eqnarray*} &&\fl\left(1+\frac{\beta }{2}\right)V_{m}^{\mathrm{j},\mathrm{n}+1}-\frac{\beta }{8}\left(V_{m}^{\mathrm{j}-1,\mathrm{n}+1}+V_{m}^{\mathrm{j}+1,\mathrm{n}+1}\right) \\ &&=\left(1-\frac{\beta }{2}-\alpha \right)V_{m}^{\mathrm{j},\mathrm{n}}+\frac{\beta }{8}\left(V_{m}^{\mathrm{j}-1,\mathrm{n}}+V_{m}^{\mathrm{j}+1,\mathrm{n}}\right) \end{eqnarray*}$$

The use of (A.18) produces the amplification factor

$$\begin{eqnarray*} &&\xi _{{\rm CN}} = \frac{1-\frac{\beta }{2}\left(\frac{1}{2}+\sin ^{2}\frac{\gamma h}{2}\right)-\alpha }{1+\frac{\beta }{2}\left(\frac{1}{2}+\sin ^{2}\frac{\gamma h}{2}\right)}. \end{eqnarray*}$$

Adding $-1+\frac{\beta }{2}\left(\frac{1}{2}+\sin ^{2}\frac{\gamma h}{2}\right)$ to the inequality −1 ⩽ ξ_CN ⩽ 1 gives

$$\begin{eqnarray*} &&-\beta \left(\frac{1}{2}+\sin ^{2}\frac{\gamma h}{2}\right)\le \alpha \le 2. \end{eqnarray*}$$

α and β are always positive, so the left inequality is always true. The right part establishes that the stability limits are set by the millisecond-scale membrane time constant τ_m = R_mC_m:

$$\begin{eqnarray*} &&\Delta t\le 2\tau _{m}. \end{eqnarray*}$$

The analytic solution for the homogeneous stimulus of a 2D cell was used to compute the reference solution in figures 2 and 3. The solution is known from [59, 18]. Potentials for a circular cell of diameter d under a homogeneous step field E starting at t = 0 are

$$\begin{eqnarray*} &&\Phi _{i}(r,\theta ,t) = -a(t)\cdot Er\cos \theta \; r<d/2\\ &&\ms \Phi _{e}(r,\theta ,t) = -Er\cos \theta -b(t)\cdot E\frac{d^{2}}{4r}\cos \theta \; r>d/2\\ &&\ms V_{m}(\theta ,t) = E{\,\rm d}\cos \theta (1-{\rm e}^{-\frac{t}{\tau }})(1-\epsilon ). \end{eqnarray*}$$

With

$$\begin{eqnarray*} &&a(t) = \frac{2\sigma _{e}}{\sigma _{i}+\sigma _{e}} \lbrace {\rm e}^{-\frac{t}{\tau }}+(1-{\rm e}^{-\frac{t}{\tau }})\epsilon \rbrace \\ &&b(t) = 1-\frac{2\sigma _{i}}{\sigma _{i}+\sigma _{e}}\lbrace {\rm e}^{-\frac{t}{\tau }}+(1-{\rm e}^{-\frac{t}{\tau }})\epsilon \rbrace \\ &&\tau = \frac{1}{\frac{1}{C_{m}R_{m}}+\frac{2\sigma _{i}\sigma _{e}}{C_{m}d(\sigma _{i}+\sigma _{e})}},\;\epsilon = \frac{\tau }{C_{m}R_{m}}. \end{eqnarray*}$$