Prospective participants were recruited through: (1) the Sexual Health Research Laboratory (SHRL) database, consisting of past participants who consented to being contacted for future research, (2) advertisements posted in the Kingston community including the university campus, (3) social media (e.g., Facebook) advertisements, (4) pamphlets in doctors' offices, and (5) emails to health care providers who may provide care for women with vulvodynia (e.g., urologists, gynecologists, pelvic health physical therapists). This study was approved by the Queen's University Health Sciences Research Ethics Board and informed consent was obtained from each participant.

Participants were screened for eligibility either over the phone or online using Qualtrics survey software. Thirty-two women between 18 and 48 years of age were recruited for this study, 16 women in each group (PVD, Control). Significant effects have been found in previous studies that utilize similar analyses in populations with and without chronic pain conditions with similar sample sizes (19, 25). Demographic information for these participants can be found in Table 1. To reduce the influence of age-related changes or hormone use confounding the results, groups were matched on age (±5 years) and hormonal contraceptive use (yes or no) (26–28). Inclusion criteria for both groups included: older than 18 years of age or under 50; no magnetic resonance imaging (MRI) contraindications (e.g., metal implants); no major brain or spinal cord injury; no use of medications that substantially affect the central nervous system (e.g., antipsychotics); and not currently pregnant or breastfeeding. In addition, women with PVD had to report idiopathic, provoked pain to the vaginal entrance during sexual and non-sexual activities involving vaginal penetration. The results of the screening process and a gynecological examination (see below) formed the basis of inclusion or exclusion of women with PVD. For Control women, the primary inclusion criterion was the absence of self-reported chronic vulvar pain. Figure 1 illustrates the study procedure and the number of participants who completed each phase.

Vulvodynia is a diagnosis of exclusion (1). Therefore, the study gynecologist ruled out factors associated with vulvar pain (e.g., infections) in women believed to have PVD. The internal and external genitals and reproductive organs were visually and manually examined by the gynecologist. The gynecologist also performed the cotton-swab test of the external genitals (i.e., labia majora and minora, vestibule, perineum)—the main gynecological test for the diagnosis of PVD (29). The cotton-swab test involved the palpation of the labia majora, inner labia minora, midline areas, and six randomly ordered sites at the vestibule (e.g., 1 o'clock, 4 o'clock, 6 o'clock, etc.). Participants were asked to rate the intensity of the pain on a numerical rating scale ranging from 0 (no pain at all) to 10 (worst pain ever felt) after each palpation. In order to be included in the PVD group, a mean average pain intensity rating of at least 1/10 during the cotton-swab test of the vestibule was required. However, all participants in this study reported pain in four or more locations. Participants with PVD reported an average pain intensity of 4.61 (SD = 2.00) across all vestibule sites.

Trait Subscale of the State-Trait Anxiety Scale (STAI). The STAI (30) is a 40-item measure of self-reported state (present) and trait (stable) anxiety. Greater scores on either 20-item subscale indicate greater feelings of anxiety. The internal consistency of the STAI has been found to be acceptable (30). In the present study, the trait subscale was administered to participants.

The Beck Depression Inventory- II (BDI-II). The BDI-II (31) consists of 21 self-report items assessing depressive symptom severity over the past 4 weeks, with greater scores indicating increased severity. The BDI-II has been shown to have excellent internal reliability and to have convergent validity with measures of depression (32).

The thermal pain testing and sham MRI sessions were completed at the Queen's MRI Facility Sham Room. The Sham Room is very similar to the actual MRI environment so that participants are able to experience undergoing an MRI, without the magnetic field. Participants were familiarized with the protocol during training on the thermal pain testing protocol and runs in the mock scanner, with the goal of reducing anxiety and providing more consistent functional magnetic resonance imaging (fMRI) results across repeated runs.

Participants were exposed to noxious heat stimulation by resting their right hand on the MRI-compatible robotic contact-heat thermal stimulator (RTS-1; Spinal Map Inc., Kingston). During the thermal pain testing session, a level of temperature for each participant that was considered “moderately painful” was determined on a 101-point scale with verbal descriptors in increments of 10 (i.e., 0 = no sensation, 10 = warm, 20 = a barely painful sensation, 30 = very weak pain, 40 = weak pain, 50 = moderate pain, 60 = slightly strong pain, 70 = strong pain, 80 = very strong pain, 90 = nearly intolerable pain, 100 = intolerable pain). If participants did not rate their pain intensity as “moderately painful” before reaching the safety limit of 51°C, then 51°C was used during trials. A thermode was raised or lowered to contact the thenar eminence (corresponding to the C6 dermatome) through a cut-out in a plexi-glass case, with the timing and temperature controlled by custom-made software on a laptop. The first three contacts lasted for 1.5 s at 45°C with onsets every 3 s. Participants rated the intensity of each individual sensation. This process was repeated with the temperatures of 46 and 47°C after 2 min of rest to avoid sensitization. The goal of this task was to train participants to rate their sensory experience. Following this task, participants rated the intensity of their sensations to 10 consecutive contacts of different temperatures (i.e., 46, 50, 44, and 48C), followed by 2 min of rest. Contacts were applied for 1.5 s, with onsets every 3 s, in order to evoke temporal summation of second pain. Repeated applications of brief stimuli can evoke temporal summation of second pain, which is mediated by C-fibers, and provides a robust BOLD response for fMRI studies of pain (25, 33–35). This protocol provided a measure of each participant's pain sensitivity and allowed us to select the appropriate temperature for stimulation during the fMRI session. The thermal pain testing and Sham session lasted ~1 h.

Imaging was conducted in a research-dedicated 3 tesla MRI system (Siemens Magnetom Trio). The imaging session took 1.5 h and included one anatomical scan and ten 4.5-min functional scans, guided by initial localizer scans. Participants were positioned with padding and blankets for comfort and to prevent movement. Posterior head, neck, and spine receiver coils were used for detecting the MRI signal, and a mirror was positioned so that participants could view the intensity rating scale and instructions on a rear-projection screen. The peripheral pulse was recorded to account for physiological noise (20, 21).

Localizer images were acquired in three planes for subsequent slice positioning. In order to obtain high quality images with good spatial fidelity, functional MR image data for BS/SC were acquired from the same region as the anatomical scans with a half-fourier single-shot fast spin-echo (HASTE) sequence with 1.5 × 1.5 × 2 mm resolution, and an echo time (TE) of 76 ms (20). Nine sagittal slices spanned from below the disc between the first and second thoracic vertebrae (T1/T2) to above the corpus collosum. This method has been demonstrated to provide optimal T2-weighted BOLD sensitivity and image quality in the spinal cord and brainstem (19–21, 25). HASTE images are less sensitive to field inhomogeneities caused by the anatomy surrounding the spinal cord. It has been shown that the majority of BS/SC fMRI studies to date have utilized variations of single-shot fast spin-echo imaging methods (36). Image data were acquired with a repetition time (TR) of 6.75 s, so each fMRI run consisted of acquiring 40 volumes to describe the BOLD time-series responses.

At the start of each acquisition participants were informed that a new run was about to begin but they did not initially know if they would receive a painful stimulus (“Pain” condition) or not (“No Pain” condition). One minute after the start of the acquisition, they were informed whether or not they would feel a heat stimulus. For runs without stimulation, data acquisition was continued for a total of 4.5 min, with no stimulus applied. In the “No Pain” condition, participants were instructed to remain still for the duration of the run. If they were to feel a stimulus, it was applied beginning 2 min after the start of the acquisition (i.e., 1 min after they were informed of the stimulus). It consisted of 10 heat contacts applied to the heel of the right thumb, with contact durations of 1.5 s and onsets every 3 s, at the temperature that elicited their moderate level of pain (Figure 2). Pain evoked by repetitive heat stimuli has been shown to result in C-fiber-mediated summation of pain (19). The stimulation period was followed by 2 min without stimulation, totaling 4.5 min of acquisition. Participants were instructed to silently rate their pain for each contact, as practiced during sham training, and at the end of each run were asked to report their ratings for the first and last contacts. If, during the MRI session, participants began to rate their pain as less or greater than moderately painful (either after leaving the sham room or across the MRI session runs), the researcher would titrate the temperature to achieve a moderately painful stimulus. Before beginning the MRI acquisition, after the 5th trial, and following the MRI acquisition participants were asked to verbally “rate your anxiety at this moment” on a 0 (I feel calm and at ease) to 10 (I feel as if I may have an anxiety attack) scale (Table 2). There were a total of five runs for each stimulus type (Pain or No Pain) presented in a pseudo-randomized, counterbalanced order with 2-min breaks to avoid sensitization of the skin.

Participants with PVD who passed the initial screening were scheduled for a gynecological examination to confirm their eligibility for inclusion in this study. Women with PVD gave informed consent at the gynecological exam and again at their MRI session, and Control women gave informed consent at their MRI session. They were asked to refrain from alcohol for at least 12 h prior to the imaging portions of the study and from caffeine for at least 6 h before the study, as these substances may affect mood, alertness, and the BOLD response. Participants were also asked to eat a regular meal at their usual mealtime prior to participation in the study.

Before the imaging session, the thermal pain testing/sham MRI session was conducted. Before entering the magnetic imaging environment, participants again confirmed that they had no contraindications for the subsequent magnetic resonance imaging. Participants were asked to change into MRI-safe clothing provided by the facility if any item of clothing possibly contained materials that were not MRI safe. They were then positioned in the MRI for scanning to commence. At the end, participants were debriefed and compensated. This study was part of a larger study that involved a separate session of brain imaging (37), saliva collection, and the administration of validated self-report measures; these additional components are not discussed in this paper.

The BS/SC data consisted of 5 runs in each of the No Pain and Pain conditions. Pre-processing was conducted with custom-written software in MATLAB called spinalfmri8 (25). First, image data were converted to Neuroimaging Informatics Technology Initiative (NIfTI) format. Sagittal slices were then co-registered to correct for bulk body motion using a 3D, non-rigid registration tool [Medical Image Registration Toolbox; (38, 39)]. These co-registration motion parameters were later used to model bulk movement and reduce noise. Next, the images were spatially normalized to fit an anatomical template previously defined based on 356 healthy participants in previous studies (23, 25). To avoid variability in T1-weighting, the first two volumes were discarded from each time series. Regions of white matter were used to model global noise. Models of noise based on the peripheral pulse, bulk motion parameters, and white matter regions, were fit to the data using a general linear model (GLM) and subtracted from the measured data.

Custom-made software written in MATLAB was used to perform SEM analyses for the purposes of assessing functional connectivity. SEM analysis, adapted for our purposes, provides measures of connectivity between regions, based on a pre-defined anatomical model of regions which are known/suspected of playing a role in pain processing, and anatomical connections between regions (40). This method allows us to model multiple sources of input signaling to each region, as opposed to using simple correlations between BOLD responses, which assumes that each region has a single dominant source of input. The additional information provided by an anatomical model also enables the directionality of interactions to be inferred. Based on previous findings of BOLD variations occurring in the spinal cord and brainstem in the periods before and after painful stimulation, and also in the absence of any sensory stimulus, pain has been shown to be continuously modulated by emotional and cognitive factors (23, 41–43). As a result, instead of modeling the influence of continuous cognitive and emotional influences, the approach to analysis used in the present paper is data driven, which has been shown to be appropriate for pain populations (44). The a priori specified model we used (Figure 3) is based on anatomical connections and regions known to be related to pain processing, and includes the thalamus (Thal), hypothalamus (Hyp), periaqueductal gray matter (PAG), locus coeruleus (LC), parabrachial nucleus (PBN), nucleus raphe magnus (NRM), nucleus gigantocellularis (NGc), nucleus tractus solitarius (NTS), dorsal reticular nucleus (DRt), and C6 right dorsal portion of the spinal cord (C6RD) (22, 45–47).

The SEM analysis method consisted of first identifying voxels in the pre-selected regions (based on the anatomical model). K-means clustering was then used to divide each region of interest into seven sub-regions (clusters) based on their BOLD responses to provide finer anatomical precision for identifying connected regions. Voxels were divided and clustered based on functional characteristics, thereby increasing spatial precision. Any clusters which were predominantly along the edge of cerebrospinal fluid (CSF) –filled spaces were removed from the analysis, because these are susceptible to effects of physiological movement. The BOLD time-course responses in each region (as a target for input signaling) was then fit to a linear combination of the BOLD time-course responses in the possible source regions. Region definitions are performed once and are the same across groups and conditions. The linear weighting factors which quantify each region's contribution as an input signal source is termed a β value. Thus, for each target region, a general linear model (GLM) was used, with the source region time-courses as the basis set, and the set of linear weighting factors (β) was calculated. This procedure was repeated for every possible combination of clusters in each source and target region and was repeated for each run and group type. We also allowed the SEM weighting factors to vary dynamically by calculating connectivity strengths using only data from specific epochs, spanning seven volumes, one within the baseline period of a run (75–120 s), one spanning the stimulation period (115–160 s), and one within the baseline period after the stimulation period (153–198 s). Volumes from the selected epoch were concatenated across runs of the same type in each person.

As mentioned above, every combination of sub-regions within each region were investigated to determine the best fit with the BOLD time-series data. The linear weighting factor value for each connection within a network component (i.e., a target cluster and a set of possible source clusters) was thus calculated multiple times with different combinations of source regions included. We calculated the goodness-of-fit by examining the proportion of the variance in the target region that was accounted for by the network component fit (i.e., R² value). We determined the probability of each fit occurring by random chance by converting R-values to Z-scores (Fisher's transform; Z = tanh⁻¹ (R)tsize-3) to estimate the significance of the fit. The fit was repeated by pruning a single source region at a time from the network component to assure all terms accounted for a significant portion of the variance in a target region, using an F-test with a corrected cutoff F-value corresponding to a p < 0.05. If any portion of the network did not account for a significant portion of the variance, it was not included in the results.

The significance of each total network component was determined by previously established probability distributions (40). The distributions depend on model parameters, and the significance thresholds accounted for family wise error rate (p_fwe < 0.05), taking into consideration the number of network combinations that were tested across the possible combinations of sub-regions. The complete best-fit network was identified for each run type for both PVD and Control women by identifying the combinations of source regions that have the highest Z-values for each network component (i.e., the set of clusters with the highest Z-scores for each network component). The best-fit model explains the BOLD signal variations in the regions of interest based on known anatomical connections.

Significant differences in connectivity between regions were identified between study conditions using paired sample t-tests. Significant differences were inferred between two conditions if the t-test culminated in a Bonferroni-corrected p < 0.05, accounting for the total number of comparisons. Connections selected for further analysis based on having significant fits (p_fwe < 0.05) in any of the four combinations of study groups and conditions. Variations in SEM weighting factors between the study groups (PVD or Control), and between conditions (Stimulation or No-Stimulation) were also investigated by means of an analysis of variance (ANOVA). The family-wise error rate correction was applied to the ANOVA. Individual SEM results (β-values) from all participants for a given connection were then analyzed with a 2 × 2 ANOVA with the study group and the study condition as the two independent variables. Results were identified with significant main effects of the group, or condition, or group x condition interactions, at p < 0.05.

Table 1 summarizes the characteristics of the participants in each study group (PVD and Control). Comparisons between groups were made using t-tests for continuous variables and Fisher's exact tests for categorical variables (Table 1). We were unable to control for menstrual cycle due to the organization of scheduling within the MRI facility; however, a Fisher's exact tests of estimated menstrual cycle phase (i.e., follicular, ovulatory, luteal, menstrual, or no longer menstruating) based on the start date of participants' most recent menstrual period revealed no significant differences between women with PVD and Control women, p = 1.00. The groups only significantly differed in their self-identified ethnicity, with the majority (13/16) of PVD participants identifying as Canadian compared to the Control participants who identified primarily as Canadian (8/16), European (4/16), and Asian (3/16). Due to time constraints, four participants in the PVD group completed the thermal pain testing training session and were familiar with the thermal pain paradigm but were unable to complete the practice runs in the mock scanner. We performed an SEM analysis of BS/SC data excluding these four participants, and there were no differences. Thus, we included all participants in the analyses.

Women with PVD reported a wide range of PVD symptom duration, from 1 to 28 years (M = 9.1, SD = 8.04; based on 10 of the 16 participants who responded to this question). The pain ratings obtained, and temperatures used, during the MRI sessions are summarized in Table 3. These values were averaged over the five trials with stimulation for each participant and then averaged within the study groups. A comparison between groups was done using a t-test. Women with PVD required lower temperatures to elicit moderate pain than Control women; however, this difference was not statistically significant.

Connectivity values across regions were compared between all four participant groups (Control Pain, Control No Pain, PVD Pain, and PVD No Pain) for each of the epochs that were analyzed (before, during, and after stimulation). Results show significant differences in weighting factors (β) between groups, and across time periods, for a number of connections. The significant differences between group and condition comparisons are summarized in Table 4, and the connectivity values for each target/source region pair, and group/condition, are listed. During the stimulation period, there were significant differences in network connectivity between the Control and PVD groups in the Pain condition, and connectivity values tended to be higher in the Control group. These include connections between the PAG and hypothalamus, from the PAG to the NGc, NTS, and DRt, and from the spinal cord to the thalamus. When comparing the No Pain conditions (Control vs. PVD) there were fewer differences between the groups, although differences in connectivity from the PAG to the NTS were found both during what would be the stimulation period, and the period after stimulation. Connectivity values during the No Pain condition tended to be higher in the PVD group.

In the Control group, comparing Pain and No Pain conditions, higher connection strengths were consistently observed in the Pain condition. This included in the period before stimulation in PAG-hypothalamus connectivity and during the stimulation period in connectivity between the PAG and hypothalamus, NGc to LC, hypothalamus to NGc, and PAG to both NGc and NTS regions. In contrast, comparing the Pain and No Pain conditions in the PVD group revealed differences only in the spinal cord to thalamus connectivity, during the stimulation period, and the connectivity strength was stronger in the No Pain condition.

The analysis of variances (ANOVAs) of connectivity strengths in the SC and BS networks in relation to the group (Control participants vs. women with PVD) and the condition (Pain vs. No Pain), as well as the group × condition interaction effects, are shown in Figure 4. The ANOVA results identified several connections for which connectivity strengths varied specifically with the group (Control or PVD), regardless of the stimulation condition (Pain or No Pain applied) in all three time periods that were analyzed (before the painful stimulus would be applied, during painful stimulation, and after the stimulus was applied). These involved connections to and from the hypothalamus, PAG, LC, NRM, NGc, NTS, DRt, and PBN regions. Additionally, during the baseline period before a stimulus would be applied, ANOVA results showed connections that varied specifically with the group (Control vs. PVD), namely from the LC to the hypothalamus, from the hypothalamus and NTS to the PAG, and from the spinal cord to the NGc and PAG, but no connections exhibited significant condition or interaction effects.

During the stimulation period, there was a main effect of study condition (Pain vs. No Pain) in connectivity from the LC and PBN to the spinal cord, and from the hypothalamus to the PAG. The connectivity strengths of several connections such as from the LC, NTS, and DRt to the spinal cord, from the hypothalamus to the PAG, from the PAG to the NTS, and from spinal cord to the NRM, NGc, and PAG also varied significantly between the Control and PVD groups. However, during this time period, connectivity strengths varied mostly with group × condition interaction effects, primarily in connections from the LC and PBN to the spinal cord, the hypothalamus to the PAG, and from the PAG and the NRM to the spinal cord.

When performing SEM, we chose regions of interest based on previous research of descending pain regulation for connectivity analyses. As a result, there could be regions relevant to pain processing that were not included in the network model. However, by including more regions to ensure all possible anatomical structures are accounted for, we would increase the number of comparisons performed and run the risk of including non-significant network connections that weaken the strength of β values for each network component. In addition, PVD is a heterogeneous patient group with a wide range of symptom severity, comorbidities, and psychosocial characteristics. It is likely that many of these factors could contribute to alterations in pain processing mechanisms, especially with the knowledge that emotional states play a significant role in pain perception, which were not controlled for in the present analyses (51). However, women with and without PVD did not significantly differ on trait anxiety or severity of depressive symptoms (Table 1), and they did not significantly differ on self-reported anxiety measured on a 0–10 scale before, during, and after the MRI acquisition (Table 2).

This study is the first to examine the BS/SC of women with PVD, and it advances our understanding of altered central processing in this condition. The results demonstrate alterations in connectivity between BS/SC regions that are relevant to pain modulation and suggest that descending pain regulatory systems play a role in the expression of PVD. In comparison to Control participants, women with PVD displayed fewer differences in connectivity between stimulation conditions, indicating an insufficient response of the pain modulation system in women with PVD. In addition, it appears that women with PVD exhibit alterations in pain processing even in the absence of provoked pain. Our results provide strong evidence that women with PVD consistently exhibit altered pain processing on a neurophysiological level, a pattern also demonstrated in an investigation of brain connectivity in women with and without PVD (37). Further research is needed to investigate whether constant engagement of descending pain regulation is a common feature of chronic pain conditions.