Model of Anisotropic Reverse Cardiac Growth in Mechanical Dyssynchrony

Based on recent single-cell experiments showing that longitudinal myocyte stretch produces both parallel and serial addition of sarcomeres, we developed an anisotropic growth constitutive model with elastic myofiber stretch as the growth stimuli to simulate long-term changes in biventricular geometry associated with alterations in cardiac electromechanics. The constitutive model is developed based on the volumetric growth framework. In the model, local growth evolutions of the myocyte’s longitudinal and transverse directions are driven by the deviations of maximum elastic myofiber stretch over a cardiac cycle from its corresponding local homeostatic set point, but with different sensitivities. Local homeostatic set point is determined from a simulation with normal activation pattern. The growth constitutive model is coupled to an electromechanics model and calibrated based on both global and local ventricular geometrical changes associated with chronic left ventricular free wall pacing found in previous animal experiments. We show that the coupled electromechanics-growth model can quantitatively reproduce the following: (1) Thinning and thickening of the ventricular wall respectively at early and late activated regions and (2) Global left ventricular dilation as measured in experiments. These findings reinforce the role of elastic myofiber stretch as a growth stimulant at both cellular level and tissue-level.

Mechanical dyssynchrony of the heart such as left bundle branch block (LBBB) is caused by the asynchronous electrical activity of the left ventricle (LV) through slow conduction across the myocytes instead of through the fast conduction Purkinje system. Besides acute changes like QRS prolongation, reduction in wall motion and changes in blood flow 1 , abnormal contraction of the dyssynchronous heart can lead to long-term ventricular remodeling if left untreated. Patients with LBBB not only have an increased risk of developing cardiac disease 2 , such as hypertension, congestive heart failure, but also have a higher mortality rate [2][3][4] . Experimental studies using animal models have shown that ventricular pacing, which induces an asynchronous electrical activation and contraction pattern, produces both global and local LV remodeling in the form of ventricular dilation and asymmetrical hypertrophy, respectively 5 .
Computational modeling is increasingly used to simulate long-term changes associated with cardiac growth. Formulated based on the volumetric growth framework in which the deformation gradient tensor is multiplicatively decomposed into an elastic and a growth component 6,7 , cardiac growth constitutive models are developed to phenomenologically describe local changes in shape and size of the myocytes in response to local alterations of cardiac mechanics (i.e., stresses) and/or kinematics (i.e., strains). These constitutive models are usually coupled with a computational cardiac mechanics model to simulate the collective geometrical changes of the myocytes associated with alterations in loading conditions [8][9][10][11] in a realistic ventricular geometry. Existing computational cardiac growth models have largely focused on describing pathologies associated with global alterations in loading conditions, such as pressure and volume overload that produce concentric and eccentric hypertrophy, respectively 9,10,12,13 . Except for one computational model 14 , little has been done to simulate long-term changes associated with alterations of the electrical conduction pattern in the heart. In that model, longitudinal and transverse growth of the myocyte are controlled by the normal strains in those respective directions. Recent experiment data 15 , however, has shown that stretching the myocyte longitudinally can produce both series and parallel addition of the sarcomeres that result in longitudinal and transverse growth, respectively. RVFW and septum, with a higher value found in the latter. By contrast, only a small pre-systolic myofiber stretch was found at the LVFW in the Normal case compared to that in the septum of the Pacing case. Consequently, myofiber stretch in the septum + RVFW and LVFW of the Pacing case deviated positively and negatively from the homeostatic value in the Normal case. This heterogeneity in myofiber stretch λ f resulted in the evolution of growth scalars θ i 's towards θ i,max in the late activated septum/RV, and θ i,max in the early-activated LVFW, leading to long-term asymmetrical geometrical changes. pressure-volume relationship. Figure 1c shows the long-term changes in the LV and right ventricle (RV) pressure-volume (PV) loops for the Pacing case. There was no immediate substantial reduction in pump function in the Pacing case (0 month). Changes were, however, noticeable at 2 months in the form of a progressive LV dilation. Specifically, in the span of 8 months, LV end-diastolic volume (EDV) increased from 104.3 ml to 118.6 ml whereas end-systolic volume (ESV) increased from 55.65 ml to 70.3 ml. This led to a rightward shift in the LV PV loop that was accompanied by a slight change in ejection fraction (EF) from 48.7% to 48.3% in 8 months. On the other hand, the simulations also show long-term changes of the RV PV loops due largely to thickening of the septum. RV EDV slightly decreased from 104.9 ml to 103.7 ml whereas RV ESV increased from 53.6 ml to 57.1 ml at 8 months.
Long-term local geometrical changes. Figure 2 shows the long-term changes in geometry of a short-axis section taken at the mid-ventricular level and a long-axis section that are typically assessed in experimental and clinical studies. As shown in the figure, long-term changes in the geometry is highly asymmetrical in the biventricular unit in the Pacing case with radial wall thickening occurring at the late-activated septum and wall thinning occurring at the early-activated LVFW. Shown quantitatively in Fig. 3, the model predicted a monotonic increase in LV EDV (13.7%), septum wall volume (60.7%) and septum wall thickness (27.5%) over an 8 months period. Conversely, RVFW and LVFW thickness decreased by (12.3%) and (20.4%) over the same period, respectively. Wall volume of the LVFW and RVFW as well as RV EDV behaves non-monotonically and were roughly unchanged after 8 months.

Discussion
Motivated by a recent experiment showing that longitudinal stretching of myocyte produces both longitudinal and transverse growth via parallel and serial sarcomeric addition 15 , we have presented an anisotropic G&R constitutive model in which the changes in lengths of the tissue in 3 orthogonal material directions are driven locally by the deviation of maximum elastic myofiber stretch (over a cardiac cycle) from its corresponding homeostatic set point value. This model is coupled to an electromechanics modeling framework to simulate the long-term effects of asynchronous activation associated with LVFW pacing. We find that the application of elastic myofiber stretch as the sole mechanical driver of G&R in the constitutive model can reproduce quantitatively both local and global features of asymmetrical remodeling found in experiments of mechanical dyssynchrony as presented below.
In terms of acute behavior, our simulations show that in the presence of electromechanics alterations induced by LVFW pacing, a pre-stretch occurs at the late activated regions (Septum + RVFW) in the beginning of systole ( Fig. 1) that resulted in a higher maximum elastic myofiber stretch compared to the homeostatic set point value (during normal activation) found in those regions. The early activated region (LVFW), on the other hand, has a lower maximum elastic myofiber stretch when compared to its corresponding homeostatic set point value. These results are consistent with observations made in animal models of asynchronous activation 18,19 and LBBB patients 20 , where abnormal stretching of the tissue at the beginning of systole (i.e., pre-systolic stretching) were found at the late activated regions.
Using the alteration of elastic myofiber stretch as a stimulant for G&R in all 3 material directions, the model predictions, after appropriate calibration of parameters, are compared favorably to local and global measurements in canine experiments with similar chronic LVFW pacing protocol 5,17 . A quantitative comparison of the experimental results with model predictions are given in Table 1. In terms of local geometrical changes, the model predicted an increase in septum wall thickness by 18.5% (cf. 23 ± 12% in the experiments 5 ) and a decrease in LVFW thickness by 19.7% (cf. 17 ± 17% in the experiments 5 ) after 6 months of pacing. In terms of global geometrical changes, the model predicts an increase in LV EDV by 9% (cf. 7.4 ± 29% in the experiments) and LVFW + Septum wall volume by 9.5% (cf. 15 ± 17% in the experiments) for the same duration. As detailed in Table 1, the chronic features of LVFW pacing predicted by the model are also found in LBBB 21 , which produces mechanical dyssynchrony via an opposite activation pattern (i.e., septum is activated first followed by the LVFW). Interestingly, our model also predicted that the RV wall thickness is changed in LVFW pacing in response to the presence of a pre-systolic stretch (Fig. 1). Unfortunately, measurements of RVFW geometrical changes in response to LVFW pacing are, to the best of our understanding, not available for comparison. While there are a small number of studies on right branch bundle block (RBBB) which has similar (left to right) activation pattern as in LVFW pacing, the presence of RBBB is either superimposed on tachypacing-induced heart failure 22 or is confounded by other diseases such as Tetralogy of Fallot 23 .
Our study shows that it is necessary to impose different forward growth rate τ g and reverse growth rate τ rg in order to reproduce asymmetrical hypertrophy. Additional simulations using the same forward and reverse growth rates produce excessive thickening in the septum compared to LVFW case and are not able to predict the observed asymmetrical hypertrophy (See Appendix). While direct measurements of reverse growth rates on isolated cardiomyocytes are, to the best of our knowledge, not available, large animal and clinical studies have suggested that the reverse growth rates are likely to be different from the forward growth rates. For example, the study by Vernooy et al. 24 found that regional LV mass does not return to normal values (100%) during biventricular pacing. Similarly, the study by Gerdes et al. 25 also suggests that cardiac hypertrophy induced by aortocaval fistula in rats may not be completely reversible upon reversal of fistula.
As with most constitutive models describing cardiac G&R, the model described here is developed based on the volumetric growth framework 6 that largely focuses on changes in the morphology of the myocytes as they undergo hypertrophy or atrophy (as opposed to turnover of extracellular matrix constituents) during G&R 11 . One of the key differences between these G&R constitutive models is in the prescription of G&R stimuli in their formulation. Goktepe 9 et al. proposed using two different types of G&R stimuli, namely, elastic myofiber stretch www.nature.com/scientificreports www.nature.com/scientificreports/ for driving growth in the myofiber direction and Mandel stress for driving growth in directions transverse to the myofiber. They show that this formulation can reproduce features that are qualitatively in agreement with eccentric and concentric hypertrophy (i.e., volume and pressure overload). It is not clear, however, if it can predict asymmetrical remodeling features found in mechanical dyssynchrony. Subsequently, Kerckhoffs et al. 10 proposed using a single type of G&R stimuli, namely, strain in their model. While the G&R stimuli in the myofiber direction is similar to Goktepe et al. 9 in their model, they proposed using normal strain components transverse to the myofiber direction as stimuli for growth in the transverse directions. They showed that this formulation is able to reproduce features associated with global remodeling in pressure and volume overload as well as local remodeling associated with LBBB 14 . A detailed comparison of existing cardiac G&R constitutive models based on equibiaxial loading can be found in Witzenburg and Holmes 26 . In contrast to the G&R models described above, elastic myofiber stretch is prescribed here as the driver for growth in both the myofiber and transverse directions, but with different sensitivity (or rates). This assumption is consistent with a recent study that observed in real time both longitudinal elongation and lateral extension occurring in the myocyte in response to longitudinal stretching 15 . With appropriate calibration, we show that the prescription of a single growth stimuli based only on elastic myofiber stretch can quantitatively reproduce the largely local G&R features that are associated with mechanical dyssynchrony, which reinforce the theory that transverse growth maybe controlled, at least to some extent, by elastic myofiber stretch that previous models have not considered.
Although the anisotropic growth model predictions are generally consistent with experimental and clinical measurements, the model presented here, nonetheless, could be further improved with time. First, we have not considered sarcomeric parallel addition under lateral stretch as found in the experiment 15 because changes in myofiber (longitudinal) stretch is the most prominent asymmetrical feature of asynchronous activation. Moreover, without any direct measurements of forward and reverse growth rates associated individually with lateral and longitudinal stretch, it is difficult to separate parallel addition of sarcomere due to changes in lateral stretch from that due to longitudinal stretch to reproduce clinically observed features of asynchronous activation. Inclusion of lateral stretch as a growth stimuli in the model can be incorporated in future studies with the availability of more experimental data. Second, we did not consider any myocardial material property changes that may www.nature.com/scientificreports www.nature.com/scientificreports/ occur during the remodeling process. This limitation is shared by growth models developed based purely on the volumetric growth framework, which cannot directly accommodate for any changes in the mechanical properties. Structural changes in the tissue such as the reorientation of the collagen fibers were also not considered in our simulations. Nonetheless, the experiments did not report any change in the myocardial collagen fraction 5 , suggesting that extracelluar matrix remodeling is limited in mechanical dyssynchrony. Third, we have also assumed a constant preload and afterload in the entire simulation. Besides ventricular G&R, the circulatory system may also adapt during the long-term 27 . Specifically, models reflecting improvements in ventricular-arterial coupling during reverse remodeling 28 are expected to better capture long-term changes in CRT. We also did not account for any electrophysiological remodeling that may occur during electromechanics alterations. Besides changes in geometry, the myocyte's electrophysiological properties such as action potential duration 29 may be altered during the remodeling process. Last, our simulations did not consider the fact that asynchronous contraction can induce mitral regurgitation (MR), which may in turn confound the G&R process by inducing volume overload. Evidence that asynchronous contraction produces MR is supported by the observation that the severity of baseline MR is decreased in CRT patients 30,31 when their mechanical dyssynchrony is corrected. The model, which we have shown is able to reproduce features associated with volume (and pressure) overload (See Appendix) using fiber stretch as the growth stimuli, can be used in future studies to evaluate such possibilities.
In summary, we have presented an anisotropic growth model that simulates asymmetric growth in mechanical dyssynchrony. We show that by using the maximum elastic fiber stretch as the sole growth stimulant, the model can reproduce long-term features that agree quantitatively with the experimental and clinical findings. The growth model may be useful for optimizing the long-term pacing effects of cardiac resynchronization therapy in future studies.

Methods
Growth constitutive model. Here we let χ κ X t e g into an elastic component F e and an incompatible growth component F g . The growth deformation gradient F g was parameterized by the scalar growth multipliers θ f , θ s and θ n in the following tensorial form where f 0 , s 0 and n 0 are local myofiber, sheet, and sheet-normal directions in the reference configuration. The evolution of the growth multipliers was described based on deviations of a prescribed stimulant s j from its homeostatic value s i,h . Specifically, the evolution of growth multipliers was described by the following ordinary differential equation 16 .
where the subscript i = f, s, n denote the association with the myofiber, sheet, and sheet-normal directions. The growth constitutive model parameters τ g,i , τ rg,i , γ g,i and γ rg,i control the rates of growth and reverse growth, which may be different. The main purpose of the rate-limiting function is to restrict the evolution of the growth multipliers θ i to lie within some prescribed limits 7 (Fig. 4). Allowing the rate limiting function k i (θ I , s) to be prescribed individually in each i direction enables a broad spectrum of anisotropic growth deformation to be produced. Based on previous experimental findings 15 , maximum elastic myofiber stretch was prescribed as the growth stimuli in all 3 material directions (i.e., s i = λ f for i = f, s, n). The myofiber stretch was defined by where C ED denotes the right Cauchy-Green deformation tensor with respect to the end-diastolic configuration that was, in turn, defined by = C F F ED ED T ED with F ED = F e (F r,ED ) −1 and F r , ED representing the map from the unloaded configuration to the end-diastolic configuration. As depicted in Fig. 4, positive deviations from homeostatic stretch in the growth constitutive model will result in the growth multiplier θ i evolving towards θ max , i , (2019) 9:12670 | https://doi.org/10.1038/s41598-019-48670-8 www.nature.com/scientificreports www.nature.com/scientificreports/ which is typically an expanded configuration. Conversely, negative deviations from homeostatic stretch will result in the growth multiplier θ i evolving towards θ min , i , which is typically a shrunk configuration. electromechanics model. Cardiac electrophysiology. Cardiac electrical activity and its propagation was modeled using the modified Fitzhugh-Nagumo and monodomain equations 32 . Specifically, the spatio-temporal evolution of cardiac action potential φ  was described in the reference configurations by 0 is the anisotropic electrical conductivity tensor (with faster conduction along myofiber direction), I s is the constant electrical stimulus for prescribing local excitation initiation and pacing or both. Variables φ and r are dimensionless and r is a recovery variable. The excitation properties of cardiac tissue were defined using: Cardiac mechanics. An active stress formulation was used to describe the mechanical behavior of the cardiac tissue. In this formulation, the cardiac tissue constitutive relation was prescribed separately for passive and active mechanical behavior. Specifically, the second Piola-Kirchkoff or PK2 stress tensor S was additively decomposed into a passive component S p and an active componentS a , i.e.: The tissue passive mechanical behavior was described using the following Fung-type 35 strain energy function based on the Green-Lagrange elastic strain tensor E e :  www.nature.com/scientificreports www.nature.com/scientificreports/ where p is a Lagrange multiplier for enforcing the kinematic constraint = = F J det 1 e . The stress tensor was derived from the augmented strain energy function bŷ p e e A phenomenological active contraction model was used to describe contraction of the tissue upon excitation. In this model, the active (fiber-directed) stress tensor was defined by: a i nit a ff max 0 0 0 Magnitude of active force T was defined by a sarcomere length-dependent contractile force generation model 36,37 : where C a0 is the peak intracellular calcium concentration, T max is a scaling factor for the tissue contractility, and EC a50 is the length-dependent calcium sensitivity given by: a a max 50 0 0 with material constant B, prescribed maximum peak intracellular calcium concentration (C a0 ) max , instantaneous sarcomere length l, and sarcomere length at which no active tension develops l 0 . The instantaneous sarcomere length was defined based on the prescribed initial length of sarcomerel s0 , and was calculated using coupling electrophysiology and mechanics. The active contraction model was modified to incorporate a spatially heterogeneous activation initiation time t init (X) that accounts for asynchronous depolarization based on the action potential φ. The function ω in Eq. (13) is given by: that explicitly defines active contraction based on time since activation t sa (X) = t current − t init (X), with t current denoting the current time in the cardiac cycle as: init defining the local initiation time that couples cardiac electrophysiology and mechanics. In Eq. (15), t 0 is the prescribed time to maximum active tension and t r is the sarcomere length-dependent active tension relaxation time that is given by t r = ml + b with parameters m and b.
computational approximation. Discretizing time using an implicit backward Euler scheme, approximate solutions of the weak formulation for the acute electromechanics problem were obtained using the finite element (FE) method, where we seek to find u ∈ H 1 (Ω), p ∈ L 2 (Ω), P LV ∈ , P RV ∈ , c 1 ∈  3 , c 2 ∈  3 , ϕ ∈ H 1 (Ω), r ∈ H 0 (Ω), such that the below equations: : : are satisfied for all test functions δu ∈ H 1 (Ω), δp ∈ L 2 (Ω), δP LV ∈ , δP RV ∈ , δc 1 ∈  3 , δc 2 ∈  3 , δϕ ∈ H 1 (Ω), δr ∈ H 0 (Ω). Domains Ω LV,cav and Ω RV,cav denote LV and RV cavities, respectively. Spring constants k spring1 and k spring2 model are associated with the robin-type boundary conditions imposed at the epicardial surface ∂Ω epi and basal surface ∂Ω b , respectively. Variables c 1 and c 2 are Lagrange multipliers 38 to constrain rigid body translation and rotation, respectively. Variables P LV and P RV are Lagrange multipliers 39 constraining the LV and RV cavity volumes to their corresponding prescribed values V LV,p and V RV,p , respectively. Cavity volumes V LV,cav (u) and V RV,cav (u) were estimated using the divergence theorem 39 . It can be shown that the computed Lagrange multipliers P LV and P RV are equivalent to the LV and RV cavity pressures, respectively. Note that P LV and P RV are different from the Lagrange multiplier p which is a scalar field, i.e., p = p(C, t) to enforce local incompressibility and depends on the spatial location of the myocardium.
We used a Galerkin FE discretization with piecewise quadratic element for the spatial approximation of displacements u(X), continuous piecewise linear elements for hydrostatic pressure p and action potential φ. We used discontinuous Galerkin elements to approximate the recovery state variable r and fiber stretch λ f . A realistic human biventricular geometry, that was reconstructed from magnetic resonance images from a previous study 40 , was used for the computational study. Helix angle defining the myofiber direction f 0 was prescribed to vary linearly across the myocardial wall from 60° at the endocardium to −60° at the epicardium using a Laplace-Dirichlet rule-based algorithm 41 . Fiber directions f 0 as well as the orthogonal sheet direction s 0 and sheet-normal direction n 0 were prescribed at the quadrature points. A FE mesh with 11988 tetrahedral elements was used to discretize the biventricular geometry and FEniCS 42 was used for computational implementation. The time interval (0, T) was discretized into subintervals[t n , t n + 1 ], n = 0, 1, …, with t n = nΔt and Δt denoting a fixed time step.
Simulating the full cardiac cycle. Complete cardiac cycles with each having a duration of 750 ms were simulated using the electromechanics model. Specifically, passive filling was simulated by incrementally applying pressure to the LV and RV endocardial surfaces until their prescribed EDV was reached. Thereafter, systole was simulated by applying a current to different locations in the biventricular unit depending on the pacing protocol. The LV and RV cavity volumes were constrained to remain constant during the isovolumic contraction phase, which terminated when the LV and RV cavity pressures exceeded the prescribed aortic and pulmonary artery pressures, respectively. The subsequent ejection phase was simulated by coupling the LV and RV cavity volumes each to a 3-element Windkessel model (Fig. 5). Ejection phase in each cavity was terminated when the cavity outflow rate became less than 0.005 ml/min. The cavity volumes were then prescribed to remain constant to describe the isovolumic relaxation phase. The isovolumic relaxation phase was terminated based on prescribed pressures for initiation of diastolic filling, which was then passively prescribed until end of diastole.

Separation of time scale between growth and elastic deformation.
To integrate growth with the electromechanics model, we adopted a similar approach to separate the time scale as described in Lee et al. 8 . Specifically, since appreciable G&R can only be observed after large number of heartbeats, the growth tensor F g was treated as a constant that is independent of time in the electromechanics model and updated after each cardiac cycle. Approximate solutions of the weak formulation for the G&R problem were obtained using the FE method where we seek to find u ∈ H 1 (Ω), p ∈ L 2 (Ω), such that the below equation G T spring1 epi is satisfied for all test functions for all δu ∈ H 1 (Ω) and δp ∈ L 2 (Ω). The above equation describing the long-term evolution differs from the Eq. (17) describing short-term evolution in the following ways: there is no active stress in the long-term model; and the short-term model was solved for cavity pressures with prescribed volumes while the long-term model is solved for cavity volumes at prescribed (zero) cavity pressures. We used a Galerkin FE discretization with piecewise quadratic FE for the spatial approximation of displacements u(X). We used discontinuous Galerkin elements to approximate the growth parameters θ i s. The algorithm used to update the geometry (long-term problem) is described follows: • Obtain maximum elastic stretch λ f by solving the short time-scale problem over a cardiac cycle using Eq. • Update F g using θ + i j 1 • Solve for new geometry using Eq. (18).
Similar to our previous approach 8 , residual stresses induced by the incompatible growth tensor were removed after each growth cycle. As a consequence of this F = F e in Eq. (17a).    Table 2.