Dynamic Pain Connectome Functional Connectivity and Oscillations Reflect Multiple Sclerosis Pain

Rachael L. Bosma; Junseok A. Kim; Joshua C. Cheng; Anton Rogachov; Kasey S. Hemington; Natalie R. Osborne; Jiwon Oh; Karen D. Davis


Pain. 2018;159(11):2267-2276. 

In This Article

Materials and Methods


Participants included 31 patients with MS (20 women; mean ± SD age = 37 ± 10 years, 11 men; mean age = 42 ± 8 years), and 31 age- and sex-matched healthy controls (HCs). Patients with MS were recruited from the St. Michael's Hospital Multiple Sclerosis clinic and HCs were recruited from the community. All participants provided informed written consent to procedures approved by the University Health Network and St. Michael's Hospital Research Ethics Boards in accordance with the Declaration of Helsinki. All patients had a diagnosis of MS and were excluded from the study if they had a painful disease unrelated/in addition to MS-related pain (eg, a diagnosis of arthritis), if they required a mobility aid, or if they had contraindications for MRI. The inclusion criteria for HCs were: (1) no previous history of chronic pain or current experience of pain on a regular basis, (2) no previous diagnosis of a psychiatric disorder, metabolic, or neurological condition, (3) no major surgery in the past 2 years, and (4) no standard contraindications for MRI.

Clinical Measures: Brief Pain Inventory and painDETECT

All patients were assessed at St. Michael's Hospital in Toronto, ON, Canada, and clinical information including scores on the Expanded Disability Status Scale, a highly valid and reliable method of assessing disability, which ranges from 0 to 10 and is scored by a neurologist in 0.5 unit increments in which higher scores indicate higher levels of disability,[39,45] disease symptom onset, date of disease diagnosis, and MS subtype were obtained. In one study session, patients underwent the neuroimaging protocol and completed a number of self-administered questionnaires including the Hospital Anxiety and Depression Scale.[58] The Hospital Anxiety and Depression Scale assesses the nonphysical symptoms of anxiety and depression, and scores above 8 are considered clinically relevant.[5] General pain and pain interference was assessed using the Brief Pain Inventory,[48] and the painDETECT[20] was used to assess NP. The Brief Pain Inventory measures 2 domains of pain: pain severity (worst pain, least pain, average pain, and pain now) and pain interference (how pain interferes with: general activity, walking, work, mood, enjoyment of life, relations with others, and sleep).[16,55] Participants rate each item on a scale from 0 (no pain or does not interfere) to 10 (pain as bad as you can imagine, completely interferes). For each patient, a composite (mean) score was calculated for pain interference measures to reflect quality of life measures. The painDETECT scale has been shown to be a valid and reliable assessment for the presence of NP features.[20,54] Scores on this scale range from −1 to 38. Scores below 13 indicate that a neuropathic component is unlikely, scores of 13 to 18 indicate that a neuropathic component may be present (and thus, these patients may have mixed neuropathic and nociceptive pain), and scores above 18 indicate that a neuropathic component is highly likely.[20] Consistent with previous studies,[53] we categorized patients with MS according to their painDETECT scores into mixed NP (scores above 12) and non-NP groups. In addition, patients with scores between 13 and 18 (mixed) also had sensory loss as determined by their mechanical detection thresholds, which lends support to classifying them as having NP features (Supplementary Figure 1, available online as supplemental digital content at http://links.lww.com/PAIN/A617). Participants were also asked about the chronicity and location of their pain (Supplementary Table 2, available online as supplemental digital content at http://links.lww.com/PAIN/A617).

Neuroimaging Acquisition

Both HCs and patients with MS underwent an MRI (3T GE) brain imaging session. Data acquisition included a T1-weighted anatomical scan (1 × 1 × 1 mm3 voxels, matrix = 256 × 256, 180 axial slices, repetition time = 7.8 seconds, echo time = 3 ms, inversion time = 450 ms) and a 9-minute, 14-second T2*-weighted rsfMRI scan (3.125 × 3.125 × 4 mm3 voxels, matrix = 64 × 64, 36 axial slices, repetition time = 2 seconds, echo time = 30 ms, flip angle = 85°, 277 volumes). For the resting-state scan, participants were instructed to "close your eyes, do not try to think about anything in particular; do not fall asleep."

Preprocessing of Resting-state Functional Magnetic Resonance Imaging Data

The preprocessing of the rsfMRI data was conducted using the fMRI expert analysis tool (FEAT) toolbox in FMRIB Software Library (FSL).[30] The following standard procedures were applied: (1) the first 4 volumes of the rsfMRI scan were removed, (2) nonbrain tissues were removed using the Brain Extract Tool (BET) within FEAT, and (3) motion correction was performed using MCFLIRT. The rsfMRI data from each participant were then linearly registered to their high-resolution anatomical image, which was skull-stripped using the opti-BET tool,[42] followed by nonlinear registration to MNI152–2 mm space using FNIRT. Scanner-related noise and physiological noise were removed by means of applying aCompCor[4,10] as described previously.[38] The aCompCor method implements principal component analysis on the thresholded (to prevent the removal of gray matter) white matter and cerebrospinal fluid maps to generate components with the greatest variance. Five components containing the most variance were removed from the resting state data as well as the 6 motion parameters derived from MCFLIRT were regressed out from each participant's rsfMRI data. Finally, spatial smoothing was applied using a 4-mm full-width at half-maximum kernel, and temporal filtering was performed in FSL to retain the signal between 0.01 and 0.1 Hz.

Definition of Seeds

To examine functional reorganization in regions of the dynamic pain connectome, we selected key regions of the salience (rTPJ) and DMN (posterior cingulate cortex [PCC] and medial prefrontal cortex [mPFC]), as well as the ascending nociceptive (postcentral gyrus; S1) and descending modulation (periaqueductal gray region [PAG]) nociceptive pathways. The location of the rTPJ seed (MNI x = 62, y = −36, x = 28) was based on previous studies that identified this region of the anterior rTPJ to highly coactivate with[32] and be highly connected to other regions involved in attention and salience.[35] The seeds in the PCC (MNI x = −8, y = −50, z = 28) and mPFC (MNI x = −2, y = 56, z = 16) were based on the coordinates from our previous study of the DMN,[24] and the seed within the PAG region (MNI x = 0, y = −32, z = −10) was according to previous studies.[17] For the connectivity analyses, we created an 8-mm spherical seed at the peak coordinates of the rTPJ, PCC, mPFC, and a 3-mm radius seed in the PAG. The S1 mask was constructed from the Harvard Oxford Cortical Structural Atlas to include all of the gray matter in the postcentral gyrus. A nonlinear transformation of each seed region from standard space to each subject's MRI space was performed and the mean time series from each was extracted.

Brain Communication: Static Functional Connectivity Analysis

The sFC was measured between the seeds of salience and DMNs (average sFC between rTPJ-PCC and rTPJ-mPFC) and between the salience-ascending (rTPJ-SI) and salience-descending (rTPJ-PAG) by means of a Pearson correlation between the average time series of each seed pair. Correlation values were then Fisher r-to-z transformed. The sFC was compared between HCs and patients with MS, as well as between pain feature-based subgroups of patients with MS (neuropathic and non-neuropathic) using unpaired sample t tests. Statistical significance was determined at P < 0.05, with false discovery rate (FDR) controlling for multiple comparisons across all P values determined from group (MS vs HC) and subgroup (mixed NP vs Non-NP) comparisons across all pairwise sFC measures (SN-DMN, SN-Asc, and SN-Des).

Flexibility in Communication: Dynamic Functional Connectivity

Dynamic FC was also measured between the same seed pairs by calculating the dynamic conditional correlation method as described previously.[12,41] Briefly, to calculate the dynamic conditional correlation, (1) each time course was prewhitened using an autoregressive and moving average (ARMA) (1,1) model, (2) a generalized autoregressive conditional heteroscedastic (GARCH) model was applied to estimate the conditional SD over time, (3) the residuals of the time series were standardized by the conditional SD, and (4) the dynamic conditional correlation was then used to calculate the time-varying correlation between the time series, using an exponentially weighted moving average that is derived from the data using maximum likelihood.[41] The SD of each dynamic conditional correlation across the time series was computed and used as the summary metric of dFC.[13] High dFC indicates greater fluctuations in the strength of connectivity between 2 regions over time. The dFC was again compared between controls and patients, and between patient subgroups using independent-samples t tests. Multiple comparisons were corrected using FDR across all P values determined from group (MS vs HC) and subgroup (mixed-NP vs Non-NP) comparisons across all pairwise dFC measures (SN-DMN, SN-Asc, and SN-Des).

BOLD Signal Variability

The standardized BOLD signal variability was calculated on a voxel-wise basis within a mask containing the key regions described above (DMN, SN, ascending nociceptive regions, and descending modulation regions) as previously reported.[49] For each voxel within the mask, we calculated the SD of the BOLD signal, then subtracted the average SD from the gray matter across the whole brain, and then divided it by the SD of the BOLD signal variability across the gray matter.[43] Previous studies of low-frequency oscillations related to pain, such as those present in regional BOLD variability time courses, have found distinct patterns of response in different frequency bands.[1,25] Therefore, we also used AFNI's 3d Bandpass function to temporally filter the rsfMRI data such that we could examine BOLD variability in the different slow-wave frequency bands (slow-5: 0.01–0.027 Hz; slow-4: 0.027–0.073 Hz; and slow-3: 0.073–0.198 Hz).[8,14,59] To compare the BOLD signal variability in each band between controls and patients with MS, we used FSL's Randomise tool for nonparametric permutation inferencing, using positive and negative contrasts with 5000 permutations and a threshold-free cluster enhancement. The group-level results are thresholded at P < 0.05 (family-wise error-corrected for multiple comparisons).

Relationship Between Brain Functional Reorganization and Pain Interference

The association between brain measures (sFC, dFC, and BOLD variability) within the dynamic pain connectome and pain interference (from the Brief Pain Inventory) was assessed by means of a Spearman correlation. Correlation analyses for both sFC (pain interference vs sFC values) and dFC (pain interference vs dFC values) were corrected for multiple comparisons using FDR (P < 0.05). Regional BOLD variability was calculated for each voxel within the dynamic pain connectome mask; therefore, in addition to assessing the univariate relationship between BOLD variability and pain interference, we also examined the multivariate association between brain patterns of BOLD variability and pain interference scores using the Pattern Recognition for Neuroimaging Toolbox (PRoNTo) (Shrouff et al., 2013) and kernel ridge regression (k-fold cross-validation, k = 5 folds). We first used permutation testing to determine which voxels demonstrated significant group BOLD variability differences between patients and controls, as described above (P < 0.05, familywise error rate corrected). Next, the voxels that demonstrated group differences were entered into the multivariate regression to determine whether these group differences in BOLD variability could be used to make generalizable inferences about pain interference. For the multivariate regression analyses, permutation testing was also implemented (5000 iterations) to determine the statistical significance.