Clinical Validation of a Multi-protein, Serum-based Assay for Disease Activity Assessments in Multiple Sclerosis

Patient samples

A total of 614 serum samples were included from two different sources: 448 retrospective samples from the Serially Unified Multicenter MS Investigation [SUMMIT] consortium; 166 prospective samples from the Rocky Mountain Multiple Sclerosis Clinic (RMMSC). SUMMIT samples were sent to Octave Bioscience, Inc. from three independent sites and studies: Comprehensive Longitudinal Investigation of Multiple Sclerosis at the Brigham and Women’s Hospital (CLIMB; n=195), American University of Beirut Medical Center Study (AMIR; n=202), and University of Massachusetts MS Center—Family Study of Demyelinating Disease (FSDD; n=51). Patient samples selected for this study were intentionally enriched for the presence of Gd+ lesions compared with the general MS population.

SUMMIT was an international multicenter, prospectively enrolled MS cohort study with standardized data structure and analysis groups that can be stratified by demographics, clinical measures, disease relapses, MRI measures, or blood sampling.²² Samples from RMMSC were collected as part of the matched serum and MRI for the Disease Activity Test Development Study. Cross-sectional (samples from a single time point) and longitudinal samples (samples from multiple time points from the same patient) from both SUMMIT and RMMSC were included in the analysis. Serum specimens were collected using standard venipuncture and processing protocols. Samples were transferred to Octave Bioscience, Inc. for analysis and stored at −65□C.

Clinical and radiographic data for biostatistical analysis

An annotated data set containing clinical and radiographical data was obtained for each serum sample. Individual data sets with additional demographic, clinical, and radiographic fields were collected; a combined data set was matched and adjusted so that the annotations were consistently labeled across the study sites and stratified, meaning that the data were split into the Train and Test subsets with a balanced distribution across demographic characteristics, sample counts per site, and Gd+ lesion counts.

Bioanalytical analysis

All serum samples were analyzed as part of a single experiment performed over several days. Assay plates contained up to 72 samples analyzed in a single well each; four serum controls, three calibrators, and a blank control, were described previously²¹ and assayed in triplicate. Analytical runs were stratified independently of the clinical Train/Test randomization, ensuring that assay plates had balanced distribution across sites, demographics, and presence/absence of or number of Gd+ lesions. Analysis was performed cross-sectionally, using the protein concentrations as predictors and algorithmic features.

Pre-processing and quality control was performed using the Olink^® Normalized Protein Expression (NPX) Manager software (Olink Proteomics, Uppsala, Sweden). Built-in quality control (three internal controls that were added into all samples and the external controls) enabled control over the technical performance of the assay at each step of the analysis. These internal controls consisted of an incubation control, extension control, and a detection control. Quality control was performed per assay run and for individual samples at each step of the analysis. Standard deviations for each internal control were established to be below a predetermined threshold (ie, 0.2 NPX) for the entire plate. Median values were calculated for the incubation and detection controls, respectively, for a sample plate. The result of each internal control was required to be within ±0.3 NPX from the plate median. If any of these internal controls deviated from this range, the sample failed quality control and was reanalyzed.

External controls consisted of serum pools with endogenous protein concentrations established at expected levels. Acceptability of a plate run was based on the percent recovery of the serum pools relative to their expected values (ie, ±3 SD). Individual samples or plates that failed the analytical quality control process were rerun.

Assessments

The primary and exploratory assessments of the study included evaluation of the association of the multi-protein and single-protein models to the radiographic (Gd+ and N/E T2 lesions) and clinical (active vs stable disease status) assessments. Radiographic annotations were derived from brain MRIs for all patients in the study and from spine MRIs when available. Association of single-protein and multi-protein models with the presence or absence of Gd+ lesions, as determined on a matched MRI performed within 60 days of the blood draw, was considered the primary endpoint. Next, the final models were evaluated relative to additional disease activity exploratory assessments (N/E T2 lesions and active vs stable disease status).

Samples were considered active if any Gd+ lesion was present, if any N/E T2 lesion was present, or if there was evidence of a clinical relapse within 30 days. Samples were otherwise considered stable (including samples missing N/E T2 lesion/clinical relapse data). This approach was used to further refine the overall Disease Activity and Disease Pathway algorithms (eg, Immunomodulation, Neuroinflammation, Myelin Biology, and Neuroaxonal Integrity). Association of the results from multi-protein and single-protein modeling with disease progression was an additional exploratory assessment, which is not reported in this manuscript. Classification models were fit to the data with proteomic results as the independent variables and presence or absence of Gd+ lesions as the dependent variable.

Statistical analysis

The study hypothesis was that a multi-protein model would significantly associate with clinical and radiographic disease activity endpoints and be superior to the highest performing, single-protein model based on the protein biomarkers included in the multi-protein model.

Inclusion criteria for the study required that information for the primary endpoint reference standard (eg, presence or absence of Gd+ lesions) was available for all patient samples. Upon completion of analytical quality control processes and any necessary reanalysis, index test data (eg, protein concentrations) were available for each sample. For secondary and exploratory endpoints, samples with missing data were excluded from the statistical analysis.

The entire data set (N=614) was split into Train and Test subsets. The Train subset was designed to optimize algorithms and included 70% of the total available samples; the Test subset established the performance specifications of the MSDA Test and comprised 30% of the total available samples. The subsets were stratified to ensure a balanced distribution across demographic characteristics that included age, sex, and disease duration; sample counts per site; and Gd+ lesion counts. Analysis was performed in the Train subset to assess unexpected differences in biomarker ranges resulting from preanalytical processing at the four sites from which the samples were obtained.

A two-layer stacked classifier using L2-penalized logistic regression models, which leveraged biological categorizations of biomarkers to calculate Disease Pathway scores, was ultimately developed. A score-based algorithm enabled the four Disease Pathway scores and an overall Disease Activity score to be derived from the probabilities from the Pathway and Disease Activity model and calculated for individual samples (Figure 1). Full details on the algorithm and model parameters are presented elsewhere.²³

• Download figure
• Open in new tab

Figure 1.

MSDA Test stacked classifier flow chart and biological characterizations model. The first layer of the model consisted of four Disease Pathway algorithms. The second layer of the model utilized the four Disease Pathway probabilities along with the individual age- and sex-adjusted biomarker concentrations as meta-features to determine an overall Disease Activity score that reflected the probability of disease activity. Thresholds were established based on the count of Gd+ lesions for the Disease Activity score, which corresponded to low (1.0–4.0), moderate (4.5–7.0), or high (7.5–10.0) levels of disease activity. MSDA, multiple sclerosis disease activity.

Single-protein models were fit using L2-penalized logistic regression with presence or absence of Gd+ lesions as the dependent variable and an intercept and the protein biomarker as independent variables.

Protein concentrations were limit of quantitation (LOQ)-imputed, log₁₀-transformed, and demographically adjusted (with age and sex, based on Ordinary Least Squares [OLS] modeling) prior to being used in the Disease Pathway, Disease Activity, and single-protein models. Sex and age were selected as demographic adjustment variables (if there was a dependence with the protein biomarkers) since they are routinely collected in a clinical setting when blood samples are taken for analysis. Previous research and development studies, samples from a cohort of healthy controls, and those from the Train subset were used to establish the biomarker-specific demographic adjustment strategy, which included removing protein concentration outliers, accounting for OLS coefficient sign consistency across the three studies and establishing statistical significance related to both age and sex.²³

Metrics for model performance including the area under the receiver operating characteristic (AUROC), sensitivity, specificity, positive predictive value (PPV), negative predictive value (NPV), accuracy, and odds ratios were used to evaluate model performance. The prevalence of Gd+ lesions was enriched in this dataset, and it is important to note that PPV, NPV, and accuracy all depend on the prevalence. An L2-penalization was used to optimize the model and minimize overfitting when training to ensure generalizability to the Test subset. The multi-protein classification model was compared with the top-performing single-protein model (demographically corrected neurofilament light polypeptide chain [NfL] for all disease activity endpoints) to assess the statistical significance of differences in AUROC performance using a boot-strapped (1000 iterations), one-sided test for significance using the pROC²⁴ package in R.²⁵

The importance of each protein biomarker included in the Disease Activity model was evaluated using mlxtend²⁶ by the mean decrease in the AUROC after permutation (repeated 1000 times) of only that protein marker in the trained Pathway and Disease Activity models, compared with the AUROC for predictions with no permutation. Proteins that displayed a larger decrease in the AUROC were more important for the Disease Activity model. NfL was identified as the most important biomarker, followed by tumor necrosis factor superfamily member 13B (TNFSF13B) (Supplementary Figure 1).

Statistical analysis was performed to characterize the algorithm at the model level (eg, prior to scoring individual samples) and at the score level. Thresholds were established for the score levels that corresponded to Disease Activity scores based on the number of Gd+ lesions. Patients with no lesions were considered to have low (1.0–4.0) disease activity, patients with 1 Gd+ lesion had moderate (4.5–7.0) disease activity, and patients with ≥2 Gd+ lesions had high (7.5–10.0) disease activity. The low versus moderate/high threshold was selected based on sensitivity. The rationale for selecting sensitivity for the low versus moderate/high threshold was that the presence of any number of Gd+ lesion(s) was an accurate and reliable indicator of active disease. The low/moderate versus high threshold was selected based on accuracy. The rationale for selecting accuracy for the low/moderate versus high threshold was that all samples within this score range were expected to have radiographic evidence of disease activity; the primary factor would therefore be optimizing the ability to distinguish a single lesion versus ≥2 Gd+ lesions.

Statistical analysis for the Train subset was performed by the Octave Bioscience Data Science Team who remained blinded to the analytical results (eg, protein concentrations) and clinical assessments in the Test subset until the algorithm was finalized. To investigate generalizability of the Disease Activity and Pathway algorithms determined using the Train subset, samples from earlier research and development studies were also evaluated prior to the Test subset analysis (analyses not reported here).

Standard protocol approvals, registration, and patient consents

The study was approved by the Mass General Brigham institutional review board (Somerville, MA, USA) and the WCG institutional review board (Puyallup, WA, USA). All patients provided written informed consent.

Data availability

Access to anonymized data not published within this article and the study protocol can be made available by request from any qualified investigator once a data-sharing agreement is in place.

Demographic and patient characteristics

A total of 614 serum samples were included from two different sources and split into Train (n=426; algorithm development) and Test (n=188; evaluation) subsets. Patient demographics and characteristics were well balanced between Train and Test subsets (Table 1).

Single-protein model evaluation

The AUROC of demographically corrected, log₁₀-transformed, LOQ-imputed individual biomarkers evaluated in the Train and Test subsets ranged from 0.436 to 0.726. NfL was the highest performing protein, with an AUROC of 0.726 for the Test subset. The biomarkers that correlated with nominal significance (p<0.05, no multiple hypothesis testing correction) with Gd+ lesion presence included NfL, cluster of differentiation 6 (CD6), C-X-C motif chemokine ligand 13 (CXCL13), interleukin 12β (IL-12β), and TNFSF13B (Supplementary Table 2).

Multi-protein model performance and optimization

Prior to incorporation of the scoring algorithm, a threshold for the multi-protein stacked classifier was chosen to provide a sensitivity of at least 0.80 for the Train subset. This threshold led to a sensitivity of 0.684 and a specificity of 0.714 in the Test subset. The AUROC of the multi-protein stacked classifier was 0.807 and 0.781 for the Train and Test subsets, respectively, based on Gd+ lesion presence (Supplementary Table 3).

Clinical validation of the multi-protein model

The multi-protein model developed on the Gd+ lesion presence endpoint using the Train subset was then applied to the Test subset. The final model performance achieved an AUROC of 0.781 relative to the Gd+ lesion presence/absence endpoint, 0.750 relative to the N/E T2 lesion presence endpoint, and 0.768 relative to the active/stable disease endpoint. In each case, the multi-protein model was found to have significantly greater (p<0.05) performance when compared with the top-performing single-protein model based on demographically corrected NfL (p<0.05).

By comparison, the single-protein model for demographically corrected NfL had an AUROC of 0.726 for Gd+ lesion presence, 0.660 for N/E T2 lesion presence, and 0.683 for active/stable disease. The multi-protein model also outperformed NfL models based on log₁₀-transformed, LOQ-imputed NfL protein concentrations with no demographic correction. The AUROCs for the single-protein, demographically uncorrected NfL models were 0.694, 0.619, and 0.645 for Gd+ lesion presence, N/E T2 lesion presence, and active versus stable disease, respectively (Figure 2).

• Download figure
• Open in new tab

Figure 2.

Summary of final performance of the multi-protein model and top performing single-protein demographically corrected (NfL) model in AUROC based on primary (A: Gd+ lesion presence) and exploratory assessments (B: N/E T2 lesion presence and C: active/stable disease status). P-values were from a boot-strapped (1000 iterations), one-sided test for significance comparing the multi-protein and single-protein models. AUROC performance and p-values are also shown for a single-protein NfL model with no demographic correction. AUROC = area under the receiver operating characteristic curve; DC = demographic correction; Gd+ = gadolinium-positive; N/E = new/enlarging; NfL = neurofilament light polypeptide chain.

Once the final model was optimized, score-level performance of the model in the Train and Test subsets was evaluated using 2 × 2 confusion matrices. Score-level precision for the Train subset is shown in Table 2. Selection of the low versus moderate/high threshold was based on sensitivity (for the comparison of no lesions vs ≥1 Gd+ lesions); selection of the low/medium versus high threshold was based on accuracy (for the comparison of 0 or 1 vs ≥2 Gd+ lesions).

For Test samples with ≥1 versus 0 Gd+ lesions, the sensitivity and NPV for the low versus moderate/high cutoff were determined to be 0.737 and 0.775, respectively. A diagnostic odds ratio demonstrated that the odds of having ≥1 Gd+ lesions among samples with a moderate/high Disease Activity score was 4.49 times the odds of having ≥1 Gd+ lesions among samples with a low Disease Activity score (Table 2). A comparison of the performance of the Disease Activity score in these samples using the multi-protein model with that of the highest performing single-protein model based on demographically corrected NfL demonstrated that the multi-protein demographically corrected model outperformed the single-protein demographically corrected NfL model (Figure 3).

• Download figure
• Open in new tab

Figure 3.

Analysis of disease activity using the A) top-performing single-protein demographically corrected (NfL) model and B) multi-protein model in patients with ≥1 Gd+ lesions. Orange bars = Test subset; Blue bars = Train subset. Gd+ = gadolinium-positive; NfL = neurofilament light polypeptide chain.

For the Test samples with either 0 or 1 Gd+ lesions when compared with those samples with ≥2 Gd+ lesions, the accuracy at the Disease Activity score to predict low/moderate versus high cutoff was determined to be 0.894. The diagnostic odds ratio demonstrated that the odds of having ≥2 Gd+ lesions among samples with a high Disease Activity score were 20.99 times the odds of having ≥2 Gd+ lesions among samples with a low/moderate Disease Activity score (Table 2).

The overall performance of the Disease Activity score and four Disease Pathway scores is shown in Figure 4. The score distribution and respective box plots with Gd+ lesions for each of the four Disease Pathway scores are shown in Supplementary Figure 2. The centering and scaling strategy for the four Disease Pathway scores resulted in sufficient correlation to the overall Disease Activity score while retaining an independent signal. Finally, a stacked bar plot of the results for the Train and Test subsets demonstrated that the calculated Disease Activity score reflected both the likelihood and severity of radiographic disease activity, based on the presence or absence of and count of Gd+ lesions. As shown, patients without Gd+ lesions had a low Disease Activity score (in blue), patients with high Gd+ lesions had a high Disease Activity score (in orange) and those with moderate disease activity had a medium range of Disease Activity score (in yellow; Figure 5).

• Download figure
• Open in new tab

Figure 4.

Correlation of the Disease Activity and Disease Pathway (A: Immunomodulation, B: Neuroinflammation, C: Myelin Biology, and D: Neuroaxonal Integrity) scores in the multi-protein model. The solid line indicates the linear regression fit between the Disease Activity and Disease Pathway scores, and the equations for the lines are given in the upper left corner of each figure. DA = Disease Activity; R² = coefficient of determination.

• Download figure
• Open in new tab

Figure 5.

Stacked bar plots of the Disease Activity score in the A) Train and B) Test subsets. Low disease activity (blue) was associated with no Gd+ lesions, moderate disease activity (yellow) was associated with ≥1 Gd+ lesions, and high disease activity (orange) was associated with ≥2 Gd+ lesions. Bars for Disease Activity scores for which there were no samples are blank. DA = Disease Activity; Gd+ = gadolinium-positive.

The distribution of individual biomarkers, and the Disease Activity and Disease Pathway scores obtained from the MSDA Test were analyzed for each of six disease-modifying therapy (DMT) categories (anti-CD20s, natalizumab, interferons, dimethyl fumarate, fingolimod, and glatiramer acetate) in 466 samples, reflecting the therapy the patient was on at the time of the blood draw. The remaining samples did not fall in one of these six categories and were categorized as either other (treated with a DMT that was not analyzed independently due to low sample count) or blank (information relating to therapy was not provided in the clinical dataset). Significant analysis of variance (ANOVA) Bonferroni corrected p-values were observed for 14 of the 18 individual biomarkers utilized in the MSDA Test algorithm (Supplementary Figure 3). Samples (n=40) from the anti-CD20 category had the lowest Disease Activity score on average (3.11 ± 1.77), followed by samples (n=129) from the natalizumab category (4.17 ± 1.63). The highest Disease Activity score on average was associated with samples (n=62) from the glatiramer acetate category (6.39 ± 1.67; Supplementary Table 4).

Classification of evidence

The MSDA Test is a multi-protein, serum-based biomarker assay designed to quantitatively measure disease activity using the protein levels of biomarkers present in the serum of patients with MS. Protein concentrations were LOQ imputed, log₁₀ transformed, and demographically adjusted for age and sex. The combination of multiple proteins was used to calculate four Disease Pathway scores and an overall Disease Activity score. The protein selection was intended to reflect the various biological pathways associated with MS pathophysiology. Using associations of the results from the multi-protein and single-protein models and the radiographic/clinical assessments, the MSDA Test was clinically validated in this study. This study provided Class II evidence to demonstrate the clinical validation of MSDA Test for disease activity assessments in MS.