The perception of acute pain is a product of neural activity in a widely distributed brain network, typically following noxious insult. In chronic pain conditions, however, altered processing within this network can create, amplify, or sustain the perception of pain independent of noxious stimulation.52 Therefore, new methods of modeling brain network activity may offer unique insights into the pathogenesis of chronic pain conditions.11,31
Brain networks can be modeled using graph theoretical tools as a set of functional interactions made up of nodes (brain regions) and edges (structural or functional connections between nodes). The organization, or topology, of a network is critical to its function because it influences the efficiency and content of information transfer among nodes.39,40
A critical component of brain network architecture is a robust hub structure, wherein hub nodes facilitate efficient information integration by occupying a highly connected and functionally central role in the network.60 Brain networks also have a higher order level of organization called “rich clubs,” in which hubs are more likely to be connected to each other than to nodes with fewer connections.7,19,59 This rich club forms a densely connected core, wherein hubs may act in concert and link different functional systems in the brain.58,59
Importantly, although hubs create efficiency, they also create network vulnerability wherein disruption of hub organization can have widespread consequences for information transfer. Hubs may be disproportionately affected in neurological disorders.8,56 A recent meta-analysis of 26 diverse brain disorders found that gray-matter alterations were more likely to be in hubs and rich-club hubs in particular.8 Alterations in optimal network structure have been reported in chronic pain disorders.2,32,36 However, how hub topology is altered in chronic pain, and its relationship to underlying neurochemistry and patient self-reported clinical pain is largely unknown. More broadly, it remains unknown how alterations in neurotransmitter levels in a clinical population may relate to changes in hub strength and rich-club membership.
Here, we used a graph theoretical approach in conjunction with resting state functional magnetic resonance imaging (fMRI) and proton magnetic resonance spectroscopy (1H-MRS) to examine global network properties, hub topology, and the relationship to clinical pain and neurochemistry in fibromyalgia (FM), a common chronic pain condition with well-established central nervous system contributions.4 In addition to widespread pain and hyperalgesia, FM patients display hypersensitivity to multiple modalities of sensory stimuli.17,25,33,37,46,65 Fibromyalgia patients consistently have altered levels of brain activity in regions related to pain and sensory processing,20,34 and differences in functional connectivity in pronociceptive and antinociceptive regions.27–29,43 Moreover, FM patients have altered levels of excitatory and inhibitory neurotransmitters within the insula and other sensory processing regions.13,23 We hypothesized that the FM brain network would have altered hub organization and rich-club membership of sensory processing regions (such as the insula), and that these differences would be related to excitatory neurochemical tone and clinical pain intensity.
2. Materials and methods
2.1. Overview of analyses
We first consider the hypothesis that global network properties are altered in chronic pain by deriving metrics of several global network properties and comparing them between the FM and healthy control (HC) cohorts. We then test group differences at the level of individual brain nodes by deriving a measure of “connectedness” to determine hub status and comparing the strength of these hubs between cohorts. We consider both global hub status and higher-order “rich-club hub” status (described below). To determine how clinical pain relates to hub metrics, we divide our FM sample into patients experiencing low, medium, and high levels of clinical pain at the time of the scan and compare these groups. Using the validation cohort, we attempt to replicate hub status of individual nodes, rich-club membership, and the association between hub metrics and clinical pain. Finally, we explore potential mediated relationships between high levels of Glx in the insula of FM patients and clinical pain ratings through measures of insula node centrality. See supplementary Table 1 for definitions of key terminology (available at http://links.lww.com/PAIN/A729) or Ref. 14 for detailed descriptions and mathematical formulations.
2.2.1. Discovery cohorts
Data from 40 female FM patients and 46 age- and sex-matched HCs were pooled from previous fMRI studies and analyzed retrospectively using a graph-theoretical approach. The data from 31 FM and 46 HCs were previously published in manuscripts focusing on functional connectivity, 1H-MRS, and treatment outcomes.22,23,25,28,42,43,49,51 Data from the remaining 9 FM patients have not been published. The University of Michigan Institutional Review Board approved the studies, and written informed consent was obtained from all participants in accordance with the Declaration of Helsinki. The major inclusion criteria for FM patients were meeting the American College of Rheumatology (ACR) 1990 criteria for FM and at least 6 months of self-reported chronic widespread pain, a score of ≥40 mm on a 100-mm pain Visual Analog Scale (VAS) at the time of consenting. Both FM and HCs were 18 to 75 years of age, female, right-handed, and capable of giving written informed consent. Major exclusion criteria for both groups were pregnant or nursing mothers, contraindications to fMRI, positive urine drug screen or history of drug or alcohol abuse within the past 2 years, body mass index greater than 36, severe psychiatric illness, concurrent autoimmune or inflammatory disease that causes pain, or systemic malignancy or infection such as HIV or hepatitis. An additional exclusion criterion for HCs was meeting the ACR 1990 criteria for FM.
2.2.2. Validation fibromyalgia cohort
A separate data set consisting of 11 female FM patients was used as a validation cohort for graph theoretical metrics of resting-state fMRI analyses. The data from these patients were included in a previously published article examining changes in resting state functional connectivity after acute pain stimulation;27 however, a cross-sectional analysis of resting state connectivity compared with HCs in that sample has not been published. Inclusion criteria for the validation cohort were having met the 1990 ACR criteria for FM, disease duration of at least 1 year, reporting continued pain presence for more than 50% of each day, willingness to forgo new medications or treatments for FM during the study, between 18 and 75 years old, right-handed, and capable of giving written informed consent. Exclusion criteria were current or past use of opioid or narcotic analgesics, history of substance abuse, concurrent autoimmune or inflammatory disease that may contribute to pain, concurrent participation in other therapeutic trials, or a severe psychiatric illness.
2.3. Demographics and clinical assessments
A standardized form was used to collect age, race, ethnicity, and current medications (9 FM patients were missing medication data, and 12 HCs were missing race and ethnicity data). Self-reported clinical pain intensity in FM patients was assessed immediately before fMRI using a VAS, with 0 being “no pain” and 10 being the “worst pain imaginable.” Depression and anxiety were assessed in a subset of FM patients (n = 28) and HCs (n = 32) using the HADS.67
2.4. Functional magnetic resonance imaging data acquisition
All participants (discovery and validation cohorts) completed a resting state fMRI scan, and 40 FM patients from the discovery cohort and 27 HCs completed a 1H-MRS sequence of the right posterior insula. Scanning was performed for all subjects on a 3.0 T General Electric (Milwaukee, WI) system with an 8-channel head coil. The MRI scanner at the University of Michigan was upgraded to a new 3.0 T GE system with identical scanning parameters during one of the studies reported on here, so data from 3 FM participants were acquired on the new scanner. Group differences in graph theoretical measures did not change when excluding these participants (data not shown).
Each participant completed a 6-minute eyes-open resting state fMRI scan as described previously.22,42,43 Data were collected using a spiral in-out gradient echo T2*-weighted blood oxygenation level–dependent pulse sequence with the following parameters: repetition time = 2000 ms, echo time = 30 ms, 180 volumes, 43 slices parallel to the anterior-posterior commissure plane, and voxel size 3.13 × 3.13 × 4.0 mm. The first 6 volumes were discarded to avoid equilibration effects. A high-resolution structural image was collected for registration purposes (spoiled gradient echo pulse sequence: repetition time = 14 ms, echo time = 5.5 ms, inversion time = 300 ms, 20° flip angle, 124 contiguous axial slices, and voxel size 1 × 1 × 1.5 mm). Cardiac and respiratory data were collected simultaneously using an infrared pulse oximeter (GE) attached to the right middle finger and a respiration belt placed around the participant's ribcage.
2.5. Functional magnetic resonance imaging data preprocessing
Data from the discovery and validation cohorts were preprocessed using identical pipelines. All data were checked for artifacts, and motion was assessed by examining 3 translations (x, y, and z) and 3 rotations (pitch, roll, and yaw). Participants were excluded if the mean absolute motion exceeded ±2 mm translation or ±1° rotation in any direction. To check for differences in motion between groups, we calculated the framewise displacement, a combined index of the 6 motion parameters, as described previously.47 Fibromyalgia and HCs did not differ according to framewise displacement (mean FM: 0.066 ± 0.04; mean HC: 0.059 ± 0.049; t = −0.768, P = 0.444). Resting state data were preprocessed using FSL (http://www.fmrib.ox.ac.uk/fsl) and SPM8 (http://www.fil.ion.ucl.ac.uk/spm/software) running on MATLAB R2014a (http://www.mathworks.co.uk/products/matlab/). Physiological correction was performed using the RETROICOR algorithm18 within the Physiological Noise Modelling toolbox in FSL. All other preprocessing steps were performed in SPM8 and included slicetiming, motion correction, coregistration of the structural and functional images, normalization to Montreal Neurological Institute space, and smoothing with an 8-mm FWHM Gaussian Kernel. Preprocessed data were entered into the Conn Toolbox,64 and a nuisance regression using the aCompCor method3 was performed with 6 subject-specific realignment parameters, the signal from white matter and CSF, and their first-order derivatives included as confounds. Finally, a temporal filter of 0.008 to 0.09 Hz was applied to focus on low-frequency fluctuations.15
2.6. Proton magnetic resonance spectroscopy data acquisition and preprocessing
During 1H-MRS, the spectroscopic voxel was placed in the right posterior insula (Supplementary Figure 1A, available at http://links.lww.com/PAIN/A729) as described previously.23,24 Single-voxel point–resolved spectroscopy spectra were acquired while participants were at rest with the following parameters: repetition time = 3000 ms, echo time = 30 ms, 90-degree flip angle, number of excitations 8, and number of averages 128, with a voxel size of 2 × 2 × 3 cm3. Shimming was optimized using auto-prescan, and the CHEmical Shift Selective (CHESS) water suppression routine was used. Spectra were analyzed offline with LCModel (Stephan Provencher, Oakville, ON, Canada).48 Raw spectra were fit with a linear combination of pure metabolite spectra using LCModel (the basis set, which was simulated by Stephan Provencher of LCModel, included the following metabolites: Ala, Asp, Cr, PCr, GABA, Glc, Gln, Glu, GPC, PCh, Ins, Lac, NAA, NAAG, Scyllo, and Tau). See Supplementary Figure 1B for a representative spectrum (available at http://links.lww.com/PAIN/A729). We have previously shown that FM patients have increased levels of glutamate + glutamine (Glx) in the posterior insula,23 and therefore, we focused on these metabolites. Glx values were calculated as concentrations rescaled using the water peak. Concentrations were corrected for CSF in each participant using voxel-based morphometry within SPM8 as previously described.23 To rule out the possibility of tissue volume biases the 1H-MRS results, we analyzed the segmented gray-/white-matter and CSF content of the spectroscopic voxel and found that it was not significantly different between groups (Supplementary Tables 2 and 3, available at http://links.lww.com/PAIN/A729). Spectra were excluded if the Glx Cramer-Rao bounds exceeded 20% or if the FWHM was greater than 0.12. One FM participant with poor quality spectra was excluded from analysis.
2.7. Graph theoretical analyses
The analysis flow is depicted in Figure 1. We defined the brain network using a set of 264 nonoverlapping nodes based on resting state and task functional connectivity meta-analyses.5,30 This set of nodes has been shown to produce reliable network topologies.6,10,54,62 The 264 nodes were entered into the Conn Toolbox as 10-mm diameter spheres30 to create Fisher z-transformed bivariate correlation (Pearson's r) matrices (264 × 264) for each participant. To exclude weak or spurious connections, matrices were thresholded beginning with connections in the top 5% and proceeding in steps of 5% up to 40% density, resulting in binary undirected graphs containing the most significant edges. Using the Brain Connectivity Toolbox,50 we calculated the following graph theoretical measures to assess global network properties: global efficiency, modularity Q score, clustering coefficient, characteristic path length, and the rich-club coefficient.
There are many graph theoretical measures used to assess hub status, many of which are based on degree (the number of connections a node has) or centrality (the relative importance of a node, or its influence on other nodes, based on network organization).14 Eigenvector centrality takes into account the connectedness of a node, in addition to the connectedness of that node's neighbors. Because this measure is more qualitative than degree, and because there is no agreed upon standard for categorizing hub nodes, we chose eigenvector centrality as our principle measure for defining hub nodes. We also calculated degree and present those results in the supplementary material (available at http://links.lww.com/PAIN/A729). Degree and eigenvector centrality were calculated using the Brain Connectivity Toolbox.50 See Ref. 14 for a detailed description and mathematical formulations of graph theoretical measures. To reduce the number of comparisons, for each metric and each node, we averaged across thresholds as previously published.1,35
2.8. Identification of hubs
We assigned hub status to a node if the eigenvector centrality was greater than 1 SD above the group mean.53,59 Analyses of between group differences on eigenvector centrality was restricted to hub regions and were examined using nonparametric permutation testing14,44 under the null hypothesis that hub nodes in FM and HCs do not differ with respect to eigenvector centrality. We randomly reassigned participants to 1 of 2 groups and calculated a 2-sample t-statistic. We repeated this procedure 10,000 times to form a randomized null distribution for each metric. We rejected the null hypothesis if the actual t-statistic was greater than or equal to the 95th percentile of the null distribution. Significance was set at P < 0.05.
2.9. Rich-club organization
We assessed rich-club organization in the average FM and HC brain networks by computing the rich-club coefficients ϕ (k) over a range of degree (k), as previously described.59 Briefly, the degree of each node was calculated for the FM and HC networks, and all nodes that had a degree ≤k were ignored. For the remaining nodes in the network, we calculated the rich-club coefficient ϕ (k) as the number of connections between the remaining nodes divided by the total number of possible connections. Random networks can also have an increasing function of ϕ (k) by chance alone; hence, the ϕ (k) is typically normalized (ϕnorm (k)) by a set of random networks. We created 1000 random networks with similar degree distribution and density, and for each level of k, calculated the average rich-club coefficient ϕrandom. We determined statistical significance of ϕ (k) at each level of k with permutation testing. Using the random networks created above, we obtained a null distribution of ϕrandom (k) values. P-values were estimated as the proportion of ϕrandom (k) that exceeded the observed ϕ (k) for FM and HC separately. This was repeated over all levels of k, and a Bonferroni correction (0.05/28) was applied, so that P < 0.001 was deemed significant at 5% density. The range of k in which ϕ (k) is significantly different from ϕrandom (k), and where the ϕnorm (k) is greater than one, is the rich-club regime. Differences between FM and HC in ϕ (k) at each level of k were tested in SPSS 24 (Armonk, NY) using independent-samples t tests. To assess how rich-club membership varied with clinical pain, we separated the FM group into age-matched tertiles (high n = 12, medium n = 16, low n = 12) based on the VAS clinical pain rating. Differences in ϕ (k) between the tertiles were examined using a 1-way analysis of variance. A Bonferroni correction was applied (0.05/28) and significance set at P < 0.001. For visualization purposes, rich-club nodes were displayed at k = 22 for FM and k = 25 for HC, the highest level of k that was significantly different from random networks for the FM and HC groups, respectively.
2.10. Examining the relationship between hub status, clinical pain, and Glx
For correlations between eigenvector centrality, clinical pain, and Glx, we performed Pearson correlations in SPSS. To test for differences in clinical variables and Glx between the FM tertiles, we performed 1-way analyses of variance in SPSS. All analyses controlled for age and significance were set at P < 0.05. To test the hypothesis that higher levels of posterior insular Glx influence clinical pain indirectly through greater eigenvector centrality in the insula, we conducted mediation analyses using MPLUS v. 8. Posterior insula Glx was used as the independent variable, and clinical pain on the VAS was used as the dependent variable. Eigenvector centrality values from 2 nodes overlapping the 1H-MRS voxel for the posterior insula (nodes 43 and 67) were standardized and averaged, and this value was used as the mediating variable. Indirect effects were evaluated by constructing 95% bias-corrected bootstrapped confidence intervals (CIs) using 20,000 resamples. All models controlled for participant age.
2.11. Validation analyses
Connectivity matrices for the FM validation cohort were created and thresholded using identical procedures as described above. We attempted to replicate the principal findings from the discovery cohort, specifically determining the hub status of sensory processing regions, group differences in eigenvector centrality of hub nodes, rich-club membership and variation with clinical pain, and finally, the association between eigenvector centrality and clinical pain. The FM validation cohort did not have 1H-MRS data, and therefore, replication of these findings could not be attempted. The FM validation cohort was compared with all 46 HC participants in the discovery data set. Because of the small sample size of the validation cohort, rich-club membership in relation to clinical pain was examined using a median split to divide patients into high and low clinical pain groups. All other analyses and statistical testing were performed in an identical manner to the discovery cohorts.
The visualization of 264 nodes on a brain surface was created in Caret.9,61 All other brain surfaces were created using BrainNet viewer.66
3.1. Clinical characteristics
There was no significant difference in age between FM and HC participants (P > 0.9; Table 1). Fibromyalgia patients reported significantly higher clinical pain (VAS) (P < 0.001) and were significantly more depressed and anxious (P < 0.001) relative to HCs. A list of current medications in the FM patients is listed in Supplementary Table 4 (available at http://links.lww.com/PAIN/A729).
3.2. Fibromyalgia and control participants had similar global network properties
Our initial aim was to characterize the global brain network properties of FM patients and controls. At every density of functional connections tested (5%-40%), FM and HC networks had small-world organization, defined as high clustering (local connectivity) and a low average path length between nodes. There were no significant differences between groups in the global efficiency, clustering coefficient, average path length, or modularity of the reconstructed brain networks (all P > 0.4, Supplementary Table 5, available at http://links.lww.com/PAIN/A729). As these measures did not differ significantly between groups, and because there is no accepted density for defining nontrivial connections, we averaged metrics across all densities (5%-40%, as previously published1,35).
3.3. Differences in hub status between fibromyalgia patients and controls
To investigate the nature or quality of edges made in the reconstructed brain network, we calculated the eigenvector centrality of each hub node. The eigenvector centrality accounts for the quantity and quality of connections by taking into account the degree of a node and the degree of that node's neighbors. We assigned hub status to nodes whose eigenvector centrality was greater than 1 SD above the mean, for each group separately (Fig. 2). The bilateral anterior insulae, mid insulae, primary motor cortex (M1), primary somatosensory cortex (S1), supplementary motor area (SMA), and superior temporal gyrus (STG) were hubs in FM. Hubs in HCs were more dominant in frontal and posterior regions, including the superior frontal gyrus (SFG) and primary, secondary, and visual association cortices (V1, V2, and V3). All hubs for each group are listed in Supplementary Table 6 (available at http://links.lww.com/PAIN/A729). To determine exactly which hubs were altered, we determined group differences in eigenvector centrality using permutation testing. Fibromyalgia patients, compared with HC, had higher eigenvector centrality in the bilateral anterior insulae, bilateral STG, 2 nodes in the left M1/operculum, 2 nodes in the right mid insula, right inferior parietal lobule, left precuneus, and the right posterior cingulate. Healthy controls had greater eigenvector centrality in the bilateral mid temporal gyrus, and medial prefrontal cortex. All significant differences are listed in Table 2. Hub regions and differences between groups were similar when measured with degree (Supplementary Fig. 2 and Supplementary Tables 7 and 8, available at http://links.lww.com/PAIN/A729). Analyses were repeated in a subset of FM patients (n = 28) and HCs (n = 32) who also had depression and anxiety measures. Most of the differences between groups (including M1, mid insula, and the bilateral anterior insulae) remain significant after controlling for depression and anxiety (all P < 0.05; data not shown).
3.4. Altered rich-club membership in fibromyalgia
Next, we examined higher-order rich-club classification in patients and controls. The rich-club curves for the group averaged networks at 5% density are displayed in Figure 3A. Both FM and HC brain networks had a rich-club organization, which differed from random networks, but the rich-club coefficients were not significantly different between groups. In FM, the rich-club regime (defined as ϕ (k) significantly different from ϕrandom (k), and ϕnorm (k) > 1) was between k = 4 to k = 6 and k = 8 to k = 22. In HC, the rich-club regime was over a range of k = 2 to k = 25. The illustrated rich-club regions for FM and HC (k = 22 for FM and k = 25 for HC; the highest level of k that was significantly different from random networks for each group) are depicted in Figures 3B and C. Although there was no significant difference in the level of rich-club organization between FM and HCs, the hubs that constituted the rich club were different between groups. Most importantly, the bilateral anterior insulae, anterior and mid cingulate cortices were members of the rich club in FM patients but not controls. Results were similar for most other network densities tested and are illustrated in Supplementary Fig. 3 (available at http://links.lww.com/PAIN/A729).
3.5. Rich-club membership is associated with clinical pain
To determine if rich-club membership was related to clinical pain, we divided the FM group into age-matched tertiles based on the level of clinical pain on the day of the scan. As expected, the tertiles had significantly different pain ratings (VAS mean ± SD, high: 7.0 ± 1.2, medium: 5.1 ± 1.5, low: 2.5 ± 1.6, F(2,37) = 29.02, P < 10−8). There were no significant differences between the tertiles in age (F(2,37) = 0.55, P = 0.58). The rich-club coefficient, ϕ (k), was not significantly different between the tertiles for any level of k; however, the hubs comprising the rich club in each tertile were markedly different (Fig. 4). The rich-club nodes in the high pain group were predominantly in S1/M1, STG, and the anterior and posterior insula. The medium pain group rich club included fewer S1/M1 regions and the anterior insula, with a general shift toward more posterior and default mode network (DMN) regions. Finally, the rich club in the low pain group primarily contained nodes in the frontoparietal, DMN, and visual networks. Results were similar across all other thresholds tested and are illustrated in Supplementary Fig. 4 (available at http://links.lww.com/PAIN/A729).
3.6. The relationship between eigenvector centrality and clinical pain
In an exploratory analysis, we correlated the eigenvector centrality of all nodes in the brain network with the level of clinical pain in FM patients. The eigenvector centrality of 14 nodes in bilateral S1/M1, 3 nodes in the bilateral STG, 2 nodes in the right posterior insula, and 1 node each in the right SMA and right supramarginal gyrus positively correlated with clinical pain (all r > 0.3, P < 0.05, Fig. 5 and Supplementary Table 9; available at http://links.lww.com/PAIN/A729). Five of these nodes (3 in bilateral M1, 1 in the supramarginal gyrus, and 1 in the right posterior insula) belonged to the rich club only in the high pain FM tertile. Again, most of these correlations remained significant after controlling for depression and anxiety in the subset of FM patients with these measures (all P < 0.05; data not shown).
3.7. Graph theory metrics mediate the relationship between glutamate within the posterior insula and clinical pain intensity
As the right posterior insula was a member of the rich club in FM patients with the highest levels of clinical pain, we sought to examine the relationship between hub strength, clinical pain, and Glx. Right posterior insula Glx (CSF-corrected values) was significantly different between the FM pain tertiles and HC (F(3,62) = 4.37, P = 0.007). Specifically, the high pain group had significantly higher Glx in the right posterior insula compared with the low pain group (mean ± SD high: 12.76 ± 1.13, low: 11.20 ± 1.42, P = 0.013) and compared with HC (mean ± SD HC: 10.97, P = 0.001).
Two nodes overlapped with the 1H-MRS voxel in the right posterior insula (Fig. 6A). The standardized and averaged eigenvector centrality of these 2 posterior insula nodes positively correlated with clinical pain (r = 0.419, P = 0.008, Fig. 6B) and Glx (r = 0.320, P = 0.05; Fig. 6B). In mediation analyses, higher posterior insula Glx was associated with greater clinical pain indirectly through increased averaged posterior insula eigenvector centrality (B = 0.132; 95% CI = 0.016-0.399). The direct effect of insula Glx on clinical pain was also significant (B = 0.387, P = 0.026). The R2 value for clinical pain in this model was 0.422. See Figure 6C for standardized estimates and CIs of each path in the model.
3.8. Validation of results in a novel fibromyalgia cohort
We repeated a majority of the analyses in a small cohort of 11 FM patients. The discovery and validation of FM patients did not differ in respect to age (mean ± SD: 39.6 ± 12.1, P = 0.87) or clinical pain intensity (mean ± SD: 5.1 ± 2.1, P = 0.76). See supplementary Table 10 for a summary of the validation results (available at http://links.lww.com/PAIN/A729). The anterior insula, M1, and S1 were not hub brain regions in this sample, and group differences in eigenvector centrality did not replicate (Supplementary Fig. 5 and Supplementary Tables 11 and 12, available at http://links.lww.com/PAIN/A729). Similar to the FM discovery sample and HCs, the validation FM brain networks possessed a rich-club organization that was significantly different from random networks. Similar to the FM discovery cohort, the FM validation rich club included the mid and dorsal anterior cingulate cortex; however, the anterior insula was not a member of the rich club (Supplementary Fig. 6, available at http://links.lww.com/PAIN/A729). The low and high pain validation groups had similar mean clinical pain intensity ratings compared with the low and high discovery tertiles (low FM discovery VAS mean ± SD: 2.5 ± 1.6; low FM validation VAS mean ± SD: 3.18 ± 1.28; high FM discovery VAS mean ± SD: 7.0 ± 1.2; and high FM validation VAS mean ± SD: 6.72 ± 0.75). Similar to the FM discovery cohort, the bilateral mid insula were members of the rich club in the high pain validation group, but not the low pain validation group (Supplementary Fig. 7A, available at http://links.lww.com/PAIN/A729). At higher thresholds, the anterior insula, mid insula, M1, and S1 were hubs only in the high pain validation group (Supplementary Figs. 7B–H, available at http://links.lww.com/PAIN/A729). Finally, we sought to replicate the correlations between eigenvector centrality and clinical pain. Eigenvector centrality in the right anterior insula (r = 0.601, P = 0.05), left posterior insula (r = 0.659, P = 0.03), and left mid insula (r = 0.644, P = 0.03) positively correlated with clinical pain intensity (Supplementary Fig. 8, available at http://links.lww.com/PAIN/A729).
In the discovery cohort, we have shown that hubs, the brain regions critical to the efficient transfer and integration of information, are altered in FM patients. Several regions linked to pain and sensory processing such as the anterior insula, STG, and M1 were hubs in FM but not in HCs. The bilateral anterior insulae, STG, right mid insula, and left M1 had more connections to other nodes with high functional importance (ie, higher eigenvector centrality) in FM compared with HCs. Furthermore, in FM, the anterior insulae, mid, and anterior cingulate were members of the higher-order hub structure (the rich club), and rich-club membership varied with the level of clinical pain, such that the somatosensory, motor, and insular cortices belonged to the rich club in patients experiencing high levels of clinical pain. Finally, this study demonstrates, for the first time to our knowledge, a neurochemical correlate of altered hub topology and its relationship to the perception of pain. Taken together, these data suggest that sensory and pain-related brain regions play an outsized role in information flow and integration in the brains of FM patients.
We attempted to replicate these findings in a small FM validation cohort. The anterior insula was not a hub or a member of the rich club in the validation group, nor was the eigenvector centrality of this region significantly higher than controls. However, the mid and anterior cingulate were members of the rich club, and rich-club membership varied with clinical pain such that the mid insulae were also rich-club hubs in the high pain validation group. Furthermore, the eigenvector centrality of the anterior, mid, and posterior insula positively correlated with clinical pain in the validation group.
Previous studies have examined network organization in other chronic pain conditions. Chronic back pain patients had increased functional connectivity between hubs in the DMN and the insula relative to HCs.57 Similarly, increased cross-network communication between the DMN and the salience network was found in patients with ankylosing spondylitis, which also correlated with the level of clinical pain.26 The most comprehensive analysis to date found altered whole-brain degree rank order in 3 chronic pain populations relative to HCs.36 However, this study did not examine hubs specifically. Instead, it assessed the variability in community membership. The insula was a member of the sensorimotor network in most HCs; however, allegiance in chronic pain conditions was split between the sensorimotor, default mode, and subcortical networks. Our finding of increased hub strength in the bilateral anterior insulae in FM is consistent with these studies.
A key attribute of brain networks is the rich-club organization, in which high-degree hubs form a densely interconnected core. The rich club may serve to integrate information because rich-club hubs are often distributed in many different resting state networks.59 In the current study, we found both FM and HC networks possessed, by global assessment, a similar level of rich-club organization. However, the specific membership of the rich club varied between groups. This is a critical finding because it suggests that information flow in FM brain networks is qualitatively, rather than quantitatively, different. The salience network, while completely absent in the rich club of HCs, was represented in the FM rich club and included the bilateral anterior insulae, mid, and anterior cingulate cortices. The rich-club visual nodes in HCs formed a stand-alone network, consistent with previous studies showing a low level of integration in the functional rich club of healthy adults.21 Interestingly, in FM, the visual nodes were integrated into the rich club. Balenzuela and colleagues reported a disruption in the community membership in chronic back pain patients such that the insula was abnormally integrated into the auditory system and dorsal visual stream.2 Given our previous findings of visual hypersensitivity and increased insular activation during an unpleasant visual stimulus in FM,25 it is possible that altered rich-club membership and integration contributes to the phenomenon of multimodal sensory sensitivity observed in many chronic pain patients. This hypothesis would also suggest that targeting FM-specific rich-club nodes could improve both pain and aversive nonpainful experiences such as sensitivity to light and sound.
Within our FM patients, rich-club membership varied as a function of clinical pain. In general, the distribution of rich-club nodes shifted from the default mode, frontoparietal, and visual networks in the low pain FM group, to the sensorimotor, cinguloopercular, and salience networks in the high pain FM group. Patients with the highest levels of clinical pain had rich-club nodes mainly in S1, M1, SMA, STG, operculum, anterior, and posterior insula. The high pain rich club had just 1 DMN node and 7 nodes in the visual network, whereas the low pain group had 13 DMN and 20 visual network regions. These results suggest that the intensity of chronic pain is related to rich-club membership.
We also found that the eigenvector centrality of the posterior insula was positively correlated with the level of clinical pain in all FM patients, and that posterior insula Glx was significantly higher in the high pain tertile compared with the low pain FM group. Importantly, we demonstrated that eigenvector centrality of nodes within the posterior insula region showing elevated Glx mediates the relationship between this neurotransmitter and clinical pain. Thus, altered neurochemistry in an individual hub could alter functional network properties that result in the experience of pain. These data are consistent with a recent preclinical study wherein either increasing glutamate or decreasing GABA in the insula led to increased hyperalgesia and allodynia in naive rats.63
The insula's role as a sensory integration region that is critical for pain perception is also supported by studies in which direct electrical stimulation of the posterior insula elicited painful sensations,38,45 and insular lesions alter pain perception.16 We previously showed that FM patients have decreased γ-aminobutyric acid (GABA), the brain's major inhibitory neurotransmitter, in the anterior insula and increased levels of posterior insula Glx relative to HCs. These alterations in GABA and Glx were associated with increased experimental pain sensitivity.13,23 Furthermore, posterior insula Glx decreased after treatments that reduced pain in FM.22,24 Together with the preclinical data,63 these findings indicate that the excitatory-inhibitory balance in the insula exerts a causal influence on pain perception.
Why might hubs be preferentially affected in chronic pain? Stam 56 hypothesized that since hubs are damaged or reorganized in many different neurological disorders, there must be a general mechanism that makes hub brain regions more vulnerable. He suggests that, in the acute phase, traffic to a failing node is rerouted to existing hubs, which become overloaded. This would be observed as an increase in the degree or centrality of hubs relative to healthy networks. In the chronic phase, new hubs seem to avoid or compensate for hub overload. This would be reflected as a decrease in the degree or centrality of “healthy” hubs and the appearance of new hubs.
In the chronic pain patients studied here, we demonstrated the presence of sensory/pain hubs (anterior/mid insula, M1, and S1) in FM and the absence of certain hubs seen in HCs (visual cortex and SFG). The lack of convincing evidence of peripheral pathology in FM4 suggests that this reorganization is not due to a barrage of nociceptive input, but rather is more likely due to alteration in excitatory neurotransmission in pronociceptive brain regions.63 The mediation analysis described here is consistent with a model in which higher Glx in the posterior insula increases the eigenvector centrality of posterior insular nodes and, consequently, clinical pain. Longitudinal analyses will be required to determine whether these effects are causal in chronic pain patients.
A methodological limitation of graph theoretical analyses is the arbitrary thresholding and binarization process, which leads to the loss of information, specifically, anticorrelations. However, binary networks retain the most significant correlations and may be less sensitive to noise and easier to interpret.50 We examined a relatively small number of female patients, and many of the patients studied here were taking medications to manage their symptoms. Finally, we have published reports of increased functional connectivity between the insula and the DMN in some of the same FM patients studied here.28,43 However, the previous studies did not use a whole-brain network approach nor did they examine the functional organization of these networks.
In conclusion, we demonstrated altered hub topology in FM and showed that disruptions in the excitatory tone within the insula altered the strength of this region as a hub and was associated with clinical pain intensity. Although cross-sectional in nature, these data suggest the possibility that disruption and reallocation of rich-club hub membership plays a causal role in the symptomatology of FM. Noninvasive brain stimulation may alter network dynamics and the neurochemistry of regions underneath and distal to the stimulating electrode.12,41,55 Future studies should explore targeting hub brain regions to personalize treatment in chronic pain patients.
Conflict of interest statement
D.J. Clauw has received research funding from Aptinyx, Cerephex and Pfizer Inc, and consulted for Abbott Pharmaceuticals, Aptinyx, Astellas Pharmaceuticals, Cerephex, Daiichi Sankyo, Pfizer Inc, Pierre Fabre, Samumed, Theravance, Tonix, and Zynerba. R.E. Harris has received research funding and consulted for Pfizer Inc. S.E. Harte has received research funding from Aptinyx, Cerephex, Forest Laboratories, Eli Lily and Merck, and served as a consultant for Pfizer, Analgesic Solutions, Longitude Capital Management, Aptinyx, and deCode Genetics. He is a member of Arbor Medical Innovations, Ann Arbor, MI. The remaining authors have no conflicts of interest to declare.
The authors thank Drs. UnCheol Lee and Hyoungkyu Kim for insightful discussions.
This work was supported in part by NIH grant R01 AT007550 to R.E. Harris; a grant from Cerephex Corporation (Los Altos, California) to D.J. Clauw, G.M. Mashour, and R.E. Harris; a grant from Pfizer Inc (Groton, Connecticut) to R.E. Harris; a pilot grant from the Functional MRI Laboratory of the University of Michigan (PF U034199); and a grant from Forest Laboratories (MD-SAV-09) to D.J. Clauw.
Appendix A. Supplemental digital content
Supplemental digital content associated with this article can be found online at http://links.lww.com/PAIN/A729.
. Achard S, Delon-Martin C, Vértes PE, Renard F, Schenck M, Schneider F, Heinrich C, Kremer S, Bullmore ET. Hubs of brain functional networks are radically reorganized in comatose patients. Proc Natl Acad Sci U S A 2012;109:20608–13.
. Balenzuela P, Chernomoretz A, Fraiman D, Cifre I, Sitges C, Montoya P, Chialvo DR. Modular organization of brain resting state networks in chronic back pain patients. Front Neuroinform 2010;4:116.
. Behzadi Y, Restom K, Liau J, Liu TT. A component based noise correction method (CompCor) for BOLD and perfusion based fMRI. Neuroimage 2007;37:90–101.
. Clauw DJ. Fibromyalgia: a clinical review. JAMA 2014;311:1547–55.
. Cohen AL, Fair DA, Dosenbach NUF, Miezin FM, Dierker D, Van Essen DC, Schlaggar BL, Petersen SE. Defining functional areas in individual human brains using resting functional connectivity MRI. Neuroimage 2008;41:45–57.
. Cole MW, Reynolds JR, Power JD, Repovs G, Anticevic A, Braver TS. Multi-task connectivity reveals flexible hubs for adaptive task control. Nat Neurosci 2013;16:1348–55.
. Colizza V, Flammini A, Serrano MA, Vespignani A. Detecting rich-club ordering in complex networks. Nat Phys 2006;2:110–15.
. Crossley NA, Mechelli A, Scott J, Carletti F, Fox PT, McGuire P, Bullmore ET. The hubs of the human connectome are generally implicated in the anatomy of brain disorders. Brain 2014;137(pt 8):2382–95.
. Dickson J, Drury H, Van Essen DC. “The surface management system” (SuMS) database: a surface-based database to aid cortical surface reconstruction, visualization and analysis. Philos Trans R Soc Lond B Biol Sci 2001;356:1277–92.
. Dosenbach NUF, Fair DA, Miezin FM, Cohen AL, Wenger KK, Dosenbach RAT, Fox MD, Snyder AZ, Vincent JL, Raichle ME, Schlaggar BL, Petersen SE. Distinct brain networks for adaptive and stable task control in humans. Proc Natl Acad Sci U S A 2007;104:11073–8.
. Farmer MA, Baliki MN, Apkarian AV. A dynamic network perspective of chronic pain. Neurosci Lett 2012;520:197–203.
. Foerster BR, Nascimento TD, DeBoer M, Bender MA, Rice IC, Truong DQ, Bikson M, Clauw DJ, Zubieta JK, Harris RE, DaSilva AF. Excitatory and inhibitory brain metabolites as targets and predictors of effective motor cortex tDCS therapy in fibromyalgia. Arthritis Rheumatol 2015;67:576–81.
. Foerster BR, Petrou M, Edden RAE, Sundgren PC, Schmidt-Wilcke T, Lowe SE, Harte SE, Clauw DJ, Harris RE. Reduced insular γ-aminobutyric acid in fibromyalgia. Arthritis Rheumatol 2012;64:579–83.
. Fornito A, Zalesky A, Bullmore E. Fundamentals of Brain Network Analysis. Amsterdam, the Netherlands; Boston, MA: Elsevier/Academic Press, 2016.
. Fox MD, Snyder AZ, Vincent JL, Corbetta M, Van Essen DC, Raichle ME. The human brain is intrinsically organized into dynamic, anticorrelated functional networks. Proc Natl Acad Sci U S A 2005;102:9673–8.
. Garcia-Larrea L, Peyron R. Pain matrices and neuropathic pain matrices: a review. PAIN 2013;154(suppl 1):S29–43.
. Geisser ME, Glass JM, Rajcevska LD, Clauw DJ, Williams DA, Kileny PR, Gracely RH. A psychophysical study of auditory and pressure sensitivity in patients with fibromyalgia and healthy controls. J Pain 2008;9:417–22.
. Glover GH, Li TQ, Ress D. Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR. Magn Reson Med 2000;44:162–7.
. Zamora-López G, Zhou C, Kurths J. Cortical hubs form a module for multisensory integration on top of the hierarchy of cortical networks. Front Neuroinform 2010;4:1.
. Gracely RH, Petzke F, Wolf JM, Clauw DJ. Functional magnetic resonance imaging evidence of augmented pain processing in fibromyalgia. Arthritis Rheum 2002;46:1333–43.
. Grayson DS, Ray S, Carpenter S, Iyer S, Dias TG, Stevens C, Nigg DA, Fair DA. Structural and functional rich club
organization of the brain in children and adults. PLoS One 2014;9:e88297.
. Harris RE, Napadow V, Huggins JP, Pauer L, Kim J, Hampson J, Sundgren PC, Foerster B, Petrou M, Schmidt-Wilcke T, Clauw DJ. Pregabalin rectifies aberrant brain chemistry, connectivity, and functional response in chronic pain patients. Anesthesiology 2013;119:1453–64.
. Harris RE, Sundgren PC, Craig AD, Kirshenbaum E, Sen A, Napadow V, Clauw DJ. Elevated insular glutamate in fibromyalgia is associated with experimental pain. Arthritis Rheum 2009;60:3146–52.
. Harris RE, Sundgren PC, Pang Y, Hsu M, Petrou M, Kim SH, McLean SA, Gracely RH, Clauw DJ. Dynamic levels of glutamate within the insula are associated with improvements in multiple pain domains in fibromyalgia. Arthritis Rheum 2008;58:903–7.
. Harte SE, Ichesco E, Hampson JP, Peltier SJ, Schmidt-Wilcke T, Clauw DJ, Harris RE. Pharmacologic attenuation of cross-modal sensory augmentation within the chronic pain insula. PAIN 2016;157:1933–45.
. Hemington KS, Wu Q, Kucyi A, Inman RD, Davis KD. Abnormal cross-network functional connectivity in chronic pain and its association with clinical symptoms. Brain Struct Funct 2015;221:4203–19.
. Ichesco E, Puiu T, Hampson JP, Kairys AE, Clauw DJ, Harte SE, Peltier SJ, Harris RE, Schmidt-Wilcke T. Altered fMRI resting-state connectivity in individuals with fibromyalgia on acute pain stimulation. Eur J Pain 2016;20:1079–89.
. Ichesco E, Schmidt-Wilcke T, Bhavsar R, Clauw DJ, Peltier SJ, Kim J, Napadow V, Hampson JP, Kairys AE, Williams DA, Harris RE. Altered resting state connectivity of the insular cortex in individuals with fibromyalgia. J Pain 2014;15:815–26.e811.
. Jensen KB, Kosek E, Petzke F, Carville S, Fransson P, Marcus H, Williams SC, Choy E, Giesecke T, Mainguy Y, Gracely R, Ingvar M. Evidence of dysfunctional pain inhibition in fibromyalgia reflected in rACC during provoked pain. PAIN 2009;144:95–100.
. Power JD, Cohen AL, Nelson SM, Wig GS, Barnes KA, Church JA, Vogel AC, Laumann TO, Miezin FM, Schlaggar BL, Petersen SE. Functional network organization of the human brain. Neuron 2011;72:665.
. Kucyi A, Davis KD. The dynamic pain connectome. Trends Neurosci 2015;38:86–95.
. Liu J, Zhao L, Li G, Xiong S, Nan J, Li J, Yuan K, von Deneen KM, Liang F, Qin W, Tian J. Hierarchical alteration of brain structural and functional networks in female migraine sufferers. PLoS One 2012;7:e51250.
. López-Solà M, Pujol J, Wager TD, Garcia-Fontanals A, Blanco-Hinojo L, Garcia-Blanco S, Poca-Dias V, Harrison BJ, Contreras-Rodríguez O, Monfort J, Garcia-Fructuoso F, Deus J. Altered functional magnetic resonance imaging responses to nonpainful sensory stimulation in fibromyalgia patients. Arthritis Rheumatol 2014;66:3200–9.
. López-Solà M, Woo CW, Pujol J, Deus J, Harrison BJ, Monfort J, Wager TD. Towards a neurophysiological signature for fibromyalgia. PAIN 2017;158:34–47.
. Lynall ME, Bassett DS, Kerwin R, McKenna PJ, Kitzbichler M, Muller U, Bullmore E. Functional connectivity and brain networks in schizophrenia. J Neurosci 2010;30:9477–87.
. Mansour A, Baria AT, Tetreault P, Vachon-Presseau E, Chang PC, Huang L, Apkarian AV, Baliki MN. Global disruption of degree rank order: a hallmark of chronic pain. Sci Rep 2016;6:34853.
. Martenson ME, Halawa OI, Tonsfeldt KJ, Maxwell CA, Hammack N, Mist SD, Pennesi ME, Bennett RM, Mauer KM, Jones KD, Heinricher MM. A possible neural mechanism for photosensitivity in chronic pain. PAIN 2016;157:868–78.
. Mazzola L, Isnard J, Peyron R, Guénot M, Mauguière F. Somatotopic organization of pain responses to direct electrical stimulation of the human insular cortex. PAIN 2009;146:99–104.
. Moon JY, Kim J, Ko TW, Kim M, Iturria-Medina Y, Choi JH, Lee J, Mashour GA, Lee U. Structure shapes dynamics and directionality in diverse brain networks: mathematical principles and empirical confirmation in three species. Sci Rep 2017;7:46606.
. Moon JY, Lee U, Blain-Moraes S, Mashour GA. General relationship of global topology, local dynamics, and directionality in large-scale brain networks. PLoS Comput Biol 2015;11:e1004225.
. Muldoon SF, Pasqualetti F, Gu S, Cieslak M, Grafton ST, Vettel JM, Bassett DS. Stimulation-based control of dynamic brain networks. PLoS Comput Biol 2016;12:e1005076.
. Napadow V, Kim J, Clauw DJ, Harris RE. Brief Report: decreased intrinsic brain connectivity is associated with reduced clinical pain in fibromyalgia. Arthritis Rheum 2012;64:2398–403.
. Napadow V, LaCount L, Park K, As-Sanie S, Clauw DJ, Harris RE. Intrinsic brain connectivity in fibromyalgia is associated with chronic pain intensity. Arthritis Rheum 2010;62:2545–55.
. Nichols TE, Holmes AP. Nonparametric permutation tests for functional neuroimaging: a primer with examples. Hum Brain Mapp 2002;15:1–25.
. Ostrowsky K, Magnin M, Ryvlin P, Isnard J, Guénot M, Mauguière F. Representation of pain and somatic sensation in the human insula: a study of responses to direct electrical cortical stimulation. Cereb Cortex 2002;12:376–85.
. Petzke F, Clauw DJ, Ambrose K, Khine A, Gracely RH. Increased pain sensitivity in fibromyalgia: effects of stimulus type and mode of presentation. PAIN 2003;105:403–13.
. Power JD, Barnes KA, Snyder AZ, Schlaggar BL, Petersen SE. Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion. Neuroimage 2012;59:2142–54.
. Provencher SW. Estimation of metabolite concentrations from localized in vivo proton NMR spectra. Magn Reson Med 1993;30:672–9.
. Puiu T, Kairys AE, Pauer L, Schmidt-Wilcke T, Ichesco E, Hampson JP, Napadow V, Clauw DJ, Harris RE. Association of alterations in gray matter volume with reduced evoked-pain connectivity following short-term administration of pregabalin in patients with fibromyalgia. Arthritis Rheumatol 2016;68:1511–21.
. Rubinov M, Sporns O. Complex network measures of brain connectivity: uses and interpretations. Neuroimage 2010;52:1059–69.
. Schmidt-Wilcke T, Ichesco E, Hampson JP, Kairys A, Peltier S, Harte S, Clauw DJ, Harris RE. Resting state connectivity correlates with drug and placebo response in fibromyalgia patients. Neuroimage Clin 2014;6:252–61.
. Sluka KA, Clauw DJ. Neurobiology of fibromyalgia and chronic widespread pain. Neuroscience 2016;338:114–29.
. Sporns O, Honey CJ, Kötter R. Identification and classification of hubs in brain networks. PLoS One 2007;2:e1049.
. Spreng RN, Sepulcre J, Turner GR, Stevens WD, Schacter DL. Intrinsic architecture underlying the relations among the default, dorsal attention, and frontoparietal control networks of the human brain. J Cogn Neurosci 2013;25:74–86.
. Stagg CJ, Best JG, Stephenson MC, O'Shea J, Wylezinska M, Kincses ZT, Morris PG, Matthews PM, Johansen-Berg H. Polarity-sensitive modulation of cortical neurotransmitters by transcranial stimulation. J Neurosci 2009;29:5202–6.
. Stam CJ. Modern network science of neurological disorders. Nat Rev Neurosci 2014;15:683–95.
. Tagliazucchi E, Balenzuela P, Fraiman D, Chialvo DR. Brain resting state is disrupted in chronic back pain patients. Neurosci Lett 2010;485:26–31.
. van den Heuvel MP, Kahn RS, Goñi J, Sporns O. High-cost, high-capacity backbone for global brain communication. Proc Natl Acad Sci U S A 2012;109:11372–7.
. van den Heuvel MP, Sporns O. Rich-club organization of the human connectome. J Neurosci 2011;31:15775–86.
. van den Heuvel MP, Sporns O. Network hubs in the human brain. Trends Cogn Sci 2013;17:683–96.
. Van Essen DC, Dickson J, Harwell J, Hanlon D, Anderson CH, Drury HA. An Integrated Software System for Surface-based Analyses of Cerebral Cortex. Journal of American Medical Informatics Association 2001;8:443–59.
. Vatansever D, Menon DK, Manktelow AE, Sahakian BJ, Stamatakis EA. Default mode dynamics for global functional integration. J Neurosci 2015;35:15254–62.
. Watson CJ. Insular balance of glutamatergic and GABAergic signaling modulates pain processing. PAIN 2016;157:2194–207.
. Whitfield-Gabrieli S, Nieto-Castanon A. Conn: a functional connectivity toolbox for correlated and anticorrelated brain networks. Brain Connect 2012;2:125–41.
. Wilbarger JL, Cook DB. Multisensory hypersensitivity in women with fibromyalgia: implications for well being and intervention. Arch Phys Med Rehabil 2011;92:653–6.
. Xia M, Wang J, He Y. BrainNet Viewer: a network visualization tool for human brain connectomics. PLoS One 2013;8:e68910.
. Zigmond AS, Snaith RP. The hospital anxiety and depression scale. Acta Psychiatr Scand 1983;67:361–70.