The role of mechano-electric feedbacks and hemodynamic coupling in scar-related ventricular tachycardia

by   Matteo Salvador, et al.

Mechano-electric feedbacks (MEFs), which model how mechanical stimuli are transduced into electrical signals, have received sparse investigation by considering electromechanical simulations in simplified scenarios. In this paper, we study the effects of different MEFs modeling choices for myocardial deformation and nonselective stretch-activated channels (SACs) in the monodomain equation. We perform numerical simulations during ventricular tachycardia (VT) by employing a biophysically detailed and anatomically accurate 3D electromechanical model for the left ventricle (LV) coupled with a 0D closed-loop model of the cardiocirculatory system. We model the electromechanical substrate responsible for scar-related VT with a distribution of infarct and peri-infarct zones. Our mathematical framework takes into account the hemodynamic effects of VT due to myocardial impairment and allows for the classification of their hemodynamic nature, which can be either stable or unstable. By combining electrophysiological, mechanical and hemodynamic models, we observe that all MEFs may alter the propagation of the action potential and the morphology of the VT. In particular, we notice that the presence of myocardial deformation in the monodomain equation may change the VT basis cycle length and the conduction velocity but do not affect the hemodynamic nature of the VT. Finally, nonselective SACs may affect wavefront stability, by possibly turning a hemodynamically stable VT into a hemodynamically unstable one and vice versa.



There are no comments yet.


page 7

page 9

page 12

page 15

page 18

page 20


A convergent numerical scheme for a model of liquid crystal dynamics subjected to an electric field

We present a convergent and constraint-preserving numerical discretizati...

Modeling flexoelectricity in soft dielectrics at finite deformation

This paper develops the equilibrium equations describing the flexoelectr...

Myocardial ischemic effects on cardiac electro-mechanical activity

In this work, we investigated the effect of varying strength of Hyperkal...

Ballooning in Spiders using Multiple Silk Threads

In this paper, three-dimensional numerical simulations of ballooning in ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Cardiac arrhythmias result from an irregular electrical activity of the human heart. Among them, ventricular tachycardia (VT), which manifests with an accelerated heart rate, is one of the most life-threatening rhythm disorders. VT may be classified as either hemodynamically stable or unstable, depending on the capability of the heart to effectively pump blood in the circulatory system. In the former case antiarrhythmic drugs are generally employed, while in the latter case cardioversion is needed

[Epstein2008]. According to the specific pathogenesis, the stability of the VT remains the same or changes over time. Moreover, it may also degenerate towards ventricular fibrillation (VF), a life-threatening condition in which the ventricular activity is fully disorganized and chaotic, leading to heart failure [Samie2001].

In the clinical framework, these pathological scenarios can be hardly ever fully investigated and predicted for all patients. For this reason, biophysically detailed computational heart models could be used to provide a deeper understanding of the hemodynamic response to VT and to characterize the electromechanical substrate leading to dangerous arrhythmias.

While electrophysiological simulations are well-established for scar-related VT identification and treatment on human ventricles [Arevalo2016, Cedilnik2018, Deng2019, Prakosa2018], patient-specific electromechanical models coupled with closed-loop cardiovascular circulation have been just recently used for the first time to enhance our knowledge on VT [Salvador2021]. Indeed, on one side, the physiological processes by which the mechanical behavior alters the electrical activity of the human heart, known as mechano-electric feedbacks (MEFs), are relevant and not fully elucidated [Kohl2004, Kohl2013, ColliFranzone2017, Timmermann2017AnIA, Taggart1999, Panfilov2010, Varela2021]. On the other hand, the identification of the hemodynamic nature of the VT has significant clinical implications [Salvador2021].

The goal of this study is threefold: (1) to combine electrophysiology, activation, mechanics and hemodynamics in several numerical simulations of scar-related VT, with the aim of uncovering the roles of myocardial deformation and the recruitment of SACs on VT stability; (2) to show that our computational model effectively reproduces both hemodynamically stable and hemodynamically unstable VT; (3) to capture relevant microscopic mechanisms during VT, such as the incomplete relaxation of sarcomeres, given the biophysical detail of our electromechanical model.

2 Mathematical models

We provide an overview of the 3D cardiac electromechanical model coupled with a 0D closed-loop model of the cardiovascular system (see [regazzoni2020model, Salvador2021]).

We depict in Figure 1 the computational domain , which represents a human left ventricle (LV) taken from the Zygote Solid 3D heart model [Zygote]. Its boundary is partitioned into epicardium , endocardium and base . We recall that .

Figure 1: Computational domain , given by a human LV.

2.1 3D-0D closed-loop electromechanical model

We consider a multiphysics and multiscale 3D-0D mathematical framework comprised of four different core models, namely cardiac electrophysiology [franzone2014mathematical, colli2018numerical, nobile2012active, levrero2020sensitivity], mechanical activation [niederer2011length, rossi2014thermodynamically, ruiz2014mathematical, land2017model, regazzoni2020biophysically], passive mechanics [guccione1991passive, guccione1991finite, ogden1997non, holzapfel2009constitutive] and blood circulation [quarteroni2016geometric, blanco20103d, kerckhoffs2007coupling, hirschvogel2017monolithic, augustin2021, arts2005adaptation, regazzoni2020model]. The volume conservation constraint  defines the coupling condition between the LV cardiac electromechanics and the remaining part of cardiovascular system. [regazzoni2020model]. All these models represent several physiological processes, ranging from the cellular level to the organ scale.

Given the computational domain and the time interval , the mathematical model reads:


in which the model unknowns are:


where is the transmembrane potential,

is the vector containing both gating and ionic variables,

are the states variables of the active force generation model, represents the tissue mechanical displacement, defines the state vector of the circulation model (including pressures, volumes and fluxes of the different compartments of the cardiocirculatory system) and and are the left and right ventricular pressure, respectively. Finally, all variables are endowed with suitable initial conditions in :


2.1.1 Electrophysiology

We model the electrical activity of the myocardium by means of Eq. (1a), that is the monodomain equation coupled with a suitable ionic model for the human ventricular action potential  [BuenoOrovio2008, franzone2014mathematical, LuoRudy1994, ten2006alternans].

In the electrophysiological model , is the surface area-to-volume ratio of cardiomyocytes, represents the transmembrane capacitance per unit area. The applied current mimics the effect of the Purkinje network [vergara2014patient, costabal2016generating, landajuela2018numerical] by triggering the action potential at specific locations of the myocardium. The reaction terms and (specified by the ionic model at hand) couple the action potential propagation to the cellular dynamics. Specifically, we use the ten Tusscher-Panfilov ionic model (TTP06), which is able to accurately describe ions dynamics across the cell membrane [ten2006alternans]. Furthermore, Eq. (1a) is equipped with homogeneous Neumann boundary conditions for at the boundary , which defines the condition of electrically isolated domain for .

represents the conductivity tensor, that reads:


where are the longitudinal, transversal and normal conductivities, respectively. The parameter takes into account the effect of scars, grey zones and non-remodeled regions. This parameter is also incorporated inside the formulation of the TTP06 model, as described in [Salvador2021].

The deformation tensor and its determinant are needed to perform the mechano-electric coupling [Trayanova2011]. Indeed, we model the so called mechano-electric feedbacks (MEFs) [Costabal2017, Hazim2021]. The geometry-mediated MEFs incorporate the effects of displacement on the cardiac tissue, while other physiological processes act at the level of single cardiomyocytes [Kohl2004, levrero2020sensitivity]. Among them, some examples are selective (e.g. -permeable) or nonselective SACs and intracellular calcium binding to sarcolemmal buffers, the latter requiring more sophisticated ventricular ionic models [Bartolucci2020, Tomek2019]. In this paper we model nonselective SACs by means of the following formulation [levrero2020sensitivity]:


where and represent the conductance of the channels and the reversal potential, respectively. SACs may alter the shape of the action potential (AP) by lengthening or shortening its duration (APD) and by generating higher resting potentials. This may induce early or delayed afterdepolarizations (EADs or DADs) and premature excitation.

Both geometry-mediated MEFs and nonselective SACs are generally known to be pro-arrhythmic in pathological scenarios, because they increase the likelihood of having extra stimuli during cardiac contraction and relaxation [Kohl2004].

In this work we consider Eq. (1a) with several degrees of complexity to assess similarities and differences in the outcomes of the electromechanical simulations during VT. We report in Tab. 1 the mathematical models and the parametrizations that we consider. In particular, the choice of three different mathematical models for the geometry-mediated MEFs, namely , and , is motivated by the numerous formulations of this type of feedback that can be found in the literature [Chapelle2009, Collin2019, levrero2020sensitivity, quarteroni2017integrated]. Indeed, we range from minimal to complete inclusion of geometry-mediated MEFs in the monodomain equation.

Model name Equation
Table 1: Modeling choices for the monodomain equation that have been used in this paper. We consider different parametrizations for in terms of and . indicates the conductivity tensor in Eq. (3) with .

2.1.2 Activation

To model how the calcium wave following the AP triggers a series of chemo-mechanical reactions within sarcomeres, resulting in the generation of an active force in the muscle, we use the mean-field version of the model proposed in [regazzoni2020biophysically], henceforth denoted by RDQ20-MF. This mathematical model is based on a biophysically accurate description of regulatory and contractile proteins and their dynamics. Thanks to suitable dimensionality reduction techniques, the dynamics of stochastic processes underlying both chemical and mechanical microscale transitions are described in only 20 ODEs (that is ). The computational cost of this model is thus comparable to that of phenomenological models [niederer_hunter_smith_2006, rice2008approximate, land2017model], while providing a biophysically detailed description that is consistent with the level of detail and mechanistic understanding of the ionic model needed for the present study.

The inputs of the RDQ20-MF model are the intracellular calcium concentration coming from , the sarcomere length (denoted by ) and its time derivative (the latter allows to account for the so-called force-velocity relationship [keener2009mathematical]). The variable is obtained as , where is the sarcomere slack length and the fourth invariant measures the tissue stretch in the fiber direction. Finally, the output of the RDQ20-MF model is the active tension generated at the microscale, that can be obtained as .

Differently from [Salvador2021], here we consider a biophysically detailed active stress model. For this reason, we do not need to put the parameter inside its parametrization. Indeed, the active force generation mechanisms are properly handled by the active stress model, which receives different intracellular calcium waves from according to the specific area of the myocardium (scar, grey zone or healthy) and provides physiological values of active tension in all cases.

2.1.3 Mechanics

We employ the momentum conservation equation reported in Eq. (0c) to model the dynamics of the displacement of the myocardium. We also consider the hyperelasticity assumption, so that the strain energy function can be differentiated with respect to the deformation tensor to obtain [ogden1997non, guccione1991finite].

In the mechanical model , represents the density of the myocardium. The Piola-Kirchhoff stress tensor is additively decomposed according to:


where the first term stands as the passive part of the tensor , while the latter as the active one; is the strain energy density function, is the active tension, provided by the activation model .

To model the passive behaviour of cardiac tissue, we employ the orthotropic Guccione constitutive law [guccione1991passive], whose strain energy function is defined as:


where the first term is the volumetric energy with the bulk modulus , which penalizes large variation of volume to enforce a weakly incompressible behavior [peng1997compressible, doll2000development], and the latter is the deviatoric energy where is the stiffness scaling parameter and the exponent reads:


in which are the entries of , i.e the Green-Lagrange strain energy tensor, being the right Cauchy-Green deformation tensor.

To model the mechanical effects of the pericardial sac [Gerbi2018monolithic, pfaller2019importance, strocchi2020simulating], we impose at the epicardial boundary a generalized Robin boundary condition by defining the tensors and where , , , are the stiffness and viscosity parameters of the epicardial tissue in the normal and tangential directions, respectively. Normal stress boundary conditions were imposed at the endocardium of the LV , in which is the pressure exerted by the blood in the LV, modeled by means of the 0D closed-loop circulation model [regazzoni2020model]. To take into account the effect of the neglected part, over the basal plane, on the ventricular domain, we set on the energy consistent boundary condition [regazzoni2020machine], where:


2.1.4 Blood circulation and coupling conditions

We model the blood circulation through the entire cardiovascular system by means of a closed-loop model recently proposed in [regazzoni2020model]. In the 0D closed-loop model, systemic and pulmonary circulations are modeled with RLC circuits, heart chambers are described by time-varying elastance elements and non-ideal diodes stand for the heart valves [regazzoni2020model].

We represent the circulation core model with a system of ODEs, which is expressed by Eq. (0d), where is a proper function (defined in [regazzoni2020model]) and includes pressures, volumes and fluxes of the different compartments composing the vascular network:

Here , , and refer to the volumes of left atrium, right atrium, left ventricle and right ventricle, respectively; , , , , , , and express pressures and flow rates of the systemic and pulmonary circulation (arterial and venous). For the complete mathematical description of the 0D circulation lumped model we refer to [regazzoni2020model]. To couple the 0D circulation model with the 3D LV model, given by , we follow the same strategy proposed in [regazzoni2020model]. In particular, we replace the time-varying elastance elements representing the LV in with its corresponding 3D electromechanical description and we introduce the coupling condition where:

in which is a vector orthogonal to the LV centerline (i.e. lying on the ventricular base) and lays inside the LV [regazzoni2020model].

Due to , in the 3D-0D coupled model (2.1), is not determined by the 0D circulation model (0d), but rather acts as Lagrange multipliers and enforces the constraint .

2.2 Reference configuration

Cardiac geometries are acquired from in vivo medical images through imaging techniques. These geometries are in principle not stress free, mainly because there is always a pressure acting on the endocardium. Therefore, one needs to estimate the unloaded (i.e. stress-free) configuration (also named reference configuration) to which the 3D-0D cardiac electromechanical model (

2.1) refers. To recover the reference configuration , starting from a geometry acquired from medical images , we apply the same procedure proposed for the left ventricle in [regazzoni2020model].

Specifically, we solve an inverse problem: find the solution of the following differential problem


obtained for and , such that we get the domain .

Finally, to properly initialize the numerical simulation, we inflate the ventricular reference configuration by solving problem (9), where we set the pressures . The value is chosen to bring the left ventricle to a defined volume , as explained in [Piersanti2021EM]. Then, the solution of the problem (9) is set as initial condition for (with ) in .

3 Numerical methods

Figure 2: Zygote LV with the idealized distribution of scars (black), grey zones (grey) and non-remodeled regions (white) over the myocardium. Volumetric view (left) and cut view (right). The first type of grey zone corresponds to , while on the second one is prescribed.

We depict the geometric model of the Zygote LV in Fig. 2. We consider an idealized distribution of ischemic regions, which is made by two scars and two different types of grey zones.

To generate the fibers distribution (field ) for our geometry, we use the Bayer-Blake-Plank-Trayanova algorithm [bayer2012novel, piersanti2021modeling] with , , , and .

We employ the strategy proposed in [Regazzoni2021] to choose proper initial conditions for the 3D electromechanical simulations, by relying on a 0D emulator. In particular, to reach a steady state for the electrophysiological variables, we trigger periodic stimuli with a period equal to , for a total duration of . On the other hand, for the 3D numerical simulations, we apply an S1-S2 stimulation protocol consisting of two gaussian stimuli, the first one applied at time and the second one applied at time . In this way the 3D electromechanical simulations are a natural continuation of the 0D initialization, and the VT can be directly induced by just using a single S1 before the S2.

We use a segregated-intergrid-staggered scheme to numerically discretize the electromechanical model [Piersanti2021EM]. Indeed, we solve , , and sequentially, by employing different time and space resolutions according to the specific core model.

We employ the Finite Element Method (FEM) for the space discretization of electrophysiological, activation and mechanical models [quarteroni2009numerical, quarteroni2019cardiovascular]. We use hexahedral meshes and finite element spaces for the space discretization of all the core models. This multiphysics problem presents different space resolutions according to the specific model at hand [regazzoni2020model, regazzoni2020numerical]. In the framework of cardiac electrophysiology, we consider a fine geometrical description (1’987’285 DOFs, 1’926’912 cells, ) to accurately capture the electric propagation due to fine-scale phenomena arising from the continuum modeling of the cellular level, especially with the aim of reproducing and properly address arrhythmias [salvador2020intergrid]. On the other hand, cardiac mechanics allows a lower space resolution (35’725 DOFs, 30’108 cells, ). This eases its numerical solution, which is computationally demanding, especially for the assembling phase [salvador2020intergrid]. Differently from [Salvador2021], in this paper we consider a significant space scales separation between electrophysiology and mechanics because we are dealing with a simple distribution of ischemic regions.

For time discretization, we use backward differentiation formula (BDF) schemes [quarteroni2009numerical]. In particular, we employ a third order BDF scheme for electrophysiology and a first order BDF scheme for activation, mechanics and circulation. These choices allow to accurately capture the fast time dynamics of the model variables without unbearable restrictions on the time step, while not introducing numerical instabilities. We treat the nonlinear terms coming from both electrophysiological and activation models in a semi-implicit fashion. Cardiac mechanics is numerically advanced in time with a fully implicit scheme. The cardiocirculatory model is solved with an explicit method [Piersanti2021EM]. To stabilize the numerical oscillations arising from the feedback of the tissue mechanics on the activation model, without having to resort to a monolithic strategy, we rely on the scheme proposed in [Regazzoni2021oscillation]. We use a smaller time step for electrophysiology () than for activation, mechanics and circulation (). We set the final time for all the numerical simulations.

The mathematical models of Sec. 2 and the numerical methods presented in this section have been implemented in lifex (, a high-performance C++ library developed within the iHEART project and based on the deal.II ( Finite Element core [dealII92].

The numerical simulations were performed on a HPC facility available at MOX for the iHEART project. The entire cluster is endowed with 8 Intel Xeon Platinum 8160 processors, for a total of 192 computational cores and a total amount of 1.5TB of available RAM.

4 Numerical results

We present electromechanical simulations to evaluate the effects of MEFs. We consider different modeling choices for the monodomain equation, as reported in Tab. 1. We start from a baseline simulation with model , in which we obtain a stable VT. Then, we compare the effects of different geometry-mediated MEFs, i.e. models , , and . We also study the impact of different parametrizations for SACs, i.e. , with respect to . Finally, we evaluate the combined effects of geometry-mediated MEFs and nonselective SACs. We report in Appendix A the values of the parameters that we use to get the numerical results that will be discussed in this paper.

4.1 Baseline simulation

Figure 11: Time evolution of the transmembrane potential (left) and displacement magnitude (right) for the Zygote LV with an idealized distribution of ischemic regions. Each picture on the right side is warped by the displacement vector . MEFs are neglected, i.e. we use model .
Figure 12: Comparison between a reference healthy PV loop in sinus rhythm (red, Appendix A, heartbeat period equal to ) and the one obtained in the baseline simulation under VT for s (light blue). We underline that, differently from [Salvador2021], here we induce a hemodynamically tolerated VT.

We depict in Fig. 11 the evolution of the transmembrane potential and the displacement magnitude for the baseline electromechanical simulation, where all MEFs are fully neglected (model in Tab. 1). We induce a sustained VT with a figure-of-eight pattern around the isthmus, which is laterally bordered by scars, that act as conduction blocks.

In Fig. 12, we compare the Pressure-Volume (PV) loop over different heartbeats for the baseline simulation under VT with a reference healthy PV loop in sinus rhythm obtained by removing scars and grey zones. We observe that the contractility increases while the stroke volume (SV) decreases. The ejection fraction (EF) remains approximately the same and we approach a steady state in which the electromechanical function is not impaired. For these reasons, we conclude that the VT is stable.

VT associated with ischemia are known to disturb the normal isovolumetric processes and to influence the end systolic/diastolic pressure volume relationship (ESPVR/EDPVR). Simultaneously, a phenomenon called incomplete relaxation may occur, especially when the VT does not leave enough time for the uncoupling of all the actin-myosin bonds between two consecutive contraction phases [Bastos2019]. The occurrence of this phenomenon is illustrated in Fig. 13, where we depict the time evolution of the minimum, maximum and average active stress in the computational domain for a reference healthy case in sinus rhythm and the baseline simulation under VT. Specifically, in the healthy case there is always a time interval between two consecutive heartbeats (precisely, during ventricular diastole) in which the active stress is virtually zero in the LV. In the VT case, instead, the cardiac muscle is never fully relaxed, thus not allowing the LV to complete its emptying. All these details are properly captured by our electromechanical model thanks to its biophysical accuracy.

Figure 13: Minimum, average and maximum active tension over time for a reference healthy case in sinus rhythm (left, Appendix A, heartbeat period equal to ) and the baseline simulation under VT (right). We see that incomplete relaxation occurs during VT.

4.2 Effects of geometry-mediated MEFs

Figure 19: Comparison among different models for geometry-mediated MEFs in terms of transmembrane potential .
Figure 20: Pointwise values of transmembrane potential , intracellular calcium concentration , sarcomere length , active tension , pressure and volume over time for , , and .
Table 2: BCL for different modeling choices in geometry-mediated MEFs. Model significantly changes BCL with respect to , , .

We consider four different modeling choices for the geometry-mediated MEFs, namely , , and in Tab. 1

, while for the moment we completely neglect the impact of SACs. Then, we perform four different electromechanical simulations by employing these four different formulations for the monodomain equation.

We illustrate in Fig. 19 the development of transmembrane potential over time. We observe minor differences in action potential propagation among , and . These differences, as we can see for , are mainly focused on the depolarization wave and occur during VT. Moreover, by looking at Tab. 2, we notice that the VT BCL is very similar among these three models. Indeed, the BCL is approximately equal to , which is long if compared to more dangerous VT and justifies a stable ventricular excitation.

On the other hand, entails major changes in VT BCL, which increases from (for model ) to , and conduction velocity, that significantly decreases.

These observations are in agreement with Fig. 20, where the electrophysiological, mechanical and hemodynamic variables retrieved in a random point of the computational domain are shifted forward for , while , and show a very similar pattern. This is also motivated by the change in the VT exit site, as we can see from Fig. 19 for . This phenomenon is particularly evident from the plot of sarcomere length over time (Fig. 20), which also presents different peak values for and . Finally, we see that wave stability is not affected by geometry-mediated MEFs. Indeed, the VT always remains hemodynamically stable in the four different cases.

4.3 Effects of SACs

Model VT type
Stable ()
, , Stable ()
, , Stable ()
, , Unstable ()
, , Stable ()
Table 3: VT classification for different SACs parametrizations. The unstable VT has a BCL that ranges from to .

In this section, we fully neglect the effects of geometry-mediated MEFs and we focus on different parametrizations for SACs in terms of and . We also compare the outcomes of with the ones of to outline similarities and differences.

We notice from Fig. 26 that both APD and wave stability are affected by SACs parametrizations. Indeed, by combining the 3D information with the pointwise evaluations of Fig. 27, we discover that there is one choice of the parameters, namely and , that converts the VT from stable to unstable. The instability derives from an extra stimulus that is completely driven by contraction, which occurs in the superior-right part of the ventricle. This extra stimulus changes the VT morphology, along with its BCL, which is not the same over time.

By considering data of Tab. 3, we observe that the VT BCL remains the same when there is no stability transition. Moreover, from the hemodynamic perspective, the onset of different types of arrhythmias is driven by the combined effects of and . Indeed, in Fig. 27 we see that given , different may change wave stability. On the other hand, in Fig. 28 we show that different may affect wave stability, with fixed a priori.

Figure 26: Comparison between model and model , for different values of ().
Figure 27: Pointwise values of transmembrane potential , intracellular calcium concentration , sarcomere length , active tension , pressure and volume over time for and with different choices of ().
Figure 28: Pointwise values of pressure and volume over time for and with different choices of ()

4.4 Combined effects of geometry-mediated MEFs and SACs

Figure 34: Comparison among models , and .
Figure 35: Pointwise values of pressure and volume over time for , (, ) and (, ).
Figure 36: Coupled effects of electrophysiology, mechanics and hemodynamics for the numerical simulation with model . The extra stimuli in the upper right part of the LV, which is driven by SACs, activate the LV electrophysiologically and mechanically. This has a direct impact on both pressure and volume transients, which in turn have an effect on the electromechanical behavior of the LV.
Model VT type
Stable ()
, , Unstable ()
, , Unstable ()
Table 4: VT classification for combinations of geometry-mediated MEFs and nonselective SACs. The unstable VT related to has a BCL that ranges from to . The unstable VT related to has a BCL that ranges from to .

We briefly evaluate the combined effects of geometry-mediated MEFs and nonselective SACs. Once SACs parametrization is fixed, we notice that switching between no formulation (i.e. ) to full formulation (i.e. ) of geometry-mediated MEFs entails significant differences. In particular, from Fig. 34 we see that triggers the extra stimuli faster than . As we can notice from Fig. 35, the VT remains unstable but both pressure and volume traces over time are very different from each other for and . From Tab. 4, we infer that the VT BCL for is lower than the one of . This potentially defines a more dangerous VT. Finally, in Fig. 36 we highlight the joint contributions of electrophysiology, mechanics and hemodynamics in the coupled model.

5 Discussion

We investigate how several types of geometry-mediated MEFs and different parametrizations for nonselective SACs affect the electric and hemodynamic stability of VT. In particular, we focus on the sustainment and the morphology of VT macro-reentrant circuits and blood supply, which is analyzed by means of PV loops.

Differently from previous studies [ColliFranzone2017, Timmermann2017AnIA, Panfilov2010], we keep into account tissue heterogeneity by introducing an idealized distribution of scars, grey zones and non-remodeled regions over the myocardium. We also consider a more sophisticated coupled mathematical model that embraces electrophysiology, activation, mechanics and cardiovascular fluid dynamics to analyze these mechano-electric couplings. To the best of our knowledge, this framework enables to study for the first time the hemodynamic effects of both geometry-mediated MEFs and nonselective SACs on VT by means of electromechanical simulations, showing strengths and weaknesses of some simplifying modeling choices that are commonly made when simulating this type of phenomenon [Chapelle2009, Collin2019, levrero2020sensitivity, quarteroni2017integrated]. The coupling between the 3D electromechanical model and the 0D circulation model permits to identify the hemodynamic nature of the VT. With our approach, we discriminate between stable and unstable VT, which might result in hemodynamically tolerated or not tolerated VT [Salvador2021].

We see that if a VT is triggered by a certain stimulation protocol and by neglecting all MEFs, the very same pacing protocol induces a VT for all possible combinations of MEFs [Salvador2021].

In our numerical simulations, we do not observe a significant impact of geometry-mediated MEFs on the induction and sustainment of the VT. Most of the modeling choices for geometry-mediated MEFs, namely , and , present very similar VT BCLs and conduction velocities, while showing a few differences in the depolarization wave [ColliFranzone2017]. On the other hand, manifests major differences in the VT BCL, which increases, and its exit site with respect to , and , as shown in [Salvador2021] for a patient-specific unstable VT. This can be justified by the simplifications introduced in , while moving towards we almost totally recover the behavior observed in . Therefore, the minimal MEFs modeling choice might lead to biased results which are in favor of less severe VT.

We observe that nonselective SACs may affect the hemodynamic nature of the VT, as they may induce EADs or DADs, which lead to ectopic foci that reactivate the LV [Hazim2021]. These extra stimuli are generally located in the regions of the myocardium in which there is a transition between scar and border zone or between border zone and non-remodeled areas, where high stretches are likely to be present [Jie2010]. According to the specific combination of and , nonselective SACs affects both APD and AP resting values [Kohl2004]. We remark that such spontaneous arrhythmias triggered by myocardial stretches cannot be assessed in electrophysiological simulations, where the mechanical behavior is neglected.

Finally, in our study we also investigate the combined effects of geometry-mediated MEFs and nonselective SACs. Significant differences between and are observed when geometry-mediated MEFs are combined with a parametrization of SACs that entails extra stimuli. Specifically, the VT BCL with is lower than the one of because the extra stimuli driven by SACs is triggered more often in the former case. This completely changes the pressure-volume dynamics of the VT, whose stability is however still not affected by the geometry-mediated MEFs.

6 Limitations

There are other types of MEFs that could be investigated in future works. Among them, an important role is certainly played by the mechanical effects mediated by fibroblasts in the extracellular matrix, buffers handling, alterations in transmembrane capacitance due to local stretch and ions selective SACs [Kohl2013]. Indeed, modeling different cellular processes in cardiac electrophysiology might be of interest to further shed light on the underlying mechanisms of arrhythmias. Nevertheless, we have provided a broad study of geometry-mediated MEFs and nonselective SACs by combining electrophysiological, mechanical and hemodynamic observations.

Furthermore, MEFs should be properly quantified in patient-specific cases. Currently, our model for nonselective SACs permits to put one single value of both and for the whole geometry, which is unlikely to happen in realistic scenarios. A precise space assessment of SACs combined with our mathematical framework could have significant clinical implications.

7 Conclusions

We studied the effects of geometry-mediated MEFs and nonselective SACs on a realistic LV geometry endowed with an idealized distribution of infarct and peri-infarct zones. We performed numerical simulations of cardiac electromechanics coupled with closed-loop cardiovascular circulation under VT.

Our electromechanical framework allows for the hemodynamic classification of the VT, which can be either stable or unstable, and permits to capture mechanically relevant indications under VT, such as the incomplete relaxation of sarcomeres. Furthermore, by combining electrophysiology, activation, mechanics and hemodynamics, we observed several differences on the morphology of the VT with respect to electrophysiological simulations. In particular, geometry-mediated MEFs do not affect VT stability but may alter the VT BCL, along with its exit site [Panfilov2010]. On the other hand, the recruitment of SACs may generate extra stimuli, which may change VT stability. These extra stimuli are driven by myocardial contraction and are induced by changes in the APD or in the resting value of the transmembrane potential [Jie2010]. We conclude that both geometry-mediated MEFs and nonselective SACs define important contributions in electromechanical models with hemodynamic coupling, especially when numerical simulations under arrhythmia are carried out.


MS, FR, SP, LD and AQ acknowledge the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 740132, iHEART - An Integrated Heart Model for the simulation of the cardiac function, P.I. Prof. A. Quarteroni). NT acknowledges the National Institutes of Health (grants R01HL142496 and R01HL142893) and a Leducq Foundation grant.


Appendix A Model parameters

We provide the full list of parameters adopted for the numerical simulations of this paper. Specifically, Tab. 5 contains the parameters related to the electrophysiological model, Tab. 6 those related to the sarcomere model RDQ20-MF, Tab. 7 the parameters of the mechanical model and, finally, Tab. 8 contains the parameters associated with the circulation model. For the TTP06 model, we adopt the parameters reported in the original paper (for endocardial cells) [ten2006alternans], with the only difference that we rescale the intracellular calcium concentration by a factor of to get more physiological values [coppini2013late].

Variable Value Unit Variable Value Unit
Conductivity tensor Applied current
Calcium rescaling
Table 5: Parameters of the electrophysiological model.
Variable Value Unit Variable Value Unit
Regulatory units steady-state Crossbridge cycling
Regulatory units kinetics
Table 6: Parameters of the sarcomere model RDQ20-MF (for the definition of the parameters, see [regazzoni2020biophysically]).
Variable Value Unit Variable Value Unit
Constitutive law Boundary conditions
3 Tissue density
Table 7: Parameters of the mechanical model.
Variable Value Unit Variable Value Unit
External circulation Cardiac chambers
0.64 0.18
0.032116 0.07
0.32 0.05
0.035684 0.07
1.2 0.06
10.0 0.55
60.0 4.0
16.0 4.0
Cardiac valves
Table 8: Parameters of the circulation model.