Machine Learning-Based Phenomapping in Patients with Heart Failure and Secondary Prevention Implantable Cardioverter-Defibrillator Implantation: A Proof-of-Concept Study

Background: Previous studies have failed to implement risk stratification in patients with heart failure (HF) who are eligible for secondary implantable cardioverter-defibrillator (ICD) implantation. We aimed to evaluate whether machine learning-based phenomapping using routinely available clinical data can identify subgroups that differ in characteristics and prognoses. Methods: A total of 389 patients with chronic HF implanted with an ICD were included, and forty-four baseline variables were collected. Phenomapping was performed using hierarchical k-means clustering based on factor analysis of mixed data (FAMD). The utility of phenomapping was validated by comparing the baseline features and outcomes of the first appropriate shock and all-cause death among the phenogroups. Results: During a median follow-up of 2.7 years for device interrogation and 5.1 years for survival status, 142 (36.5%) first appropriate shocks and 113 (29.0%) all-cause deaths occurred. The first 12 principal components extracted using the FAMD, explaining 60.5% of the total variability, were left for phenomapping. Three mutually exclusive phenogroups were identified. Phenogroup 1 comprised the oldest patients with ischemic cardiomyopathy; had the highest proportion of diabetes mellitus, hypertension, and hyperlipidemia; and had the most favorable cardiac structure and function among the phenogroups. Phenogroup 2 included the youngest patients, mostly those with non-ischemic cardiomyopathy, who had intermediate heart dimensions and function, and the fewest comorbidities. Phenogroup 3 had the worst HF progression. Kaplan–Meier curves revealed significant differences in the first appropriate shock (p = 0.002) and all-cause death (p < 0.001) across the phenogroups. After adjusting for medications in Cox regression, phenogroups 2 and 3 displayed a graded increase in appropriate shock risk (hazard ratio [HR] 1.54, 95% confidence interval [CI] 1.03–2.28, p = 0.033; HR 2.21, 95% CI 1.42–3.43, p < 0.001, respectively; p for trend <0.001) compared to phenogroup 1. Regarding mortality risk, phenogroup 3 was associated with an increased risk (HR 2.25, 95% CI 1.45–3.49, p < 0.001). In contrast, phenogroup 2 had a risk (p = 0.124) comparable with phenogroup 1. Conclusions: Machine-learning-based phenomapping can identify distinct phenotype subgroups in patients with clinically heterogeneous HF with secondary prophylactic ICD therapy. This novel strategy may aid personalized medicine for these patients.


Introduction
Evidence-based clinical guidelines recommend implantable cardioverter-defibrillator (ICD) for primary and secondary prevention of ventricular arrhythmias (VA) and sudden cardiac death (SCD) [1,2].In fact, most primary prevention ICD recipients derive no benefit in the subsequent follow-up [3][4][5].Therefore, previous studies have focused on identifying risk factors, constructing risk assessment tools, and finding new parameters and modalities to facilitate risk stratification in these patients [4,[6][7][8].In contrast, the benefit of ICD therapy in secondary prevention is more certain due to the established substrate for VA [9].Consequently, little attention has been paid to secondary prophylactic ICD implantations.Risk stratification in these patients is imperative to understand the heterogeneity of disease development, progression, and prognosis [10].
Despite efforts being made, previous studies could not predict VA recurrence in patients undergoing secondary prevention [11,12].The individual risk of VA recurrence is high, irrespective of other factors [13].Thus, in this circumstance, instead of identifying underlying risk factors, segregating patients into different subgroups based on their similarities, called phenomapping, might provide a more feasible solution [14].In theory, patients within a phenogroup share similar baseline features and long-term prognoses, whereas patients within different phenogroups differ in these aspects.Through phenomapping, personalized treatment and management, including programming settings, medication optimization, and VA ablation, can be refined.
This study aimed to first classify patients into different phenogroups using unsupervised machine learning clustering based on principal component analysis and then validate the clinical utility of phenogroups by comparing the baseline characteristics and outcomes (including all-cause death and the first appropriate shock) among the established phenogroups.

Study Population
We retrospectively included all patients with heart failure (HF) who underwent secondary prophylactic singleor dual-chamber ICD implantation at our center between January 2010 and December 2020 (n = 756).The exclusion criteria were as follows: patients with cardiac channelopathies, hypertrophic cardiomyopathy, and congenital heart disease (n = 132); patients with heart failure with preserved ejection fraction (n = 172); and patients without visit after ICD implantation (n = 63).The study was performed in accordance with the Declaration of Helsinki and approved by the Fuwai Hospital.In addition, written informed consent was obtained from all the patients.

Outcome Assessment
All-cause of death and first appropriate ICD shock on a VA event were primary outcomes.Patient survival status was obtained from medical records, death certificates, or phone contact with their family members.Additionally, ICD discharge information was extracted using device interrogation.The patients were required to have planned interrogation at least yearly and unplanned interrogation after sensible device therapy.The ICD programming settings were at the discretion of the treating physicians.No standard protocols were requested; nonetheless, shock therapies were generally set to be delivered after ATP could not terminate the VA events.Patients were censored at the time of their successful VA ablation, and patient survival was further ascertained after their last device interrogations.

Data Collection and Processing
Forty-four baseline variables, including demographics, physical examination, cardiovascular risk factors, laboratory tests, electrocardiogram, echocardiographic parameters, comorbidities, and medications, were collected.All the variables had a missing rate of <3.9% (Supplementary Table 1).Missing values were imputed using random forests.Data standardization was applied to eliminate different scales of variables, and the original form was used for interpretation.

Phenotypic Clustering
To determine the natural intrinsic connection within ICD recipients, all baseline characteristics, excluding the medications (thirty-four variables in total), were used for phenomapping.Hierarchical k-means clustering was chosen for phenotypic clustering [15,16], which includes the two most widely used unsupervised machine learning algorithms: hierarchical clustering and k-means clustering.Hierarchical clustering was first computed to determine the optimal number of clusters and their respective centroids.Subsequently, k-means clustering was performed using these centroids as the initial cluster centers.Both steps combined to make the clustering results more reliable and robust.Therefore, this is also called k-meansconsolidated hierarchical clustering.It has been successfully applied in clustering patients with dilated cardiomyopathy [17].
Specifically, in this study, as it is not directly applicable to mixed data types (continuous and categorical), the factor analysis of mixed data (FAMD) was a preliminary step to clustering [15,16].To obliterate redundant information contained in the raw variables, only principal components extracted using the FAMD with eigenvalues >1 were retained.Based on the selected principal components, hierarchical clustering was performed using the Euclidean distance and Ward's criterion.The optimal number of clusters (1-10) was determined by the Elbow method, using the total within the sum of squares as metrics.Finally, k-means consolidation was performed, leading to the final phenomapping.
The representative variables in each cluster were analyzed using the v test based on the hypergeometric distribution.Briefly, a positive v-test statistic indicates the overrepresentation of a variable, whereas a negative statistic indicates the underrepresentation of a variable.The baseline characteristics and clinical outcomes of the phenogroups were compared to validate the clinical utility of phenomapping.

Statistics
Continuous data were expressed as mean ± standard deviation or median with 25th and 75th percentiles, as ap-propriate; categorical data were presented as numbers and percentages.Baseline characteristics between phenogroups were compared using ANOVA or Kruskal-Wallis test for continuous variables and χ 2 or Fisher test for categorical variables.Cumulative incidences were estimated using the Kaplan-Meier method and compared using the log-rank test.Unadjusted and adjusted relationships between clusters and outcomes were assessed using Cox proportional hazard regression and described using hazard ratios (HR) and 95% confidence intervals (95% CI).Medications were included in the adjusted models.The backward stepwise selection was applied using the Akaike information criterion to obtain the most parsimonious model.Analyses were conducted using R (version 4.1.2,R Foundation for Statistical Computing, Vienna, Austria), mainly through the "miss-Forest", "FactoMineR", and "factoextra" packages.Statistical significance was set at p < 0.05.

Baseline Characteristics
A total of 389 patients with HF with reduced or mildly reduced ejection fraction were identified.Table 1 summarizes the baseline characteristics of the study population during the ICD implantation.The mean age was 60.5 ± 12.5 years.The patients were predominantly male (83.5%), and approximately half had ischemic cardiomyopathy (52.4%).The mean left ventricular ejection fraction (LVEF) was 35.5 ± 8.4%, and 46.3% of patients had New York Heart Association (NYHA) Class III or IV.

Principal Component Analysis and Phenomapping
The scree plot obtained using FAMD analysis is illustrated in Supplementary Fig. 1.The first 12 principal components with eigenvalues >1 and a cumulative explained variance of 60.5% were used for clustering.The first and second principal components accounted for 9.9% and 9.1% of the variation, respectively.The variables that contributed the most to these two dimensions, as illustrated in Supplementary Fig. 2, were cardiomyopathy etiology and renal function parameters.Phenomapping was based on the extracted components.The Elbow method presented in Fig. 1 suggests that the best number of clusters was three in hierarchical clustering.Subsequently, the k-means consolidation was performed.As illustrated in Fig. 2, the patients were stratified into three distinct phenogroups.Phenogroups 1, 2, and 3 included 163, 148, and 78 patients, respectively.

Baseline Characteristics among Phenogroups
The most representative features of each phenogroup are illustrated in Fig. 3.By adding the group differences presented in Table 1, these three clinical phenogroups were well defined.Phenogroup 1 was the oldest (65.2 ± 10.5 years), and nearly all (95.1%) had ischemic cardiomyopathy.These patients were the most likely to be smokers (64.4%) and had the most favorable cardiac dimensions (left ventricular end-diastolic diameter: 62.0 ± 7.6 mm; left atrial diameter: 42.0 ± 5.5 mm; right ventricular diameter: 21.9 ± 2.9 mm) and function (LVEF: 38.3 ± 7.5%) among the phenogroups.Consistent with these findings, these patients had the highest prevalence of cardiovascular and metabolic comorbidities (diabetes mellitus, hypertension, and hyperlipidemia), the lowest levels of Nterminal pro-brain natriuretic peptide (NT-proBNP) [805.0 (329.0-1787.0)pmol/L], and were the least symptomatic Values are presented as the mean ± standard deviation, median (interquartile range), or frequency (%).ACEI/ARB/ARNI, angiotensin-converting enzyme inhibitor/angiotensin receptor blocker/angiotensin receptor-neprilysin inhibitor; BP, blood pressure; BUN, blood urea nitrogen; CLBBB, complete left bundle branch block; CRBBB, complete right bundle branch block; hs-CRP, highsensitivity C-reactive protein; ICD, implantable cardioverter-defibrillator; IVS, interventricular septum thickness; LAD, left atrial diameter; LVEDD, left ventricular end-diastolic diameter; LVEF, left ventricular ejection fraction; MRA, mineralocorticoid receptor antagonist; NT-proBNP, N-terminal pro-brain natriuretic peptide; NYHA, New York Heart Association; PVC, premature ventricular contractions; RVD, right ventricular diameter.1. as assessed using the NYHA (NYHA III/IV, 31.3%).Conversely, phenogroup 2 included the youngest patients (54.3 ± 12.4 years) and was predominately composed of nonischemic cardiomyopathy (86.5%).These patients had intermediate cardiac structure and function; however, the overall best health condition was reflected by the highest hemoglobin, the lowest creatinine and blood urea nitrogen levels, and the lowest burden of comorbidities across phenogroups.Phenogroup 3 exhibited the most severe HF progression, as evidenced by the largest cardiac dimensions, lowest LVEF, and highest rates of atrial fibrillation.In addition, NT-proBNP levels increased approximately threefold compared to the other groups.Moreover, these patients had the highest renal function parameters, indicating the worst cardiac and renal function status among the phenogroups.Based on these features, individuals in phenogroup 3 were more often prescribed diuretics, digoxins, and anticoagulants.However, beta-blocker use was disproportionately low in these patients.

Discussion
In this follow-up study of secondary ICD recipients with HF, we identified three phenotypically and prognostically distinct phenogroups using principal component analysis and unsupervised phenomapping.Phenogroup 1 comprised the oldest patients with ischemic cardiomyopathy and had the highest diabetes, hypertension, and hyperlipidemia burden.Phenogroup 2 was characterized by the youngest age, non-ischemic cardiomyopathy, and lowest burden of comorbidities.Phenogroup 3 had the worst HF progression.In addition, in phenogroups 1-3, heart structure and function deteriorated.In accordance with these findings, the risk of appropriate shock doubled from phenogroups 1 to 3.Not surprisingly, no VA recurrencefree group was observed.The death risk also doubled in phenogroup 3 compared to that in phenogroup 1.These results offer a novel perspective on predicting the longterm outcomes of different HF subgroups that received secondary ICD implantation and might serve as a triage for these patients.
Prior studies using traditional statistical analysis failed to construct risk assessment tools for secondary ICD recipients [11,12].This may be largely explained by the complex pathophysiological interactions underlying VA [13,18].Therefore, instead of finding a specific mode between baseline features and prognosis, unsupervised machine learning phenomapping by classifying patients into several clusters based on their similarities may represent a promising alternative [14].Patients within a cluster, also called a phenogroup, share similar characteristics and prognostic outcomes, whereas patients from different clusters differ [19][20][21].Phenomapping has been utilized extensively in cardiovascular medicine, including the categorization of HF comorbidities [22], classification of patients with HF with preserved ejection fraction [20,21], dilated cardiomyopathy [17], and identification of cardiac resynchronization therapy responders [19,23].To date, there is no gold standard for selecting an optimal algorithm for phenomapping.However, in this study, k-means consolidated hierarchical clustering based on principal components was chosen because of its reliability and effectiveness [15][16][17].Notably, FAMD was used first to extract principal components from mixed data, which also contributed to reducing the redundancy of the raw variables.These factors led to robust clustering, as demonstrated by the unique clinical features and prognoses observed across these phenogroups.
Assessing the risk of VA recurrence and ICD shock is essential to optimize the management of patients with ICD [10,13].On the one hand, although life-saving, the ICD shock is related to a higher risk of death, decreased quality of life, and psychiatric disorders, such as anxiety and depression [1,2].On the other hand, intensified therapies should not be encouraged indiscriminately, as antiarrhythmic medications have adverse side effects and may increase mortality [24], while VA ablation in ICD recipients has been inconsistent with VA recurrence and mortality [25,26].Therefore, it is strongly advised to identify high-risk individuals and take more aggressive actions in these patients, including personalized device programming to reduce appropriate yet unnecessary ICD shocks, innovative drugs (such as sodiumc-glucose cotransporter-2 inhibitors [27]), and catheter ablation to reduce VA recurrence.Aging, male sex [12], VT as a presenting arrhythmia [11,12,28], incomplete revascularization [29,30], and cardiac remodeling [11,29,31] are risk factors for VA occurrence in secondary prevention; however, no risk models have been developed.In contrast, we identified three phenogroups with a gradient increase in shock risk parallel with progressive cardiac remodeling and function, which might aid in real-world risk stratification.This finding is also consistent with previous findings that deteriorated left ventricular dimension and function contribute to the increased risk of VA [11,29,31,32].
In addition to identifying the risk of appropriate shock, this phenomapping can also classify the risk of all-cause death.Phenogroup 3 had a higher mortality risk than phenogroup 1, which is consistent with HF worsening.Conversely, phenogroup 2 had a lower risk of death than phenogroup 1, although the difference was not significant.At first glance, this is obscure.However, individuals in phenogroup 2 were much younger (11 years younger) and had a better overall health status, except for heart conditions.Thus, the rate of non-cardiac death might be lower in phenogroup 2 than that of patients in phenogroup 1, which ultimately contributes to a lower risk of all-cause mortality.Notably, ICD is merely designed to terminate the VA.Nevertheless, it can neither reverse pump failure nor prevent non-cardiac deaths.Furthermore, data have revealed that patients with comorbid hypertension [33], renal dysfunction [34], and diabetes [35,36] derive less or no benefit from ICD implantation.Therefore, it is imperative to minimize non-cardiac deaths among ICD recipients to achieve a mortality benefit.Specifically, for patients in phenogroup 1, optimized control of diabetes, hypertension, hyperlipidemia, and renal function protection should be implemented.
The risk of cardiac death, including pump failure death and SCD, increases as HF worsens [32,37].Therefore, it is reasonable to identify individuals with a higher likelihood of SCD than pump failure to maximize ICD implantation benefits [38][39][40].However, defining the extent of HF progression that meets this condition is challenging.Therefore, phenomapping patients with varying degrees of comorbidities and HF is a feasible alternative strategy.For instance, phenogroup 2 had a lower risk of death but a higher risk of appropriate shock than phenogroup 1.In this case, phenogroup 2 may derive greater benefits from ICD therapy.Notably, as our patients exclusively comprised secondary prevention indications, no subgroup of patients could be deemed as having no need for ICD implantation.Except for identifying those most likely or unlikely to benefit from ICD, phenomapping is also paramount in the historical period when subcutaneous ICD (S-ICD) is gaining momentum.The S-ICD is designed to overcome some transvenous ICD-related complications, such as vascular injury, lead failure, and systemic infections; however, it is not free from device-related complications and inappropriate shocks [41].Moreover, there remains a need for a certain population implanted with S-ICD to upgrade to a transvenous ICD due to the need for anti-bradycardia pacing, anti-tachycardia pacing, or cardiac resynchronization therapy [42].Phenomapping techniques can identify patients suitable for S-ICD implantation.
Our study has some limitations.It was conducted exclusively in a tertiary hospital; therefore, the conclusions drawn not be directly applicable to external datasets.Therefore, validation of phenomapping is required.Nonetheless, as these phenogroups are well characterized and resemble real-world situations, phenomapping can be inferred to be robust.Second, no standard protocol for device programming was used, and parameter tuning might have occurred during the follow-up.Nevertheless, shock therapies were generally triggered after ATP failed to terminate the VA.Therefore, an appropriate shock could be considered a reliable surrogate for lethal VA.Third, due to the difficulties of adjudication, we did not distinguish between deaths caused by HF progression, SCD, or non-cardiac causes.However, determining the cause of death would undoubtedly provide a more comprehensive perspective for interpreting this phenomapping.Fourth, the patients were undertreated with guideline-directed medical therapy for HF at baseline.Nevertheless, given that phenomapping is solely based on baseline features, except for medications, and computed by similarities between individual patients, outcome differences between phenogroups were unlikely to change significantly even after modifying medications.Finally, only the routine clinical variables were included.Genetic testing, inflammatory markers, computed tomography, and cardiac magnetic resonance parameters were unavailable.This hinders the possibility of diving into a deeper understanding of the phenomapping.

Conclusions
As a proof-of-concept study, we have demonstrated that unsupervised machine learning phenomapping based on principal components is a promising way to identify subgroups of patients with HF implanted with ICD, among whom clinical characteristics and prognosis substantially differed.This novel approach has the potential to facilitate risk stratification and guide the individualized management of patients with ICD.

Fig. 1 .
Fig. 1.Building the hierarchical clustering.(A) The Elbow method to identify the ideal number of clusters.The optimal number of clusters was three, as illustrated.(B) The dendrogram construction.

Fig. 2 .
Fig. 2. Visualization of the phenomapping after k-means consolidation.The X and Y axes represent the first and second extracted principal components, respectively.

Fig. 3 .
Fig. 3. Characteristics of the three phenogroups by their major representative features.A higher v-test statistic means a more remarkable representation of a feature compared to the overall population.Positive and negative values represent overrepresentation and underrepresentation, respectively.Abbreviations as in Table1.

Fig. 4 .
Fig. 4. Kaplan-Meier curves for the clinical outcomes stratified by phenogroups.(A) The first appropriate shock.(B) The all-cause death.

Table 2 . Association between phenogroups and clinical outcomes.
⋆ Phenogroup 1 was used as the reference.The adjusted analysis included the following medication prescriptions: ACEI/ARB/ARNI, amiodarone, beta-blockers, calcium channel blockers, diuretics, MRA, digoxin, statin, antiplatelet, and anticoagulants.HR, hazard ratio; CI, confidence interval.The other abbreviations are listed in Table1.