Noninvasive vagus nerve stimulation alters neural response and physiological autonomic tone 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 the nVNS group (p < .0005) 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 treatmentresistant 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 FDAapproved 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]. Moreover cervical nVNS is known to result in characteristic evoked potentials when measured with EEG that match evoked potentials produced by implanted vagal nerve stimulators [18]. Collectively, transcutaneous cervical nVNS results in vagal activation that affects pain transmission and experience.

Pain autonomic responses and vagus 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 [19,20]. Noxious thermal (painful) stimulation activates a sympathetic response, as measured by an increase in galvanic skin response (GSR); [21-23], with a dose response relationship to increasing thermal stimulus magnitude [24]. 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 [25][26][27][28][29]. 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 [30][31][32][33][34][35]. 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) [32]. 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 [34] [32] [35]. Adding to preclinical work, multiple translational clinical studies also show similar antinociceptive effects of acute [10,[36][37][38][39] and chronic VNS [40].
Recent fMRI studies have revealed that nVNS affects brain areas important in pain processing (e.g. the medial thalamus, dorsal ACC, IC, and PFC; [41][42][43][44], thus highlighting a potential supraspinal vagal influence on pain perception. Only a single small pilot study (n = 20) has evaluated the neural effects of transcutaneous VNS using auricular "Arnold's nerve" stimulation on experimental pain [37]. The results did not show a difference between groups, but a post-hoc analysis of "responders", i.e. subjects (n = 12) with increased pain threshold post-nVNS, showed decreased activation during the application of pain stimuli in the left dorsoposterior insula, ACC, ventromedial PFC, caudate nucleus, and hypothalamus [37]. Notably, this study performed continuous transcutaneous auricular VNS during the noxious thermal challenge, possibly confounding the results as emerging literature shows pronociceptive effects during actual VNS, while the antinociceptive effects occur post-VNS [45,46]. Taken together, the evidence accumulated to date suggests that VNS alters clinical pain perception, but that VNS must be carefully timed to produce antinociceptive effects.

Study objectives
The objective of this study was to gain a richer understanding of post-nVNS effects on sensory discriminative neurocircuits, affective pain neurocircuits, and the peripheral autonomic response to noxious thermal stimuli. Our goal was to determine the extent of post-nVNS neural effects on pain-related brain activation and autonomic tone. Taken together, this knowledge could guide and improve the efficacious use of nVNS in pain-disease states.

Participants
Thirty male and female subjects (age range, 18-54) were recruited through the Altman Clinical and Translation Research Institute at the University of California, San Diego Health System. Screening, exclusion, and inclusion criteria are found in Supplementary Information (S1 File). All participants were right-handed and provided written, informed consent to participate in the study. The Institutional Review Board at the University of California, San Diego Health Systems approved this study (UCSD IRB project # 150202).

Intervention
Subjects were randomized to receive either nVNS (n = 15) or sham (n = 15) treatment ( Fig  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 4.5 V. With sham stimulation, subjects generally experience greater discomfort (because of stimulation of neck muscles), with maximal tolerable amplitude typically only 2-8 V[44, 47]. Therefore, we employed a maximum sham voltage of 4.5 V to the far lateral neck position [44,47]. Prior work suggests that use of the far lateral neck position with low voltage does not result in stimulation of the vagus nerve [44,47]. Mourdoukoutas and colleague's [11], recent work shows (with finite element modeling) that active 20 V nVNS positioned directly over the carotid artery results in electric field penetrance that activates the vagus nerve. Based on this modeling we chose the maximum setting of 24 V known to activate the vagus nerve in the cervical neck. Both nVNS and sham stimulation were carried out 9.5 minutes prior to the noxious thermal stimulus paradigm (Fig 1B).

Thermal stimulus task
The thermal heat threshold and thermal heat tolerance were obtained prior to the MRI scan, as previously described [48] (S1 File). During the MRI scan, noxious thermal stimulation up to a (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. (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. fMRI, and GSR acquisition were taken 9.5 to 16.8 minutes after nVNS or sham treatment.  temperature of 49.8˚C was applied for 5 seconds via a fMRI-compatible thermode (probe size  3 x 3 cm; TSA-II, NeuroSensory Analyzer, MEDOC Advanced Medical Systems, Rimat Yishai,  Israel) attached via Velcro strap, to the left lower extremity (left anteromedial lower leg, anterior to the medial gastrocnemius) 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 (Fig 1B). To tabulate the subjective pain report, we used the numerical pain rating scale (NPRS), 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 [49]. Ten seconds after each noxious thermal stimulus ended(noxious thermal stimulus application from 1.5-6.5 sec), subjects visualized a projector screen that displayed NPRS, at which point they were asked to rate their pain intensity) numbered 0 to 10 (where 0 = no pain, and 10 = most intense pain possible). On the visualized screen a cursor was slowly moved across the NPRS scale from 0-10(left to right over approximately a 10 second time interval). Prior to scanning, the subject was instructed to raise the right thumb when the cursor indicated the pain level experienced. A video camera visualized the subject's right thumb as it was raised (when the cursor passed under the specific number indicating their (i.e., #5/ 10) numerical pain rating). The subject's pain report (number when thumb was raised) was then recorded in the source document. Pain ratings were provided 10 seconds after termination of the pain stimulus and thus could be separated in the slow event related design.

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 GSR 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 [50].

Statistical analysis
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 betweengroup (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 [50].
MRI preprocessing. Structural and functional image processing and analysis were completed using analysis of functional neuroimages (AFNI) software [51] 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 halfmaximum). 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.
Group differences in the time course of Blood Oxygen Level-Dependent (BOLD) responses over the entire course of the pain experience were measured over the 5 noxious thermal applications. Time-course data were modeled using AFNI's 3dDeconvolve TENT function. The TENT function is a linear interpolation of the hemodynamic response function over time described as piecewise linear splines. A group (nVNS or sham) × time, and (nVNS or sham) × time × GSR linear mixed-effects analysis (LME) using AFNI's 3dLME was conducted to compare time-course data from nVNS vs sham. Effects of interest included (group × time) and (group × time × GSR) interactions, in which all were fixed effects without covariates. Group and GSR were handled as between subject factors and time was a within subject factor. Multivoxel multiple comparisons were performed by Monte Carlo simulations (using AFNI 3dClustSim modeled with 3-perameter modeling noise) to reduce the potential for false positive results. A per-voxel threshold of p < .001, a cluster-wise threshold of p < .001, and a minimum number of 14 voxels per cluster were used. The Montreal Neurological Institute (MNI) atlas was used to identify clusters. Brainstem nuclei localizations in the group × time × GSR LME were compared with graphical representations of brainstem nuclei from the Duvernoy atlas [52] and compared to prior grey and white matter brainstem maps by Beissner and colleagues [53].

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 1). 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.

Pain and physiologic measures
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 2). Heat thresholds, measured using the method of limits, were similar across groups (nVNS, 41.2˚-C ± 2.8˚C; vs sham, 41.9˚C ± 2.0˚C; p = .935), as was heat tolerance, also, measured using the method of limits,(nVNS, 49.0˚C ± 1.4˚C; vs sham, 48.71˚C ± 1.2˚C; p = 0.467; Table 2).

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 ( Fig 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 (S1 Fig). Both groups reported NPRS scores that were lower with the fifth thermal stimulus (decrease in NPRS, -0.678, ± 0.209; t = -3.241; p = .002) compared with the first stimulus. We then compared the change in mean pain report (NPRS) across each of the successive noxious thermal stimuli (T1-T5) between groups. In contrast to the nVNS group, subjects who underwent sham stimulation showed an increase in NPRS with each of the successive noxious thermal stimuli from the second to the fourth (T2-T4) (this change in pain score for each of the successive noxious thermal stimuli (T2-T4) was calculated as a slope, i.e. sham slope; 0.150 ± 0.122) vs the decrease in NPRS with each of the successive noxious thermal stimuli observed for nVNS (T2-4), (nVNS slope; -0.233 ± 0.122; p = .0301) (S1 Fig). 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.  nVNS alters neural and physiological responses to a noxious thermal challenge

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. Mixed-model 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; Fig  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; Fig 2A; S1 Table). 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 ( Fig 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; Fig 2B).

Imaging results
Group differences during the application of thermal stimuli. During the application of a noxious thermal stimulus, 21 regions met cluster thresholds in group × time LME analyses (i.e., nVNS vs sham × time). Examination of this interaction indicated that regions in the left insula, right cerebellum/declive, and right cuneus had large clusters of greater activation (sham > nVNS). Additional regions important in the processing of thermal stimuli included the left somatosensory cortex, bilateral mediodorsal thalamus, right dorsal anterior cingulate gyrus, left supramarginal gyrus, and right medial frontal gyrus (orbitofrontal cortex [OFC]; Table 3). A TENT function analysis showed significantly greater activation during the application of noxious thermal heat stimuli in the sham group in the SI (Fig 3A and Fig 3B), SII (Fig 3C-3E), left dorsoposterior insula (Fig 3F and Fig 3G), and bilateral mediodorsal thalamus, as well as in the dorsal anterior cingulate (area 24; Fig 3H and 3J), and right medial frontal gyrus (OFC; Fig 3K and Fig 3L). 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). Mixedmodel 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. Red circles = nVNS group. Blue circles = sham group. https://doi.org/10.1371/journal.pone.0201212.g002 nVNS alters neural and physiological responses to a noxious thermal challenge

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 ( Fig 4A and Fig 4B), cerebellum/medullary brainstem ( Fig 4C  and Fig 4D), and left occipital gyrus (Table 4). 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 [52] and compared with prior grey and white matter brainstem maps by Biessner and colleagues [53]. 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 (Fig 4C  and Fig 4D).

Discussion
The effects of VNS on the central and peripheral neural circuits involved in pain and autonomic physiology are not well elucidated. In this study nVNS treatment (when compared to nVNS alters neural and physiological responses to a noxious thermal challenge sham) resulted in reduced responses in highly relevant pain-processing nodes. There was a significant alteration of autonomic tone, as determined by a decrease in sympathetic activity (measured with GSR) and attenuated activity in brainstem nuclei known to contribute to pain-mediated autonomic responses. These results provide preliminary evidence of significant nVNS modulation of central and peripheral autonomic neural circuits relevant to pain perception. nVNS alters neural and physiological responses to a noxious thermal challenge

Post-Non-invasive vagus nerve stimulation with noxious thermal stimuli; neural effects on bilateral somatosensory cortex 1 (SI) and somatosensory sortex 2(SII) (i.e., Lateral Pain Pathway)
Compared with subjects in the sham group, group × time LME analysis showed that subjects in the nVNS group had decreased neural activation of SI and SII, the medial dorsal thalamus, ACC, IC, and OFC-all brain regions associated with the processing of painful stimuli. Metaanalysis of human data from fMRI, EEG, magnetoencephalography (MEG), and positron emission tomography (PET) studies has shown that the commonest regions found to be active during an acute pain experience [19] are the SI and SII, thalamus, ACC, IC, and PFC, (comparable to areas that show decreased activity with nVNS in this study). Analysis of the group x time interaction showed a decrease in responses of the bilateral SI and SII somatosensory cortex, suggesting that nVNS mediates this signaling during the application of a thermal stimulus. These nVNS-mediated response changes in the SI somatosensory cortex strip match bilateral somatotopy to the lower leg, consistent with the placement of the Peltier heat probe. It is generally believed that somatosensory stimuli are processed primarily or preferentially by the hemisphere that is contralateral to the point of stimulation. However, evidence from clinical studies in patients with brain lesions, and from brain-imaging studies of noxious painful stimuli have called this theory into question [54]. 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 [55-60]; and, in some instances, SI [55, [61][62][63]. 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 painprocessing 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.

Post-Non-invasive vagus nerve stimulation with thermal stimuli; neural effects on left dorsoposterior insula
In addition to nVNS-mediated bilateral SI and SII responses, our analysis showed a unilateral decrease in the left dorsoposterior insula. The dorsoposterior insula exhibits an anterior-toposterior somatotopic organization in response to innocuous or noxious/painful stimuli as measured with fMRI [64][65][66][67][68]. Various painful stimuli, including hypertonic saline injection [65], thermal stimuli [66], and laser stimuli [69], have consistently reproduced this anteroposterior somatotopy within the dorsoposterior insula; specifically, rostral targets (head/neck) localizing more anteriorly whereas caudal targets (leg) localizing posteriorly [70]. Dorsoposterior insular stroke results in discrete thermoanesthesia and analgesia that equivalently mapped anteroposterior somatotopy, further supporting the idea that the dorsoposterior IC plays a critical role in the pain experience [71][72][73][74][75]. Neuroanatomical data have demonstrated that the lamina I spino-thalamo-cortical pathway convey both nociceptive and interoceptive information mapped to the viscerosensory cortex in the posterior and mid-insular cortex, which is then represented in the anterior insula [76][77][78]. Surgically implanted vagal nerve stimulators (FDA-approved for treatment of resistant depression and epilepsy) consistently [79][80][81][82][83][84] modulate insular cortex activity, thus pointing to the insula as a possible neuromodulatory target for nVNS. Moreover, while insular activity is known to increase during acute VNS [79][80][81], recent work has shown a resultant decrease in insular activity at 10 to 15 minutes post-nVNS [43]. In our cohort, there was a significant left dorsoposterior insula decrease in activity 10 to 17 minutes post-VNS that further support the temporally dependent dose-response effects of VNS. As a whole, the observed changes in the response to pain in the SI, SII, and left dorsoposterior insula with nVNS infer possible nVNS-mediated changes in neuronal firing patterns, either through direct brainstem effects, afferent cortical, or afferent cortical-to-efferent brainstem/ spinal cord effects on nociceptive signaling.

Post-Non-invasive vagus nerve stimulation with thermal stimuli; neural 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 ventroposteriormedial thalamic nuclei) to the SI and SII, known to relate the sensory-discriminative aspects of pain spinal pathways to limbic structures, the medial thalamic nuclei provide inputs to emotion-related brain areas, including the insula, ACC, amygdala, PFC, and other regions important in processing the affective-motivational dimension of the unpleasant pain experience [85]. In our study, the nVNS group showed a decreased response in the bilateral mediodorsal thalamus and dorsal ACC (Brodmann area 24) during the application of thermal stimulation, (with the group x time interaction). Prior clinical work shows that the mediodorsal thalamus is important in antinociceptive regulation [86], the processing of emotions [87], affective pain processing (pain unpleasantness) [66,88,89] [90][91][92], [87], thought to occur through mediodorsal thalamic connections with dorsal ACC (area 24). 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 report an "unpleasant" feeling with the application of painful stimulus, suggesting in vivo separation of the affective and sensory discriminative pain pathways [89]. 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 with previous studies [86,93,94]. Based on this remarkable (but preliminary) finding in a future study we will measure nVNS effects on affective pain (i.e. pain unpleasantness).

Post-Non-invasive vagus nerve stimulation with thermal stimuli; neural 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 [95,96]. The group × time analysis in the current study showed decreases in the right OFC response, suggesting that nVNS mediates this signaling during nociceptive stimulation (Fig 3K  and Fig 3L). 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 [97,98]. Multiple pain-imaging 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 [66,99,100]. 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 [101,102], and it has been postulated that the valuation of pain is context-sensitive, as classified by the OFC [103]. Activity in the ventromedial cortex and the OFC has repeatedly been shown to be modulated by acute [80,[104][105][106][107] and chronic VNS [80,104,108]. 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 including: 1) RVLM, 2) Rt, 3) NAmb, 4) DMNX and 5) the RVM, (all found superior to the obex at the level of the olive spanning to lower medulla) [109]. In this study, medulla/brainstem clusters from sham and nVNS groups were separated into high and low mean GSRs. Only the sham treatment group showed a high GSR, demonstrated by greater activity in the medulla/brainstem, compared with other groups (sham low, nVNS high, nVNS low). At the level of the medulla, where this interaction is found (i.e., superior to the obex at the level of the olive from the lower pons, spanning to the lower medulla) multiple afferent fibers enter the brainstem, including the vagus nerve, and the glossopharyngeal, hypoglossal, and accessory nerves, as well as multiple nuclei and tracts (i.e., DMNX, NTS, NAmb, RVM, RVLM, and the Rt). In particular, the Rt is proposed to be primarily a pronociceptive center that integrates multiple excitatory and inhibitory functions important in nociceptive processing [110]. The premotor nuclei (i.e., NAmB and DMNX) are critical in autonomic response patterns evoked by physiological and sensory stimuli [111] that culminate in efferent parasympathetic outflow and play a crucial role in parasympathetic reflexes, accepting input from the NTS that is the principal nucleus for incoming afferent signals from the vagus nerve [112]. The RVM is intricately involved in areas of endogenous pain modulation in the brain, conveying descending pain modulatory influences from the PAG to neurons located in the dorsal horn of the spinal cord. The ON and OFF cells of the RVM increase or decrease activity during the application of painful stimuli, respectively [113], with notable effects on descending pain-inhibitory circuits [114]. In sum, decreased activity in the medulla found in this study can be seen in reduced autonomic tone (reduced GSR in the lower GSR sham group), or through vagal nerve stimulation (in both nVNS groups, regardless of GSR (high or low)) suggesting that the default regulation of GSR can be decoupled through nVNS. We postulate that, even with increased GSR output (high GSR in the nVNS group), nVNS inhibits the response of the central nervous system to pain (in part) by blunting the response in key nuclei in the medulla that relay autonomic responses. Support for the relationship between the nVNS neural response and physiological response stems from altered autonomic sympathetic output (i.e., time to peak GSR and decrease in GSR slope). Future study is planned to examine this interaction in disease states such as Posttraumatic Stress Disorder or Major Depressive Disorder where dysfunctional emotional regulation and dysregulated autonomic output coincide.

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 [21][22][23], and the subsequent emotional regulation of aversive stimuli [115,116]. Loggia and colleagues demonstrated the existence of a dose-response relationship between the magnitude of a thermal stimulus and the time to peak GSR [24]. 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.
We did not detect a statistically significant difference in subjective reports of pain for each thermal stimulus with a near maximal noxious thermal stimulus (S1 Fig). However there was a significant difference in the change in response between subjects who underwent nVNS vs sham stimulation across thermal stimuli (T2-T4) in reports of pain, as measured by the NPRS. The group that underwent sham stimulation showed a progressive increase in NPRS (across thermal stimuli T2-4), whereas the nVNS group demonstrated a significant decrease in NPRS (across thermal stimuli T2-4). To further characterize the effects of nVNS, additional work is needed that carefully measures affective pain, such as unpleasantness and catastrophizing, associated with the application of noxious thermal stimuli.
Non-invasive vagus nerve stimulation potential temporal dependent effects on brain & pain. Henry and colleagues [79] 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 [43,44,117]. In this study, we showed that subjects in the nVNS group had nVNS-mediated activity decreases in the dorsoposterior insula, low medullary brainstem, medial thalamic, and ACC compared with subjects in the sham group (occurring 10 to 17 minutes after nVNS treatment). Similar to our post-nVNS effects observed on the low medullary brainstem, Frangos and colleagues also show post-cervical transcutaneous VNS effects in this time frame, (13-15 minutes after stimulation), in posterior insula, lower medullary brainstem and medial thalamic/ACC deactivation at rest [44], that provide a convergence of preliminary evidence supporting a temporal nVNS dose-response curve [44]. In line with the aforementioned post-VNS neural effects, emerging clinical literature also demonstrate post-VNS antinociceptive effects [40,118] while pronociceptive effects during VNS have also been reported [45,46]. Both prior literature and this study suggest that the temporally dependent neural effects (i.e. during vs post-stimulation) of VNS may be critical to clinically relevant pro-or anti-nociceptive effects of VNS treatment, and therefore should be taken into account in future clinical study designs. Moreover, future studies are planned to determine the temporal dose response curve on affective pain 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 statistically significant difference in subjective reports of pain for each thermal stimulus with a near maximal noxious thermal stimulus (S1 Fig). We chose a near maximal noxious thermal stimulus to ensure clear autonomic responses (GSR). Our own work [48] 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) [119][120][121] [122]. 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 and decreased spatial resolution. Although others have shown a similar pain and autonomic tone interaction at the same medullary brainstem level (Sclocco and colleagues [109]) future study is planned in larger cohorts and with high spatial resolution multiband imaging to confirm this interaction at this brainstem level.

Conclusion
We examined the neural effects of nVNS during a noxious thermal stimulus challenge, in the context of autonomic responses. We demonstrated 3 major findings; first, nVNS activity not only reduces peak responses to thermal stimuli in the SI, SII, medial thalamus, dorsal anterior cingulate (area 24), dorsoposterior insula, and OFC, which are important nodes in sensory discriminative pain, affective emotional pain, and interoception pathways, but also changes temporal dynamic responses within these nodes. Second, nVNS alters autonomic responses to noxious thermal stimuli, as measured by GSR, and therefore affects critical autonomic pain networks. Third, even with a higher GSR response being provoked by the application of noxious thermal stimulus, nVNS decreased the central nervous system response by blunting the usual reactions in key nuclei in the medulla that relay autonomic responses. These significant findings may improve effectual nVNS that, if tuned with careful dose-response curves in mind, could translate into efficacious targeted effects on pain and autonomic neural circuits.
Supporting information S1 Fig. nVNS versus sham numerical pain rating with maximal noxious thermal stimuli. After either nVNS or sham stimulation, 5 successive noxious thermal stimuli were applied (up to 49.8˚C) for 5 seconds each (T1-T5). Mean pain, as 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. (DOCX) S1 Table. 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. There was no difference in GSR absolute mean change between groups for all time points examined. (DOCX) S2 Table. 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. (DOCX) S3 Table. Within-group comparisons for the time to peak GSR and absolute mean GSR for the nVNS group. In the nVNS group, the time to peak GSR did not change between each successively applied noxious thermal stimulus. The mean GSR measured after each noxious thermal stimulus did not increase from T4 to T5, or from T1 to T2. (DOCX) S1 File. Inclusion and exclusion criteria. Heat tolerance and threshold measurements. Correlations between autonomic tone and pain reports. This file contains information on inclusion exclusion criteria, a description of heat tolerance and threshold measurements and finally pertinent correlations between autonomic tone and pain reports. (DOCX) S1 Data Set. Demographic and pain data sets. This data set contains demographic information and pain rating data sets.