This single-center observational study integrated three data modalities, more specifically clinical, MRI, and motor EP data, collected retrospectively at Noorderhart, Rehabilitation and MS Centre in Pelt, Belgium. To construct a combined dataset, all visits from the separate modalities were matched for each patient using unique patient identifiers. The final dataset included data collected during routine clinical care from June 2011 until May 2017. Clinical data were extracted from the MS data entry portal iMed (iMed, ^© 2022 MSBase Foundation, Australia) and included the following features: age, gender, and MS course. We adhered to the Transparent Reporting of a multivariable prediction model for Individual Prognosis Or Diagnosis (TRIPOD) (44) guidelines for reporting our study (see Appendix 10). This study has been approved by the Medical Ethics Committee of Hasselt University (CME2019/046). Due to the retrospective study design and the pseudonymization of data, obtaining patient informed consent was not required. All our preprocessing and model code are available at https://github.com/UHasselt-BiomedicalDataSciences/Radiomics_Epomics.

Pseudonymized T2-weighted FLAIR MRI scans were acquired using a Philips Achieva 1.5T scanner with three distinct protocols (Protocols A, B, and C), depending on the acquisition date. Further details regarding the acquisition protocols are summarized in Appendix 1. A standardized pre-processing pipeline was established to prepare the MRI data for downstream analysis. First, all MRI images were denoised using adaptive non-local means filtering and corrected for bias field inhomogeneities with the N4 algorithm. Given that Protocols A and B consist of multiple low-resolution FLAIR images per session, a super resolution reconstruction (SRR) approach called Perceptual Super-Resolution in Multiple Sclerosis (45) was applied to enhance the through-plane resolution of multi-slice structural MRIs containing MS lesions. This technique harmonizes the spatial resolution across scans, which is crucial for downstream radiomics analysis and segmentation tasks (45).

To achieve whole-brain segmentation, we applied Sequence Adaptive Multimodal SEGmentation (SAM-SEG) to all FLAIR protocols (46). SAMSEG segmented 41 anatomical brain structures (see Appendix 6) and is adaptive to different MRI contrasts and scanner types. Among the segmented structures were normal-appearing white matter (NAWM), grey matter (GM), thalamus, and cerebrospinal fluid (CSF). Additionally, white matter lesions (WML) were segmented using the lesion prediction algorithm (47) as implemented in the lesion segmentation tool (LST) toolbox¹, version 1.2.3 for SPM8². Beyond segmentation, anatomical features and estimated lesion volumes obtained from SAMSEG and LST were normalized by the intracranial volume for comparability. Furthermore, intensity normalization was performed using adaptive histogram matching on all skull-stripped FLAIR images, as previously described by Khan et al. (24).

Adding to the anatomical features and lesion volumes, radiomics features were extracted using PyRa-diomics 2.20 from multiple brain structures, including WML, NAWM, GM, CSF, and thalamus. The extracted radiomics features belonged to six feature classes: shape-based features (48), first-order statistics (FO), grey level co-occurrence matrix (GLCM) (49), grey level run length matrix (GLRLM) (50), grey level size zone matrix (GLSZM) (51), and grey level dependence matrix (GLDM) (52). Shape-based features capture the geometric properties of a region of interest (ROI), FO detail the distribution of intensities, while the remaining four feature types (GLCM, GLRLM, GLSZM, and GLDM) are collectively referred to as textural features. Grey level features were computed by discretizing the images into 50 intensity bins, following the recommendations of the Image Biomarker Standardization Initiative (53). A total of 40 anatomical features, including lesion volume, along with 14 shape-based, 18 FO, and 68 textural features per ROI, were used as input to our ML pipeline. The details of the number of radiomics features extracted per class are provided in Appendix 2.

We used motor EP data of which an anonymized version was previously published by Yperman et al. (28). In this dataset, MEPs were recorded bilaterally from the abductor pollicis brevis (APB) and abductor hallucis (AH) muscles following TMS of the corresponding motor cortex areas using two acquisition systems (devices A and B). Stimulation was delivered using a Magstim 2002 or Bistim device (The Magstim Company Ltd., Whitland, UK) with a 9 cm round coil at maximal stimulator output (2.2 T). Signals were recorded for 100 ms post-stimulation and digitized at either 20 kHz (device A) or 19.2 kHz (device B), with all 20 kHz recordings down-sampled to 19.2 kHz for consistency. Device A applied a 0.6 Hz – 10 kHz band-pass filter, whereas Device B employed a 100 Hz high-pass filter for noise reduction.

Electromyographic activity was recorded using three surface electrodes per limb: for upper limbs, electrodes were placed on the APB muscle, the proximal phalanx of the thumb, and the dorsum of the hand (ground); for lower limbs, on the AH muscle, the big toe, and the dorsum of the foot (ground). Stimulation intensity started at 45% (upper limbs) and 50% (lower limbs) of maximal stimulator output, increasing in 5% increments until MEP amplitude reached 1 mV or no further increase was observed. Trials affected by artefacts or poor quality were excluded. MEPs were recorded across multiple excitation levels, producing a series of EPTS values. Following clinical expert recommendations, only the EPTS corresponding to the maximal PPA was retained for analysis as the most representative and reliable for predictive modelling.

As model input, we incorporated the key motor EPTS features identified by Yperman et al. (37). Specifically, the original study employed a time series feature extraction by using highly comparative time series analysis (HCTSA) (54), followed by a multi-step feature selection pipeline that included mutual information-based filtering, hierarchical clustering to remove redundant features, and the Boruta algorithm (55) to identify the most relevant predictors for MS disability worsening. The most relevant feature 201 for the APB muscle (i.e., EPTS Sliding Window Feature) was found to characterize how quickly the time series returns to an average baseline following an initial peak. For the AH muscle, an autoregressive modelling-based feature (i.e., EPTS Mean Absolute Lag-1 Autocorrelation of Prediction Errors) demonstrated strong predictive value. For this feature, a high value means the EP still carries a clear, repeating pattern that the simple autoregressive model could not remove, whereas a low value means that, after accounting for the obvious trends, the signal exhibits a behavior similar to a random background noise. For comparability with prior literature, we also incorporated established electrophysiological biomarkers such as latency and PPA.

The prediction outcome of disability worsening was defined following established criteria (56) based on changes in the EDSS between two different time points: the baseline measurement t0, corresponding to the MRI acquisition, and the closest EDSS evaluation two to three years later t2y. Worsening was defined using the thresholds outlined in Equation 1. Consequently, each participant was assigned a binary outcome: worsened or stable (non-worsened). A single follow-up EDSS evaluation at t2y was considered sufficient for determining disability worsening, rather than requiring a subsequent confirmatory assessment. This choice increased the number of usable training instances.

Patients were included in the final cohort if they had: (1) an MRI scan at baseline (t₀), (2) a motor EP and EDSS assessment performed within six months before to three months after t₀, and (3) an EDSS evaluation between two and three years after the baseline (t_2y). If multiple EDSS or EP assessments were available within the baseline window, the measurement closest to the MRI date was selected. Likewise, when several EDSS assessments were available beyond two years, the one closest to the two-year mark was used for t2y. A complete case analysis was performed, excluding episodes with missing data on any required modality. Figure 1 illustrates the inclusion criteria for the final cohort. A two-year prediction horizon was selected following the approach of De Brouwer et al. (20).

The clinical data export from iMed (iMed, ^©2022 MSBase Foundation, Australia) initially included 1,025 unique patients and 20,456 clinical visits. Following the exclusion of visits lacking an EDSS follow-up at least two years later, 821 patients with a total of 14,940 visits remained. Within this subset, 332 patients had available EPTS data (1,741 visits), and 144 patients had corresponding MRI data (605 visits). As some individuals had multiple assessments meeting the inclusion criteria, their longitudinal records contributed several, occasionally overlapping, clinical episodes. After applying the full inclusion criteria, the final cohort comprised 127 unique patients contributing 424 eligible episodes. The distribution of the number of episodes contributed per patient is shown in Appendix 9. All statistical and ML analyses were conducted using the episode-level data. A flowchart detailing the inclusion process is shown in Figure 2.

In summary, for the predictive modelling, the input consisted of baseline clinical data (EDSS, age, and gender), MRI-derived anatomical and radiomics features, and motor EPTS features, while a binary outcome indicated whether disability worsening had occurred two to three years after the baseline.

At the patient level, individuals in the final cohort had a mean age of 42.3 years (± 11.9 SD) and a mean EDSS score of 2.5 (± 1.6 SD) at baseline (t₀). The cohort was predominantly female (75.6%), with the majority diagnosed with relapsing-remitting MS (RRMS, n = 105), followed by secondary-progressive MS (SPMS, n = 17), primary-progressive MS (PPMS, n = 4), and clinically isolated syndrome (CIS, n = 3). A total of 16 patients (12.6%) were classified as worsened, meaning they experienced EDSS-based disability worsening in at least one episode during follow-up. The cohort represented a heavily treated population, with 76.4% receiving moderate to high DMTs. All patients met the McDonald criteria applicable at the time of diagnosis (57). Summary statistics at the episode level are presented in Table 1. Overall, 5.2% of the total episodes showed EDSS worsening. The remaining values reflect the underlying patient-level characteristics.

For the statistical analysis, two-tailed Mann-Whitney U-tests were used to assess differences in all features between the worsening and stable patients prior to feature selection. This non-parametric test was selected because it does not assume normality of the data. To account for multiple comparisons, the Benjamini–Hochberg procedure was applied to control the false discovery rate, and the Bonferroni correction was used to reduce the risk of type I errors. Statistical significance was defined as p< α = 0.05. All analyses were performed using Python 3.10.12, pandas 2.2.2 (58) and NumPy 1.26.4 (59). Statistical analyses were carried out using pingouin 0.5.4 (60) and SciPy 1.14.0 (61).

To assess the viability and potential benefits of integrating motor EPTS and MRI features, an extensive ML analysis was conducted. The objective was to evaluate the predictive performance of both modalities, MRI radiomics and motor EPTS, separately, as well as in combination, with or without the inclusion of clinical variables. The entire ML pipeline is illustrated in Figure 3. A total of 96 different model setups were explored, systematically varying five key parameters:

A key aspect of our methodology was the use of a repeated stratified patient 3-fold cross-validation scheme, which ensures that each patient’s data is exclusively used in either training or testing, thereby minimizing overfitting. Given the low prevalence of disability worsening in the dataset (22 worsening episodes; 5.2%), a single fixed hold-out test split would have yielded only a very limited number of positive class instances in the test set, resulting in unstable and high-variance performance estimates. Accordingly, three folds were constructed with approximately equal numbers of worsening episodes (stratified), while assigning all episodes from a given patient to a single fold to prevent within-patient information leakage. We then repeat this process 20 times, analogous to the framework proposed by Yperman et al. (37). This ensures there are enough worsening and stable patients in each test set to obtain valid metrics. The repetitions ensure that different divisions of the dataset are explored. In the following sections, we provide a detailed breakdown of our ML framework within this cross-validation scheme. A visual overview is shown in Figure 3.

To identify features that differed significantly between stable and worsening groups on an episode level, we performed a Mann-Whitney U test, applying Bonferroni correction. After correction, only one feature remained statistically significant, specifically GLSZM Grey Level Non-Uniformity of the NAWM (p-value = 0.035), as shown with a red outline in Figure 4. This feature quantifies signal heterogeneity in the NAWM, with higher values observed in patients who experienced disability worsening. If, instead, the less strict Benjamini-Hochberg correction was used, 23 features were identified as significant. Boxplots for these 23 features are also displayed in Figure 4. These included 8 textural radiomics features (e.g., GLRM Run Length Non-Uniformity of the WML), 10 shape-based or anatomical features (e.g., Lesion Volume, Surface Area of the WML), 2 First-order statistics features (e.g., Total Energy of the WML), and three EPTS-derived features (e.g., EPTS Latency (APB, Right)). Most of these features showed significantly higher median values in the worsening group.

To statistically assess whether the combination of different data modalities resulted in a significant performance advantage over multiple metrics, we employed a critical difference diagram (Figure 5). Statistical significance was first determined using a Friedman test (p = 1.9 · 10⁻⁹), followed by a post-hoc Wilcoxon test with Holm correction to assess pairwise differences between models (67). The figure shows that the combined feature sets perform better with statistical significance, compared to the modalities alone. No significance can be determined when clinical features are present together with the combined feature set, nor are there any significant differences between the modalities alone.

Figure 6 presents a comparative analysis of model performance across various evaluation metrics (detailed in Appendix 4). The boxplots illustrate the distribution of performance scores across cross-validation folds, with slight improvements observed in models incorporating both motor EPTS and MRI data. Each modality and combination plotted represents the data with the inclusion and exclusion of clinical data, as no significant differences were observed in Figure 5. These models exhibited statistically significant superior median performance compared to the other modalities alone in AP, specificity, and Brier score. Additionally, the combination performs better with statistical significance than MRI radiomics features alone when considering AUROC.

Since there are many trade-offs concerning the metrics reported by the different models, the Brier score is chosen to determine the best-performing model. The model with the lowest Brier score (0.062) is an LGBM model (hyperparameters in Appendix 7). Figure 8 presents the SHAP summary plot. Appendix 8 contains the numerical SHAP values corresponding to the features in the plot. The majority of high-impact features in the model were derived from MRI radiomics. These included textural features such as GLSZM Grey Level Non Uniformity, which reflects heterogeneity within NAWM and lesion areas. Higher values were associated with an increased likelihood of disability worsening, suggesting that greater heterogeneity in brain tissue is linked to disease progression. Shape-based MRI features also played a significant role, particularly those related to brain atrophy, such as the Major and Least Axis Length of the NAWM. In these cases, lower values were associated with a higher probability of worsening, indicating that reduced structural integrity in non-lesional tissue contributes to the model’s predictions.

In addition to MRI-based features, several EPTS features emerged as important predictors. Notably, Sliding Window–based measures ranked highly in the SHAP analysis, with increased values generally corresponding to a greater risk of disability worsening. Furthermore, the PPA of EPTS signals was selected as a relevant feature, where lower values were linked to a higher probability of disability worsening.

Among clinical features, gender was the only variable retained by the model. SHAP values indicated that male gender was associated with a lower predicted risk of disability worsening compared to female gender.

In our cross-validation framework, all configurations and metrics evaluated demonstrated that combining high-dimensional MRI and EPTS features outperformed models that relied solely on one modality. This finding confirms the additional value of multimodal approaches. Structural and functional biomarkers from distinct data modalities, imaging and time series, appear to provide complementary information about disease processes that may not be captured by any single assessment. Moreover, we utilized high dimensional features instead of conventional measures, revealing subtle structural changes in seemingly normal tissue and functional disturbances that might otherwise go undetected (24, 37).

Consistent with previous research on brain texture analysis in MS (70–73), texture-based features of WMLs, most prominently GLSZM Gray Level Non-Uniformity and GLDM Dependence Non-Uniformity, were found to be associated with subsequent worsening of disability (74). The former feature quantifies the dispersion of signal intensities within a lesion, while the latter captures the variability in the spatial grouping of similar intensities. Higher values for either feature indicate greater heterogeneity and structural irregularity.

Similarly, increased textural complexity in the NAWM, exemplified by GLSZM Grey Level Non-Uniformity, also predicted EDSS worsening in both the statistical and SHAP interpretability analyses (73, 74). Together, these findings support the idea that diffuse microstructural abnormalities in both WML and NAWM reflect underlying processes such as demyelination, inflammation, and axonal loss (75, 76). In line with this, longitudinal quantitative MRI data have shown that reductions in myelin and axon volume fractions in the NAWM precede the development of more destructive lesions, such as T1-weighted hypointense black holes, further highlighting the prognostic relevance of early microstructural disruption (77).

Shape-based features of the NAWM, specifically the Least Axis Length and Major Axis Length, also appeared in SHAP analyses. These features measure the smallest and largest axes, respectively, of a three-dimensional ellipsoid that encompasses the ROI in the brain (78). Higher values in these shape measures may indicate swelling or more extensive structural alterations in the NAWM. This suggests that atrophy or morphological distortions in non-lesional tissue can have significant prognostic implications (79).

Furthermore, certain anatomical features, such as atrophy in the cerebellum or putamen, have been identified as predictors of worsening disability, supporting the established correlation between damage to subcortical and cerebellar structures and declines in motor and cognitive functions (80–83). In contrast, other anatomical features related to atrophy, such as the volumes of white matter, cerebrospinal fluid, and lateral ventricles, did not effectively distinguish between patients with stable disease and those experiencing disability worsening, which is inconsistent with earlier research (84). This inconsistency may arise from patients exhibiting significant atrophy being categorized as stable due to their already elevated baseline EDSS scores. Consequently, the two-year follow-up period may not have been long enough to capture meaningful EDSS worsening in this subgroup (see Appendix 3). Additionally, factors like variability in imaging between centers, reliance on reconstructed FLAIR MRI, and class imbalance may also contribute to the differences observed between our findings and previous studies.

Most of the features selected by our ML model were derived from MRI data, likely due to the larger number of radiomics-based inputs. However, certain EPTS features emerged among the top-ranking selections. Notably, the Sliding Window-type features, which were previously identified as highly relevant from a comprehensive pool of 7,700 motor EPTS features by Yperman et al. (37), were confirmed again in our study. Yperman et al. explain that this feature is calculated by moving a window equal to half the length of the EPTS across the signal in 25% increments, resulting in three overlapping segments. For each window, the mean is computed, and the SD of these means is divided by the SD of the full EPTS. Essentially, this metric quantifies the variability of the temporal signal after the initial stimulation peak, providing an indirect measure of how quickly and consistently the signal returns to baseline. While its precise physiological interpretation remains unclear, the repeated selection of this feature across different modalities suggests that it carries meaningful prognostic information. This insight may be particularly relevant in the context of a well-treated and clinically stable cohort, where traditional clinical and imaging markers might lack sensitivity to detect subtle progression. In such cases, the prognostic value of EPTS features could lie in their ability to reflect subclinical activity that is otherwise undetectable. Further research is needed to explore the underlying physiological significance of this feature.

In addition to the Sliding Window-type features, our model also identified PPA in the APB muscle as a key predictor, which contrasts with the findings from Yperman et al. (37), where incorporating amplitude-based features did not consistently enhance model performance. One possible explanation for this discrepancy is that, in our multimodal dataset, PPA gained importance through synergistic interactions with a wide range of MRI radiomics features, thereby increasing its overall predictive contribution.

In our SHAP analysis, gender emerged as the only relevant clinical feature. Notably, being male was associated with a lower predicted risk of worsening disability. This finding contradicts the common belief that men typically experience worse outcomes in MS (85, 86). One possible explanation for this result is that men constituted a much smaller percentage of our cohort, reflecting the typical female-to-male ratio seen in MS (87, 88). Also, fewer men experienced EDSS worsening or exhibited MRI or EPTS abnormalities detected by our model (6 men compared to 16 women). As a result, the model may have learned that male gender was less correlated with the subtle functional and structural deficits associated with MS progression. Another contributing factor may be the use of synthetic oversampling. Given the small number of worsening episodes, ADASYN-generated samples may outnumber real positive cases, increasing the risk that the models might be biased. Moreover, the observed effect may reflect center-specific clinical decision-making, whereby male patients might have received more intensive treatment based on the well-known prognostic factors, potentially leading to a lower observed rate of progression in this group. Thus, the apparent protective effect of being male in our analysis could be due to cohort-specific biases rather than a genuine difference in disease progression between sexes.

While the exact physiological mechanisms underlying high-dimensional motor EPTS features are not yet fully understood, our findings demonstrate that these functional measures carry substantial prognostic information. EPTS features consistently ranked among the most informative predictors, even when evaluated alongside a large set of MRI radiomics features, indicating that electrophysiological markers capture disease-related information that is distinct from, and complementary to, structural imaging. This complementarity is particularly relevant in the context of progression independent of relapse and MRI activity (PIR(M)A) (89–91), where patients may continue to accumulate disability despite stable conventional MRI and absence of relapses, posing the increasingly important clinical question of how to identify such patients at an earlier stage.

Interestingly, while the model combining MRI and EPTS data generally achieved the best overall balance across most performance metrics, EPTS-only models demonstrated the highest AUROC and sensitivity. This suggests that EPTS alone may be more effective at identifying those few individuals who are likely to experience disability worsening, even in an imbalanced dataset where predicting the minority class was essential. High sensitivity is particularly important in a prognostic context, as it helps minimize the chances of giving false reassurance to patients who are actually at risk. Since we did not analyze the temporal sequence of changes across modalities in this study, future longitudinal research should investigate whether alterations in EPTS occur before structural imaging abnormalities. At the same time, MRI radiomics features remained essential in our best-performing models, underscoring the established role of MRI imaging in the diagnosis and monitoring of MS.

From a clinical perspective, incorporating MEP-based measures into routine evaluations for MS could be beneficial, given their relatively low cost and general feasibility. Brief MEP sessions impose a minimal additional burden and may uncover functional impairments in PwMS who appear clinically stable or radiologically unchanged, making them suitable for interim monitoring of patients between scheduled MRI examinations. In this context, multimodal prediction frameworks could, once externally validated, be applied during routine follow-up to flag patients at increased risk of disability worsening despite apparent clinical and MRI stability. Such risk stratification could motivate closer monitoring, adapted follow-up intensity, or additional functional assessments, and may support earlier consideration of treatment adaptation, particularly in light of emerging therapeutic strategies targeting progression. Early identification of individuals at risk for worsening disability, therefore, has the potential to inform both monitoring strategies and treatment decision-making in a manner aligned with a patient’s prognostic profile (9, 92).

Finally, the integration of MRI radiomics and EPTS highlights the utility of ML approaches for combining high-dimensional structural and functional data. Such multimodal quantitative patterns, which arise from complex mathematical transformations of the input data, are difficult to detect solely through visual inspection or conventional summary measures. ML methods can therefore augment clinical interpretation by capturing non-linear relationships and interactions that are not readily apparent to the clinician’s eye. At the same time, the features contributing to these predictions should be viewed as mathematical descriptors of complex patterns rather than as directly interpretable physiological markers. Consequently, the present work represents a proof of concept, and our findings provide quantitative support for multimodal follow-up strategies in MS and motivate future validation of multimodal, high-dimensional ML approaches in larger, multicenter cohorts.

This study has several limitations that must be considered. First, the small sample size (127 patients) and low proportion of disability progression events (5.2%) limited statistical power and introduced substantial class imbalance, and within-patient dependencies. This phenomenon reflects the proactive treatment environment of a tertiary MS center, where fewer patients accumulate measurable disability. To counter this limitation, we employed a repeated stratified patient K-fold cross-validation, ensuring that episodes from one patient never appeared in both training and test sets. While this method introduces some bias due to its repetitive nature, it effectively addresses the inherent variance of our small cohort. Additionally, ADASYN oversampling was applied to synthesize additional minority-class examples, enhancing our analysis.

Although these steps partially mitigate overfitting and imbalance-related bias, we acknowledge that they do not fully compensate for the small sample size and skewed distribution of outcomes. The large whiskers observed in our boxplots reflect this variability across validation folds and highlight the fragility of performance estimates in a limited dataset. Future work should aim to validate these findings using larger, multicenter datasets, thereby enhancing statistical power and the robustness of performance estimates. Nevertheless, we evaluated model performance using a variety of metrics, including both discrimination (e.g., AUROC, AP) and calibration measures (e.g., Brier score), to examine the trade-off between competing objectives. Across nearly all of these metrics and their combinations, the multimodal framework encompassing both MRI radiomics and EPTS features, with or without clinical variables, consistently outperformed single-modality models.

Another limitation is the exclusion of DMT use as a predictive feature, which may have introduced confounding, as different treatment regimens significantly impact progression risk. This omission was due to incomplete or heterogeneous treatment data in this retrospective cohort. Nevertheless, the high proportion of patients receiving moderate-to-high efficacy therapies reflects contemporary clinical practice and underscores that measurable disability worsening can still occur despite effective treatment. Future models should incorporate detailed longitudinal treatment data, such as DMT type, timing, and duration, to enhance clinical relevance. Moreover, our data-driven approach may not sufficiently account for existing biases in the dataset, such as the notable gender imbalance inherent to MS prevalence, which could distort model outputs.

Third, disability worsening was defined using a binary outcome based on a single EDSS assessment approximately two years after the baseline assessment. Confirmed disability progression could not be reliably assessed in the present retrospective cohort due to limited visit density. As a result, this outcome definition may capture transient EDSS fluctuations and does not fully exploit the granularity of the continuous EDSS scale, and may be less sensitive in patients with higher baseline EDSS scores. Nevertheless, this pragmatic choice was made to ensure a sufficient number of events and comparability with prior real-world prognostic studies. Future work should incorporate denser longitudinal follow-up to enable confirmed progression endpoints over a longer period and continuous modelling of EDSS trajectories to avoid bias in outcome ascertainment.

Fourth, conducting a hyperparameter search prior to the cross-validation could introduce slight bias, as hyperparameters may have been exposed to testing sets, potentially inflating observed performance metrics. However, because this process was uniformly applied across all configurations, the comparative power of the study (EPTS vs. MRI and their combination) remains intact.

Fifth, the use of multiple MRI acquisition protocols and different EP devices for data collection inevitably introduces heterogeneity in feature extraction. We sought to mitigate this through MRI harmonization and through stringent standardization of EP acquisition and processing at the center of data collection, thereby minimizing variability across devices. Nevertheless, the presence of multiple MRI acquisition protocols necessitated the use of a SRR approach using PRETTIER (45), as radiomics analysis requires three-dimensional MRI data (24). While our results indicate that radiomics features derived from the NAWM contributed more strongly to the predictive models, this should not be interpreted as evidence that lesion-based features are of lesser biological relevance. Rather, it is plausible that bias introduced during image reconstruction and subsequent automated lesion segmentation affected lesion representation, thereby reducing the apparent importance of lesion radiomics features.

Sixth, an important limitation of this study is the absence of external validation. All data were acquired within a single center and a single country, using a Philips 1.5T scanner with three distinct acquisition protocols. Independent external validation in larger, multicenter cohorts is required to confirm the robustness and transportability of the results. Although MRI reconstruction using SRR was applied to reduce technical variability, the generalizability of the proposed models to other centers, scanners, and patient populations cannot be assessed. Accordingly, our findings should be interpreted as proof-of-concept evidence supporting the potential value of multimodal data integration rather than as a clinically deployable prediction tool. We further note that harmonization approaches similar to those applied here have been successfully used in multicenter neuroimaging studies and may facilitate future external validation across scanners and sites (24). However, such validation is beyond the scope of the present dataset.

Lastly, the complex preprocessing required for obtaining both MRI radiomics and EPTS features poses challenges for generalizability and limits the clinical applicability of our findings. The limited interpretability of high-dimensional features, alongside their derivation at a group level, presents obstacles to effectively translating model outputs into individual clinical decisions. Nonetheless, these features have the potential to capture subtle structural and functional abnormalities that clinicians often perceive during examinations but may struggle to quantify objectively.