Immediately prior to imaging, participants underwent a 45-minute stimulus and paradigm training session in a “sham” MRI lab within the Queen's University MRI Facility. The purpose of training was to familiarize participants with the study paradigm, including pain rating scales, the noxious thermal stimulus and timing of stimulation, as well as to ease anxiety and practice laying still to reduce overall bulk motion while in the MRI system. This practice also enabled us to confirm that participants could identify a noxious stimulus as being painful and an innocuous stimulus as being non-painful, using the pain rating scales. Participants were trained to use validated 100-point numerical pain intensity and unpleasantness rating scales (NPS), with verbal descriptors at intervals of 10 (Figure 1A) (41–43). Participants were encouraged to rate in increments of 5 (i.e., using values such as 45, 50, 55, etc.) to simplify the rating process, and the researcher continuously checked each rating with the participant to ensure that they were accurately implementing the scales. They were informed that pain intensity describes more of the discriminative aspect of pain whereas unpleasantness describes the emotional/affective component of perceived pain. The ratio of each participants' pain intensity rating to the temperature needed to elicit that pain rating was used as a “normalized pain score”. A higher pain score may indicate that participants who experienced a particular pain rating at a lower temperature are more sensitive to noxious stimuli (perceiving it as being more painful) than those that reported the same pain rating but required a higher temperature to produce it. This method was used to standardize our pain measures as each participant was individually calibrated to an experimental stimulus temperature which would produce a target pain rating of approximately 50 on a 0–100 scale, as described below.

To elicit acute experimental pain, noxious thermal stimulation was applied with an MRI-compatible robotic contact-heat thermal stimulator (RTS-2), which raised and lowered an aluminum thermode to and from the participants' skin via a pneumatic piston. This device was custom-made within our lab. The stimulus was applied to the thenar eminence of the right hand, corresponding to the sixth cervical segment of the spinal cord. The timing and duration of heat-contacts, along with thermode temperature, were under precise control by custom-made software in MATLAB® (Mathworks Inc., Natick, MA). Each test consisted of ten 1.5 s heat contacts with onsets every 3 s, over the span of 30 s to elicit sustained behavioural and neural responses and to avoid habituation of receptors. Participants experienced a variety of temperatures presented in the same order (ranging from 45 °C to 52 °C) and, once familiarized with the stimulus, were individually calibrated to a temperature corresponding to a tolerable average pain rating of 50 intensity units (“Moderate Pain”, Figure 1A). Participants were kept blinded to this objective, as well as to the temperatures used during the tests, to avoid any potential response bias. The upper limit of 52 °C was set to avoid causing damage to the skin. Additionally, participants were instructed to remove their hand from the RTS-2 if their pain ratings ever exceeded 70 intensity units in order to avoid causing distress or very strong pain. Participants were also informed that there would be two conditions applied while they were in the MRI system, a painful condition at the calibrated temperature (“Noxious”), and a warm, non-painful condition at a standard 40°C for all participants, similar to body temperature (“Innocuous”). Once participants were calibrated, they moved on to the next stage of training.

A mock-up of the MRI scanner (sham-MRI) was used to train participants on the stimulation paradigm and timing that they would experience in the MRI, and to familiarize them with the confined environment. This process was also intended to reduce variations in the data that may be caused by anxiety and bulk motion across repeated fMRI acquisitions. Participants were positioned supine in the sham-MRI with a mirror over their eyes to view a rear projection screen showing the pain intensity and unpleasantness scales, and the RTS-2 under their right hand. A simulated version of the Pain fMRI protocol was carried out at the calibrated temperature. Beginning with a baseline period of 120 s of no contact, the stimulus then made 10 contacts to the skin over 30 s, followed by 120 s of baseline (270 s total, Figure 1B). This method of stimulation was chosen to enable comparisons with our previous studies, and because it produces temporal summation of pain while avoiding habituation of nociceptors in the skin, and thus evokes robust BOLD responses (15, 44–46). Participants were instructed to silently rate the intensity and unpleasantness of each contact as they felt them, and to remember only the highest ratings on both scales. The peak ratings of pain intensity and unpleasantness were recorded, and the calibration temperature was confirmed or adjusted based on these ratings. While the task of mentally rating the intensity and unpleasantness of each contact may have influenced the pain or sensations that were experienced, this procedure was common to runs with noxious and innocuous stimuli and is essential to enable comparisons of the participant's experience to the measures obtained with fMRI.

Functional MRI was carried out on a Siemens 3 tesla MRI system (Siemens Magnetom Prisma, Erlangen, Germany). Participants were positioned head-first and supine with foam supports under their knees and surrounding their head to minimize bulk motion during scanning. The peripheral pulse was recorded from all participants with a sensor attached to their left index finger, and they were provided with a squeeze-ball to signal the experimenter in the event of an emergency, or if they did not wish to continue the study. The RTS-2 was positioned at their side, under the palm of the right hand, and foam earplugs were provided to dampen the loud noise from repeated scans. A 64-channel head and neck coil array was used to obtain images of the brain and brainstem, and a mirror positioned above the participants' eyes allowed them to view a rear projection screen which displayed a presentation of prompts and the pain rating scales during each run which were synchronized with the timing of each scan.

Before the onset of an fMRI run, participants were informed that scanning would commence, and at 60 s after the start of the run they were told which condition to expect for that run (i.e., Pain or Innocuous), at 120 s the participant was informed of the impending stimuli and the pain rating scales were displayed during and after the thermal contacts, and finally, at the end of the run, the participant was instructed to remember their peak ratings and wait for the experimenter to ask for them over the two-way intercom (Figure 1B). After setup, participants were instructed to remain as still as possible and wait for the experimenter to announce the beginning of each scan.

Each imaging session began with the acquisition of “localizer” images in three planes to provide a reference for subsequent slice positioning. A sagittal, T1-weighted anatomical scan was acquired to aid in normalization of functional data using a 3D MPRAGE sequence with 1 × 1 × 1 mm³ resolution, a repetition time (TR) of 1,760 msec, echo time (TE) of 2.2 msec, inversion time of 900 msec, and flip angle of 8°. A gradient-echo imaging method with echo-planar spatial encoding (GE-EPI) was used with a flip angle of 90°, and BOLD contrast. The 3D volume spanned from the top of the first cervical vertebra to the apex of the skull and was imaged with a TE of 30 msec for optimal T₂*-weighted BOLD sensitivity in the brain. The TR was set at 3 s per volume, and 90 volumes were recorded to produce a time-series spanning 270 s (4.5 min). Data were acquired in 49 contiguous axial slices, 3 mm thick, with a 192 × 192 mm field of view, and a 64 × 64 matrix, resulting in 3 mm isotropic resolution, with an anterior/posterior phase-encoding direction.

Multiple runs of each condition (Pain and Innocuous) were acquired in a randomly interleaved order and participants were informed of which condition to expect after 60 s of scanning (60 s prior to stimulation). The stimulation paradigm followed the sham-MRI run exactly, beginning with 120 s of baseline scanning, followed by 30 s of 10 heat contacts, and then 120 s of baseline. Participants provided their peak pain intensity and unpleasantness ratings at the end of each run, and these ratings were recorded. In between each run, the MRI operator confirmed that the participant was comfortable and alert before continuing. We aimed to acquire 10 total runs for each participant, 5 Pain and 5 Innocuous, creating a complete dataset with 450 volumes/condition/participant, however 3 participants ended the scanning session with 9 total runs due to technical issues and time-constraints.

To confirm that noxious and innocuous stimulation produced distinct behavioural responses, intensity and unpleasantness ratings and normalized scores were investigated across study conditions using non-parametric Wilcoxon signed-rank tests at a significance threshold of p < 0.05. The relationship between questionnaire scores, pain ratings, and normalized pain scores were tested across all individuals using Spearman's rho correlations, with significance inferred at a threshold of p < 0.05. Although the sample size is small, we were interested in understanding whether individuals' personal characteristics impact their experience of acute experimental pain.

Functional MRI data were pre-processed using Statistical Parametric Mapping software (SPM-12, The Wellcome Centre for Human Neuroimaging, UCL Queen Square Institute of Neurology, London, UK) in MATLAB (version 2021A, MathWorks, Natick, MA, USA). Pre-processing steps included conversion to NIfTI format, co-alignment to correct for bulk motion, slice-timing correction, and spatial normalization to pre-defined anatomical templates from the Montreal Neurological Institute (MNI). The first three time points were excluded to avoid periods of variable T1-weighting. No participants were found to have excessive motion during any of the fMRI runs (>1 mm translation, or >2°of rotation). Images were re-sized to 2 mm cubic voxels prior to normalization for compatibility with the MNI template, and data were “cleaned” to reduce noise by fitting and subtracting signal variations corresponding to the motion parameters determined during co-alignment. Normalizing the data to a common reference also served the purpose of eliminating any effects of participants changing positions between runs, and we were therefore able to combine data across repeated fMRI runs in each person.

Subsequent data analyses focused on characterizing temporal BOLD responses and relationships between regions known or suspected to be involved in pain processing, emotion, and autonomic regulation. To achieve this, we modified a network model that we had previously used in our lab to suit the purposes of this analysis (45, 47). We aimed to identify the relationships between study conditions, individual pain ratings and scores, the timing of the stimulation paradigm (i.e., before and during stimulus application), and personal characteristics (questionnaire scores). For the purposes of prior studies, we had created a combined anatomical template and anatomical region map that spans the brain, brainstem, and spinal cord (47, 48). For this study, the relevant reference images consisted of the MNI152 template, included in SPM12, and anatomical regions maps from the CONN15e software (49). Brainstem regions not included in the CONN15e region map were supplemented based on examples and anatomical descriptions (50–55), and freely shared atlases as described by Pauli et al. (56), Keren et al. (57), and Harvard atlases (http://www.med.harvard.edu/AANLIB/).

Structural equation modeling (SEM) is a family of statistical techniques which are used to identify patterns of correlation/covariance among a set of BOLD responses within and across ROIs, and to explain as much variance as possible. This method requires a pre-defined model of directional anatomical connections across the brain and brainstem, based on known neuroanatomy between ROIs (Figure 2), but does not require assumptions about the timing of BOLD signal variations. The network model used in the following analyses includes: brain regions—pre-frontal cortex (PFC), primary somatosensory cortex (SI), secondary somatosensory cortex (SII), anterior cingulate cortex (ACC), posterior cingulate cortex (PCC), insular cortex (IC), thalamus (Thal), amygdala (Amg), hippocampus (Hipp), and nucleus accumbens (NAc); midbrain regions—hypothalamus (Hyp) and periaqueductal gray matter (PAG); pontine regions—locus coeruleus (LC), and parabrachial nucleus (PBN) (52). These areas were chosen to cover a comprehensive array of centres for somatosensation, pain processing and perception, cognitive and emotional processing, and autonomic homeostatic regulation (3, 4, 11, 52, 58–61). Some existing anatomical connections were pruned from this network model to limit the number of comparisons and highlight the most important regions involved in nociceptive regulation and modulation. The same network model was used for all subsequent analyses. Data were averaged across voxels within sub-regions to reduce the number of statistical comparisons to be made and to increase the signal-to-noise ratio over that of single-voxel analyses. Each ROI was functionally divided into seven sub-regions based on time-series characteristics using k-means clustering. Once defined, identical sub-regions were used across the group for both study conditions. This process limits potential bias when dividing each ROI into sub-regions as it assumes that each ROI can have more than one function (15, 62–65). This method has successfully identified robust networks of connectivity across the brain, brainstem, and spinal cord in our previously published work (15, 47, 48, 63–67). For the present study, this method provides the means to investigate coordination of activity across the network of regions, and also enables identification of the sub-regions that best fit the model network. These sub-regions are used as regions-of-interest for subsequent analyses of BOLD responses to the different study conditions.

SEM was carried out by means of a general linear model to calculate linear weighting factors (β) which indicate the contribution of each connection to the overall network model, using the time-series data across participants, separately for each condition. If region A receives input signaling from regions B and C, and the BOLD signal time-series responses in these regions are S_A, S_B, and S_C respectively, then: SA=βABSB+βACSC+eA; where e_A is the residual signal variation that is not explained by the fit (68). The weighting factors were calculated separately for each network component, consisting of a region receiving input (target) with multiple regions providing input (source). However, in order to allow for dynamic variations in connectivity between time periods, the SEM method was applied using data only from selected time periods spanning 45 s (15 volumes/run) before stimulation, and spanning the stimulation period, for each run in each participant. Data from separate runs were combined for each participant. The fitting procedure determines the b values that minimize the residual signal variance. Networks were investigated for every combination of anatomical sub-regions of each ROI to identify the sub-regions that resulted in the best fits to the measured BOLD responses.

The significance of connectivity values (β) was determined based on their average values across the group, and the estimated standard errors. Significance was inferred at a family-wise-error Bonferroni corrected p_FWE < 0.05 which accounted for the total number of network combinations that were tested across combinations of anatomical sub-regions. The network consists of 2–5 source regions for each target region, resulting in between 49 and 16,807 possible network combinations, depending on the target. This statistical approach been validated previously (68). With this process, connections with significant β values were identified and used for subsequent second-level analyses. Differences in connectivity across conditions were then tested using paired, two-tailed t-tests at a family wise error corrected p_FWE < 0.05, accounting for the 51 distinct connections in the network. Connectivity weighting factors (β) were also analyzed across individuals by computing correlations of β-values with individual pain ratings and scores. Ratings from the pain condition were used as a measure of a personal characteristic to test relationships with connectivity values. R² values from these correlations were converted to Z-scores using the Fisher Z-transform with the number of participants. The significance was estimated based on a normal distribution and was inferred at a family-wise error corrected p_FWE < 0.05, again accounting for the total number of possible network connections, depending on the target region.

In order to investigate details of BOLD responses in specific regions-of-interest, a Bayesian regression technique was used to characterize variations in BOLD signal intensity at each time point (i.e., volume) across participants in relation to pain ratings and stimulation temperatures. This analysis was used to identify consistent features of BOLD responses in specific regions which were dependent on individual pain responses. Bayesian regression was applied to the BOLD time-course responses for each sub-region, in each individual, for both conditions (Noxious and Innocuous).

Bayesian regression was applied to each point in the time series, using pain ratings and temperatures as independent variables. The pain ratings and temperatures were first centered so that the average values across all participants were equal to zero, and scaled so that the largest differences from the average were equal to one. BOLD responses in each run were also expressed as the percent signal change from the average for the run. The data were then fit to approximate the consistent BOLD responses (S_BOLD) at the average pain and temperature ratings (S₀), plus linear estimates of the BOLD variations with pain ratings (S_p) and temperature (S_t) (15): SBOLD=S0+painSp+temperatureSt. The fitting process therefore enables us to estimate BOLD response patterns (S₀) independent of individual differences in sensitivity or the stimulation temperature used, as well as to identify how the BOLD responses varied systematically across participants with different pain responses. The expected BOLD response for a region can thus be identified at the average stimulation temperature, as being S₀ + S_p at the highest pain rating, and S_0—S_p at the lowest pain rating. Comparisons of activity were done by linear regression between average time-courses across conditions, for one sub-region at a time, and ranking the slopes of these fits based on their difference from a slope equal to 1 which would indicate perfect correspondence. From this, we chose specific examples of the most similar and least similar activity across conditions, and differences in average time-course activity were determined using point-by-point t-tests, corrected for multiple comparisons at a p_FWE < 0.05.

Significant connectivity was found in the periods before and during stimulation across the network in both conditions (i.e., |β| > 0), however no significant differences in network connectivity between conditions were found at the group level.

Differences were found across conditions in terms of which connections were correlated with pain intensity and unpleasantness ratings in the periods before and during stimulation (Table 3). In the noxious condition, connections with β values that were correlated with pain intensity ratings included Amygdala → Thalamus and PFC → LC before stimulation, and SII → IC during stimulation. Connections with β values that were correlated with pain unpleasantness included NAc → Thalamus before, and PFC → LC during stimulation. In the innocuous condition, correlations between connectivity and pain intensity ratings were seen in connections between PBN → LC before, and NAc → Thalamus during stimulation. Only one connection between PAG → Hypothalamus was correlated with pain unpleasantness ratings during stimulation in the innocuous condition.

Numerous connections were identified as having strong relationships between the difference in connectivity values across conditions (Δβ) and normalized pain scores. These relationships identified effects of individual differences across pain measures and neural connectivity in the brain and brainstem. Examples above a threshold of R = 0.7 are shown in Figure 3 for connections between the PCC → Hippocampus, PFC → Amygdala, and ACC → IC in the periods before and during stimulation. These connections all show strong relationships between normalized pain scores and connectivity across the two conditions.

The results of Bayesian regression analyses demonstrated BOLD time-course responses at the average pain intensity rating and temperature for each sub-region and condition. Figure 4 shows BOLD responses from specific sub-regions identified as being the most similar across conditions (left panel), and the least similar (right panel). The strongest agreement across conditions (most similar) were seen primarily in sensory and association regions including SI, SII, ACC, PCC, IC, thalamus, and hippocampus, while the weakest agreements were seen primarily in cognitive, affective, and autonomic regions including the PFC, hypothalamus, amygdala, PBN, and LC. Time-courses were also tested point-by-point for significant differences across conditions, and these are indicated on each plot (*). Note that the BOLD responses and the regions do not depend in any way on the SEM analysis, and the Bayesian regression serves only to provide a linear fit estimate of the time-series values at the average pain rating and stimulus temperature for the group.