Noninvasive Vagus Nerve Stimulation Alters the Neural and Physiological Response to Noxious Thermal Challenge

The mechanisms by which noninvasive vagal nerve stimulation (nVNS) affect central and peripheral neural circuits that subserve pain and autonomic physiology are not clear, and thus remain an area of intense investigation. Effects of nVNS vs sham stimulation on subject responses to five noxious thermal stimuli (applied to left lower extremity), were measured in 30 healthy subjects (n=15 sham and n=15 nVNS), with fMRI and physiological galvanic skin response (GSR). With repeated noxious thermal stimuli a group × time analysis showed a significantly (p < .001) decreased response with nVNS in bilateral primary and secondary somatosensory cortices (SI and SII), left dorsoposterior insular cortex, bilateral paracentral lobule, bilateral medial dorsal thalamus, right anterior cingulate cortex, and right orbitofrontal cortex. A group × time × GSR analysis showed a significantly decreased response in nVNS group (p < .0005) in bilaterally in SI, lower and mid medullary brainstem, and inferior occipital cortex. Finally, nVNS treatment showed decreased activity in pronociceptive brainstem nuclei (e.g. the reticular nucleus and rostral ventromedial medulla) and key autonomic integration nuclei (e.g. the rostroventrolateral medulla, nucleus ambiguous, and dorsal motor nucleus of the vagus nerve). In aggregate, noninvasive vagal nerve stimulation reduced the physiological response to noxious thermal stimuli and impacted neural circuits important for pain processing and autonomic output.


Noninvasive vagus nerve stimulation
Afferent and efferent vagus nerve signaling are critical mediators of physiological homeostasis, modulating heart rate, gastrointestinal tract motility and secretion, pancreatic endocrine and exocrine secretion, hepatic glucose production, and other skeletal and visceral functions that together make the vagus nerve the principle nerve of the parasympathetic nervous system (1). Vagal fibers can be activated with exogenous electrical stimulation carried out with surgically implanted vagus nerve stimulation (sVNS) devices (implanted around the vagus nerve in the carotid sheath). Surgically implanted vagus nerve stimulation is approved by the United States Food and Drug Administration (FDA) for the treatment of epilepsy (2) and for treatment-resistant major depression (TRMD); (3)(4)(5). However, cervical sVNS can result in complications, including hoarseness, dyspnea, nausea, and postoperative pain (6,7).
Noninvasive techniques for VNS have beneficial effects in treating epilepsy, depression, and pain. Treatment includes the use of devices that activate the auricular branch (termed Arnold's nerve) of the vagus nerve (8)(9)(10) and the cervical vagus nerve (found within the carotid sheath) (11). Cervical transcutaneous noninvasive vagus nerve stimulation (nVNS) has shown promising therapeutic effects in the treatment of acute and chronic migraine headaches (12)(13)(14), and acute and chronic cluster headaches (15), and is now FDA-approved to treat both episodic cluster (14) and acute migraine headaches (7,16,17). Recent work has shown that, with finite element modeling of cervical nVNS, the electrical field significantly penetrates the human neck and is sufficient to activate the cervical vagus nerve (11). Thus, growing evidence supports the notion that transcutaneous cervical nVNS results in vagal activation that affects pain transmission and experience.

Pain autonomic responses and vagal nerve stimulation
Pain is a multimodal experience represented by a broad network of cortical and subcortical structures, including the primary (SI) and secondary somatosensory (SII) cortices, bilateral insular cortex (IC), anterior cingulate cortex (ACC), prefrontal cortex (PFC), thalamus, and brainstem nuclei (18,19). Noxious thermal (painful) stimulation activates a sympathetic response, as measured by an increase in galvanic skin response (GSR); (20)(21)(22), with a dose response relationship to increasing thermal stimulus magnitude (23). Prior work has identified pain-mediated increased activation of the IC, amygdala, ACC, and PFC that correlates with pain-evoked sympathetic activity (i.e. GSR), and together offer a baseline construct for the neural basis of this autonomic pain dimension (24)(25)(26)(27)(28). In the present study, we used functional magnetic resonance imaging (fMRI) and primary physiological outcomes (GSR) to test the hypothesis that nVNS may alter typical cortical and subcortical neural and physiological autonomic responses to aversive noxious thermal stimuli more than to sham treatment. Prior literature supports antinociceptive effects of vagal nerve stimulation in preclinical pain models (29)(30)(31)(32)(33)(34). The antinociceptive effects of VNS are postulated to depend on afferent signaling to the nucleus tractus solitarius (NTS), nucleus raphe magnus (NRM), and locus coeruleus (LC) (31). Based on this work, it has been proposed that vagal afferent inputs to NTS, NRM, and LC result in a summative signal (including activation of descending noradrenergic, serotonergic, and spinal opiodergic tracts) that inhibits dorsal horn neurons (33) (31) (34).

Participants
Thirty male and female subjects (age range,

Intervention
Subjects were randomized to receive either nVNS (n=15) or sham (n=15) treatment (Figure 1a). A pair of nonferromagnetic stainless-steel surface electrodes (1-cm diameter) were placed on the subject and secured with an adjustable Velcro strap collar.
The 2 devices were identical in appearance and subjects were blinded to specific intervention. Application of the device was made to either the right anterior cervical area (overlying the carotid artery) for active nVNS, or the right lateral cervical area (posterior to sternocleidomastoid) for the sham treatment. Surface electrodes were connected to the battery-powered stimulation unit by a 6-m shielded, grounded cable. Both the sham and nVNS devices delivered 1-ms duration bursts of 5 sinusoidal wave pulses at 5000 Hz with a repetition rate of 25 Hz, and a continuous train duration of 2 minutes. In both the nVNS and sham treatments, a computational fixed, initial 30-second ramp-up period was followed by 90 seconds of peak stimulation. In the nVNS treatment, the voltage was increased to 24 V, whereas in the sham stimulation it was increased to 9 V. Sham low-voltage stimulation applied to the neck far lateral to the sternocleidomastoid produces a slight tingling sensation (activating cutaneous afferents) but does not penetrate deep below the skin surface or result in muscle activation (11). Both nVNS and sham stimulation were carried out 9.5 minutes prior to the noxious thermal stimulus paradigm (Figure 1b).

Figure 1. Study Design (a).
Subjects were screened and randomized to either the sham treatment or nVNS group. Sham stimulation was carried out posteriolateral to the sternocleidomastoid. In the nVNS group, stimulation occurred anteromedial to the sternocleidomastoid and lateral to the trachea. In both the nVNS and sham treatments, a computational fixed, initial 30-second ramp-up period was followed by 90 seconds of peak stimulation. In the nVNS treatment, the voltage was increased to 24 V, whereas in the sham stimulation it was increased to 9 V. Experimental design (b). Subjects were allowed to rest for 5 minutes before undergoing 2 minutes of nVNS (electrodes placed over carotid) or sham stimulation (electrodes placed far lateral to the sternocleidomastoid). Subjects then rested for an additional 5 minutes. Nine and a half minutes after either nVNS or sham stimulation, 5 successive noxious thermal stimuli were applied in bouts of 5 seconds each, up to 49.8°C. Each heat stimulus began 110 seconds after the start of the previous one. Measurements were taken 9.5 to 16.8 minutes after nVNS or sham.

Thermal stimulus task
The thermal heat threshold and thermal heat tolerance were obtained prior to the MRI scan, as previously described (46)  in all participants. Five noxious thermal stimuli were successively applied for 5 seconds each, with a 105-second interval between each application. The total duration of the task was 9 minutes and 15 seconds (Figure 1b). Ten seconds after each thermal stimulus ended, each subject was asked to rate their pain intensity on the numerical pain rating scale (NPRS). In response to each thermal heat stimulus, subjects indicated the appropriate pain-intensity score with a cursor pointing to the NPRS number 0 to 10 (where 0 = no pain, and 10 = most intense pain possible). The NPRS is a validated pain-intensity score, with a test-retest reliability of 0.71 to 0.99 that is highly correlated with the numerical pain rating scale and McGill Pain Questionnaire (47).

Galvanic skin response
We used the BioPac MP150 Psychophysiological Monitoring System (BioPac System Inc., Santa Barbara, CA) to measure psychophysiological reactivity at rest and during the noxious thermal stimulus pain paradigm. The GSR was recorded using 2 electrodes positioned on the volar pads of the distal phalanx of the middle and ring fingers of the right hand, and was sampled with a frequency of 1000 Hz. The mean GSR (in microsiemens) prior to the application of each (#1-#5) noxious thermal heat stimulus (baseline GSR) was compared to the peak GSR response after the application of noxious thermal stimulus for each trial (#1-#5). The slope of GSR from baseline to peak was calculated (microsiemens/s). Additionally, the time (in seconds) from baseline (prior to each noxious thermal stimulus) to the peak GSR response (each post-noxious thermal stimulus) was measured and compared within and between groups. The mean GSR response was defined as the average GSR (over 25 seconds) obtained after the peak GSR was reached. Data analysis, including sample selection and artifact removal, was carried out with AcqKnowledge software (version 4.42, BioPac System Inc.) and the R statistical programming language, version 3.4.3 (48). TR=8.1 s, TE=3.17 ms, 1 × 1-mm in-plane resolution) were acquired to permit activation localization and spatial normalization.

Group demographics of GSR analyses
Group differences in questionnaires and demographic analyses were calculated with Mann-Whitney U tests. BIOPAC system measurements of GSR were incorporated into a mixed-model regression to evaluate within-and between-group (nVNS vs sham) changes in GSR with each noxious thermal stimulus (from baseline, i.e. prior to each (#1-#5) noxious thermal stimulus to after the noxious thermal stimulus has been applied (#1-#5).
The within-and between-group GSR post-thermal noxious stimulus mean value (microsiemens), time to peak (seconds), and slope from the baseline GSR to the peak (microsiemens/seconds) were compared. All statistical calculations were performed using the R statistical programming language, version 3.4.3 (48).

MRI preprocessing
Structural and functional image processing and analysis were completed using analysis of functional neuroimages (AFNI) software (49) and R statistical packages. Echo planar images were slice-time and motion-corrected and aligned to high-resolution anatomic images in AFNI. Volumes with >20% voxels marked as outliers using 3dToutcount were censored and dropped from the analysis. For all group data points in the LME analyses 1.5 % data censor were identified as outlier. Percentage Outlier voxels in the time series were interpolated using 3dDespike. Functional data were aligned to standard space, resampled to 4-mm isotropic voxels, and smoothed with a Gaussian spatial filter (to 6 mm full width at half-maximum). Hemodynamics of the pain experience were modeled using line interpolation (3dDeconvolve/3dREMLfit modeled with TENT) for the span from the initiation of thermal heat stimulus and the following 15 seconds as modeled by 5 regressors overtime. These regressors were reconstructed to form a time series with 11 data points 1.5 seconds apart, which was used in subsequent analysis.

Participant demographics and psychiatric assessments
The mean age between the nVNS (24.7 ± 3.7 years) and sham group (30.7 ± 10.3 year) was not statistically different, as determined by a Mann-Whitney U test (p = .349).
Subjects did not report having elevated anxiety, depression, or posttraumatic stress disorder (PTSD), as measured by the Beck Anxiety Index (BAI), Beck Depression Inventory 2 (BDI-2), or the PTSD Check List-Civilian version (PCL-C). Accordingly, no significant difference in mean scores between groups was noted for these measures.
There were no significant differences in gender or race between the sham and nVNS groups. Two subjects failed the initial screen and were excluded from the study; one had a preexisting arrhythmia disorder (Wolf-Parkinson-White syndrome) and the other had braces ( Table I). The total sample used for analysis (after exclusion of the 2 subjects who failed screening) was 15 subjects in each of the VNS and sham groups.

Baseline pain measures
Subject responses to the baseline MPQ, measured at rest prior to thermal threshold or tolerance testing, were not different between the groups (Table I)

Pain reports during the fMRI task as measured by the NPRS
During the MRI task, 5 successive noxious thermal stimuli were administered based on thermal tolerance measures, up to 49.8°C (Figure 1b). The pain intensity score, measured as the mean NPRS score reported during the noxious thermal stimulus paradigm, was similar between the groups for each application of thermal stimulus  Figure 1). One subject in the sham group was unable to complete the fifth 5-second noxious thermal stimulus due to discomfort. No other adverse events occurred during the study.

Galvanic skin response
The GSR was recorded with each noxious thermal stimulus. The time from the onset of the each of noxious thermal stimuli to the peak GSR was measured in seconds. Mixedmodel regression analyses conducted across all noxious thermal stimuli (T1-5) and between groups (nVNS vs sham) showed a significantly shorter time to peak in the nVNS group (p = .020; Figure 2a). Post-hoc comparisons between groups (with a 2-sample t test) revealed that subjects who underwent nVNS had a shorter time to peak GSR compared with sham subjects during the application of noxious thermal stimuli T1 and T2 (p < .05). Similar trends also approached significance for T3 and T4 (p < .09; Figure   2a; Supplementary Table I). We then measured the GSR slope (in microsiemens) from the baseline GSR (prior to the application of each noxious thermal stimulus) to the peak GSR (accompanying each noxious thermal stimulus) and compared how this slope changed with each of the noxious thermal stimuli (T1-5). This GSR slope decreased equally in both groups for T1 to T3 (Figure 2b). But in contrast to the nVNS group, which had an average decrease in slope (-0.0461 microsiemens/second) for T3 to T5, the sham group showed an increase in the average slope to peak GSR from T3 to T5 (0.049 microsiemens/seconds), with a significant between-group difference observed (group x time interaction, -0.09508; p = .0412; Figure 2b).

Figure 2. nVNS vs Sham Autonomic Measures of Sympathetic Tone Galvanic Skin
Response (GSR) with noxious thermal challenge. (A) The time to peak galvanic skin response (GSR) measured in seconds after the application of each of the noxious thermal stimuli was significantly reduced in the nVNS group for noxious thermal stimuli 1 and 2 (T1 and T2) (**p < .05) compared with the sham group, and approached significance for T3 and T4 ( δ p < .09). Mixed-model regression showed that the combined (T1-T5) time to peak GSR in the nVNS group was significantly shorter compared with the sham group (p < .02). (B) The GSR slope (in microsiemens) from the baseline GSR (prior to the application of each noxious thermal stimulus) to the peak GSR (accompanying each noxious thermal stimulus) was measured in each group. The slope from the baseline GSR to the peak response decreased in both groups with each successively applied noxious thermal stimulus from T1 to T3. However, whereas the nVNS group showed a negative average slope to peak GSR of -0.0461 from T3 to T5, the sham group showed a positive average slope to peak GSR of 0.049 from T3 to T5. The between-group difference (group x time interaction = -0.09508) for T3 to T5 was significant at *p < .05.
Within-group analysis conducted using a Mann-Whitney U test showed that the mean GSR (measured for each of the successive noxious thermal stimuli) was successively lower in the sham group after the application of the noxious thermal stimulus

Group differences during the application of thermal stimuli
There were no between-group differences in BOLD responses 1.5 seconds before the application of noxious thermal stimuli. During the application of a noxious thermal stimulus, 21 regions met cluster thresholds in group × time LME analyses (i.e., nVNS vs   followed by an increase in OFC response (greater than sham) that was most evident at the 10 to 15 second mark.

Imaging results with LME analysis
To better understand the relationships between neural and autonomic measures during thermal stimuli, the GSR mean, measured from the peak after thermal stimulus for 15 seconds, was incorporated into a group (nVNS vs sham) × linear time x GSR LME analysis using AFNI's 3dLME to compare time-course data from the nVNS and sham groups. The group × time × GSR interaction showed that 3 regions met cluster thresholds; the postcentral gyrus/somatosensory cortex (Figure 4a, b), cerebellum/medullary brainstem (Figure 4c, d), and left occipital gyrus (Table IV). At the medullary level (i.e., level of the olive from the lower pons, spanning to the lower medulla) multiple afferent fibers enter the brainstem, including vagus, glossopharyngeal, hypoglossal, and accessory nerves that synapse on multiple brainstem nuclei (i.e., nucleus ambiguous (NAmb), dorsal motor nucleus of the vagus nerve (DMNX) and nucleus tractus solitarius (NTS)). Other brainstem nuclei important for pain processing (i.e., the rostral ventromedial medulla (RVM), rostral ventrolateral medulla (RVLM), and nucleus reticularis (Rt)) are also found at this level. Brainstem nuclei localizations were compared with graphical representations of brainstem nuclei from the Duvernoy atlas (50) and compared with prior grey and white matter brainstem maps by Biessner and colleagues (51). Subjects in the sham and nVNS groups were separated by median into high and low mean GSR categories, and the group × time × GSR interaction in the areas corresponding to the above nuclei (within medulla/brainstem) were examined. During the application of noxious thermal stimuli, subjects who underwent sham treatment and showed a high GSR demonstrated greater activity in the medulla/brainstem, compared with other groups (Figure 4c, d).

Discussion
The clinical studies in patients with brain lesions, and from brain-imaging studies of noxious painful stimuli have called this theory into question (52). Well-established brain regions that show bilateral activation upon the application of painful stimuli include the ACC, PFC, SII, insula, thalamus, and inferior parietal lobe (53)(54)(55)(56)(57)(58); and, in some instances, SI (53,(59)(60)(61). It is likely that nVNS-mediated bilateral decreases in SI represent modulation of cortical context and or anticipatory neurocircuits. We postulate that the observed effects of nVNS on bilateral pain-processing pathways may represent bilateral nVNS afferent signaling effects; possible afferent to bilateral efferent effects on the thermal (and possibly nociceptive) signaling pathways of the spinal cord; or direct disruption of normal bilateral thermal and nociceptive afferent neural firing patterns that either independently or collectively change the temporal dynamics of pain processing.

Effects on Left Dorsoposterior Insula:
In addition to nVNS-mediated bilateral SI and SII responses, our analysis showed

Effects on Bilateral Mediodorsal Thalamus and Anterior Cingulate Cortex (Area 24) (i.e. Medial Pain Pathway):
Beside lateral thalamic nuclei projections (i.e., ventroposterior-lateral and ventroposterior-medial thalamic nuclei) to the SI and SII, known to relate the sensorydiscriminative aspects of pain spinal pathways to limbic structures, the medial thalamic nuclei provide inputs to emotion-related brain areas, including the insula, ACC, amygdala, . In an illustrative case study, a patient with a somatosensory cortex stroke that spared the dorsal ACC (area 24) and thalamus (including mediodorsal thalamus) reported usual contralateral limb analgesia to painful stimuli, but the patient continued to reported an "unpleasant" feeling with the application of painful stimulus, suggesting in vivo separation of the affective and sensory discriminative pain pathways (87). We observed mediodorsal thalamus and dorsal anterior cingulate deactivation in the nVNS group, which likely indicates a key mechanism of the effect of nVNS on the medial affective pain pathway, in agreement previous studies (84,91,92). Based on this remarkable (but preliminary) finding in a future study we will measure nVNS effects on affective pain (i.e. pain unpleasantness).

Effects on Right Orbitofrontal Cortex:
In addition to the medial dorsal thalamic connections to ACC, there are known medial thalamic projections to the PFC, ventromedial-prefrontal, and orbitofrontal (OFC) cortices (93,94). The group × time analysis in the current study showed decreases in the right OFC response, suggesting that nVNS mediates this signaling during nociceptive stimulation (Figure 3 k, l).
Prior clinical work has also demonstrated involvement of the prefrontal and frontal cortical regions in reflecting the emotional, cognitive, and interoceptive components of pain conditions, negative emotions, response conflicts, decision-making, and appraisal of unfavorable personal outcomes (95,96). Multiple painimaging studies have found that the frontal cortical regions are critical for controlling functional interactions among key brain loci that produce changes in the perceptual correlates of pain, independent of changes in nociceptive inputs (64,97,98). Manipulating the cognitive aspects of pain, such as reappraisal, control, and coping, produce neural changes in the brain thought to be important in top-down processing. The lateral OFC expresses a contextual modulation of response that is widely implicated in emotional regulation and decision-making behaviors (99,100), and it has been postulated that the valuation of pain is context-sensitive, as classified by the OFC (101). Activity in the ventromedial cortex and the OFC has repeatedly been shown to be modulated by acute (78,(102)(103)(104)(105) and chronic VNS (78,102,106). In this cohort, we showed initial decrease in OFC activation during nociceptive thermal stimulation followed by an increase in OFC response (greater than sham) that was most evident at the post-thermal stimulation 10 to 15 second mark. This interesting finding suggests a decrease in the OFC affective appraisal of pain (0-6 seconds) followed by a subsequent late hemodynamic response increase that may reflect a resultant increase in pain-coping behavior. The observed nVNS-mediated decrease in the response of the OFC during the application of maximal noxious thermal stimulation is consistent with the results of prior VNS-treatment imaging studies, and suggests that the effect of nVNS on the OFC likely plays a role in the processing of painful and aversive stimuli.

Non-invasive Vagus Nerve Stimulation; Combined Neural Effects and Physiological Measures (GSR):
The group × time × GSR analysis highlighted differential interactions among nVNS, GSR, and the temporal dynamic of pain responses in the cerebellum, medulla/brainstem nuclei, bilateral SI, and a right occipital gyrus cluster. In addition to cortical nodes, the mid and lower medullary brainstem have been shown to be important sites that demonstrate an interaction between sympathetic output and pain, with decreases in sympathetic output (as measured with cardiac vagal tone) shown to correlate with brain stem nuclei

Non-invasive Vagus Nerve Stimulation; Autonomic Measures and pain report
Time to peak GSR (i.e., time from GSR measured immediately prior to each 5 second noxious thermal stimuli to peak post-noxious thermal stimuli) in the nVNS group was more rapid than in the sham group, indicating changes in the temporal dynamics of pain processing and subsequent sympathetic output. The temporal dynamic of GSR during the application of a thermal stimulus is an important component of autonomic responsivity (20)(21)(22), and the subsequent emotional regulation of aversive stimuli (113,114). Loggia and colleagues demonstrated the existence of a dose-response relationship between the magnitude of a thermal stimulus and the time to peak GSR (23). Specifically, their study showed that the greater the impact of a stressor (increased thermal temperature), the greater the rise in GSR, thus resulting in a longer time to peak response.
In addition to the longer time to peak observed in the sham group, significant differences between the nVNS and sham groups in the slope of the GSR rise from baseline (prior to each noxious thermal stimuli to peak after noxious stimuli) for the latter thermal stimuli (T3-T5) were observed. In particular, the slope of the GSR response decreased across the length of the task in the nVNS group, whereas the slope of the response in the sham group increased. Taken together the longer time to peak and increase in GSR slope in the sham group compared to the nVNS group further suggest nVNS alters sympathetic output, possibly to due to the brainstem and cortical effects described.
Together with group differences in GSR (time to peak and slope) there was also a significant difference in the change in response between subjects who underwent nVNS

Non-Invasive Vagus Nerve Stimulation Potential Temporal Dependent Effects on Brain & Pain
Henry and colleagues (77) first argued (2002), that neural effects which occur during VNS are very different from those that occur after VNS, while others continue to confirm this phenomenon (42,43,115). In this study, we showed that subjects in the nVNS group had nVNS-mediated activity decreases in the dorsoposterior insula, low medullary brainstem, processing that may also be of clinical import to guide efficacious use of VNS for clinical comorbid pain and psychiatric syndromes.

Limitations
Our work has some important limitations. The study was carried out in healthy control subjects. Because, as a pilot study, we involved only 15 subjects per group, the small sample size may not adequately represent a larger population. Therefore, our results as described here should be considered preliminary. However, the positive findings observed in this small cohort of healthy control subjects were robust and significant, warranting further investigation of the effects of cervical transcutaneous nVNS on the brain in a larger cohort of healthy control subjects, and in subjects who may experience a greater magnitude of affective pain subtypes, that may include Posttraumatic Stress Disorder or Major Depressive Disorder. Our study found significant neural alterations in the temporal dynamics of noxious thermal-stimuli processing known to be important in affective pain processing, and group differences in changes in the subjective pain report across the thermal stimuli (T2-T4). But we did not detect a difference in subjective reports of pain for each thermal stimulus with a near maximal noxious thermal stimulus (Supplementary Figure 1). We chose a near maximal noxious thermal stimulus to ensure clear autonomic responses (GSR). Our own work (46) as well as that of other studies has described maximal noxious stimuli that result in maximal reports of pain and, therefore, blunting of group differences in mean reports of pain (i.e. a ceiling effect on pain report) (117-119) (120). This phenomenon also could have occurred in this study. Future studies that measure affective pain (such as pain unpleasantness and catastrophizing) using maximal and submaximal noxious thermal stimuli are now needed to further characterize the antinociceptive effects of nVNS, as measured by reports of pain. While we correct for motion artifact at the brainstem level with the (Group x time x GSR) interaction, this area can be artifact-prone due to motion.
Although others have shown a similar pain and autonomic tone interaction at the same medullary brainstem level (Sclocco and colleagues (107)) future study is planned in larger cohorts to confirm this interaction at this brainstem level.

Conclusion
We    to the peak response decreased in both groups with each successively applied noxious thermal stimulus from T1 to T3. However, whereas the nVNS group showed a negative average slope to peak GSR of -0.0461 from T3 to T5, the sham group showed a positive average slope to peak GSR of 0.049 from T3 to T5. The between-group difference (group x time interaction = -0.09508) for T3 to T5 was significant at *p < .05.

Inclusion and exclusion criteria
All subjects were instructed to refrain from taking any over-the-counter analgesics or anti-

Heat tolerance and threshold measurements
Thermal thresholds and heat tolerance were calculated using the previously described method of limits (Yarnitsky & Sprecher, 1994) by taking the average of 5 thermal stimuli successively applied with a fMRI-compatible MEDOC probe to the right lower extremity (anterior shin) at an increasing slope of 1°C/s, from 32°C up to a maximum of 50°C.

Correlations between autonomic tone and pain reports
Autonomic measures taken during the application of thermal stimuli (time to peak GSR, mean GSR) and pain score were correlated within each group to better understand the relationships among these factors. In the nVNS group, mixed-model regression for the change in mean GSR with covariates for thermal stimuli showed significant negative correlations between mean GSR and thermal stimuli T3 (-0.955 ± 0.271; t = -3.528; p < .001), T4 (-1.429 ± 0.85; t = -5.106; p < .001), and T5 (-1.593 ± 0.295; t = -5.394; p < .001), that support the difference from baseline GSR to the peak slope (ie, a decrease in slope from T3-T5 in the nVNS group only). Lastly, mixed-model regression results (for the nVNS group only) showed that the time to peak GSR covaried with thermal stimulus trial and pain score, with a significant effect of pain score detected (1.035 ± 0.449; t = 2.305; p = 0.025). reported by subjects using the numerical pain rating scale (NPRS) after each noxious thermal stimulus did not differ between the sham and nVNS groups. Both groups had lower NPRS scores at T5 compared with T1 (NPRS decreased by -0.678 ± 0.209; t = -3.241; p = .002). In contrast to findings for the nVNS group, subjects who underwent sham stimulation had a positive slope in NPRS scores across thermal stimuli (i.e. the change in NPRS score with successive noxious thermal stimuli T1-T5) for T2 to T4 that was significantly different (slope in the sham group, 0.150 ± 0.122; vs the slope in the nVNS group, -0.233 ± 0.122; p = .0301) and also approached significance from T1 to T4 (sham group, 0.010 ± 0.847; vs nVNS group, -0.203 ± 0.847; p =.0785). Red circles = nVNS group. Blue circles = sham group. *p < .05; δ p < .08.

Supplementary
Supplementary Table I. Between-group comparisons for the time to peak and absolute mean GSR. The nVNS group showed significant decreases in the time to peak GSR for T1 and T2.
Supplementary Table II. Within-group comparisons for the time to peak GSR and absolute mean GSR for the sham stimulation group. In the sham group, the time to peak GSR increased from T1 to T4 and T5. The mean GSR measured after the application of noxious thermal stimuli consistently increased from T2 to T5 and from T1 to T4.