A central mechanism enhances pain perception of noxious thermal stimulus changes

Pain perception temporarily exaggerates abrupt thermal stimulus changes revealing a mechanism for nociceptive temporal contrast enhancement (TCE). Although the mechanism is unknown, a non-linear model with perceptual feedback accurately simulates the phenomenon. Here we test if a mechanism in the central nervous system underlies thermal TCE. Our model successfully predicted an optimal stimulus, incorporating a transient temperature offset (step-up/step-down), with maximal TCE, resulting in psychophysically verified large decrements in pain response (“offset-analgesia”; mean analgesia: 85%, n = 20 subjects). Next, this stimulus was delivered using two thermodes, one delivering the longer duration baseline temperature pulse and the other superimposing a short higher temperature pulse. The two stimuli were applied simultaneously either near or far on the same arm, or on opposite arms. Spatial separation across multiple peripheral receptive fields ensures the composite stimulus timecourse is first reconstituted in the central nervous system. Following ipsilateral stimulus cessation on the high temperature thermode, but before cessation of the low temperature stimulus properties of TCE were observed both for individual subjects and in group-mean responses. This demonstrates a central integration mechanism is sufficient to evoke painful thermal TCE, an essential step in transforming transient afferent nociceptive signals into a stable pain perception.

Nevertheless, there are reasons to suspect some other CNS mechanism underlies nociceptive TCE. Across modalities contrast enhancement constructs behaviorally relevant perceptual features late in sensory processing, allowing for signal convergence and modulation by feedback projections 16,17 , and, in the temporal domain, for longer information retention 6,18,19 . Additionally, asymmetric spatial interactions have been observed between time varying painful stimuli and simultaneous but remote noxious stimulation 7 . This at least reveals events in the CNS can interfere with pain dynamics because peripheral nociceptive neurons have receptive fields on the order of centimeters and are unaffected by remote stimulus changes 20,21 . Consequently, we hypothesize a central mechanism underlies nociceptive TCE.
Responses to step-wise noxious stimuli cannot be reduced to a linear combination, or superposition, of responses to constituent boxcar stimuli, so this TCE is a nonlinear response. Here we leverage this property to probe for a central mechanism using psychophysics, theoretical modeling, and painful thermal stimuli. An optimal step-wise low-high-low stimulus paradigm is identified by exploring a recurrent feedback model trained on a preliminary dataset 22 . The complex stimulus is then fractured in two: a short high intensity boxcar stimulus bounded by a longer low intensity boxcar stimulus. Delivered simultaneously, their composite reproduces the low-high-low stimulus profile. A nonlinear response to the low-high-low stimulus is confirmed by comparison with responses to individual boxcars. We then bypass the ability of primary afferent fibers to detect changes in the overall stimulus profile by simultaneously delivering each boxcar stimulus to remote locations using two thermodes, and test if the response matches a superposition of constituent stimulus responses or is nonlinear. Nonlinear integration across peripheral receptive fields is sufficient to demonstrate a central mechanism for nociceptive TCE.

Results
Model exploration and stimulus design. A nonlinear stimulus-response function with perceptual feedback was fit to evoked pain ratings averaged across 11 naïve subjects (Eq. 1, Fig. 1A). The model fit performs well in capturing pain response features, sharing a greater proportion of variance with VAS ratings than does the stimulus (pain vs. stimulus: r 2 = 0.656; pain vs. model fit: r 2 = 0.855; forced timeseries data is not independent or identically distributed, but Pearson correlations serve as a similarity measure). The fit also shows a prominent nonlinearity following stimulus offsets which resolves within 15 s, matching the pain response during the first low-high-low stimulus epoch, but not later epochs, consistent with published reports on offset analgesia 7 . Model fitting was performed blind to model exploration or subsequent testing.
Model exploration was used to design a boxcar stimulus and a low-high-low complex stimulus which would evoke prominent and distinct nonlinearities in pain response. To find a candidate boxcar stimulus we first examined step functions that crossed the pain threshold (θ; initial temperature, T i ; final temperature, T f ; T i < θ < T f ; Fig. 2A, top left). Peak response was predicted 14-15 s after stimulus onset regardless of T i and with magnitude proportional to T f . Time to peak and subsequent adaptation dynamics were independent of T f . For the complex stimulus we examined positive and negative steps taken from suprathreshold holding stimuli (θ < T i < T f and θ < T f < T i ). These also showed approximately temperature invariant dynamics. They were distinguished by a shorter 3-4 s time to peak (or trough), followed by prominent damped oscillations which stabilized within about 40 s after stimulus onset ( Fig. 2A, top right, bottom left).
In typical fashion 23 predicted ratios of pain intensity were proportional to ratios of stimulus intensity, but with a time and T i dependent proportionality constant ( Fig. 2A, insert). In other words, In the range of experimentally tractable stimulus intensities (T i ∈ [45 °C, 49 °C]), the dependence on T i is trivial ( Fig. 2A, insert, shaded gray area). Thus, as a practical consequence of this relationship, an approximately 25% decrease in temperature relative to threshold (~2.5 °C for 45-49 °C stimuli) was predicted to transiently abolish pain. These predictions were unambiguously confirmed by preliminary testing in one subject (Pearson correlation, VAS vs. stimulus: r 2 = 0.570; VAS vs. model: r 2 = 0.755; Fig. 2A).
We applied these lessons in the design of a new low-high-low complex stimulus (15 s-15 s-45 s at T 1 -T 2 -T 1 , T 1 > θ, T 2 = T 1 + 2 °C) and also designed corresponding decompositions of this stimulus into two constituent boxcars (15 s-15 s-45 s at either T 1 -T 1 -T 1 or baseline-T 2 -baseline, hereon 'T 1 stimulus' and 'T 2 stimulus' , respectively, baseline ≪θ; Fig. 3). 15 s stimulation periods ensure peak response is achieved, the final 45 s period resolves Figure 1. A nonlinear feedback model (red) captures mean response of 11 subjects (black) to a painful thermal stimulus sequence (gray). Stimuli are similar to those previously used to study "offset analgesia" and include a 1 °C × 5 s up/down step. Pearson correlations between pain vs. stimulus and pain vs. model prediction are shown.
Complex stimulus paradigm maximizes nonlinear response. Nonlinear dynamics were tested using boxcar and complex stimuli in a novel group of 20 subjects. Results are consistent with model predictions (Fig. 3, Table 1). For boxcar stimuli, mean T 1 stimulus response showed a time to peak of 15.3 s (bootstrapped CI .95 = [11.5, 18.3], n = 50000 samples; times taken from start of response which lags the stimulus), and the T 2 stimulus response also follow this trajectory only truncated by the shorter stimulus duration (time to peak pain: 12.7 s, CI .95 = [9.3, 13.2]). The complex stimulus is a linear combination of T 1 and T 2 stimuli, but the mean pain The model stepresponse was explored outside the training space, and three different types of step responses were identified. Responses to threshold crossing steps peak at ~15 s while steps from supra-threshold 'holding' stimuli (T i ) peak after ~7.5 s. Oscillations follow which stabilize after 45 s. Finally, positive and negative supra-threshold steps show symmetric responses. For supra-threshold steps, relative change in pain is proportional to relative change in stimulus with a time and T i dependent proportionality constant (insert; color coordinated slopes are plotted for two time points marked on supra-threshold step response plots; shadows illustrate T i ∈ [45 °C, 49 °C], T i = 47 °C featured). Response dynamics are approximately invariant with T i and representative step-response functions are shown. (B) Model predictions were qualitatively evaluated in a single subject using a variety of threshold crossing and supra-threshold steps. Results confirm presence of nonlinear damped oscillations in all three types of step-responses, confirms dependence of nonlinearity amplitude on step size, and that 2.5 °C steps transiently abolish pain (epochs 3-5, 7). response is not a linear combination of evoked responses to boxcar stimuli. Although T 2 stimulation occurs on top of a preexisting T 1 stimulation, the T 2 evoked response features a shorter time to peak pain (t 2 maximum: 5.4 s, CI .95 = [5.0, 8.2]) than occurs with the T 2 boxcar stimulus (bootstrap p < 0.001). The model predicts this kind of TCE for supra-threshold step responses (3.6 s for supra-threshold vs. 14.6 s for threshold crossing steps). Additionally, the complex stimulus intensity returns to T 1 after the 2 °C stimulus offset, at which point another prominent nonlinearity is visible in the form of an "offset analgesia" type transient downward deflection (t 3 minimum predicted at 4.5 s; observed at 9.8 s, CI .95 = [6.6, 10.8]; predicted amplitude: 0.96, observed amplitude: 0.85 ± 0.067, mean ± SE; VAS amplitude calculated at t 3 minimum with respect to mean of t 1 max and t 3 max). This demonstrates the intrinsic contrast enhancing nonlinearities predicted by the model can be captured in the group level pain responses.
Subject level responses to dual thermode stimuli are nonlinear. Simultaneous dual thermode stimulation tested if central mechanisms are sufficient to evoke nonlinear contrast enhancing (TCE) pain dynamics. Only simple boxcar stimuli were used, but stimulus epochs featured coincident delivery of the T 1 and T 2 stimuli at different sites (Fig. 4). By design T 1 stimulation began 15 s before -and ended 45 s after-T 2 stimulation. Thus, the superposition of these two stimuli reproduces the optimal timeseries of the single thermode complex stimulus (Figs 3 and 4). Critically, spatial separation between thermodes during dual thermode stimulation exceeded single peripheral receptive field sizes. In order to potentially disambiguate the level of the neuraxis at which the prospective CNS mechanism might act, stimuli were delivered at three distances: near and far on the same arm, or on opposite arms (Fig. 4).
In order to distinguish monotonic site specific sensitization and habituation from intrinsic nonlinear dynamics (TCE), responses were first evaluated according to a simple on-off-on criterion. At the single subject level, many responses reflect transient absence of pain during the final 45 s of stimulus epochs (e.g. Fig. 5). However, if a pain free period interrupts distinct periods of unambiguous T 1 evoked pain, then neither habituation (a continuous decay process) nor sensitization (a continuous rise process) offers adequate explanation. Instead, a nonlinear TCE response provides the most parsimonious explanation (e.g. Fig. 5, gray dashed boxes).
This rudimentary on-off-on criterion for nonlinearity demonstrates the dual thermode stimuli evoke nonlinear dynamics which can be discerned within individual subject responses. The complex single thermode stimulus acts as a positive control, and evokes significantly more frequent nonlinear responses (69% non-linear, n = 52) than the simple T 1 boxcar stimulus which in turn serves as a negative control for random VAS variance (4% non-linear, n = 55; χ 1 2 = 48.5, p < 1e-05; Cochran-Mantel-Haenszel Chi-squared test, accounting for stratification by stimulus intensity which was adjusted according to pain tolerance). The greatest incidence of nonlinearities among dual thermode stimuli occurs with delivery of both to the same forearm, either near each other or far apart (50% non-linear, n = 32,  on the opposite arm shows less frequent nonlinear responses, but still more than from the single thermode T 1 stimulus (19% nonlinear, χ 1 2 = 4.44, p = 0.035, n = 27; all tests taken with respect to the negative control and significant after Holm-Sidak correction for multiple comparisons; Table 2). Thus, dual thermode stimuli can produce nonlinear stimulus response dynamics, clearly discernable at the level of individual subjects, although there may be a distance or laterality dependence.

Group level dual thermode responses match the nonlinear complex stimulus response dynamics.
Dual thermode stimuli evoke nonlinear responses, but single thermode complex stimuli show rich and predictable dynamics not captured by the simple on-off-on criterion. For instance, the psychophysical model exploration and single thermode stimuli highlight characteristic maxima and minima of group mean TCE responses (Figs 2 and 3), while individual subject responses are discontinuous, resembling a random jump process with indistinct maxima and minima (e.g. Fig. 5), obscuring this precise temporal information. Employing group mean responses thus allows for further comparison between the precise dynamics of different stimulus responses.

Figure 5.
Stimulus and pain ratings from a representative subject exhibiting nonlinear pain for single complex stimuli and for dual stimuli configurations (near, far, and opposite arm). Each stimulation sequence consisted of two dual-thermode stimuli (red, blue, gray ratings), and three single thermode stimuli (black ratings). Dual thermode stimuli were delivered at three different distances from one another: near or far on the same arm or on opposite arms, and the order of these was randomized across subjects. The disproportionate diminution of perceived pain following a small offset shows that this subject rates all dual and single complex stimuli nonlinearly, in contrast to a linear sum of the simple stimuli. A more stringent but less ambiguous on-off-on criterion for nonlinear dynamics is also satisfied in many of these instances (gray boxes, one example and one counter example explicitly labeled): despite invariance of stimulus intensity there is pain before T 2 stimulation, no pain immediately after T 2 stimulation, and then pain comes back before the T 1 stimulus ends. Thermode roles switch between the first and second dual thermode stimulus epochs.

Complex
Near Using group mean responses, time of (a) t 1 peak pain, (b) t 3 minimum pain, and (c) t 3 peak pain were identified for the complex stimulus, and used as reference points for evaluating TCE in dual thermode stimulus responses (Fig. 6A). Differences in stimulus intensity between subjects did not affect reported pain or moderate TCE (p > 0.4 for T 1 and TCE*T 1 interactions, Table 3), so these terms were dropped from the model. However, the TCE reference points show significantly different pain ratings in near dual thermode stimuli, with pain at (a) and (c) significantly higher than (b) (planned comparison mean (a, c) -(b): 7.86 ± 2.75 VAS difference, mean ± SE., 47% pain reduction, p = 0.01). The same pattern occurs with far dual thermode stimuli responses (planned comparison: 7.15 ± 1.99 VAS difference, 50% pain reduction, p = 0.002), but not for opposite arm dual thermode stimuli (planned comparison: −0.39 ± 2.29 difference, 2% pain increase, p = 0.87; Table 4). Similarly, timing of t 3 extrema observed with near and far stimuli compare favorably with extrema of the single thermode complex stimulus (bootstrapped contrast of local minima timing, near vs. complex: p = 0.24, far vs. complex p = 0.18, opposite vs. complex: p < 0.001; Table 1). Therefore, during a stimulus of invariant intensity, nonlinear dynamics are evoked across near and far dual thermode stimuli at the same timepoints which characterize the single thermode complex stimulus response.  ) and (c) in dual thermode stimuli, but pain responses at (a) and (c) were significantly higher than (b) after near (p = 0.010) or far T 2 stimulation (p = 0.002; planned contrast). No differences were found after T 2 stimulation of the opposing arm (p = 0.87). (B) In near and far dual thermode stimuli the nonlinear response dynamics are the same as for the single thermode complex stimulus. Pain return dynamics were specifically examined between (b) and (c) (top panel, insert). Group averaged responses to near or far dual thermode stimuli significantly correlate with responses to the complex but not T 1 stimulus while responses to opposite arm stimulation correlates with neither (left). Correlations between responses to near, far, and single thermode complex stimulation support the conclusion that the same process underlies all three nonlinear responses (multiple regression with ARIMA errors; near vs. complex: t = 4.54, p < 1e-05; far vs complex t = 1.99, p = 0.047; df = 437). Dynamical equivalence is best illustrated with variance normalized (z-scored) timeseries (right). *p < 0.05. NS = non significant.
Scientific RepoRts | 7: 3894 | DOI:10.1038/s41598-017-04009-9 Although local maxima and minima reflect characteristic stimulus dynamics and edge enhancement, trajectories between extrema also characterize pain dynamics. However, most group mean stimulus responses were too abrupt to robustly model. Changes in pain in response to stimulus changes affect the statistical properties of   Table 4. Final mixed effects model shows a statistically significant nonlinear response for near and far dual thermode stimuli, but not the opposite arm stimulus.
the VAS timeseries (e.g. mean and variance) in a complicated way. These changes effectively introduce heterogeneity in the probability distribution of the timeseries from one stimulation period to the next (i.e. timeseries structure is "non-stationary"). This problem is compounded by the relatively short duration between stimulus changes (e.g. t 1 and t 2 pain ratings only span 214 frames each), which limits the number of effective degrees of freedom available to model these statistical properties. This confounds linear timeseries modelling of complete stimulus responses. Nevertheless, the pain return during t 3 is sustained during a period of constant stimulation, which simplifies modelling. Dual thermode responses between timepoints (b) and (c) were therefore modeled as a linear combination of two independent variables, each representing a competing hypothesis: (1) the single thermode complex stimulus response (nonlinear hypothesis) and (2) the single thermode T 1 stimulus response (linear hypothesis). Although competing, these hypotheses are not necessarily mutually exclusive, and combining both together in a single model recognizes the possibility that responses to dual thermode stimuli may be a hybrid of single thermode stimulus-response patterns. Regression models of near and far dual thermode responses support the nonlinear hypothesis (multiple regression with ARIMA(1, 1, 0) errors; near β complex = 0.167, t = 4.54, p < 1e-05; far β complex = 0.083, t = 2.03, p = 0.043; df = 437). The psychophysical model predicts the overall shape of the nonlinear trajectory should be the same regardless of how much pain is reported (i.e. dynamics show amplitude invariance) ( Fig. 2A), a property clearly illustrated when comparing amplitude (i.e. variance) normalized responses to the complex stimulus with responses to ipsilateral dual stimuli. Conversely, neither linear nor nonlinear hypothesis offered a consistent explanation when modeling responses to the opposite arm stimulus, suggesting it may have unique dynamical properties (Fig. 6B, Table 5).

Discussion
Temporal contrast enhancement in pain has thus far mainly been explored in the context of "offset analgesia" and a single stimulus pattern. Here we formally generalize this process to an entire class of stimuli with a nonlinear recurrent feedback model of temporal integration, and identify an invariance relationship underlying responses to step-wise dynamic stimuli. This enables simple and accurate construction of a highly effective novel stimulus paradigm for studying TCE in pain. Using simultaneous stimulation of different skin sites, we show the CNS is sufficient to perform nonlinear spatiotemporal integration of pain, and that resulting TCE pain dynamics are shared with nonlinear temporal integration of single stimuli.
Model extrapolations predict three distinct nonlinear stimulus-response motifs, two of which are contrast enhancing, and all are confirmed by single thermode stimuli. TCE responses to incremental stimulus increases and decreases are symmetric, peak around five seconds after stimulus onset and rapidly rebound. These dynamics enhance saliency of stimulus changes, and are absent in responses to initial stimulus onset. Critically, the stimulus features of the model training dataset are mainly 5 s long, too short to fully resolve stimulus-response motifs which show characteristic extrema and relaxation dynamics over 15-45 s of constant stimulation. Nevertheless, modeling successfully predicts the magnitude of downward deflecting nonlinearities, and peak and rebound dynamics of their upward deflecting counterparts. Although each of these motifs has been shown in some form before 8,9 , this is the first time to our knowledge that they are systematically characterized and unified under a single theoretical framework. This shows that TCE pain dynamics behave as if there is an underlying recurrent feedback process.
The model's ability to successfully extrapolate stimulus-responses outside of its training space lends credence to its use for mechanistic insight and as an interpretive lens, but limitations must also be recognized. Because our usage intent was primarily instrumental other models were not considered for comparison, but alternatives have been proposed which successfully predict response dynamics like adaptation and habituation 24 , processes our  model was not designed to accommodate. Perhaps as a result our model overestimates the rate of pain rebound for the complex stimulus. Additionally, the predicted response symmetry between stimulus increases and decreases contradicts recent findings using painful laser stimuli 9 , which suggest more abrupt contrast enhancement for stimulus increments. However, it should be noted our model breaks symmetry in a similar manner as stimuli deviate from idealized steps, e.g. by using finite stimulus ramp rates. Any mechanistic insights gained from this model must undoubtedly be interpreted in the context of supplemental mechanisms affecting stimulus response properties. In fact, continuous rating of subjective states could be influenced by many factors such as attention, expectations, prior experience, and even trust in the experimenter. We imposed minimal control on such factors. Most importantly, the instructions of what to rate are critical, and subtle changes in the latter may lead to different outcomes. Finally, our model parameters were derived based on data collected mainly in males, while dual thermode experiments were conducted primarily in females. Although the robust extrapolation of model outcomes to the dual thermode results suggests that sex may not be an important contributor to pain TCE mechanisms, this remains to be systematically studied, as sex is an important factor in pain psychophysics. The small dynamical range of our thermodes limited the number of subjects that we could study; it remains unclear whether this introduced bias in our model fit. Further studies are needed to clarify the limits of the model's extrapolations and dependence on stimulation apparatus and modality, task conditions and demographic composition. The role of peripheral nociceptive neurons in TCE has not been systematically explored. Different classes of first order afferents contribute to stimulus dynamics by conveying unique stimulus information depending on response thresholds, latencies and adaptation dynamics. For instance, low threshold nociceptive C-fibers begin responding at 39-41 °C 25 , are unmyelinated and relatively slow to transmit afferent signals 26 and rapidly adapt to stimulation 27 , properties conducive to detecting slow stimulus changes. Type-II Aδ fibers also rapidly adapt to stimulation but are much faster acting 28 providing a second complementary method through which to detect fast stimulus changes 29 , while type-I Aδ fibers 27 are slow adapting and might signal persistent stimulation. However, activity of first order afferents does not correlate with experimental tonic pain sensations 27 . Although we cannot a priori preclude that mechanisms in the periphery contribute to TCE dynamics, pain sensations appear more likely to be indirectly constructed from stimulation of the sensory epithelium, just as in vision, where a crisp and stable sensory representation is indirectly constructed from the radically different retinal image that is subject to constant saccadic jitter and optic distortions by ocular tissues and fluids.
Nevertheless, this study's primary finding shows central mechanisms are sufficient for temporal contrast enhancement in thermal pain perception. When complex stimuli are deconstructed into constituent boxcars and delivered to the same arm either near or far apart, the subject's response is not merely a superposition of responses to the constituent stimuli. Instead, a characteristic nonlinear, group mean response is observed with unambiguous subject level incidence. Receptive fields of first order nociceptive afferents are on the order of centimeters, while distances between thermodes exceed these and are up to an order of magnitude larger in the 'far' configuration. Therefore, subjects must integrate information from different peripheral afferent receptive fields to experience a TCE pain response. This entails nonlinear spatiotemporal integration of these stimuli at a higher level in the CNS.
Although we don't distinguish between pain ratings and pain experience, there are undoubtedly differences due to reporting errors and the interpretive decisions inherent in self report. As we illustrate with single subject data, self report is discontinuous and is frequently a succession of abruptly changing steps, while (subjectively speaking) experience is continuous. Much of this error variance cancels out when averaging across subjects. We mitigate the risk of an interpretive confound by using an optimal stimulus paradigm. Identifying a nonlinear response depends on the temporal contrast of experienced pain, and our optimized stimulus paradigm ensures this contrast is maximized and highly salient, as shown by single thermode stimuli (85% pain decrease) and ipsilateral dual stimuli (50% pain decrease). This improves on all existing studies of offset analgesia since gross changes in pain intensity are much easier to identify and report than subtle ones.
Our experimental design precludes peripheral temporal integration of dual thermode stimuli, but the results do not provide further insight regarding the level of the neuraxis at which integration occurs. It may occur anywhere there are neurons with receptive fields large enough to circumscribe both thermodes, but candidates already exist at the level of the spinal dorsal horn (e.g. wide dynamic range neurons 30 ). Thermodes were positioned on different sides of the radial midline of the arm, so it is likely different dermatomes were stimulated and spatial integration occurred across spinal segments. It was further hoped responses to opposite arm stimuli might identify a supraspinal mechanism with a bilateral receptive field. While the incidence of complete transient abolition of opposite arm pain remains statistically significant, the group dynamics during the prospective pain return period (t 3 ) do not match the boxcar stimulus (as required by the superposition principle), the complex stimulus, nor a linear combination of the two. An additional nonlinear spatiotemporal integration process may yet be involved. Therefore, our findings do not distinguish between prospective segmental, intersegmental or supraspinal mechanisms.
Inhibition of pain after painful stimulation of a remote site might be achieved by spatial filtering acting in parallel with temporal contrast enhancement. For instance, in CPM one painful stimulus suppresses responses to additional painful stimuli elsewhere on the body. Although CPM and TCE in pain are mediated by distinct mechanisms 15 , dual thermode stimuli are both spatially and temporally integrated. However, transient pain reduction during our dual thermode stimuli follows the same dynamics as a single thermode complex stimulus, supporting the conclusion that nonlinear responses to dual thermode stimuli are the result of a mechanism shared with nonlinear responses to single thermode stimuli rather than any additional mechanism of spatial filtering.
The construction of sensory representation through successive levels of spatiotemporal filtering and amplification is a fundamental principle underlying the organization of the nervous system. By matching dynamics elicited by spatially localized or dispersed time-varying noxious stimuli we demonstrate the existence of a CNS mechanism for temporal contrast enhancement of thermal painful stimulus changes in pain perception.

Methods
Subjects. Data were collected from 20 individuals recruited from the Northwestern University community.
All provided informed consent for all experiments, were right handed, at least 18 years old (age: 26 ± 3 years, mean ± SD; 6 males), and had no history of chronic pain, depression, dementia or other psychiatric conditions. A model training dataset was also obtained from 11 subjects recruited from the Chicago area (age: 23 ± 2 years, mean ± SD; 8 males). An additional 9 subjects were excluded from the training dataset due to equipment malfunction (3 subjects) or due to failure of a subset of training stimuli to evoke pain (6 subjects). All experiments were approved by and conducted in accordance with the ethical principles set forth by the Northwestern University institutional review board.
Equipment. Stimuli were delivered by two Medoc PATHWAY ATS thermodes (30 × 30mm) controlled by two separate PCs, but were calibrated using a common external reference thermometer to protect against unintended stimulus mismatches. Temperature and event related data was sampled at 200 Hz. Real-time pain ratings were obtained using a finger device with VAS feedback, and sampled at 14.3 Hz using in house LabView software. Visual feedback was rendered as a moveable yellow bar alongside a number scale (0-100) 31 . The LabVIEW software also controlled a digital input/output device (National Instruments NI USB-6501) which delivered a 0.4 Hz heartbeat (synchronization) signal to each Medoc device during trial runs. These reference signals synchronously triggered stimulus delivery and facilitated alignment of stimuli and rating data during preprocessing.
Experiments were conducted in a purpose built space consisting of two adjacent rooms to minimize contact between experimenter and subject. The experimenter, PCs and Medoc systems were in one, while the subject, stimulating thermodes, feedback monitor and finger span device were in the other. Subjects were provided with and instructed in the use of an emergency stop button which would terminate stimuli if they became intolerable.
Modeling. Pain was modeled using a nonlinear second order differential equation which has been previously described 22 .
The model was fit to a training dataset which was not used further in this study. Five stimuli were delivered consisting of 5s-5s-15s or a 5s-5s-30s stimulation, following a T 1 -T 2 -T 1 temperature profile, where T 1 ∈ [44 °C, 48 °C] and T 2 = T 1 + 1 °C. Stimulus timing profiles and 15 s or 45 s interstimulus intervals were pseudorandomly selected and preprogrammed into the stimulator, but delivered identically to all participants (Fig. 1). Mean stimulus-response data was used to obtain the model fit with a random walk over the parameter space. The obtained model fit was used throughout the rest of this study: (θ, α, β, γ, λ) = (37.1913, 2.4932, 36.7552, 0.0204, 0.0169).
Experimental Task. Subjects were instructed to rate their overall pain from thermal stimulation using the finger device. The experimental task at times involved rating pain during delivery of two simultaneous stimuli of different intensity, but our a priori hypothesis only requires examining the periods flanking this period. Rather than introducing superfluous task demands by asking the subject to attend to a specific thermode, instructions were deliberately vague with respect to the source of the stimulus. Thus subjects were allowed to freely choose how to reconcile the two potentially competing stimuli in their rating.
Before experimental trials, subject pain thresholds and tolerance were evaluated. Six 10 s thermal stimuli were delivered at differing intensities to the right volar forearm. Stimuli differed in 1 °C increments, and each was used twice in random order. The procedure was repeated until (1) all stimuli used were painful but tolerated, (2) subject rated stimuli of equal intensity the same (±10 VAS units). Between repetitions stimuli were adjusted as needed to accommodate tolerances and pain thresholds of subjects. The least tolerant individuals completed training using 44 °C, 45 °C and 46 °C stimuli, while those with highest pain threshold completed training with 47 °C, 48 °C and 49 °C stimuli.
Experimental stimulation involved 2 noxious stimulus intensities (T 1 and T 2 , where T 1 + 2 °C = T 2 ; baseline = T 1 − 16 °C), used in 3 stimulation sequences, each with six stimulation epochs (Fig. 5). Overall, these thermal pulses were either a complex, low-high-low sequence as used in "offset-analgesia" studies, or its component low and high pulses, presented either on the same skin site, or at different skin locations. Sequences had identical stimulus programming: 15 s of T 2 , a dual thermode stimulus (75 s T 1 , 15 s T 2 ), a second dual thermode stimulus (temperatures swapped), a 15s-15s-45s T 1 -T 2 -T 1 complex stimulus, a 75 s T 1 stimulus and finally another 15 s T 2 stimulus. The dual thermode stimuli were timed such that their superposition would match the profile of the complex stimulus. Relative stimulus intensities were fixed, but absolute intensities were adjusted for pain thresholds and tolerances (T 1 = [44.5 °C, 45 °C, 46 °C, 47 °C], N = [1, 5, 7, 7], resp.). Spatial configurations of thermodes (near, far or opposite arms, Fig. 4) were changed between sequences, and the order of sequences was randomized across subjects. Stimuli were delivered to the left ventral forearm, except for the opposite arm stimulus, which was delivered to a site on the right ventral forearm unaffected by the training task. 6 °C/s temperature ramp rates were used. Data analysis. Preprocessing. Matlab was used to resample thermode data to match the finger device sampling rate using heartbeat signals to determine appropriate timestamps. Stimulus-response lag was determined by the maximal cross correlation of stimulus and response for single thermode boxcar stimuli (average lag: 3.36 ± 1.6 s, mean ± SD). Finally, onset and offset times for responses to specific stimulus features of interest Scientific RepoRts | 7: 3894 | DOI:10.1038/s41598-017-04009-9 (i.e. t 1 , t 2 and t 3 ) were projected by discounting the temperature rise and fall rates and incorporating the average stimulus-response lag.
Hierarchical Linear Modelling. Inference on mean TCE was drawn using hierarchical linear models in order to control for between subject covariates (e.g. stimulus intensity) while also modelling variability within subject. Pain intensity (VAS) served as the dependent variable, measured at three timepoints per stimulus epoch (Fig. 6A, insert). Several considerations informed model selection. First, the study's a priori hypothesis calls for a non-monotonic response characterized by two periods of high pain interrupted by a period of lower pain. This suggests a specific planned contrast (a and c vs. b in Fig. 6A). Second, the relationship between stimulus intensity and TCE (if any) is unknown and is modeled. Third, previous studies have shown site specific and nonspecific sensitization depend on stimulus intensity 24 , but neither affect magnitude of offset analgesia 32 (a special case of TCE), therefore sensitization and its interaction with stimulus intensity, but not with TCE, are modeled. Finally, within subject stimulus intensity is constant across all timepoints of interest (a, b, c), therefore stimulus intensity and its higher order effects are only modeled at the group level. At the subject level the only random effects modeled are TCE and sensitization. The full model is summarized by the following equations, where X is the higher level design matrix, and Z is a block diagonal design matrix composed of Z i , the subject level design matrices. S is the stimulus number, since we delivered two, and models sensitization effects. TCE ab-c is the planned contrast. The model was fit using restricted maximum likelihood estimation as implemented in the LME4 package in R. The lmerTest package was used for inference with Satterthwaite's approximation for effective second level degrees of freedom. Because TCE is tested in three different thermode configurations (near, far and opposite arms), Holm-Sidak correction for multiple comparisons is employed to control the 0.05 false positive rate. The higher level design matrix fits are reported in (Table 3). Although plausible a priori, many of these higher order parameter estimates are not statistically significant. Stepwise regression strategies would argue for their removal, but such methods can be problematic in hierarchical models due to the effect of higher level terms on lower level parameter estimates. Terms that are present at the first level of the model are not considered for removal due to the difficulty of correctly quantifying the effective reduction in model complexity, and such terms are also retained at the second level, but Bayes information criteria is used to determine which temperature terms to remove, since these are only present at the second level (Table 4). Nevertheless, because the correct approach to such model selection problems remains an open question in the statistical literature 33 , and different approaches might have justified a different final model selection we also present the preliminary models.
Dynamics of nonlinearity. Timeseries data was analyzed using an autoregressive integrated moving average (ARIMA) model 34 ARIMA models capture linear relationships between points within timeseries (u(t)), producing a modified error term (ε) against which the statistical significance of parameter estimates between timeseries (B) can be evaluated (L i , lag operator; ϕ i , autocorrelation weight of u(t) vs u(t-i); θ i , moving average weights).
To determine the error model order (p, d, q) we regress independent variables from the dependent variables and inspect the autocorrelation and partial autocorrelation plots of the residuals. These revealed a parsimonious solution with timeseries dominated by lag 1 autocorrelation at the level of the first difference and no moving average component, an ARIMA(1, 1, 0) model. Analysis of Bayesian information criteria for models with autoregressive and moving average terms within 2 orders of an ARIMA(1, 1, 0) model supports this conclusion. Thus all effect sizes for parameters of interest are evaluated while controlling for ARIMA(1, 1, 0) errors. This yields more conservative inferences than would result from a standard regression model. Data availability. Raw and preprocessed data used in this study are available from the corresponding author upon request.