Solvers for the cardiac bidomain equations

Abstract

The bidomain equations are widely used for the simulation of electrical activity in cardiac tissue. They are especially important for accurately modeling extracellular stimulation, as evidenced by their prediction of virtual electrode polarization before experimental verification. However, solution of the equations is computationally expensive due to the fine spatial and temporal discretization needed. This limits the size and duration of the problem which can be modeled. Regardless of the specific form into which they are cast, the computational bottleneck becomes the repeated solution of a large, linear system. The purpose of this review is to give an overview of the equations and the methods by which they have been solved. Of particular note are recent developments in multigrid methods, which have proven to be the most efficient.

Fig. 1

Transmembrane polarization induced by extracellular stimulation. A) Unipolar anodal stimulation of a thin 3D sheet with equal anisotropy ratios. Only changes in transmembrane voltage of one polarity (anodal hyperpolarization in this case) are observed. These decay quickly within a few space constants from the stimulation site, with the space constant being a function of fiber orientation (indicated by the arrow). The location of the stimulus is indicated by a battery pole. The ellipsoidal region of stimulus-induced hyperpolarization under the positive pole is referred to as the anode (AN). B) With unequal anisotropy ratios, the polarization pattern is strikingly different. Changes in transmembrane voltage of both polarities are observed. Compared to (A), the region of anodal hyperpolarization extends much farther away from the stimulus site in the direction transverse to the fiber orientation. Tear-drop shaped virtual cathodes (VCs) form along the fiber direction. This peculiar pattern of hyperpolarization, consisting of an AN and VCs is referred to as a “dogbone”. It should be noted that the setups used in A and B were identical except that in (A) the transverse interstitial conductivity, g_et, was reduced so that the ratio of longitudinal to transverse conductivity was the same in the intracellular and interstitial spaces. C) The transmembrane voltage pattern induced in rabbit ventricles by application of a 1 volt battery to the right ventricular free wall with the anode and cathode (CA) separated by 1 cm is shown. Battery terminals are indicated by the “+” and “−” signs. Virtual electrodes, both VCs and Virtual Anodes (VAs) are seen as areas of depolarization and hyperpolarization adjacent to the anode and cathode (CA) respectively. The model contained 862,000 extracellular nodes and 547,000 myocardial nodes with the membrane ionic current modelled by a modified Beeler-Reuter representation(Skouibine et al., 1999). The time instant shown is 50 ms after the application of the battery, by which time a depolarization has emanated from the site and propagated across the entire myocardium.

Fig. 2

Bidomain representation of cardiac tissue in 2D. Intra- and extracellular (interstitial) domains are represented by the gray and orange planes respectively. Within each domain, conductivities are anisotropic as indicated by the different resistances in each direction. Each point in the intracellular domain has a potential associated with it ɸ_i, and a corresponding point in the extracellular domain with the potential ɸ_e. The voltage between the points is denoted V_m, and the points are linked by transmembrane current flow, I_m. Each green cylinder represents a section of membrane which is described by a model-dependent nonlinear current-voltage relationship.

Fig. 3

A) MRI scan of a rabbit heart (B₀ =11.7 T, resolution ≈ 25 μm) shown in longitudinal view. B) The same image stack shown in axial cross-section. A substack, indicated by the red solid lines, was extracted containing a wedge of the left ventricular (LV) freewall. C) LV wedge substack after segmentation removing blood. D) A multimaterial finite element method (FEM) mesh was created from the substack by an octree-based mesh generation technique that generated boundary-fitted, locally refined, hex-dominant hybrid (i.e., consisting of hexahedra, prisms, pyramids and tetrahdera) FEM meshes. The myocardial volume in (D) was discretized at an average resolution of 100 μm. The surfaces of the model are smooth (detailed inset) to avoid artifactual currents when large extracellular stimuli are applied. E) The LV wedge is immersed in a cuboid bath (1.13×1.93×1.0 cm³). The element size in the bath increased with distance from the myocardium allowing a significant reduction in overall degrees of freedom for the bidomain problem. This is particularly important when large volume conductors surround the heart (for instance, when a torso has to be modelled.) F) To demonstrate the computational feasibility of the approach, the wedge was paced by applying an external electric field of 6 V/cm via two plate electrodes located at the bath boundaries next to endocardial and epicardial surfaces. Extracellular potential, Φ_e, in the bath and transmembrane voltage, V_m, in the wedge are shown 5 ms after the end of stimulation. The model contained 1.31 million extracellular nodes and 1.02 million myocardial nodes. Simulating a pacing pulse with a basic cycle length of 300 ms took about 3 hours using 128 processors.

Fig. 4

One iteration of geometric multigrid PCG. To solve the system at each level, one iteration of PCG with ILU preconditioning was performed in parallel. For the coarsest grid, the system was solved on one machine by a direct solver. Prolongation operators, p_n, transfer values from a coarser grid to a finer grid, while restriction operators, $p_n^T$, transfer quantities from the fine to coarse grids.

Fig. 5

Setup used to benchmark sequential simulation runs: A sustained anatomical reentry was induced with an S1-S2 pacing protocol. Shown are polarization patterns of V_m for selected instants. Red arrows indicate the conduction pathways of the activation wavefronts. The left upper panel shows the bath geometry and the location of the reference electrode.