In total, out of 141 MS patients that were retrospectively screened for this project, 111 MS patients that underwent a comprehensive and detailed clinical assessment were finally included in the analysis ( Figure 1 ). The selected cohort included MS patients with MRI (T2-hyperintense LV and GMV), blood (sNfL and sGFAP), and OCT (RNFL and GCIPL) measurements at the outpatient clinic of the Department of Neurology, at the University Medical Center Mainz (Germany) ( Table 1 ). All included patients had relapsing-remitting multiple sclerosis (RRMS) as diagnosed according to the 2017 revised McDonald diagnostic criteria (15). The mean (± standard deviation) disease duration of all patients at study inclusion was 3.15 ± 4.26 years. All diagnostic baseline measurements were performed within 6 months of study inclusion. An experienced neurologist clinically assessed patients and their Expanded Disability Status Scale (EDSS) score at study entry and follow up visit (3.74 ± 1.25 years), along with clinical relapse history over the study period. EDSS progression was defined as an increase of ≥ 1 point in the EDSS score for a baseline score of ≥ 1.5 or a 1.5 point increase for a baseline score of 0 (16). A clinical relapse was defined as a monophasic clinical episode with new neurological symptoms, lasting more than 24 h and in the absence of fever or infection (15). The annualized relapse rates (ARR) were calculated by dividing the total number of all observed relapses by the total number of patient-years. All measurements were performed at least 30 days after a clinical relapse and/or a high-dose corticosteroid treatment.

Serum samples were collected by attending physicians at the University Medical Center Mainz. Samples were processed at room temperature within 2 hours. Serum samples were spun at 2000xg at room temperature for 10 minutes, aliquoted in polypropylene tubes and stored at −80°C. sNfL and sGFAP concentrations were measured as previously described (10, 14). In brief, sNfL and sGFAP levels were determined using the highly sensitive single molecule array (SiMoA) technology (17). Samples were measured in duplicates by SiMoA HD-1 (Quanterix, USA) using NF-Light Advantage kits according to the manufacturer’s instructions. The mean inter-assay and intra-assay coefficient of variation was less than 10%. Measurements were performed in a blinded fashion without information about clinical data.

MRI data acquisition was performed as previously described (14). In brief, structural MRI was performed on a 3-Tesla MRI scanner (Magnetom Tim Trio, Siemens, Germany) with a 32-channel receive-only head coil. In all patients, imaging was performed using a sagittal 3D T1-weighted magnetization-prepared rapid gradient echo (MP-RAGE) sequence (TE/TI/TR = 2.52/900/1900 ms, flip angle = 9°, field of view = 256 × 256 mm2, matrix size = 256 × 256, slab thickness = 192 mm, voxel size = 1 × 1 × 1 mm³) and a sagittal 3D T2-weighted fluid-attenuated inversion recovery (FLAIR) sequence (TE/TI/TR = 388/1800/5000 ms, echo-train length = 848, field of view = 256 × 256 mm², matrix size = 256 × 256, slab thickness = 192 mm, voxel size = 1 × 1 × 1 mm³). A clinician scientist blinded to the patient data excluded major anatomical abnormalities based on the subject’s T1-weighted and FLAIR images of the whole brain.

The quantification of WM (white matter) volume, lesion volume and GMV was performed as previously described (14). Using voxel-based morphometry (VBM) analysis in the Statistical Parametric Mapping (SPM8) software, the GM and WM volumes were calculated. The volumes of WM lesions were assessed using the cross-sectional lesion growth algorithm of the lesion segmentation toolbox (18) included in the SPM8 software. 3D FLAIR images were co-registered to 3D T1-weighted images and bias corrected. After partial volume estimation, lesion segmentation was performed with 20 different initial threshold values for the lesion growth algorithm (18). By comparing manually and automatically estimated lesion maps, the optimal threshold (ĸ value, dependent on image contrast) was determined, and average values were calculated for each patient. A uniform ĸ value of 0.1 was applied in all patients in order to automatically estimate lesion volume and filling of 3D T1-weighted images. Subsequently, the filled 3D T1-weighted images and the native 3D T1-weighted images were segmented into GM, WM, and CSF and then normalized to the Montreal Neurological Institute (MNI) space. The quality of the segmentations was visually inspected to increase reliability.

The analysis was performed as previously described (19, 20). In brief, the Advised Protocol for OCT Study Terminology and Elements (APOSTEL) recommendations were followed (21) including a quality control for the raw OCT scans complying with the OSCAR-IB criteria (22). MS patients with accompanying diseases potentially affecting the optic nerve or other ocular disease were excluded in advance. Hence, none of the patients had a history of glaucoma, retinopathy or other neurological disorders (besides RRMS). An experienced operator performed OCT image acquisition following a unified standard acquisition protocol using a spectral domain OCT (Heidelberg Spectralis, Heidelberg Engineering, Germany) with Heidelberg Eye Explorer software (HEYEX, version 1.10.2.0). The measurements were acquired in a shaded room at ambient light without pupillary dilation. Intra-retinal layers of the macula were gauged by a standardized scan comprising 61 vertical or horizontal B-scans while focusing on the fovea at a scanning angle of 30° × 25° and a resolution of 768 × 496 pixels. Automatic real time was set to nine at high-speed scanning mode. Confocal scanning laser ophthalmoscopy was performed in parallel and revealed no evidence of pathology. No further fundoscopic imaging was carried out. To account for inter-eye within-patient dependencies, we calculated the mean of both eyes in patients with no history of optic neuritis; in patients with a history of unilateral optic neuritis, we only used the OCT scan of the non-affected eye. Hence, the main statistical analysis was performed at a per-patient level. All B-scans were automatically segmented (followed by manual correction by a trained rater) using segmentation beta-software (Spectralis Viewing Module version 6.9.5.0) of the Heidelberg Eye Explorer (version 1.10.2.0) provided by the manufacturer. The segmentation lines were the following retinal layers: RNFL, GCIPL, inner nuclear layer, outer plexiform layer and outer nuclear layer. The mean volume of the individual retinal layers was computed in an area of a radius of 3.45 mm around the fovea including the fovea using the Early Treatment of Diabetic Retinopathy Study (ETDRS) grid. Lastly, RNFL and GCIPL were finally selected as primary estimate for neuroaxonal damage of the retina, as both have been associated with brain atrophy and disability worsening (23, 24).

Statistical analysis was performed using SPSS 23 (SPSS, Chicago, IL, USA), MedCalc (Version 20.115) and GraphPad Prism 9 software. Summary statistics are presented as mean ± standard deviation (SD), or median (25^th and 75^th percentile), or number (percentage), where applicable. To create a combined variable for each biomarker combination, a binary logistic regression model for each combination (corrected for sex, age, disease duration and disease-modifying treatment) was estimated in order to get the predicted probability from each model. Then, we used this probability as the test variable in the subsequent receiver operating characteristic (ROC) procedure (14).

A ROC analysis was performed to calculate the predictive discriminating values for each biomarker and the combinations. This statistical method is preferentially used to make a series of discriminations into two different states based on a specific diagnostic variable. Here, the presence or absence of relapses or EDSS worsening, served as binary classifiers. Every value of that discriminating variable is used as a cut-off with calculation of the corresponding sensitivity and specificity.

The analysis was performed as previously described (25) using the SEM toolbox for MATLAB (version 13a; Mathworks, Natick, MA, USA). SEM represents a statistical technique that is used to test and estimate structural relationships between variables in a model. By structural, we mean that we incorporate causal assumptions as part of the model. Hence, SEM represents a multivariate technique that is able test complex relationships among multiple variables simultaneously, and estimate the strength and direction of these relationships. In our model, we explored the association between multimodal biomarkers and the clinical outcomes (clinical relapses and EDSS progression). We used the Maximum likelihood method of estimation to fit the models. In order to adjust the models for a large sample size, we used the Root Mean Square Error of Approximation (RMSEA) index, which improves precision without increasing bias (26). The RMSEA index estimates lack of fit in a model compared to a perfect model and therefore should be low. In all models, the Invariant under a Constant Scaling (ICS) and ICS factor (ICSF) criteria should be close to zero, indicating that models were appropriate for analysis. Finally, based on the Akaike Information Criterion (AIC) the quality of each model relative to other models was estimated, with smaller values signifying a better fit of the model. The strength of associations between the variables in the models was quantified by standardized coefficients (s), ranging from 0 (no association) to 1 (very strong association). To correct for potential confounders the models were adjusted for sex, age, disease duration and disease-modifying treatment (DMT). P-values less than 0.05 were considered statistically significant.

All demographics and clinical characteristics of the investigated cohort are summarized in Table 1 . In total, 141 early MS patients with baseline MRI and OCT were selected. Thirty patients were excluded from the final analysis because either there was no serum sample available or they were lost to clinical follow-up ( Figure 1 ). The mean follow-up time in our longitudinal cohort of 111 patients was 3.74 ± 1.25 years. The mean age ± SD was 34.8 ± 9.67 years; 79 patients (71.0%) were female and 32 (29.0%) were male. The mean disease duration at study inclusion was 3.15 ± 4.26 years. All patients had a relapsing-remitting disease course (RRMS) according to the 2017 revised McDonald criteria (15). At the time of inclusion, 18 patients (16%) were not receiving any DMT, 69 (62%) were receiving a mild to moderate efficacy DMT, and 24 (22%) were receiving a high efficacy DMT. The median baseline disability, quantified with EDSS, was 1.0 (25th and 75th percentile: 0.0−2.0). Overall, 46 patients (41.4%) experienced EDSS progression during the observation period. The mean ARR was 0.21 ± 0.33; 33 (30%) patients had a history of optic neuritis. The results from blood biomarker, MRI, and OCT measurements are also summarized in Table 1 .

An overall ROC analysis was performed to determine the predictive discriminating value of the individual and combined measures to distinguish MS patients with and without disease activity (determined through the presence or absence of relapses during this time) and with and without disability progression (determined through the presence or absence of EDSS worsening over four years). Resulting values with AUC, standard error, 95% confidence interval and p-values are presented in detail in Figures 2A and 3A .

In general, none of the individual biomarkers were able to predict the occurrence of clinical relapses within the 4-year observation period (AUC-range: 0.523 – 0.602). All p-values for testing AUC = 0.5 vs. AUC ≠ 0.5 were greater than 0.05 and were hence not significantly different from a random classifier ( Figure 2B ). Only LV showed a trend towards significance (AUC = 0.602; p = 0.060). In the ROC analysis based on the presence or absence of EDSS progression, GMV was the only single biomarker to show significant predictive capability for EDSS progression on its own (AUC = 0.614, SE = 0.054; p = 0.035), whereas all other single biomarkers did not (AUC-range = 0.502 - 0.596) ( Figure 3B ).

When we combined biomarkers within their respective modality, MRI markers (LV + GMV) were able to predict both relapses (AUC = 0.631, SE = 0.054; p = 0.015) and EDSS progression over the four-year period (AUC = 0.621, SE = 0.055; p = 0.026). Combined blood biomarkers (sNfL + sGFAP) were only able to predict EDSS progression (AUC = 0.632, SE = 0.059; p = 0.025), while combined OCT measures (RNFL + GCIPL) were unable to predict either clinical relapses (AUC = 0.599, SE = 0.054; p = 0.069) or EDSS progression (AUC = 0.507, SE = 0.058, p = 0.906) ( Figures 2C , 3C ).

However, all combinations of two biomarker modalities significantly predicted clinical relapses (AUC range = 0.636 – 0.643) and EDSS progression (AUC range = 0.631 – 0.699) ( Figures 2D , 3D ). The best prediction for EDSS progression using two modalities was achieved with a combination of MRI and blood biomarkers (AUC = 0.699, SE = 0.055; p < 0.001).

Most notably, the combination of all six biomarkers achieved the highest AUC for discriminating MS patients with clinical relapse activity from those without (AUC = 0.678, SE = 0.057; p = 0.002) and for discriminating progressive from non-progressive MS patients (AUC = 0.706, SE = 0.055; p < 0.001) ( Figures 2E , 3E ). Overall, these results demonstrate that the predictive capability of single biomarkers remains limited except for GMV, whereas combining multimodal biomarkers stepwise improves their accuracy in prediction of both relapse activity and disease progression within early multiple sclerosis.

In order to create a prediction model analyzing complex relationships among multiple variables, we next applied SEM to assess the causal relationship of the most promising biomarker combinations determined in the ROC approach, namely MRI (LV + GMV) and blood (sNfL + sGFAP) biomarkers. In addition to the ROC analysis, SEM allows us to test a model for its compatibility with the data in its entirety simultaneously. In the predictive modeling approach, the RMSEA index for the models was below 0.03 and the AIC comparing the models varied between 0.006 and 0.019. The obtained fit indices in the SEM analysis implied a good fit of the constructed models to the observed data, providing robust relations between the variables. Within the SEM model quantifying the pathways, the input variables (GMV, sNfL, sGFAP and LV) predicted both ARR and EDSS progression. Our model with resultant standardized coefficients (s) identified that GMV (s = 0.58; p < 0.01) and sNfL (s = 0.63; p < 0.01) significantly predict ARR and EDSS progression through lesion volume as mediator (ARR [s = 0.59; p < 0.01] and EDSS [s = 0.73; p < 0.001]) ( Figure 4 ). Taken together, LV mediates the path between GMV and sNfL on the one side, and ARR and EDSS progression on the other side.