Study participants

We recruited 50 patients with episodic migraine and 50 age- and sex-matched healthy controls from an academic headache clinic between August 2017 and July 2018 [30]. The inclusion criteria for patients were as follows: 1) age 18–50 years, 2) not taking preventive medications, and 3) premenopausal status in female patients. The exclusion criteria were as follows: 1) chronic migraine, medication-overuse headaches, chronic pain disorders other than migraine, and psychiatric disorders such as bipolar affective disorder or schizophrenia; 2) contraindications for 3 T MRI, including use of a tissue expander, pacemaker, non-detachable metal objects, orthodontic devices, electrical leads, or implants in the body; 3) pregnancy; 4) claustrophobia requiring sedation during MR scanning; 5) inability to report headache or complete the headache diary due to cognitive decline; and 6) disagreement with the study procedures. Controls were recruited using the same inclusion and exclusion criteria after confirming that they had no history of memorable headache episodes in their lifetime or during the last year and no family history of migraine. Migraine was diagnosed by a headache specialist (MJL) according to the ICHD 3rd version [9]. The study was approved by the Samsung Medical Center Institutional Review Board, and all participants provided written informed consent.

Study scheme

At baseline, the patients completed headache diaries and clinical evaluations and underwent structural and functional MRIs following the protocol described below. Headache status within 2 days of MRI scanning was obtained. Patients were considered ictal if they had headaches of any intensity on the day of scanning, interictal if they were headache-free for 2 days before and after scanning, and peri-ictal if they were headache-free on the day of scanning, but developed headaches within two days of scanning. The same participants underwent multimodal MRI scanning and clinical evaluation using the same protocol one year after the baseline scan. This study was registered at ClinicalTrials.gov (NCT03487978).

MRI acquisition

Participants underwent T1-weighted and resting-state functional MRI (rs-fMRI) using a Philips Ingenia 3 T scanner. The acquisition parameters for T1-weighted imaging were as follows: repetition time (TR) = 9.9 ms; echo time (TE) = 4.6 ms; flip angle = 8°; field of view (FOV) = 240 mm × 240 mm; the number of slices = 180 (reconstructed to 360); slice thickness = 1 mm (reconstructed to 0.5 mm); and in-plane resolution = 1 mm × 1 mm (reconstructed to 0.5 mm × 0.5 mm). The rs-fMRI data were acquired as follows: TR = 3,000 ms; TE = 30 ms; flip angle = 90°; FOV = 220 mm × 220 mm; number of slices = 33; slice thickness = 4 mm; in-plane resolution = 3 mm × 3 mm (reconstructed to 1.25 mm × 1.25 mm); and number of volumes = 200. During the scan, the participants were instructed to keep their eyes open and focus on a fixation cross to prevent falling asleep.

Data preprocessing

The imaging data were preprocessed using the fusion of a neuroimaging preprocessing (FuNP) surface-based pipeline incorporating FSL, AFNI, FreeSurfer, ANT, and Workbench (https://github.com/CAMIN-neuro/FuNP/) [31,32,33,34,35]. The gradient nonlinearity and B0 distortions were corrected from the T1-weighted images, non-brain tissues were removed, and intensity was normalized. Pial and white surfaces were generated, and the two surfaces were averaged to obtain a mid-thickness surface [33, 36,37,38]. For the rs-fMRI data, preprocessing involved the following steps: the first four volumes were discarded to allow for magnetic field saturation, and slice timing was corrected. Volumes with large head motion (i.e., framewise displacement > 0.5 mm) were discarded [39], and motion correction was performed. Skull removal and intensity normalization were performed, and nuisance variables, including head motion, white matter, cerebrospinal fluid signals, cardiac and respiration signals, and large vein-related contributions, were removed using the FMRIB’s independent component analysis-based X-noiseifier (ICA-FIX) [40]. The rs-fMRI data were registered onto preprocessed T1-weighted structural data and subsequently to the Montreal Neurological Institute standard space. The band-pass filter with a frequency range between 0.009 and 0.08 Hz was applied. Finally, the preprocessed volumetric rs-fMRI data were mapped to the cortical surface using a cortical ribbon-constrained volume-to-surface mapping algorithm, and spatial smoothing with a full width at half maximum of 5 mm was applied.

Multimodal imaging features

We calculated the morphological features from preprocessed T1-weighted MRI and the functional features from rs-fMRI (Fig. 1).

Fig. 1

The migraine diagnostic model constructed using multimodal MRI features. a Schema of the research pipeline to construct a multimodal MRI-based migraine diagnosis model using data at baseline and validating its performance using data at follow-up. b The concepts of the features, including curvature, sulcal depth, cortical thickness, functional gradient, degree centrality, and betweenness centrality, are shown. Mean values of the features at baseline are shown on brain surfaces

Morphological features

The morphological features included cortical thickness, curvature, and sulcal depth. Cortical thickness assesses the thickness of gray matter by measuring the distance between the pial and white surfaces. Curvature was defined as the degree of folding of the cortical surface, while the sulcal depth was the vertical distance between the sulci and pial surfaces. Vertex-level data were mapped to brain regions defined using the Shaefer atlas with 200, 300, and 400 parcels [41] to assess the robustness of the model across spatial granularities.

Functional features

From the rs-fMRI data, a functional connectivity matrix was constructed by calculating the linear correlations of the time series between different brain regions defined by the Schaefer atlas containing 200, 300, and 400 parcels [41]. The correlation coefficients were calculated using Fisher’s r-to-z transformations. The nodal degree and betweenness centrality were computed after leaving the top 10% of the elements per node. Degree centrality was defined as the sum of all edge weights connected to a given node, while betweenness centrality was defined as the number of shortest paths between any two nodes running through that node. Additionally, we generated a low-dimensional representation of functional connectivity (gradient) by applying a nonlinear dimensionality-reduction technique (diffusion map embedding) [42, 43]. The functional gradient represents the hierarchy of the brain because it shows an axis that distinguishes between the low-level sensory and higher-order default mode regions [42]. First, we generated a group-level template gradient from the averaged functional connectivity matrix and aligned the individual gradients to the template gradients using Procrustes rotation [43, 44].

Feature selection

Because the size of each imaging feature was equal to the number of brain regions (200, 300, or 400), the combination of \(N\) imaging features resulted in \(N \times\) number of brain regions. When constructing machine learning models, too many features commonly result in overfitting problems. Thus, a feature selection process is required to identify a subset of features to build a classification model to distinguish patients with migraine from healthy controls. We utilized the elastic net regularization method, which combines the properties of the least absolute shrinkage and selection operator (LASSO) and ridge regression models [45]. It adjusts the weights of each feature by penalizing the coefficients of the model with the L1 (LASSO) and L2 (ridge) norms, as follows:

$$\left(\widehat_},\widehat}}\right)=\underset_,}}}\left(\sum_^_-_-}}^}}}}_}\boldsymbol}^+\lambda \sum_^\left(\frac_^+\alpha \left|_\right|\right)\right)$$

where \(}}_}}\) is the feature set of \(^\) subject and \(_\) is the dependent variable indicating the disease state. \(_\) denotes the regression coefficient of \(^\) feature, \(\lambda\) determines the overall strength of the regularization, while \(\alpha\) determines the relative weight of the ridge and LASSO penalties (i.e., L1 ratio). We determined the optimal \(\uplambda\) using greedy search algorithms by varying the L1 ratio as . A hundred \(\uplambda\) values were generated at each L1 ratio, and the optimal \(\uplambda\) showing the best performance was selected. The feature selection procedure was performed using the baseline data, and the selected features were applied to the follow-up data.

Diagnostic model

We constructed a multimodal MRI-based diagnostic model to distinguish patients with migraine from healthy controls using a random forest classifier with selected features. Given the constraints of small sample size and the high dimensionality of imaging features, the random forest classifier, an ensemble-based machine-learning approach that aggregates decisions from multiple decision trees, is an appropriate option as it effectively captures feature interactions in high-dimensional data [46]. It is more suitable than deep learning or boosting methods, which are prone to overfitting and unstable training with limited data. In this study, a hundred trees and the maximum depth at which all leaves become pure nodes were considered. The model was trained using all the subjects at baseline. We conducted a classification task with a five-fold cross-validation. Specifically, 100 participants were randomly divided into five subsets of equal size (i.e., 20 participants per subset), ensuring that the proportions of patients and controls were balanced. Four of the five folds were used as training data, while the remaining fold served as the test set. We repeated this process 100 times with different training and test datasets to mitigate potential subject selection bias. To assess the stability of the trained model across different migraine phases, we applied the model constructed using all data to patients in the ictal/peri-ictal phase and matched controls, as well as patients in the interictal phase and matched controls. In addition, the model built using all participants at baseline was applied to the follow-up data. The classification performance was assessed based on accuracy, sensitivity, specificity, F1 score, and area under the curve (AUC). The model accuracy was assessed for all designed L1 ratios and parcellation scales. We further evaluated the classification performance by configuring the features in all possible combinations of morphological (i.e., cortical thickness, curvature, and sulcal depth) and functional measures (i.e., degree centrality, betweenness centrality, and gradient).

Between-group differences in the features

To evaluate how each feature differed between the groups, we quantified the selected probability of the feature set that showed the highest average classification accuracy across baseline and follow-up, interictal and ictal/peri-ictal phases, L1 ratios, and parcellation scales. We subsequently compared the features of the brain regions selected during the feature selection process between patients with migraine and healthy controls using two-sample t-tests. We visualized the selected probabilities and t-statistics according to seven intrinsic functional networks: visual, somatomotor, dorsal attention, ventral attention, limbic, frontoparietal, and default mode [47]. Additionally, we calculated the Shapley values from the random forest classifier across various parameter settings to identify the features contributing to the classification between patients with migraine and healthy controls.

Clinical associations

The clinical applicability of our findings was demonstrated by examining the relationship between the features with the highest classification performance (i.e., curvature, sulcal depth, thickness, and gradient) and clinical measures, including headache frequency and disease duration. Specifically, we predicted the clinical scores by applying a linear regression model to the top 20 features that showed high correlations with each clinical score.

Statistical analysis

Data are presented as the mean ± standard deviation (SD), unless otherwise specified. All analyses and visualizations were performed using Python v3.10.9. Elastic net regularization and random forest classifiers were conducted using scikit-learn v1.3.2, and a two-sample t-test was conducted using SciPy v1.10.0.