Unaltered intrinsic functional brain architecture in young women with primary dysmenorrhea

Primary dysmenorrhea (PDM), painful menstruation without organic causes, is the most prevalent gynecological problem in women of reproductive age. Dysmenorrhea later in life often co-occurs with many chronic functional pain disorders, and chronic functional pain disorders exhibit altered large-scale connectedness between distributed brain regions. It is unknown whether the young PDM females exhibit alterations in the global and local connectivity properties of brain functional networks. Fifty-seven otherwise healthy young PDM females and 62 age- and education-matched control females participated in the present resting-state functional magnetic resonance imaging study. We used graph theoretical network analysis to investigate the global and regional network metrics and modular structure of the resting-state brain functional networks in young PDM females. The functional network was constructed by the interregional functional connectivity among parcellated brain regions. The global and regional network metrics and modular structure of the resting-state brain functional networks were not altered in young PDM females at our detection threshold (medium to large effect size differences [Cohen’s d ≥ 0.52]). It is plausible that the absence of significant changes in the intrinsic functional brain architecture allows young PDM females to maintain normal psychosocial outcomes during the pain-free follicular phase.


Results
Demographic data and psychological assessments. There were no significant between-group differences regarding age (PDM: 23.1 ± 2.27 years of age, control: 23.7 ± 2.40 years of age, P = 0.147), age at menarche (PDM: 12.2 ± 1.19 years of age, control: 12.2 ± 1.11 years of age, P = 0.811), years of menstruating (PDM: 10.9 ± 2.53 years, control: 11.5 ± 2.69 years, P = 0.194), or average duration of one menstrual cycle (PDM: 29.3 ± 1.41 days, control: 29.5 ± 1.19 days, P = 0.525). The PDM group had a long history of menstrual pain (8.8 ± 2.75 years), with the pain lasting approximately 1 to 3 days during one menstrual cycle (2.0 ± 0.84 days). The current menstrual pain experience, as assessed by the pain rating index (29.3 ± 12.70) and present pain intensity (2.6 ± 1.01) of McGill Pain Questionnaire, confirmed that the PDM group experienced moderate to severe menstrual pain. The PDM females reported significantly higher state anxiety, trait anxiety, Beck Anxiety Inventory, and Pain Catastrophizing Scale scores during both the menstrual phase (MENS phase) and periovulatory phase (POV phase) ( Table 1).
Serum gonadal hormone measurements. No significant between-group differences were found in the concentrations of estradiol, progesterone, and testosterone during either the MENS phase or POV phase (Table 1).
Global network metrics of the resting-state functional network. Since the stringent statistical correction for multiple comparisons revealed negative findings, we therefore adopted a liberal uncorrected approach for the following analyses in order to unravel possible subtle effects. For the respective weighted and binary networks (covariates of gonadal hormones and psychological status adjusted), no main effect of group, menstrual cycle phase, or interaction between them was noted for any of the global network metrics (mean clustering coefficient, characteristic path length, global efficiency, and local efficiency of the network) with cost values ranging from 0.03 to 0.40 (all P > 0.05, uncorrected for multiple comparisons of 38 cost levels) ( Fig. 1: weighted network; Fig. 2: binary network). For the definitions of mean clustering coefficient, characteristic path length, global efficiency, and local efficiency of the network, please refer to the Methods section (see Global network metrics of the resting-state functional network). Small-worldness of brain functional network was confirmed for each phase of both the PDM and control groups. The general efficiency and small-worldness property of the resting-state brain functional network were not altered by acute menstrual pain in young PDM females.
SCIeNTIfIC REPORTs | (2018) 8:12971 | DOI: 10.1038/s41598-018-30827-6 Modular structure. For the full range of network costs (0.03-0.40) (covariates of gonadal hormones and psychological status adjusted), no main effect of group, menstrual cycle phase, or interaction between them was noted for the modularity and the number of partitioned modules in the respective weighted and binary networks (all P > 0.05, uncorrected for multiple comparisons of 38 cost levels) (Fig. 3: weighted network; Fig. 4: binary network). By assessing normalized mutual information (NMI) at a specific cost level, no significant differences were found in the similarity of modular partitions for the between-group comparisons during each of the MENS and POV phases (Fig. 5) (all P > 0.05, uncorrected for multiple comparisons of 38 cost levels). Moreover, no significant differences were found in the similarity of modular partitions for the between-phase comparisons in the respective PDM and control groups (Fig. 6) (all P > 0.05, uncorrected for multiple comparisons of 38 cost levels). For the definitions of modularity and NMI, please refer to the Methods section (see Modular structure). The modular structures of the PDM and control groups were found to be similar within the same respective phase.

Modular assignment of specific nodes.
For the full range of network costs (0.03-0.40) in the respective weighted and binary networks, no significant differences were found in the modular assignment of specific nodes for the between-group comparisons during each of the MENS and POV phases (all P > 0.05 for 90  . For the definition of node degree, please refer to the Methods section (see Regional network metrics of the resting-state functional network). The local connectivity property of brain functional networks in terms of node-level connectivity analyses was not altered by acute menstrual pain in young PDM females.

Discussion
In this study, we examined the global and regional network metrics and modular structure of brain functional networks in young PDM females using networks constructed by parcellated cerebral regions according to the AAL and Harvard-Oxford cortical and subcortical probabilistic atlases. No significant between-group differences were found for the metrics of global and local network efficiency of information transfer between nodes. It indicates that the young PDM females may retain the integrity of global and local connectivity properties of brain functional networks despite the presence of maladaptive neuroplasticity of DPMS 11 and DMN 10,12,13 . Although the young PDM females in our study exhibited higher levels of depressive mood and anxiety compared to the controls (Table 1), these levels did not reach the degree of clinical severity, such as an anxiety disorder that is mandatory for psychiatric intervention. The average Beck Depression Inventory score during the MENS phase in the PDM group was conventionally interpreted as minimal depression 32 , and the average Beck Anxiety Inventory score as mild anxiety 33 . of the binary network, no significant differences were found in the global network metrics among the PDM MENS , CON MENS , PDM POV , and CON POV by conducting linear mixed models. C p , clustering coefficient; L p , characteristic path length; CON MENS , menstrual phase of the control group; CON POV , periovulatory phase of the control group; PDM MENS , menstrual phase of the primary dysmenorrhea group; PDM POV , periovulatory phase of the primary dysmenorrhea group.
In terms of unaltered global and regional network metrics and modular structure of brain functional networks in young PDM females, the negative findings support that PDM is essentially dissimilar to major neurological 34 and psychiatric diseases [35][36][37][38] . This is particularly important because otherwise healthy PDM females usually do not exhibit overt cognitive, affective and psychosocial liability and disability, despite the fact that the PDM females may exhibit minimal or mild depressive and anxious symptoms during the MENS phase 3,8 . In this study, we deliberately excluded the females with severe premenstrual syndrome or suspected premenstrual dysphoric disorder, which features a combination of cyclical cognitive and affective debilitations and physical symptoms and is associated with significant functional impairments for daily living 39,40 . Premenstrual dysphoric  disorder has been regarded as a psychiatric disease 40 , and psychiatric diseases may exhibit disrupted topological organization of large-scale structural and functional brain networks 36,41 . The absence of overt psychological and psychiatric disturbance is succinctly reflected by the unaltered global and local connectivity properties of brain functional networks in young PDM females. It is documented that the small-world architecture and optimal topological organization of human brain networks are prominently altered in depression 36,37 , social anxiety disorder 38 , Alzheimer's disease 34,[42][43][44] , and schizophrenia 35,41,45,46 , and the graph metrics can be correlated with the symptom severity of diseases 34,37,42 . These neurological and psychiatric diseases not only target distributed brain regions but also disrupt the topological configuration of large-scale brain networks 24,26 .
Although previous studies reported alterations in local connectivity of the brain in young PDM females by means of regional homogeneity 12,47 and seed-based functional connectivity analyses [10][11][12][13] , the current study showed unaltered local connectivity property in the context of regional network metrics by using graph theoretical network analysis. The discrepancies may be inherent in the different methodologies used. Graph theoretical network analysis quantifies the functional connectivity profiles among node neighborhoods in the topological space and provides the information of connection pattern on a whole-brain scale 24 . Seed-based approach compares the region-to-region functional connectivity in the network, and can only provide the information that is limited to the connection pattern of the selected seed region 24 . Regional homogeneity analysis measures the extent of functional synchronizations among the neighboring voxels and provides the information of connection pattern on a local spatial scale 48 .
The abnormal functional connectivity and networking of the brain in chronic pain disorders cast light on the neurophysiological mechanisms underlying pain chronification and the different facets of associated neurocognitive conditions 18,[49][50][51][52] . We speculate that the alterations in functional connectivity in young and otherwise healthy PDM females can be initially limited to certain systems of the brain (e.g., maladaptive functionality of the DPMS 11 and the concomitant adaptive reorganization of cross-network connectivities 12 ), without the disruption of topology of brain functional networks. Dysmenorrhea can be associated with high rates of comorbidity with chronic functional pain disorders later in life 1 , and chronic functional pain disorders prominently exhibit not only hypo-connectivity of the DPMS with the DMN 18,51,53 but also alteration of topology in the whole-brain networking [18][19][20][21][22] . In addition, the extent of alterations of functional connectivity in migraine, a disease that is often comorbid with dysmenorrhea, proceeds from local systems to topological organization of whole-brain networking as the disease duration increased 22 . Chronification of pain involves large-scale, dynamic reorganization of brain connectivities since the transition from acute to chronic pain parallels the transition of brain activation from nociceptive to emotional and reward circuits 54 . We have also recently reported that genetic factors, such as the BDNF Val66Met 9 and OPRM1 A118G 55 polymorphisms, may affect, in the genotype-specific process, the functional connectivity dynamics of the DPMS in young PDM females. Such brain resilience that vary with the genotypes contributes to individual pain differences and susceptibility, which in turn may have an influence on the vulnerability for the subsequent development of chronic pain disorders. It is tempting to propose that the cumulative aberrancy of the dysfunctional DPMS 11 and the repetitive overstimulation of menstrual pain may eventually lead to the alteration of topology of brain architecture in later life of vulnerable PDM females through the chronification processes and genetic attributions. Longitudinal follow-up is needed to investigate the possibly progressive alterations in functional brain architecture of PDM females, and therefore to improve the understanding of the co-occurrence of chronic functional pain disorders.
The correlations between the regional structural changes in gray matter volume and topology of resting-state brain functional networks in PDM females are of interest for further investigation. One may speculate that the patterns of regional structural and functional alterations in the brain shown in previous PDM studies 5-7,10-13,56 represent a disordered configuration and will translate into abnormal interregional functional connectivity in terms of large-scale networking. The cingulate cortex, a brain region integrating pain, negative emotions and cognition 57 , exhibits trait-related changes in glucose metabolism and state-and trait-related changes in gray matter volume in PDM females [5][6][7] . Medial brain motor areas, including the supplementary motor area, exhibit altered functional connectivity in chronic pain conditions that are often comorbid with dysmenorrhea, including painful bladder syndrome 58 and localized provoked vulvodynia 59 . However, the modular assignments and regional network metrics of anterior and posterior cingulate gyri and supplementary motor area parsed based on the AAL and Harvard-Oxford cortical and subcortical probabilistic atlases were not altered in this study. The interactions between functional networks (functional connectivity) and their structural neural substrates of PDM remain to be explored, although functional and structural brain networks are closely related and may share common topological characteristics 24,60,61 . Different brain imaging modalities assess distinct but complementary types of brain connectivity. Integrated models of structural and functional connectivity constructed from multimodal brain imaging data can further shed light on the effect of PDM on topological network architecture.
By using power analysis (G*Power software; http://www.gpower.hhu.de/), the sample size of our study (n = 119 in total and ~60 in each group) may achieve statistical power (80%) to detect the medium effect size larger than 0.52 for the between-group/within-phase comparison (independent t test) and 0.37 for the within-group/between-phase comparison (paired t test). However, statistical significance in our study cannot be achieved should the effect size is smaller (e.g., underpowered concern in our statistical construct). Possible existence of more subtle changes cannot be precluded. Future studies of larger sample sizes are needed to detect possible statistical differences under smaller effect size. We had investigated the global network metrics and modular structure across a wide range of cost levels, using liberal, uncorrected statistical threshold for 38 cost levels to better detect the differences. Thus, it can be reasonable to conclude that there would be no significant differences in the graph metrics between the PDM and control groups.
Several points for further consideration should be addressed in our study. Firstly, the different atlas parcellation schemes may influence the results of topological network characteristics [62][63][64] , and the parcellation scheme of the cerebral cortex for graph theoretical network analysis requires further optimization and standardization 27,65,66 . The AAL and Harvard-Oxford cortical and subcortical probabilistic atlases are the parcellation schemes that are commonly used in graph theoretical network analysis of fMRI studies 36,67 . Secondly, it should be noted that the functional connectivity network in this study was based on the 90 cerebral regions of the AAL atlas, excluding the cerebellar regions. Although our previous PDM study showed trait-related changes of gray matter volume in the right cerebellar tonsil 6 , the specific role of the cerebellum in pain processing and modulation is not completely understood 68 . Notwithstanding, technical challenges for image preprocessing of the cerebellum (e.g., normalization procedure to the atlas) arise due to the small size and functional heterogeneity of the cerebellum 69 . Finally, the young PDM females in our study group experienced moderate to severe menstrual pain with a wide spectrum of individual pain differences. Owing to the sample size, we did not perform subgroup analyses according to the pain severity (i.e., moderate and severe PDM). Future studies of larger sample sizes to address the relationship between the topological profiles and the clinical characteristics (e.g., pain severity) can be important to disclose the intricacy between imaging endophenotypes and clinical phenotypes of PDM.

Conclusions
Our current study indicates that, despite regional dysfunction and structural alterations, the overall integrity of intrinsic functional brain architecture in young PDM females, in terms of global and regional network metrics and modular structure, showed no differences compared to females without PDM at our detection thresholds (Cohen's d between-group ≥ 0.52 and Cohen's d within-group ≥ 0.37). Although we cannot completely preclude small to medium effect size differences of functional brain architecture owing to the underpowered concern, it is plausible that the absence of significant changes in the functional brain architecture allows young PDM females to maintain normal psychosocial outcomes during the pain-free follicular phase.

Methods
Subjects. The subjects of this study (59 PDM and 68 control females) were a subset of the participants (smaller in the sample size) from our previous genetic association and behavioral study of PDM 8 who had completed the whole study protocols and were eligible for neuroimaging analyses in this study. Two PDM and 6 control females were excluded for further analyses owing to the significant head motion (translation >2 mm or rotation >2°) during the MRI scan. Eventually, 57 otherwise healthy females with PDM (age, 23.1 ± 2.27 years) and 62 education-matched, healthy control females (age, 23.7 ± 2.40 years) were recruited for the present study (see Table 1 for demographic data). The inclusion criteria for subjects with PDM were as follows: 1) 20-30-year-old Taiwanese (Asian) females; 2) a regular menstrual cycle of approximately 27-32 days; 3) a history of menstrual pain longer than 6 months; 4) average menstrual pain under regular treatment with a rating higher than 4 on a verbal numerical scale (VNS, 0 = not at all, 10 = the worst imaginable pain) in the last 6 months; and 5) right-handedness, as confirmed by the Edinburgh Handedness Inventory 70 . All PDM females were clinically examined and diagnosed in the gynecology clinic by the same gynecologist (H.T.C.). All subjects in the PDM group received pelvic ultrasonography to exclude secondary dysmenorrhea caused by organic pelvic diseases, such as endometriosis or adenomyosis. The inclusion criteria for the healthy control females were similar to those for the PDM group, except that the subjects of the control group had no pain whatsoever during menses (VNS = 0). The exclusion criteria for all of the subjects were as follows: 1) using oral contraceptives, hormonal supplements, Chinese herbal medicine, or any centrally acting medication (e.g., opioid, anti-epileptics) within 6 months prior to the study; 2) pathological pituitary gland disease; 3) organic pelvic disease; 4) any psychiatric or neurological disorders (e.g., premenstrual dysphoric disorder); 5) any head injury with loss of consciousness or brain surgery; 6) immediate plans for pregnancy or a positive pregnancy test; 7) a history of childbirth; and 8) having a metal/pacemaker implant, claustrophobia, or any contraindications in relation to MRI. No analgesics had been taken by the subjects within 24 hours before the study. Some of these subjects have been analyzed in the previously published neuroimaging studies 9,11,12 . The overall study was conducted in accordance with the Declaration of Helsinki and was approved by the Institutional Review Board of Taipei Veterans General Hospital. All subjects provided written informed consent. Experimental design. MRI scanning was individually scheduled according to the first day of menstruation for each subject. All of the subjects received psychological assessments, blood sampling for gonadal hormone assays, and brain MRI scanning (T1 and resting-state fMRI images) during the MENS phase (days 1-3 of the menstrual cycle) and POV phase (days 12-16 of the menstrual cycle). Ovulation was confirmed using a urinary luteinizing hormone test (Han Chiun Proper LH Rapid Test).

Psychological assessments.
To evaluate the psychological status throughout the menstrual cycle, all of the subjects in the two groups completed self-reported psychological measurements, including the Spielberger State-Trait Anxiety Inventory 71 , the Beck Anxiety Inventory 72 , the Beck Depression Inventory 73 , and the Pain Catastrophizing Scale 74 , during the respective MENS and POV phases. The McGill Pain Questionnaire was completed by the PDM females at the inception stage and the MENS phase during the experiment to assess their respective overall and current experiences of menstrual pain.
Serum gonadal hormone measurements. The sera extracted from the blood samples drawn during the respective MENS and POV phases were stored for batch analysis using commercialized assays (UniCel DxC 800 Synchron Clinical Systems, Beckman Coulter, Inc., Brea, CA). The total serum concentrations were assayed using a chemiluminescence immunoassay technique for estradiol and progesterone and a radioimmunoassay technique for testosterone.
Our previous report revealed: 1) there were significant main effects of group and menstrual cycle phase as well as the interaction between them on the psychological measurements and 2) there was a significant main effect of menstrual cycle phase, but no main effect of group or interaction between group and menstrual cycle phase, on the serum gonadal hormone levels 8 . For the purpose of this neuroimaging subset study, we only conducted the two-sample t-test to examine the between-group differences of the psychological and gonadal hormone measurements in each phase using SPSS Statistics 20.0 (SPSS Inc., Chicago, IL), without testing the overall interaction effects between group and menstrual cycle phase. This aim was to ascertain that the imaging data was not confounded by hormonal differences between the groups. The results were considered significant when P < 0.05. For more comprehensive statistical analyses and results of psychological and gonadal hormone measurements, please refer to our previous report 8 . Image preprocessing. All EPI images were preprocessed using Statistical Parametric Mapping software (SPM8, Wellcome Trust Centre for Neuroimaging, University College London, London, United Kingdom, http:// www.fil.ion.ucl.ac.uk/spm) in MATLAB (The MathWorks, Inc., Natick, MA, USA) with the following procedures: the correction of slice timing, the realignment for head motion correction (6-parameter rigid body transformation), and the spatial normalization. The time course of head motion for each subject was obtained by estimating the translation and rotation in each axis for the 200 consecutive EPI volumes. Because head motion has a significant influence on intrinsic functional connectivity measurements 75 , we excluded subjects with significant head motion (translation >2 mm or rotation >2°) of any volume from further analysis. The EPI images were spatially normalized using the SPM's standard EPI template in Montreal Neurological Institute (MNI) space and re-sampled to an isotropic 2 × 2 × 2 mm 3 voxel size. For the time series of blood oxygen level-dependent (BOLD) signals in each voxel of the normalized image, the effects of head motion (the 6 motion parameters estimated from rigid-body realignment) and the signals of white matter regions and ventricular system were removed by linear regression, and the shift of BOLD signals was detrended. The band-pass filter was set at 0.01-0.08 Hz for the resulting time series to extract the low-frequency oscillating components that contributed to resting-state functional connectivity 76 . Global signal regression was not performed due to the concern of potential distortion on correlation patterns 77,78 . Construction of connectivity network. Brain functional networks, composed of nodes (parcellated brain regions) and edges (inter-nodal functional connectivity), were constructed from resting-state fMRIs.

Image acquisition.
Definition of node. The node of the network was defined as the parcellated brain regions according to the AAL atlas (consisting of 90 cerebral and 26 cerebellar anatomical regions) 31 . For the nomenclature of the 90 parcellated cerebral regions (45 areas per hemisphere) used in network construction, please refer to our previous report 79 . The resulting templates of the 90 nodes were then respectively co-registered with the preprocessed fMRI images of each subject. The mean time series of each node was calculated by averaging the time series of each voxel in that node. Definition of edge. The edge of the network was defined as the degree of correlation (Pearson correlation coefficient, r) between the mean time series of each pair of nodes 22,24,26,27,79 . To obtain a better normality of the correlation, the 90-by-90 'r-value matrix' was converted to a 'z-score matrix' using Fisher's r-to-z transformation. The weight of the edge was defined as the absolute value of the z-score matrix.
The z-score matrix was then thresholded at different levels of sparsity using the minimum spanning tree method, followed by global thresholding 80 . The minimum spanning tree method aims to force the connectedness of sparse or low-density graphs and to prevent the emergence of fragmented networks with isolated islands, which are neuroscientifically unexplainable and unjustified. For each node, the edge of the highest weight was retained to keep the node connected with at least one neighboring node (i.e., local thresholding). The step of local thresholding was iterated until all nodes were able to connect with one another, giving rise to a fully connected backbone network without the emergence of isolated islands. The edges of the strongest weight, ranked by all the edges, were sequentially added into the backbone network. This 'growing' method was iterated until the number of edges meets with the assigned network sparsity (i.e., global thresholding). To investigate the effect of different network sparsity on graph metrics, weighted and binary networks were constructed with cost values ranging from 0.03 to 0.40, at increments of 0.01. For a weighted network, the edge weight exceeding the threshold was retained, otherwise it was discarded. For a binary network, the edge weight exceeding the threshold was set to '1' , otherwise it was set to '0' . Binary networks contain only the information whether an edge (connection) between two nodes is present or not according to the preset threshold. Weighted networks contain not only the information whether an edge is present or not, but also the connection strength between two nodes. Binary networks are succinct to characterize and compute the organizational characteristics 29 . Weighted networks can, compared to the binary networks, retain more biologically relevant information and provide different but complementary information of network organization 29,81 . We used both binary and weighted network constructions to cross-validate the results.
Network analysis software. The global and regional network metrics and modularity of brain functional networks were computed using the Brain Connectivity Toolbox 29 (https://sites.google.com/site/bctnet/) in MATLAB. Between-group differences of network community structures (i.e., similarity analysis), in terms of modular partitions and modular assignment of specific nodes, were examined according to the approach reported by Alexander-Bloch et al. 82 . We adopted the R codes provided by the authors and computed the NMI and Pearson's phi for the similarity analyses of modular partitions and modular assignment of specific nodes, respectively. For the codes to perform the similarity analyses please refer to the Supplementary information. The similarity comparisons were carried out by means of permutation procedure for the planned contrast. We compared the average within-group similarity of the actual data to that of the permuted data where the group memberships are randomized. The P value is defined as the number of instances that the permuted within-group similarity exceeds the actual within-group similarity, divided by the number of permutations. Significant group differences are indicated should the within-group similarity of the actual data always exceed that of the permuted data 82 . For the details of computational processes, please refer to the respective sections (Modular structure & Modular assignment of specific nodes).
Global network metrics of the resting-state functional network. The mean clustering coefficient 83,84 , characteristic path length 83 , global efficiency 85 , and local efficiency 85 of the network were assessed for each of the MENS and POV phases in the PDM and control groups. The mean clustering coefficient of the network represents the tendency of network nodes to form local clusters (i.e., local connectedness), and the characteristic path length indicates how well the network nodes are interconnected (global connectedness) 34 . Small-worldness reflects the combination of a high clustering coefficient and a short characteristic path length 83 . A high clustering coefficient characterizes an efficient local information processing (high local efficiency) and a short characteristic path length features a high level of global communication efficiency (high global efficiency) 28,86 . Modular structure. The modularity 87,88 , the number of partitioned modules, as well as the similarity between two modular partitions were assessed. The modular structure is defined as a subset of nodes that are more densely interconnected among them than with other nodes within the network, indicating the formation of sub-networks 24 . Modularity is a measure of degree to which a network can be subdivided into smaller sub-networks 26 .
The global network metrics (mean clustering coefficient, characteristic path length, global efficiency, and local efficiency) and modular structure (modularity, number of partitioned modules, and similarity of modular partitions) of brain functional networks were assessed on the respective weighted and binary networks across a range of network cost (0.03-0.40). The statistical analyses for the mean clustering coefficient, characteristic path length, global efficiency, local efficiency, modularity, and number of partitioned modules of the network at a specific cost level were conducted using SPSS Statistics 20.0. Linear mixed models were used to examine the fixed effects of group (PDM vs. control) and menstrual cycle phase (MENS phase vs. POV phase), as well as the interaction between them. Subject was entered as a random effect. As menstrual cycle phase (i.e., gonadal hormones) may have an influence on the resting-state functional connectivity between brain regions 89 and PDM females exhibit minimal or mild depressive and anxious symptoms during the MENS phase (see Results), the hormonal fluctuations (estradiol, progesterone, and testosterone levels) and psychological assessments (the Spielberger State-Trait Anxiety Inventory, the Beck Anxiety Inventory, and the Beck Depression Inventory scores) were entered as covariates in the statistical model. The variance-covariance matrix was assumed to be unstructured since there were no presumptions on the correlation pattern of graph metrics between menstrual cycle phases. Significance was initially set at a threshold of corrected P < 0.05 (FDR correction for multiple comparisons of 38 cost levels), and no significant differences were found. Hence, we further lowered the significance level at a liberal threshold of SCIeNTIfIC REPORTs | (2018) 8:12971 | DOI:10.1038/s41598-018-30827-6 uncorrected P < 0.05 (uncorrected for multiple comparisons of 38 cost levels) to better the detection of possible differences.
Similarity between two modular partitions (quantified by the NMI 82,90 ) was compared by a permutation procedure 82 . The NMI, a similarity index, can be used to evaluate different modular partition algorithms and to estimate the similarity of modular structures across subjects 90 . The NMI value ranges from 0 to 1, where 1 denotes that the two modular partitions are identical and 0 that they are totally independent 82 . For the equation of NMI, please refer to previous reports 82,90 . Between-group comparisons for each phase and between-phase comparisons for each group were computed to address the state and trait effects of dysmenorrhea on the similarity of modular structure of brain functional networks 11,12 . The alterations of modular partitions during the MENS phase in the between-group comparison (MENS phase of PDM vs. MENS phase of control) and in the between-phase comparison for the PDM group (MENS phase of PDM vs. POV phase of PDM) were regarded as state-related changes. The alterations during the POV phase in the between-group comparison (POV phase of PDM vs. POV phase of control) were regarded as trait-related changes. At a specific cost level, the average pairwise NMIs during the same menstrual cycle phase were separately calculated across subjects within the control and PDM groups. We tested the hypothesis that the average within-group pairwise similarity is higher than the average between-group pairwise similarity, which exhibits a genuine difference in similarity of modular partitions between the control and PDM groups during the same menstrual cycle phase 82,91 . The permutation procedure to compare the similarity of modular partitions was performed by varying the group membership, with 10000 permutations 82 . Significance was set at a stringent threshold of corrected P < 0.05 (FDR correction for multiple comparisons of 38 cost levels) initially, and no significant differences were found. Hence, we further lowered the significance level at a liberal threshold of uncorrected P < 0.05 (uncorrected for multiple comparisons of 38 cost levels) to unravel possible subtle differences.

Modular assignment of specific nodes.
For the analyses of modular assignment of specific nodes, we analyzed the respective weighted and binary networks with cost values ranging from 0.03 to 0.40. To investigate the difference in modular assignment of a specific node of interest (NOI) at a specific cost level during the same menstrual cycle phase, the similarity of module labels of two subjects, in terms of the network community of the NOI, was calculated. For a given NOI, for each subject, all of the other nodes were labeled '1' if they shared the same module with the NOI, and '0' if not. The similarity of the module labels between each pair of subjects was calculated using Pearson's phi, a method of Pearson correlation for dichotomous variables. For the given NOI, a higher phi value indicates a higher similarity of modular assignment between subjects. We tested the hypothesis that the average within-group pairwise similarity is higher than the average between-group pairwise similarity, which exhibits a genuine difference in similarity of modular assignment of specific nodes between the control and PDM groups during the same menstrual cycle phase 82 . The permutation procedure to compare the similarity was performed by varying the group membership, with 10000 permutations 82 . The permutation tests were performed for every one of the 90 AAL nodes at each of the 38 cost levels. Under the premise of equivalent statistical stringency for the global network metrics and modular assignment of specific nodes, FDR correction was only applied to correct for multiple comparisons of 90 AAL nodes, but not for multiple comparisons of 38 cost levels.
Regional network metrics of the resting-state functional network. For each node, the regional network metrics (local efficiency, clustering coefficient, and degree) were assessed on the respective weighted and binary networks across a range of network cost (0.03-0.10, at increments of 0.01) for each of the MENS and POV phases in the PDM and control groups. Low cost levels were selected for network construction to avoid spurious or statistically non-significant connections between nodes 28,29 . Node degree is defined as the number of links connected to a node 29 . The node with a high degree is interacting with many other nodes in the brain functional network 29 .
Data were analyzed using SPSS Statistics 20.0. At a specific cost level, linear mixed models were used to examine the fixed effects of group (PDM vs. control) and menstrual cycle phase (MENS phase vs. POV phase), as well as the interaction between them on the regional network metrics of each node. The hormonal fluctuations (estradiol, progesterone, and testosterone levels) and psychological assessments (the Spielberger State-Trait Anxiety Inventory, the Beck Anxiety Inventory, and the Beck Depression Inventory scores) were entered as covariates in the statistical model. The variance-covariance matrix was assumed to be unstructured. FDR correction was applied to correct for multiple comparisons of 90 AAL nodes, but not for multiple comparisons of 8 cost levels.

Robustness of methodological variation.
To validate our observations using results obtained from different methodological approaches, we analyzed the difference in the global and regional network metrics and modular structure using the networks constructed by parcellated brain regions according to the Harvard-Oxford cortical and subcortical probabilistic atlases. The atlases are composed of 48 cortical and 7 subcortical brain regions in each hemisphere, and a total of 110 nodes (55 areas per hemisphere) were used in the network construction.