Electrical excitation and wave propagation in the myocardium were simulated as part of our electro-mechanical simulation framework (Gerach et al., 2021). The electrophysiological part is based on the verified framework acCELLerate (Seemann et al., 2010; Niederer et al., 2011). Tissue parameters were set such that conduction velocities (CV) of 800 mm/s along myocyte orientation and 550 mm/s in transverse directions were reached. Further information of the underlying monodomain formulation and numerical details are included in the Supplementary Materials (section 1.1.1).

For atrial myocytes, the cell model by Courtemanche et al. (1998) (CRN) was used to calculate the membrane currents, with alterations as described in van Oosterom and Jacquemet (2009), as the model would not be able to reach limit cycle (“steady state”) otherwise. Ventricular myocyte membrane currents were calculated using the O'Hara–Virág–Varró–Rudy (OVVR) model (O'Hara et al., 2011), with alterations to the sodium current as published by Dutta et al. (2017). The alterations were necessary to achieve physiological CV in tissue at the mesh resolution used in our model. Since the OVVR model has an epi and an endo variant of myocyte properties, the ventricular walls were subdivided into a layer where the endo variant (accounting for 20 % of the wall thickness), and a layer where the epi variant was assigned (remaining; cf. Figure 2). Further, an apico-basal heterogeneity in the slow delayed rectifying potassium channel conductance g_Ks was introduced, causing a shortening of action potential duration (APD) toward the apex by approximately 17%. An apico-basal APD gradient (with strong contribution of I_Ks) was observed by Szentadrassy et al. (2005) and shown by Keller et al. (2011) to result in physiological T-Wave morphology in the simulated ECG. Before initiating whole-heart simulations, all cell models were paced in a single cell environment until limit cycle was reached.

Excitation in the whole-heart simulations was initiated using a transmembrane stimulus current in the sinus node area of the atria and, then after assuming an atrio-ventricular node delay of 120 ms, the ventricles were stimulated using a pattern mimicking the activation via Purkinje–muscular junctions. This activation pattern was calculated using the optimization scheme published by Kahlmann et al. (2017), with the recorded ECG signals from Keller et al. (2012), i.e., from the same volunteer the geometry was based on, as a target. Here, the node density, maximal tree height, CV, and activation time of an activation sequence, mimicking the Purkinje network, are iteratively optimized. In each iteration, ventricular local activation times are estimated based on the current Purkinje tree using a fast marching algorithm. Combining these activation times with predefined action potential snippets yields an estimate of V_m distribution, from which an ECG is calculated using a lead-field approach. This ECG estimate is compared with the target ECG.

Based on the measured heart rate in the ECG, the stimulus cycle length was set to 850 ms, resulting in a heart rate of approximately 70 bpm. The instances of the cell models that describe the currents at each node of the geometry were initialized in an isolated and static mode. Thus, the impact of strong coupling was not present during initialization. Therefore, the electrophysiological side of the fully-coupled model took another 3 beats until it stabilized, determined by comparing the resulting ECG of two consecutive heart cycles (root mean square error < 0.1 mV), after being coupled to the mechanics model.

The forward problem of the ECG (i.e., the calculation of body surface potential maps (BSPM) from transmembrane voltages computed with the monodomain equation) was solved as published by Keller et al. (2012) with few minor alterations: Since excitation propagation was already calculated on a tetrahedral mesh, no further interpolation between the meshes of different types was needed. The deformed mesh of the heart was inserted into the torso mesh by remeshing the non-cardiac tissue layer surrounding the heart as described in section 2.1.

ECG leads were derived from the resulting BSPM as traces of the extracellular potentials over time at fixed locations on the body surface. Wilson's central terminal was used as reference potential for BSPM voltages. Numerical details on the solution of the forward problem can be found in the Supplementary Materials (section 1.1.2).

In order to calculate the force generated by the myocytes during contraction, the tension model by Rice et al. (2008) was coupled to the electrophysiological models. For all models, the temperature was set to a physiological value of 37°C. As the tension model itself has a stretch and stretch/shortening-velocity dependency, it was coupled to the mechanical part of the simulation, which is described later. Therefore, myocyte length changes and their derivatives (i.e., velocity) were interpolated using quadratic shape functions from the mechanical to the electrophysiological mesh for each time step (Gerach et al., 2021).

Parameters of the ventricular tension model were set as published by Timmermann et al. (2017) under consideration of the formulation for strong coupling. This means, that, additional to the previous stated coupling, the length and velocity dependent modulation of troponin C binding was calculated within the tension model, replacing the original formulation in the electrophysiological model. As the original formulation for troponin C binding in the electrophysiological model is neither length nor velocity dependent, this can lead to an alteration of the resulting AP, depending on contraction. From a physiological point of view, strong coupling is a more accurate representation of the underlying mechanisms. But as calcium dynamics in electrophysiological models are often quite delicate, such modification can disturb their balance and result in unphysiological behavior.

In terms of the atrial tension model, the original rat species parameters (Rice et al., 2008) were applied, as proposed by Brocklehurst et al. (2015). Yet, strong coupling led to instabilities in the underlying electrophysiological model (CRN), manifesting as the model being unresponsive to pacing after 2 fully coupled beats. Therefore, strong coupling was not applied to atrial cells.

Passive: To describe its passive stress-strain relationship, the cardiac tissue is assumed to be nearly incompressible, anisotropic, and hyperelastic. There are numerous published material laws describing these properties, with the one by Usyk et al. (2000) being used for the ventricles as well as atria in this work, with the parameters being set according to Gurev et al. (2015). The choice for this specific material law was based on its behavior during our proposed initialization scheme. Other materials model would have needed to be re-parameterized, similar to as proposed by Kovacheva et al. (2020). For all other types of tissue, the energy function was characterized by using a neo-Hookean material (Gerach et al., 2021). See Supplementary Materials (section 1.2.1) for the respective energy functions and used parameters.

Active: The stress calculated due to active tension (T_k) as described in 2.2.3 was linearly interpolated from the nodes of the electrophysiological mesh to the Gauss integration points of the quadratic elements of the mechanical mesh. For this, a constant mapping matrix between the two meshes was determined prior to any deformation. The tension was then added to the first element of the second Piola–Kirchhoff stress tensor (S), resulting in a contraction in main myocyte direction (see also Gerach et al., 2021).

Surrounding Tissue: As the heart is embedded in the pericardial sac, its influence on the motion of the heart has to be taken into account as well. As shown by Fritz et al. (2014) and Strocchi et al. (2020b) the inclusion of these influences is crucial, in order to be able to simulate a physiological contraction. Within our framework, this is done by modeling the interaction between the inner (i.e., heart-facing) surface of the non-cardiac tissue and the epicardial surface of the heart as a friction-less contact problem (Fritz et al., 2014; Pfaller et al., 2019; Gerach et al., 2021), yet different approaches, such as using spatially varying Robin boundary conditions, do exist (e.g., Strocchi et al., 2020b). For the friction-less contact problem, two surfaces based on the used geometry had to be defined: One being the heart-facing surface of the tissue surrounding the heart and the other being the epicardial surface of the heart (defined manually). As these surfaces are biologically equivalent to the pericardium and pleural surface, which form a lubricated laterally-movable contact surface, they will be hereinafter referred to as the pericardium. The gap function (g_N) governs the acting force on the respective aforementioned surfaces in normal direction. The formulation as defined by Fritz et al. (2014) in Equation (24), was replaced with a non-linear continuous function, as it resulted in a better convergence for our used geometry. A more extensive description, as well as comparison between the old and new gap function can be found in the Supplementary Materials (section 1.2.2).

Circulation: To model effects of the circulatory system on cardiac mechanical load, a fully coupled closed-loop lumped model representing the blood flow through the heart chambers, the lungs, and the rest of the body was used as described by Gerach et al. (2021). Therefore, for our given geometry four distinct closed cavities were defined (inner surfaces of each of the four chambers of the heart), with the openings in the atria being closed via triangulation without adding additional nodes. All other distinct parts of the circulatory system are concentrated and represented by a series of diodes, resistors, and capacitances. The finite element mechanics model and the lumped circulation model are then coupled with the goal of iterative adaption of the respective cavity pressures, such that the in- and outflows match one another. A more detailed description can be found in the Supplementary Materials (section 1.2.3).

In order to damp oscillations, the parameters of the Rayleigh damping matrix (α_RM+β_RK) as proposed by Fritz et al. (2014) in Equation (11) were adjusted. The coefficient α_R was set to 100 s⁻¹ (Fritz et al., 2014; Gerach et al., 2021) and β_R to 0.04 s. It should be noted that K (i.e., the isotropic tangential stiffness matrix of the passive stresses) and the resulting damping are highly dependent on the state of contraction and the used material law, hence, making a direct comparison to previously published parameter values difficult. The values chosen for damping here were selected such that oscillation introduced by the external stresses were damped sufficiently. Especially during isovolumetric relaxation, the problem to solve becomes rather stiff. Here the interactions during the iterative coupling of the mechanics and circulation problems can manifest as oscillations which, without damping, would build up and result in diverging solutions. The subsequent deformation was then calculated using a Newmark-beta scheme as described in the supplement of Gerach et al. (2021), where the resulting non-linear system was solved using distributed supernodal LU decomposition (Li et al., 1999; Li and Demmel, 2003), provided in the PETSc software package (Abhyankar et al., 2018), as the direct solver such that the residual norm fulfilled ||r|| < 10⁻⁸.

The simulation of one heart beat took approximately 20 h using all cores on a second generation 64-Core AMD EPYC processor (7H12), varying depending on if the time step had to be reduced momentarily during some part. But, this also included the time needed to save the results. As the deformation of both meshes needs to be saved, this can take up a significant amount of time—especially due to the higher resolution electrophysiological mesh.

Our 3D model of the heart was based on MRI data obtained in diastatic state, meaning in the presence of respective stresses. For initialization, the geometry of an unloaded state had to be determined. Thus, the goal was to define a mesh state which, when end-diastatic pressure is applied, converges toward the initially obtained mesh.

Different from Gerach et al. (2021), we used an approach aiming to find said unloaded state by applying a force to the first and fifth element of the second Piola–Kirchhoff stress tensor S. Looking at the previously defined myocyte coordinate system this meant adding a force in main (f) and in sheet (s) myocyte direction, thus inverting the stretching effect of the diastatic pressure:

Of course this requires that transmural and the sheet-normal myocyte direction align—if they do not, a rotation matrix would need to be defined. The applied force T_init was then iteratively optimized, so that the residual left ventricular reference volume (V₀) was compliant with the data from Klotz et al. (2006):

with V_diastatic being the volume of the segmented diastatic geometry, and P_diastatic the diastatic pressure). After the residual volume was reached, end-diastatic pressure was applied and the factor ε of interaction between the heart and its surrounding tissue was linearly increased within a time span of 50 ms from 0 kPa to 20 kPa, see Supplementary Figure S1. The ventricular volumes were then given another 250 ms to converge toward their end-diastatic volume.

By applying an initial unloading stress of 30 kPa according to Equation (2), the ventricular volume decreased from 142 to 146 mL (as segmented in the end-diastatic state) to 84 mL/89 mL (left/right), which is in concordance with Equation (3) (approximate target residual volume for the left ventricle: 78 mL; right ventricle: 81 mL). Succeeding the unloading stress, diastatic ventricular pressures (determined from previous simulations as convergence values) of 7.3 mmHg for the left side and 6.5 mmHg for the right side of the heart were applied. This resulted in initialized diastatic volumes of ≈152 mL (left ventricle) and ≈156 mL (right ventricle). No further adaptation of stiffness parameters was deemed necessary as these volumes only slightly exceeded the expected (segmented) values.

It should be noted that without interaction with the pericardium, the left ventricular volume showed similar characteristics as with interaction. Yet, the atrial and right ventricular volumes surpassed their segmented diastatic volume (exceeding a volume of 180 mL in case of the right ventricle) illustrating the importance of modeling the pericardial interaction. Additionally, it should be pointed out that the initialized end-diastatic volumes further changed during the course of multiple heart beats, converging toward a left ventricular end-diastatic volume of 136 and 168 mL after three heart beats.

Figure 5 shows the resulting deformation and Figure 6 the corresponding pressure and volume traces. Additionally, in more detail, a complete cardiac cycle can be seen in Supplementary Video 1. The simulated heart beat shows all major characteristics expected from a physiological circulatory point of view. First, during the atrial contraction phase (0–150 ms), the atria contract, leading to an increase in ventricular volume by ≈11/21% (left/right) (Figure 6). This is followed by the isovolumetric contraction phase of the ventricles, marking the beginning of systole (225–250 ms). Here, the pressure increases rapidly until it surpasses the “diastolic” aortic or pulmonary artery pressure (85/13 mmHg respectively). Subsequently, during the ejection phase (250–350 ms), the ventricles contract, exhibiting noticeable wall thickening, and eject 58/49% of their end-diastolic content. Once the individual ventricular pressure drops below the “systolic” arterial pressure (125/20 mmHg) of the respective outflow vessel, but is still above the atrial pressure, isovolumetric relaxation takes place (350–400 ms), indicating the end of the systole and the beginning of the diastole. After ventricular pressure drops below the atrial pressure, the ventricles relax back to their initial end-diastatic state (400–800 ms).

To assess the occurring displacement during the heart cycle, both the rotational movement and the atrioventricular plane displacement (AVPD) were further assessed. Rotation was measured at 5 points of the endocardial side of the left ventricular free wall (Figure 7). The overall twisting motion is hereby in concordance with published data from Sengupta et al. (2008) in terms of characteristic and amplitude during the different contractile phases. During atrial contraction, there is a small clockwise (viewing from base toward apex) rotation in the base and a minor counterclockwise rotation in the apex. During isovolumetric contraction, this rotation then starts to flip to a rapid counterclockwise rotation in the base and a clockwise rotation in the apex. The peak twisting angle of −5° in the base and 9° in the apex is then reached amid the ejection phase (−5° and 10° in the measured data). Subsequently, during relaxation, the ventricle rotates back in two phases (fast/slow), a characteristic which can also be observed in the data measured by Sengupta et al. (2008).

AVPD was extracted at three points: at the base of the left ventricular free wall, at the base of the right ventricular free wall, and at the base of the septal wall, see Figure 8. Positive displacement represents a movement of the valve toward the apex and negative toward the top of the atria. Due to the atrial contraction, the valve plane first moves upwards by –6/–7/–15 mm (LV/septum/RV) to then move downwards during ventricular ejection by 14/11/27 mm. The resulting traces are comparable to reported data from Maffessanti et al. (2013), Fritz et al. (2014), and Pfaller et al. (2019). Unfortunately, as the original Cine-MRI data-set from Keller et al. (2011) was recorded solely in short-axis view, we could not validate against the volunteer's AVPD.