The chronic heart failure evolutions: Different fates and routes
An accurate and individualized prognostic assessment of patients suffering from chronic heart failure (HF) in order to identify those at higher risk of complications and ominous outcome is still an unmet need. A better patients' stratification will allow identification of candidates who need treatment upgrading and tailoring, personalization of follow‐up and optimized resource allocation. Currently, clinical definitions classify HF according to left ventricle ejection fraction and HF aetiology while several parameters are used to define HF severity including plasma level of B‐type natriuretic peptides, exercise derived data, haemoglobin concentration, iron deficiency or overload, kidney function, presence of HF co‐morbidities and, more recently, proteomic derived data to assess inflammation, redox state, and coagulation. , , , , Regardless, at present, none of the several HF prognostic tools alone or in combination consents to precisely define a reliable HF evolution model. , , , , Moreover, on top of prognosis, it is important to properly phenotyping HF patients to understand who of the HF patients' plethora will negatively evolve toward a higher risk condition and possibly how fast this evolution will be. With this regard, of major importance is the identification of specific conditions or parameters changes leading to or being a marker of a change in HF evolution trajectory. Indeed, characterizing subgroups of patients at variable risk and HF evolution is crucial to improve the decision‐making process in HF. In this context, the application of unsupervised learning methodologies could help the identification of specific patient phenotypes with predictable HF evolution and possibly with a predictable speed of this evolution. Among them, topological data analysis (TDA) is a robust and effective unsupervised methodology, able to preserve the intrinsic characteristics of data and the mutual relations among observations, depicting complex data in a graph‐based representation. Indeed, building topological models as networks, TDA allows complex diseases to be inspected in a continuous space, where subjects can ‘fluctuate’ over the graph, sharing, at the same time, more than one adjacent node of the network. Of note, TDA has proved to be able to identify specific ‘communities’ of patients revealing similar clinical characteristics, for example, those diagnosed with diabetes, coronary artery disease, or aortic stenosis. , , , The purpose of the present work was a multi‐step analysis of a sizable population of HF patients studied in a holistic mode including demographic, clinical, therapy, laboratory, cardiac ultrasound, and exercise data with a precise predefined follow‐up. At first, we aimed to identify HF patients' phenotypes (communities) who share several characteristics. Second, we wanted to identify the possible evolution pathways of each HF patients' community and, finally, to assess the risk of adverse events in each HF patients' community. 

The present study considers the MECKI registry population, enrolled between 1993 and 2022 from 26 Italian sites for a bulk of 7948 patients with an average follow‐up of 4.09 [1.78–7.47] and unfixed HF severity. In the present analysis, a minimum 2‐year follow‐up was considered since the original MECKI score refers to a 2‐year prognosis. Inclusion criteria for the MECKI study are the following: previous or present HF symptoms (NYHA functional class I–IV) and history or presence of left ventricular ejection fraction (LVEF) < 40%, unchanged HF medications for at least 3 months, ability to perform a CPET, and no major cardiovascular treatment or intervention scheduled. Exclusion criteria were history of pulmonary embolism, moderate‐to‐severe aortic and mitral stenosis, pericardial disease, severe obstructive lung disease, exercise‐induced angina and significant electrocardiographic (ECG) alterations, or presence of any clinical co‐morbidity interfering with exercise performance. A requirement for inclusion in the registry was a CPET performed with the following modalities: (1) a ramp exercise on an electronically braked cycle ergometer with a protocol set, to reach peak exercise in 8 to 12 min ; (2) symptom limited tests, regardless of the respiratory exchange ratio (RER) reached. Ventilation and respiratory gases were collected breath‐by‐breath and analyzed following a standard technique. Oxygen uptake at peak exercise (peakVO 2 ) was calculated as the 20‐s average of the highest recorded VO 2 , while ventilation vs. CO 2 production slope (VE/VCO 2 slope) was calculated as the slope of the linear relationship between VE and VCO 2 from one minute after the beginning of loaded exercise to the end of the isocapnic buffering period. Peak VO 2 expressed as % of predicted value (VO 2 %) was calculated according to Hansen et al. In addition to CPET‐related variables, the MECKI registry collects echocardiographic, ECG, therapy, and blood chemistry data at enrolment, as well as records vital status and causes of events during follow‐up. Patient follow‐up and data management procedures were performed as previously described. In brief, follow‐up was carried out according to the local HF programme, and it ended with the last clinical evaluation or with patients' death, heart transplantation, and left ventricular assist device (LVAD) implantation. If a patient died after hospital discharge, medical records of the event and the reported cause of death were collected. For prognosis evaluation, the endpoint was the composite of cardiovascular death, urgent heart transplantation, or LVAD implantation. The study was approved by the local ethics committee (protocol number: CCM04_21 PA). The present multi‐steps study consists in a data pre‐processing phase to consolidate the MECKI score dataset, in an unsupervised analysis aimed at identification of the study cohorts by TDA, communities' characterization and trajectory inference and, finally, in a supervised analysis for risk assessment and results validation ( Figure ). Topological data analysis TDA needs some parameters to be properly assessed and set. The first 2 dimensions extracted from the principal component analysis (PCA) have been chosen as lenses, while the number of TDA bins (NB) and the bin overlapping ratio (BO) have been selected via a ‘grid search’ approach, ranging NB between , , , , , , , , , , , , and BO between [0.2–0.5]. To ensure a scale‐free configuration of the resulted network, the lowest Pearson's R correlation index between node rank (k) and the frequency of node rank (p(k)) has been chose as evaluation metrics. The best parameters' set was: NB = 28 and OB = 0.2 (see details on the Supplementary File). TDA was performed by the ‘PIUMA’ R package ( https://bioconductor.org/packages/devel/bioc/html/PIUMA.html ). Networks have been built and enriched by Cytoscape; each node represents a group of patients with similar characteristics, while edges thickness depicts the number of samples shared by nodes (Jaccard Index). Edges with a Jaccard Index > 0.1 were considered; on the other hand, nodes collecting <2 patients or without any connections were used as a validation cohort (see below). The Girvan–Newman algorithm implemented in the ‘clustermaker2’ Cytoscape plug‐in, in which edges are iteratively removed based on the value of their betweenness, was used to identify communities (i.e., clusters) in the network. Of note the cluster labels were randomly assigned. Trajectory analysis Trajectory analysis allows deriving paths, that is, trajectories, that describe the evolution of samples traversing through different states and, at the same time, estimating a temporal unit, called pseudo‐time, which describes the progression from a state to the next one; this method was firstly proposed to model single cell RNA‐seq experiments. Of note, during their HF journey, patients may progressively increase or reduce their pseudo‐time in case of HF worsening or improvement, respectively. Here, we assumed that each cluster is a disease state. Network nodes were projected onto a 2D space by the ‘Uniform Manifold Approximation and Projection’ (UMAP) dimension reduction method, implemented in the ‘umap’ (v.0.2.10) R package. Then, the ‘slingshot’ R package was exploited to identify the clusters' centroids, assess clusters' relationships, determine trajectories and, then, infer pseudo‐times. Overall, this analysis allows inferring'state‐by‐state’, the sequence and the progression of a disease, iteratively assessing the lowest difference between all adjacent clusters' characteristics. Survival analysis Survival analyses were performed, exploiting the ‘survival’ (v. 3.2.13) and ‘survminer’ (v. 0.4.9) R packages. Kaplan–Meier analysis was implemented to generate time‐to‐event curves up to five‐year follow‐up in patients with and without reaching the study endpoint, recorded in the MECKI registry. All patients included in the study (study cohort and internal + external validation cohorts, Figure ) were considered in the survival analysis. A P ‐value < 0.05 was considered statistically significant. Internal and external validation The cohort for the internal validation was built with patients not resulted to be connected to any nodes of the network generated by TDA ( n = 527). On the other hand, the cohort for the external validation was gathered by collecting data from European and Chinese centres, joining the MECKI consortium activities on December 2022 ( n = 777). Data of the latter patients were recently published. For both the validation cohorts, patients were assigned to a specific cluster based on their features' similarity with the clusters' hallmarks, identified by the TDA on the study cohort; specifically, for each cluster, a coefficient of variation was calculated for each feature, using the averaged values estimated in the study cohort as reference. The cluster with the highest number of coefficients of variation < 5% was assigned to the patient. We performed two validation approaches: firstly, a Pearson's correlation analysis was implemented to compare event frequency occurrence at 5 years in each cluster of the study cohort with those in each validation cohorts, to assess the reproducibility of TDA procedure in different settings. Second, we exploited an event‐free survival analysis (Kaplan–Meier) to evaluate if the survival trend of each cluster were similar between the study cohort and the validation cohort; for this second approach, we join the internal and external validation cohort to allow a more reliable 5‐year Kaplan Meyer analysis. Statistical analysis For continuous variables, the Anderson–Darling test has been performed to assess the hypothesis of data normality. Continuous variables, not normally distributed were presented as median and interquartile range, while categorical data are reported as frequencies and percentages. Group comparisons for continuous and categorical variables were performed by Mann–Whitney U test and by chi‐square ( χ 2 ) test, respectively. On the hand, normally distributed variables were reported as mean ± standard deviation and compared by unpaired t‐test. The relationship between variables has been assessed by Pearson's R correlation coefficient. A P ‐value < 0.05 was considered statistically significant. All plots were generated by the ‘ggplot2’ R package. The data imputation was performed by the ‘randomForest’ R package. 

TDA needs some parameters to be properly assessed and set. The first 2 dimensions extracted from the principal component analysis (PCA) have been chosen as lenses, while the number of TDA bins (NB) and the bin overlapping ratio (BO) have been selected via a ‘grid search’ approach, ranging NB between , , , , , , , , , , , , and BO between [0.2–0.5]. To ensure a scale‐free configuration of the resulted network, the lowest Pearson's R correlation index between node rank (k) and the frequency of node rank (p(k)) has been chose as evaluation metrics. The best parameters' set was: NB = 28 and OB = 0.2 (see details on the Supplementary File). TDA was performed by the ‘PIUMA’ R package ( https://bioconductor.org/packages/devel/bioc/html/PIUMA.html ). Networks have been built and enriched by Cytoscape; each node represents a group of patients with similar characteristics, while edges thickness depicts the number of samples shared by nodes (Jaccard Index). Edges with a Jaccard Index > 0.1 were considered; on the other hand, nodes collecting <2 patients or without any connections were used as a validation cohort (see below). The Girvan–Newman algorithm implemented in the ‘clustermaker2’ Cytoscape plug‐in, in which edges are iteratively removed based on the value of their betweenness, was used to identify communities (i.e., clusters) in the network. Of note the cluster labels were randomly assigned. 

Trajectory analysis allows deriving paths, that is, trajectories, that describe the evolution of samples traversing through different states and, at the same time, estimating a temporal unit, called pseudo‐time, which describes the progression from a state to the next one; this method was firstly proposed to model single cell RNA‐seq experiments. Of note, during their HF journey, patients may progressively increase or reduce their pseudo‐time in case of HF worsening or improvement, respectively. Here, we assumed that each cluster is a disease state. Network nodes were projected onto a 2D space by the ‘Uniform Manifold Approximation and Projection’ (UMAP) dimension reduction method, implemented in the ‘umap’ (v.0.2.10) R package. Then, the ‘slingshot’ R package was exploited to identify the clusters' centroids, assess clusters' relationships, determine trajectories and, then, infer pseudo‐times. Overall, this analysis allows inferring'state‐by‐state’, the sequence and the progression of a disease, iteratively assessing the lowest difference between all adjacent clusters' characteristics. 

Survival analyses were performed, exploiting the ‘survival’ (v. 3.2.13) and ‘survminer’ (v. 0.4.9) R packages. Kaplan–Meier analysis was implemented to generate time‐to‐event curves up to five‐year follow‐up in patients with and without reaching the study endpoint, recorded in the MECKI registry. All patients included in the study (study cohort and internal + external validation cohorts, Figure ) were considered in the survival analysis. A P ‐value < 0.05 was considered statistically significant. 

The cohort for the internal validation was built with patients not resulted to be connected to any nodes of the network generated by TDA ( n = 527). On the other hand, the cohort for the external validation was gathered by collecting data from European and Chinese centres, joining the MECKI consortium activities on December 2022 ( n = 777). Data of the latter patients were recently published. For both the validation cohorts, patients were assigned to a specific cluster based on their features' similarity with the clusters' hallmarks, identified by the TDA on the study cohort; specifically, for each cluster, a coefficient of variation was calculated for each feature, using the averaged values estimated in the study cohort as reference. The cluster with the highest number of coefficients of variation < 5% was assigned to the patient. We performed two validation approaches: firstly, a Pearson's correlation analysis was implemented to compare event frequency occurrence at 5 years in each cluster of the study cohort with those in each validation cohorts, to assess the reproducibility of TDA procedure in different settings. Second, we exploited an event‐free survival analysis (Kaplan–Meier) to evaluate if the survival trend of each cluster were similar between the study cohort and the validation cohort; for this second approach, we join the internal and external validation cohort to allow a more reliable 5‐year Kaplan Meyer analysis. 

For continuous variables, the Anderson–Darling test has been performed to assess the hypothesis of data normality. Continuous variables, not normally distributed were presented as median and interquartile range, while categorical data are reported as frequencies and percentages. Group comparisons for continuous and categorical variables were performed by Mann–Whitney U test and by chi‐square ( χ 2 ) test, respectively. On the hand, normally distributed variables were reported as mean ± standard deviation and compared by unpaired t‐test. The relationship between variables has been assessed by Pearson's R correlation coefficient. A P ‐value < 0.05 was considered statistically significant. All plots were generated by the ‘ggplot2’ R package. The data imputation was performed by the ‘randomForest’ R package. 

The analysis of the present study consists of three consecutive steps: ‘step 1’ was dedicated to the selection of suitable patients and features (data pre‐processing), ‘step 2’ to unsupervised analysis for community identification, characterization and trajectory inference, and ‘step 3’ to a supervised analysis for risk assessment based on evaluated communities and findings validation ( Figure ). Step 1: Baseline characteristics of study population Before implementing TDA, we pre‐processed the original MECKI registry dataset (7948 samples and 51 clinical features) as follows: (1) features with >10% missing values were removed; (2) samples with >5 missing values or <2‐year follow‐up without a MECKI outcome (non‐cardiovascular death) were discarded; (3) data imputation was performed on the remaining dataset (4876 samples and 43 features) ( Figure ). Imputation was needed for 1, 2, 3, 4, and 5 parameters in 989, 429, 286, 190, and 84 patients, respectively. In 2898 patients, all 43 clinical features were available and no imputation was needed. Three hundred and ninety‐nine (8%) and 789 (22%) patients reached the study endpoint at 2‐ and 5‐year follow‐up, respectively. At 5‐year 1235 cases had <5 years follow‐up without events. Both at 2‐ and 5‐year follow‐up patients with events were more likely male, older, with atrial fibrillation, lower blood pressure and LVEF, larger LV volumes, longer QRS duration, worst kidney function, lower haemoglobin concentration and cardiopulmonary exercise performance. As regards therapy, patients with events have a lower prescription rate of ACE/AT1/ARNI blockers but higher MRA, amiodarone and diuretic prescription ( Table ). Step 2: Identification and characterization of patients' communities Based on the 43 clinical features ( Table ), TDA generated a network‐based representation of the cohort with 629 nodes linked by 1088 edges, aggregating 4349 samples. Conversely 175 nodes, collecting 527 patients who, albeit with similar HF characteristics, were not closely linked to the graph; these patients were chosen as the cohort for the internal validation. Upon the network, 19 patients' communities (i.e., cluster of nodes) were, then, automatically identified ( Figure ). Interestingly, distinct sets of categorical and continuous features are predominantly associated to specific clusters ( Figure most relevant variables, Table all variables and ). For instance, cluster ‘cl.01’ and ‘cl.02’, albeit with a similar NYHA class (2.80 ± 0.40 and 2.87 ± 0.26, respectively), are differently characterized. The cluster ‘cl.01’ grouped patients with the highest frequency of exercise oscillatory ventilation (average frequency = 65%), and with the lowest values of VO 2 peak % (32.8 ± 6.94), Na + (137 ± 1.73), haemoglobin (11.9 ± 0. ±0.9 g/dL), and MDRD (43.7 ± 13.6 mL/min/1.73m 2 ). On the other hand, patients belonging to the cluster ‘cl.02’ showed the lowest values of LVEF (19.4 ± 2.01), K + (4.27 ± 0.50 meq/L) and peak exercise respiratory rate (29.1 ± 3.55). Inference of the evolution states of heart failure At the enrolment, HF severity of patients varied from moderate to severe. This allows hypothesizing that the clusters identified by TDA could represent specific stages of HF evolution, from the lowest to the worst severity. In addition, a trajectory analysis is able to identify all the routes connecting one or more starting clusters to several endpoint clusters. This necessarily implies that the pathways (i.e., the disease evolution) in some clusters can follow at least two different directions. We referred to those clusters as ‘bifurcation points. Here, by the trajectory analysis, a path was identified ( Figure ), starting from cluster ‘cl.16’ (average pseudo‐time: 2.8 ± 2.6). In the HF evolution path, 3 bifurcation points (red squares) at clusters ‘cl.13’ (pseudo‐time: 8.1 ± 1.9), ‘cl.09’ (pseudo‐time: 10.1 ± 1.3), and ‘cl.07’ (pseudo‐time: 12.8 ± 1.4) were identified. Finally, 4 stopping points (yellow diamonds) at clusters ‘cl.15’ (pseudo‐time: 9.1 ± 0.9), ‘cl.06’ (pseudo‐time: 14 ± 3.4), ‘cl.01’ (pseudo‐time: 17 ± 0.8), and ‘cl.02’ (pseudo‐time: 18 ± 0.8, Figure and Figure ) ended the HF evolution trajectory. The stopping points can be considered specific end‐stages of the disease and are associated to peculiar characteristics ( Table ) and patients’ prognosis. Indeed, the 2‐year study endpoint prevalence was 56%, 36%, 9%, and 2%, the 5‐year study endpoint prevalence was 73.5%, 60%, 25%, and 6%, while average MECKI score at enrolment was 20%, 13%, 11%, and 10% for cluster ‘cl.01’, ‘cl.02’, ‘cl.06’, and ‘cl.15’, respectively. Differently, the bifurcation points represent key states to determine distinct disease evolutions; therefore, we identified the features that are significantly different comparing the samples belonging to the bifurcation point with those in the adjacent clusters ( Table ). In Table , the bifurcation cluster is reported in the middle while on the left and right the less and more adverse evolution, respectively. For example, for patients in cluster ‘cl13’ a negative HF evolution, that is, evolution to cluster ‘cl.09’, is foreseen for patients who develop atrial fibrillation, increase NYHA class, reduces kidney function and peak VO 2 with a simultaneous increase of the VE/VCO 2 slope while a more favourable evolution, that is, to cluster 12, is foreseen for patients who had a higher QRS duration, underwent a pacemaker therapy and had a more frequent MRA treatment. Similarly, at the second bifurcation, that is, ‘cl.09’, a negative HF evolution is expected in patients who reduced BMI, haemoglobin, kidney function and increase resting heart rate and frequency of exercise induced oscillatory breathing; differently a less negative evolution is suggested by the presence of reduced digitalis, MRA therapy and ICD, higher LVEF and short QRS. At the third bifurcation a negative HF evolution is foreseen in patients with increase in heart rate and NYHA class as well as a decrease in K + , NA + level, and VO 2 peak % as these patients will follow the route through the ‘cl.01’ end‐stage; differently, patients following the other route (‘cl.02’ end‐stage) showed a significant increase of haemoglobin levels ( P < 0.001). Step 3: Survival analysis and MECKI study endpoint (cardiovascular death, urgent heart transplant, or left ventricular assist device) by communities So far, we identified HF patients' communities in an unsupervised way, that is, without knowing any disease outcomes. However, the MECKI registry recorded also a prolonged follow‐up. Therefore, we performed a survival analysis, stratifying data by communities. Interestingly, ‘cl.01’ and ‘cl.02’, which were also identified as end‐stages in the trajectory analysis, showed the lowest survival rate at 2 (44.1% and 62.5%) and 5 year (20% and 35.2%, respectively). On the other hand, patients belonging to ‘cl.11’, ‘cl.12’, ‘cl.13’, ‘cl.14’, ‘cl.15’, ‘cl.16’, ‘cl.17’, ‘cl.18’, and ‘cl.19’ (the first stages in trajectory analysis) showed a survival rate higher than 90% at a 5‐year follow‐up ( Figure ). In addition, we confirmed that each cluster can be considered as a specific state of heart failure progression, exploiting the blood concentrations of B‐Type Natriuretic Peptide (BNP), a well‐known biomarker of Heart Failure severity, available at the baseline for 2840 patients (58%) and therefore not included in the initial analysis. A Pearson's correlation analysis highlighted a relevant negative relationship ( R = −0.96, P < 0.001) between the average concentrations of BNP with the event‐free survival rates at 5 years of each cluster ( Figure ). It is acknowledged that we have not considered NTproBNP but BNP only. Indeed, NTproBNP was available only in 1085 cases. Moreover, we have not included in the present subanalysis patients on ARNI treatment due to the direct influence of ARNI on BNP value. Validation To confirm that reliability of our finding we performed an internal and an external validation study. To do so, we inferred the cluster membership for the 527 patients of the internal validation cohort and 777 of the external validation cohort ( Figure ). Interestingly, the percentage of events in each cluster for both the validation cohorts is highly correlated to that of the study cohort, being the Pearson's correlation index R = 0.94 ( P < 0.001, Figure ) and R = 0.84 ( P < 0.001, Figure ), when comparing the internal and the external validation cohort, respectively. Of note, the survival analysis was performed grouping patients according the 19 clusters and considering the HF population of the two validation cohorts combined ( n = 1304). The results were comparable with those of the original study population ( Figure and ); indeed, comparing the Kaplan–Meier curves trend for each cluster in both cohorts, we confirmed that there no significant difference between any cluster ( P > 0.05), except cluster ‘cl.06’ ( Figure ). Looking deeper into the difference between the ‘cl.06’ in both cohorts, we confirmed that the high majority of variables are not statistically different. However, we noted in the validations cohort a significant increase of the average rest heart rate (74.7 vs. 68.1 b.p.m.), and of the atrial fibrillation patients' frequency (32% compare to 4% in study population). Moreover, survival differences between the two ‘cl.06’ cohorts are mainly related to the first year ( Figure ). 

Before implementing TDA, we pre‐processed the original MECKI registry dataset (7948 samples and 51 clinical features) as follows: (1) features with >10% missing values were removed; (2) samples with >5 missing values or <2‐year follow‐up without a MECKI outcome (non‐cardiovascular death) were discarded; (3) data imputation was performed on the remaining dataset (4876 samples and 43 features) ( Figure ). Imputation was needed for 1, 2, 3, 4, and 5 parameters in 989, 429, 286, 190, and 84 patients, respectively. In 2898 patients, all 43 clinical features were available and no imputation was needed. Three hundred and ninety‐nine (8%) and 789 (22%) patients reached the study endpoint at 2‐ and 5‐year follow‐up, respectively. At 5‐year 1235 cases had <5 years follow‐up without events. Both at 2‐ and 5‐year follow‐up patients with events were more likely male, older, with atrial fibrillation, lower blood pressure and LVEF, larger LV volumes, longer QRS duration, worst kidney function, lower haemoglobin concentration and cardiopulmonary exercise performance. As regards therapy, patients with events have a lower prescription rate of ACE/AT1/ARNI blockers but higher MRA, amiodarone and diuretic prescription ( Table ). 

Based on the 43 clinical features ( Table ), TDA generated a network‐based representation of the cohort with 629 nodes linked by 1088 edges, aggregating 4349 samples. Conversely 175 nodes, collecting 527 patients who, albeit with similar HF characteristics, were not closely linked to the graph; these patients were chosen as the cohort for the internal validation. Upon the network, 19 patients' communities (i.e., cluster of nodes) were, then, automatically identified ( Figure ). Interestingly, distinct sets of categorical and continuous features are predominantly associated to specific clusters ( Figure most relevant variables, Table all variables and ). For instance, cluster ‘cl.01’ and ‘cl.02’, albeit with a similar NYHA class (2.80 ± 0.40 and 2.87 ± 0.26, respectively), are differently characterized. The cluster ‘cl.01’ grouped patients with the highest frequency of exercise oscillatory ventilation (average frequency = 65%), and with the lowest values of VO 2 peak % (32.8 ± 6.94), Na + (137 ± 1.73), haemoglobin (11.9 ± 0. ±0.9 g/dL), and MDRD (43.7 ± 13.6 mL/min/1.73m 2 ). On the other hand, patients belonging to the cluster ‘cl.02’ showed the lowest values of LVEF (19.4 ± 2.01), K + (4.27 ± 0.50 meq/L) and peak exercise respiratory rate (29.1 ± 3.55). 

At the enrolment, HF severity of patients varied from moderate to severe. This allows hypothesizing that the clusters identified by TDA could represent specific stages of HF evolution, from the lowest to the worst severity. In addition, a trajectory analysis is able to identify all the routes connecting one or more starting clusters to several endpoint clusters. This necessarily implies that the pathways (i.e., the disease evolution) in some clusters can follow at least two different directions. We referred to those clusters as ‘bifurcation points. Here, by the trajectory analysis, a path was identified ( Figure ), starting from cluster ‘cl.16’ (average pseudo‐time: 2.8 ± 2.6). In the HF evolution path, 3 bifurcation points (red squares) at clusters ‘cl.13’ (pseudo‐time: 8.1 ± 1.9), ‘cl.09’ (pseudo‐time: 10.1 ± 1.3), and ‘cl.07’ (pseudo‐time: 12.8 ± 1.4) were identified. Finally, 4 stopping points (yellow diamonds) at clusters ‘cl.15’ (pseudo‐time: 9.1 ± 0.9), ‘cl.06’ (pseudo‐time: 14 ± 3.4), ‘cl.01’ (pseudo‐time: 17 ± 0.8), and ‘cl.02’ (pseudo‐time: 18 ± 0.8, Figure and Figure ) ended the HF evolution trajectory. The stopping points can be considered specific end‐stages of the disease and are associated to peculiar characteristics ( Table ) and patients’ prognosis. Indeed, the 2‐year study endpoint prevalence was 56%, 36%, 9%, and 2%, the 5‐year study endpoint prevalence was 73.5%, 60%, 25%, and 6%, while average MECKI score at enrolment was 20%, 13%, 11%, and 10% for cluster ‘cl.01’, ‘cl.02’, ‘cl.06’, and ‘cl.15’, respectively. Differently, the bifurcation points represent key states to determine distinct disease evolutions; therefore, we identified the features that are significantly different comparing the samples belonging to the bifurcation point with those in the adjacent clusters ( Table ). In Table , the bifurcation cluster is reported in the middle while on the left and right the less and more adverse evolution, respectively. For example, for patients in cluster ‘cl13’ a negative HF evolution, that is, evolution to cluster ‘cl.09’, is foreseen for patients who develop atrial fibrillation, increase NYHA class, reduces kidney function and peak VO 2 with a simultaneous increase of the VE/VCO 2 slope while a more favourable evolution, that is, to cluster 12, is foreseen for patients who had a higher QRS duration, underwent a pacemaker therapy and had a more frequent MRA treatment. Similarly, at the second bifurcation, that is, ‘cl.09’, a negative HF evolution is expected in patients who reduced BMI, haemoglobin, kidney function and increase resting heart rate and frequency of exercise induced oscillatory breathing; differently a less negative evolution is suggested by the presence of reduced digitalis, MRA therapy and ICD, higher LVEF and short QRS. At the third bifurcation a negative HF evolution is foreseen in patients with increase in heart rate and NYHA class as well as a decrease in K + , NA + level, and VO 2 peak % as these patients will follow the route through the ‘cl.01’ end‐stage; differently, patients following the other route (‘cl.02’ end‐stage) showed a significant increase of haemoglobin levels ( P < 0.001). 

So far, we identified HF patients' communities in an unsupervised way, that is, without knowing any disease outcomes. However, the MECKI registry recorded also a prolonged follow‐up. Therefore, we performed a survival analysis, stratifying data by communities. Interestingly, ‘cl.01’ and ‘cl.02’, which were also identified as end‐stages in the trajectory analysis, showed the lowest survival rate at 2 (44.1% and 62.5%) and 5 year (20% and 35.2%, respectively). On the other hand, patients belonging to ‘cl.11’, ‘cl.12’, ‘cl.13’, ‘cl.14’, ‘cl.15’, ‘cl.16’, ‘cl.17’, ‘cl.18’, and ‘cl.19’ (the first stages in trajectory analysis) showed a survival rate higher than 90% at a 5‐year follow‐up ( Figure ). In addition, we confirmed that each cluster can be considered as a specific state of heart failure progression, exploiting the blood concentrations of B‐Type Natriuretic Peptide (BNP), a well‐known biomarker of Heart Failure severity, available at the baseline for 2840 patients (58%) and therefore not included in the initial analysis. A Pearson's correlation analysis highlighted a relevant negative relationship ( R = −0.96, P < 0.001) between the average concentrations of BNP with the event‐free survival rates at 5 years of each cluster ( Figure ). It is acknowledged that we have not considered NTproBNP but BNP only. Indeed, NTproBNP was available only in 1085 cases. Moreover, we have not included in the present subanalysis patients on ARNI treatment due to the direct influence of ARNI on BNP value. Validation To confirm that reliability of our finding we performed an internal and an external validation study. To do so, we inferred the cluster membership for the 527 patients of the internal validation cohort and 777 of the external validation cohort ( Figure ). Interestingly, the percentage of events in each cluster for both the validation cohorts is highly correlated to that of the study cohort, being the Pearson's correlation index R = 0.94 ( P < 0.001, Figure ) and R = 0.84 ( P < 0.001, Figure ), when comparing the internal and the external validation cohort, respectively. Of note, the survival analysis was performed grouping patients according the 19 clusters and considering the HF population of the two validation cohorts combined ( n = 1304). The results were comparable with those of the original study population ( Figure and ); indeed, comparing the Kaplan–Meier curves trend for each cluster in both cohorts, we confirmed that there no significant difference between any cluster ( P > 0.05), except cluster ‘cl.06’ ( Figure ). Looking deeper into the difference between the ‘cl.06’ in both cohorts, we confirmed that the high majority of variables are not statistically different. However, we noted in the validations cohort a significant increase of the average rest heart rate (74.7 vs. 68.1 b.p.m.), and of the atrial fibrillation patients' frequency (32% compare to 4% in study population). Moreover, survival differences between the two ‘cl.06’ cohorts are mainly related to the first year ( Figure ). 

To confirm that reliability of our finding we performed an internal and an external validation study. To do so, we inferred the cluster membership for the 527 patients of the internal validation cohort and 777 of the external validation cohort ( Figure ). Interestingly, the percentage of events in each cluster for both the validation cohorts is highly correlated to that of the study cohort, being the Pearson's correlation index R = 0.94 ( P < 0.001, Figure ) and R = 0.84 ( P < 0.001, Figure ), when comparing the internal and the external validation cohort, respectively. Of note, the survival analysis was performed grouping patients according the 19 clusters and considering the HF population of the two validation cohorts combined ( n = 1304). The results were comparable with those of the original study population ( Figure and ); indeed, comparing the Kaplan–Meier curves trend for each cluster in both cohorts, we confirmed that there no significant difference between any cluster ( P > 0.05), except cluster ‘cl.06’ ( Figure ). Looking deeper into the difference between the ‘cl.06’ in both cohorts, we confirmed that the high majority of variables are not statistically different. However, we noted in the validations cohort a significant increase of the average rest heart rate (74.7 vs. 68.1 b.p.m.), and of the atrial fibrillation patients' frequency (32% compare to 4% in study population). Moreover, survival differences between the two ‘cl.06’ cohorts are mainly related to the first year ( Figure ). 

In the present study, utilizing a dedicated analytical approach, we showed that HF progression is characterized by a continuous but heterogeneous evolution process. The HF evolution leads to four terminal stages differentiated as for characteristics and prognosis. During HF progression, patients meet 3 bifurcation points, which determine HF evolution and outcome. Accordingly, special attention needs to be dedicated to patients when their HF trajectory reaches a bifurcation point since the subsequent evolution leads to a different outcome. Of note, albeit HF is usually seen as a progressive disease internodes trajectories are potentially bidirectional with a reverse direction in case of HF recovery. The study population was sizable with several variables derived from anthropometric, clinical, laboratory, therapy, cardiac ultrasound, and cardiopulmonary exercise testing data. Each of the analyzed variable has a recognized and frequently independent prognostic role in chronic HF. The study patients' population was recruited by a large number of HF centres. Patients had mostly moderate HF based on survival, NYHA class, exercise performance, cardiac volumes and LVEF, kidney function, haemoglobin levels, blood pressure and incidence of atrial fibrillation. Of note, to be eligible patients were required to have a history of HF with reduced LVEF, but this condition could not be present at the time the CPET was performed. Treatment was optimized in accordance with HF guidelines in place at study time. As a matter of fact, beta‐blockers, angiotensin converting enzyme inhibitors/angiotensin receptor blockers/angiotensin receptor‐neprilysin inhibitor (ARNI), and mineralcorticoid receptor antagonist were prescribed in a relevant fraction of patients as were ICD and CRT devices. We required the existence of >2 years follow‐up as in the original MECKI algorithm and assessed patients' prognosis up to 5 years. Indeed, for a HF population with mostly moderate HF, a 2‐ and 5‐year timeframe seems to us a reasonable length of time, being a shorter interval more appropriate for patients with more severe HF and a longer one not particularly relevant due to several possible intercurrent confounding events. Consequently, albeit the continuous updating of the MECKI database, we evaluated patients recruited up to 2020, while patients recruited afterwards were not considered in the present analysis. We recognize that this period was arbitrary and that potential prognostic benefit of new HF drugs such as gliflozins and ARNI was therefore not considered or limited, respectively. Accordingly, a further re‐evaluation of HF trends will be needed after the widespread and more prolonged introduction of those drugs as suggested by 2021 ESC guidelines. In the era of precision medicine, a relevant unmet need is still the accurate identification of homogeneous subgroups of patients in heterogeneous cohorts ; this is, particularly evident for the heart failure. Today, many unsupervised artificial intelligence techniques could help face this issue ; one of these, the TDA, has proved to be an effective in identify sub‐phenotypes of patients in cardiovascular settings. Indeed, the TDA has been recently used for phenotyping patients with diabetes and aortic stenosis as well as to describe the evolution of coronary artery diseases. , , In the present study, we applied an unsupervised topology‐based approach to identify HF patients' clusters with peculiar characteristics and estimate HF evolution; indeed, exploiting the pseudo‐time estimation, an established technique to estimate the cell fate evolution in single cell RNA‐Seq data, we inferred the continuous progression of HF disease states, assuming a cluster as a disease stage. The pseudo time analysis allows to evaluate in a comparable fashion the evolution speed of the different HF settings. Overall and most importantly, we showed that HF progression is characterized by a continuous but heterogeneous evolution process in which many bidirectional trajectories, both toward a disease progression but also toward recovery, are conceivable. Specifically, the TDA allowed the identification of 19 groups of patients with peculiar clinical characteristics and the identification of possible path of HF evolution. Indeed, it is important to identify the cluster to which the patient belong both to precisely define his/her prognosis but also to understand in which direction the disease is evolving. In addition, we described that the HF progression leads to four terminal clusters associated with different characteristics and prognosis. For instance, comparison of ‘cl.01’ and cl.02, showed albeit a similar NYHA class a different phenotype with highest frequency of exercise oscillatory ventilation, lowest values of haemoglobin, peak VO 2 and kidney function in the former and a lower LVEF in the latter suggesting a mainly ‘peripheral’ phenotype in ‘cl.01’ and ‘central’ in ‘cl.02’. Moreover, during HF progression, patients encounter 3 bifurcation points, which determine distinct HF evolution and outcome. Accordingly, special attention needs to be dedicated to patients when their HF trajectory reaches a bifurcation point since the subsequent evolution leads to a different outcome. Therefore, the identification of a cluster which has multiple evolution paths is of relevant clinical impact. Of note, albeit HF is usually seen as a progressive disease, internodes trajectories are potentially bidirectional with a reverse direction in case of HF recovery. The present cluster analysis was validated considering an internal cohort of 527 patients. These patients were recruited in the same time period of the original study population and a relevant concordance between events was observed. A second validation cohort ( n = 777) was obtained from the MECKI international registry and again a relevant concordance between events was observed. Of note, the two implemented validation analysis confirmed the reproducibility of our results. All together, these findings imply that the phenotyping of patients in 19 clusters is reliable as regards events. Clearly further validations are needed particularly since a few new, effective, disease modifying drugs have been introduced as new is the therapeutic approach of the latest guidelines. On this regards it is of note that, applying the present analysis to data derived from a recently published paper on MECKI score reassessment where patients were evaluated twice, we ascertained that either patients, who improved (23.7% of tested patients) or worsened (40% of tested patients) the MECKI score, followed the above described HF trajectories ( Figure ). Moreover, in another experimental settings, where patients, not included in the main study and with a short follow‐up, were evaluated before and after ARNI ( n = 47) or dapagliflozin ( n = 19) treatment, we observed that all patients who improved their condition (43.4%), followed backward the HF trajectories ( Figure ). Again, these findings albeit very preliminary, reinforce the reliability of the present analysis. Consequently, we believe that the application of our or similar HF evolution models will allow a more comprehensive patients phenotyping, a more appropriate patients management and consequently resource allocation. Indeed, the patients at the bifurcation points seems to us those who need a stricter follow up being at a point where the evolution path may change. The present investigation suffers of a few limitations, which need to be acknowledged. First, we analysed only patients with history of low LVEF. Hence, patients with preserved and mildly reduced LVEF have not been tested systematically. Second, the number of female HF patients, 903 (19.0%) and 679 (18.9%) at 2‐ and 5‐year follow‐up, respectively was relatively small and not tested independently. Similarly, the great majority of our patients were Caucasian, so that we do not know how the present model performs in other races. Third, albeit we analyzed several variables some of those known to have an impact on HF prognosis such as sleep disorders, genetics and family history, prevalence of catheter ablation for atrial fibrillation or cardiac rehabilitation were not considered. On this regard, BNP, probably the most used HF prognostic and diagnostic biomarker was not included in the present analysis because available only in 2840 out of 4348 cases. Indeed, in the early cases BNP measurement was not routinely performed in all centers and, in the most recent cases only NT‐proBNP was available. Of note, in each cluster BNP value was strongly associated with the 5‐year event free survival rate ( Figure ) supporting the consistency of the implemented unsupervised approach. Fourth, the effect of treatment changes or of events needs to be assessed in larger populations. Indeed, the few cases we analysed ( Figure . B) are just a tentative and should be seen as a work in progress and a very preliminary analysis which needs to confirmed. Fifth, patients' recruitment period was large, expanding from 1993 to 2022, a period during which significant advancement in treatment modalities and diagnostic features happened. However, we assessed the impact of time of patients recruitment on our findings enriching the TDA network using the time of enrollment, and found no influence. Furthermore, it must be recognized that the present HF evolution model needs to be tested and validated by an independent revision. As a matter of facts, a sequential assessment of clinical data was not available so that the real heart failure progression paths still need to be evaluated. Indeed, although we believe that the present approach represents a relevant step forward in the evaluation of HF management it must be recognized that a perfect prognostic algorithm is still unreached and a continuous update involving more patients with different HF phenotypes and therapies is mandatory. In other words, HF prognosis evaluation is an endless story which remember to me a famous poem … espera luz espera, Y corro ansiosio loco. Espera luz espera, Y cuando voy a legar a su lado, se oscurece, fria, … (Juan Ramòn Jimènez). Research in the context This study, leveraging TDA on a cohort of 4876 chronic heart failure (HF) patients, identified 19 distinct HF phenotypes with unique progression trajectories. Crucially, bifurcation points and stopping points were showed, signifying pivotal junctures and end‐stage conditions, respectively. Survival rates varied significantly among clusters, offering a nuanced understanding of HF prognosis. The robust validation of findings in internal and external cohorts reinforces the clinical applicability of our approach, providing a foundation for personalized HF patient management. Our findings redefine the landscape of HF prognostication, offering a novel framework for tailoring assessments to individual patients. The identified bifurcation points highlight critical decision junctures, necessitating focused clinical attention. Looking forward, this work paves the way for future research into targeted interventions at these key points, aiming to improve outcomes and reshape the management of chronic HF. The identified phenotypes serve as a foundation for developing precision medicine approaches in HF care, emphasizing the potential for more effective and individualized treatment strategies. Moreover, the described HF evolution pathways may be reverted in case of therapeutic intervention. The latter needs to be validated by further larger studies. 

This study, leveraging TDA on a cohort of 4876 chronic heart failure (HF) patients, identified 19 distinct HF phenotypes with unique progression trajectories. Crucially, bifurcation points and stopping points were showed, signifying pivotal junctures and end‐stage conditions, respectively. Survival rates varied significantly among clusters, offering a nuanced understanding of HF prognosis. The robust validation of findings in internal and external cohorts reinforces the clinical applicability of our approach, providing a foundation for personalized HF patient management. Our findings redefine the landscape of HF prognostication, offering a novel framework for tailoring assessments to individual patients. The identified bifurcation points highlight critical decision junctures, necessitating focused clinical attention. Looking forward, this work paves the way for future research into targeted interventions at these key points, aiming to improve outcomes and reshape the management of chronic HF. The identified phenotypes serve as a foundation for developing precision medicine approaches in HF care, emphasizing the potential for more effective and individualized treatment strategies. Moreover, the described HF evolution pathways may be reverted in case of therapeutic intervention. The latter needs to be validated by further larger studies. 

The study was supported by Italian Ministry of Health (Ricerca Corrente; CUP=B43C24000090001). 

None declared. 

Figure S1. Evolution states of HF by pseudo‐time analysis. As in Figure 4A the dot plot shows the UMAP representation of TDA network nodes with pseudo time analysis. For symbols see Figure 4A. Figure S2. BNP analylis. Linear correlation between BNP (log) vs. 5 year survival defined as the composity of cardiovasualr death, urgent heart transplant or LVAD implantation. Figure S3. Survival analysis. Kaplan Meyer curves in each cluster comparing the study population and the 2 validation cohorts combined. The study endpoint was the combination of cardiovascular death, urgent heart transplant or LVAD. Figure S4. Tentative application of TDA in clinical settings. A) A first tentative application of our finding in a cohort of patients who repeating CPET after 1–2 years from the first examination. As expected, most of samples worsened (40%, red arrows) during time, following the route we identified in trajectory analysis, while a few of them improved their condition (23.7%; blue arrows). B) A second tentative application in a small cohort of samples who repeated CPET after 6–9 months treated with ARNI or dapagliflozin. About half (56.6%) of samples did not change their condition while all the other improved their status (43.4%; blue arrows). In both diagrams, each node represent the cluster centroid as in Figure 4A, while arrows could have specific direction whether patients improve (blue arrows) or worsen (red arrow) their condition. Table S1: feature distributions for each cluster. Data S1. Supporting Information. 

