ABSTRACT
Background Up to 30-50% of selected patients with chronic heart failure do not respond to cardiac resynchronization therapy (CRT). Optimization of pacing lead placement in ventricles remains a challenge.
Objective We utilize a machine learning (ML) classifier to predict the position of an optimal left ventricular (LV) pacing site maximizing the probability of CRT response for a certain patient.
Materials and Methods Retrospective data from 57 patients with implanted CRT devices were utilized. Positive response to CRT was defined by a 10% improvement in the LV ejection fraction in a year after implantation. For each patient, a personalized model of ventricular activation and ECG was developed based on MRI and CT images. The total ventricular activation time, QRS duration and electrical dyssynchrony indices during intrinsic rhythm and biventricular (BiV) pacing with clinical pacing lead position (ref-LP) were computed and used in combination with clinical data to train the ML algorithm. We built a logistic regression classifier predicting CRT response with a high ROC AUC=0.84 and an average accuracy of 0.77. It generates a ML-score estimating the probability of CRT response. ML-scores were computed from model-driven features for varying LV pacing sites. Then Bayesian optimization was used to interpolate the ML-score over the available LV surface and an optimal LV lead position for BiV pacing that maximizes ML-score (ML-LP) was defined.
Results The optimal LV pacing site position increased the average ML-score by 17% in the patient cohort. Moreover, 11 out of 34 (20%) non-responders classified as true negative at ref-LP were re-classified as positive at ML-LP. In a patient group (n=14, 25% of the cohort) with LV pacing lead deployed in close proximity to the optimal position, the ratio of responders to non-responders was three times higher than in the entire cohort.
Conclusion We have developed a new technique based on simulations and ML to define an optimal position for LV lead for BiV pacing maximizing the ML-score of CRT response on the available LV epicardial surface. This technique demonstrates a high potential for the improvement of CRT outcome with guided lead implantation.
1 INTRODUCTION
In addition to being an optimal medical treatment, cardiac resynchronization therapy (CRT) is an effective therapy for patients with chronic heart failure (CHF). According to the current 2021 patient selection guidelines (1), CRT is recommended for selected CHF patients in sinus rhythm with a reduced left ventricular ejection fraction (LV EF) <=35% and a long QRS duration (QRSd>=130 mc) reflecting left bundle branch block (LBBB). CRT delivers biventricular (BiV) pacing to correct electromechanical dyssynchrony in order to increase cardiac output.
Despite the well-documented CRT benefits for improving patient outcomes and reducing patient hospitalizations and mortality in the general population of CRT recipients, it still remains ineffective in 30-50% of cases depending on CRT response definition (2).
The problem of CHF patient stratification for CRT implantation and optimization of the procedure and device setting up is addressed in a number of clinical trials and simulation studies and still calls for the community attention (3, 4, 5). Pacing lead configuration for CRT was shown to be an essential determinant of patient improvement (4, 6). Recently, quadripole and multipole left ventricle (LV) leads have become available, which allow an operator to test different LV pacing configurations and choose the most appropriate setting (7, 8, 9). In addition to BiV pacing, His-Purkinje conduction system pacing is also considered as an alternative or additional possibility for optimizing pacing effects (10, 11).
Generally, there are two steps in CRT optimization that are being discussed. The first concerns the intra-operative guidance of LV lead implantation. In the standard procedure, a LV lead is introduced into the coronary sinus vein at an appropriate position by intra-operative X-ray imaging guidance based on empiric rules developed by practice but essentially dependent on the operator’s experience and skills. Much effort has been devoted to improving lead implantation. Several criteria of intra-procedural optimization have been proposed. Overall, clinical trials suggest avoiding the apical and anterior regions for lead positioning (4). In electrically guided implantation, the narrowing of QRS duration (QRSd) is used to predict the best electrical synchronization (12), but this approach has not shown credible evidence. New techniques of Body Surface ECG Mapping for non-invasive assessment of the ventricular activation are also being developed to improve electrical synchronization (13, 14). Novel cardiac imaging approaches demonstrate the benefits of pre-operational myocardial fibrosis and scar area assessment to avoid proximity to this area when implanting the LV lead (15, 16, 17). Frequently, the late activation time (LAT) area is considered as a target for LV lead implantation (18), assuming that pacing from this area provides more effective electrical synchronization.
In mechanically guided implantation, the area of late mechanical activation derived from modern echocardiography imaging data was suggested as a target for LV lead placement (19). However, although a recent systematic study and a meta analysis of the randomized clinical studies using this approach confirmed improvements in the NYHA HF class in the patients, no improvements were found in the LV remodeling and ejection fraction (LV EF) (20). Recent studies using intra-operative evaluation of hemodynamic features showed that an LV lead location maximizing both LV dP/dtmax and stroke work has potential for CRT response (21, 3, 8) improvement. Although great efforts have been made to improve LV lead implantation, none of the above approaches is included in the current guidelines for the procedure.
The other, more feasible step of LV pacing optimization is post-procedural choice of active pacing sites among the two, four or multiple sites that are available depending on the LV lead configuration (22, 12, 23). Many clinicians consider QRSd shortening or a QRS area reduction as a criterion for pacing pole selection or used echocardiography to control LV EF improvement during device programming. In today’s practice, CRT devices are conventionally programmed to pace from the LV active pole with the longest electrical delay (QLV delay). This parameter is measured in the sinus rhythm as the time interval from the ECG Q wave onset to the moment when the electrical wave arrives at the pole of the LV sensing lead (24). Two ongoing trials are also investigating whether a QLV targeted approach is beneficial compared to standard LV lead implantation (DANISH-CRT NCT03280862, ENHANCE CRT NCT01983293). In some studies, interventricular delay (RV-LV delay), as measured by the difference in activation time between the right ventricular (RV) and LV leads in sinus rhythm or at RV pacing is also used to optimize LV lead placement (25, 26). Nevertheless, although the results of numerous clinical trials have demonstrated various potential ways of lead placement optimization, no image-based approach has yet been recommended to guide the implantation strategy in routine practice. This emphasizes the need for further validation of useful imaging-based approaches and the development of new CRT optimization strategies.
Although it is generally accepted that the problem of CRT procedure optimization is essentially multifactorial, most of the approaches used to optimize ventricular lead placement have been based on some single features characterizing the resynchronization of the ventricular activation or contraction.
Modern machine learning approaches open up new perspectives on the analysis of a variety of patient-specific data that can be obtained using current cardiac imaging techniques of different modality. Recently, the use of unsupervised ML analysis enabled researchers to elicit four phenotypes in a patient population showing different probabilities of improvement due to CRT (27). Using supervised ML algorithms, a ML calculator based on a minimal set of conventional preoperative clinical data was developed to predict LV EF improvement with a high accuracy of about 0.7, which was higher than what other classifiers were able to attain (28).
Predictive models have been developed to estimate mortality or hospitalization risks from baseline clinical parameters (29, 30, 31), to assess improvements in LV EF based on baseline indices and analysis of medical records (32), and to stratify patients by an unsupervised learning approach implementing electrocardiography data (27, 33).
In current modeling studies, detailed cardiac anatomical models are also used to identify significant features related to CRT improvement (34, 35). In a recent article by Rodero et al. (36), an optimal pacing design with quadripolar LV lead in terms of minimizing total ventricular activation time (TAT) as a measure of the electrical dyssynchrony in ventricles was analyzed using personalized electrical ventricular models. The authors compared the effects of single and multi-site pacing with personalized and population-based optimal lead designs in case of infarct absence. Finally, they concluded that a single optimal lead design is sufficient to obtain near-optimal results across most patients. Moreover, modeling predicts a decreased effect of pacing on TAT reduction due to postprocedural ventricular reverse remodeling, indicating the need to re-tailor optimal lead design in postoperative follow-up. The above study shows the potential of virtual clinical trials as a tool for exploring new pacing lead placement designs.
In another recent study, Lee and co-authors (37) showed that simulations from personalized computational cardiac models and image based measured mechanical indices are predictive of an optimal LV lead placement for a positive acute hemodynamic response (37). Using machine-learning (ML) techniques, these authors showed that the combination of the RV-LV electrical delay and mechanical regional time to peak contraction can predict a positive response with an accuracy of about 70%.
In our recent proof-of-concept study, we utilized a combination of retrospective clinical data in patients who had undergone CRT, personalized cardiac modeling and an ML approach in order to develop a new technique for predicting the response of conventional BiV pacing prior to the CRT procedure (38).
We developed the ML classifier for CRT response based on a hybrid combination of clinical and model-derived data on ventricular geometry and electrical activation at both intrinsic LBBB pattern and BiV pacing (see Fig. 1 for a pipeline of ML-classifier development based on the hybrid data). A total of 7 selected model-driven features of myocardial activation and clinical biomarkers are fed to the ML classifier for calculating a ML score which allows one to estimate the probability of a more than 10% LV EF improvement and thus to predict potential CRT responders or non-responders. The ML classifiers on the hybrid dataset outperformed classifiers built upon the clinical data, showing a higher accuracy over 0.8 with sensitivity and specificity higher than 0.7. Moreover, the ML score showed a positive correlation with the percentage of LV EF improvement in our patient cohort, suggesting the possibility of a quantitative prediction of the CRT outcome.
Schematic outline of the ML model development. I. Building and calculation of a personalized electrophysiological ventricular model: 1. Processing of the CT imaging data. 2. Segmentation of the finite element meshes of the torso, lungs and ventricles; 2*. Personalization of the ventricular model: A. Rule-based generation of myocardial fibers. B. Assignment of the scar/fibrosis area in the ventricles (shown in back) and computing of the ventricular activation map at the baseline LBBB pattern and BiV pacing with clinical lead position. 3. Calculation of ECG signals from the ventricular activation map. II. Development of a supervised machine learning classifier: creation of a dataset contacting combination of the clinical data and simulated features from the electrophysiological model from each of the 57 patients labeled into responders and non-responders, supervised training of a ML classifier and calculation of ML-scores of CRT response.
As a hypothesis for the current work, we suggested that such ML classifier on hybrid data could be used to predict an optimal LV lead position for guiding lead implantation. The idea of the strategy is as follows (Fig. 2). A personalized ventricular model for a given CRT candidate is constructed using imaging data and is then used to calculate ML scores for various LV pacing sites located on the epicardial surface of different LV segments with the exception of carefully defined scarring regions. Then Gaussian process regression is applied to the ML score array for predicting the location of a pacing site at the LV surface that would maximize the ML score of CRT response for the patient. This location can be then used to target LV lead implantation.
Schematic outline of an algorithm for finding the optimal LV pacing site position on the epicardial LV surface of a personalized ventricular model. The pipeline includes three major steps: 1. Precalculation of the ventricular activation maps from the personalized model at BiV pacing from the centers of available segments of the AHA LV model with the exception of the septal and postinfarction scar segments (marked as dark gray at the AHA LV scheme on the left panel). 2. Iterative Bayesian Optimization procedure to interpolate ML-scores on the LV surface, which includes: computation of simulated features which in combination with patient clinical data are fed to the LR classifier for creation of an initial array of the ML-scores for interpolation on the LV surface. Gaussian Process (GP) regression model trained on the current ML-score array for estimating a GP acquisition function L(μ, σ) and predicting ML-score values on the entire available LV surface (see two color maps on the LV surface with shades of red for ML-score>0.5 and shades of blue for ML-score<0.5); finding a target point candidate, in which the maximum ML-score of the acquisition function is approached; calculation of a new ventricular activation map and simulated features at BiV pacing with the LV site located at the current candidate point. The simulated features in the next iteration step are fed again to the LR classifier to generate the ML-score which is then added to the ML-score array for retraining the GP regression for further interpolation of Ml-score on the LV surface. 3. If two iterations of the algorithm predict the same candidate point, the algorithm is considered to converge and the last point with maximal ML-score value provides an optimal LV pacing site. Resulting ML-score map is shown on the LV surface of the personalized LV model and on the LV AHA segment scheme.
The effects of ML-based optimal LV lead position on ventricular activation outlined in this study are compared with other strategies for lead placement based on the LAT area and TAT minimization.
2 METHODS
2.1 ML classifier of CRT response
In the present study, we used the same pipeline to develop a supervised classifier of CRT response as previously described (38) (Fig. 1). The ML classifier was trained and tested on a hybrid dataset consisting of clinical data from patients who had undergone CRT and simulated data from personalized computational models of cardiac electrophysiology.
2.1.1 Clinical data
The study involved the clinical data from 57 CHF patients. All patients were on optimal drug treatment following CRT device implantation at Almazov National Medical Research Centre between August 2016 and August 2019. The participants signed approved informed consent forms. The study protocol was approved by the Institutional Ethical Committee.
The criteria for patient inclusion in the study and the complete list of clinical data used to perform feature importance analysis for machine learning classifier development are presented in the Supplementary Materials (sec. Clinical data description).
Patients were evaluated before CRT device implantation and during the follow-up period of 12 months after implantation. The baseline clinical data of the patient cohort are presented in Table S1.
In addition to the standard protocol of patient evaluation for CRT implantation, we also acquired data from 12-lead ECG and echocardiography recordings prior and after device implantation. Moreover, we collected multichannel ECG with a special system (Amycard, EP Solutions SA). This tool helped to record a maximum of 224 unipolar ECG signals at various pacing configurations (intrinsic rhythm, RV or LV pacing, BiV pacing) during the follow-up in a year after CRT implantation. After recording multichannel ECG, we performed computer tomography (CT) to visualize the torso with applied electrodes and the heart. The series captured with a scanner (Somatom Definition 128, Siemens Healthcare, Germany) were imported into special Wave program version 2.14 (Amycard, EP Solutions SA) to reconstruct the 3-dimensional geometry of the torso and heart. Finally, epi/endo ventricle models were manually built with marked active poles of RV and LV leads deployed for BiV pacing simulations. Data from magnetic resonance imaging (MAGNETOM Trio A Tim 3 T, Siemens AG or INGENIA 1.5 T, Philips) with contrast (Gadovist or Magnevist) before CRT implantation were also used to detect scar/fibrosis area in the myocardium and to incorporate these data into a personalized ventricular myocardial model.
Responders and non-responders
Patient data were annotated into responder and non-responder groups according to LV EF improvement by more than 10%. Table S2 shows clinical data in the groups, indicating significant differences in the echocardiography indices.
2.1.2 Simulated data
Ventricular anatomy models
Based on the segmentation of CT imaging data, finite element models were constructed for the torso, lungs and RV-LV ventricles for each of the 57 patients (Fig. 1 I. 1.). A rule-based approach was used to simulate myocardial fibers architecture (39). MRI data on scarring and fibrosis areas in the myocardium were accounted for in the LV model using expert annotation of these areas within the 17-segment American Heart Association (AHA) model of LV (see (38) for more detail). The scar regions were simulated as non-conducting and non-excitable areas and the conductivity of fibrosis regions was decreased by 50%.
Myocardial electrical activation models
As in the previous work (38), we used an Eikonal model (40) to calculate electrical activation times at each point on the ventricular mesh. Cardiac tissue was simulated as an anisotropic medium with conductivities resulting in an excitation velocity ratio of 4:1 along vs across the myocardial fibers. The Eikonal model is currently widely used; it allows one to simulate the evolution of the cardiac excitation wavefront (41, 42, 43, 44, 45).
ECG calculation was performed using the Lead Field method proposed by Pezzuto et al. (42, 46). ECG signals were computed according to the standard 12-lead ECG definition and the lead-field approach allowed us to reduce calculation time more than 100x times.
Model parameter personification
Each patient-specific model assumed a uniform conductivity in the myocardial tissue across the entire ventricles. Then an optimization problem was solved for each of these models to estimate the global conductivity parameter, and minimize discrepancy between simulated and clinically measured QRS complex durations from 12-lead ECGs recorded in the patient. We used the L-BFGS-B algorithm to handle optimization in the model and the method proposed in (47) for automatic QRS onset and offset detection.
Pacing protocols
We simulated two pacing protocols – LBBB activation pattern and BiV pacing, and calculated model-derived features to be used for developing an ML classifier of CRT response.
To simulate the LBBB activation pattern, the RV endocardial surface was annotated and a Purkinje network was generated using the model proposed by Costabal et al. (48). The Purkinje system was isolated from the working myocardium and connected to it only at the ends of the Purkinje fibers. Activation started at the His node and spread throughout the conduction system with an excitation velocity of 3 mm/ms before approaching Purkinje-myocardial junction points. This activation map was then applied to initiate activation within the ventricular myocardium according to the Eikonal model.
Furthermore, we used the simulated LBBB ventricular activation map to define the area of LAT in every patient model. This area was then used as one of the optimal LV pacing site positions, and the effects of this pacing optimization were compared against other designs (see next section below).
For BiV pacing we used referent RV and LV pacing lead locations derived from CT images. A zero time delay was set between the RV and LV pacing sites as programmed in patients.
For both the LBBB and BiV protocols, the global conductivity parameters were fitted to the 12-lead ECG data recorded from the patients.
Simulated features used for developing ML-classifiers
Patient-specific models allowed us to identify several clinically important features affecting ventricular activation and the geometric properties of ventricles.
The first group of model-derived indices was defined from the ventricular anatomy models based on CT and MRI data, coupled with electrophysiological model simulations. We estimated the volume of postinfarction scar and non-ischemic fibrosis and their size relative to the myocardial tissue volume (MTV). Knowing the exact location of RV and LV pacing leads, we measured the time delay in the activation of the LV electrode later than the RV electrode (RV-LV delay) in LBBB. We also calculated the physiological distances between the RV and LV pacing sites (RV-LV distance), and the distances from the LV lead to the scar area (Scar-LV distance) and to the area of LAT (LAT-LV distance) under intrinsic rhythm, by solving an isotropic Eikonal equation, to define the distance from a certain point on the ventricular surface to a specific area. The latter distances mimic ones that can be directly measured from CT and MRI images by a ruler.
The second group of model-derived indices were calculated for LBBB and the BiV pattern of myocardial activation. We simulated activation maps and 12-lead ECG signals. Final calculations included the following biomarkers for both pacing modes: total ventricular 95% activation time (TAT95), computed as the time interval for 95% ventricular activation; maximum QRS complex duration; time delay between the total LV and RV activation time (ATRV LV); relative difference between the mean activation time of the LV free wall and septum (mATST LV = (LV latmean − STmean)/TAT, where LV latmean (ms) is the average activation time of the LV free wall, STmean is the average activation time of the septum, and TAT is the total ventricular activation time. Changes in indices at BiV pacing compared to LBBB in either absolute values or normalized to the LBBB values were also used for developing the ML classifier of CRT response.
A summary of simulated data is presented in Tables S1, S2 in Supplementary materials.
2.1.3 Machine learning model
We built ML classifiers based on the hybrid dataset containing the clinical and model-derived indices described above. The complete list of the clinical and simulated features fed to ML algorithms is shown in Figure S3. To train ML classifiers, the dataset was labeled as responders and responders according to the definition of response to CRT as an increase of more than 10% in LV EF. (38).
At the preprocessing step, non-categorical data were normalized by substracting the mean and dividing by standard deviation. Highly correlated features were also removed from the dataset by threshold > 0.85. The next step was to select a set of most significant features for CRT response prediction. Feature selection was done using a Leave-One-Out cross-validation technique: for each patient, we formed a training dataset from which we excluded that patient’s features and then predicted the ML Score value for that patient.
The logistic regression (LR) classification algorithm was trained on the hybrid dataset, and 7 most significant features with the highest LR weights were chosen as follows: 3 pre-operational clinical features: left ventricle ejection fraction before CRT (LV EF, %), body mass index (BMI, dimensionless), and end-diastolic diameter of LV (EDD, mm). Also, 4 model-derived features were taken into account: distance from LV pacing site to postinfarction scar area (Scar-LV distance, mm), total activation time 95% myocardium in the LBBB activation pattern (TAT95, ms), and inter-ventricular dyssynchrony indices in LBBB and at BiV pacing (ATRV LV,LBBB, ATRV LV,BiV, ms).
Finally, these 7 features were used to train the final LR classification model. The features and corresponding logistic regression weights are listed in Supplementary Table S3. The LR classifier generates an ML-score based on the clinical and simulated features of the patient, which gives an estimate of the probability that the patient would respond positively to CRT. These ML-scores then were combined into one set to build a receiver operation characteristic (ROC) curve, and the area under the ROC curve (AUC) was calculated (Fig. S1 in Supplementary Materials). A threshold value of the ML-score=0.5 was adopted to predict either positive (ML-score ≥ 0.5) or negative (ML-score<0.5) response to CRT in our patient cohort.
2.2 Optimization of LV pacing site position based on the ML-score
In step of the study, we used the final LR classifier to predict an optimal LV pacing lead position which maximizes the ML-score throughout the LV epicardial surface for each personalized model of the cohort (see Fig. 2 for the pipeline employed for searching an optimal LV lead position). Septal regions and scarring area were excluded from possible LV pacing site locations. First, we varied the position of the LV pacing site between the centers of LV AHA segments on the epicardial surface with the RV pacing site located at the clinical position derived from the CT scan.
For each LV lead position (up to 12 positions, 10 per model on average), we computed the electrophysiological model and extracted model-derived features from the simulations and fed them into the LR classifier to generate the ML-score. At the end of this step, an initial distribution map of the ML-score on the LV epicardial surface was generated. The small number of points with defined ML-scores did not allow us to accurately predict the optimal LV pacing site with maximum ML-score. Therefore, we used a Bayesian optimization method to interpolate the ML-score on the entire LV surface accessible for pacing. This method involves building a regression model and its iterative refinement before converging at the optimal solution.
Bayesian Optimization
Bayesian optimization is a derivative-free global optimization method that requires only a model evaluation function.
The iterative process of ML-score interpolation was performed using Gaussian process regression (GP regression) model (49). We used the current ML-score set (an initial pre-calculated ML-score vector in the first iteration step) to train GP regression and to predict the ML-score at every point on the LV epicardial surface (at every node of the mesh). Then we calculated the so called acquisition function: L(μ, σ) = μ + 2σ, where μ is an expected ML-score value predicted by GP regression and σ is a standard deviation of GP at this point (GP uncertainty value) (Fig. 2, step 2). After that, maximum L(μ, σ) was defined throughout the LV nodes, and this position was further used to calculate the electrophysiological model at BiV pacing with this LV pacing site and to compute the corresponding ML-score according to the LR classifier fed with new simulated features.
The Bayesian optimization method thus strikes a balance between finding points that allow one to refine the GP regression model (points with large uncertainty, i.e., large σ), and finding points where the value of the regression function is maximum (points with maximum μ).
In the next iteration step, GP regression was re-trained with the addition of the data from the new point on the LV surface and the algorithm was repeated. The optimal solution was considered to be found if the last two iterations of the Bayesian optimization algorithm predicted the same point. Finally, we obtained an LV epicardial surface map of ML-score values, predicting areas of LV pacing with either positive (ML-score ≥ 0.5) or negative (ML-score < 0.5) expectation of response to CRT and suggested the optimal position of LV pacing site maximizing the ML-score among all available LV surface positions. This map can guide LV lead placement during CRT implantation.
2.3 Optimization of LV pacing site based on LAT area or TAT reduction
In addition to the ML-score based optimization of pacing lead position in the personalized models, we also used LAT area defined in LBBB for LV pacing site location, as suggested in several clinical studies (50, 51). Another approach to LV pacing site optimization in our models was based on the minimization of TAT95, which is also considered as a target for LV lead positioning (52). The latter approach was implemented in our personalized models using an iterative procedure similar to the one we used for ML-score optimization. To this end, we generated an initial set of simulated TAT95 with BiV pacing under variation in the LV pacing site position between ventricular segments and further used Bayesian optimization of TAT95 over the available LV surface. As a result, we found an LV pacing site position with minimal TAT95 on the available LV surface in each personalized model of our cohort.
The effects of each of the approaches to LV pacing site optimization were compared with the effects at the clinical LV lead position and with each other.
2.4 Software
Cardiac electrophysiology was simulated using an in-house software based on the FENICS library (for solving PDE problems) (53) and VTK (for working with meshes). The scikit-learn library was employed for the machine learning: classifier development, statistical modelling, feature selection, cross validation, and ROC-AUC calculation, and the Pyro (54) library for GP regression and Bayesian optimization.
2.5 Statistics
Detailed analysis was performed using the IBM SPSS Statistics 23.0.0.0 software package (USA). For qualitative data, the frequency and percentage of total patients in the cohort were calculated. Quantitative data are presented as mean ± standard deviation. Comparisons between two dependent groups were made using Wilcoxon’s test for quantitative data and McNemar’s test for qualitative data. Comparisons between dependent groups (ref-LP and opt-LPs vs LBBB; opt-LPs vs ref-LP) were made using nonparametric Friedman’s two-way ANOVA, followed by a pairwise comparison adjusted for multiple comparisons. Comparison between two independent groups (non-responders vs responders) was carried out using Mann-Whitney test for quantitative data and Pearson’s chi-square test for qualitative data. Feature dependence was assessed using Spearman rank correlation test. The critical level of statistical significance was taken equal to 0.05.
3 RESULTS
3.1 Hybrid dataset of clinical data before and after CRT device implantation and model simulations in LBBB and at BiV pacing
Retrospective clinical and imaging data from 57 patients who had undergone CRT implantation were collected and analyzed. Each patient-specific RV and LV pacing lead positions as deployed during the implantation procedure were derived from the CT-scan and used further as a reference pacing lead position (ref-LP). RV electrodes were placed at a standard apical position in all the patients. Figure S2 shows the distribution of the referent LV pacing sites between the segments according to the 17 segment AHA LV model. It can be seen that in 50 (88 %) out of the 57 cases the LV lead was placed in the lateral wall, mostly in the mid- and basal segments. The LV lead was delivered to the inferior segment in only one case and to the anterior segments in 5 cases. In 2 participants, an apical LV lead position was observed.
The location and transmurality of myocardial postinfarction scar and fibrosis were defined from MRI image description (see Methods for detail). Figure S2 shows the distribution of the segments with scar and fibrosis between the 17 AHA LV segments in the patient cohort. In 10 (18%) out of the 57 patients the segments with implanted LV electrodes were concordant.
A population of personalized ventricular electrophysiology models was built for the patient cohort. The ref-LP was used to evaluate the effects of BiV pacing in the personalized models. A summary of the statistics for clinical data, CT/MRI derived data and model-driven biomarkers in the patient cohort is presented in Supplementary Materials S1.
The entire patient cohort showed an average positive response to BiV pacing, revealing itself in a decrease of -23 ± 14% in QRSd, a reduction in the end-diastolic and end-systolic volume (EDV and ESD, correspondingly) by -18 ± 31 and -24 ± 36%, and in LV EF improvement by 9 8% as compared with pre-implantation data. In consistency with the clinical data, average positive changes in the simulated features of myocardial activation in response to BiV pacing manifested themselves as a decrease in QRSd by -23 ± 13%,in TAT95 by -32 ± 17%, and in all computed indices of the inter- and intraventricular dyssynchrony by about 100% see Table 1.
Statistics of simulated total activation time 95% (TAT95) and QRSd at baseline LBBB activation and BiV pacing with different LV lead positions in the total patient cohort and groups of responders and non-responders.
The hybrid dataset combining clinical and model-driven data in our patient cohort was then classified into responders (n=23/40%) and non-responders (n=34/60%) to CRT according to the LV EF improvement of more or less than 10% (referred to below as the EF10 criterion). The two groups demonstrated significant distinctions in several indices associated with CRT response. Table S2 in Supplementary Materials compares the clinical and model-derived variables of the groups.
According to the classification criteria, the average LV EF improvement by 17 ± 5% in the responders group is well above 3 ± 5% for non-responders. The average reduction in EDV and ESV is about 5 times larger in responders. At the same time, in consistency with the clinical data, no difference in the relative decrease in both simulated TAT and QRSd at BiV pacing was found between the two sub-populations of models.
Regarding the CT/MRI derived geometry indices, we found that in 2 (9%) responders and in 8 (24%) non-responders the segments with implanted LV electrodes were concordant with scar area. Accordingly, a shorter distance from the LV pacing site to the scar/fibrosis zone in the non-responders group (26±24 mm in non-responders vs 39±22 mm in responders) was revealed, suggesting less effective pacing of the normal tissue in non-responders. No other indices showed any significant difference between the groups (see Table S2 in Supplementary Materials).
3.2 ML-classifier of CRT response built on a hybrid dataset
As we showed in our previous study (38), none of the individual biomarkers in the intrinsic LBBB activation pattern derived from either clinical or CT/MRI data allowed us to classify the responders and non-responders groups with a sufficient accuracy. Similarly, none of the simulated electrophysiological biomarkers enabled us to distinguish between the groups both in the LBBB mode of activation and at BiV pacing (see Table S2 in Supplementary Materials) and, thus, consider them as individual classifying features.
This prompted us to use a combination of the clinical data recorded prior to the operation and MRI/CT derived biomarkers together with simulated features from personalized models of ventricular excitation in the LBBB and BiV pacing activation modes to build an ML-classifier, which improved significantly the accuracy of CRT response predictions in our patient cohort. The features fed into the feature selection algorithms when developing CRT response classifiers are listed in Figure S3 in the Supplementary Materials in the descending order of feature importance.
We used the hybrid dataset containing 57 data entries for every patient from our cohort. The supervised classifiers were trained using an EF10 criterion (ΔEF>10%) of CRT response. In this article, we developed a Logistic Regression (LR) classifier using Leave-One-Out and five-fold cross-validation and 3 different feature selection methods inside the cross-validation loop to train the classifier. Seven most important features were selected for developing the final LR classifier used further below in this article (see Table S3 for the LR independent variables and coefficients). The feature sub-set contains 3 pre-operational clinical features: LV EF, %, body mass index (BMI, dimensionless), end-diastolic LV diameter (EDD, mm); and 4 model-driven features: distance from LV pacing site to scar area (Scar-LV distance, mm), total activation time of 95% of myocardium in the LBBB activation pattern (TAT95LBBB, ms), and inter-ventricular dyssynchrony indices in LBBB and at BiV pacing (ATRV LV,LBBB, ATRV LV,BiV, ms).
The best supervised LR classifier we developed features a high ROC AUC of 0.84 with a total accuracy of 77%, sensitivity of 65% and specificity of 85% (see the complete list of the classifier characteristics in Table S3 in Supplementary Materials). Based on the hybrid clinical and model data, the LR classifier calculates an ML-score of CRT response as an estimate of the probability of a higher than 10% LV EF improvement for a patient. The threshold ML-score classifying patient data into the responders or non-responders group was found to be 0.5. This value will be used hereafter to classify test data in various LV pacing configurations.
The ML-score generated by the LR classifier allowed as to distinguish between the groups (0.64 ± 0.30 vs 0.25 ± 0.22, p<0.01, see Table S2). Moreover, the ML-score correlates with post-operative improvement in the EF (r=0.57, p<0.001).
These findings convinced us to use the ML-score from the LR classifier to evaluate the effects of pacing lead position on CRT response prediction and to suggest an ML-based strategy for lead placement optimization.
3.3 Optimal pacing lead position based on the model simulations and ML prediction of CRT response
In the previous section, we trained our LR classifier using post-operative CT data on precise clinical location of RV and LV lead pacing sites. However, when we use the classifier to compute an ML-score for selecting CRT candidate prior to the procedure, we do not know the exact lead location. In this section, we describe a technique that uses ML classifier to define an optimal LV lead position (opt-LP) maximizing the ML-score of CRT response for a given patient.
A general scheme of the technique for pacing site optimization based on ML-score is shown in Figure 2. A personalized ventricular model for a given CRT candidate is constructed using imaging data and employed to calculate model-driven features in the LBBB activation pattern and at BiV pacing with various LV pacing sites located on the epicardial surface in the centers of LV segments. Septal segments are excluded from LV stimulation area tested as unavailable for the transvenous approach. Neither are scarring regions used for stimulation as non-excitable. The selected clinical data and model-driven features for every tested pacing configuration are then fed into the LR classifier we developed to compute ML-scores. At the next step, Gaussian regression is applied to the ML-score array on the LV surface for predicting the location of the pacing site that maximizes the ML-score of CRT response for the patient. If the best possible ML-score corresponding to the opt-LP is higher than the threshold value of 0.5, this patient is classified as a responder. In this case, the opt-LP could be used to guide targeted LV lead implantation.
Hereinafter in the text, we refer to the BiV pacing configuration using the CT-derived spots of RV and LV lead tips as the reference lead positions (ref-LP) for each patient-specific model. The simulated characteristics of ventricular electrical activity for different opt-LP criteria are compared with model results for the original LBBB activation pattern, BiV pacing from ref-LP, and between each other.
3.3.1 ML-score based optimal LV pacing site versus reference lead position
Figure 3 shows three examples of optimal LV pacing sites in personalized ventricular models. Two-color maps of the ML-score value are shown on the LV epicardial surface of the personalized models and on the LV AHA segment schemes. Red shades show ML-scores>0.5 (see the color scale in the Figure) in desirable LV segments with pacing sites predicting a positive response to CRT. In contrast, shades of blue show ML-scores<0.5 in LV segments unwanted for LV pacing. Blue and red dots on the map show the locations of the clinical and optimal LV pacing sites, respectively.
Examples of optimal LV pacing sites in personalized ventricular models. Two color maps of the ML-value are shown on the LV surface of personalized models and on the LV AHA segment schemes. Dark gray at the AHA LV scheme marks segments containing postinfarction scar, which are excluded from pacing. Shades of red show ML-scores>0.5 and shades of blue show ML-scores<0.5. Blue and red dots show locations of the clinical and optimal LV pacing sites. From left to right are shown examples of the ML-score map in the clinical responder (patient #2), non-responder (patient #1), and non-responder (true negative at the ref-LP) predicted as positive to CRT response at the optimal ML-based lead position (patient #7).
The left panel in Figure 3 demonstrates the ML-score map in a clinical responder (patient #2) with LV EF improvement of 12%, which is higher than 10%. Here, the referent and optimal pacing sites are located in adjacent LV segments and the maximum ML-score of 0.95 at ML-LP is slightly above the referent value of 0.94. So, this patient is predicted as a true positive for CRT response (ML-score>0.5) at both the ref-LP and ML-LP. It is interesting that the patient is predicted as positive to CRT with any available LV site located at the lateral wall. The LR prediction is in line with the great extent of LV EF improvement in this patient.
The center panel in Figure 3 shows the ML-score map for a clinical non-responder (patient #1). Here, the ML-scores at both the referent and optimal LV pacing sites are blue colored (0.14 and 0.27 < 0.5, respectively). Moreover, the overall map of ML-scores on the entire available LV surface is blue colored, predicting a low possible response to CRT in this patient. Correspondingly, this patient has a large postinfarction scar spreading over half of the LV segments, and LV EF improvement in this patient is 7%, which is less than 10%.
The right panel in Figure 3 shows the ML-score map for a clinical non-responder (patient #7). The patient with an LV EF of 6% was classified by the LR predictive model as a true negative at the clinical LV pacing lead location (ML-score=0.38, see the referent pacing lead located in the blue color area on the ML-score map). At the same time, our algorithm predicts a narrow red area at the basal inferior segments where the patient is predicted as a positive for CRT response, particularly with the optimal ML-based LV lead position (ML-score=0.77, see optimal pacing site located in the red color area on the ML-score map). Thus, our simulations suggest that this patient could possibly improve with the ML based optimal pacing lead placement (ML-LP).
Figure 4 compares the distributions of the ML-scores for the reference and optimized LV pacing sites. The average ML-score for the ML-score based opt-LP (ML-LP) is higher as compared to that for ref-LP (0.58 ± 0.30 vs 0.41± 0.31, p<0.01). In particular, the ML-score is increased in 88% of the patients (51 out of 57, 19 responders and 32 non-responders). Importantly, our model predicts a much higher increase (almost double) in ML-score for ML-LP in the non-responders group (see Table 2).
Dependence of CRT response characteristics on the LV pacing lead position. Top panel: ML-score at baseline LBBB activation and at different LV pacing lead positions. Bottom panel: TAT at LBBB and different LV pacing sites. Comparison of dependent groups was performed using the Friedman’s test, followed by a pairwise comparison adjusted for multiple comparison.
Statistics of the ML-score and distance from LV pacing site to the scar area (Scar-LV distance) at different pacing Lv lead positions in the total patient cohort and groups of responders and non-responders
In a majority of models (52 (91%) out of the 57 models), BiV pacing from ML-LP reduces considerably all simulated features characterizing ventricular activation as compared with the LBBB models. The average TAT (104 ± 26 ms) and QRSd (146 ± 27 ms) at ML-LP are significantly shorter (p<0.01) than those at LBBB activation (150 ± 26 and 189 ± 24, respectively). At the same time, no significant difference in average TAT and QRSd was found between ML-LP and ref-LP (Table 1).
Both the inter- and intra-ventricular dyssynchrony indices in the LBBB activation pattern have positive average values reflecting significantly later activation of LV versus RV, and LV lateral wall versus septum (Table S5). Both indices reduce several times at ML-LP as compared with LBBB. However, no difference between ML-LP and ref-LP was found in the inter-ventricular dyssynchrony index ATRV LV, showing a delay in total activation between LV and RV. In contrast, the intra-ventricular dyssynchrony index mATST LV, representing the difference in total activation time between the lateral LV wall and septal part of LV, is slightly higher for ML-LP. The average positive mATST LV (0.09 ± 0.06) at ML-LP suggests later activation of the LV lateral wall as compared to the septum, while the negative index (−0.07 ± 0.10) at ref-LP reflects later activation of the septum.
3.3.2 Optimal ML-score based LV lead position versus other optimized LV lead placement
In this subsection, we compare the ML-LP results with other opt-LPs based on different criteria used to optimize lead placement in clinical studies.
First, we defined the latest electrical activation time (LAT) area on the LV subepicardial surface in the LBBB activation pattern. Then, we used these LAT spots for LV pacing (referred to hereafter as LAT lead position, LAT-LP) in patient-specific models, calculated the ventricular electrical activation features and fed them together with clinical data into the LR classifier to predict ML-scores of CRT response under LAT-LP pacing for the patients.
In addition, we defined a model feature based opt-LP minimizing a certain single model-derived feature characterizing ventricular activation. We tested opt-LP based on either minimal QRSd or TAT, or LV electrical dyssynchrony indices derived from model simulations. Here, we again used a similar Gaussian regression approach for predicting the optimal feature value on the LV surface as employed for ML-score optimization. Note that the same reference RV pacing site was used for each of the tested opt-LP, unless otherwise specified. We then compared the effects of optimized LV pacing positions between ML-LP, LAT-LP and minimum TAT based opt-LP (TAT-LP).
A summary of the model-derived features characterizing ventricular activation and ML-scores generated by LR classifier of CRT response using model-derived features of different opt-LP is presented in Tables 1, 2, S5 and S6. First of all, we compared ML-scores predicting the probability of LV EF improvement with different pacing lead configurations based on the model simulations.
Figure 4 compares the distributions of ML-scores generated by the LR classifier of CRT response using the features simulated for different opt-LPs at BiV pacing. The only maximum ML-score based optimization of LV pacing site predicted an increase in the average ML-score in our patient group as compared with ref-LP and other tested lead configurations (Table 2). Unexpectedly, no difference in the average ML-score was found between optimized LAT-LP, TAT-LP and ref-LP in the entire population. Each optimized LV lead position predicted a higher ML-score in the responders group versus non-responders. However, only ML-LP pacing caused an increase in ML-scores in each group, showing a much higher increase in the non-responder group (Table 2).
As defined above, the ML-score=0.5 is a threshold separating potential responders and non-responders according to our LR classifier trained on the ref-LP data. For each opt-LP, we found cases where our LR predictive model classified a patient into the opposite group as compared to ref-LP. In other words, we found patients who were classified as negative (potential non-responders) with ref-LP but re-classified as positive (potential responders) with opt-LP and vice versa. Figure 5 shows such transitions from the group of ML-score<0.5 to the group of ML-score>0.5 and back due to opt-LP pacing in the model. It can be seen that in the case of ML-LP pacing (left panel), there are only upwards transitions from the underside group of potential non-responders with an ML-score <0.5 for ref-LP to upside group of potential responders with an ML-score >0.5. There are 11 such transitions, which are shown in Table S7 in more detail. Here, 5 responders classified by the LR classifier as false negative in ref-LP move upward into the positive group in ML-LP (see +5 in the top left cell coming up from the bottom left cell). Moreover, 6 non-responders truly classified as negative in ref-LP are re-classified as positive in ML-LP (see +6 in the top right cell coming up from the bottom right cell). In total, according to the LR classifier the ratio of positive to negative CRT responses with optimized ML-based pacing lead placement rose considerably to 31-to-26 (54-to-46%) versus 20-to-37 (35-to-65%) in ref-LP.
Transitions from negative to positive prediction of CRT response and reverse transitions when switching from ref-LP to opt-LP. Top: ML-LP; Middle: LAT-LP; Bottom: TAT-LP.
In contrast to ML-LP pacing, LAT-LP and TAT-LP pacing did not show such promising predictions (Fig. 5, Table S7). There are far fewer positive transitions from negative to positive predictions (4 for LAT-LP and 3 for TAT-LP) with these opt-LP. Moreover, there is a number of reverse transitions from positive to negative prediction in opt-LP (3 for LAT-LP and 6 for TAT-LP), making one to anticipate an unlikely CRT response as compared to ref-LP. In total, the ratio of positive to negative predictions in our population for LAT-LP and TAT-LP is not improved as compared with that in ref-LP pacing (Table S7).
3.3.3 Effects of different opt-LPs pacing on simulated features
To explain the difference in the distributions of ML-scores depending on pacing lead configuration, we compared the model-derived features characterizing ventricular activation in different opt-LPs. While BiV pacing at every optimized LV lead position reduces ventricular dyssynchrony as compared to the LBBB activation pattern, only opt-LP minimizing LV TAT causes a significant reduction in the average TAT and QRSd in the entire population as compared to ref-LP, ML-LP, and LAT-LP (Table 1, Fig. 4). No difference in both TAT and QRSd was found between the latter three opt-LPs (Table 1). No correlation between TAT and LV EF improvement (r=-0.20, p=0.134) was found for ref-LP pacing. So, the shortest TAT could unlikely predict with confidence the best response in the patients.
Similarly, there were no effects of pacing configuration on TAT and QRSd in the responders and non-responders groups except for TAT-LP, which shortened these features in each group. No significant difference between the groups was found independently of the BiV pacing configuration.
The inter-ventricular dyssynchrony index ATRV LV reduced under BiV pacing by about 100% independently of LV pacing site optimization approach in the entire population, and no difference in change in the index was observed between the responder and non-responder groups (Table S5). Similarly, the intra-ventricular electrical dyssynchrony index mATST LV significantly reduced under BiV pacing as compared to the LBBB activation pattern, and there was not much difference in the index between the opt-LPs and between responders and non-responders.
Thus, the peculiarities of these model features under different opt-LPs could not explain the differences in the ML-score we found.
As described above, the distance from the LV pacing site to postinfarction scar area (Scar-LV distance) was selected as one of the significant model-driven features affecting LR classifier accuracy, and ML-score value. The average Scar-LV distance is significantly higher for ML-LP as compared to ref-LP, LAT-LP and TAT-LP, contributing to the higher ML-scores for this stimulation pattern. This conclusion is supported by the correlation (r=0.673, p=0.000) found between the improvement in the ML-score and the extension of the Scar-LV distance when switching the pacing lead configuration from ref-LP to ML-LP.
Comparing the responder and non-responder groups, we observe that the distance from LV lead to scar is shorter for the non-responder group for ref-LP, while in ML-LP the distance is significantly increased in the non-responder group, blurring the difference between the groups (Table 2). This contributes to a stronger improvement of the ML-score in the non-responder group under ML-LP pacing. Finally, no correlation between maximal ML-score and the Scar-LV distance was found for ML-LP. In contrast to ML-LP, the Scar-LV distance reduced for LAT-LP as compared to ref-LP (Table 2). This could have been expected, since the electrical propagation should slow down closer to the scarring areas. Hence, the LAT area could be close to the scar. Unexpectedly, the Scar-LV distance reduced with TAT-LP as well, compared to ref-LP (Table 2).
For both LAT-LP and TAT-LP, the reduction in Scar-LV distance is larger in the responders group so that the overall effect of such opt-LP on Scar-LV distance in the entire group is not significant.
4 DISCUSSION
4.1 ML-score based optimal LV lead position
In this study, we have developed an LR classifier of CRT response predicting a positive long-term LV EF improvement of more than 10%. This classifier is based on pre-operative clinical data in combination with simulated features from personalized ventricular anatomy and electrophysiology model computing ventricular activation in the LBBB pattern and under BiV pacing. The classifier was trained and tested on the data obtained from BiV pacing with reference RV and LV lead locations as that delivered in the clinic. The precise RV/LV pacing lead sites were determined from the post-operative CT scans and used in personalized ventricular models for BiV pacing. This was performed to exclude the effects of uncertainty in lead position on the ML prediction results. The classifier showed a great accuracy of more than 0.8 with a sensitivity and specificity over 0.7. This was the first step in the validation of our new technique suggesting its high potential for CRT response prediction.
In actual practice, such an ML classifier should be used before the procedure when the lead position is yet uncertain. Moreover, a good prognostic model should help to first decide whether a patient should be selected for CRT and, on top of that, help with procedure planning once the patient has been selected. The main advantage of using personalized computational models is the possibility of computing the characteristics of ventricular activation from any accessible pacing sites and to predict an optimal lead placement against certain optimization criteria prior to the procedure. ML classifier allows one to estimate a probability of CRT response for various opt-LPs and help with decision making.
In this article, we propose a novel approach to using a ML classifier directly for choosing LV lead position on the LV surface. The patient assessment algorithm involves building a personalized ventricular model, performing simulations in the LBBB activation mode and at multiple LV pacing site positions for BiV pacing. It then involves calculating the ML-score using the LR classifier based on clinical data and model-derived characteristics for each BiV pacing configuration, and finally, applying Gaussian regression with respect to the ML-score throughout the entire LV surface (Fig. 2). Thus, this approach provides a CRT operator with a surface ML-score map that predicts areas of LV pacing lead location with either positive (ML-score>0.5) or negative (ML-score<0.5) predictions of LV EF improvement. Moreover, this approach suggests a target position for optimal LV lead placement that maximizes the ML-score, predicting the highest possible probability of response to CRT (Fig. 2). If the optimal ML-score is high enough (higher than 0.5 in the case of our LR classifier), the patient could be considered a candidate for CRT. The ML-based optimal LV lead position (ML-LP) could serve as a target for lead implantation (see examples of ML-score maps in patient specific models in Figure 3).
Existing lead optimization strategies are predominantly based on a single pre- or intra-operative feature (either e.g. LAT, or QRSd, or TAT, or LAT, or dP/dtmax, etc.). In contrast to these strategies, our ML-based lead optimization accounts for several significant features related to the CRT response. This approach includes realistic model-driven data on the size and location of the myocardial damage area, distance between the LV pacing site and the scarring area, and simulated features of ventricular activation at BiV pacing from certain LV pacing sites.
In fact, we did know the response to CRT in our patients with just one reference pacing lead configuration (ref-LP). Therefore, we first compared the results of optimal BiV pacing design with those of ref-LP, and the ML-score predicted by our LR classifier was used as an estimate of the probability of CRT response depending on lead configuration.
The most essential result of our study is that the maximal ML-score based opt-LP (ML-LP) provides the highest ML-score in our patient population among the pacing configurations we tested (Fig. 4, Table 2). This result was to be expected but was not obvious. Indeed, we were sure that variations in LV pacing site location throughout the available LV surface (with the exception of the septal segments and scarring area) would affect the ML-score essentially. The ML-score varied more than 10-fold across the LV surface, and the range was much broader in the non-responders group of patients. This allowed our algorithm to choose a maximal ML-score value for each patient, which exceeded the reference value by 17±14% in 89% of our patients. Moreover, in the non-responders group the maximal ML-score increased by 19% in 32 out of 34 patients. This high potential of ML-based optimization of pacing lead placement is clearly visualized in Fig. 5. It shows a great number of transitions among the patients classified at ref-LP into the group of negative expectation for CRT response (ML-score < 0.5) to the group of positive expectation (ML-score > 0.5) at ML-LP. It is especially promising that 6 true non-responders in ref-LP are predicted as positive for CRT response in ML-LP pacing.
4.2 Validation of the optimal ML-LP approach
Working with retrospective clinical data in the present study, we had no possibility to verify our predictions of optimal LV lead position with maximum ML-score in a prospective group of patients. That is why we chose to start working with clinical data having the referent LV pacing lead (ref-LP) deployed during the procedure close to the optimal position of the pacing site (ML-LP) predicted from our simulations. The proximity of the referent pacing site to an optimal position was defined as 25% of the shortest distance from the distribution of the distances in our population of 57 patients (lower quartile).
There are 14 models in this suboptimal group (with a suboptimal referent LV lead position). The distance from the referent to the optimal LV pacing site position varies from 0 to 26 mm. In 8 models the optimal and referent LV sites are located in the same LV segment of the 17-segment AHA model, while in the other 6 models, - in neighboring segments. There are 9 responders and 5 non-responders in the suboptimal group. Our ML-classifier predicts 11 positive responses with an ML-score>0.5 for ML-LP and 3 negative responses with an ML-score<0.5. From the 11 individuals predicted as positive in ML-LP, 9 cases contain all true-positive responders, and 2 cases are false-negative non-responders. Considering data for the 2 false-negative individuals more accurately, we found that patient #25 was predicted as positive for both ref-LP and ML-LP with rather a high ML-score of about 0.8. An optimal LV pacing site was predicted at the basal-anterior segment #7, while the referent location was in the neighboring mid-anterior segment #1. The patient was initially labeled as a non-responder according to the EF10 criterion because his LV EF showed an improvement of less than 10% at 2% only. However, his AHA HF functional class improved from 3 to 2, and other echocardiography indices used to evaluate CRT response indicated significant LV reverse remodeling in this patient, particularly a reduction both in EDV by -22% and in ESV by -25%. Thus, this patient would have been surely labeled as responder if other response criteria had been used. So, in this case the ML-score based positive prediction should be considered as true-positive in terms of CRT response. The second patient #45 was also labeled as a non-responder as his LV EF improvement was 8%, which is less than 10% but close to this threshold. At the same time, this patient demonstrated also an improvement in the functional class from 3 to 2, and an essential reduction in EDV by -24% while ESV reduction of -3% was not as much pronounced. Based on the clinical data and simulations with ref-LP, the patient was classified as true negative with an ML-score of 0.39 (< 0.5). However, the maximum ML-score of 0.55 generated from the simulations at the optimal pacing site position predicts this patient as positive, although the ML-score is also close to the threshold. The optimal LV pacing site was suggested to be in the basal-anterolateral segment #6 which is the same as the referent position. Nevertheless, the optimal LV lead position was predicted to be 34 mm away from the scar, while the suboptimal referent position was 17 mm closer to the scarring area. This close proximity would possibly not allow the patient to improve as much as could be expected from the optimal ML-LP.
If we consider the 3 non-responders from the suboptimal group, each is classified as true-negative with both ref-LP and ML-LP (ML-score < 0.5), and the ML-score at the optimal lead position is not much different from the reference, suggesting a small chance of improvement. CRT should not possibly be recommended for these patients, as it would unlikely be effective.
Thus, in the patient group with suboptimal reference LV lead position the ratio of patients with positive to negative prediction of CRT response is 11 (79%) to 3 (21%), which is much higher than that in the entire group of patients. The accuracy of the LR classifier in the group is also higher (85%). In addition, the ML-score based predictions of CRT response with optimal LV lead position are supported by the difference in the range of LV EF improvements in the positive and negative patients with a median of 14[11.20] for responders versus 2[1.10] for non-responders within the group with suboptimal BiV pacing.
By contrast, in the 25% group of patients with a long distance between the reference and optimal pacing sites (higher quartile of the distance distribution in our population, ranging from 79 to 117 mm) the ratio of 3 responders to 11 non-responders is much lower than in the entire group. This also indicates that pacing lead placement away from the optimal position is associated with a highly likely negative prognosis of CRT response. Note that in a majority of non-responders from this group the pacing lead is placed close to the scarring area, which could fortunately be avoided in the case of optimal lead positioning. The general trend is that the closer the referent pacing site is to the scarring area, the farther away is the optimal pacing site from the reference position and from the scar. Moreover, even in this far-from-optimal group, 3 non-responders were classified as positive at the optimal pacing site position suggesting a chance for a better outcome against what empiric lead implantation provided.
The above results evidence a high potential of our ML-score based optimal LV lead placement for stratifying CRT candidates and guiding lead implantation.
4.3 Which LV lead position is better?
To check the ML-score based optimization approach for other advantages, we compared ML-LP with the other two opt-LPs based on model-derived feature optimization. In the first approach the personalized models were computed for the LBBB activation pattern, and LAT area was derived from the simulated activation map on the LV surface. The LAT sites were then used as LAT-LP to pace LV in the models. This approach is similar to that automatically implemented in CRT device programming, particularly for quadripolar LV electrodes. The maximal Q-LV delay used for choosing the LV active pole reflects the latest activation time among the electrode tips/rings suggesting the most effective re-synchronization of the ventricular activation. We implemented our personalized models to calculate the interventricular RV-LV delay as a time interval between activation of the simulated LV and RV pacing sites in the LBBB pattern. This feature is often used instead of Q-LV, especially in case of RV pacing.
Comparing the different lead optimization designs, the RV-LV delay was found to be expectedly longer for LAT-LP (147 ± 42 ms) as compared to ref-LP and optimized ML-LP, TAT-LP (103 ± 40, 96 ± 47, 80 ± 37, p<0.01, respectively). Meanwhile, no pairwise difference in the average RV-LV delay was found between the latter three pacing configurations. Moreover, we found no significant difference in the average RV-LV delay at ref-LP between the responder and non-responder groups in our patient population (see Table S6). Neither was there any correlation between the RV-LV delay and LV EF improvement (r=-0.14, p=0.314). The lack of correlation was also confirmed between the geometry distance from RV to LV pacing site and EF improvement. Furthermore, no statistically significant correlation was observed between the RV-LV delay and corresponding ML-score predicting the probability of EF improvement in a certain pacing lead configuration (r=0.02, p=0.904).
There is no consensus on the role of RV-LV delay for optimal lead configuration. There are clinical studies showing RV-LV delay to be predictive of CRT response (18). In a recent simulation study (37) using personalized ventricular models, the authors showed the RV-LV delay as a significant simulated feature the value of which above 60% predicts solely an acute hemodynamic response to BiV pacing. In our investigation, the RV-LV delay did not demonstrate the ability to classify the patient cohort into responders and non-responders with applicable accuracy (ROCAUC=0.623, p=0.118). No significant difference in the LV EF improvement in the groups with an RV-LV delay of more or less than 60% (11 ± 10% versus 8 ± 7%, p=0.085) was revealed in our cohort. The cause of the inconsistency could be an acute response to BiV pacing used to train the LR classifier based on the RV-LV delay in (37). By contrast, we used clinical data on long-term EF improvement in a year after the procedure where the impact of the RV-LV delay could have been less important. This hypothesis is supported by our recent data from clinical observations (55). We compared two groups of patients with quadripolar LV leads: one group featured an optimal choice of LV pacing pole with maximal RV-LV delay while in the other group the maximal RV-LV parameter could not be set for BiV pacing. We showed a faster improvement in the first group during the first 3-6 months after implantation. Our finding in that study coincided with the simulation data from Lee et al. (37)). On the other hand, in that study we observed no difference in the LV EF improvement or ESV reduction in 12 months after operation between the groups. That finding is in agreement with our model predictions. In addition, our simulation data are in line with the results of a recent ENHANCE-CRT study that showed no discernible difference between Q-LV delay optimized localization in comparison to conventional anatomical lead location in the long-term response of non-LBBB patients (56).
Another tested opt-LP approach is minimum TAT based optimization of LV lead position (TAT-LP), which is frequently addressed in clinical and in-silico studies. TAT is often used as a measure of ventricular dyssynchrony and its reduction via BiV pacing or by other pacing settings is considered as a target for stimulation design (36, 52). Many clinical trials report the effects of pacing on TAT as an estimate of its quality. In clinical practice, direct assessment of both LAT and TAT is complicated and requires invasive electrophysiological mapping to be performed. Currently, noninvasive body surface ECG mapping is also used to solve a reverse problem of electrophysiology to assess the activation pattern in the ventricles (57, 58, 59). However, it still requires manipulations on patients and can be performed for CRT patients only after device implantation. By contrast, personalized cardiac models present another useful tool for noninvasive prediction of LAT in LBBB and TAT for various pacing configurations prior to the procedure.
It is noteworthy that both simulated TAT and QRSd showed a significant reduction with each BiV pacing setting as compared to an inherent LBBB activation pattern (Table 1, Fig. 4). Comparing the three optimal pacing settings tested, we found essentially shorter mean values of TAT and QRSd for TAT-LP minimizing TAT across the LV pacing site on the surface of the LV. At the same time, neither TAT nor QRSd showed any significant pairwise intergroup difference between ML-LP, LAT-LP, and ref-LP.
Analyzing our clinical data, we found no correlation between LV EF improvement and either absolute QRSd at BiV pacing or its change relative to the value at sinus rhythm (−0.14<r<0, p>0.318). Consequently, no correlation was found between simulated absolute values of TAT or QRSd at BiV pacing or their relative change against the LBBB pattern (ΔTAT, ΔQRSd) and ML-score (0<r<0.18, p>0.321, not shown). Similar results were obtained for inter- and intra-ventricular dyssynchrony indices, which highly reduced by about 100% at every BiV pacing setting against LBBB activation mode, but did not change to any considerable degree in the optimal pacing configuration as compared to the ref-LP (Table S5).
Our simulation results are consistent with the simulation results published in a recent article (37). Lee et al. also did not find that changes either in QRSd or in the bulk (10–90%) ventricular activation time (risetime), or in the time of LV activation at BiV pacing in comparison to RV pacing were predictive of a more than 10% acute hemodynamic response to BiV pacing.
Much more unexpectedly, our ML-LP pacing predicts higher ML-scores as compared to the other two optimal pacing lead configurations LAT-LP and TAT-LP. Moreover, the average ML-scores for LAT-LP and TAT-LP do not differ from that for ref-LP and between each other (Fig. 4, Table 2). A similar tendency can be seen in the responders and non-responders groups. Surprisingly, in several models LAT-LP and TAT-LP pacing even caused a decrease in the ML-score and in transition from the group of positive to negative CRT response prediction (Fig. 5).
Thus, neither the longest RV-LV delay for LAT-LP nor the narrowest TAT and QRSd for TAT-LP account for an essential difference in the ML-scores versus ref-LP and between each other. In contrast, ML-LP provides higher ML-scores despite higher mean TAT and QRSd values compared to TAT-LP or shorter RV-LV delay compared to LAT-LP. This result suggests that a uniparametric strategy (based on either TAT or LAT) for targeting LV lead placement cannot guarantee the best possible effect of pacing in terms of CRT response prognosis.
4.4 On the regional distribution of optimal LV lead position
Several clinical trials recommend avoiding apical and anterior regions for LV pacing, where possible (4). Figure S2 compares the distribution of segments with LV pacing sites depending on LV optimization approach. It can bee seen that in ref-LP the lateral segments with LV pacing sites are more frequent (50 out of the 57 cases) in our population.
In LAT-LP, the segment distribution is very similar to the 52 cases of LV pacing lead located in the lateral wall, but inferior lateral segments are more frequent than anterior lateral ones as compared to the ref-LP pacing sites. On the contrary, in TAT-LP pacing, lateral segments are less frequent (16 cases), while in the majority of 21 cases (16 non-responders and 5 responders), LV pacing sites are located in anterior segments against 5 cases (3 non-responders and 2 responders) for ref-LP. Thus, the prognosis of the number of potentially negative patients (all of them are true non-responders in ref-LP) paced from anterior segments is significantly increased at minimum TAT optimization as compared to ref-LP, suggesting worse CRT response expectations in the entire population.
In ML-based optimization, lateral segments are more frequent (36 cases) as well as in ref-LP, but both anterior (11 cases) and inferior (10 cases) segments are also representative in terms of the maximum ML-score. Analyzing the 11 cases with chosen optimal anterior segments for ML-LP, we found 4 responderss which were predicted as true positive (ML-score > 0.5) at ML-LP, and their referent lead position was close to the optimal one (in the same or neighboring segments with a small distance from the optimal to referent site). The rest 7 non-responders were still predicted as negative at ML-LP (ML-score > 0.5), despite an increase in their ML-score at the optimal pacing site as compared to the referent lead position. Thus, in these cases even an optimal lead configuration with the highest possible ML-score is very much unlikely to improve CRT response in the patients. So, our model predictions are in line with clinical observations showing a small fraction of anterior segments among positive responses even for optimal lead positions.
4.5 On the role of scarring area for optimal LV lead position
In our study, the extent of LV myocardial damage (both absolute and relative to the surviving myocardium volume) was not selected as one of the most important features by the ML feature selection algorithms and thus was not included as an independent variable in our predictive LR classifier. At the same time, distance from LV pacing site to scarring area (Scar-LV distance) was selected as the third most important feature for CRT response prediction. It was selected by every feature selection algorithm we tested from the total hybrid dataset analyzed, and was used as one of the 7 most significant input variables for building the predictive LR model (see Fig. S3). This distance is the only model-driven feature that distinguishes responders from non-responders in our population at the referent LV lead position (see Table 2), although, no correlation was found between Scar-LV distance and LV EF improvement (r=0.18, p=0.211).
At the same time, we discovered a low positive correlation between ML score and Scar-LV distance at the referent lead position suggesting that it has a role to play in the integrative estimate of CRT response prediction (r=0.419, p=0.000). As a support for this hypothesis, we revealed strong positive correlations between the change in ML-score and the change in Scar-LV distance for the optimal LV lead position against the referent one (r=0.673 for ML-LP; 0.855 for LAT-LP; and 0.881 at TAT-LP; p<0.01). Here, the absence of difference in average ML-score between LAT-LP, TAT-LP and ref-LP was consistent with no difference in Scar-LV distance between the pacing settings. By contrast, the higher Scar-LV distances were associated with the maximum ML-scores in our patients (57 ± 21 mm in ML-LP versus 32 ± 24 mm in ref-LP, Table. 2).
Our findings are consistent with the results of clinical studies which assessed the significance of postinfarction scar size for CRT response. Marsan and co-authors (60) performed MRI for CRT candidates to derive LV mechanical dyssynchrony and the extent of scar tissue to predict CRT response. Higher LV dyssynchronies were strongly associated with echocardiographic response to CRT, while the total extent of scar correlates with non-response. Importantly, a univariable logistic regression analysis showed that the presence of a match between the LV lead position and a transmural scar was also significantly associated with non-response to CRT. The location of scar in the posterolateral region of the LV, which is empirically thought to be a target site for LV lead implantation, was associated with lower response rates following CRT (17). Pezel and co-authors (16), found no difference in the presence and extent of scar between CRT responders and non-responders. However, in non-responders, the LV lead was more often over an akinetic/dyskinetic area suggesting the presence of tissue lesions, a fibrotic area, or an area with myocardial thickness < 6 mm. By contrast, the extent of scar core and gray zone was automatically quantified using cardiac MRI analysis (15) and the highest percentage of CRT response was observed in patients with low focal scar values and high QRS area before operation. Such area was calculated using vector-cardiography. Lee and co-authors used LV wall thickness <5 mm as an anatomical index for scar, and this feature was not shown to be predictive of an acute hemodynamic response to pacing (37). However, in every patient from their cohort the LV had a wall thickness > 5 mm throughout. So, the lack of MRI information regarding scarring in the cardiac tissue was mentioned as one of the limitations of the data they used. In our models we accounted for such data and showed that this is essential for model predictions.
4.6 On the great potential of using personalized ventricular models for CRT response prediction
In recent simulation studies, personalized cardiac models were used to reveal model-derived features correlating with CRT response (61, 62, 63, 64, 65). In two recent papers, model simulations were demonstrated to be predictive of the LV pacing site optimization (36, 37). In our present study, we have developed an ML-based technique, using both clinical and simulated features. This technique provides an LV surface map predicting areas of positive and negative response and indicating the best possible place for LV lead guidance with a very high probability of CRT response.
Further prospective analysis of indices suggested by simulation studies could enable the clinicians to test acute responses at fewer pacing sites (if any) intra-operatively and during follow up. Thus, it will provide some benefits like reducing time, costs and risks to patients and enhancing the chances of a positive outcome. The results of simulation studies demonstrate the potential of virtual clinical trials as a tool for exploring new approaches for pacing lead placement optimization .
5 LIMITATIONS
There are several limitations in our study. First, the ventricular geometry in our personalized models was derived from CT images obtained after CRT device implantation, not before it. Despite the supposed difference in the ventricular geometry our simulated ECGs in the LBBB mode have a high correlation with pre-operative clinical ECGs (r=0.84, p<0.05), thus demonstrating the effect of ventricular geometry as being secondary. Second, the CRT response definition we used was based on an LV EF improvement that has a low-to-moderate correlation with ventricular reverse remodeling in our patient cohort. It has been recently shown in a simulation study by Rodero and co-authors (36) that cardiac reverse remodeling after CRT implantation can reduce the effect of ventricular pacing. According to this hypothesis, if our predictions classify a patient as positive with post-operative ventricular geometry, it is all the more reasonable to expect to see this participant as positive with preoperative geometry. However, more false negative predictions could be expected, and an additional prospective study with pre-operative ventricular geometry will have to be performed to prove our approach.
Next, we have shown a high importance of the distance from the LV pacing site to the myocardial damage area in ML predictions of optimal LV pacing position. In this study, we simulated LV scarring area based on the labeling of damaged LV segments performed by the expert that analyzed the MRI scans. A more accurate segmentation of the raw MRI data should be used to compare model predictions with a different accuracy of scarring area simulation. We think more objective information on the scar and fibrosis morphology may improve the predictive models of CRT response as well as optimal pacing site location.
In this proof-of-concept study, we used our ML-based approach to optimize BiV pacing configuration. Currently, there are data emerging on the different pacing modality enhancing effects of pacing in a certain patient depending on the ischemic or non-ischemic origin of CHF, LBBB or non-LBBB activation, the levels of conduction system block, etc. New techniques for His-Purkinje and/or conduction system pacing, selective LV, RV or BiV, or multipole pacing, epicardial versus endocardial pacing and their possible combinations together with AV and VV delay optimization for the optimal ventricular fusion create a challenge for the best choice. Only computational modeling provides a tool to test every possible combination and suggest ones that help to optimize patient outcome. The technology we have developed is capable of solving such complex problems.
In the present study, CRT response prediction involved simulated characteristics of ventricular activation and ECG derived from electrophysiological models. However, the synchronization of ventricular contraction and subsequent improvement in the mechanical performance of the ventricles is the main goal of the therapy. Recent studies have shown the predictive power of the mechanical indices that could be measured from CT or echocardiography images and accounted for in the predictive models of optimal pacing designs (37). Moreover, electromechanical models of cardiac activity (such as reported by (65, 64, 62, 63) and being developed by our team) could help perform direct simulations of LV EF, dP/dtmax changes and other mechanical biomarkers of CRT response which can further improve ML based optimization of CRT procedure.
Last but not least, in this study we had a limited data sample from 57 patients. However, to the best of our knowledge, this is the largest model population used in simulation studies. Our predictive classifiers based on hybrid clinical data from and computational models have demonstrated high performance with accuracies higher than those achieved with classifiers developed on the basis of clinical data from a thousand of patients. Still, a prospective study using our technique is needed to validate the approach and confirm its usefulness for patient stratification and optimal lead guidance.
6 CONCLUSIONS
We have developed a new technique combing personalized heart modeling and supervised ML to predict optimal LV pacing lead position in CHF candidates for CRT. We suggest an optimal LV pacing site based on ML-scores from an LR classifier of positive CRT response (LV EF improvement > 10%). In a patient group with suboptimal LV lead position deployed in close proximity to predicted optimal pacing site with maximal ML-score, the number of positive responses is two times higher than negative responses. We showed the distance from the LV pacing site to scarring area to be an important feature for predicting optimal lead location.
This novel approach has great potential clinical implications for patient care improvement. With an ML classifier on hybrid data created and thoroughly validated, the range of generated ML scores at any pacing site throughout the accessible LV surface would classify this patient as a potential responder or non-responder to the therapy, thus supporting individual selection for CRT. At the same time, the best pacing site location predicted from model simulations and corresponding ML scores could be used for guiding lead deployment during the CRT procedure and optimizing the outcome for the patient.
Data Availability
All data produced in the present study are available upon reasonable request to the authors
CONFLICT OF INTEREST STATEMENT
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could conflict of interest
AUTHOR CONTRIBUTIONS
OS, AD, TC, AB, and SK participated in conceptualization and design of the study, data collection, processing, and analysis. SZ, TL, VL and DL performed patient selection and treatment, clinical data collection and preprocessing. SZ, AB, and AD performed segmentation CT and MRI images and construction of anatomical models. AD, SK, AB performed personalized computational simulations of myocardial electrical activation. AD implemented algorithms of the pacing site optimization. SK, AD, TC performed machine learning and built classifiers. TC contributed to the statistical analysis and data interpretation. OS supervised all stages of study execution and data analysis, and was a major contributor to writing the manuscript. All authors contributed to manuscript preparation and approved the final version of the manuscript.
FUNDING
This work was supported by Russian Science Foundation grant No. 19-14-00134.