Participants

Between November 2021 and April 2022, a total of fifty-seven right-handed individuals were enrolled in the present study, including 27 patients with EM—20 without aura (MWoA) and 7 with aura (MWA)—recruited consecutively from the Headache Clinic of First Affiliated Hospital of University of Science and Technology of China (USTC; Anhui Provincial Hospital), and 30 HCs recruited from the local community via advertisements. EM (with and without aura) was diagnosed according to the International Classification of Headache Disorders-III (ICHD-III) [22]. No migraine preventive medication was used by the participants in the past 3 months. The inclusion criteria for patients and controls included 18–55 years of age, right-handedness, and Han ethnicity. The exclusion criteria were as follows: (I) presence of other neurological diseases; (II) a history of significant physical or psychiatric illnesses; (III) a history of head injury with loss of consciousness; and (IV) pregnancy or any contraindications for MRI. To avoid measuring imaging changes associated with acute migraine symptoms, all patients were scanned during an interictal period, at least 72 h after and 24 h prior to a migraine event. The study procedures were approved by the Ethical Committee of First Affiliated Hospital of USTC and complied with the Declaration of Helsinki. Written informed consent was obtained from all participants before study entry.

Clinical assessment

Demographic information of the participants (including age, sex, years of education) was recorded. Migraine family history, migraine duration, the Headache Impact Test-6 (HIT-6) [23], and a visual analogue scale (VAS) [24] were used to assess the impact of migraine. The 14-item Hamilton Rating Scale for Anxiety (HAMA) [25] and the Beck Depression Inventory, 2nd edition (BDI-II) [26] were applied to assess the anxiety and depression status of the patients. The Montreal Cognitive Assessment (MoCA) [27] was used to evaluate cognitive function. The Big Five Inventory–60 items (BFI) measures personality traits [28]. Several subjects did not undergo the whole clinical assessment due to personal reasons, e.g., low education level or insufficient time.

MRI acquisition

All the subjects underwent MRI scans on a GE 3.0 T MR system (DISCOVERY MR750, GE Healthcare, Milwaukee, WI, USA) with a 24-channel head coil at the MRI Center of The First Affiliated Hospital of University of USTC. Earplugs were used to reduce scanner noise, and tight but comfortable foam padding was used to minimize head motion. Before the scanning, all subjects were instructed to keep their eyes closed, relax, move as little as possible, think of nothing in particular, and not fall asleep during the scans. High-resolution, three-dimensional (3D), T1-weighted structural images were acquired using a brain volume (BRAVO) sequence with the following parameters: repetition time (TR) = 8.5 ms; echo time (TE) = 3.2 ms; flip angle (FA) = 12°; field of view (FOV) = 256 mm × 256 mm; matrix = 256 × 256; slice thickness = 1 mm, no gap; 144 axial slices; and acquisition time = 240 s. Resting-state BOLD data were acquired using a gradient-echo single-shot echo planar imaging (GRE-SS-EPI) sequence with the following parameters: TR = 2,000 ms; TE = 30 ms; FA = 90°; FOV = 240 mm × 240 mm; matrix = 64 × 64; slice thickness = 4 mm without gap; 36 interleaved axial slices; 240 volumes; and acquisition time = 480 s. DTI data were acquired using a spin‒echo single-shot echo planar imaging (SE-SS-EPI) sequence with the following parameters: TR = 5, 260 ms; TE = 99 ms; FA = 90°; FOV = 220 mm × 220 mm; matrix = 128 × 128; slice thickness = 5 mm; slice gap = 1 mm; 19 axial slices; 25 diffusion gradient directions (b = 1000 s/mm2) plus five b = 0 reference images; and acquisition time = 142 s. In addition, conventional MRI examination was underwent to exclude the subjects with cerebral infarction, malacia, or occupying lesions.

Definition of thalamic subregions

The thalamic subregions were defined according to a connectivity-based parcellation study using multimodal neuroimaging techniques [29]. In each hemisphere, the thalamus was parcellated into the medial prefrontal thalamus (mPFtha), premotor thalamus (mPMtha), sensory thalamus (Stha), rostral temporal thalamus (rTtha), posterior parietal thalamus (PPtha), occipital thalamus (Otha), caudal temporal thalamus (cTtha), and lateral prefrontal thalamus (lPFtha). Thus, we defined a total of 16 regions of interest (ROIs) for the bilateral parts of the thalamus (Fig. 1).

Fig. 1

Illustration of subregions of the bilateral thalami. Abbreviations: mPFtha, medial prefrontal thalamus; mPMtha, premotor thalamus; Stha, sensory thalamus; rTtha, rostral temporal thalamus; PPtha, posterior parietal thalamus; Otha, occipital thalamus; cTtha, caudal temporal thalamus; lPFtha, lateral prefrontal thalamus; L, left; R, right

fMRI data preprocessing

The functional data were preprocessed and analysed using the Statistical Parametric Mapping software (SPM12; http://www.fil.ion.ucl.ac.uk/spm) and the Data Processing and Analysis for Brain Imaging (DPABI_ V3.1_180801; http://rfmri.org/dpabi) [30] in MATLAB R2016b (MathWorks, Inc.) as follows: (1) Removal of the first 10 volumes of the resting-state functional images; (2) Slice timing correction; (3) Head motion correction; (4) Regression of several nuisance covariates (linear drift, estimated motion parameters based on the Friston-24 model, spike volumes with FD > 0.5, white matter signal, and cerebrospinal fluid signal) from the data; (5) Data detrending and bandpass-filtering from 0.01 to 0.1 Hz; (6) Spatial normalization using DARTEL; and (7) Data smoothing with a Gaussian kernel of 6 × 6 × 6 mm3 full-width at half-maximum (FWHM). We also calculated FD, which indexes the volume-to-volume changes in head position.

The rsFC analysis was processed using DPABI software (V3.1_180801). For each individual, Pearson's correlation coefficients between the mean time courses of each thalamic subregion and those of each voxel in other parts of the brain were computed. Then, the correlation coefficients were converted into Fisher's z values to improve normality. For each group, individuals' z values were then entered into a random-effect one-sample t test in a voxelwise manner to identify brain regions that showed significant positive correlations with each ROI. Finally, a 2-sample t test was performed within the positive rsFC mask to quantitatively test group differences in the rsFC of each ROI. Multiple comparisons for these analyses were corrected using a cluster-level familywise error (FWE) method with a corrected threshold of P < 0.05.

GMV calculation

Structural scans were processed using CAT12 (CAT12, http://www.neuro.uni-jena.de/cat/) for SPM12 in MATLAB R2016b for VBM analysis). VBM includes spatial normalization, segmentation and smoothing. In brief, each participant's original T1 image was spatially normalized and segmented into grey and white matter and cerebrospinal fluid (CSF). After data preprocessing, the modulated normalized GMV was smoothed using a 6 mm FWHM Gaussian kernel.

The GMV of each thalamic subregion was extracted and compared between the two groups using the two-sample t test. Multiple comparisons were corrected using the Bonferroni method with a significance threshold of P < 0.05/16 = 0.003 (16 thalamic subregions). Moreover, the GMV of the whole thalamus was also compared between the two groups, and a p value < 0.05 was considered significant.

DTI analysis

The DTI datasets were preprocessed with the FMRIB Software Library (FSL v6.0.1, https://fsl.fmrib.ox.ac.uk/fsl/fslwiki) [31, 32]. The FSL Diffusion Toolbox (FDT) was used to correct eddy current distortions and head motion. The brain extraction tool (BET) was used to create brain masks from the b0 images. An automated quality control framework was used to assess the diffusion MRI data [33]. AD (axial diffusivity), FA, MD, and RD (radial diffusivity) were calculated by using the FSL toolbox DTIFIT. These images were then coregistered to each subject’s T1-weighted images using the FLIRT linear registration tool, yielding the normalized FA, AD, RD, and MD. Finally, the mean FA, AD, RD, and MD values of each thalamic subregion were extracted and compared between the two groups using the two-sample t test. Multiple comparisons were corrected using the Bonferroni method with a significance threshold of P < 0.05/16 = 0.003 (16 thalamic subregions).

Correlations between imaging and clinical parameters

To determine whether thalamic rsFC, GMV, and DTI abnormalities of the thalamic subregions with significant intergroup differences were associated with illness duration and symptom severity (HIT-6, VAS, HAMA, BDI-II, MoCA, and BFI scores), we calculated partial correlations (two‐tailed) in the migraine group after controlling for age, sex, education, TIV, and FD to explore the association between the mean values extracted from each significantly different region and clinical parameters. A p value < 0.05 was considered significant.

Sample size calculation

As our primary goal was to detect differences in the functional connectivity of each thalamic subregion between EMs and HCs, the sample size was calculated based on pilot data from 20 subjects, ten for each group. For the EM group, the rsFC between L-mPMtha and Frontal_Mid_L was 0.13 ± 0.119, between L-rTtha and Frontal_Sup_L was 0.18 ± 0.062, and between L-PPtha and Precuneus_R was 0.25 ± 0.142. For the HC group, the rsFC between L-mPMtha and Frontal_Mid_L was 0.31 ± 0.101, between L-rTtha and Frontal_Sup_L was 0.40 ± 0.135, and between L-PPtha and Precuneus_R was 0.45 ± 0.106. To achieve a desired power of 90% with a significance level of 5%, the required sample size was 10 subjects for each group as calculated by PASS software (https://www.ncss.com/software/pass). For a more conservative estimate, we decided to complete the enrolment when 27 EM and 30 HC subjects had been included in the current study.