Anatomically accurate high resolution modeling of human whole heart electromechanics: A strongly scalable algebraic multigrid solver method for nonlinear deformation

Abstract

Electromechanical (EM) models of the heart have been used successfully to study fundamental mechanisms underlying a heart beat in health and disease. However, in all modeling studies reported so far numerous simplifications were made in terms of representing biophysical details of cellular function and its heterogeneity, gross anatomy and tissue microstructure, as well as the bidirectional coupling between electrophysiology (EP) and tissue distension. One limiting factor is the employed spatial discretization methods which are not sufficiently flexible to accommodate complex geometries or resolve heterogeneities, but, even more importantly, the limited efficiency of the prevailing solver techniques which are not sufficiently scalable to deal with the incurring increase in degrees of freedom (DOF) when modeling cardiac electromechanics at high spatio-temporal resolution. This study reports on the development of a novel methodology for solving the nonlinear equation of finite elasticity using human whole organ models of cardiac electromechanics, discretized at a high para-cellular resolution. Three patient-specific, anatomically accurate, whole heart EM models were reconstructed from magnetic resonance (MR) scans at resolutions of 220 μm, 440 μm and 880 μm, yielding meshes of approximately 184.6, 24.4 and 3.7 million tetrahedral elements and 95.9, 13.2 and 2.1 million displacement DOF, respectively. The same mesh was used for discretizing the governing equations of both electrophysiology (EP) and nonlinear elasticity. A novel algebraic multigrid (AMG) preconditioner for an iterative Krylov solver was developed to deal with the resulting computational load. The AMG preconditioner was designed under the primary objective of achieving favorable strong scaling characteristics for both setup and solution runtimes, as this is key for exploiting current high performance computing hardware. Benchmark results using the 220 μm, 440 μm and 880 μm meshes demonstrate efficient scaling up to 1024, 4096 and 8192 compute cores which allowed the simulation of a single heart beat in 44.3, 87.8 and 235.3 minutes, respectively. The efficiency of the method allows fast simulation cycles without compromising anatomical or biophysical detail.

Keywords: Algebraic Multigrid; Cardiac Electromechanics; Parallel Computing; Whole Heart Model.

Figure 1

Comparison of V_m, TnC_L and $\left[{\text{Ca}}_{\mathrm{i}}^{2+}\right]$ single cell traces between the GPB model (blue) and the GPB+LN model at steady state with λ = 1 (red).

Figure 2

Steady state traces of V_m, isometric tension, TnC_L and $\left[{\text{Ca}}_{\mathrm{i}}^{2+}\right]$ in a single cell pacing experiment under varying stretch ratios.

Figure 3

Mechanical boundary conditions: Homogeneous Dirichlet boundary conditions were enforced at the termini of the meshed pulmonary veins and venae cavae, and at the bottom of the soft material attached to the apex.

Figure 4

The CM (left), MM (middle) and FM (right) models of a four chamber heart geometry, discretized with average resolutions of 880 μm, 440 μm and 220 μm respectively. Insets illustrate level of detail in geometry resolution and smoothness of organ surfaces.

Figure 5

Electrical activation patterns for the CM, MM and FM models, shown from anterior and posterior views.

Figure 6

Electromechanical simulation of a human heartbeat. Top row shows displacement ∥u∥ and bottom row transmembrane voltage V_m at instants 0 ms, 40 ms, 80 ms, 150 ms and 350 ms after the stimulus delivered to the endocardia.

Figure 7

Comparison between different multigrid redistribution strategies used in the ptAMG preconditioner. The time needed to solve the linear systems in each newton step of the MM model during a 20 ms simulation was measured. The variable g denotes the group size of the AMG redistribution as described in Section 4.2. With g = 1 redistribution is turned off, while starting with g = 2, multigrid redistribution is activated.

Figure 8

Comparison of the average number of linear solver iterations (left panel) per linear solver step and the number of nonlinear Newton iterations (right panel) over a full heart beat for simulation runs using N_p = 512.

Figure 9

Strong scaling profile of an electromechanical model using a monodomain (left) or bidomain (right) EP model. Red, blue and green traces refer to CM, MM and FM simulation runs respectively, over a time frame of 500 ms. On the left hand side the assembly time for mechanics (33–34) and the solve time (i.e. preconditioner setup and solving) for all linearized systems of equations (32) is given. In addition, we show the total solving time of a monodomain electromechanical simulation. Depicted in the right plot is the computational time for all elliptic systems (38), the computational time for all ODE systems (39–40) and the computational time for all parabolic systems (42). In these three cases the computational time includes assembling, preconditioner setup and solving times. Solid lines correspond to the total solving time of a bidomain electromechanical simulation.

Figure 10

Effects of length dependence calcium binding affinity of troponin C (TnC). Shown are the distribution of intracellular calcium [Ca²⁺]_i and stretch ratio λ at t = 130 ms, and active stress S_a at t = 190 ms in the strong coupling (TnC(λ)) and weak coupling (TnC(λ_ref)) scenario.

Figure 11

Volume changes of the left ventricle, V_LV(t), in milliliters.