2.1. Derivation cohort (EAST-AFNET 4)

The early treatment of AF for stroke prevention trial (EAST-AFNET 4) randomized patients with recently diagnosed AF and stroke risk factors to systematic early rhythm control or usual care, including symptom-based rhythm control.¹⁴ All patients were followed up for a median of 5.1 years. The primary outcome of the trial was a composite of stroke, cardiovascular death, and unplanned hospitalization for heart failure or acute coronary syndrome.¹⁴ Details of the EAST-AFNET 4 biosample study collecting a baseline blood sample from 1586 patients enrolled in the EAST-AFNET 4 trial have been published.¹⁶ In brief, all consenting patients provided a blood sample at baseline. Samples were shipped to the core biostorage facility at UKE Hamburg, spun, shock frozen, and stored at −8°C. EAST-AFNET 4, and its biomolecule substudy were approved at all study sites. Written informed consent was obtained from all patients. This study complied with the Declaration of Helsinki.

2.2. Validation cohort (BBC-AF subcohort)

Details of the BBC-AF cohort have been described before.¹⁵ In brief, consecutive patients eligible for recruitment had electrocardiogram (ECG)-diagnosed AF or presented with at least two cardiovascular conditions (congestive heart failure, hypertension, diabetes, prior stroke, or vascular disease) to a large teaching hospital (Sandwell and West Birmingham NHS Trust). Patients who did not have a diagnosis of AF underwent 7-day ambulatory ECG monitoring to rule out undiagnosed ECG-documented AF. For this analysis, only patients with ECG-documented AF were included. All patients underwent a detailed interview, physical examination, 12-lead ECG, and transthoracic echocardiography at the time of recruitment. Follow-up data for events were collected by assessing local hospital records corroborated against Hospital Episode Statistics data, general practitioner records, and mortality data from NHS Digital, at 2.5 years after the final patient was recruited.¹⁷ The follow-up duration was calculated as the time between the baseline assessment date and an event or the date of record review, where no events were documented. This study complied with the Declaration of Helsinki, was approved by the National Research Ethics Service Committee (IRAS ID 97753), and was sponsored by the University of Birmingham. All patients provided written, informed consent.

2.3. Selection of biomolecules and their quantification

The methodological approach taken for this work tried to combine the multifaceted information of several circulating biomolecules and the reliability of high-precision assays. A group of scientists within the EU-funded CATCH-ME consortium searched literature and scouted meetings for disease mechanisms related to AF and to cardiovascular events in patients with AF. One summary of this initial effort has been published.¹ In a next step, biomolecules that could potentially reflect these disease processes (called ‘health modifiers’ in Fabritz et al.¹) were identified based on a literature and patent search enriched with knowledge available in the EU Horizon 2020 CATCH-ME consortium. A modified Delphi expert consensus process was conducted to identify biomolecules representing these disease mechanisms available for high-precision quantification. Thirteen biomolecules were identified (Table ​Table11): angiopoietin 2 (ANGPT2), bone morphogenetic protein 10 (BMP10), cancer antigen 125, C-reactive protein (CRP), D-dimer, endothelial-specific molecule 1, fatty acid–binding protein 3 (FABP3), fibroblast growth factor 23 (FGF23), growth differentiation factor 15 (GDF15), insulin-like growth factor–binding protein 7 (IGFBP7), interleukin-6 (IL-6), N-terminal pro-B-type natriuretic peptide (NT-proBNP), and cardiac troponin T (TnT).

Biomolecules were centrally quantified using pre-commercial and commercial high-throughput, high-precision platforms (Roche, Penzberg, Germany). The biomolecule quantification was provided as an in-kind contribution by Roche to the EU Horizon 2020 CATCH-ME consortium. Absolute protein concentrations were centrally quantified in ethylenediaminetetraacetic acid (EDTA) plasma. Run controls and calibrators were measured twice each run, and lab staff involved were blinded to clinical status and data. Blood samples were shipped to and quantifications were conducted at the Roche biomolecule research facility in Penzberg, Germany. This is the first analysis of the biomolecules in the EAST-AFNET 4 trial substudy.

2.4. Data pre-processing and clustering

Biomolecule concentrations were 1% winsorized and Blom-transformed¹⁸ and each patient was assigned into one quintile for each biomolecule. These quintiles were used to cluster patients using polytomous variable latent-class analysis (poLCA). K-means clustering was used as sensitivity analysis. Unsupervised models were developed using LCA, available in the package poLCA¹⁹ in R (https://www.r-project.org/). LCA was performed on 13 biomarker variables. Patients with any missing biomolecule concentrations were excluded. The models were created in a bootstrapping fashion by repeating 100 times with data resampling with replacement. Each data resampling was performed with a fixed initialization seed to ensure reproducibility. Models between 2 and 10 clusters were assessed. The Bayesian information criterion (BIC) was used to assess the best number of clusters by penalizing models with too many parameters. The number of clusters with the lower BIC score was counted over all bootstraps. This most frequent number of clusters was used to create the final model using the original data, without resampling. We also fitted k-means clustering models against the same Blom-transformed biomolecules without conversion of those into categorical variables. To assess the optimal number of cluster groups, we followed the exact same approach as for poLCA models. As base model instance, we made use of SciKit-learns algorithm (https://scikit-learn.org/stable/index.html) with Lloyd algorithm,²⁰k-means++²¹ as initial cluster centroid selection strategy and let the algorithm run for 10 iterations with different centroid seeds to fit the model against our data set.

2.5. Phenotypic description

Clusters were formed agnostic to any clinical data. Patients’ characteristics were summarized for each cluster. Body mass index (BMI) was categorized into obese (BMI ≥ 30) and non-obese. The estimated glomerular filtration rate (eGFR) was calculated using the Chronic Kidney Disease Epidemiology Collaboration equation^22,23 and categorized into normal kidney function (eGFR ≥ 60 mL/min) and chronic kidney disease. The left ventricular ejection fraction was categorized into groups of ≥50 and <50%. Differences between categorical variables were calculated using generalized logistic mixed models with the study centre as random effect. For continuous variables, linear mixed models with study centre as random effect were used for normally distributed and non-normally distributed variables. P-values resulted from Analysis of Deviance Table (Type II Wald χ² tests). The R packages lme4 and car were used for this analysis. In the non-multicentric BBC-AF cohort, t-test and χ² test were applied for quantification of differences between continuous and categorical variables.

2.6. Commonality analysis

Biomarkers are biologically secreted and reabsorbed reflecting different disease processes. They are excreted or shed based on common pathways (e.g. secretion via the kidney), creating collinearity in the data. The presence of multicollinearity, as can be demonstrated using a correlation matrix or calculating the variance inflation factor, complicates the interpretation of regression model outcomes, since both unique and shared variances are contributing to an effect on the outcome. Commonality analysis allows an exploration of relationships between biomarkers by decomposing the R² of the regression model or respectively the pseudo R² of the binomial regression model to quantify unique and shared variances of each biomarker in explaining the outcome. The analysis returns 2^k − 1 commonality coefficients (k = number of variables entered).

2.7. Dominance analysis

As commonality analysis is one way to assess the relative importance of predictors (p; the 13 biomolecules) on an outcome (the cluster group), dominance analysis (DA) is a different approach that makes two distinct contributions. First, it measures the relative importance of a predictor in a pairwise fashion and secondly, it does this in the context of all $2^{(p-2)}$ models that contain any subset of the remaining predictors. In a refined version of DA by Azen and Budescu,²⁴ the concepts of complete, conditional, and general dominance were introduced. We employed the Dominance-Analysis Python package²⁵ (Python 3: https://www.python.org/) that refines this concept further by introducing individual dominance, average partial dominance, interactional dominance, and total dominance. Finally, the percentage relative importance can be calculated from the mean of those four dominance measures for each predictor variable.

2.8. Change in clustering after removal of biomolecules

To generate a more global understanding of the contribution of each biomolecule to the clustering, the patient clustering was repeated with reduced feature sets by removing one, two, three, four, or five of the biomolecules. For each possible biomolecule combination of the reduced feature sets, we estimated the optimal number of clusters, predicted the clusters, and used those to partition the patients again. We computed the adjusted rand index (ARI) for those new partitions in comparison with the original ones derived from the poLCA model fitted against 13 biomolecules and predicting four clusters. The biomolecule clusters were used as predictors in Cox proportional hazards (PH) model instances, with the first primary composite outcome as an event of interest to obtain hazard ratios and c index.

2.9. Risk of cardiovascular events

Cox PH models were fitted using cluster membership as the predictor to predict a composite outcome of cardiovascular death, heart failure hospitalization, stroke or systemic embolism, and acute coronary syndrome. To infer hazard ratios, we used models with the centre as shared frailty term and the R package survival. For sensitivity analyses, we added age, sex, and randomization group as confounding variables to the models. To compare the unsupervised cluster assignment with existing predictors, separate risk prediction models were built using (i) CHA₂DS₂VASc score (congestive heart failure, hypertension, age, diabetes, stroke, vascular disease, sex category), (ii) ABC stroke¹² and bleeding¹³ scores, (iii) discretized¹⁷ Troponin T, and (iv) NT-proBNP quartiles. For the ABC scores, published criteria^12,13 were computed. The ABC (age, biomarker, clinical history)-stroke and ABC-bleeding risk scores incorporate clinical variables and cardiovascular biomarkers to estimate risk of stroke or systemic embolic events and bleeding, respectively, in patients with AF.

In the BBC-AF validation cohort, there are 68 missing values for the first primary composite outcome and 59 missing values for the first primary safety outcome for either the event-status information or the time-to-event information. We dropped those participants from the primary analysis, but added a best- and worst-case scenario analysis. For the best-case scenario, all missing event values were set to one (occurrence of an event) and for the best-case scenario all event values were set to zero (no occurrence of an event). We imputed all missing time-to-event data by the median censoring time.

To calculate the area under the receiver operating characteristic curve (AUC), we fitted unadjusted Cox PH models without centre as frailty as this information is missing in the validation data set, and our aim was to measure the discriminatory power of inter-cohort generalizable models. We used the Python lifelines²⁶ and Sklearn packages for this analysis. The ABC scores^12,13 and genomic risk scores were used as continuous variables, and all other predictors were discretized as proposed in the literature (details in legend to Supplementary material online, Figure S5).

3.2. Risk of outcome events in each cluster

Each cluster had a distinct risk of primary and safety outcomes (Figures ​Figures33A and ​and33C). Patients in the highest risk cluster had a five-fold higher rate of primary outcomes than patients in the lowest risk cluster in the derivation (Figure ​Figure33A) and validation (Figure ​Figure33B) data sets. Each component of the composite outcome moved in the same direction as the composite for the primary outcomes (Table ​Table44, Figures​33A and​33B) and for the safety outcome (Table ​Table55, Figures ​Figures33C and ​33D). The clustering using k-means resulted in a similar risk gradient across clusters (see Supplementary material online, Figure S3). Early rhythm control was effective across all biomolecule clusters (P_interaction = 0.63, Supplementary material online, Table S5). The clustering model outperformed other risk scores for most of the tested outcomes (see Supplementary material online, Figures S5–S8). For the first primary composite outcome, the poLCA cluster model yielded an AUC 0.76 [95% confidence interval (CI): 0.72–0.79], the next best predictive model relied on the ABC bleeding score and yielded an AUC 0.74 (95% CI: 0.70–0.78) in the validation data set. Hazard ratios for the clustering were higher or similar to hazard ratios obtained by applying established risk prediction models using clinical features, combinations of clinical features and biomolecules, or a single biomolecule (see Supplementary material online, Table S5).

Open in a separate window

Figure 3

Risk of cardiovascular events for each biomolecule cluster. (A) Aalen–Johansen curves for cluster groups from poLCA clustering model for the first primary outcome in EAST-AFNET 4, a composite of all-cause mortality, stroke, or unplanned hospitalization for heart failure or acute coronary syndrome. (B) Aalen–Johansen curves for cluster groups from poLCA clustering model for the first primary outcome in BBC-AF, a composite of all-cause mortality, stroke, or unplanned hospitalization for heart failure or acute coronary syndrome. Administrative censoring has been applied to events occurring after number of at-risk patients dropped below five. (C) Aalen–Johansen curves for cluster groups from poLCA clustering model for the safety outcome in EAST-AFNET 4. (D) Aalen–Johansen curves for cluster groups from poLCA clustering model for the safety outcome in BBC-AF. Administrative censoring has been applied to events occurring after number of at-risk patients dropped below five.

3.3. Biomolecule combinations are required to assign patients to cluster groups

To estimate the relevance of each biomolecule for the assignment to a patient cluster, we computed the ARI and the C statistic for each combination of biomolecules as further post hoc analyses. Removal of five or more biomolecules consistently yielded ARIs of 0.55 or less, indicating that almost half, or more, of the patients were no longer assigned to their original cluster (Figure ​Figure44A). Figure ​Figure44B plots the ARI and the C statistic for each possible model with fewer biomolecules. While risk prediction remains reasonable with only two to three biomolecules (C statistic estimates 0.67–0.69), assignment to cluster requires more information.

Open in a separate window

Figure 4

Effects of removing biomolecules on clusters. (A) Plot of ARIs for each simulated clustering process when using less biomolecules. Each colour indicates one number of biomolecules (2–12) used for clustering. Each dot represents one simulated set of clusters. ARIs are shown for each model in ascending order from right to left. (B) Plot of the ARI, a measure of the assignment of a patient to a biomolecule-derived cluster, and of the corresponding c index, for each virtual model using poLCA clustering with a reduced number of biomolecules (2–12). Each dot represents one clustering model. The colour indicates the number of biomolecules used. Models relying on eight or less biomolecules (yellow, orange, and red colours) consistently yield ARIs below 0.55. Only models using 7–12 biomolecules achieve correct assignment of patients to biomolecule clusters. The c index, a summarized measure of the accuracy of risk prediction, changes only marginally (x axis). (C) Importance of each biomolecule included in this study based on the effect of its removal from the clustering process. For each number of biomolecules (2–12), the five clusters with the lowest rand indices were selected. The missing biomolecules were listed and ranked by number of clusters lacking that biomolecule. This list provides an estimate of the importance of each biomolecule in the clustering process.

To estimate the relevance of each biomolecule for the assignment of patients to the risk clusters, another in silico exercise was performed: For each number of biomolecules, the five models with the lowest ARIs were selected. The missing biomolecules were listed and counted. Biomolecules whose removal often lead to a low ARI were considered relevant for the clustering process. Figure ​Figure44C provides a list of all biomolecules sorted by the number of relevant removals in this exercise.

4.2. Comparison with other risk estimation scores and limitations of the C statistic

As expected, the biomolecule-based clustering process evaluated here provides better risk estimation than the CHA₂DS₂VASc score.¹² Its C statistic was better than or comparable with other proposed risk scores, including the ABC stroke and bleeding scores^12,13 (Figure ​Figure44B). In view of the summative nature of the C statistic, this may not come as a surprise. The UK Biobank provided the first insights into the added value of multiple biological measurements for risk prediction.⁴² Previous work on biomarkers tested the predictive value of a biomolecule when added to clinical characteristics,⁴³ or for single biomolecules on their own,^10,44 prior to combining biomolecules into scores.^12,13 Classical statistical methods, including forward and backward selection, did not identify these biomolecules, probably due to a different handling of shared and common information.

4.3. Comparison with proteomic methods

Proteomic technologies now enable the quantification of thousands of proteins from a small sample of plasma or blood. Earlier iterations of these technologies contributed to the discovery of AF-related biomolecules quantified here, e.g. FGF23,¹⁵ while RNA sequencing contributed to the discovery of BMP10⁹ as a biomolecule of interest in patients at risk of AF. Others used proteomic analyses in all-encompassing analyses of circulating proteins related to cardiovascular events in patients with AF.¹¹ Such proteomic analyses are hypothesis free but necessarily highlight proteins with large concentration ranges that can be quantified with high precision using proteomic technologies. These proteomic methods will be extremely helpful in identifying additional proteins related to AF. Such work is likely to confirm, refine, and extend the present findings. The present analysis pre-selected 13 biomolecules hypothesized to reflect different modifiable disease processes that can be quantified at high precision. These biomolecules were used to identify groups of patients who share pathophysiological patterns based on these biomolecules. This method identifies groups of patients with different predominant disease mechanisms. It may be useful alongside continued hypothesis-free research aiming to identify additional disease mechanisms leading to age-related diseases⁴⁵ and chronic cardiovascular diseases.¹¹