Sensitivity of functional connectivity to periaqueductal gray localization, with implications for identifying disease-related changes in chronic visceral pain: A MAPP Research Network neuroimaging study.

Previous studies examining the resting-state functional connectivity of the periaqueductal gray (PAG) in chronic visceral pain have localized PAG coordinates derived from BOLD responses to provoked acute pain. These coordinates appear to be several millimeters anterior of the anatomical location of the PAG. Therefore, we aimed to determine whether measures of PAG functional connectivity are sensitive to the localization technique, and if the localization approach has an impact on detecting disease-related differences in chronic visceral pain patients. We examined structural and resting-state functional MRI (rs-fMRI) images from 209 participants in the Multidisciplinary Approach to the Study of Chronic Pelvic Pain (MAPP) Research Network study. We applied three different localization techniques to define a region-of-interest (ROI) for the PAG: 1) a ROI previously-published as a Montreal Neurological Institute (MNI) coordinate surrounded by a 3 mm radius sphere (MNI-sphere), 2) a ROI that was hand-traced over the PAG in a MNI template brain (MNI-trace), and 3) a ROI that was hand-drawn over the PAG in structural images from 30 individual participants (participant-trace). We compared the correlation among the rs-fMRI signals from these PAG ROIs, as well as the functional connectivity of these ROIs with the whole brain. First, we found important non-uniformities in brainstem rs-fMRI signals, as rs-fMRI signals from the MNI-trace ROI were significantly more similar to the participant-trace ROI than to the MNI-sphere ROI. We then found that choice of ROI also impacts whole-brain functional connectivity, as measures of PAG functional connectivity throughout the brain were more similar between MNI-trace and participant-trace compared to MNI-sphere and participant-trace. Finally, we found that ROI choice impacts detection of disease-related differences, as functional connectivity differences between pelvic pain patients and healthy controls were much more apparent using the MNI-trace ROI compared to the MNI-sphere ROI. These results indicate that the ROI used to localize the PAG is critical, especially when examining brain functional connectivity changes in chronic visceral pain patients.

Introduction

The periaqueductal gray (PAG) is considered one of the most important centers of activity in the brain, integrating and controlling behavior crucial for survival (Koutsikou et al., 2017, Koutsikou et al., 2015, Linnman et al., 2012, An et al., 1998, Coulombe et al., 2016, Kong et al., 2010, Wei et al., 2016). Functions specific to the PAG in humans have been attributed to autonomic control (Linnman et al., 2012, An et al., 1998, Coulombe et al., 2016, Benarroch, 2012, Öngür and Price, 2000, Fowler et al., 2008), such as with micturition and cardiovascular regulation, but also relate to functions as diverse as fear and anxiety (Linnman et al., 2012, Mobbs et al., 2007, Silva and McNaughton, 2019), sexual behavior (Ferris et al., 2004), and antinociception (Bandler and Shipley, 1994). The PAG connects to the medial prefrontal cortical regions in Macaque monkeys (An et al., 1998, Öngür and Price, 2000) and governs descending pain modulation pathways in humans (Stamford, 1995, Tracey et al., 2002). These connections support an anatomical basis for top down control of spinal and sensorimotor circuits that could be particularly important in chronic pain disorders.

Defining brain regions and their relationships noninvasively through neuroimaging, such as with functional connectivity, has gained traction in recent literature (Coulombe et al., 2016, Tillman et al., 2018, Cohen et al., 2008, Eickhoff et al., 2018, Ezra et al., 2015). This approach is fundamental to understanding the brain function that may underlie human behavior (Eickhoff et al., 2018, Genon et al., 2018). Much work has been done to define the PAG in animals and humans invasively (Koutsikou et al., 2015, An et al., 1998, Öngür and Price, 2000, Bandler and Shipley, 1994, Bittar et al., 2005, Herrero et al., 1991, Richardson and Akil, 1977), but a benefit to neuroimaging is that it can measure human behavior in vivo, noninvasively (Linnman et al., 2012) (see Linnman et al. for a review on neuroimaging PAG). Anatomically, the PAG, a small, curved, hollow partial cylinder approximately 10–14 mm long and 4–5 mm in external diameter encircling the central aqueduct in the midbrain (n.b. the nuclei ventral to the aqueduct are distinct from the PAG), can be identified visually using high-resolution structural magnetic resonance sequences (Eippert et al., 2009 Aug, Tracey and Iannetti, 2006, Satpute et al., 2013). However, the cytoarchitectonic pattern and the behavioral role of the four longitudinal columns within the PAG, the dorsomedial, dorsolateral, lateral, and ventrolateral, well-established in invasive animal studies, remains unclear in humans (Linnman et al., 2012, Coulombe et al., 2016, Ezra et al., 2015, Tracey and Iannetti, 2006). Recent work has attempted to define these subregions of the PAG using data-driven rs-fMRI connectivity methods (Coulombe et al., 2016). Another approach activates the region through known functional stimuli, pain, fear, or autonomic regulation changes. The challenge with this type of measurement is in anatomically localizing the region of activation, optimally performed with additional landmarks and across a range of experimental conditions (Genon et al., 2018, Tracey and Iannetti, 2006, Poldrack et al., 2011, Costafreda, 2009). At this juncture, functional magnetic resonance imaging (fMRI) still provides marginal spatial resolution and significant scanner and physiological artifacts (Linnman et al., 2012, Eippert et al., 2009 Aug, Tracey and Iannetti, 2006, Satpute et al., 2013) when identifying brainstem regions.

Many functional neuroimaging studies have used PAG coordinates as a seed or region of interest (ROI) analysis based on previously established functional connectivity measures. However, a number of these studies (Kong et al., 2010, Wei et al., 2016, Wei et al., 2017, Harper et al., 2018, Chen et al., 2017, Yu et al., 2014, Li et al., 2016, Liu et al., 2017, Linnman et al., 2012) have modeled the PAG as a sphere, despite its curved cylindrical shape, around coordinates located several millimeters anterior to the anatomical location of the PAG. These coordinates were localized in a low and high heat pain fMRI contrast (Kong et al., 2010) and validated in a follow-up study analyzing the functional connectivity of the ventrolateral region of the PAG (Kong et al., 2010). Though the ventrolateral region is justified in this study, others have used these coordinates to represent the whole PAG. The discrepancy between the functional and the anatomical location of the PAG leads us to further establish benchmarks that define the true PAG and discern approaches that render functional neuroimaging more reproducible and robust.

We defined three different localization techniques for the PAG as a whole, including coordinates from the literature, and compared each ROI with various resting-state neuroimaging approaches. These three ROIs were labeled the MNI-sphere, coordinates from literature surrounded by a 3 mm radius sphere, the MNI-trace, a hand-traced ROI of the PAG from a MNI template, and the participant-trace, another hand-traced ROI in the structural space of individual participants (this trace was assumed to represent the optimal rendering of an individual’s true PAG region, or “gold-standard”). We measured differences in extracted time series signals of resting-state functional neuroimaging data from healthy controls, analyzing the signal directly from the PAG ROI defined in these three ways, and compared the functional connectivity throughout the brain across methods. For this study, we also demonstrated that the location of the ROI affects the differences seen in voxelwise whole-brain functional connectivity between urologic chronic pelvic pain syndrome (UCPPS) patients and healthy controls. We hypothesized that the differences that emerged between healthy control subjects and patients would differ based on the anatomic region identified as the ‘PAG’. UCPPS patients were used in previously published cross-sectional functional neuroimaging studies with involvement of the PAG (Kleinhans et al., 2016), but this is the first study to directly compare the whole-brain resting-state connectivity of the PAG in such a population.

Methods3

Participants

Our initial cohort consisted of 318 participants from the first phase of the multi-site Multidisciplinary Approach to the Study of Chronic Pelvic Pain (MAPP) Research Network study4, selected based on a clinical diagnosis of UCPPS or as a healthy control without a history of UCPPS. Participants were recruited from five collection sites: Northwestern University, University of California, Los Angeles, University of Michigan, University of Alabama at Birmingham, and Stanford University (Alger et al., 2016, Clemens et al., 2014, Bagarinao et al., 2014). All participants provided informed consent according to the Declaration of Helsinki and the institutional review board approved the collection at each site (Kutch et al., 2017). Selection criteria were based on whether neuroimaging data met specific quality standards as regards to motion and dataset balancing protocols.

The first step in our analysis examined PAG location in healthy controls. These participants were selected with the following goals in mind: 1) equalizing number of males and females, 2) spanning a wide age range, and 3) distributing across the five sites equally to account for any differences in acquisition technique. Due to the time-intensive work of hand tracing the PAG voxel by voxel for each participant, we selected 15 healthy controls (Table 1, n = 9 female, n = 6 male) after excluding three participants (n = 1 female, n = 2 males) due to image preprocessing issues.

Table 1

Demographic characteristics of participant datasets. Key: NW, Northwestern University, UCLA, University of California, Los Angeles, UAB, University of Alabama at Birmingham.

Demographic characteristics of participants
	age (y)	sex (male/female)	site
Healthy control datasets
Test (n = 15)	37.6 ± 19.1	6/9	NW, UCLA, Michigan, UAB, Stanford
Validation (n = 15)	37.1 ± 11.4	7/8	UCLA
Patients vs. healthy controls datasets
Patients (n = 100)	39.2 ± 13.3	34/66	NW, UCLA, Michigan, UAB, Stanford
Healthy controls (n = 109)	36.7 ± 12.2	34/75	NW, UCLA, Michigan, UAB, Stanford

An independent dataset (Table 1, n = 8 female, n = 7 male) with all participants from the University of California, Los Angeles, was matched to our first dataset for validation. For clarification, these two datasets are labeled as test, which is our first dataset, and validation. For our validation dataset, we followed the same procedure as our test dataset, but the selection was made from only one site to verify that the results from our test dataset were driven by the localization technique and not site differences. For example, if site differences drive our findings, the test and validation datasets would differ significantly since they originated from different site groupings. Participants from the validation dataset were selected only after removing participants included in the test dataset from the University of California, Los Angeles. The validation dataset was also matched in size and balanced by the number of males and females to the test dataset. Age is known to be related to volumetric variability of the brainstem (Lambert et al., 2013). Maximum deviation in age, an iterative selection of participants with the largest age difference, was calculated in the selection of participants in both test and validation datasets for males and females separately. The age distribution was not significantly different in our test and validation datasets (p = 0.92).

In our secondary analysis, our initial 318 MAPP cohort was reduced to 209 (Table 1, n = 100 patients, n = 109 healthy controls), based on the conditions described by Landis and colleagues (Landis et al., 2014) characterizing UCPPS patients. Furthermore, only natural history was recorded in MAPP (Clemens et al., 2014) and we did not control for treatment factors within patients. Inclusion criteria included 1) UCPPS or healthy control participants and 2) participants that passed image preprocessing quality control. This resulted in the exclusion of 68 positive control patients with non-urologic associated syndromes (Landis et al., 2014) and 41 participants based on image preprocessing failures. Focus was placed only on neuroimaging measures and no symptom data were utilized in study analysis.

Image collection

All 3D T1-weighted structural MRI and resting-state fMRI data were collected and quality controlled by MAPP on 3 T scanners (Alger et al., 2016). A high resolution structural T1-weighted magnetization-prepared rapid gradient-echo (MP-RAGE) pulse sequence was used to collect the data from each participant from all sites except Stanford University, which used an inversion-recovery fast spoiled gradient echo (IR-FSPGR) sequence. The University of California, Los Angeles collected slices with interleaved sequence, all other sites collected slices sequentially. Parameters for structural MRI collection included a repetition time (TR) = 2200 ms, echo time (TE) = 3.26 ms, slice thickness = 1 mm, 176 slices, 256x256 acquisition matrices, and voxel resolution = 1x1x1 mm (Alger et al., 2016, Kutch et al., 2017, Kutch et al., 2015). The participants closed their eyes for 10 min during the rs-fMRI, with 40-slice whole-brain volume, slice thickness = 4 mm, TR = 2000 ms, TE = 28 ms, flip angle = 77°, and FOV = 220. Further details are published in previous MAPP Research Network studies (Alger et al., 2016, Glover et al., 2012).

Image processing

Resting-state fMRI data were preprocessed with a software package developed by the Analysis of Functional NeuroImages (AFNI, , Version AFNI_18.1.18, (Cox, 1996) with parameters adjusted according to Plitt and colleagues and Drysdale and colleagues (Jo et al., 2013, Plitt et al., 2015, Drysdale et al., 2017 Jan). One of the benefits for this pipeline is its stringent correction for distance-dependent motion artifacts and use of local white matter regressor (AFNI’s ANATICOR). First, FreeSurfer () is used to segment the T1-weighted image and generate the white matter and ventricle masks. The rs-fMRI image was preprocessed by removing the first four volumes, despiked, slice-time corrected, co-registered to the T1-weighted image output from FreeSurfer, automatically resampled to 3 mm isotropic voxels, normalized to reflect the percent signal change, non-linearly registered to the MNI template scan, and spatially smoothed with a 6 mm full width half-maximum Gaussian kernel (Plitt et al., 2015, Drysdale et al., 2017 Jan). Nuisance variables included the average ventricular time series, local average white matter time series (AFNI’s ANATICOR), twelve parameter estimates for head motion, including their first derivative (Jo et al., 2013). The residual BOLD time series from the rs-fMRI images are simultaneously bandpass filtered at 0.01–0.1 Hz and nuisance regressed (AFNI’s 3dTproject). The time series is censored by removing volumes exceeding the Euclidean norm of the motion derivative at 0.3 mm or volumes with outliers, defined as time points a certain distance from the trend of the time series (see AFNI’s 3dToutcount), whose voxels exceed 10% of the masked voxels. The median number of volumes removed was 2 across participants, with a maximum number of 84 (Plitt et al., 2015, Drysdale et al., 2017 Jan). Alignment between the normalized rs-fMRI and normalized T1-weighted image was verified across the combined 30 participants from the test and validation datasets (see Supplemental Fig. 1).

PAG ROI definition

Three different localization techniques were used to define PAG ROIs. As seen in Fig. 1.A, these ROIs are defined as MNI-sphere, MNI-trace, and participant-trace. The MNI-sphere was a 3 mm radius sphere with a location previously published as MNI coordinates for the PAG, (MNI: 4 –26 −14). The MNI coordinates for the PAG can be found in Wei et al. (Wei et al., 2016). The MNI-trace ROI was created from a 1 mm MNI152 standard-space T1-weighted average structural template image and the participant-trace ROI was identified on each participant’s structural T1-weighted image in 1x1x1 mm grid-space using a combination of AFNI and FreeSurfer’s Freeview.

Fig. 1

Analysis of the BOLD rs-fMRI time series signal from the three different localization techniques of the PAG. A. Top left: An 3D illustration of the brain highlighting the brainstem. Top right: The average normalized rs-fMRI across participants in the healthy control dataset. Bottom: The three different regions of interest (ROIs) for the PAG: a previously published Montreal Neurological Institute (MNI) set of coordinates surrounded by a 3 mm radius sphere (MNI-sphere in green, shown top right in average image), a hand-traced ROI in a MNI template brain (MNI-trace in blue, shown top right in average image), and a hand-drawn ROI in each of the structural images from healthy control participants (participant-trace in red). B. Example rs-fMRI time series signals comparing the MNI-sphere ROI and the MNI-trace ROI to the participant-trace ROI. C. The differences in connectivity across all healthy controls (test dataset, n = 15, p < 0.001 and validation dataset, n = 15, p < 0.0001). The connectivity is given by the standardized correlation of the signal extracted from the MNI-sphere ROI to the participant-trace ROI and from the MNI-trace ROI to the participant-trace ROI. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The T1-weighted MPRAGE image allowed for optimal differentiation of grey and white matter to identify the PAG. Using AFNI and FreeSurfer’s Freeview to trace the PAG, the outline of the PAG ROI was identified using visual differentiation between grey and white matter. Since brainstem anatomy textbooks vary slightly in their demarcation of the caudal and rostral limits of the PAG, we consulted several individuals who work with the PAG on a regular basis and whose work has required its identification from different perspectives. These experts included brainstem anatomists, a neuroradiologist with expertise in the brainstem, and a physiologist who focuses on PAG function. A consensus was reached by all experts on the exact anatomic boundaries of the PAG on the anatomic images. Literature which informed our anatomic localization of the PAG included Linnman et al., and Satpute et al. (Linnman et al., 2012, Satpute et al., 2013), as well as standard anatomy textbooks. The caudal margin was considered as the location where the cerebral aqueduct no longer appeared as a circle, due its opening into the fourth ventricle. The rostral border of the PAG was defined as the point at which the cerebral aqueduct had definitively formed in the axial image from the third ventricle. This can be delineated using a combination of axial and sagittal images to determine the inferior border of the third ventricle. The general morphology of the PAG resembles a three-dimensional kite, which tapers as the cerebral aqueduct enters the fourth ventricle. We did not fix a set PAG length due to variations in brain size, which may affect the length of the cerebral aqueduct and the PAG, but rather allowed the border definitions to provide this information, intrinsic to the brain-shape of the individual participant. In instances of uncertainty about the definitive border of the PAG, a more conservative approach was taken: voxels that could not be identified as either PAG or surrounding white matter were not included in the PAG ROI mask. The cerebral aqueduct was then outlined and the area of the PAG was filled in. After the PAG ROI mask was created, the edges were further refined on a separate day. In total, each mask took approximately 1.5 h to complete. This process was completed for the PAG ROI of the MNI template, or the MNI-trace, as well as the 30 individual masks, or the participant-traces, from MAPP.

To assess intra-rater and inter-rater reliability, two separate manual tracings of the PAG were completed by a single rater for 10 participants and by two independent raters for 4 participants in a separate dataset. All these separate sets of tracings were performed in different random orders and the raters were blinded to the order. The degree of spatial overlap between the two traced ROIs for each participant was measured using the Dice similarity coefficient (Dice, 1945), which ranges for 0 (no overlap) to 1 (perfect agreement). The mean intra-rater Dice coefficient (+standard deviation) was 0.83 + 0.05 (range: 0.76 – 0.91), and the mean inter-rater Dice coefficient was 0.80 + 0.02 (range: 0.78–0.83), which both suggest strong agreement.

PAG ROI signal comparison

To compare the three different PAG localization techniques, we looked for non-uniformities in the extracted rs-fMRI brainstem signals from each ROI in healthy controls separately in our test (n = 15) and validation datasets (n = 15). A combination of FMRIB Software Library (FSL, ) and AFNI functions were used to transform PAG ROIs into the participants’ preprocessed rs-fMRI 3x3x3 mm grid-space and then to extract signals from the residual BOLD time series. The signal from the MNI-sphere ROI was directly extracted, while the MNI-trace ROI was output to a coordinate text file based in MNI space (AFNI’s 3dmaskdump). To transform the participant-trace ROI into the correct image space for signal extraction, a similar technique to the registration of the white matter mask to standardized participant specific space was applied. With this procedure, FSL’s FLIRT was used to register the output of the anatomical image from FreeSurfer and AFNI’s 3dNWarpApply, a nonlinear 3D warping technique, into the preprocessed rs-frmi image space. Extracted signals were averaged for further analysis.

We assumed the participant-trace ROI to be the “truest” measure of the anatomical location of the PAG because the ROI, a structure that is easily identifiable visually, is traced in anatomical space specific to the participant. The MNI-sphere and MNI-trace extracted signals were therefore each compared to the participant-trace signal using Pearson’s correlation (Fig. 1.B). Differences between the resulting correlation coefficients were determined after standardization using Fisher z-transformation with a paired t-test.

Whole-brain connectivity

Whole-brain connectivity analysis allowed us to identify quantifiable differences throughout the brain between the MNI-sphere and MNI-trace ROIs with reference to the participant-trace ROI. For this analysis, voxelwise signals were extracted and averaged in each of the Power 264 atlas regions (Power et al., 2011) from the residual BOLD time series output file from AFNI in healthy controls. Individual voxelwise signals were also extracted without averaging from a standard template brain-mask covering the entire brain, including areas of the brain not covered in the Power 264 atlas.

Connectivity error was first assessed in a test dataset (n = 15) and then validated in a matched dataset (n = 15). Connectivity error was defined as the standardized correlation value from each signal extracted from the Power 264 atlas regions to the participant-trace signal (z1, the “gold standard”) subtracted from the standardized correlation values from each signal extracted from the Power 264 atlas regions to the MNI-sphere signal (z2) or the MNI-trace signal (z3) (Fig. 2.A). The interquartile range (IQR) across participants and distance of the median (D) of the participants from zero were measures used to ascertain connectivity error. For example, in Fig. 2.B, smaller differences (i.e. smaller IQR and D) are seen in the connectivity error of the MNI-trace as compared to the MNI-sphere for an extracted signal (MNI: x  = 44, y = -53, z = 47) from the Power 264 atlas. In this case the connectivity of the MNI-trace ROI shows better performance in estimating PAG connectivity compared to the connectivity of the MNI-sphere ROI across participants for this particular region, as the MNI-trace ROI has a smaller connectivity error in relation to our reference trace, participant-trace ROI. A larger connectivity error (i.e. larger IQR and D), would indicate worse performance. In the voxelwise analysis, the same comparison was applied to individual voxels using the IQR and D for each extracted voxel instead of each Power 264 atlas region.

Fig. 2

Whole-brain functional connectivity analysis of three different localization techniques of the PAG. A. A measure of functional connectivity from each PAG ROI to 264 regions throughout the brain defined by the Power 264 atlas. An example of rs-fMRI BOLD residual time series signal from an atlas region and its connectivity measure to each PAG ROI. The connectivity error of the example region and the MNI-sphere ROI to the participant-trace ROI (z2-z1) and the connectivity error of the example region and the MNI-trace ROI to the participant-trace ROI (z3-z1) are calculated per healthy control participant. B. The distribution of the connectivity error as calculated in A (n = 15, test dataset). For this example region, smaller differences are seen in the MNI-trace ROI to the participant-trace ROI as compared to the MNI-sphere ROI to the participant-trace ROI. Signals from A come from the participant circled in purple. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Cluster-based analysis in UCPPS patients versus healthy controls

Due to the PAG’s association with chronic pain (Coulombe et al., 2016, Harper et al., 2018, Yu et al., 2014, Kutch et al., 2015, Segerdahl et al., 2018) and other chronic conditions (Linnman et al., 2012, Hamani et al., 2011, Harricharan et al., 2016), our objective was to lay the foundation for optimal localization of the PAG in UCPPS patients. Data from 100 UCPPS patients and 109 healthy controls was used for this analysis. A participant-trace is not practical for widespread use with more than a few participants. Our analysis was therefore focused on the two PAG ROIs traced in standard MNI space: the MNI-sphere ROI and the MNI-trace ROI. Connectivity was calculated by following a voxelwise correlation procedure established in AFNI (). Next, to control for site related differences, the participant data was demeaned for each MAPP Research Network site and the global mean for all sites was added. We first established healthy control connectivity and defined clusters from each PAG ROI, and subsequently used this connectivity as a backdrop against which to compare our findings in patients with chronic pelvic pain. Specific clusters were defined by their significant connection to each PAG ROI across healthy control participants (AFNI’s 3dClustSim) (p < 0.00001, , u = 6). The use of the most recent update of 3dClustSim for cluster analysis is taken into account after results in regard to false positives (Eklund et al., 2016). Each set of clusters was binarized and saved as an image file generated in AFNI’s interactive clustering interface accessed via AFNI’s GUI and based on parameters from 3dClustSim. The image file saved contains multiple masks created from clusters in the healthy control dataset (n = 109). The image file was then used to extract the average signals (AFNI’s 3dROIstats) from both patient and healthy control participants’ voxelwise whole-brain connectivity to one of the two PAG ROIs, respectively. AFNI’s CA_ML_18_MNIA atlas (Cox, 1996, Eickhoff et al., 2005) was queried (AFNI’s whereami) for cluster labels by centering on the voxel peak coordinates and location was checked with FSL’s Harvard-Oxford cortical and subcortical structural probabilistic atlases (). UCPPS patients and healthy control z-scores were statistically compared with a two-sample t-test for each cluster and significance was determined after p-values were corrected for multiple comparisons (false discovery rate using the Benjamini and Hochberg procedure) across each set of clusters.

Independent component analysis signal correction for potential physiological effects

Finally, the impact of physiological noise on whole brain connectivity differences needs to be addressed considering the localization of our PAG traces within the brainstem. In the absence of physiological recordings, we have implemented this correction using FSL’s ICA-FIX in an attempt to auto-classify signals and remove potential noise confounds due to physiological noise at the individual subject level (Salimi-Khorshidi et al., 2014, Griffanti et al., 2014). FSL-FIX trains features based on either the spatial map, time series, or frequency spectrum and provides a score of the likelihood of whether the feature is a signal (Bijsterbosch, J., Smith, S., Beckmann, C. An Introduction to Resting State fMRI Functional Connectivity. Oxford, New York: Oxford University Press;, 2017). A threshold is then used on these scores to classify which features are labeled as noise. In our analysis we followed this protocol, beginning with individually preprocessed, but not nuisance regressed rs-fMRI images. One of the key aspects in following FSL-FIX’s protocol is to train components classified as noise in the analyzed dataset. However, in order to maintain an automated preprocessing pipeline with as little experimenter input as necessary, we used the training data provided in FSL-FIX, “Standard.RData” (3.5 mm isotropic resolution, TR = 3 s), which was the closest match to our dataset’s parameters and has been implemented successfully in previous literature (Boyacioğlu et al., 2015, Marchitelli et al., 2016). For the reason that we did not create a unique training set specific to our dataset, we examined the impact of signal to noise threshold at a range of values suggested by the literature (Coulombe et al., 2016, Clemens et al., 2014, Descarries et al., 1982, Salimi-Khorshidi et al., 2014, Marchitelli et al., 2016). Following independent component analysis noise correction, resulting datasets were nuisance regressed (AFNI’s 3dTproject). The same methods as above were completed by 1) comparing the three different PAG localization techniques in healthy controls (test and validation) and 2) the cluster-based analysis in UCPPS patients versus healthy controls. These two analyses were chosen to evaluate potential differences in the direct PAG signal and to evaluate the sensitivity of healthy control versus patient differences, respectively.

Results

PAG ROI BOLD signal non-uniformities in healthy control participants

Our results included a total of 30 healthy control participants in our test (n = 15) and validation (n = 15) datasets (Table 1). To verify proper alignment between PAG ROIs traced in anatomical space and the extracted signal from functional space, alignment was inspected visually and quantified using the Dice coefficient. We found that the alignment between the normalized rs-fMRI and anatomical T1-weighted MRI was consistent across the healthy control participants in the test and validation datasets, showing no gross misalignment across participants (Supplemental Fig. 1). To test BOLD signal uniformities from each PAG ROI, the residual time series signal from the MNI-sphere and the MNI-trace ROIs was correlated to the signal from the participant-trace ROI, our best measure of the true PAG location. Fig. 1.C shows significant differences in the standardized correlation coefficient of the BOLD signal across participants in a paired t-test (test: p < 0.001, validation: p < 0.0001). This relationship demonstrates that the MNI-trace signal is significantly more similar to the participant-trace signal than the MNI-sphere signal (test: MNI-sphere z-score = 0.51 ± 0.17; MNI-trace z-score = 0.72 ± 0.04; validated: MNI-sphere z-score = 0.50 ± 0.13; MNI-trace z-score = 0.69 ± 0.08).

Whole-brain functional connectivity differences in healthy control participants

We evaluated whether the MNI-sphere and the MNI-trace signals performed better or worse by their connectivity error, or the standardized connectivity difference, to the participant-trace signal. Using the Power 264 atlas regions, 73% of regions (out of 264) showed better performance, or regions with smaller IQR and D when compared to the MNI-sphere connectivity error, in the test dataset for the MNI-trace connectivity error (Fig. 3.A.i). This was replicated in the validation dataset with 69% of regions and in the overlap of the test and validation datasets with 52% of regions with better performance (Fig. 3.A.i). Worse performance, or regions with larger IQR and D, was only seen in 3% of regions in the test dataset, 2% of regions in the validation dataset, and no regions in the overlap (Fig. 3.B.i) for the MNI-trace compared to the MNI sphere connectivity error. No particular network or pattern of regions emerged, though clearly the MNI-sphere and the MNI-trace connectivities differ in major ways throughout the brain.

Fig. 3

Performance of the PAG ROIs on whole-brain functional connectivity. A. Brain regions whose interquartile range (IQR) of the connectivity error was smaller and median closer to zero (D) in the MNI-trace ROI to participant-trace ROI compared to the MNI-sphere ROI to participant-trace ROI. B. Power 264 atlas regions whose interquartile range of the connectivity error were greater and median further to zero in their connectivity to the MNI-trace ROI to participant-trace ROI compared to the MNI-sphere ROI to participant-trace ROI. The comparison was performed at a for each Power 264 atlas region (A.i and B.i) and at each voxel within a brain-masked covering the entire brain (A.ii and B.ii). The analysis included the test dataset (n = 15), validation dataset (n = 15), and overlap of the two datasets. Key: SMM, somatosensory/motor (hand), SMMm, somatosensory/motor (mouth), COTC, cingulo-opercular task-control, Aud, auditory, DMN, default-mode, MR, memory retrieval, Vis, visual, FPTC, frontoparietal task-control, SN, salience, subC, subcortical, VAN, ventral attention, DAN, dorsal attention, CBL, cerebellum, Other, other regions.

Whole-brain voxelwise analysis showed similar results. Better performance of the MNI-trace connectivity error was seen with 72% of brain masked voxels (out of 73409) in the test dataset, 68% of voxels in the validation dataset, and 49% in the overlap (Fig. 3.A.ii). Worse performance was found in 2% of voxels in the test dataset, 2% of voxels in the validation, and less than 0.1% in the overlap (Fig. 3.B.ii). These results indicate smaller connectivity errors from the MNI-trace signal throughout the brain when compared to the MNI-sphere signal across healthy control participants.

Cluster-based functional connectivity differences in UPPS patients versus healthy controls

Our final goal examined not only whole-brain connectivity differences in choice of PAG localization, but also tested whether the choice would alter connectivity in determining patient versus healthy control differences in the MNI-sphere and the MNI-trace signals. The presence of both ROI’s in MNI space removed any restriction on participant numbers. Whole-brain functional connectivity differences in healthy controls suggest that the MNI-trace ROI was more similar, though not equivalent, to our participant-trace ROI. After neuroimaging preprocessing and quality controls, 209 participants (17.3% reduction from our initial number) were included in our final dataset (n = 100 patients, n = 109 healthy controls).

Following connectivity analysis and corrections for site effects we found significantly connected clusters in healthy controls from the MNI-sphere and the MNI-trace connectivity to the entire brain. We used these clusters to extract the average signal from each participant’s standardized functional connectivity data in both the patients and healthy controls. Table 2 shows all cluster labels defined using voxel peak coordinates except for the Right Posterior Cingulate for each PAG ROI. Peak coordinates for this cluster were not located in any atlas region and were defined based on neighboring voxels in the cluster. Also, we did not subdivide the largest cluster for either of the ROIs (MNI-sphere: 24,130 voxels, MNI-trace: 13,876 voxels) because our cluster defining primary threshold was already very low (p < 0.00001) (Woo et al., 2014) and we were able to acquire a number of regions throughout the brain for comparison. Collectively, there were smaller p-values from the group difference measured in patients versus healthy control clusters associated with the MNI-trace connectivity than those associated with the MNI-sphere connectivity (Table 2). The only clusters with significant differences after multiple corrections (p < 0.01, FDR q < 0.05) were those significantly connected to the MNI-trace signal, defined as the left rostral part of the inferior frontal gyrus (pars orbitalis), or the orbitofrontal area, and the left inferior parietal lobule, or the second somatosensory cortex (Fig. 4). Fig. 4.A shows that there were no clusters significantly connected to the MNI-sphere signal with significant differences in patients and healthy controls even in the inferior frontal gyrus (pars orbitalis) (Table 2, p = 0.42, FDR q = 0.51), a coinciding cluster significantly connected to the MNI-trace signal. To emphasize this difference in choice of PAG ROI, we compared the differences in connectivity between UCPPS patients and healthy controls within the orbitofrontal area (Supplemental Fig. 2). Significant differences (p < 0.01) are only seen in the connectivity to the MNI-trace ROI.

Table 2

The average connectivity of significant clusters found connected to the MNI-Sphere ROI and MNI-trace ROI in UCPPS patients and healthy controls. Cluster labels, Montreal Neurological Institute (MNI) for the peak voxel coordinates, total number of voxels in each cluster, average UCPPS patient z-score and standard deviation, average healthy control z-score and standard deviation, and p-value, uncorrected and false discovery rate corrected. AFNI’s CA_ML_18_MNIA atlas (Kutch et al., 2015, Hamani et al., 2011) was used for cluster labels based on their peak voxel coordinates (except for the Right posterior cingulate, whose label was defined using neighboring voxels within the cluster). Key: PAG, periaqueductal gray, SMA, supplementary motor area.

Average connectivity of significant clusters in UCPPS patients and healthy controls
	X (mm)	Y (mm)	Z (mm)	Volume (3 mm2)	patients (z-score)	healthy controls (z-score)	p-value/FDR corrected p-value
Cluster name
MNI-sphere
PAG ROI/midline	3	−27	−15	24,130	0.12 ± 0.07	0.13 ± 0.07	0.57/0.74
Right medial temporal pole	36	18	−42	44	0.08 ± 0.12	0.09 ± 0.12	0.35/0.74
Left medial temporal pole	−30	12	−42	40	0.09 ± 0.15	0.09 ± 0.13	0.68/0.74
Right precentral gyrus	42	−15	42	21	0.08 ± 0.15	0.09 ± 0.17	0.37/0.74
Right posterior cingulate	21	−39	21	20	−0.07 ± 0.14	−0.09 ± 0.15	0.73/0.74
Right postcentral gyrus	21	−33	60	18	0.08 ± 0.17	0.10 ± 0.17	0.74/0.74
Left calcarine gyrus	−15	−81	12	17	0.08 ± 0.16	0.09 ± 0.19	0.15/0.67
Left postcentral gyrus	−42	−21	39	11	0.07 ± 0.16	0.10 ± 0.18	0.14/0.67
Left inferior frontal gyrus (pars orbitalis)	−30	39	−15	8	0.08 ± 0.18	0.08 ± 0.15	0.53/0.74
MNI-trace
PAG ROI/midline	0	−33	−9	13,876	0.12 ± 0.09	0.14 ± 0.09	0.14/0.31
Right inferior frontal gyrus (pars triangularis)	51	24	30	210	0.06 ± 0.11	0.10 ± 0.12	0.03*/0.13
Left inferior frontal gyrus (pars opercularis)	−54	15	33	189	0.06 ± 0.11	0.09 ± 0.13	0.02*/0.13
Left middle frontal gyrus	−24	30	48	158	0.08 ± 0.13	0.10 ± 0.15	0.43/0.57
Right superior temporal gyrus	51	−18	6	151	0.07 ± 0.15	0.10 ± 0.13	0.19/0.36
Right precentral gyrus	45	−15	42	76	0.07 ± 0.14	0.10 ± 0.15	0.09/0.25
Left postcentral gyrus	−42	−21	42	63	0.09 ± 0.13	0.11 ± 0.16	0.26/0.50
Right posterior cingulate	12	−36	18	29	−0.07 ± 0.16	−0.11 ± 0.12	0.05*/0.17
Right superior parietal lobule	45	−48	57	26	0.05 ± 0.16	0.09 ± 0.16	0.07/0.23
Left precentral gyrus	−42	−6	51	18	0.08 ± 0.16	0.09 ± 0.15	0.58/0.64
Left superior frontal gyrus	−24	63	18	17	0.09 ± 0.15	0.09 ± 0.16	0.86/0.86
Right precentral gyrus	21	−30	63	16	0.09 ± 0.20	0.10 ± 0.18	0.71/0.75
Left SMA	0	18	63	13	0.09 ± 0.18	0.11 ± 0.20	0.57/0.64
Right SMA	6	24	66	11	0.05 ± 0.19	0.09 ± 0.17	0.10/0.25
Left inferior frontal gyrus (pars orbitalis)	−33	33	−9	10	0.03 ± 0.16	0.09 ± 0.14	0.004**/0.04*
Left inferior parietal lobule	−54	−39	36	7	0.02 ± 0.16	0.08 ± 0.15	0.004**/0.04*
Right superior frontal gyrus	24	63	12	6	0.07 ± 0.20	0.09 ± 0.16	0.40/0.57
Right superior frontal gyrus	18	21	69	6	0.07 ± 0.17	0.09 ± 0.17	0.41/0.57
Right SMA	6	9	75	6	0.09 ± 0.22	0.12 ± 0.22	0.45/0.57

Fig. 4

Functional connectivity differences in UPPS patients versus healthy controls. A. Selection of clusters based on significant voxelwise connectivity in healthy controls to the MNI-sphere ROI (HC: p < 0.00001, , uc = 6). B. Selection of clusters based on significant voxelwise connectivity in healthy controls to the MNI-trace ROI HC: p < 0.00001, , uc = 6). Uncorrected significant clusters identified between patients and healthy controls are circled in purple (HC vs. UCPPS: p < 0.05) and corrected significant clusters are identified with an asterisk (HC vs. UCPPS: p < 0.01, FDR q < 0.05). Key: HC, healthy control, UCPPS, patients. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Head motion

Head motion is known to be a confounding factor in the estimation of BOLD functional connectivity (Drysdale et al., 2017 Jan, Power et al., 2012, Van Dijk et al., 2012, Satterthwaite et al., 2012). Even small amounts of movement are known to bias functional connectivity measures, especially those regions found in the brainstem, and have been shown to increase relative to the distance between two brain regions of interest (Power et al., 2012). Our neuroimaging preprocessing measures attempt to correct these issues limiting framewise displacement by censoring volumes related to movement greater than 0.3 mm, setting this parameter in our AFNI pipeline according to previous publications (Plitt et al., 2015, Drysdale et al., 2017 Jan).

However, in our analysis, motion differences may contribute to differences between UCPPS patients and healthy controls. We used the same time series file for censoring movement volumes in our AFNI pipeline based on the Euclidean norm of the derivatives of the motion parameters for patients (n = 98, two files removed due to incomplete motion times series) and healthy control (n = 109). After comparing average movement across all time series per participant, we found no significant differences between these groups in overall head motion (p = 0.64).

Independent component analysis signal correction effect

To control for potential physiological effects from our PAG signal we implemented the auto-classification FSL-FIX at multiple signal to noise thresholds. First, we compared the three different PAG localization techniques across FIX thresholds 5, 40 and 75 in both test and validation datasets (Fig. 5.A). The higher the threshold indicates the inclusion of more potential noise components that are removed in the preprocessed rs-fMRI image. Our results show that significance was maintained after implementation between the standardized correlation coefficient of the MNI-sphere and MNI-trace ROIs to the participant-trace ROI. In fact, the p-value decreases at higher FIX thresholds (FIX Thr.40: p < 0.0001, FIX Thr.75: p < 0.0001). The average standardized correlation coefficient of the MNI-Sphere ROI to participant-trace ROI consistently remains lower in value compared to the MNI-trace ROI coefficient values.

Fig. 5

Independent component analysis (ICA) correction for potential physiological noise effects without ICA application and across FIX thresholds 5, 40 and 75. A. The differences in connectivity across healthy controls in test (n = 15) and validation (n = 15) datasets (No ICA: p < 0.001, FIX Thr.5: p < 0.001, FIX Thr.40: p < 0.0001, FIX Thr.75: p < 0.0001) B. A comparison of UCPPS patient versus healthy control connectivity differences to the MNI-trace ROI in significant clusters, left inferior frontal gyrus (IF) and left inferior parietal lobule (IP), from Fig. 4 (IF HC vs. UCPPS: No ICA: p < 0.01, FIX Thr.5: p < 0.01, IP HC vs. UCPPS: p < 0.01, FIX Thr.5: p < 0.01, IF HC No ICA vs IF HC FIX Thr.75: p < 0.01, IF HC FIX Thr.5 vs IF HC FIX Thr.75: p < 0.001). C. ICA noise component contribution at the MNI-trace ROI, IF, and IP, calculated from the maximum value at each location and averaged across participants. Threshold values are based on the number of participants with a signal above a voxel threshold of 0.1. Key: HC, healthy control, UCPPS, patients, IF, left inferior frontal gyrus, IP, left inferior parietal lobule.

Likewise, in our cluster-based analysis, in which we determine patient versus healthy control connectivity differences to the MNI-trace ROI, qualitative differences are maintained in significant clusters from our previous analysis without independent component analysis application (Fig. 5.B). These clusters, defined as the inferior frontal gyrus cluster (IF) and inferior parietal lobule cluster (IP), have consistently higher connectivity to the MNI-trace ROI in healthy controls than in patients across signal to noise thresholds. This same quality is not necessarily maintained in other clusters connected to the MNI-trace ROI. However, as more noise components are removed with higher thresholds, there is a loss in the significant differences between the groups.

Based on these changes in connectivity, we were interested in key noise components that influenced connectivity from the MNI-trace PAG ROI. The contribution from these potential noise components was calculated at the MNI-trace PAG ROI, inferior frontal gyrus cluster, and inferior parietal lobule cluster (Fig. 5.C). The maximal value at each location was used to determine the noise component averaged across participants. The voxel threshold across components is based on the normalized voxel z-scored values for the chosen maximum component (calculated by dividing the z-scores for that component by the maximum average value at the location, e.g. MNI-trace PAG ROI). Fig. 5.C shows the number of participants as a percentage (out of n = 209) whose signal is above a voxel threshold of 0.1. For example, at the location of the MNI-trace PAG ROI, there is an increased number of participants with a voxel threshold of 0.1 or greater around the area of the brainstem and cerebellum. There were fewer noise components removed at a FIX threshold of 5 across all participants than at higher FIX thresholds, 40 and 75.

Discussion

The goal of this study was to determine whether the type of localization methodology affected PAG functional connectivity and impacted the detection of differences in chronic pain patients compared to healthy controls. Using the time-intensive, hand-drawn, individual participant-trace ROI as a reference, we found significant non-uniformities in the extracted brainstem rs-fMRI signals from our MNI-sphere ROI and MNI-trace ROI in healthy controls. ROI localization methodology impacts functional connectivity findings throughout the brain using either region-based or voxel-based analysis. When tested in a larger population of healthy controls and UCPPS patients, we found significant differences in connectivity in regions previously examined in pain patients (Kong et al., 2010, Wei et al., 2016, Chen et al., 2017, Napadow et al., 2010 Aug, May, 2008, Apkarian et al., 2004). To the best of our knowledge, this is the first study to show that such a comparison depends on PAG localization methodology. How the PAG is localized as a whole, impacts findings of connectivity to other regions of the brain and should be considered when analyzing disease-states associated with PAG connectivity.

This study establishes that the PAG localization methodology frequently reported in literature (Kong et al., 2010, Wei et al., 2016, Wei et al., 2017, Harper et al., 2018, Chen et al., 2017, Yu et al., 2014, Li et al., 2016, Liu et al., 2017, Linnman et al., 2012), as representing the entire PAG region and labeled in this study as the MNI-sphere, does not adequately represent the anatomical PAG. This finding leads us to ask what brainstem portion the MNI-sphere may represent. The MNI coordinates reported in literature (MNI: 4 –26 −14) align with another study (Mayberg et al., 2000) whose peak voxel coordinates from a positron emission measurement relate to improvement in depression symptoms after treatment of a selective serotonin reuptake inhibitor. The dorsal raphe nucleus (DRN), a region in the brain with the largest population of serotonergic neurons (Silva and McNaughton, 2019, Hornung, 2003, Michelsen et al., 2008), sits very close to the ventrolateral PAG in the brainstem (Hornung, 2003, Michelsen et al., 2008, Naidich et al., 2009). The DRN is also involved in descending pain modulation (Stamford, 1995, Babalian et al., 2019) and, although it is histologically different from the PAG (Silva and McNaughton, 2019, Descarries et al., 1982), produces analgesia when stimulated (Stamford, 1995, Liebeskind et al., 1973, Oliveras et al., 1974); similarly, the ventrolateral PAG has been stimulated to evoke opioid mediated analgesia (Linnman et al., 2012, Kong et al., 2010, Bandler and Shipley, 1994, Tracey et al., 2002). The MNI-sphere in the literature could therefore be identifying the DRN. Nevertheless, fMRI coordinates reported in literature for the DRN region are ambiguous (Napadow et al., 2005, Pedroni et al., 2011, Beliveau et al., 2015), often cited as part of the anatomical location of the PAG and it is unclear if these regions are functionally distinct (Bandler and Shipley, 1994, Stamford, 1995, Liebeskind et al., 1973).

Clearly the optimal ROI would be defined both anatomically and functionally from within each participant (Poldrack et al., 2011, Nieto-Castañón and Fedorenko, 2012). Though not interchangeable, we demonstrated that the participant-trace and MNI-trace ROIs were similar in their connectivity patterns. Caution should still be taken when considering the use of the MNI-trace ROI as a PAG ROI because the trace is derived from an average of 152 T1-weighted structural images. For example, if the study of interest includes variability between participants, such as size and shape variations of the PAG, there is a likelihood of mixed functional activity patterns (Eickhoff et al., 2018). Our results showed no Power 264 atlas regions and less than 0.1% of voxels with consistent error differences of the whole-brain functional connectivity of the MNI-trace ROI compared to the MNI-sphere ROI in our test and validation datasets (Fig. 3). This was indicative of a better estimation of the MNI-trace ROI to our participant-trace ROI using participant variation as a measure of differences (interquartile range and distance of the median from 0). Even so, the method for defining the “gold standard” remains complicated for extracting the BOLD time-series signal from the entire PAG.

Characterizing a region in the brainstem using fMRI has been met with difficulty (Tracey and Iannetti, 2006, Brooks and Faull, 2013, Sclocco et al., 2018). Recent efforts focus on defining a brain region first, before testing its functionality with different experimental paradigms (Genon et al., 2018). Traditionally brain regions were defined by stimulation of the region of interest or lesion-deficient approaches, such as with the functional segregation of the PAG in humans (Linnman et al., 2012, Coulombe et al., 2016, Genon et al., 2018, Bittar et al., 2005, Richardson and Akil, 1977). Heat stimulation tasks are the main choice for eliciting a PAG response for measuring pain-induced activation with fMRI or positron emission tomography (Linnman et al., 2012). However, one behavioral outcome may be related to multiple neuronal behaviors (Genon et al., 2018, Price and Friston, 2002) and the behavior elicited to identify a region is not specific to what that region actually does. We defined the PAG by expertly hand tracing the grey matter voxels in the brainstem of a T1-weighted structural MRI in patient specific space, labelled participant-trace, or from the MNI152 1 mm standard space, labelled MNI-trace. This process can be time consuming and details about distinct subregions are unidentifiable (Tracey and Iannetti, 2006, Sclocco et al., 2018). An alternative to our method is the creation of a probabilistic grey matter map for the study group of interest encompassing the entire PAG (Brooks et al., 2017). For the definition of specific functional subgroups within the PAG, a combination of imaging approaches and atlas references can be used to maximize spatial resolution and compare intrinsic functional connectivity within the region (Coulombe et al., 2016, Tillman et al., 2018). A normalization approach optimized for the brainstem as to reduce further sources of potential error with PAG definition (Diedrichsen et al., 2011) is an additional consideration. Overall, the integration of different approaches by pooling data from multiple sources, meta-analysis (Linnman et al., 2012, Poldrack et al., 2011, Costafreda, 2009), and established databases, such as Neurosynth (Yarkoni et al., 2011) and BrainMaps (Mikula et al., 2007), rather than focusing only on single study activations to define regions, provides a more reliable brain definition.

The purpose of this study was not to provide a specific method for analyzing the PAG, but to highlight the erroneous application of published PAG coordinates in chronic visceral pain. Nonetheless, the MNI-trace ROI provides a potential representation of the PAG. Anatomical PAG connections to the cortical and subcortical parts of the brain are well established. Connections to the parts of the medial prefrontal, anterior cingulate, dorsal medial, orbital, posterior cingulate, temporal, ventral insula, and amygdalar regions of the brain were revealed by anterograde and retrograde tracers in Macaques (An et al., 1998). Human studies have shown similar connections using deep brain stimulation (DBS) and diffusion tractography depending on the seeded region; though not all connections fully correlate to animal studies (Benarroch, 2012, Ezra et al., 2015) and differences exist across species (Ezra et al., 2015). It is important to keep in mind the BOLD signal provides only a surrogate marker for the actual neuronal brain responses (Logothetis, 2008) and since functional connectivity depends entirely on the BOLD signal, it cannot demonstrate actual anatomical connections (Linnman et al., 2012). In addition, evaluating functional connectivity throughout the entire brain produces a high family-wise error rate causing some important connections to be overlooked (Poldrack et al., 2011), especially when dealing with small volume regions. By limiting our analysis between patients and healthy controls to predefined regions connected to our ROIs within our healthy control dataset, we were able to distinguish important differences in two clusters of the MNI-trace ROI associated with PAG connectivity that the MNI-sphere ROI failed to identify. While this approach may bias the connectivity differences towards healthy controls, most cluster differences do not pass multiple comparisons and only two in the MNI-trace ROI: inferior frontal and inferior parietal, were significant (assessed in with both one and two-tailed, two sample t-tests).

The first cluster was in the orbitofrontal cortex, a region known for emotionally driven cognitive states and its impairment in chronic pain (Chen et al., 2017, May, 2008, Apkarian et al., 2004, Babalian et al., 2019). Connectivity from the orbitofrontal cortex to the PAG is established in humans in pain-related studies (Coulombe et al., 2017, Wager et al., 2007) and was shown within both the MNI-sphere and MNI-trace ROI cluster-based connections from our voxelwise analysis (Table 2); however, significant differences between UCPPS patients and healthy controls were only seen in the MNI-trace ROI. By focusing on the orbitofrontal area, we show that patient versus healthy control differences are not simply due to healthy control bias from cluster definition and the resulting differences between PAG ROIs are a feature of ROI location (Supplemental Fig. 2).

The second cluster occurred in the inferior parietal lobule. Despite this region not showing a direct connection in animal studies to the PAG (Öngür and Price, 2000), there is minor evidence of its linkage in rat studies (Herrero et al., 1991) and direct linkage to the dorsolateral prefrontal cortex and its association with visceromotor control (Öngür and Price, 2000, Quintana and Fuster, 1999). In patients with chronic inflammation, the inferior parietal lobule has recently been found to be associated with the magnitude of inflammatory response (Schrepf et al., 2018). Additionally, a fibromyalgia study (Napadow et al., 2010 Aug) between patients and healthy controls showed resting-state network connectivity differences in the secondary somatosensory cortex of the default node network (DMN), a region located very close to the parietal area reported in our study. Another study (Kucyi et al., 2013) found that the secondary somatosensory cortex and the temporoparietal junction were activated during an attention to pain fMRI task. Though these studies validate our patient and healthy control comparison for the MNI-trace connectivity, further studies are needed to understand its relationship as a functional measure for the PAG.

It is important to note that these connectivity findings fit well with a larger view of the PAG’s role in chronic pain states as an orchestrator of the acute threat response (Benarroch, 2012, Mobbs et al., 2007). Prior studies comparing irritable bowel syndrome with ulcerative colitis have demonstrated that frontal connectivity with PAG may determine whether the PAG is “set to sound the alarm” (Mayer et al., 2005) which would generate both an autonomic sympathetic response and a conscious pain experience. A volumetrically larger PAG primarily characterizes women with endometriosis who do not have pain, compared to those who do (As-Sanie et al., 2012), confirming that PAG function and its trophic state may play an important role in chronic pelvic pain. It is also known that vagal parasympathetic influence is reduced in both chronic pain and inflammatory states (Williams et al., 2015, Koenig et al., 2016), and the vagus nerve plays a major role in the control of inflammation (Pavlov and Tracey, 2006). The PAG may carry out its primary function of orchestrating a response to an acute threat by linking cortical processes focused on the readiness of bodily systems such as the inflammatory response, with descending vagal and sympathetic influences which carry out these tasks. Further, the PAG could play a central role in a chronic pelvic pain syndrome, if the systems designed to reset it after the acute threat response fail to do so, resulting in a state of “chronic threat response”, and consequent exhaustion of resources designed only for an acute need.

Unlike the MNI-sphere ROI, our MNI-trace and participant-trace ROIs surround the cerebral aqueduct which carries cerebrospinal fluid (CSF) between the third and fourth ventricles (Bandler et al., 1991). One of the challenges in fMRI of the brainstem is spatial resolution at 3 T, which shows the brainstem as a mostly homogeneous structure (Sclocco et al., 2018). The anatomical size of the PAG is 4–5 mm around the aqueduct (Linnman et al., 2012, Sclocco et al., 2018), therefore the MNI-trace and participant-trace likely include artifacts from the CSF due to cardiac and respiratory effects (Brooks and Faull, 2013, Sclocco et al., 2018). The PAG’s association with autonomic control suggests that these movement artifacts could also be temporally correlated to behavior (Linnman et al., 2012, Benarroch, 2012, Ezra et al., 2015, Tracey and Iannetti, 2006, Napadow et al., 2005, Roberta et al., 2016, Chang et al., 2013). The significant differences in UCPPS patients and healthy controls in our results (Fig. 4) could be due to autonomic functional differences or coupled physiological noise effects. We recognize that the p-values across a number of the clusters in MNI-trace ROI appear smaller than in the MNI-sphere ROI, which may be due to proximity of this region to the cerebral aqueduct and the effect of cardiac pulsation on the CSF (Salimi-Khorshidi et al., 2014, Brooks and Faull, 2013, Griffanti et al., 2017, Beissner, 2015). The variance in global signal, which was not corrected for in our pipeline because of distant dependent bias (Jo et al., 2013), is linked to both fluctuations in cardiac rate and respiratory activity (Liu et al., 2017), and could account for the relationship in patient versus control differences in the MNI-trace ROI connectivity to the rest of the brain.

Even though we have not directly controlled for physiological measures such as heart rate and breathing, we utilized an independent component analysis approach in our preprocessing pipeline to address the potential effect. At increasing FIX thresholds, we show that the signal differences from ROI placement exist regardless of ICA correction (Fig. 5.A). However, significant differences between patients and healthy controls were lost at higher signal to noise FIX thresholds (Fig. 5.B). This was likely due to potential noise around the aqueduct, caused by fluctuations in cardiac and respiratory activity on the CSF, and arterial blood flow from the posterior cerebral artery (Salimi-Khorshidi et al., 2014) near the PAG (Fig. 5.C). The spread across the PAG at higher FIX thresholds in Fig. 5.C could also be enhanced because of spatial smoothing (Beissner, 2015). In lower spatial resolution datasets, such as in our study, signal can spread to nearby veins in brainstem and cerebellar areas (Griffanti et al., 2017) and should not be removed from the data. Nevertheless, we also observe a maintenance in the qualitative relationship between the connectivity in healthy controls and patients (a quality not maintained in other clusters connected to the MNI-trace ROI) and a valid significant difference in connectivity at a more conservative FIX threshold. When using automatic classification methods without the input from expert training on noise component selection, more conservative FIX thresholds are recommended (Salimi-Khorshidi et al., 2014) so as not to remove signal. We demonstrate that FSL-FIX’s independent component analysis may have an important impact in interpreting the relationship of the PAG to other areas of the brain if physiological measures have not been recorded directly. This automated technique removed spurious correlation artifacts that affect the PAG connectivity and may constitute potential physiological noise. In moving forward with further analysis on the PAG region, use of FSL-FIX is strongly recommended. According to FSL-FIX documentation, more selective approaches in defining ICA components should also be considered, such as creating a training set for selecting heart rate and breathing components by several trained experts as defined by Salimi-Khorshidi and colleagues (Salimi-Khorshidi et al., 2014).

In addition to potential physiological noise artifact we reduced motion through our AFNI pipeline (AFNI’s ANATICOR) which accounts for local artifacts rather than using the global signal (Cox, 1996, Jo et al., 2013). Even so, we have found no significant differences in overall head motion in our patient and healthy control populations (p = 0.64). Given that the PAG’s location is so close to such a physiologically dynamic structure affected by heart rate and breathing, and dealing with signals derived from this activity, the aqueduct must be considered carefully, and these findings show that caution should be taken when analyzing PAG functional connectivity. Many prior studies did not address these potential confounds (Kong et al., 2010, Wei et al., 2016, Wei et al., 2017, Chen et al., 2017, Yu et al., 2014).

This analysis of PAG localization demonstrates the dependence of disease-state differences on ROI localization methodology in functional neuroimaging studies, particularly in an area such as the brainstem. Improved imaging modalities will likely lead to better understanding of both anatomical morphology and functional connectivity of the PAG. For example, enhanced resolution of the GE 7 T MR950 MRI at the Center for Imaging Research at the Medical College of Wisconsin improved anatomic identification of the PAG borders and facilitated identification of subregions within the PAG during aversive image viewing (Satpute et al., 2013). However, this particular imaging device suffers from blurring in areas of heterogeneous density weakening its capacity for functional imaging particularly in the frontal areas where significant fluid-bone-air interfaces occur. Ultimately, it will be critical to corroborate PAG localization methods anatomically when performing neuroimaging meta-analyses in order to account for activation peaks found outside the PAG region (which weaken spatial specificity), to align processing procedures, and to limit mislabeling (Linnman et al., 2012).

In summary, this study evaluates an anatomical PAG localization method directly from structural MRI images. Future PAG studies may benefit from this foundational work. Based on this study, we recommend anatomic rather than functional localization procedures, either through hand-traced participant specific ROI, or a PAG traced in average space, depending on the research problem. As a note of caution, we have supported use of the MNI-trace ROI based on our connectivity results corroborated by prior literature and our similarity measures to the participant-trace, differences between UCPPS patients and healthy controls alone do not determine choice of ROI. Interest in the PAG will likely increase in reference to chronic pelvic pain and bladder pain, as the PAG has been shown to activate in response to bladder filling and during micturition (Benarroch, 2012, Fowler et al., 2008, Blok et al., 1997). The MAPP research group has shown that painful bladder filling or urgency is associated with a possible subset of UCPPS patients (Kleinhans et al., 2016, Henry et al., 2015). The present MNI-trace ROI can be used to compare changes in PAG activity in UCPPS patients and healthy control participants and their relationship to pelvic pain. The MNI-trace ROI may also be useful in predicting longitudinal symptom evolution in UCPPS patients, since whole-brain functional connectivity patterns of the MNI-trace ROI signal are similar to our reference participant-trace ROI across defined Power 264 atlas regions. Connectivity from the MNI-trace ROI may link to individual differences in pain symptoms in the larger UCPPS population. Given these applications, there is still a need for further study to establish the reproducibility of the functional connectivity of the PAG in humans.

Conclusion

To determine the optimal localization of the PAG ROI in fMRI, we compared the rs-fMRI functional connectivity in three different traces of the PAG: 1) traditional MNI-sphere, 2) MNI-trace, and 3) participant-trace. The direct fMRI time series signal from the MNI-trace ROI paralleled the participant-trace ROI, our presumptive reference measure of the PAG signal, far more closely than the traditional MNI-sphere. The same finding emerged from further investigation using whole-brain functional connectivity to these ROIs and likewise, found that the connectivity to the MNI-trace ROI was more similar to the participant-trace ROI. In our second analysis, we compared whole-brain voxelwise connectivity to the MNI-sphere and MNI-trace ROI, respectively in patients versus healthy controls. Regions established to be connected to these ROIs in healthy controls only showed patient versus healthy control differences in the MNI-trace ROI and not in the traditional MNI-sphere ROI. Our results confirm that the PAG ROI, MNI-sphere, often reported in literature as the PAG is not equivalent to the anatomical trace of the PAG (both the MNI-trace and the participant-trace), and may reflect another region in the brainstem related to pain, perhaps the dorsal raphe nucleus. This study demonstrates that care and validation are required when defining the PAG, especially when probing for disease-state differences. Assumptions about localization may alter our understanding of neuronal processing and influence prospective clinical conclusions. The connectivity findings reported here support a specific hypothesis linking the PAG’s role as orchestrator of the response to threat and a potential neuronal infrastructure for the development and maintenance of a chronic pelvic pain syndrome and its autonomic co-morbidities (Chelimsky et al., 2019).

CRediT authorship contribution statement

Sonja Fenske: Data curation, Investigation, Methodology, Validation, Formal analysis, Software, Writing - original draft. Douglas Bierer: Data curation, Methodology, Writing - review & editing. Gisela Chelimsky: Resources, Supervision. Lisa Conant: Formal analysis. Candida Ustine: Formal analysis. Ke Yan: Formal analysis. Thomas Chelimsky: Supervision, Conceptualization, Resources, Methodology, Writing - review & editing. Jason J. Kutch: Supervision, Conceptualization, Data curation, Methodology, Resources, Writing - review & editing, Project administration, Funding acquisition.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.