### Description of computational model

The mathematical model used in this study consisted of four components: (1) an electrical conduction model of the ventricular tissue, (2) the Purkinje fiber that transmits electrical impulses to the ventricular tissue, (3) a mechanical contraction model of the ventricles, and (4) a lumped-parameter model of the cardiovascular system [28]. Figure 1a shows a schematic of the complete electromechanical model of the electrical and mechanical components and their coupling.

To construct a computational model of the cardiovascular system, we used a 3D finite-element model of a failing human ventricle that was developed in previous studies [27, 29,30,31] in combination with a lumped-parameter model of the circulatory system [28]. The model was constructed with a ventricular geometry based on MR and DTMR imaging of the heart [27]. According to the imaging data, the ventricles were segmented as explained in previous studies, and a finite-element mesh was generated using the segmented images [32, 33]. Fiber and sheet structural information was also assigned into the mesh. Subsequently, two-dimensional Purkinje fibers, as described by Berenfeld and Jalife [34], were mapped onto the 3D surface of the endocardium in the ventricles [35].

### Electrical model

For the electrophysiological simulation of the ventricles, we used the HyperMesh software to generate finite tetrahedral linear elements (214,319 nodes and 1,061,379 elements). The mesh represents a realistic heart that consists of the epicardium, mid-myocardium, endocardium, and Purkinje fibers. The electrical component of the model simulated the propagation of the AP in the ventricular tissue by solving an electrical conduction equation. This equation represents a continuum of the current flow through cardiomyocytes that are electrically connected. We used the membrane dynamic model proposed by Tusscher et al. [36], which had been validated. The following partial differential equations for reaction–diffusion, which were proposed by Vigmond et al. were used for the electrical conduction in the 3D ventricular tissue [18]:

$$\nabla \,.\,\widetilde{\sigma }\nabla V_{m} = \beta I_{m}$$

(1)

$$I_{m} = C_{m} \frac{{\partial V_{m} }}{\partial t} + I_{ion} (V_{m} ,\upsilon ) - I_{trans} ,$$

(2)

where \(\tilde{\sigma }\,\) represents the intracellular conductivity, *β* represents the surface-to-volume ratio of the cardiac cells, *C*_{m} represents the capacitance of the cells per unit of surface area, *V*_{m} represents the membrane potential, and *I*_{trans} represents the current density of the transmembrane stimulus. Additionally, *I*_{ion} represents the total transmembrane ionic current and is given as [36].

$$\begin{aligned} I_{ion} & = I_{Na} + I_{K1} + I_{to} + I_{Kr} + I_{Ks} + I_{CaL} \hfill \\ & \quad + I_{NaCa} + I_{NaK} + I_{pCa} + I_{pK} + I_{bCa} + I_{bNa} , \hfill \\ \end{aligned}$$

(3)

where *I*_{Na}, *I*_{K1}, *I*_{to}, *I*_{Kr}, *I*_{Ks}, *I*_{Ca,L}, *I*_{NaCa}, *I*_{NaK}, *I*_{pCa}, *I*_{pK}, *I*_{bCa}, and *I*_{bNa} represent the rapid inward Na^{+} current, inward rectifier K^{+} current, transient outward K^{+} current, rapid delayed rectifier K^{+} current, slow delayed rectifier K^{+} current, L-type Ca^{2+} current, Na^{+}/Ca^{2+} exchanger current, Na^{+}/K^{+} pump current, plateau Ca^{2+} current, plateau K^{+} current, background Ca^{2+} current, and background Na^{+} current, respectively.

To model the excitation–contraction (EC) coupling condition, the electrical component was coupled with a mechanical component, as shown in Fig. 1a. Cardiac EC coupling occurs during cellular depolarization in the electrical component, initiating the release of Ca^{2+} from the sarcoplasmic reticulum [37]. Subsequently, the cooperative bindings of Ca to troponin C and cross-bridge cycling occur. The cross-bridge cycling forms the basis for protein movement contractility and the development of active cellular tension, thereby resulting in the deformation of the ventricles. The Ca^{2+} transient response, obtained by the electrical component, serves as an input to the contractile myofilament dynamic model in the mechanical component (see Fig. 1a).

### Mechanical model

The mathematical model of the mechanical contraction of cardiac tissue was based on continuum mechanics [38, 39], with the assumption that the myocardium is hyperelastic and almost incompressible. Additionally, the passive mechanical properties are described as follows [40]:

$$W = \frac{C}{2}\left( {e^{Q} - 1} \right)$$

(4)

$$Q = b_{1} E_{ff}^{2} + b_{2} \left( {E_{rr}^{2} + E_{cc}^{2} + E_{rc}^{2} } \right) + 2b_{3} \left( {E_{fr}^{2} + E_{fc}^{2} } \right),$$

(5)

where *W* represents the strain energy function, and E_{ij}, the Green–Lagrange strain, represents the local fiber coordinate system. Information on the fiber orientation and laminar sheet was used to determine the orthotropic electrical conductivity and passive mechanical properties of the ventricular myocardium.

In this study, the differential equations of the cross-bridge model of muscle contraction were based on Rice et al. [41]. By assuming an isosarcometric simulation, *dSL/dt* was set as 0, and *SL* (the sarcomere length) was fixed at the initial value *SL*_{0}. When the sarcomere contracts, *SL* is computed as follows:

$$\frac{d}{dt}SL = \frac{{Integral_{Force} + (SL_{0} - SL) \times viscosity}}{mass},$$

(6)

where the viscosity and mass (especially in the mechanical component) are shown in Fig. 1a. *Integral*_{Force} represents the normalized forces and is calculated as

$$Integral_{Force} = \int\limits_{0}^{1} {\left( {F_{active} (x) + F_{passive} (x) - F_{preload} - F_{afterload} (x)} \right)dt} ,$$

(7)

where *F*_{active} (*x*) is an active force, *F*_{passive} (*x*) is a passive force, and *F*_{preload} is a constant force (an applied force that induces an initial *SL*). In the case of *F*_{afterload}, two conditions can be used: an isotonic contraction or a fixed-muscle length (isometric) contraction. For the isotonic contraction, the afterload force is fixed after the release. Conversely, for the fixed-muscle length contraction, the afterload is calculated as a series elastic element (as shown in the mechanical component in Fig. 1a), as follows:

$$F_{afterload} (x) = KSE \times (x - SL_{0} ),$$

(8)

where *x* represents the *SL*, and *KSE* represents the stiffness (*Force*/µm).

To simulate the hemodynamic responses, the electromechanical model was coupled with the lumped-parameter model of the circulatory system proposed by Kerckhoffs et al. [28]. This lumped-parameter model is shown in the mechanical component in Fig. 1a and is expressed as follows:

$$- R_{SA} \dot{Q}_{SA} + \frac{1}{{c_{SA} }}Q_{SA} = V_{SV}$$

(9)

$$- R_{SV} \dot{Q}_{SV} + \frac{1}{{c_{SV} }}Q_{SV} = V_{RA}$$

(10)

$$- R_{RA} \dot{Q}_{RA} + \frac{1}{{c_{RA} }}Q_{RA} = V_{RV}$$

(11)

$$- R_{RV} \dot{Q}_{RV} + \frac{1}{{c_{RV} }}Q_{RV} = V_{PA}$$

(12)

$$- R_{PA} \dot{Q}_{PA} + \frac{1}{{c_{PA} }}Q_{PA} = V_{PV}$$

(13)

$$- R_{PV} \dot{Q}_{PV} + \frac{1}{{c_{PV} }}Q_{PV} = V_{LA}$$

(14)

$$- R_{LA} \dot{Q}_{LA} + \frac{1}{{c_{LA} }}Q_{LA} = V_{LV}$$

(15)

$$- R_{LV} \dot{Q}_{LV} + \frac{1}{{c_{LV} }}Q_{LV} = V_{SA} ,$$

(16)

where *R* represents the resistance, *Q* represents the flux, *C* represents the compliance, *V* represents the volume, *SA* represents the systemic artery, *SV* represents the systemic vein, *RA* represents the right atrium, *RV* represents the right ventricle, *PA* represents the pulmonary artery, *PV* represents the pulmonary vein, *LA* represents the left atrium, and *LV* represents the left ventricle.

### Heterogeneous myocardium

The heart wall consists of the epicardium (outermost heart layer), mid-myocardium (cardiac muscle tissue in the middle layer of the ventricles), and endocardium (innermost ventricular surface) [2, 3]. In the present study, we used three different 3D ventricular models—called Models 1, 2, and 3—as shown in Fig. 1b. The models differ with regard to the thicknesses of each ventricular layer. Model 1 exhibits the thickest mid-myocardium and the thinnest epicardium and endocardium. Model 3 exhibits the thinnest M cell layer and the thickest epicardium and endocardium. Model 2 exhibits intermediate-thickness layers of the heart walls.

At a position 1 cm from the base (see Fig. 1b), the proportions of the epicardium and endocardium layers of Model 1 were 12% in the right ventricle (RV) wall and 9% in the left ventricle (LV) wall. The total percentage of the endocardium layer was 22% in the septum wall (11% in each side of the septum wall—right and left). In Model 2, the proportions of the epicardium and endocardium layers were 23.5% in the RV wall and 18% in the LV wall. The endocardium layer in Model 2 corresponded to 19% for each side of the septum wall. In Model 3, the proportions of the epicardium and endocardium layers were 41% in the RV wall, 22% in the LV wall, and 26% for each side of the septum wall.

### Simulation protocol

We investigated two scenarios: (1) a normal sinus rhythm, representing the control group, and (2) reentry, mimicking a pathological condition. We observed and compared the electromechanical responses of three different ventricular models (Models 1, 2, and 3) during the normal sinus rhythm and reentry, which was associated with a tachyarrhythmia condition. For each scenario, we performed two simulations: electrical and mechanical. For the normal sinus rhythm, a basic cycle length (BCL) of 600 ms was used for the electrical simulation. The sinus rhythm was performed for 3 s to obtain a steady-state condition. However, only the data for the last cycle (2.4 to 3 s) were used to simulate the ventricular electromechanics. The electrical impulse was initiated via Purkinje fibers with an electrical conduction velocity (CV) of 200 cm/s. Additionally, a myocardial CV (MCV) of 70 cm/s was used for the normal condition [6].

For the electrical simulation with the normal sinus rhythm, the electrical activation time (EAT) and electrical deactivation time (EDT) were calculated. The EAT and EDT are the times required for tissue depolarization and cellular repolarization, respectively. In the simulation, the EAT was the time when the membrane potential of the ventricles exceeded − 30 mV, and the EDT was the time when the membrane potential of the ventricles decreased below − 75 mV.

The reentry scenario was initiated using an S1–S2 protocol that was introduced by Tusscher et al. [36]. Three S1 stimuli were applied in the apex of the ventricles with a BCL of 600 ms. The stimuli generated waves that propagated in all directions. After the refractory tail of the wave passed through half of the length of the medium, an S2 stimulus was applied parallel to the S1 stimulus at three-quarters of the medium length. This produced a second wavefront with a curly tip and generated a reentrant spiral wave.

In the present study, the ventricular electromechanical simulations of the reentry scenario were conducted over 10 s, and two different initialization methods were used to examine the reentry arrhythmogenesis and the instability of the reentrant waves in the electrical simulation. In the first initialization method, the electrical simulation was initiated via the S1–S2 protocol as described previously, and the MCV was 43 cm/s. In the second method, the initial condition was from the outputs that were derived from the previous electrical simulation using the S1–S2 protocol that already generated reentry. The MCV of this method was 70 cm/s.

The results of the electrical simulation for both the sinus rhythm and the reentry conditions were coupled with the intracellular Ca^{2+} transient and were used as inputs for the mechanical simulation. Subsequently, the results of the mechanical simulation were used to analyze the ventricular mechanical responses, including the pressure, volume, total adenosine triphosphate (ATP) consumption rate, ejection fraction (EF), stroke volume (SV), stroke work (SW), and cardiac output (CO).

The SV represents the amount of blood pumped from the LV and is obtained by subtracting the end systolic volume (ESV) from the end diastolic volume (EDV). The EF represents the amount of blood ejected from the ventricles with each heartbeat and is obtained by dividing the SV by the EDV and multiplying the result by 100%. The SW corresponds to the work or pressure of blood processed in one cycle during ventricular contraction and is obtained by integrating the enclosed area of the pressure–volume (P–V) curve (Fig. 4b). The CO represents the amount of blood pushed out by the heart per minute.