OUP user menu

Length-dependent tension in the failing heart and the efficacy of cardiac resynchronization therapy

Steven A. Niederer, Gernot Plank, Phani Chinchapatnam, Matthew Ginks, Pablo Lamata, Kawal S. Rhode, Christopher A. Rinaldi, Reza Razavi, Nicolas P. Smith
DOI: http://dx.doi.org/10.1093/cvr/cvq318 336-343 First published online: 14 October 2010


Aims Cardiac resynchronization therapy (CRT) has emerged as one of the few effective and safe treatments for heart failure. However, identifying patients that will benefit from CRT remains controversial. The dependence of CRT efficacy on organ and cellular scale mechanisms was investigated in a patient-specific computer model to identify novel patient selection criteria.

Methods and results A biophysically based patient-specific coupled electromechanics heart model has been developed which links the cellular and sub-cellular mechanisms which regulate cardiac function to the whole organ function observed clinically before and after CRT. A sensitivity analysis of the model identified lack of length dependence of tension regulation within the sarcomere as a significant contributor to the efficacy of CRT. Further simulation analysis demonstrated that in the whole heart, length-dependent tension development is key not only for the beat-to-beat regulation of stroke volume (Frank–Starling mechanism), but also the homogenization of tension development and strain.

Conclusions In individuals with effective Frank–Starling mechanism, the length dependence of tension facilitates the homogenization of stress and strain. This can result in synchronous contraction despite asynchronous electrical activation. In these individuals, synchronizing electrical activation through CRT may have minimal benefit.

  • Computer modelling
  • Heart failure
  • Contractile function
  • Cardiac resynchronization therapy
  • Ventricular function

1. Introduction

Heart failure (HF) is a costly, progressive disease characterized by a reduced quality of life and high mortality.1 Across a range of interventions, cardiac resynchronization therapy (CRT) has emerged as one of the few effective24 and safe5 treatments for HF. The goal of this intervention is to provide synchronous atriobiventricular stimulation, through pacemaker implantation, in order to resynchronize the electrical activation and mechanical contraction of the heart. However, despite improving life expectancy, HF symptoms, and quality of life within a significant patient cohort, identifying patients that will benefit from CRT still remains problematic, a situation which is exacerbated by the cost and invasiveness of the procedure.

The current selection criteria for CRT is limited to patients that are medically refractory, have symptomatic New York Heart Association (NYHA) Class III or IV limitation, electrocardiogram (ECG) QRS duration ≥130 ms, ejection fraction ≤35%, and left ventricular (LV) diameter ≥55 mm.4 Approximately 70% of patients who meet these criteria have been shown to benefit from CRT, and the criteria have the advantage that they can readily be evaluated non-invasively. However, a limitation of this selection method is evidenced by its inability to differentiate patients who show no signs of improvement following the procedure, so-called ‘non-responders'.6

This patient selection challenge has prompted multiple clinical trials to identify improved metrics in order to identify CRT responders. Initial studies from single centres showed promising results for echocardiographically derived metrics; however, multicentre trials have systematically demonstrated the inadequacy of these markers in identifying CRT responders compared with existing selection criteria.7 Furthermore, combinations of quantitative metrics proved of limited diagnostic value resulting in significant proportions (∼27%) of false-negatives.8 Hence, the current selection criteria based on ejection fraction, NYHA class, QRS duration, and ventricle diameter, although lacking precision, continue to direct patient selection.

A novel approach to address these issues and to improve our understanding of the mechanisms underpinning the CRT response is through the integration of information across the full range of clinically available diagnostic and medical imaging modalities. However, although this information has the capacity to provide a rich description of the electrical and mechanical state of the patient's heart, its inherent complexity also requires new approaches for its analysis. Specifically, parallel development of an analysis framework is required, to facilitate the integration of data and to exploit the full value of the combined information content that they represent.

In this study, we develop and apply such a framework in the form of a mathematical model which is constrained by invariable physical laws. The model links sub-cellular mechanisms through to whole organ function and is comprehensively validated using patient data collected pre- and post-CRT. Using this approach, we analyse the physiological mechanisms that regulate CRT response to identify novel, and potentially patient-specific, selection criteria on which to base treatment.

2. Methods

A comprehensive description of the model development and the diagnostic acquisition methods is given in the Supplementary material online. The development of the model is summarized in Figure 1. The model was developed in six stages. (1) The model geometry was fitted to magnetic resonance image (MRI) data. (2) A generic fibre field was then described based on a priori knowledge. (3) The patient's scar tissue was mapped onto the geometry using late enhancement gadolinium MRI. (4) The electrophysiology of the patient was fitted to EnSite™ endocardial activation data. (5) The mechanics of the heart were fitted to LV pressure and volume transients. (6) The model was validated against cine MRI deformation, EnSite™ activation maps during pacing, and changes in pressure transients following pacing.

Figure 1

Flow diagram showing the six steps in the model development and validation process.

2.1 Patient data

The study conforms to the principles outlined in the Declaration of Helsinki and was carried out as part of a local ethics committee-approved protocol with informed consent being obtained from the patient. We present data on a 60-year-old female with NYHA Class III heart failure despite optimal medical treatment. There was significant LV systolic dysfunction with an LV ejection fraction of 25%. The surface ECG demonstrated significant conduction disease with a left bundle branch block morphology and a QRS duration of 154 ms. Cardiac MRI had shown a small area of sub-endocardial posterolateral scar and severe LV dysfunction. A detailed description of the invasive study procedure is provided in the Supplementary material online.

2.2 Model geometry

The model geometry was fitted to MRI images with 1.5625 mm × 1.5625 mm × 10.0 mm resolution in the x-, y- and z-axes, respectively, with the z-axis running from apex to base. The LV and right ventricular (RV) endocardium surfaces and epicardium surfaces were segmented and fitted by three bicubic Hermite surfaces (LV/RV endocardium and epicardium). The resulting surfaces were then combined into a tricubic Hermite mesh describing the patient's heart geometry (Figure 2).

Figure 2

Schematic representation showing the fibre orientation in the heart model. The cubes show the transmural variation of fibre angle, as shown by the blue arrows and sheet angle, as shown by the pink planes at (A) lateral LV free wall base, (B) lateral LV free wall apex, (C) posterior LV wall, (D) anterior LV wall, (E) septum, and (F) RV free wall.

2.3 Model fibres

The fibre orientation cannot presently be obtained in vivo from patients. Although the fibre orientation is known to change in HF, the degree of change and the consistency of these changes are not known. Hence a generic, healthy fibre field was determined from animal studies911 augmented with human data12 where available. This resulted in the definition of fibre and sheet orientation that varied transmurally, between apex and base, and between the LV, RV, and septum. The fibre orientation and model geometry are summarized in Figure 2.

2.4 Model scar geometry

Late enhancement gadolinium MRI studies identified regions of myocardial infarction in the heart. The data were used to separate the heart into two regions: viable and scar. The scar tissue was described by tricubic Hermite volumes and were used to describe the scar regions in both the mechanical and electrical models.

2.5 Modelling electrical activation

Electrical activation was simulated using the monodomain equations. Scar tissue was simulated as isotropic low-conduction regions. The timing and location of stimulation for baseline conditions were determined from the baseline EnSite™ activation data.13 The conduction parameters in the scar and tissue were fitted to the EnSite™ activation maps.

2.6 Modelling mechanical deformation

Mechanical deformation was simulated using the finite elasticity equations to ensure that mass, momentum, and work are conserved. The passive material properties of the heart were modelled using a transversely isotropic material law, aligned to the microstructure of the myocardium. Scar was simulated as an isotropic region with increased stiffness. The parameters for the material constitutive equation were fitted using the pressure and volume relationship during atrial contraction, when the myocardium is assumed to be quiescent.

A Windkessel model provided the pressure–volume relationship for the heart model during ejection. The model parameters are fitted to the pressure–volume relationship during ejection. During isovolumetric contraction and relaxation, the volume of the LV and RV cavities were held constant. The diastolic period was of limited interest for this study and was approximated by defining a simple pressure transient that was linearly scaled over the duration of diastole.

Active tension was simulated by a simple model which was a function of electrical activation time. The model has length-dependent peak tension and length-dependent rate of activation, which share a common nonlinear length-dependent function. Model parameters correspond to the rate of relaxation, tension transient duration, rate of activation, length-dependent rate of activation scalar, peak tension, minimum length of tension generation, and degree of length dependence. The model parameters were fitted to the baseline pressure transient.

2.7 Software and hardware

The geometry was created and visualized using the interactive 3D model development environment CMGUI (http://www.cmiss.org/cmgui). The electrophysiology mesh was created using the medical image meshing software package Tarantula (http://www.meshing.at). Electrophysiology was solved using the cardiac reaction diffusion simulation software CARP (http://carp.meduni-graz.at). Mechanics were solved using large deformation functionality of the simulation software package CMISS (http://www.cmiss.org). Simulations were performed using 128–1024 cores on Hector (www.hector.ac.uk) or using 4–128 cores on the Oxford Supercomputing Centre (www.osc.ox.ac.uk) machines.

3. Results

3.1 Validation patient-specific heart model

Figure 3 summarizes the fitting and validation of the mechanical, electrical, and coupled model. The conduction properties of the electrophysiological model were adjusted to match baseline simulations with the corresponding EnSite™ endocardial activation maps (Figure 3A). The electrophysiological model was then validated by comparing model and EnSite™ activation maps in the presence of LV endocardial pacing (Figure 3D). Quantitative comparisons between EnSite™ data and the model results were not possible because of the discrepancies in the endocardial geometries from MRI and EnSite™, shown as the overlaid elliptical mesh and spherical colour-mapped surface, respectively, in the top rows of Figure 3A and D. However, a qualitative comparison shows that the model provides an accurate representation of the available data and that it is then able to predict the activation pattern following pacing.

Figure 3

(A) Comparison of simulated and measured baseline/sinus rhythm endocardial activation times. Top panel shows the EnStie™ activation maps on the geometry derived from the EnSite™ system with the MRI-derived endocardium geometry represented by the transparent mesh. Bottom panel shows the activation times calculated by the model on the LV endocardium with conduction parameters adjusted to match EnStie™ data. (B) The LV pressure transient measured by catheter (yellow line) or simulated by the model (red line). (C) The normalized volume for each MRI time point (yellow squares) and the normalized volume transient predicted by the model (red line). (D) Comparison of the EnSite™ activation maps and model simulations for LV pacing. (E) The model and (F) the catheter dP/dt transients for sinus rhythm (red line) and LV pacing (yellow line). (G) Comparison of the wall motion predicted by the model (red line) with the cine MRI.

The active contraction model was fitted to the percentage change in maximum rate of LV cavity pressure development (dP/dtmax) upon LV pacing and the pressure transient (Figure 3B), using the steepest decent method and parameter sweeps. The final coupled electromechanics model captured the volume transient (Figure 3C), with the difference in ejection fraction of only 3.9% between the model and the measurements, which lies within the range of inter-operator variability of cine MRI analysis.14 Figure 3E and F demonstrates similar changes in the dP/dt transient between sinus rhythm (red) and LV pacing (yellow) for the model and the clinical data, respectively. A qualitative comparison is made due to variation in the baseline pressure transient over the course of the procedure and the differences in the pacing frequency and heart rate between paced and the baseline transients. The mechanics model was further validated by comparing the model wall motion directly to cine MRI data. Figure 3G shows that the model is in close agreement with the measured wall motion (Figure 3G). When combined, the validation of the electrical and coupled model provided confidence that the model is operating with physiologically plausible boundary conditions and parameters.

3.2 Factors contributing to CRT efficacy

The proposed model replicates baseline and pacing conditions in the heart. The model spans from cellular excitation and tension generation through to whole organ mechanics and electrical activation. To predict the properties of the heart that impact the outcome of CRT, the sensitivity of the efficacy of CRT, defined as the percentage change in dP/dtmax upon LV pacing, to 16 of the model parameters was calculated. The sensitivity is equal to the derivative of the normalized efficacy of CRT with respect to the normalized value of each parameter. This means that if a percentage change in a parameter causes a greater percentage change in the efficacy of CRT, the sensitivity will be higher than 1. The 16 parameters impact the passive material properties, active contraction, electrical conduction, and boundary conditions in the model. Table 1 summarizes the factors and the sensitivity of efficacy of CRT to each of the parameters.

View this table:
Table 1

Sensitivity of efficacy of CRT to model parameters

DffElectrical conduction in the fibre direction0.608
DxxElectrical conduction in the cross-fibre direction0.282
C0Scalar of passive material stiffness0.065
αScalar of passive material anisotropy0.043
RaWindkessel model aortic resistance0.077
RcWindkessel model circulatory resistance0.021
CcWindkessel model circulatory compliance0.117
ApAortic pressure0.373
EDPEnd diastolic pressure0.149
troBase active tension development rate0.633
a4Scalar for length-dependent rate of tension development0.240
a6Level of length-dependent tension1.640
a7Stretch at which active tension reaches zero3.226
tdRate of relaxation0.021
ToPeak tension0.339
tmaxDuration of active tension transient0.058

This analysis reveals two significant parameters in determining the efficacy of pacing in the model, corresponding to the length dependence of tension15 (a6) and the minimum length of tension generation (a7). As the length of minimum tension generation was regulated by sarcomere structure, and is likely not to change during HF, this result highlighted the length dependence of tension as a fundamental mechanism underpinning CRT efficacy.

3.3 Length dependence of tension and CRT

To further investigate how this cellular-scale mechanism alters whole organ function, the regional stress, strain, and rate of tension development were calculated. Figure 4H shows that, within the model, the efficacy of CRT increased with decreasing the effect of the length-dependent tension mechanism. To further elucidate this observation, the average active tension and strain in regions of early and late activation were analysed.

Figure 4

Active tension, rate of tension development, and regional stretch in the presence of baseline, amplified, and attenuated length-dependent tension regulation. (A) The three regions: (1) septum (red), (2) posterior lateral wall (yellow), and (3) anterior lateral wall (blue). All plot lines colours correspond to these regions. For all line plots, dashed and solid lines correspond to paced and baseline stimulation. (BD) The active tension, (EG) the rate of change of active tension, and (IK) the fibre stretch. (H) The percentage change in dP/dtmat between the baseline and LV paced case in the presence of −50, 100, and +50% length-dependent tension regulation.

Figure 4BD shows the development of active tension in the septum (early activation, region 1, Figure 4A) and LV free wall (late activation, regions 2 and 3, Figure 4A), for the baseline (sinus rhythm) and LV-paced cases, with a6 (length-dependent tension) at 50, 100, and 150%, respectively. Similarly, Figure 4IK shows the fibre stretch in the fibre direction for these two regions for different levels of the length-dependent tension. Figure 4H shows the impact of varying the length dependence of tension effect on the efficacy of CRT.

4. Discussion

Through a sensitivity analysis, we have identified the length dependence of tension as a significant contributor to the efficacy of CRT that potentially may vary between patients. Furthermore, the model analysis enables us to demonstrate that effective length-dependent modulation of tension causes the homogenization of tension generation despite an asynchronous activation field, such that when length-dependent regulation is present, limited mechanical benefit for pump function can be provided by CRT.

4.1 Length-dependent tension

Mechanical synchrony is the result of the combined effects of electrical activation times and the regulation of tension development. As tension development is secondary to electrical activation, it can either amplify or dampen any asynchrony induced by the electrical activation. The simulations, in matching clinical measurements, have asynchronous electrical activation times, with the septum activating early and the LV free wall late (Figure 3A). The personalized human model developed in this study demonstrates that the length dependence of tension is a major regulator of mechanical synchrony in the presence of electrical asynchrony. This is shown by the decreasing synchrony in strain (Figure 4IK), tension development (Figure 4BD), and both the rate of tension development and the time that the rate of tension development peaks (tpeak) (Figure 4EG) with decreasing length-dependent tension regulation.

Figure 4E shows that in the presence of effective length-dependent tension regulation, tpeak can be synchronized (as seen by the converging peaks of the solid lines in Figure 4E), despite asynchronous electrical activation. In contrast, attenuated length-dependent tension regulation results in a dispersion of tpeak with asynchronous electrical activation (Figure 4G). The introduction of a pacing site in the LV free wall results in a significant increase in electrical synchrony (Figure 3A and D). In the attenuated length-dependent tension case, synchronizing activation also synchronizes the tpeak (Figure 4G); however, in the effective length-dependent tension case, tpeak was already synchronized, prior to pacing, so there is limited improvement in tpeak synchrony (Figure 4E). The significant and nominal change in the dispersion of tpeak with attenuated or effective length-dependent tension, respectively, is manifested in the significant differences in the change in dP/dtmax upon pacing between these cases (Figure 4H).

The change in dP/dtmax, a common metric for CRT responders,1618 was used as a metric for the efficacy of CRT in model simulations. As the value of dP/dtmax is intrinsically a function of the integrated rate of active tension development in the myocardium, it is no surprise that significant decrease in the dispersion of tpeak causes significant improvement in dP/dtmax. The dispersion of tpeak is thus a potential indicator that a patient will improve upon receiving CRT. However, the majority of clinically derived indices for identifying CRT responders are based on deformation or electrical activation.3,4,6,7 These measures of cardiac function are readily available in the clinic but they are unable to characterize local active tension development, which cannot be measured directly, limiting the use of a dispersion of tpeak metric to identifying CRT responders.

The Frank–Starling mechanism can be taken as an organ-scale surrogate for tissue-scale length-dependent tension regulation. However, studies in muscle preparations from failing human hearts are inconsistent, with the length dependence of tension generation reported as absent,19,20 attenuated,21,22 or unchanged.2325 It is likely that this variation can partially be explained by diversity in the HF patient population, with some patients having effective length-dependent tension development, whereas others have attenuated or absent length-dependent tension regulation. If the length dependence of tension varies across the HF patient population, then identifying patients with attenuated length-dependent tension could facilitate CRT selection.

A central difficulty of proposing such an identifying metric is that the effects of the Frank–Starling mechanism are not routinely measured in vivo in HF patients, with recent studies using excised hearts,23 papillary muscles,20 and trabeculae.25,26 However, the efficacy of the Frank–Starling mechanism is characterized in patients suffering from hypotension by altering the preload through a fluid challenge27 or passive leg raising28 and recording of the subsequent change in cardiac output. These could be performed with a non-invasive cardiac MRI study and phase contrast aortic through-plane cine velocity and therefore flow measurement (which has been validated in a number of studies for its accuracy and reproducibility29) of cardiac output before and immediately after an intravenous fluid challenge or leg elevation while in the MRI scanner. The model prediction that for individuals with attenuated length-dependent tension regulation, the effectiveness of CRT will be improved, combined with a viable clinical measurement of the effectiveness of the Frank–Starling mechanism, offers a potentially promising new metric to assist in identifying responders to CRT.

4.2 Generalization

The current model requires a broad range of clinical data to constrain and validate the model. To maintain consistency, the data have been collected from a single patient. However, from the detailed analysis of this specific model, we can distil the role of length-dependent tension as a general phenomenon by considering the case of a strand of cardiac muscle that is stimulated at one end, electrical activation propagates along the muscle, resulting in dyssynchronous early and late-contracting regions. Figure 5 shows idealized schematic model tension transients in each of these regions in the presence of length-dependent tension (solid lines) and in the absence of length-dependent tension (dashed lines).

Figure 5

Schematic diagram showing the development of tension in the early-activated and late-activated regions of a papillary muscle in the presence of (solid line) and in the absence of length-dependent tension generation (dashed line). The blue arrows highlight changes in the two transients.

In the presence of length-dependent tension, the early-contracting region (red lines) will initially shorten as the late-activating region is quiescent, thus elongating the late-activating region and reducing the rate of tension development in the early-activating region. When the late-activating region contracts, it will do so from a longer length, increasing its rate of tension development, as indicated by arrows between the solid yellow and dashed lines in Figure 5. When the late-contracting region shortens, it will elongate or reduce the rate of shortening of the early-contracting region increasing its rate of tension development, thus tension development in the early-activated region is delayed compared with the case where there is no length-dependent tension regulation, as indicated by the arrows between the solid and dashed red lines in Figure 5. In this way, the rate of tension development and the time of peak tension will become homogenous across the muscle despite heterogeneous activation. In the presence of effective length-dependent tension development, the homogenization of electrical activation has limited benefit as the development of tension is already synchronous prior to CRT.

4.3 Limitations

The model proposed in this work inevitably provides only a simplified representation of the patient's heart. There are many physiological aspects of the heart that are necessarily absent or simplified within the model as they cannot be readily characterized (for more details, see Supplementary material online), even with the extensive data sets available (Figure 1). For this reason, both the model results and conclusions we have reached must be considered within the context of a number of assumptions.

The clinical data and model results only characterize the acute changes that occur immediately following pacing and may not reflect the long-term outcome of these patients. However, patients with effective length-dependent tension will see nominal changes in cardiac function and hence are unlikely to improve over time. Thus, although this is an acute study, it is likely that the conclusions have long-term implications.

The model has identified a single contributing factor to the efficacy of CRT. However, the complexity and multiple confounding factors that regulate heart function during HF mean that the length dependence of tension alone is unlikely to be the sole indicator of CRT efficacy. Further studies will be needed to provide a comprehensive suite of metrics to improve on the sensitivity and accuracy of CRT patient selection.

4.4 Conclusion

Comprehensively validated computational models provide a promising method to rationally integrate data from multiple diagnostic and imaging modalities for developing patient-specific cardiac models to extract the maximum value from these data sets.

Developing and applying such a model, we have identified the length dependence of tension as a significant regulator of CRT efficacy. The model predicts that a patient with dyssynchronous electrical activation but effective length-dependent tension regulation at the cellular scale is less likely to respond to CRT treatment than a patient with attenuated or no length dependence of tension.

Supplementary material

Supplementary material is available at Cardiovascular Research online.

Conflicts of interest: M.G. receives a St Jude Medical educational grant and C.A.R. receives St Jude educational support and is on advisory boards for St Jude Medical and Medtronic.


This work was supported by United Kingdom Engineering and Physical Sciences Research Council for support through grants EP/F043929/1 and EP/F059361/1.


View Abstract