In silico screening method for non-responders to cardiac resynchronization therapy in patients with heart failure: a pilot study

Abstract

Background

Cardiac resynchronization therapy (CRT) is an effective treatment option for patients with heart failure (HF) and left ventricular (LV) dyssynchrony. However, the problem of some patients not responding to CRT remains unresolved. This study aimed to propose a novel in silico method for CRT simulation.

Methods

Three-dimensional heart geometry was constructed from computed tomography images. The finite element method was used to elucidate the electric wave propagation in the heart. The electric excitation and mechanical contraction were coupled with vascular hemodynamics by the lumped parameter model. The model parameters for three-dimensional (3D) heart and vascular mechanics were estimated by matching computed variables with measured physiological parameters. CRT effects were simulated in a patient with HF and left bundle branch block (LBBB). LV end-diastolic (LVEDV) and end-systolic volumes (LVESV), LV ejection fraction (LVEF), and CRT responsiveness measured from the in silico simulation model were compared with those from clinical observation. A CRT responder was defined as absolute increase in LVEF ≥ 5% or relative increase in LVEF ≥ 15%.

Results

A 68-year-old female with nonischemic HF and LBBB was retrospectively included. The in silico CRT simulation modeling revealed that changes in LVEDV, LVESV, and LVEF by CRT were from 174 to 173 mL, 116 to 104 mL, and 33 to 40%, respectively. Absolute and relative ΔLVEF were 7% and 18%, respectively, signifying a CRT responder. In clinical observation, echocardiography showed that changes in LVEDV, LVESV, and LVEF by CRT were from 162 to 119 mL, 114 to 69 mL, and 29 to 42%, respectively. Absolute and relative ΔLVESV were 13% and 31%, respectively, also signifying a CRT responder. CRT responsiveness from the in silico CRT simulation model was concordant with that in the clinical observation.

Conclusion

This in silico CRT simulation method is a feasible technique to screen for CRT non-responders in patients with HF and LBBB.

INTRODUCTION

Cardiac resynchronization therapy (CRT) is one of the established treatment options for patients with heart failure (HF) and left ventricular (LV) dyssynchrony [1]. Patients with CRT benefit from the effects of electric and mechanical synchronization of their ventricular muscle and LV function is ultimately improved. However, CRT is not effective in approximately 1/3 of the patients and this CRT non-responder issue has not been resolved to date [2, 3]. Therefore, for achieving a more ubiquitous performance of CRT in patients with HF, patient-specific simulation methods for screening CRT non-responders need to be established. In view of the many difficulties associated with performing an experimental or clinical study using a living heart, computer simulation can be an alternative method to understand the physiological mechanisms of CRT. Hitherto, computer simulation of cardiac electrophysiology and mechanics has been considered a useful tool for the investigation of the pathophysiology of heart disease and for the development of novel therapeutic clinical techniques [4]. For example, a multi-scale heart modeling was developed by integrating cell, tissue, and organ scales to investigate the electrophysiological behavior of the heart [5]. However, most of the computational studies on the mechanism of heart disease have been performed not on a patient-specific heart model, but on a typical heart model. For CRT simulation, a patient-specific heart model is essential. In this study, we propose a patient-specific cardiac electromechanical model of pre- and post-CRT treatment to screen for CRT non-responders in patients with HF and left bundle branch block (LBBB).

Methods

Electromechanical heart model

To construct an integrated CRT simulation model, we combined the 3-dimensional (3D) image-based electromechanical model of a patient with LBBB with the lumped model of the circulatory system and electric pacing device [6]. A schematic diagram of the integrated model is shown in Fig. 1. To simulate the electric and mechanical behavior of the heart during cardiac tissue excitation and contraction, we used the electromechanical model that has two dynamic components, electrical and mechanical, as described previously [7]. Physiologically, as an electric wave propagates through the heart, the depolarization of each myocyte initiates the release of calcium (Ca²⁺) from the stores in the sarcoplasmic reticulum, which is followed by the binding of Ca²⁺ to troponin C and cross-bridge cycling. The cross-bridge cycling forms the basis for sliding of myosin (thick) filaments relative to actin (thin) filaments and the development of active tension in the cells, resulting in contraction of the ventricles.

Fig. 1

The schema of the finite element ventricular electromechanical model coupled with the circulatory and CRT models. C_a = systemic artery compliance, C_pa = pulmonary artery compliance, C_pv = pulmonary vein compliance, CRT = cardiac resynchronization therapy, C_v = systemic vein compliance, ECG = electrocardiography, R_a = systemic artery resistance, R_lo = left ventricular outflow resistance, R_pa = pulmonary artery resistance, R_pv = pulmonary vein resistance, R_v = systemic vein resistance

The electrical component of the model simulates the propagation of a wave of transmembrane potential by solving the monodomain equations on the electrical mesh as shown in Eq. 1.

$$ \frac{{\partial V_{m} }}{\partial t} = - \left( {I_{{{\text{ion}}}} + I_{{{\text{app}}}} } \right)\frac{1}{{C_{m} }} + \left( {\nabla \cdot D\nabla V_{m} } \right) $$

(1)

Here, t, V_m, I_ion, I_app, C_m, and D are time, membrane voltage, current across the cellular membrane, stimulation current, membrane capacitance, and electric diffusion coefficient, respectively. To solve the equation, we spatially discretized the equation by using the finite element method and time derivative in the equation was approximated by forward Euler method. More detailed explanation for the solution of the equation is shown in our previous paper [8].

A calcium transient serves as an input to the cell myofilament model representing the generation of active tension within each myocyte, where a set of ordinary differential equations and algebraic equations describes Ca²⁺ binding to troponin C, cooperativity between the regulatory proteins, and cross-bridge cycling. Contraction of the ventricles results from the active tension generation represented by the model of myofilament dynamics. Contraction is described by the equations of passive cardiac mechanics, with the myocardium being an orthotropic, hyperelastic, and nearly incompressible material with passive properties defined by an exponential strain–energy function as in our previous paper [9]. Simultaneous solution of the myofilament model equations with those representing passive cardiac mechanics on the mechanical mesh constitutes the simulation of cardiac contraction. The interaction between the electric model and mechanical one is depicted in Fig. 2. Also circulatory dynamics represented by vascular resistances (hereafter ‘R’) and capacitances (hereafter ‘C’) was combined with the heart mechanics.

Fig. 2

The schema of the cardiac electromechanical model. Here, the cellular excitation–contraction is incorporated into the tissue mechanical model and vascular hemodynamics implemented in the lumped parameter model interacting with cardiac electromechanical contraction. C_a = systemic artery compliance, C_pa = pulmonary artery compliance, C_pv = pulmonary vein compliance, CRT = cardiac resynchronization therapy, C_v = systemic vein compliance, R_a = systemic artery resistance, R_lo = left ventricular outflow resistance, R_pa = pulmonary artery resistance, R_pv = pulmonary vein resistance, R_v = systemic vein resistance

Patient-specific heart model for in silico CRT simulation

The numerical methods for solving the equations of the electromechanical model have been described previously [10]. To establish LBBB model in the heart model, the left part of the Purkinje fiber networks was blocked as shown in Fig. 3. Here, the network was from our previous work [11]. A patient-specific 3D ventricular geometry was reconstructed from the computed tomography (CT) images by using a commercial segmentation software package (Aquarius iNtuition Version 4.4.11 TeraRecon, Inc, San Mateo, CA, USA). To estimate the parameters in the model, the clinically measured values of the physiological data were used. The procedure for the estimation is as follows:

Fig. 3

Electric potentials in Purkinje fiber networks. Here, color indicates electric potentials. Red means electrically activated whereas blue means resting state. a. In normal ventricular conduction model, most parts of the ventricles activate almost simultaneously. b. In the LBBB model, most of the left part of the Purkinje fiber network remains electrically inactivated. LBBB = left bundle branch block

Step 1. Measure the physiological data of a patient in pre-CRT state: systolic (SBP) and diastolic blood pressures (DBP), heart rate (HR), and cardiac output (CO).

Step 2. Construct the ventricle model with 3D geometry from the CT images of the patient.

Step 3. Add the Purkinje fiber network to the 3D ventricular geometry.

Step 4. Assume the resistance (R) and capacitance (C) values of the vascular system.

Step 5. Couple the circulation model having the assumed parameters of R and C with the 3D ventricular mechanics model.

Step 6. Perform the simulation and obtain the computed SBP, DBP and CO of the patient.

Step 7. Calculate the errors of the computed SBP, DBP and CO compared with the clinically measured ones.

Step 8. If the errors are greater than 5%, go to Step 2 with newly assumed R and C parameters. If the errors are less than 5%, go to Step 9 since these parameters satisfy the physiological state of the patient.

Step 9. Based on the parameters, we perform CRT simulation by virtually pacing two areas on the left and right ventricular surface.

Step 10. From the computed results of CRT simulation, we determine whether the patient is a responder or non-responder.

Using the algorithm, we can match the model parameters in the simulation with clinically measured physiological properties. Then, CRT simulations based on the parameters are performed to assess whether the patient is a non-responder to CRT.

Comparison of changes in LV volumes in the in silico CRT simulation model with those from clinical observation

The study protocol was approved by the institutional ethics board of Severance Hospital, Yonsei University College of Medicine (IRB number: 1–2019-0077 and 1–2020-0008). A 68-year-old female with nonischemic HF with LBBB was retrospectively included (Table 1). Electrocardiography (ECG), echocardiography, cardiac CT and magnetic resonance (MRI) were performed before implantation of CRT device. Because LV function had not improved despite optimal medical therapy for HF for 3 months, a CRT-defibrillator was implanted in this patient. Electrocardiography and echocardiography were performed 6 months after CRT for assessing CRT responsiveness. LV end-diastolic (LVEDV) and end-systolic volumes (LVESV), and LV ejection fraction (LVEF) were measured by echocardiography before and 6 months after CRT implantation. LVEF was defined as (LVEDV − LVESV)/LVEDV × 100. A CRT responder was defined as absolute increase in LVEF ≥ 5% and relative increase in LVEF ≥ 15% at 6 months after CRT [12]. Absolute ΔLVEF was defined as LVEF_post-CRT − LVEF_pre-CRT. Relative ΔLVEF was defined as (LVEF_post-CRT − LVEF_pre-CRT)/LVEF_pre-CRT × 100. LVEDV, LVESV, LVEF, absolute and relative ΔLVEF, and CRT responsiveness measured from the in silico simulation model were compared with those found in clinical observation.

Table 1 Baseline characteristics, ECG, and echocardiographic data of the patient

Results

The patient-specific CRT model was constructed from the CT images of the patient. For the validation of the LBBB model, we compared the computed ECG data with the clinically measured values. As shown in Fig. 4, the computed ECGs show the typical pattern of LBBB [14]. The baseline characteristics, ECG, and echocardiographic data of the patients are shown in Table 1. Based on the measured SBP, DBP, and HR, we iteratively computed the 3D heart model to find the simulation parameters (R and C values) that can reproduce the measured CO. Figure 5 shows an activation map of both ventricles at 0, 80, and 140 ms under the LBBB and CRT conditions. In the case of LBBB, the right ventricle (RV) was activated before the LV, which shows that the model of LBBB was correctly implemented. When CRT was performed in the model, however, the RV and LV were activated almost simultaneously, which demonstrates that this patient model responded to the CRT protocol. Figure 6 shows the changes of the arterial and LV pressures as well as the LV volume with time under the LBBB and CRT conditions. CRT increased SBP from 120 to 124 mmHg and DBP from 68 to 69 mmHg. The in silico simulation modeling revealed that changes in LVEDV, LVESV, and LVEF by CRT were from 174 to 173 mL, 116 to 104 mL, and 33 to 40%, respectively. From in silico modeling, absolute and relative ΔLVEF were 7% and 18%, respectively, and we took that to mean a CRT responder (Table 1). In clinical observation, echocardiography showed that changes in LVEDV, LVESV, and LVEF by CRT were from 162 to 119 mL, 114 to 69 mL, and 29 to 42%, respectively. Absolute and relative ΔLVESV were 13% and 31%, respectively, and that means that the patient was a CRT responder (Table 1). Therefore, CRT responsiveness from the in silico CRT simulation model was concordant with that in the clinical observation.

Fig. 4

Computed ECG in the patient in the LBBB (a) and CRT pacing (b) conditions. CRT = cardiac resynchronization therapy, ECG = electrocardiogram, LBBB = left bundle branch block

Fig. 5

The activation maps of ventricles under the LBBB (a) and CRT (b) conditions. While the RV was activated before the LV in the LBBB condition, both ventricles were activated almost simultaneously in the CRT condition. CRT = cardiac resynchronization therapy, LBBB = left bundle branch block, LV = left ventricle, RV = right ventricle, V_m = membrane voltage

Fig. 6

Changes in pressure and LV volume with time under the LBBB (a) and CRT (b) conditions. CRT increased systolic LV pressure from 120 to 124 mmHg and decreased LV end-systolic volume from 116 to 104 mL, that resulted in increased LVEF from 33 to 40%. CRT = cardiac resynchronization therapy, LBBB = left bundle branch block, LV = left ventricle, LVEF = left ventricular ejection fraction

Discussion

Main findings

In this study, we constructed the patient-specific three-dimensional electromechanical heart model with LBBB, simulated CRT on the model, and compared the cardiac parameters at pre- and post-CRT to test the applicability of the model to CRT simulation. The simulated ECGs were similar to those reported in the literature under both pre- and post-CRT conditions. The simulated electrical wave propagations in the ventricles exhibited the sequence of activation typically observed under both LBBB and CRT conditions. An improvement in LVEF was observed after CRT simulation.

Prior CRT simulation models

Because CRT is not effective in about one-third of the patients with HF and LBBB [12, 13], screening for CRT non-responders before CRT implantation is an important strategy. Computer simulation can be a useful tool for testing the effect of CRT for an individual patient before treatment. The factors affecting the efficacy of CRT include location of the LV and RV leads, atrioventricular and interventricular delay and underlying pathology [15]. Moreover, the optimal location of ventricular pacing varies in each individual patient [16, 17]. The pacing locations and time delays can be tested using computer simulation. Isotani et al. [18] tested four different lead positions and examined the change in the maximum time derivative of ventricular pressure (%ΔdP/dt_max) for each location. They hypothesized that patients with insignificant gain in %ΔdP/dt_max even at the best lead position would be non-responders. Placing the ventricular leads on the scar region would have minimal effect on the efficacy of resynchronization. The locations and size of the scar can also be examined by constructing a 3D heart model using images from modalities such as MRI. The accuracy of the heart model for an individual patient would be critical in the interpretation of the simulation results. Okada et al. [19] constructed a heart model and simulated CRT, but they modeled the Purkinje fiber network as a thin layer on the endocardial surface. Our 3D electromechanical heart model has all the necessary components for individualized CRT simulation such as the Purkinje fiber network, fiber orientation, mechanical contraction/relaxation, and systemic circulation. The CRT simulation performed on one patient model in this study exhibited the possibility of utilizing this model for screening non-responders to CRT before treatment. Appropriate simulation protocols or biomarkers would have to be developed for the screening. Okada et al. [20] observed that dP/dt_max showed the best correlation with clinically observed improvement in LVEF. They also observed that the change in the duration of the QRS complex did not correlate with the clinical improvement in LVEF.

Criteria of CRT responders

It still remains controversial how and when to define CRT response. Several clinical, echocardiographic, and hemodynamic criteria for CRT response have been suggested. The echocardiographic CRT response criteria were the following: decrease in LVESV > 15%, decrease in LVEDV > 15%, absolute increase in LVEF ≥ 5%, relative increase in LVEF ≥ 15%, and decrease in functional mitral regurgitation ≥ 1 grade. “Decrease in LVESV > 15% at 6 months after CRT implantation” is generally considered as the clinically useful criterion for CRT response. However, “Decrease in LVESV > 15%” might not be appropriate as a CRT response criterion in the in silico CRT simulation model, because “decrease in LV volume” is one of the markers of long-term reverse remodeling of the LV and long-term reverse remodeling is not taken into account in the in silico CRT simulation model. In contrast, the in silico CRT simulation model adequately reflects the acute changes in LV contractility induced by CRT. The present results show that ΔLVEF in the in silico CRT simulation model was well correlated to that in the clinical observation. Therefore, absolute and relative ΔLVEF may be useful CRT response criteria in the in silico CRT simulation model.

The strength of this in silico CRT simulation model

Unlike the previous CRT simulation methods, our approach reflects the physiological measurement data in the model. Using this model, the CRT simulations for one patient were conducted, and the outcomes of the patient were analyzed. Furthermore, we demonstrated that this computational model can be used as a screening tool for CRT non-responders before CRT device implantation.

Study limitations

The individual anatomy of the Purkinje fiber networks cannot be visualized. The vascular resistance cannot be noninvasively measured. Therefore, the anatomy of the Purkinje fiber networks and vascular resistance were assumed in this in silico CRT simulation model based on the data from the general population. [11] This in silico CRT simulation model can reflect acute changes in LV volume and motion induced by CRT. Long-term reverse remodeling of the LV is not considered in this model. The present study was a pilot study and limited to only one patient. Because the present stud was retrospective, bias might be involved. Validating our model for a large number of patients is necessary for our model to be used as a screening tool for response to CRT.

Conclusion

This in silico CRT simulation method of a 3D patient-specific heart is feasible as a screen for CRT non-responders in patients with HF and LBBB.