The simulations are performed on a truncated ellipsoidal geometry with a rule-based myofiber field, which represents an idealized LV (Bovendeerd et al., 2009). The wall and cavity volumes are set to 136 and 44 mL. In the ellipsoid, two infarct geometries are considered, where the infarct core extends from 14% to 51% of the ventricle height (Figure 2A). A 100%-transmural core enclosed by a 9 mm thick BZ is implemented to ensure the induction of VTs, as it is known that larger BZs give higher risk for VT. Secondly, an endocardial 67%-transmural core with a 9 mm thick BZ is implemented to investigate whether the extension has the same effect for a non-transmural infarct. The total infarct volumes, including core and BZ, are 15 mL and 11 mL respectively.

The model of Janssens et al. is applied to simulate the mechanical behavior of a chronic infarct ventricle (Bovendeerd et al., 2009; Janssens et al., 2023). The material model and parameter values are copied from this study as well. A brief summary of their methods will be given in this section.

The myocardium is modeled as a non-linear, transversely isotropic, nearly-incompressible material. The Cauchy stress tensor is given by σ = f p a s   σ p a s + f a c t   σ a c t t , l   e ⇀ f e ⇀ f . (2.1)

This tensor is composed from two parts to account for both passive, σ p a s , and contractile active, σ a c t , myocardial stresses. σ a c t depends on the time t elapsed since cardiac activation, and on the length l of the contractile element. This active component is only present in parallel with the myofiber orientation e ⇀ f . To simulate a chronic infarct, the passive stress σ p a s is increased by changing the multiplication factor f p a s from 1 to 10 in the infarct core to represent stiff fibrosis. Furthermore, the active component is eliminated by setting the multiplication factor f a c t to 0 to represent the loss of contractile myocytes in this area. To simulate the structural BZ, which is a mixture of infarct and healthy tissue, both f p a s and f a c t values are linearly interpolated from the infarct core to the surrounding healthy tissue, where they are set to 1. The infarct size and geometry differ from those in Janssens et al., and are described in Sec. 2.1.

The ventricle is coupled to a 0D closed-loop lumped parameter model, which represents the circulatory system. Finally, the equation of conservation of momentum is solved: ∇ → ∙ σ = 0 → . (2.2)

Rigid body motion is suppressed by restricting the in-plane solution of basal nodes to its nullspace, which suppresses the average rotation about the long axis as well, and by constraining the base in the out-of-plane direction. The endocardial surface is subject to a uniform pressure. During isovolumetric contraction and relaxation, this pressure is determined such that mechanical equilibrium is obtained at a constant LV volume. The pressure during filling and ejection is computed from the interaction with the circulatory system. The simulated cardiac cycle is initiated by homogeneous activation of the myocardium. Strains are computed by using the LV geometry at t = 0 as reference, and fiber-strain amplitudes are defined as the difference between the local maximum and minimum logarithmic fiber-strain as experienced during the computed cardiac cycle.

Unstructured finite-element tetrahedral meshes are generated with an average resolution of 3 mm and 500 μm for the mechanics and electrophysiology model respectively (Figure 2F). Mechanics simulations are performed in FEniCS with timesteps of 2 ms (fenicsproject.org). The displacement field is represented by quadratic Lagrangian basis functions, corresponding to a system of 106,852 degrees of freedom. Electrophysiology simulations are performed in the CARP framework with timesteps of 20 μs (opencarp.org). The membrane voltage potential is approximated by linear Lagrangian basis functions, corresponding to 108,2784 degrees of freedom. To incorporate mechanical results into the electrophysiology simulation, fiber-strain amplitudes are linearly interpolated from the lower resolution mechanics mesh to the higher resolution electrophysiology mesh. From these interpolated strain-amplitudes, CV_T values are obtained through Eq. 2.4 and discretized in 9 levels from which the center values are imposed to the extended electrophysiology model using the methods of Costa et al. (Costa et al., 2013).

The effect of strain-induced CV_T remodeling on the electrophysiological simulation outcomes is explored for the two infarct cases as depicted in Figure 2A. For each case, one mechanics cardiac cycle was simulated. Fiber-strain amplitudes from these simulations were used in the extended model to modify the distribution of CV_T. Virtual stress test simulations from Sec. 2.2.3 are performed for both the reference and extended model to assess the effect of strain-dependent CV_T.

To investigate the sensitivity of the results to the parameter settings of the novel strain amplitude-CV_T relation, the values of α (steepness) and ε C V T , r e f (shift) were varied with respect to the baseline values as presented in Table 1 α was set to 50%, 70%, 100%, 150% and 250% compared to the default value. ε C V T , r e f was set to 122%, 111%, 100%, 89%, and 78%. This results in even distributions of strain amplitude-CV_T curves as shown in the top row of Figure 7, which are applied to the transmural infarct case. For datasets with overall higher measured strains (large ε C V T , r e f ), typically a lower sensitivity to variations in strain (small α ) is expected. Therefore, when both parameters are varied, only smaller values for α are considered for larger values in ε C V T , r e f .

The simulation results are analyzed in terms of the spatial repolarization time gradients (RTG) map. This map is created for the last S1 stimulus in the simulation, where repolarization time is defined as the moment at which the membrane voltage repolarized to −70 mV. From these RTG maps, RTG _vol and RTG _mean metrics are extracted where RTG_vol [ml] represents the tissue volume with gradients ≥10 ms/mm, which are critical for conduction block (Restivo et al., 1990). RTG_mean [ms/mm] indicates the mean RTG value within this volume with RTGs above 10 ms/mm. These metrics are computed for all performed experiments, such that the effects of mechanical feedback and parameter variability on RTG can be explored. The RTG metrics during baseline rhythm are an indication for VT risk (Cluitmans et al., 2021), and only require S1 stimuli. For the sensitivity test, these metrics are computed for one pacing site only (site A in Figure 2E) to save computational costs.

For the virtual stress test simulations, all six pacing sites are used and S2 stimuli are applied as well. The results are captured in vulnerable windows. These windows define for each pacing site and S1-S2 interval whether VT was induced or not. Here, VT is defined as electrical activity sustained for at least 800 ms after S2. In case of occurrence of VT, exit points are defined as the location with earliest second activation time, counted from the start of S2 sustained for at least 800 ms after S2. In case of occurrence of VT, exit points are defined as the location with earliest second activation time, counted from the start of S2.

To include the effect of mechanical tissue load, we introduced several modifications to our reference electrophysiology methods. First, in the reference model, CV_T is altered in the structural BZ based on measurements in canine infarcts (Yao et al., 2003). Supported by cultured cell experiments, we proposed a phenomenological model in which CV_T is altered by recurring strain amplitudes, on top of the effect of structural abnormalities. We found a strong positive correlation (r² = 0.952) between experimental CV_T (Poelzing et al., 2004) and simulated myofiber strain amplitude, as they were both lower at the epicardium and higher in the last 60% towards the endocardium. We fitted the model on this correlation, by setting the central CV_T to 0.3 m/s as is done conventionally, around which CV_T varies due to the effects of strain. This adaptation results in different CV_T values for infarct and healthy regions as strain modulates around the stiff infarct.

Second, the effect of fibrotic density on CV_T is modelled as a continuous linear relation, instead of the discrete 3-level relation assumed in conventional models. In this way, the highest fibrotic density naturally results in no electrical activity in the non-viable infarct core. Additionally, the spatial CV_T distribution is more heterogeneous as it results from the combined effects of fibrosis and fiber-strain amplitudes. These amplitudes are, for example, not constant over the wall depth, and are modulated nearby the stiff infarct. The CV_T distribution is therefore discretized in 10 levels instead of the aforementioned 3 (core, structural BZ, and remote) in the reference model.

Thirdly, in order to follow the Poelzing data used in the relationship between experimental CV_T and strain, we introduced in our model the transmural APD gradient present in those experimental canine data too. The presence of such transmural APD gradients in human data has been a topic of discussion (Boukens et al., 2015). For completeness, we have performed simulations without the transmural APD gradient too (Supplementary Appendix A). Here we implemented a homogeneous APD (287 ms when paced with 600 ms intervals) over the wall depth. For these simulations, RTGs were lower than for the simulations with transmural APD gradient, and no VTs could be induced. This might be caused by the smaller difference in APD between remote and border zone tissue in this model. Moreover, the occurrence of VT is a binary ‘all-or-nothing’ response, sensitive to various modeling choices. Introducing the effect of strain on CV_T might result in VTs in more critical cases. However, the impact of mechanics on electrophysiology substrate is similar under both conditions (with and without transmural APD gradient), enhancing RTGs on the border of the infarct area.

Fourth, we apply a weak coupling between the model of mechanics and electrophysiology. This is motivated from the long time-scale of conduction remodeling observed in cultured cell experiments. This remodeling took a few hours of applied pulsatile stretch, and persisted several hours after the stretch was removed. We assumed that recurring ‘pulsatile’ strains are mimicked during hours of sinus rhythm, which will induce remodeling of CV_T. These effects would persist long enough to remain in the short timeframe during which stress test pacing is applied. Additionally, the period of pacing is too short to induce further remodeling. It was therefore sufficient to only simulate mechanics during sinus rhythm, from which the results are used to define a static distribution of CV_T. The numerical uncoupling of mechanics and electrophysiology saves computational costs, and easily allows usage of different simulation environments and mesh compositions for both domains. The downside is that it overlooks instantaneous processes, like the (dynamic) response of strain-sensitive ion channels. However, this was not needed for the focus of our study on tissue remodeling (permanent) over time, in search for underlying mechanisms of VT evolution.

Our results show a significant effect of mechanical load on the simulation outcome of cardiac electrophysiology. For the healthy LV, transmural variety in fiber strains resulted in heterogeneity in CV_T over the ventricular wall depth (Figure 4). The largest strains and, therefore, the highest CV_T are found at the endocardium, which is in agreement with experimental results (Poelzing et al., 2004). For the infarcted LV, elevated or reduced fiber-strain amplitudes caused further CV_T heterogeneity at the infarct area. It should be noted that this resulted in CV_T slowing in remote regions adjacent to fibrotic areas, opposed to conventional models where conduction slowing is limited to structural regions only. Including this functional BZ within regions of impaired mechanics might be a first step in modeling continuous infarct remodeling and expansion.

A second effect of the extension is an increase in tissue volume with severe conduction slowing and in RTG_vol, which are both substrate metrics. The increase in RTG_vol was localized at the structural BZ boundary (Figure 5). We think the lowered CV_T adjacent to the structural BZ facilitates cell-to-cell uncoupling on the BZ boundary. This enhances the expression of differences in APD between the structural BZ (longer) and adjacent remote tissue (shorter), increasing RTG on the boundary in between. The elevated RTG caused the lines of conduction block to be longer, resulting in larger reentrant circuits. This, in combination with the lowered CV_T in the remote part of reentrant pathways, can prolong reentrant cycle time and more easily result in sustained VT. Those effects might be even more pronounced if cardiomyocyte and APD heterogeneity after infarction is caused by mechanical load as well, as suggested in recent work (Amoni et al., 2023).

LV ejection fraction reduced from 61% in the healthy case to 55% and 56% in the two infarct cases, and did not indicate the difference in VT inducibility. It may show that the LV ejection fraction as current clinical marker is a too global indicator for mechanical function and VT risk (Deng et al., 2016). Our hypothesized strain-induced remodeling might be a way to include local mechanical dysfunction in the assessment of VT risk.

The hypothesized strain-induced remodeling of CV_T might help detect critical regions in the modelled heart. Our extended model could replicate the VT morphologies and their exit sites as identified by the reference model. This is desirable as in previous studies the exit sites as predicted by conventional models and observed retro- and prospectively matched quite well for several swine and human cases (Prakosa et al., 2018). On top of that, novel induced VTs were observed which can lead to identification of additional ablation targets.

The methods presented in this study might be a first step in modeling the effects of pathology induced changes in mechanics on VT risk, which is of interest as post-infarct patients can develop heart failure or dilation. This can alter strain, and a correlation between those diseases, heterogeneous connexin distribution, dispersed impulse conduction, and VT risk is shown (Peters et al., 1997; Kitamura et al., 2002 van Rijen et al., 2004; Boulaksil et al., 2010). However, further experiments should be performed to validate our results and the implemented strain-induced remodeling of CV_T, as it remains a theoretical hypothesis unproved.

In this study, the phenomenological model proposed captures the effects of strain on CV_T, implicitly capturing the role of connexin. Connexin proteins are often assumed to be the link between strain and conductivity, because strain affects connexin expression (Shanker et al., 2005), which in turn affects myocardial conductance (van Rijen et al., 2004; Boulaksil et al., 2010). In future work, explicitly including connexin in a more detailed model might help to obtain additional insight into the process of CV_T remodeling over time.

Both our mechanics and electrophysiology models are partly based on animal data. For example, the reference electrophysiology model from literature already included conduction slowing and a prolonged APD in the border zone, which they based on canine data. In our model extension, we had to determine the CV_T-strain relationship from canine CV_T data too. Although we believe our current approach is sufficiently relevant to human conditions to study the potential modulating effects of mechanical strain on conduction velocity, this species mismatch should preferably be addressed in future studies.

Strain data from endo-to epicardium in literature is divers (Ashikaga et al., 2007; Bogaert and Rademakers, 2001; Costa et al., 2013; Takayama et al., 2002), which results in uncertainty in the accuracy of the relationship between strain and CV_T proposed. Therefore, we investigated the effect of parameter variation in the extended model. In particular, we focused on parameters that alter the steepness and shift of the relation proposed (Figure 7). We found the definition of the relation for low strain amplitudes (˂0.12) to be critical. Relational changes for this range affected RTG results the most, because those low amplitudes are found in and around the BZ, where they determine the amount of conduction slowing. This finding stresses that further research is needed to investigate whether the proposed relation correctly describes strain-induced remodeling of CV_T, especially for lower strain amplitudes. Further improvements may be obtained with personalized models, which allows for validation by linking patient-specific strain-imaging data.

In the model of mechanics, chronic infarct areas were implemented by increasing the stiffness as a result of fibrotic tissue, and by reducing contractility due to a loss of myocytes. However, contraction was initiated simultaneously throughout the LV, neglecting asynchronous electrical activation. As a next step, this can be included which would allow simulation of a feed forward cycle to describe post-infarct tissue evolution. Altered strains around the infarct may affect the electrical conductivity, which will affect mechanical activation in turn and therefore further modulate the strain. This cycle can then be repeated. Such a feed forward cycle is not yet included.

Finally, in this work, only one mechanics cardiac cycle was simulated which is not enough to reach hemodynamic steady state. To investigate the consequences, we performed a convergence study in which we ran the model for transmural infarct for 5 subsequent cycles. On top of that, both the mechanics and electrophysiology mesh was refined, and the number of CV_T discretization levels was doubled. Strain amplitudes remained similar in remote areas while low strain amplitudes in the BZ region (±0.15) further reduced by up to 0.02. The spatial CV_T distribution remained similar and CV_T in the BZ reduced maximally with 0.06 m/s. The lowered amplitudes and CV_T further enhance the effect of the mechanics extension, and the effect on the critical RTG substrate remained the same (Supplementary Appendix A).