Imaging data processing

All T1 imaging data will be processed and examined for voxel-based morphological (VBM) analyses using the VBM8 toolbox (http://dbm.neuro.uni-jena.de/vbm/) in Statistical Parametric Mapping (SPM) (SPM8, Wellcome Department of Imaging Neuroscience Group, London, UK; http://www.fil.ion.ucl.ac.uk/spm) running in MATLAB 2014a platform. We will perform an optimized VBM protocol in which a brain-tissue-only template will be used other than a whole-brain template. Brain Extraction Tool (integrated in MRIcro; http://www.mricro.com) was used to remove non-brain region images [25]. After removal of nonbrain regions, images will be segmented into GM, white matter (WM), and cerebrospinal fluid (CSF) by using the Diffeomorphic Anatomical Registration Through Exponentiated Lie algebra (DARTEL) algorithm [26]. Afterward, the normalized images would be averaged by using SPM mean function. To create the final template, all average images would be smoothed using an isotropic Gaussian kernel with a full width at half maximum (FWHM) of 8 mm.

All DKI images will be processed using a combination of FMRIB’s Software Library (FSL) [27] and in-house image processing tools developed in MATLAB. The diffusion dataset will be modulated to get potential 3D head motion and eddy current distortion using FSL eddy correct. The toolbox implement in MTALAB will be applied to deal with diffusional kurtosis tensors. Region of interest (ROI) of headache will be drawn by hand. DKI parameters of migraine pain-ROIs, such as mean kurtosis (MK), mean diffusivity (MD) and fractional anisotropy (FA), will be measured. All feature maps will be spatially normalized to the standard MNI space by using the transformation fields derived from tissue segmentation of structural images and resampled to 3 mm isotropic voxels.

RS-BOLD data will be preprocessed in the toolbox of MATLAB (Data Processing Assistant for Resting-State fMRI, DPARSF; http://www.restfmri.net/forum/DPARSF). SPM8, RS-fMRI data analysis toolkit (REST1.6; http://www.restfmri.net) and graph-theoretical network analysis (GRATNA; https://www.nitrc.org/projects/gretna) will also be selected to deal with the images. The first 10 time points will be discarded to ensure stable magnetization and allow the participants to adapt to the EPI scanning environment. After that, slice timing and head motion will be conducted, and scans with head motion of translation > 3.0 mm or rotation > 3° will be excluded. The functional images will then be spatially normalized to Montreal Neurological Institute (MNI) space using the transformation fields derived from tissue segmentation of individual structural images and resampled to 3 × 3 × 3 mm³. The resulting images will be smoothed with 8 mm FWHM isotropic Gaussian kernel. Linear trends will be removed from the image time series, and data will be band-pass filtered at 0.01–0.08 Hz. Finally, nuisance signals such as 24-parameter head motion profiles, white matter, and cerebrospinal fluid signals will be regressed out from each voxel’s time series to exclude non-neuronal sources [28].

Amplitude of low frequency fluctuation (ALFF), regional homogeneity (ReHo) and functional connectivity (FC) values will be estimated in the DPARSFA software that is developed to extract abnormal regions between two groups. Small-world properties and network efficiency will be calculated in the GRETNA software to compare the topological characteristics between two groups.

For the ASL fMRI data, corresponding CBF images were obtained using an automated image postprocessing tool embedded in the GE healthcare MR-750 system. Subsequently, the CBF images will be spatially normalized to the standard MNI space by using the transformation fields derived from tissue segmentation of structural images and resampled to 3 mm isotropic voxels. The resulting images will be transformed to z scores using Fisher’s transformation approach and then will be smoothed with 8 mm FWHM isotropic Gaussian kernel.