Multiple sclerosis (MS) is an inflammatory disorder of the brain and spinal cord in which focal lymphocytic infiltration leads to damage of myelin and axons (Compston and Coles, 2008). MS affects approximately 2.5 million people worldwide, with an onset age commonly between 20 and 40 years, and an incidence more than twice as high in women compared to men (McAlpine and Compston, 2005).

The majority of MS patients (85%) usually experience a first attack defined as Clinically Isolated Syndrome (CIS), and will develop a relapsing-remitting (RR) form (Miller et al., 2012). Two thirds of the RR patients will develop a secondary progressive (SP) form, while the other third will follow a benign course (Scalfari et al., 2010). The rest of MS patients (15%) will start directly with a primary progressive (PP) form.

The criteria to diagnose MS forms was originally formulated by McDonald et al. (2001) and revised by Polman et al. (2005, 2011). They all rely on using conventional magnetic resonance imaging techniques (MRI) such as T1-weighted, gadolinium-enhanced T1-weighted MRI, as well as T2-weighted and FLAIR, due to a high sensitivity for visualizing MS lesions. Conventional MRI is also used for quantifying lesion load (LL), a marker of inflammation process but only a moderate predictor of MS evolution (Filippi et al., 1994).

More recently, advanced magnetic resonance techniques such as ¹H-Magnetic Resonance Spectroscopic Imaging (MRSI), Diffusion Tensor Imaging (DTI), and Magnetization Transfer Imaging (MTI) have been shown (Rovira et al., 2013) to provide a better characterization of the normal appearing white matter (NAWM) and thus a better understanding of the pathological mechanisms of MS. MTI metrics reflect the demyelination and remyelination processes and have been shown to predict the evolution of MS lesions. DTI metrics are very sensitive to the MS pathology and have been shown to be mainly affected by myelin loss and decreased neuronal integrity. MRS metrics provide high MS pathological specificity as well as high sensitivity to biochemical changes. Decrease of N-acetyl-aspartate (NAA) was observed in both chronic lesions and NAWM, reflecting a neuronal integrity loss (Rovira et al., 2013). Choline (Cho) and Creatine (Cre) contents were found to be increased in WM lesions and in NAWM, indicating the presence of severe demyelination and cell proliferation in relation with inflammatory processes (Tartaglia et al., 2002; Sajja et al., 2009).

Therefore, in this study we investigate the added value of magnetic resonance metabolic features (NAA/Cho, NAA/Cre, Cho/Cre) combined with routinely collected clinical MS data [e.g., patient age, disease duration (DD), Expanded Disability Status Scale (EDSS)] and lesion load values (LL). To this purpose, we build multiple binary classifiers to automatically discriminate between different clinical forms of MS patients, by training each classifier on combinations of clinical data, lesion loads, and metabolic features.

Figure 1 shows boxplots comparing MR metabolic features (Figures 1A–C) and lesion loads (Figure 1D) extracted from HC and each MS course. Boxplots are drawn using default style in MatLab, meaning the middle line inside the box represents the median value, the vertical limits are the 25th and 75th percentiles (q₁ and q₃), each whisker covers 1.5 the interquartile range (i.e., q₃−q₁), and the crosses outside the whiskers represent outliers. Supplementary Figures 1–4 show the MS data points in various 2-D feature spaces.

Using the previously described (Section 2.6) linear mixed-effects models we found that the fixed effect MS course is statistically significant in the evolution of NAA/Cho, NAA/Cre, Cho/Cre, and LL, with corresponding p-values of: 3.4 × 10⁻⁶, 2 × 10⁻⁴, 2.3 × 10⁻², and 2.6 × 10⁻⁴. Table 3 provides adjusted p-values for multiple comparisons between the MS groups.

Table 4 shows F1-scores after training LDA using only metabolic ratios, as clinical data and lesion loads were not available for healthy controls. Corresponding BAR, sensitivity and specificity values of this table can be found in Table A1 in Appendix. If F1-scores are missing, then the classifier assigned all data points to the negative class (second MS group).

Surprisingly, the F1-scores for separating HC from any MS course are very low, and the same holds true for separating very early MS form (CIS) and the most likely MS evolution, RR and RR+SP. In contrast, for RR vs. PP we find that all three metabolic ratios have F1-scores higher than 75, with a maximum of 78 for NAA/Cre. For RR vs. SP the F1-scores are lower, with a maximum of 69 after combining all metabolic features.

Table 5 shows F1-scores of classification tasks involving only MS patients. Training was done on seven different combinations of features to evaluate the classification power of clinical data, lesion loads, and metabolic features. Corresponding BAR, sensitivity, and specificity values can be found in Appendix in Tables A2–A4, respectively. If F1-scores are missing, then the classifier assigned all data points to the negative class (second MS group).

The highest F1-scores for CIS vs. RR and CIS vs. RR+SP, respectively 71 and 72, were achieved by SVM-rbf trained on clinical data and lesion loads. Training any classifier only on metabolic features yielded very low F1-scores.

The highest F1-score for RR vs. PP (85) was achieved by LDA using patient age, disease age, and EDSS. Adding all spectroscopic information maintained the F1-score at 85, while adding lesion load lowered the F1-score at 79. LDA outperformed SVM-rbf and RF in all RR vs. PP cases, always achieving an F1-score higher than 70.

The highest value for RR vs. SP (87) was first achieved after training SVM-rbf on clinical and metabolic features, but also with LDA trained on all features combined (clinical data, lesion loads, and metabolic features). SVM-rbf outperfomed LDA in the majority RR vs. SP cases, but only with 1–2%.

In this paper, we present results for nine binary classification problems using clinical data, lesion loads and metabolic features extracted from MS patients and healthy controls. We focused on metabolic features as numerous studies showed significant metabolic alterations in MS patients of different MS forms. It has been demonstrated that metabolic abnormalities in MS patients are not restricted to lesions alone (Husted et al., 1994; Doyle et al., 1995; Narayanan et al., 1997; Fu et al., 1998; Narayana et al., 1998; Sarchielli et al., 1999; He et al., 2005) and NAWM tissue is well known to be altered in MS (Narayana, 2005; De Stefano and Filippi, 2007). Concentrations of NAA in NAWM were shown to be significantly lower in MS patients (Bitsch et al., 1999; Suhy et al., 2000; Bjartmar et al., 2001; Inglese et al., 2003; Tiberio et al., 2006; Wattjes et al., 2007, 2008). Concentrations of Cho and Cre in NAWM were shown to be significantly higher in MS patients (Narayana et al., 1998; Tourbah et al., 1999; Suhy et al., 2000; Tartaglia et al., 2002; Inglese et al., 2003). Concentrations of NAA/Cre in NAWM were shown to be significantly lower in MS patients (Leary et al., 1999; Narayana et al., 2004). Multiple studies also report significant differences between metabolite concentrations in lesions vs. NAWM of HC: lower NAA and increased Cho and Cre (Wolinsky et al., 1990; Larsson et al., 1991; Davie et al., 1994, 1997; Narayana et al., 1998; Arnold et al., 2000; He et al., 2005).

Our findings are in agreement with these previous reports as decreased NAA and increased Cho and Cre contents were measured in NAWM and lesions of MS patients. After building linear mixed-effects models to properly analyze the statistical difference between the four clinical courses, we observed significant differences at the disease starting points of all MS courses using four response variables, namely the lesion load, NAA/Cre, NAA/Cho, and Cho/Cre ratios. A cross-sectional study (Hannoun et al., 2012) based on a subset of our MRSI data found statistical differences in the NAA/Cre and NAA/Cho ratios between HC and RR, PP, SP, and RR+PP+SP patients. To our knowledge, there is only one study that reports sensitivity and specificity values for classifying healthy controls from MS patients based on spectroscopic features. Inglese et al. (2003) show that absolute values of choline in NAWM can differentiate 9 controls and 10 out of 11 RR patients.

Other MS classification studies are Muthuraman et al. (2016) and Kocevar et al. (2016), both based on diffusion features. The first one reports a classification accuracy of 97% between 20 CIS and 33 RR patients. The second one analyzes classification tasks based on DTI data from a cross-sectional subset of our database. They found very high F1-scores (91.8% for both HC-CIS and CIS-RR) after training SVM-rbf on six global brain connectivity metrics. For RR vs. PP their maximum F1-score was 75.6%, which is lower than our results based on metabolic features, while for RR vs. SP, their maximum F1-score was 85.5%, which is comparable to our results. It is also worth mentioning that they did not use any clinical data, which might improve their results.

In this study, we analyzed the added value of combining standard clinical data with quantitative magnetic resonance features. To this purpose, we trained linear and non-linear classifiers only on advanced MR features, and then only on clinical data. Afterwards we train the classifiers on clinical data combined with lesion loads and metabolic features.

Although, MS patients are expected to have significantly different WM metabolism compared to healthy controls, this difference was not reflected in the metabolic average obtained from “good quality” voxels (Supplementary Figures 1A,B). This result is not entirely surprising, considering that we averaged over a high number of voxels, and the subtle lesion information could be lost in the average. However, we can visually see in Supplementary Figures 1C,D that the two progressive MS courses tend to have lower NAA/Cho and NAA/Cre ratios than healthy controls.

CIS and RR patients' distribution over the NAA/Cho and NAA/Cre feature space do not differ much, as seen in Supplementary Figure 2A. Disease duration interval for RR patients is much larger than for CIS patients, as most of CIS patients have a disease duration lower than 5 years, which can be seen in Supplementary Figure 3A. Because RR patients have more relapses than CIS patients, the number of lesions will be higher and the lesion volume as well, while EDSS scores are mainly in the same range, as seen in Supplementary Figure 4A. BAR values in Table A2 show a maximum of 85, when combining patient age, disease duration, EDSS, and lesion load. However, the corresponding maximum F1-score of 71 is much lower because the dataset is unbalanced (61 CIS vs. 214 RR), heavily influencing the classifier's precision. In this case the F1-score reflects better than BAR the difficulty of discriminating CIS from RR forms.

CIS and SP patients' distribution over different features is visible in Supplementary Figures 2B, 3B, 4B and it is clear that these two are the least and most advanced forms of MS. Because RR patients will eventually evolve into SP forms during their lifetime, we grouped together RR and SP patients for a separate classification task versus CIS patients. BAR values in Table A2 show a maximum of 92, when combining patient age, disease duration, EDSS and lesion load. The same discussion as for CIS vs. RR apply: the corresponding maximum F1-score is only 72 because the dataset is very unbalanced (61 CIS vs. 410 RR+SP) and the precision will be very low.

RR and PP patients can be discriminated using only EDSS by visually inspecting Supplementary Figure 4C. Training a linear classifier on clinical data (patient age, disease duration, and EDSS) gives the maximum F1-score of 85. Adding the 3 metabolic features keeps the score at 85, while adding lesion load information lowers the score to 79. This drop in the F1-score suggests that lesion load is not useful in differentiating RR from PP patients. Indeed, these two MS forms have the closest lesion load averages (16.7 and 20.8 ml), as shown in Table 1. In contrast, the clinical status of RR and PP patients are very different, as reflected by the EDSS values of 2 for RR and 4 for PP. Moreover, training LDA on individual metabolic features always provided higher F1-scores than lesion load, therefore we can conclude that for RR vs. PP, metabolic features have a higher discrimination power than LL. BAR values in Table A2 are also closer to the F1-scores in Table 5 because the dataset is more balanced compared to previous cases.

RR and SP patients can also be discriminated using only EDSS by visually inspecting Supplementary Figure 4D. Our results showed that EDSS is very important in differentiating RR patients from primary or secondary progressive patients. We also report consistent higher F1-scores for classifiers trained only on lesion load compared to classifiers trained only on metabolic features. Furthermore, it is clearly visible in Table 4 that we obtain higher F1-scores for this classification task using multiple features, compared to the rest of 8 tasks. These findings suggest that in the future it might be possible to build a decision support system using clinical data combined with lesion loads and metabolic features.

However, this study suffers from a few limitations caused by the low scanning frequency of only 1.5 Tesla. Firstly, it is known that the sensitivity of lesion load segmentation is improved by scanning at higher frequencies (Sicotte et al., 2003). Therefore, our LL values may not reflect entirely the pathological changes inside the brain. Secondly, the signal to noise ratio of MRSI is proportional to the scanning frequency, meaning our metabolites' quantification is not entirely accurate. In order to obtain true metabolites values, we would have to measure T1 and T2 relaxation times of water for each patient, which would be impossible in clinical practice. Moreover, spectroscopic signal scales can differ from patient to patient, resulting in large metabolite variations. To overcome some of these limitations, we use as features all three metabolite ratios (NAA/Cho, NAA/Cre, Cho/Cre). By doing so, we expect to retain the most valuable information.

When comparing classification tasks from a computational point of view, LDA is clearly the winner as the training period last only 3 h using a computer with 8 threads. Training both SVM-rbf and RF took around 20 days in total and it was done using 60 threads, meaning LDA is approximately 600 times faster than SVM-rbf or RF. Also, the maximum F1-scores for RR vs. PP and RR vs. SP were obtained with LDA and SVM-rbf, suggesting that a linear classifier performs equally good as a non-linear classifier in these cases.

This study is a proof of concept that investigates the added value of MR metabolites combined with clinical data and lesion loads, in classifying MS patients and healthy controls. Clinical data is routinely collected by doctors, lesion load is a known marker of neurodegeneration, while MR metabolites have been shown to provide high specificity of MS pathology. In order to better understand the underlying MS pathological mechanisms, we used three different machine learning methods, one linear and two non-linear, and had a strict quality control for extracting metabolic features. Despite all our efforts, averaging metabolite ratios over “good quality” voxels provides only moderate biomarkers for discriminating MS groups (i.e., RR vs. PP). In general, combining patient age, disease duration, EDSS, and averaged metabolic ratios, leads to the highest classification results. We believe extracting metabolic information from specific brain sub-regions of the MRSI grid (e.g., NAWM) should provide a more detailed view of MS pathology and help the classification tasks. Therefore, further investigations about the MS patients' evolution will be done in the future on sub-regions metabolite quantification, DTI-based brain connectivity metrics, patient treatment, and multi-class classification.

In this paper, we performed nine binary classification tasks and report F1-scores and BAR values after learning linear and non-linear classifiers on combinations of clinical data, lesion loads, and metabolic features. We presented a simple method to compute metabolic features by averaging metabolite ratios over “good quality” voxels of a MRSI grid. Using linear mixed-effects models we found that the MS course is statistically significant in the evolution of four response variables: Lesion Load, NAA/Cre, NAA/Cho, and Cho/Cre ratios. Our results showed that the best classifier for discriminating CIS from RR or RR+SP is SVM-rbf trained on clinical data and lesion loads. We also showed that discriminating RR from PP or SP with high accuracy is possible when training LDA on clinical data. For RR vs. PP, adding metabolic features will not change the results, while for RR vs. SP, adding metabolic features and lesion loads will slightly improve the results.

AI, GK, CS, DS, SV, and DS designed the methodological part of the study. FD and DS designed the clinical part of the study and collected the data. AI analyzed the data. All authors contributed to the manuscript.