Although, classical models represent the big picture of the electrochemical wave traversing cardiac tissue well, they may fail to reveal the finer details of cardiac conduction. For example, it is well-established that local perturbations to the conduction velocity may be arrhythmogenic; in particular, slowed conduction will increase the risk of arrhythmias [33]. It is therefore essential to understand how various remodelings of the heart affect the conduction velocity. Individual perturbations of the size and shape of the cardiac cells clearly affect the conduction velocity (e.g., [34]), but such changes are very hard to represent in a classical homogenized model, since a detailed representation of the individual cells in the tissue is needed. Furthermore, local density distributions of ion channels on the cell membrane will affect local conduction properties and such effects are also very hard, if even possible, to represent in the classical models.

The electrical conduction of the heart is believed to depend on direct cell-to-cell contact realized in terms of gap junctions (e.g., [35–37]). These connections are reduced under heart failure, resulting in impaired conduction velocity that may in turn increase the probability of arrhythmias (e.g., [38, 37]). However, even when conduction through gap junctions is significantly reduced, electrical signals are still conducted (e.g., [39]). This conduction is believed to rely on ephaptic coupling between neighboring cells via the extracellular space. The effect depends on the shape and size of the extracellular space and is thus not directly amenable to analysis via the homogenized bidomain or monodomain models.

The number of cardiomyocytes in the human ventricles can be estimated to be around 8 billion (e.g., [40]), and the number is close to 4 million for the mouse heart (e.g., [41]). In both cases, a homogenized model may be justified by the large number of cells involved. However, for experimental setups with monolayers of cardiac cells, the number of cells is much lower (hundreds or a few thousands) and the validity of the homogenized continuum approach becomes questionable. The EMI model, on the other hand, is very well suited, since it represents every individual cell. The ability to faithfully simulate monolayers of cardiac cells has become very important since it has become possible to simultaneously measure the transmembrane potential and the intracellular calcium concentration (e.g., [42]). Therefore, at least in principle, the inversion of spatial models of monolayers may be applied to characterize properties of single cells using monolayer experiments. This is particularly important because of the development of human induced pluripotent stem cells (hi-PSC). Based on skin samples, such cells can be used to derive cardiac cells with certain properties identical to a patient's cardiac cells. Therefore, this technology is believed to have great potential in the development of personalized drugs for rare diseases (e.g., [43–45]).

Twenty-five years ago, the best mathematical model of cardiac tissue was solved using 257 computational nodes [12, 46]. At that time, an accurate representation of cardiac tissue in terms of the representation of individual cells was inconceivable for reasons of both storage and computing time. This has changed dramatically; in recent computational studies, 29 million computational nodes were used to represent cardiac tissue [27, 47]. The computational mesh size in these simulations was about 59 μm, which should be compared with 100 μm, the typical length of a cardiac cell. This means that current simulators of the electrophysiology of the heart are, at least in principle, able to resolve features at the individual cell level.

The main purpose of the present paper is to assess the computational challenges of the EMI modeling framework. We will show how the model's complexity increases as the number of cardiac cells in the simulations increases and how the complexity of the membrane model affects the overall CPU demands.Furthermore, we will demonstrate that the EMI framework opens the possibility of simulating local properties of the cell that are hard to represent in homogenized models.

We introduce an operator splitting scheme for the EMI model and propose and compare three numerical schemes for the discretization of the resulting partial differential equations (PDEs): one finite difference-based (FDM) and two finite element-based (FEM) schemes of various degrees of complexity, computational cost, and accuracy. We compare the three schemes numerically in terms of convergence rates and computational cost. Moreover, to illustrate the distinctive capabilities of the EMI model, we present new results for simulating monolayers of cardiac cells with spatially heterogeneous distributions of ionic channels across the cell membrane.

Our results demonstrate that the EMI approach is computationally feasible: We can solve systems relevant for simulating monolayers of cardiac cells with sufficient resolution. Moreover, we show, using numerical computations based on the FDM code, that the computational effort per cell is bounded independently of the number of cardiac cells, and thus that the effort increases at most linearly with the number of cells.

In the next section, we will present the EMI model and three numerical methods used to solve the model. Next, we will discuss the numerical accuracy of the solutions, show convergence under mesh refinements, and assess the methods' CPU demands. To illustrate the ability to model local properties of individual cells, we present an example showing the difference in the conduction velocity of cells with uniform and non-uniform distributions of sodium channels. In the final sections, the results will be summarized and discussed.

We will use the EMI model to simulate collections of cardiac cells. However, to present the model, it is sufficient to consider the case of two coupled cells.

We assume that the complete computational domain consists of intracellular spaces Ωik, with k = 1, 2 in the case of two cells, that are connected by gap junctions Γ_1,2 and surrounded by a connected extracellular space Ω_e. The membrane is defined to be the intersection between each intracellular domain Ωik and the extracellular domain and is denoted by Γ_k, while the remaining boundary of the extracellular domain is denoted by ∂Ω_e. Figure 1 illustrates a two-dimensional (2D) version of this setup, showing two connected cells surrounded by extracellular space. In our computations (except in the first simple test case with an analytical solution) all cells are 3D and the cells can be connected in one-, two-, or three-dimensional collections. In one-dimensional strands of cells, the cell coupling is as illustrated in Figure 1; for two and three-dimensional collections of cells, the coupling in the y- and z-directions are similar to the x-coupling illustrated in the figure.

For the case illustrated in Figure 1, the EMI model can be formulated as follows: Find the extracellular potential u_e defined over Ω_e, the intracellular potentials uik defined over Ωik, and the transmembrane potentials v^k defined over Γ_k for k = 1, 2 and w defined over Γ_1,2 satisfying

for k = 1, 2. Here, n_e is the normal pointing out from Ω_e and nik is the (outward) normal pointing out from Ωik for k = 1, 2; σ_i and σ_e are the intracellular and extracellular conductivities, respectively; Iionk represents the ionic current density, which typically depends on additional state variables such as ionic concentrations; and I_gap represents the gap junction current density. In terms of units, the potentials u_e, uik, v^k, and w are given in mV; the current densities Imk, Iionk, I_1,2, and I_gap are given in μA/cm²; the conductivities σ_i and σ_e are given in mS/cm; the capacitances C_m and C_1,2 are given in μF/cm²; length is given in cm; and time is given in ms. In the following, we will refer to (1–9) as the EMI model. For brevity, we will write u_i in place of uik, v in place of v^k, I_ion in place of Iionk and Γ in place of Γ_k for k = 1, 2 when context allows.

In our computations, we will consider both a passive and an active model for the dynamics on the cell membrane between the intracellular and extracellular spaces. In the passive model, I_ion is given by the linear model

where R_m represents the resistance of the passive membrane (in kΩcm²) and v_rest denotes the resting potential of the membrane. In the active model, we let I_ion be represented by the action potential (AP) model of Grandi et al. [48]. In this case, Equation (6) is replaced by a system of the form

where v represents the membrane potential and s represents a collection of additional state variables introduced in the AP model. Furthermore, I_ion represents the sum of the ionic current densities across the membrane through a number of different types of ion channels, pumps, and exchangers and F(v, s) represents the ordinary differential equations (ODEs) describing the dynamics of the additional state variables. The Grandi model is implemented by defining a membrane potential v and a set of state variables s for each of the membrane nodes of the mesh. We let all state variables of the Grandi model, including the intracellular ionic concentrations, be defined only on the mesh nodes located on the cell membrane, and we allow the value of these variables to vary for different membrane nodes located on the same cell. The values of the state variables are updated in each time step using an operator splitting scheme described below. Intracellular and extracellular gradients of the ionic concentrations are ignored (see comment in section 4).

Finally, we represent the gap junction between neighboring cells by a passive membrane:

where R_gap represents the resistance of the passive membrane (in kΩcm²). A discussion of the modeling of the gap-junctions is given in Hogues et al. [1] where a boundary element method is used to solve a model similar to the system (1–9).

The ionic current density I_ion entering the EMI model through (6) typically introduce a significant number of additional states [e.g., as in (11)]. For this reason, we consider an operator splitting approach to solve the EMI model defined by (1–9).

The system (1–9) is solved by first applying given initial conditions for v and w. Then, for each time step n, we assume that the solutions vⁿ⁻¹ and wⁿ⁻¹ are known for t = t_n−1 on Γ and Γ_1,2 respectively. We then find the solutions at t = t_n using a two-step (first-order) operator splitting procedure, but note that a three-step (second-order) operator splitting could equally well be used (e.g., [11]).

In the first step, we update the solutions for the membrane potential by solving a system of ODEs of the form (11) and (12) over Γ with I_m set equal to zero. In the following numerical experiments, the ODE system (11) and (12) is solved by taking m forward Euler steps of size Δt^* = Δt/m for each global time step, though any other suitable ODE scheme could be used.

In the second (PDE) step of the operator splitting procedure, we solve the linear system arising from an implicit discretization in time and space of (1–9) with Iion1 and Iion2 set to zero. For the discretization in time of (6) and (9), we use an implicit Euler scheme using the solution from the first (ODE) step of the operator splitting scheme as the previous state.

When a linear model for I_ion is considered, the first (ODE) step of the splitting scheme is redundant and thus omitted, and Iionk for k = 1, 2 is kept in the PDE step, altering the linear system to be solved.

We propose and compare three different approaches for the spatial discretization of the PDE step in this study, each presented in the following sections. For the numerical experiments, the finite difference method (FDM) was implemented directly in MATLAB, while the finite element methods (FEMs) were implemented using the FEniCS finite element library [49, 50]. All computations were run on a Dell PowerEdge R430 with dual Intel Xeon processors (E5-2623 v4 2.60 GHz) and 12 x 32 GB RDIMM; each processor runs four kernels with two threads each.

A common problem in scientific computing is to solve a linear PDE defined on a certain geometry. After applying some sort of discretization characterized by a mesh parameter h, the remaining problem is to solve a linear system of algebraic equations. The linear solution process is usually said to be order optimal provided that the number of floating point operations (FLOPs) required to solve the problem grows linearly in the number of unknowns as h decreases. For self-adjoint, linear PDEs, optimal solvers are well understood (e.g., see the review papers [59, 60] for the theory of saddle point problems). In simulating cardiac tissue, optimal solvers exist for both the monodomain model and the bidomain model (e.g., [11, 61, 62]).

Feynman [63] suggested an alternative, but related, definition of order optimality: Suppose a numerical method is used to simulate a small space–time volume of a physical process and the mesh is refined to convergence. Then computational complexity should only grow linearly as the space–time volume is increased. For our application, this definition is very well suited; we consider a single cell surrounded by an extracellular space, and we carry out numerical simulations to find the mesh resolution in time and space necessary to obtain convergence. Then we define a numerical solution as being order optimal provided that the CPU efforts only increase linearly in the number of biological cells in the computation.

In the first unitless test problem we consider a 2D domain consisting of an extracellular domain and a single cell. In the remaining simulations, we consider 3D domains consisting of a number of connected cells and the surrounding extracellular space. The coupled cells are organized as a single layer where the cells are connected to each other in a grid in the x and y directions by gap junctions. The shape and size of the cells and the extracellular domain will be specified for each simulation below. We primarily consider cells of the shape illustrated in Figure 4, where each part of the intracellular domain, Ω_O, Ω_W, Ω_E, Ω_S, and Ω_N, is shaped as a rectangular cuboid.

The parameter values used in the simulations are given in Table 1 unless otherwise specified. Moreover, we use the initial condition w = 0 in all the simulations of connected cells. When the Grandi model is used to model I_ion, we mainly use the default initial conditions of the Grandi model for v and the remaining state variables. When a passive model is used for I_ion, we primarily use the initial condition v = v_rest.

As mentioned in the section 1, simulation of the electrophysiology of cardiac tissue is usually based on homogenized models such as the monodomain model or the bidomain model. The motivation for this is certainly that it requires considerably less computing power than the EMI approach considered here. Therefore, it is very important to understand the computational complexity of the EMI model to appreciate the applications in which this approach can be used.

The somewhat clunky geometry of the cells used above does not reflect reality very well. Indeed, cardiac cells have cylindrical shapes, but such shapes are inconvenient to address using FDMs, and we therefore apply the FEM. Figure 6 shows the membrane potential and surrounding extracellular potential for a simulation of two connected cylinders using the parameters given in Tables 1, 10. We note that the FEM is well suited for handling cylindrical geometry, and we expect that the method can also be used to handle the even more complex geometries that will arise when the T-tubules of ventricular cells (e.g., [71]) are incorporated in the model.

As mentioned in the section 1, it is difficult to represent a non-uniform distribution of ion channels along the cell membrane using classical homogenized models. This is an important shortcoming of the classical methods, because a non-uniform distribution of sodium channels is believed to affect the conduction velocity. In the EMI modeling framework, the representation of non-uniform distributions of ion channels is straightforward.

Figure 7 shows the solutions of two simulations of a collection of 30 × 5 cells with different distributions of the sodium channel conductance, g_Na. In the left panel the value of g_Na is constant on the entire membrane, while, in the right panel, g_Na is zero over a large part of the membrane and only non-zero on Ω_W and Ω_E, that is, near the ends of the cell in the x direction. The mean value of g_Na over the cells is the same in the two simulations. We observe that the conduction velocity is increased for the case with a varying value of g_Na compared to the case with a constant value. The conduction velocities reported in the figure are computed from the 10th and 20th cells in the third row in the y direction, and are defined as the distance between the cell centers divided by the time between each of the two cell centers reaches a membrane potential of v = 0 mV.

Figure 8 shows a more detailed view of the two cells at the center of the domain. Here, we observe that the conduction velocity across the first part of the cell appears to be higher for the case with a constant value of g_Na than for the varying case, but that the traveling wave moves faster across the gap junction for the case with a varying g_Na than for the constant case, leading to an overall increased conduction velocity for the varying case.

This effect is studied more closely in Figure 9, which shows the activation times and conduction velocities along the length of two cells in a similar pair of simulations. The gap junction between the two cells is located at x = 548μm, and we observe that there is a delay of about 0.1 ms between when the membrane on each side of the gap junction is activated. The delay appears to be slightly longer for the case with a constant g_Na compared to the varying case. We also observe that, overall, the wave uses less time to activate the two cells for the case with a varying g_Na than for the constant case, consistent with the results of Figures 7, 8. In addition, we observe that the shape of the activation curve is different in the two cases. In the case with a constant g_Na distribution, the activation curve becomes steeper toward the end of the cells, corresponding to a decreasing conduction velocity along the cell lengths. For the varying case, however, the activation curve flattens out toward the end of the cells, corresponding to an increased conduction velocity toward the cell ends.

The dynamics of the cell membrane are absolutely critical for the functioning of the cell and have been subject to intense studies for centuries. A wide variety of models are available through the open CellML library [72]. In our computations, we have used the ventricular cell model by Grandi et al. [48]. That model consists of a system of 39 ODEs defined on every computational node of the cell membrane and is believed to provide an accurate representation of the cell's action potential. From a computational point of view, we could have used numerous other models (e.g., [73–80]) with comparable complexity of the membrane computations. Common to all these models is that the ion currents are represented using models where the ion channel density must be specified. When the models are used as part of a monodomain or bidomain model, the channel density is most conveniently treated as a constant for each cell, but in the EMI model, the ion channel density associated with any of the currents can be specified as a function of space on the cell membrane. We noted in Figure 7 that a non-uniform distribution of sodium channels significantly affects the tissue's conduction velocity.

The numerical accuracy of the discretizations considered was assessed using a single-cell problem (24–30) and the method of manufactured solutions. As seen in Tables 2–4, all the discretizations provide converging numerical solutions. However, taking the L^∞ norms of the computed potentials for comparison, there are considerable differences in the convergence properties of the methods.

The convergence of the FDM discretization is linear and this method is the least accurate. The first order convergence of the FDM is to be expected, since all internal boundary conditions are approximated using first-order finite differences. Compared to the FDM, the mortar FEM yields solutions with considerably smaller error and the observed rates are about 1.5. It should be noted that, on a given grid, the methods lead to identical numbers of unknowns. The H(div) FEM is the most accurate among the methods considered with errors much smaller than those of the mortar FEM. As noted in section 3.3, the H(div) FEM leads to larger linear systems than the other two methods do (see Tables 7, 9). Finally, let us note that the manufactured solution employed was particularly simple and thus the numerical results obtained may not be universally valid.

The CPU efforts needed to solve the system (1–9) using the FDM or the FEM are given in Table 7 (FDM, passive membrane dynamics), Table 8 (FDM, membrane dynamics given by the Grandi model), and Table 9 (H(div) FEM, passive membrane dynamics). We observe that, for the FDM, the CPU efforts per cell seem to be bounded independently of the number of cardiac cells. This result implies that the solver is optimal in the sense defined above. When the Grandi model is used for the membrane dynamics, the CPU efforts per time step per cell are around 1.5 s. This enables us to simulate 16,384 cells, which defines a linear system with over 117 million unknowns. Since the CPU efforts per cell seem to be bounded independent of the number of cells, the CPU efforts will not explode as more cells are added to the computations. With proper parallelization strategies, it should be possible to simulate huge numbers of cells. In fact, the mouse heart, with around 4 million cells, may be within reach; with a computer that is 1,000 times faster (for large problems) than the one we used, it would require about one week to perform 100 time steps for 4 million cardiac cells.

We observe from Tables 8, 9 that the FDM method in the current implementation is significantly faster than the FEM code even though the FDM code is implemented in Matlab. It has proven to be difficult to derive optimal preconditioners to be used for the FEM, but we hope to be able to improve this part of the code in the future.

In the present report, we have used very simple geometries to represent the cells. We have assumed that the geometries are simple rectangular cuboids or have cylindrical shapes. However, real cardiac cells have much more complex geometries and future work will investigate the effect of the geometry on the solutions. Of particular importance is the effect of T-tubules in ventricular cells and how they change during illness (e.g., [71, 81]). The diameter of the T-tubules ranges from 20 to 450 nm (see [82]) and therefore a very fine computational mesh would be needed to represent these invaginations. Presently, we have run the FDM code with spatial resolution of 5 nm (in the case of only three coupled cells), so including T-tubules is within reach for the FDM code but not for the FEM code.

The focus of the present report has been to show that it is possible to simulate the electrical potential of cardiac cells based on explicit representation of the cells. We have focused on models that represent the membrane dynamics in terms of interchange over the cell membrane and we have ignored the spatial gradients of the ionic concentrations away from the cell membrane. Certainly, this is a major simplification; the intracellular concentration of calcium is essential and can be modeled using partial differential equations defined in the intracellular domain, see e.g. [83–86]. In Nivala et al. [83], a model based on Calcium released units (CRUs) is presented and the number of CRUs for a single cell used in the computations is typically ~20,000. In our model, a cardiac cell with a volume of 30 pL and a typical mesh length of 1 μm would consist of 30,000 computational blocks within each cell. A reasonable representation of the CRUs in the EMI model is therefore within reach.

As mentioned above, the conduction velocity is essential for the stability of the electrochemical wave underpinning the rhythmic contraction of the cardiac muscle. In numerical computations, the distribution of ion channels is usually assumed to be constant, but experimental evidence suggests that the ion channel density is non-uniform along the cell membrane. For instance, the density of sodium channels is believed to be much higher closer to the intercalated discs separating individual cells (e.g., [87]). The difference between uniform and non-uniform distributions of sodium channels was addressed in Figures 7–9. We observed that the conduction velocity was significantly lower in the case of a constant distribution of the sodium channels compared to the case of a non-uniform distribution. Interestingly, we also observed (Figure 9) that for the uniform case, the conduction velocity decreased along the cell, whereas it increased in the case of non-uniform distribution. Again, such effects are difficult to observe in the classical models (monodomain, bidomain). This effect deserves closer scrutiny and the EMI model provides a suitable framework for such studies.

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.