RETRO is a multicentre, randomized, open, prospective, controlled parallel-group study. Details of the study are described in the original publication (15). The objective of the study is to evaluate tapering or discontinuation of DMARDs in patients with RA. All enrolled patients fulfilled the 2010 American College of Rheumatology (ACR) criteria for RA (19). Patients had to have sustained clinical remission defined by the Disease Activity score (DAS28 < 2.6) criteria for at least 6 months (20). Ethics committee of the Friedrich-Alexander-University of Erlangen-Nuremberg approval was granted.

Patients were randomized to one of three arms: Arm 1 continued with existing DMARD regimen at full dose for 12 months, arm 2 reduced the dose of all DMARDs by 50%, while arm 3 reduced the dose of all DMARDs by 50% in the first 6 months, then discontinued all medications. Relapse of disease was defined as a DAS28-ESR score greater than 2.6. Participants were assessed for clinical disease activity every 3 months until month 12.

Age and sex were recorded in all patients. Disease duration, tender joint count (68), swollen joint count (66), patient visual analogue scale (VAS) for pain and patient global were assessed and recorded. C-reactive protein (CRP), ESR, Rheumatoid Factor (RF), ACPA, DAS28-ESR and Health Assessment Questionnaire (HAQ-DI) were recorded.

SOMAscan is a proteomics assay that measures 1307 proteins using an aptamer library. This high-throughput proteomics assay has been used in recent publications to study the aging proteome (21, 22) along with other human diseases (23). 130 baseline serum samples were available from the RETRO study. Briefly, a library of aptamers were incubated with serum, and those that bind are isolated and hybridized to DNA microarray for detection. The identity and relative concentration of the detected proteins are revealed by localization and fluorescence intensity. Protein quantification is reported as relative fluorescence units (RFU), an arbitrary value. In general, agreement between aptamer and antibody-based assays is high (24). Further details regarding the SOMAscan assay are available (18).

Descriptive results ( Table 1 ) are stated in means and standard deviation. SOMAmer protein expression RFU values for the study patients were transformed into a log2 scale for differential analysis ( Supplemental File 2 ). Batch effect was removed in our SOMAmer data using internal controls within each plate to adjust proteomic intensity as per standard SOMAscan protocols. Batch effect was assessed between plates and determined to require no further correction. Data was loaded and analyzed in the R (v3.5.3) environment unless otherwise stated. Missing clinical data was imputed using multiple imputation by chained equations (MICE) (25). Differential analysis between groups was undertaken using linear modeling with the package LIMMA ( 26). GO pathways analysis was performed using clusterprofiler (27). Graphs were generated using the ggplot2 package. Correlation analyses were performed for select proteins using Pearson correlation. Multi-dimensional scaling (MDS) was used for dimension reduction on all SOMAscan proteins.

The 200 most variable proteins measured by coefficient of variation were used to determine optimal number of clusters ranging from k = 2 to 10 and identify sample clusters using the R package Consensusclusterplus. We used 80% protein resampling and 80% patient resampling and selected Pearson as our distance function. Multinomial logistic regression implemented in R package glmnet was used to identify clinical variables that are independently associated with cluster assignment. Sliding window analysis of DAS28 scores and protein expression was performed using a previously published algorithm, DE-SWAN (28). Briefly, this algorithm analyzed serum protein expression across quintiles of DAS28 scores using linear modeling, while adjusting for baseline demographic factors, in this case age and sex. A protein expression score was developed on proteins that correlated with DAS28 which were identified by DE-SWAN. We filtered the proteomic data on the 34 score members, scaled the data by the mean and standard deviation, and multiplied by 1 (positively associated with DAS28) or -1 (negatively associated with DAS28) for each protein. The final score was the mean expression of all 34 proteins for each patient. We randomly generated 5000 data sets with 34 randomly selected proteins in each set to evaluate the significance of the association score.

We applied two supervised machine learning (ML) techniques to develop algorithms to classify flare or remission based on serum proteomics. The first approach we used is XGBoost (Extreme Gradient Boosting), which employs a regularization term to overcome the overfitting (29). The second approach is the LASSO model (30), which was used as a baseline to compare its performance with that of XGBoost. Data was loaded into Python, and samples were randomly split into a training (n = 104, 80% of the samples) and test (n = 26) set. The training set was used to train and tune the parameters in the two models and the test set was used evaluate the models’ performance, which were measured by the area under of the curve (AUC) of receiver operating characteristic (ROC), accuracy, sensitivity and specificity. To increase the interpretability of the XGBoost model to predict the flare status of a given sample, we used SHAP values (Shapley Additive Explanation) (31). A higher SHAP value of a given feature in the model represents its strong influence on the model output. The final model parameters we used in the XGBoost are as follows: learning_rate = 0.01, max_depth = 3, subsample = 0.6, colsample_bytree = 0.7, n_estimators = 100, gamma = 0.0, reg_alpha = 0.5, the parameters used in the LASSO model is as follow: cost=1.17 and max_iterations = 5000. We used 5-fold cross-validation to get the optimal hyperparameters.

Baseline characteristics for 130 patients enrolled in the RETRO study are found in Table 1 . Overall, the group had maintained clinical remission for 16.6 (16.2) months and mean disease duration of over 6 years. 57.7% of the patients were ACPA positive, while 40.0% required biologics to achieve remission. 76.6% of patients had achieved the most stringent definition of remission [ACR/EULAR remission (32)]. After 12 months of follow-up, 62.3% of the overall population remained in clinical remission (50% in those undergoing withdrawal).

Given the paucity of data aimed at understanding subclinical disease activity in RA patients who achieve remission, we sought to explore underlying heterogeneity using serum proteomics in this established cohort. We quantified over 1300 serum proteins from 130 RETRO patients at their baseline visit, all of whom were in stable clinical remission (DAS28 < 2.6). We hypothesized that despite the clinical similarities amongst individuals within this cohort, proteomic differences may provide important insights by identifying sub-clusters of patients. We applied consensus clustering to assign individuals to one of the 4 clusters ( Figure S1 ) and clustered scaled protein expression by hierarchical clustering, which can be seen in Figure 1 . MDS analysis revealed separation of the hierarchical clusters ( Figure S2 ).

Baseline characteristics split by cluster are listed in Table 1 . We found no differences in sex, age, biologic use, or serological status across our 4 proteomic clusters. Cluster 4 had significantly lower DAS28 scores compared to the remaining clusters. With respects to future flare, Cluster 2 trended toward lower rates relative to the remaining clusters, however this did not reach statistical significance (16.7% vs 39.8%, p = 0.21). To assess this association by multinomial regression, we assigned Cluster 2 as the reference cluster and found that that Cluster 3 had higher odds of flare (OR 5.6, 0.97 to 33.06, p = 0.05), relative to Cluster 2 with similar trends observed for Cluster 1 and Cluster 4 ( Table S1 ). Indeed, no clear distinction between individuals who developed future flare was observed in the MDS plot ( Figure S2 ). Overall, these results suggest that global proteomic clusters within a clinically homogenous cohort can be identified but are associated with current clinical status rather than future outcomes.

We next aimed to use the serum proteome to identify biomarkers associated with future disease flare in stable RA, given that clustering did not clearly associate with risk of flare. We identified DEP’s between these groups ( Table S2 ) and observed upregulation in Ectodysplasin A receptor (EDAR, FC = 1.20) and Serine peptidase inhibitor (SPINT2, FC 1.1), and downregulation of Fractalkine (CX3CL1, FC 0.95) and Ephrin type-B receptor 2 (EPHB2, FC 0.95) in individuals who eventually went on to flare ( Figure 2A ). However, after adjustment for multiple comparisons, none of the differentially expressed proteins reached statistical significance. This suggests that although subtle differences exist in the serum proteome between individuals who experience future flare, it’s unlikely that singular biomarkers accurately predict this outcome in this population.

To test this hypothesis, we explored the use of two machine learning algorithms, LASSO and XGBoost, to build predictive models that classify future flare using baseline serum proteomics. We generated two models, both of which were validated on a test cohort (20% of total cohort). The LASSO model achieved 69.2% accuracy, with an AUC of 0.58, along with sensitivity of 0.5 and specificity of 0.78 based on the test cohort (Features in Table S3 ). XGboost delivered a model with higher specificity (0.78) than sensitivity (0.63) and an overall accuracy of 73.1% with an AUC 0.80. Therefore, we found that the XGboost model outperformed the LASSO model by the metric area under the curve ( Figure 2B , AUC, 0.80 vs 0.58), accuracy (73.1% vs 69.2%) and sensitivity (0.63 vs. 0.5).

To interpret this XGBoost model, we explored the impact of essential features in terms of SHAP values on the classifier’s output for a single prediction which are shown in Figure 2C . We identified Interleukin 17F (IL17F) and Myc proto-oncogene protein (MYC) expression as indicators of future flare, while Roundabout homolog 3 (ROBO3), Synuclein alpha (SNCA), complement factor B (CFB) and vascular endothelial growth factor A (VEGF-A) expression were indicators of sustained remission ( Figure 2C and Figure S3 ). Given the small number of proteins that derived our boosted model, we next explored whether there were any functional links between these proteins. We developed gene concept plots to identify potential protein interactions and found that SNCA, MYC, VEGFA and PLAUR were connected by a single pathway, endopeptidase mediated apoptosis. We then analyzed IL17F restricted networks, given its conflicting role in RA (33–36), and that its expression was associated with future flare in our model. We found that IL17F interacted with VEGFA though growth factor function, and with SNCA through common effects on the inflammatory response ( Figure 2D ). IL17F independently regulated GM-CSF production, a key driver of RA disease activity through the recruitment of neutrophils (37). Overall, this network analysis suggests that cellular apoptosis and GM-CSF production may be associated with future disease flares in RA patients who are otherwise stable.

In our hierarchical clustering, we observed a lower mean DAS28 score in Cluster 4 compared to the remaining 3 clusters ( Figure 3A ). RA patients who achieve DAS28 defined remission often have residual disease activity, however, since this population is not typically the focus of translational studies, little is known regarding biomarkers that are reflective of ongoing disease activity. Indeed, we found that several protein members correlated with DAS28 score ( Figure S4 ), including Integrin alpha 2B (ITGA2B), Bactericidal permeability-increasing protein (BPI) and chemokine ligand 2 (CXCL2, Pearson R, all p value < 0.01). Further, complement proteins (C3, C4A, C1S) were all negatively associated with DAS28 score, suggesting activation and consumption of complement proteins ( Figure 3B ) were indicators of disease. There was no indication that these parameters varied based on ACPA status ( Figure S5 ).

To further explore the relationship between disease activity and serum biomarkers, we developed a sliding window model (SWAN) which examined protein variability across DAS28 quintiles, after controlling for Age and Sex. Across DAS28 a total of 34 proteins varied significantly with disease activity ( Figure 4A ). We used these 34 protein members to annotate a meta-protein expression score (the mean expression profiles of the 34 proteins), which correlated with DAS28 (R = 0.45, p < 0.001, Figure S6 ). To test the robustness of this finding, we sampled 5000 random sets with 34 proteins in each set and correlated their mean expression with DAS28 scores. We found a range of -0.37 – 0.27, associated with a low probability (0) that the correlation of 0.45 would occur by chance ( Figure S7 ). This protein disease activity score was significantly lower in Cluster 4 compared to the remaining 3 clusters ( Figure 4B ), concordant with their lower DAS28 scores. Gene concept plots revealed that these 34 proteins interacted through nodes that included humoral immunity, apoptosis, and complement activation ( Figure 4C ). Overall, these data suggest that stable RA disease activity is marked by complement consumption and humoral immune responses, which is reflected in a serum protein signature that is detectable in patients with stable RA.