Neuroimaging features of whole‐brain functional connectivity predict attack frequency of migraine

Abstract Migraine is a chronic neurological disorder characterized by attacks of moderate or severe headache accompanying functionally and structurally maladaptive changes in brain. As the headache days/month is often measured by patient self‐report and tends to be overestimated than actually experienced, the possibility of using neuroimaging data to predict migraine attack frequency is of great interest. To identify neuroimaging features that could objectively evaluate patients' headache days, a total of 179 migraineurs were recruited from two data center with one dataset used as the training/test cohort and the other used as the validating cohort. The guidelines for controlled trials of prophylactic treatment of chronic migraine in adults were used to identify the frequency of attacks and migraineurs were divided into low (MOl) and high (MOh) subgroups. Whole‐brain functional connectivity was used to build multivariate logistic regression models with model iteration optimization to identify MOl and MOh. The best model accurately discriminated MOh from MOl with AUC of 0.91 (95%CI [0.86, 0.95]) in the training/test cohort and 0.79 in the validating cohort. The discriminative features were mainly located within the limbic lobe, frontal lobe, and temporal lobe. Permutation tests analysis demonstrated that the classification performance of these features was significantly better than chance. Furthermore, the indicator of functional connectivity had a higher odds ratio than behavioral variables with implementing a holistic regression analysis. The current findings suggested that the migraine attack frequency could be distinguished by using machine‐learning algorithms, and highlighted the role of brain functional connectivity in revealing underlying migraine‐related neurobiology.

that occur >15 days/month) (Dodick, 2018). Without effective treatments, patients with episodic migraine would experience more headache attacks and even develop chronic migraine Lipton & Silberstein, 2015). Since attack frequency is a risk factor for headache progression, it is possible that preventive medication may reduce the risk of progression to chronic migraine (Dodick, 2018;Lipton et al., 2007). The headache days per month is often evaluated in the clinical contexts by patient self-report, which may rely on patient memory Niere & Jerak, 2004). However, memories may be enhanced to various extents, causing the measurement of headache days to become unreliable or inaccurate (i.e., recall headache days may be higher than the actual headache days) (Berger et al., 2018;Haywood et al., 2018). An objective measurement of headache days per month might be helpful in the decision-making process of providing preventive medication for migraine patients.
Migraine is mainly a disorder of brain function that involves the brain regions of sensory discrimination of pain, affective-emotional processing, cognitive processing, and pain modulation (Schwedt, Chiang, Chong, & Dodick, 2015), which serve as an integrated network in its development and maintenance . Functional magnetic resonance imaging (fMRI) can be used to non-invasively acquire the data of brain activities, which would be useful to explore the central mechanisms underlying migraine (Fox & Raichle, 2007). Synchronous fluctuations in the blood oxygen-level-dependent (BOLD) signal measured by fMRI are deemed to reflect functional connectivity, or functional communication between brain regions (Rosenberg et al., 2016;Schwedt et al., 2015). Compared with the analysis using only several predefined regions or networks of interest, whole-brain functional connectivity analysis can ensure the optimal use of the pattern information of the brain (Zeng et al., 2012). Neuroimaging studies have consistently shown the abnormal brain functional network organization in patients with migraine (Liu et al., 2015;Maleki et al., 2012;Schwedt et al., 2015;Yuan et al., 2013), and these connection patterns may sharp the cognition and behavior to individual differences in migraine vulnerability (Kucyi & Davis, 2015). Hence, it is possible to locate the neuroimaging features in the brain functional network, which could be used to objectively evaluate the headache days per month (Fox & Raichle, 2007). Furthermore, the machine-learning-based predictive model does not require months of journaling to determine prognosis, and could help predict headache days for clinical settings without the presence of patient-reported ratings.
In the current study, to identify neuroimaging features that could objectively evaluate patients' attack frequency of migraine, we built a multivariate logistic regression model by using the brain functional connectivity data, in which a total of 179 migraineurs without aura (MO) were recruited from two data centers with one dataset used as the training/test cohort (N = 151) and the other used as the validating cohort (N = 28).

| MATERIALS AND METHODS
All research procedures were approved by the West China Hospital Subcommittee on Human Studies and the Medical Ethics Committee of the Affiliated Hospital of Chengdu University of Traditional Chinese Medicine and were conducted in accordance with the Declaration of Helsinki. All participants gave written, informed consent after the experimental procedures had been fully explained.

| Participants
Inclusion criteria for the migraine patients were according to the International Classification of Headache Disorders third edition (beta version) (2013): (a) migraine attacks last 4-72 hours (untreated or unsuccessfully treated); (b) featuring at least two of the following characteristics: unilateral location, pulsating quality, moderate to severe pain intensity and aggravation by causing avoidance of routine physical activity; and (c) there is nausea and/or vomiting, or photophobia and phonophobia during migraine. Exclusion criteria were: (a) any physical illness such as a brain tumor, hepatitis, or epilepsy as assessed according to clinical evaluations and medical records; (b) existence of other comorbid chronic pain conditions (e.g., tension type headache, fibromyalgia, etc.); (c) existence of a neurological disease or psychiatric disorder; (d) pregnancy; (e) use of prescription medications within the last month; (f) alcohol, nicotine, or drug abuse; and (g) claustrophobia.
A total of 151 MO patients (age: 28.69 ± 0.83 years, 39 males and 112 females) were recruited from the West China Hospital, and these patients were used as the training/test cohort. Another 28 MO patients (age: 31.43 ± 1.72 years, 8 males and 20 females) were recruited from the Chengdu University of Traditional Chinese Medicine and the surrounding community, and were used as the validating cohort. All patients were right-handed and were evaluated by a neurologist. Patients were instructed to complete a daily headache diary every evening before retiring. The diaries were returned to the researchers at the end of a 4-week period. Telephone reminders were given if diaries were not received within 2 days of the end of the four-week period. During the 4 weeks before the MRI scans, patients carefully rated the average headache intensity (0-10 scale, 10 being the most intense pain imaginable), headache days and migraine attack duration (hours) with the migraine diary. The Zung Self-Rating Anxiety Scale (SAS) and Zung Self-Rating Depression Scale (SDS) were used to quantify anxiety/depression-related symptoms of the patients. Analysis of migraine diaries was performed by two blinded evaluators.
The first five volumes were discarded to eliminate nonequilibrium effects of magnetization and to allow participants to become familiar with the scanning circumstances. Functional volumes were then slice timecorrected and realigned, and normalized to the Montreal Neurological Institute template brain, and smoothed with a 6-mm 3 isotropic Gaussian kernel. Individuals with an estimated maximum displacement in any direction larger than 1 mm or head rotation larger than 1 were discarded from the study, and no data were excluded under this criterion. A band-pass filter (0.01 Hz < f < 0.1 Hz) was applied to remove the effects of lowfrequency drift and high frequency physiological noise. Finally, regression of nuisance covariates including Friston 24 head movement parameters, cerebrospinal fluid signals and white matter signals from the fMRI data were carried out. The data used to support the findings of this study are available from the corresponding author upon request.

| Head-motion calculations
To test whether our observations would be held when considering the effects of head motion, we eliminated volumes from each subject's resting fMRI time series that were associated with sudden head motion.
An index of framewise displacement (FD) was applied to mark the volumes that tended to behave as burst noise. This would result in temporal masks for our data and similar approaches have been used in several previous rs-fMRI studies (Liu et al., 2015;Power, Barnes, Snyder, Schlaggar, & Petersen, 2012). For each subject, the flag volume was censored if its derivative values were above 0.5 (Power et al., 2012).

| Functional connectivity measures
The Human Brainnetome Atlas (Fan et al., 2016) containing 210 cortical and 36 subcortical regions of interest (ROIs) was used to create fMRI time course correlation matrices (i.e., functional connectivity matrices) for the participants' processed echo-planar image time series. Time series were extracted and averaged within each region.
For each participant, Pearson correlation coefficients were calculated between the average time series of each region in the atlas to form a 246 × 246 symmetrical matrix (30,135 unique connections). All correlation values were converted using a Fisher Z-transformation to remove the effects of global levels of correlation (Plitt, Barnes, Wallace, Kenworthy, & Martin, 2015). The GRETNA toolbox (https:// www.nitrc.org/projects/gretna), implemented in Matlab (MathWorks, Inc.), was used for the construction of functional connectivity matrices (Wang et al., 2015). The 30,135 unique connections for each participant were used as the initial features. Considering the effects of age and sex, a linear regression model was applied to obtain an age-and sex-adjusted feature for each connection.

| Grouping threshold definition
According to guidelines for controlled trials of prophylactic treatment of chronic migraine in adults, preventive medications should be given when migraine attacks are frequent (Dodick, 2018;Silberstein et al., 2008). Specifically, the migraineur who had greater than or equal to 8 headache days per month should take preventive medications, which may reduce the risk of progression to chronic migraine (Dodick, 2018;Silberstein et al., 2008). Thus, the frequency of 8 headache days/month was used as the grouping threshold in our study.

| Multivariate logistic regression model
Following the general methodology developed by Vallières et al (Vallières, Freeman, Skamene, & El Naqa, 2015), an adapted method was applied to find the migraine attack frequency-related neuroimaging features. The schematic overview can be seen in Figure 1. Specifically, the training/test cohort was divided into MOl (MO with a low attack frequency, that is, headache days < 8 per month) group and MOh (MO with a high attack frequency, that is, headache days ≥ 8 per month) group (Table 1). Multivariate models were built for the initial features of the training/test cohort and modeled outcome using imbalanced-adjusted logistic regression. The initial feature set first underwent feature dimension reduction using a ReliefF algorithm (Robnik-Šikonja & Kononenko, 2003), which reduced the feature dimension from 30,135 to 10,000. For short, every feature was assigned with a weight value indicating its relevance to the group label, and the first 10,000 features were remained. Then, Gain algorithmbased feature selection was performed using 100 bootstrap training samples to yield a reduced feature set of 25 outcome-related but lowredundancy features (Vallières et al., 2015). After feature selection, stepwise forward model construction was conducted by maximizing the 0.632+ bootstrap area under the receiver operating characteristic curve (AUC) metric in 1000 bootstrap training and testing samples (training sample: test sample≈2:1) to obtain logistic regression models combining 1 to 10 features (corresponding to model order: 1 to 10) (Vallières et al., 2015). For the outcome, the model providing the combination of functional connection features with the best parsimonious properties (i.e., the minimum model order and the maximum AUC) was chosen. Finally, the classification performance of the chosen model (i.e., best model) was estimated using the average AUC, sensitivity, specificity, and accuracy obtained in 1000 bootstrap testing samples.
To further evaluate the generalizability of the best model, we used a validating cohort (28 MO, containing 18 MOl (headache days <8) and 10 MOh (headache days ≥ 8)). The same functional connections were extracted from the validating cohort, and then fed into the logistic regression model. The AUC, sensitivity, specificity, and accuracy were computed.

| Robustness of the features
To test whether the classification performance of the true features of the best model was higher than chance level and whether counterpart random features of functional connectivity could distinguish MOl from MOh, the random connectivity matrices with the same degree distribution as the real network were generated. Specifically, the feature matrix of the training/test cohort was extracted from the best model, and the corresponding label vector (i.e., 0 for the MOl group and 1 for the MOh group) was randomly shuffled. Then, the corresponding training and test sets of the best model were used to train and test the model, respectively. As for the random features, the same dimension random connectivity matrix (i.e., 246 × 246 matrix) for each participant was first generated by using the GRETNA toolbox to ensure the same degree distribution as the real network; the random features that were located within the same location of the connection matrix were extracted to conduct permutation tests. Then, the coefficients of the logistic regression model were computed using the shuffled training sample, and the classification performance was tested on the shuffled test sample. The null distribution of the classification AUC for the training/test cohort was built.

| Additional validation
To evaluate the influence of an imbalanced grouping, we adjusted the grouping threshold to 5 headache days/month (the median of F I G U R E 1 Schematic overview of the multivariate logistic regression analysis. Whole-brain functional connectivity matrix was computed for each subject using a 246-region cortical and subcortical atlas, for a total of 30,135 unique features. The training/test cohort (N = 151) was used to build a multivariate logistic regression model with bootstrap sampling and ReliefF and Gain algorithms-based feature selection procedures. The area under the receiver operating characteristic curve (AUC) was used to select the best model that had a high AUC and few features. Then, the same features were extracted for the validating cohort (N = 28) according to the feature index of the best model. Finally, the group label for the unseen subjects was calculated using the regression coefficients of the best model headache days) of the training/test cohort to insure a balanced grouping. Specifically, the training/test cohort was divided into a MOl (MO with a lower attack frequency, that is, headache days < 5) group and a MOh (MO with a higher attack frequency, that is, headache days ≥ 5) group (Table S1). Then, the feature selection, model selection, and model generalizability evaluation procedures described in the Multivariate logistic regression model part were repeated. Note that the grouping threshold of headache days was also adjusted to 5 for the validating cohort (i.e., 28 MO, containing 15 MOl (headache days <5) and 13 MOh (headache days ≥ 5)).

| Development of a holistic evaluation model
To explore the relationship between migraine attack frequency and demographical, behavioral and neuroimaging indicators, an individual logistic regression analysis was performed. According the best classification model, a neuroimaging score was calculated for each participant via a linear combination of selected functional connectivity features weighting by its corresponding regression coefficients (Huang et al., 2016). Then, an individual logistic regression analysis was conducted with the following candidate indicators: age, sex, disease duration, average duration of a migraine attack, average pain intensity, SAS, SDS, and neuroimaging score. A backward stepwise selection was adopted and the likelihood ratio test with Akaike's information criterion was used to determine the stopping rule (Huang et al., 2016).

| Statistical analysis
As for the demographic and headache information, group comparison was carried out using a two-sample t-test and chi-square test in SPSS 20.0 (SPSS Inc., Chicago, Illinois) with a statistical power p < .05.
As for permutation tests, after 5,000 permutations, the significance value p for each metric (i.e., AUC, sensitivity, specificity, and accuracy) was computed by dividing the number of times that showed a higher value than the actual value derived from the nonpermuted model by the total number of permutations, and p < .05 was accepted (Plitt et al., 2015).

| Demographic and headache information
A total of 179 MO (i.e., the training/test cohort (N = 151) and the validating cohort (N = 28)) were enrolled in this study. For the training/ test cohort (both imbalanced and balanced grouping), there were no significant differences between the MOl and MOh groups with regard to age, sex, education, disease duration, average duration of a migraine attack, and average pain intensity (p > .05, Table 1 and S1).
For the validating cohort, the age was 31.43 ± 1.72 years (mean ± SE), and the headache days during past month were 7.82 ± 1.25 days (mean ± SE) ( Table 2).

| Functional connectivity-based classification of migraine attack frequency
The following results are for the training/test cohort using the imbalanced grouping (i.e., grouping threshold = 8 headache days/month), unless otherwise specified. No significant between group difference T A B L E 1 Demographic and headache information for training/ test cohort with an imbalanced grouping was found for the mean FD. Functional connectivity matrix for each participant was measured before the model construction. After the feature selection and model selection procedures, we found that the best model that could accurately discriminate MOh from MOl contained 8 features, which were mainly located within the limbic lobe, frontal lobe, and temporal lobe ( Figure 2,  Figure 3).

| Classification with balanced grouping
Considering the influence of an imbalanced grouping, another grouping threshold of headache days was adopted to obtain a balanced grouping for the training/test cohort (i.e., 151 MO, containing 77 MOl (headache days < 5) and 74 MOh (headache days ≥ 5), Table S1). F I G U R E 3 Classification AUC (indicated by a vertical red line) and corresponding null distribution with 5,000 random permutations. (a) The AUC was significantly greater than chance level (p < .0001) using the true features. (b) The AUC wasn't significantly greater than chance level (p = .1676) using the random features Abbreviations: CI, confidence interval; SAS, self-rating anxiety scale; SDS, self-rating distress scale. and 46.67%, respectively, which seemingly failed to distinguish MOh from MOl in the unseen dataset.

| Holistic evaluation
We finally investigated the relationship between migraine attack frequency and demographical, behavioral and neuroimaging indicators from an aspect of holistic evaluation, and we found that the functional connectivity had a higher odds ratio than other variables (Table 5).
Additionally, the β value corresponding to the neuroimaging score was positive (β = .666, p < .001, Table 5), indicating that the aberrant brain functional connectivity patterns were an influence factor for frequent migraine attacks.

| DISCUSSION
In this study, whole-brain functional connectivity was used to build a multivariate logistic regression model to find migraine attack frequency (i.e., headache days per month) related neuroimaging features.
Our model successfully predicted migraine attack frequency in independent datasets, demonstrating that the neuroimaging feature-based model is generalizable for individual disease characterization in migraineurs. Furthermore, an individual logistic regression analysis using the combined variables (i.e., age, sex, disease duration, average duration of a migraine attack, average pain intensity, SAS, SDS, and neuroimaging score) indicated that the aberrant brain functional connectivity patterns were an influence factor for frequent migraine attacks.

| Reliable identification of lower or higher frequency of migraine attacks
Preventive medications are of use to reduce the frequency and severity of attacks in people with frequent migraine (Dodick, 2018), and one of the criteria for considering or offering prevention was based on headache frequency (Lipton et al., 2007). Several studies characterize the extent of measurement error arising from rounding in headache frequency reporting , and pointed out that recall error in describing and quantifying pain was sufficiently common to distort the resulting population estimates in headache chronification research  and social factors such as anger, anger-expression, anxiety, and depression have been used to predict the frequency of migraine attacks, but only obtained an accuracy of 69.8% (Bernardy et al., 2007). Due to a lack of sensitive predictive features, these studies using behavioral indicators to predict headache attacks seemingly failed to achieve good performance.
Borsook and colleagues proposed an allostatic model to understand the effect of high-frequency headache attacks, and pointed out that repeated migraine attacks and physical or psychological stressors may cause progressive functional and structural changes in brain networks (Borsook, Maleki, Becerra, & McEwen, 2012). In our previous study, we investigated how the between-group differences of functional connections in MO were organized along with the changing trend by using resting-state fMRI, and found that the presence of chronic headache altered the functional connectivity from the local central nervous system due to a disruption in whole-brain networks with increased disease duration (Liu et al., 2015). In the current study, we extended our previous findings and observed that the variability of headache days in migraineurs could be predicted from whole-brain functional connectivity.

| Consideration of the grouping threshold
According to the guidelines for controlled trials of prophylactic treatment of chronic migraine in adults (Dodick, 2018;Silberstein et al., 2008), a frequency of 8 headache days per month was used as the grouping threshold to subdivide our training/test cohort into a MOl group (102 participants) and MOh group (49 participants). For better learning from imbalanced data, prior probabilities for each group were set equally with an imbalance-adjusted bootstrap approach (Vallières et al., 2015;Zhou et al., 2017

| Migraine attack frequency classificationrelated features
In the present study, by using a multivariate analysis approach, discriminative functional connections were found to be related to many brain areas (i.e., the temporal lobe, limbic lobe, prefrontal lobe, and insula) that belong to different functional brain systems, rather than with activity in dedicated "pain centers" within the brain. Previous studies have repeatedly shown that multiple brain regions are involved in sensory, emotional, and evaluative aspects of pain processing in healthy individuals (Martucci & Mackey, 2018). These regions of the brain include the primary and secondary somatosensory cortices, primary and supplementary motor cortices, anterior cingulate cortex, prefrontal and parietal cortices, and limbic system (Martucci & Mackey, 2018). Additionally, abnormal brain activation of the limbic lobe, frontal lobe and temporal lobe was also observed in chronic pain, including migraine (Schwedt & Dodick, 2009). Insights from positronemission tomography have observed that the limbic regions were activated during migraine attacks Matharu et al., 2004). It suggested that the limbic regions may modulate the affective components of pain in migraineurs, and may play an important role in the association between chronic migraine and psychiatric disturbances. Kucyi and Davis (2015) pointed out that the prefrontal cortex may be act as an "exchange hinge" for distributing pain-related information and relate to painattention interactions in chronic pain. Recently, our group found that baseline gray matter volume of the medial prefrontal cortex and its functional connectivity could predict a future placebo response in an 8-week sham acupuncture treatment for migraine Liu, Ma, et al., 2017). However, there is little specificity or congruence of these regions across studies in migraine processing.
Looking at the brain as an integrative complex system, these discriminative features were optimized to jointly predict the patients' status, suggesting that migraine attacks may be associated with an abnormal integrated network configuration, rather than in one or more isolated brain circuits. It highlighted the importance of datadriven analysis methods, which do not restrain the features to a priori nodes of interest. Our findings were built upon previous literature suggesting that multivariate analysis approaches in classifying functional connectivity data offered a valuable technique in understanding network-level differences in migraine attack frequency-related neurobiology and extended previous findings by demonstrating the heterogeneity in resting-state networks of individuals with migraine with different headache frequency (Vallières et al., 2015;Zhou et al., 2017). Our results also demonstrated good utility for distinguishing different patient groups, which may have potential to further understand the functional mechanisms contributing to the development and maintenance of migraine.
It should be noted that our results only indicated that the observed network features could identify the lower and higher migraine attacks migraineurs. It was not necessarily implied that the neural activity from the limbic lobe, frontal lobe and temporal lobe was derived corresponding to the neural activity generating migraine.
As our current neuroimaging investigation is not a longitudinal study, it cannot be determined whether the observed functional connection pattern is preexisting to the continual migraine attacks, or whether the continual migraine attacks result in the abnormal functional connection pattern within the brain. Nonetheless, these discriminative features were tested on different data sets (i.e., the training/test cohort and validating cohort) by combining the bootstrap and permu-

CONFLICT OF INTEREST
The authors declare no potential conflict of interest.

DATA AVAILABILITY STATEMENT
The data and code that support the results of this study are available from the corresponding author upon request.