Headache specialists at the Headache Clinic of the Taipei Veterans General Hospital prospectively surveyed patients aged 20–60 years with newly-diagnosed migraine from May 2011 to Jan 2017. The diagnosis of migraine followed the International Classification of Headache Disorders-3 criteria . We included patients with high frequency episodic migraine and CM, and excluded patients with severe depression, i.e., Hospital Anxiety and Depression Scale-Depression (HADS-D) score ≥ 15, or those with comorbid medication overuse, i.e., using abortive headache medication for more than 10 or 15 days per month . We recruited age- and sex-matched HC for comparison.
All participants were right-handed, denied any history of systemic or major neurological diseases, and presented with normal neurologic examinations. We encrypted information that could potentially expose individual identity.
Standard protocol approvals, registrations, and patient consents
All participants completed informed consent forms after receiving a complete explanation of the study. The Institutional Review Board of our hospital approved the study protocol.
All participants filled out a semi-structured questionnaire at their first visit to obtain demographic information and headache profiles. We defined headache frequency as the average number of headache days per month in the last 3 months and duration of headache history as the duration in years between the first migraine episode and the date of the first visit to our headache clinic. We evaluated mood based on severity of depression and anxiety using the HADS  and functional disability caused by migraine using the Migraine Disability Assessment Scale (MIDAS) .
Each participant underwent scheduled MRI during the interictal period, defined as absence of acute migraine attack within 2 days prior and subsequent to the date of image acquisition. We rescheduled the scanning session if there was an acute migraine episode during this period or use of analgesics, triptans, or ergots for any reason within 48 hours prior to scanning.
After the MRI scanning, all included patients were treated at our headache clinic by headache specialists according to their clinical experience. Two years after the first visit, physician interviewed the patients by telephone to assess their migraine status and headache frequency within the last 3 months. Based on the comparison of current headache frequency to that at first visit, we defined good outcome as ≥ 50% reduction in baseline headache frequency and poor outcome as lower than 50% reduction in headache frequency or as frequency increase. Medical records were also reviewed to assess the use of migraine medication in these patients. Among the 56 patients included in our analysis, 48 were treated with both migraine prophylactic (including topiramate, propranolol, metoprolol, flunarizine, valproic acid, tricyclic acid, or amitriptyline) and abortive medications (sumatriptan or non- steroidal anti-inflammatory drugs), while the other 8 patients were treated with abortive medications only.
MRI data acquisition
We used the same 3.0 T GE Discovery MR750 scanner (General Electric Healthcare, Milwaukee, WI, USA) at Taipei Veterans General Hospital to acquire all data with a standard eight-channel phase array head coil.
We acquired the T1-weighted anatomical scans with two different acquisition pulse sequences including: 1) an inversion recovery prepared fast spoiled gradient-recalled echo sequence (IR-FSPGR) and 2) an IR-FSPGR-brain volume imaging (BRAVO) sequence.
Of note, 32 of 56 patients with HFM and 29 of 37 HC in this study were also included in a dataset that was recently used to analyze hippocampal volume changes in patients with migraine ; these patients underwent the first acquisition pulse sequence, whereas the rest and the HC underwent the second acquisition pulse sequence.
The detailed imaging parameters of each protocol were as follows: repetition time/echo time/inversion time = 9.4 (IR-FSPGR sequence) or 9.2 (IR-FSPGR-BRAVO sequence) / 4.0 (IR-FSPGR sequence) or 3.7 (IR-FSPGR-BRAVO sequence) / 450 ms; flip angle = 12°, matrix size = 256 x 256, field of view = 256 x 256 mm2, number of excitations = 1, slice thickness = 1 mm without inter-slice gap and interpolation, and 172 (IR-FSPGR sequence) or 168 (IR-FSPGR-BRAVO sequence) axial contiguous slices. An experienced neuroradiologist visually inspected all MRI scans to exclude any organic brain disorders; no participant was excluded for brain abnormalities. Before subsequent image processing, we reoriented all T1-weighed scans using a center-of-mass approach to minimize the position difference during the data acquisition.
Calculation of gray matter volume information
To estimate individual voxel-wise gray matter volume (GMV) maps, we processed T1-weighted scans using the voxel based morphometry (VBM)  pipeline with Statistical Parametric Mapping 12 (SPM12, version 7487, Wellcome Institute of Neurology, University College London, UK) under the MATLAB environment (version R2015b; Mathworks, Natick, MA). We corrected individual T1-weighted scans for intensity inhomogeneities, segmented them into GM, white matter (WM), and cerebrospinal fluid (CSF); and initially rigid-aligned them to the Montreal Neurological Institute (MNI) space using the SPM12 “Segment” function. To improve tissue classification accuracy of the subcortical areas, which are highly involved in the pathophysiology of migraine, we incorporated the enhanced tissue probability maps for subcortical regions to the above tissue segmentation procedure . We then used the Diffeomorphic Anatomical Registration through Exponentiated Lie algebra (DARTEL) toolbox to generate study-specific tissue templates, by iteratively registering the ridge-aligned GM and WM segments of all study participants, and to further warp individual tissue segments to the constructed templates . We modulated the individual MNI-space GM tissue segments with the corresponding DARTEL flow field to ensure that the following statistical analyses would be more sensitive to local GMV changes. Finally, we smoothed the modulated GMV maps with an isotropic Gaussian filter (full width at half maximum = 8 mm) and we further excluded voxels with GM probability lower than 0.2. We set the final spatial resolution of all GMV maps to 1.5 mm3. We estimated the individual global tissue volume and total intracranial volume (TIV = GM + WM + CSF volumes) in native T1 space to adjust for the effect of global brain size in the following statistical analyses.
Minimization of the influence of different acquisition protocols using data harmonization modeling
In this study, we used the ComBat harmonization approach to reduce the potential influence of different image acquisition protocols in GMV measurements . ComBat was originally designed to correct for “batch effects” in genomic studies, and recent multi-site neuroimaging studies adapted this approach to remove unwanted non-biological variability while preserving meaningful associations between image variables and covariate of interests [21-23]. This approach uses a multivariate linear mixed effects regression with terms for biological variables and imaging protocols to model quantitative measurements (for example: voxel-wise GMV maps in the current study). In more detail, we included “group” as a covariate of interest to preserve potential biological trends in the data and simultaneously corrected for the effect of different acquisition protocols.
1. Analyses of demographic and clinical data
The descriptive data in the demographic and clinical profiles are presented as mean ± standard deviation or numbers and percentages. We used the chi-square test to test for differences in categorical data. We used Student’s t test to compare the means of normally-distributed continuous variables, and the Mann–Whitney U test to compare non-normally-distributed variables, i.e., the MIDAS scores.
2. Analyses of voxel-wise imaging-based investigations
We used SPM12 with appropriate statistical models to perform the following voxel-wise statistical analyses. We adopted the cluster-extent statistical approach with the updated “3dFWHMx” and “3dClustSim” functions (available in the Analysis of Functional Neuroimages software, version 19.1.18; 10,000 Monte Carlo simulations with explicit GM mask) to correct for multiple comparisons across the whole-brain voxels. An initial voxel-level P-value < 0.005 with 257 extended voxels was considered statistically significant at cluster-level family-wise error (FWE) rate-corrected P-value < 0.05 for all voxel-wise statistical analyses (including VBM and structural covariance (SC) network analysis). The details of the statistical models are listed below.
2.1 Identification of GMV alterations among patients with HFM with different outcomes and HC
To identify the GMV difference among patients with HFM with different outcome status and the HC, a single-factor three-level (HFM with good outcome, HFM with poor outcome, and HC) analysis of covariance (ANCOVA) design was employed, with age, sex, TIV, and HADS scores entered as nuisance variables. The following contrasts were tested: HFM vs. HC, HFM with good outcome vs. HC, HFM with poor outcome vs. HC, and HFM with good outcome vs. HFM with poor outcome. We extracted, averaged, and correlated the GMV of the clusters with a significant between-group effect with headache profiles using partial Pearson’s correlation to investigate the clinical relevance. We entered age, sex, TIV, and HADS scores as nuisance variables in the correlation analysis.
2.2 Exploration of the changes in the SC network in patients with different outcomes
SC is a volumetric correlation measurement between two brain regions. SC network analysis was recently proposed as a surrogate approach to characterizing structural connectivity profiles between distinct anatomic brain regions . We generated SC networks using hypothesis-driven, seed-based correlation analysis  or data-driven, independent component analysis. SC networks reflect the shared covariance of brain morphologic features across the study participants and provide a quantitative means to studying cortical morphometric organization. In addition, the topology of the SC network is highly concordant with gene expression patterns and recapitulates intrinsic functional network architecture [26, 27].
We used three cerebellar regions with significant GMV alterations between the outcome groups as predefined regions of interest (ROIs) for the SC analyses. We constructed three multiple regression models for the respective ROIs to explore potential differences in the SC network between the outcome groups . We integrated a group main effect term (good vs. poor outcome), a mean ROI volume main effect, and a group x mean ROI volume interaction term into the statistical model and included age, sex, TIV, and HADS scores to adjust for potential nuisance effects. The changes in SC integrity of the predefined ROI with identified clusters between the outcome groups could be identified by evaluating the statistical significance of the corresponding interaction term of the constructed models.
3. Investigation of SC integrity and long-term headache frequencies
To investigate if the changes in SC integrity between the predefined ROIs and identified clusters were correlated to 2-year headache frequencies, we applied a recently-proposed approach to obtain a single measure that could quantify the integrity of SC for each individual . We correlated the SC measures for each patient with their 2-year headache frequencies using partial Pearson’s correlation after controlling for age, sex, and HADS scores.
4. Predictive values of clinical profiles and neuroimaging data for headache outcomes
We constructed two logistic regression models to assess the predictive values of clinical profiles only and clinical profiles with neuroimaging results for headache outcomes, respectively. We estimated the predictive values of the regression models using the area under the receiver operating characteristic curve (AUC).
We performed all statistical analyses using SPSS version 21.0 for Windows (IBM Corp., Armonk, NY), and a P value < 0.05 was regarded as significant.