The virtual multiple sclerosis patient

Abstract

Multiple sclerosis (MS) diagnosis typically involves assessing clinical symptoms, MRI findings, and ruling out alternative explanations. While myelin damage broadly affects conduction speeds, traditional tests focus on specific white-matter tracts, which may not reflect overall impairment accurately. In this study, we integrate diffusion tensor immaging (DTI) and magnetoencephalography (MEG) data into individualized virtual brain models to estimate conduction velocities for MS patients and controls. Using Bayesian inference, we demonstrated a causal link between empirical spectral changes and inferred slower conduction velocities in patients. Remarkably, these velocities proved superior predictors of clinical disability compared to structural damage. Our findings underscore a nuanced relationship between conduction delays and large-scale brain dynamics, suggesting that individualized velocity alterations at the whole-brain level contribute causatively to clinical outcomes in MS.

Keywords: health sciences; neuroscience.

Graphical abstract

Figure 1

Overall pipeline The patient undergoes non-invasive brain imaging (MRI, DTI). Based on these images, the subject-specific connectome is constructed. (A) The subject-specific tractography is then used to generate synthetic data by placing a mathematical model of neural activity at each region, (B) and incorporating global couplings and conduction velocities randomly drawn from physiologically plausible ranges as the prior. Subsequently, summary statistics of spectral power are extracted from the synthetic data. Probabilistic machine learning techniques are then employed to learn an invertible function (C) capable of relating parameters to spectral features, (D) along with their associated uncertainty. By providing the source-reconstructed magnetoencephalography data acquired from both MS patients and control subjects, (E and F) the spectral properties, serving as the data feature for model inversion, are calculated. (G) Then, posterior distributions for parameters (global coupling and conduction velocities) are approximated. (H) This leads to an efficient estimation of the most likely velocities given the observed summary statistics, namely spectral features. (I) Finally, based on the estimated velocities, demographic features, clinical data, (J) the total volume of lesions and global atrophy, (K) we built a multilinear model to predict individual clinical disability, (L) as quantified by the Expanded Disability Status Scale (EDSS).

Figure 2

Empirical data features Data features extracted from empirical MEG data for control and MS groups. (A) Median power spectral density plotted for controls (red) and MS patients (blue) in the 1–30 Hz range. Shaded region represents the alpha band (8–13 Hz). (B) The peak frequency shows non-significant changes in the control group relative to the patient group (p = 0.62). (C) The amplitude of peak frequency illustrates a significant decrease from the control to the patient group (p < 0.01). (D) The area under the PSD shows a significant decrease from the control to the patient group (p = 0.02).

Figure 3

Model output (A‒C) Heatmaps indicate the median frequency, the median peak alpha PSD and the total alpha power, respectively, for the range $G\in \left[800,1800\right]$ and $V\in \left[1,20\right]$. (D) Simulation output represented as the median spectrogram for $G$ = 1200, $V$ = 16 m/s (indicated as ∗ in heatmaps) and (E) The corresponding median PSD (blue), for a given simulation run showing a peak in the alpha band (smoothed average in red). All PSDs reported in log-scale.

Figure 4

Estimated parameters for individuals Estimated posterior distribution for two brain parameters from the PSD of the empirical MEG data, pooled over control and MS patient groups. (A) The global coupling strength $G$ shows non-significant changes (p = 0.87). (B) The velocity parameter $V$ significantly decreases in the control group relative to the patient group (p < 0.01).

Figure 5

Diagnostics of the inference process for aggregated across individuals (A and B) Observed (red) and predicted (blue) PSD of MEG data and corresponding time-series median average over brain regions, respectively. (C and D) Inferred posterior distributions for the global coupling strength $G$, and the velocity $V$, respectively, given PSD features (amplitude, median frequency and total power). Increasing the number of simulations for training steps yields progressively tighter posteriors and, thus, a more accurate estimate. (E) Joint posterior distribution between parameters $G$ and $V$ estimated from 20000 simulations (correlation = 0.75). The ground truth parameters are shown in red, the high-probability parameters in yellow, and the low-probability ones in blue. (F) Sensitivity analysis using the estimated posterior, indicating stronger model sensitivity to $V$ than to $G$ (the Eigenvalues for $G$ and $V$ are 1.2e-05 and 8.2e-05, respectively).

Figure 6

Prediction of clinical outcome Multilinear regression model with leave-one-out cross-validation (LOOCV) performed to test the capacity of the estimated speed to predict the EDSS scores in MS patients. (A and D) Variance explained by the additive model including five variables (i.e., gender, age, disease duration, lesion load, and estimated speed). Adding the estimated conduction velocities significantly increased the predictive power in both classical multilinear (panel A, p = 0.028) and cross-validated (panel D, p = 0.0417) models. Data are represented as mean ± SD. (B and E) Empirical and predicted EDSS scores. (C and F) standardized residuals of the model.