Ultra-low frequency neural entrainment to pain

Nervous systems exploit regularities in the sensory environment to predict sensory input and adjust behavior, and thereby maximize fitness. Entrainment of neural oscillations allows retaining temporal regularities of sensory information, a prerequisite for prediction. Entrainment has been extensively described at the frequencies of periodic inputs most commonly present in visual and auditory landscapes (e.g. >1 Hz). An open question is whether neural entrainment also occurs for regularities at much longer timescales. Here we exploited the fact that the temporal dynamics of thermal stimuli in natural environment can unfold very slowly. We show that ultra-low frequency neural oscillations preserved a long-lasting trace of sensory information through neural entrainment to periodic thermo-nociceptive input as low as 0.1 Hz. Importantly, revealing the functional significance of this phenomenon, both power and phase of the entrainment predicted individual pain sensitivity. In contrast, periodic auditory input at the same ultra-low frequency did not entrain ultra-low frequency oscillations. These results demonstrate that a functionally-significant neural entrainment can occur at temporal scales far longer than those commonly explored. The non-supramodal nature of our results suggests that ultra-low frequency entrainment might be tuned to the temporal scale of the statistical regularities characteristic of different sensory modalities.


Introduction 45
The sensory environment is dynamic in nature, with its temporal structures unfolding across multiple 46 timescales. Time is therefore an indispensable aspect of sensory experiences. The ability of the brain to track 47 and predict the temporal dynamics of sensory inputs allows an organism to take appropriate actions to meet 48 the changing environmental demands. What neural processes may be responsible for these functions, 49 however, is still an open question. Accumulating evidence supports a theory that neural codes of temporal 50 information build on brain oscillations (1-3). Taking situations involving rhythmic sensory inputs as an 51 example, the brain may adapt to the external rhythm through entrainment of ongoing neural oscillations at 52 the corresponding frequency (2, 4-6). Thus, neural entrainment constitutes a flexible mechanism through 53 which the brain adjusts the power or the phase of ongoing oscillations as a function of sensory input, with 54 consequences on brain dynamics that can persist after the sensory input has ceased (7, 8). 55 Temporal dynamics are ubiquitous in sensory domains, including pain (9-13). Despite this fact, most 56 neuroimaging studies investigating the neural mechanisms of pain were conducted using transient painful 57 stimulation (14-16). This approach poses two main problems. First, the neural responses to transient painful 58 stimulation are dominated by supramodal neural activities -i.e. activities associated with the detection of 59 salient environmental events regardless of their modality (17-19) -which limits the usefulness of this 60 approach in identifying nociceptive-specific brain processing (16). Second, the presentation of a single 61 intense painful stimulus does not capture the rich and often long-lasting dynamics of pain perception, 62 leaving the critical question of how the brain processes dynamic pain information largely unanswered. There 63 is, therefore, a growing consensus in the field that a shift is needed from measuring brain responses elicited 64 by transient painful stimulation to more naturalistic approaches that allow the capture of the temporal 65 dynamics of pain (15,16). intensity (22,25). In other words, the slow temporal dynamics of both nociceptive input and pain seem to 74 be reflected in neural oscillations at frequencies several orders of magnitude higher than that of ongoing 75 pain. 76 The above considerations lead to an outstanding question about the role of low-frequency neural 77 oscillations: can the slow dynamics of pain be represented in neural oscillations at the same timescales? This 78 is physiologically plausible, for several reasons. First, a 1:1 phase locking between the rhythm of external 79 sensory inputs and neural oscillations is theoretically possible (26-28), and has been in fact repeatedly 80 observed in auditory and visual domains in the form of neural entrainment (7, 29-33). In these cases, the 81 neural oscillations, and specifically their phase structures, may serve as a substrate of temporal 82 representation and prediction of the incoming input (1, 2). Second, even when the power of higher 83 frequency oscillations is modulated at the lower frequency of the sensory input, oscillations at the same 84 frequency of the input can still be functionally important: they may coordinate rapidly-changing and local 85 neural processing, usually reflected by high-frequency brain activities, with brain activities operating at 86 slower and more behaviorally-relevant timescales (1, 34, 35). 87 Therefore, in the current study, we investigated whether the brain encodes long-timescale dynamics of 88 painful input through entrainment of ultra-low frequency neural oscillations. Specifically, we investigated (i) 89 whether and how amplitude and phase modulations contribute to such entrainment, (ii) whether this 90 entrainment is supramodal, and (iii) whether the strength of entrainment reflects the variability of 91 perceptual sensitivity across participants. 92 To address these questions, we recorded high-density EEG (128 electrodes) from participants receiving 93 continuous thermo-nociceptive and continuous loud auditory stimuli oscillating at 0.1 Hz. Painful stimuli 94 were delivered over the hand dorsum using a feedback-controlled laser stimulator with high temporal and 95 temperature precision. Participants were requested to rate the perceived stimulus intensity on a visual 96 analogue scale (VAS). To control for the confounding effect of the rating task, in two additional conditions 97 participants received painful and auditory stimuli but were not required to rate them. Finally, to control for 98 stimulus-intensity effects, we included an additional condition in which painful stimuli of lower intensity 99 were delivered. Thus, there were five conditions in total. 100 Our results provide clear evidence that the long-timescale dynamics of nociceptive input are encoded by 101 neural oscillations at the same ultra-low frequency of the input. The observations that this encoding was 102 dependent on the phase of ongoing neural oscillations, and particularly that it outlasted the stimulus input, 103 imply a true neural entrainment mechanism. This ultra-low frequency entrainment was not supramodal, as 104 it was robust during nociceptive stimulation but not present during auditory stimulation of similar intensity. 105 Remarkably, the strength of neural entrainment to the nociceptive input was predictive of pain sensitivity 106 across individuals. 107 108

Stimulus input and perceptual ratings have similar temporal profiles 110
Participants' continuous rating of perceived intensity roughly followed the temporal pattern of the rhythmic 111 stimulation. In both pain (Fig 1A) and auditory (Fig 1B) rating time series, we observed three cycles whose 112 period was similar to that of the stimulus. We formally compared the peak amplitude and latency of ratings 113 across conditions and cycles. The results (Fig 1C and Supplementary Information) showed that auditory 114 ratings peaked earlier than pain ratings in all three cycles, and that pain ratings were decreased and delayed 115 in the last two cycles compared to the first one (Fig 1C). 116 Low-frequency nociceptive input enhances neural activity at the frequency of the stimulus 117 Next, we examined whether brain activities encoded the low-frequency rhythmic stimulation. When 118 inspecting the time-domain EEG responses to the nociceptive input, we observed a clear oscillatory pattern 119 reminiscent of that of the stimulus at central electrodes (Fig 2A). Importantly, these EEG oscillations were 120 not visible in response to the auditory stimulation. To quantify the frequency contents of these EEG 121 responses, we transformed single-trial EEG signals into the frequency domain and then averaged the 122 resultant power spectra across trials, for each subject and condition (Fig S1). Note that we performed the 123 frequency decomposition at trial level rather than on single-subject average waveforms because in the 124 former case it is possible to detect an increase of power regardless of whether neural activity is phase-125 locked across trials . To reveal the frequency of power increases, we subtracted the average power of  126  neighboring frequencies from each given frequency, yielding background-subtracted power (BSP), as  127 previously recommended (20, 36, 37). We found strong evidence for a power enhancement only at 0.1 Hz 128 for all three conditions entailing nociceptive stimulation (Fig 2B; all t29>7.802, P<0.0001, Cohen's d>1.424). 129 This effect was maximal in central scalp regions (Fig 2B; also see Fig S2 for detailed results at other  130 frequencies and scalp positions). This power enhancement could be consistently detected across single trials 131 in the majority of subjects (Fig 2D). 132 We observed strong evidence that the BSP at 0.  (Fig 2C), indicating that the power enhancement was 136 not dependent on the rating task. 137 These results showed an enhancement of neural activities specifically at the frequency of rhythmic painful 138 stimulation. Importantly, this effect was neither supramodal (i.e., it was not present in audition), nor 139 dependent on whether the participants were performing a rating task. Also, the observation of an increase 140 of BSP calculated using an extremely narrow frequency window of 0.06 Hz (i.e. from -0.03 Hz to +0.03 Hz 141 with respect to 0.1 Hz) strongly indicates that the observed power increase is not consequent to a general 142 autonomic response, but has instead a neural origin. 143

Phase reorganization of neural activity at the frequency of the nociceptive stimulus 144
The finding that neural activity and stimulus profile have the same peak frequency does not necessarily 145 indicate a stable phase relationship between the two (38). To test for such a phase relationship, we 146 quantified the EEG phase locking across trials using inter-trial phase coherence (ITPC), separately for each 147 subject and condition. Since EEG trials were time aligned to the onset of the rhythmic stimulation, it follows 148 that here ITPC also quantified the consistency of the phase relationship between stimulus and EEG data 149 within each individual. We observed that during painful, but not auditory, stimulation there was a clear peak 150 of ITPC at 0.1 Hz. This was observed in both the mean ITPC across subjects (Fig 3A) and the percentage of 151 subjects with significant ITPC (Fig 3B; ITPC significance was determined in each subject using the Rayleigh's 152 test for circular uniformity (39)). This phase-locking effect at the frequency of the nociceptive stimulation 153 was maximal in central scalp regions (Fig 3A-B). To test for the group-level consistency of this effect, we 154 compared the group-mean ITPC and the percentage of subjects with significant ITPC to randomized data 155 (see Materials and Methods). In all pain conditions, both ITPC measures at 0.1 Hz in central regions were 156 consistently greater than chance (all P<0.001) (Fig 3A-B oscillations. Similar to the power increase, this phase-locking effect was also not supramodal, and not 186 related to whether the participants were performing a rating task. 187

Pain sensitivity across individuals is reflected in the strength of neural entrainment 188
Perceived pain was stronger in participants with a stronger neural entrainment. Specifically, we correlated 189 the pain rating time series to the indices of power enhancement (BSP) and phase locking (ITPC) at 0.1 Hz in 190 central scalp regions. We observed clear positive correlations in the time interval around each rating peak 191 (Fig 4). Remarkably, the across-subject variability in perceived pain intensity was also reflected in the phase 192 relationship between the entrained oscillations and the stimulus (Fig 5). To evaluate this relationship, we 193 fitted a single-cycle cosine to the subject-mean peak rating as a function of the phase of 0.1-Hz oscillation in 194 central scalp regions. In subjects who rated the nociceptive stimulus as more painful the phase of neural 195 activity at 0.1 Hz was closer to that of the stimulus input. This is an important finding given that almost all 196 nociceptive-evoked neural responses fail to track pain sensitivity across subjects, in both human and animal 197 studies (42). 198 Importantly, all these relationships (i.e. the correlation between BSP/ITPC and pain ratings, and the 199 relationship between phase and pain ratings) existed not only (i) within the conditions entailing nociceptive were observed in the conditions entailing auditory stimulation (Fig 4G-H; Fig 5G-H). 206 To test whether the variability in stimulus temperature contributed to the above results, we also tested for 207 an across-subjects relationship between stimulus temperature and pain ratings, as well as between stimulus 208 temperature and neural entrainment. We did not observe such relationships (see Table S1 for detailed 209 statistical results). 210 Taken together, these results show that the strength of entrainment could predict how sensitive each 211 participant is to painful stimulations. 212

Stimulus-induced neural oscillations outlast nociceptive input 213
A remarkable observation was that neural oscillations around 0.1 Hz continued for at least 10 seconds after 214 the end of the rhythmic nociceptive input (Figs 2A and 6). The observation that the scalp topographies at 215 the peak of the additional cycle resembled those of the previous cycles (Fig 6) further indicates self-216 sustained activity of the same underlying neural process. These observations provide additional evidence for 217 an actual neural entrainment of ongoing EEG activities to the nociceptive input. They also provide further 218 evidence that the observed ultra-low frequency oscillations were not consequent to autonomic responses. 219 Finally, the amplitude of this additional oscillation after the end of rhythmic nociceptive stimulation was 220 correlated with pain ratings across subjects (Fig 6). 221

222
Discussion 223 Here we aimed to identify neural activity present during tonic sensory stimuli that produce slowly 224 fluctuating sensations. We delivered intensity-matched auditory and painful stimuli fluctuating at 0.1 Hz and 225 observed that only painful stimuli resulted in both a power enhancement and a phase locking of brain 226 activity at the same frequency of the stimulus (Figs 2 and 3). Thus, this brain response was not supramodal, 227 and possibly selective for the somatosensory system. This response likely reflected a true neural 228 entrainment, since (i) the degree of phase locking depended on the phase of ongoing brain oscillations 229 occurring before the onset of the rhythmic input (Fig 3D), and (ii) the stimulus-induced brain oscillations 230 outlasted the rhythmic input (Fig 6). Importantly, this neural entrainment to the rhythmic painful input was 231 not due to the rating task, as it was also present when participants did not have to rate the painfulness of 232 the stimuli (Figs 2 and 3). Finally, the strength of the neural entrainment was correlated with pain sensitivity 233 across individuals (Figs 4 and 5), a relationship that persisted even in the neural activity outlasting the 234 rhythmic stimulus (Fig 6). 235 These findings show that the brain encodes long-timescale dynamics of nociceptive input through 236 entrainment of ultra-low frequency neural oscillations. This work not only represents a step towards 237 analyzing brain processes in more clinically-relevant models of long-lasting and dynamic pain (15, 16), but 238 also sheds new light on the functional significance of neural oscillations at frequencies well below the 239 traditional boundaries of human EEG rhythms. 240

Neural entrainment or evoked responses? 241
The power enhancement (Fig 2) and phase-locking (Fig 3) of neural activities observed at the frequency of 242 nociceptive input could be explained by either repeated evoked responses or an entrainment of ongoing 243 neural oscillations to the rhythmic stimulus (6, 7, 40, 41, 43-47). Two lines of evidence in our data support an 244 entrainment mechanism. 245 First, we found clear evidence that the degree of phase locking depended on the phase of ongoing neural 246 oscillations occurring before stimulus onset (Fig 3D). Indeed, such phase interactions are predicted by 247 theoretical models and experimental data of entrainment of neural oscillators: although rhythmic stimuli 248 drive brain activities, the effectiveness of this process, reflected in the phase reorganization, also depends 249 on the state of the neural oscillators (33, 38, 48, 49). 250 Second, as depicted in Fig 6,  regularities, which could facilitate sensory processing of incoming events (1, 2). 261

Neural entrainment to rhythmic painful input is independent of rating task 262
We observed neural entrainment to the rhythmic painful input regardless of whether the participants were 263 required to continuously rate pain intensity (Figs 2, 3, 6, S2, and S3). 264 Continuous perceptual ratings are commonly used to investigate the neural correlates of percepts occurring 265 at long timescales (e.g., tonic experimental pain in healthy volunteers (20, 22, 25) or spontaneous 266 fluctuations of pain in chronic pain patients (12, 13)). Such continuous ratings are typically obtained with a 267 finger-span device (12, 13, 22, 25) or a slider (20). While such ratings provide valuable information on the 268 dynamic fluctuations of perception, they heavily confound the analysis of neural data due to the 269 superposition of brain activities related to the motor and cognitive activity related to the rating task (53). A 270 strategy to control for this confound is to ask participants to continuously rate the perceived intensity of 271 stimuli belonging to another sensory modality (e.g., vision) (12, 13, 22, 25, 54): the brain activity measured 272 during painful but not visual stimulation should then reflect pain-selective activity. This paradigm assesses 273 whether the rating task is sufficient to explain the brain response sampled during painful stimulation. 274 However, it cannot resolve whether the same brain response remains when no rating task is performed -275 that is, whether the rating task is necessary for observing the brain response. To effectively address this 276 issue, we included in our experiments additional control conditions in which participants had to rate the 277 intensity of neither the auditory nor the painful stimuli. 278 Thus, our results demonstrate that rating-related brain activities were not necessary for the observed 279 entrainment of brain oscillations to the rhythmic painful input. Indeed, the power enhancement, the phase 280 locking, as well as the continuation of the entrained neural oscillations after stimulus offset were also 281 present in the pain condition not involving the rating task (Figs 2, 3, 6, S2, and S3). 282 While we did not observe strong evidence for an effect of rating task at the stimulation frequency of 0.1 Hz 283 for either power or phase, the task of rating auditory stimuli enhanced the power and the phase locking at 284 0.2 and 0.3 Hz (Figs 2B, 3A-B and S3). In contrast, rating of painful stimuli only slightly enhanced the phase 285 locking at 0.2 and 0.3 Hz (Fig S3). This frequency and modality pattern suggests that the rating effect reflects 286 a different mechanism than the neural entrainment to the rhythmic input (44). The fact that the increases 287 and decreases were more regular in auditory than in pain ratings (Fig 1) might explain why the effect of 288 rating auditory stimuli was more evident. 289 Is the observed neural entrainment modality-specific? 290 As we discussed above, the entrainment of ultra-low frequency brain oscillations to stimulus input was 291 clearly not supramodal, given that (i) the power increase of EEG signal at stimulus frequency was only 292 consistently observed during nociceptive stimulation but not during loud auditory stimulation (Fig 2), and (ii) 293 only the EEG signals during nociceptive stimulation showed a predominant peak of phase locking at stimulus 294 frequency (Fig 3). Importantly, these findings by no means imply that the entrainment we observed during 295 nociceptive stimulation was modality-specific. Indeed, strictly speaking, demonstrating the pain specificity of 296 a neural response is virtually impossible, as it would require testing all stimuli that could in principle elicit 297 that response, and show that that response only occurs when pain is experienced (16). Instead, the current 298 findings provide evidence that entrainment at the ultra-low frequency used in this study occurs 299 preferentially in response to somatosensory input. That the observed entrainment is preferential to pain 300 would require testing whether it occurs less strongly during non-nociceptive somatosensory stimulation. 301 The lack of entrainment to the auditory stimuli is not trivial, since there is a considerable amount of 302 evidence for entrainment of neural oscillations to rhythmic auditory stimulation (4, 7, 30, 31, 33, 46). A 303 possible explanation is that the frequency of the delivered auditory stimulation is substantially lower than 304 the timescale of dynamics for which the auditory neural circuits are optimized. Indeed, the temporal 305 structures of speech and music largely occur in the subsecond range (3, 55). Accordingly, previous evidence 306 for auditory entrainment is primarily observed during stimulation at delta (1-4 Hz) and theta (4-8 Hz) 307 frequencies (7, 8, 30, 33, 56, 57). 308 These previous observations, together with the current results, suggest a different frequency preference for 309 oscillatory entrainment across sensory systems (58). Why would the nociceptive system respond to lower 310 ranges of stimulus frequencies than the auditory system? Given that the temporal dynamics of fluctuations 311 of spontaneous pain and somatosensory detection occur at very low frequencies (11, 59), the brain 312 representations of long-timescale variability might be more essential and relevant to the processing of pain 313 information. These considerations suggest that neural entrainment is tuned to the temporal scale of the 314 statistical regularities characteristic of different sensory modalities, a hypothesis that would require further 315 testing. 316

The strength of neural entrainment reflects pain sensitivity across individuals 317
We observed a clear relationship between perceived pain intensity and the strength of neural entrainment 318 to the nociceptive input. Individuals with higher pain sensitivity had greater power enhancement and phase 319 locking at the frequency of the nociceptive input (Fig 4), as well as a phase of the entrained oscillations 320 closer to the phase of the nociceptive input (Fig 5). This is a remarkable result, for two reasons. First, even 321 the commonly-observed within-subject correlation between nociceptive-evoked responses and subjective 322 pain ratings have been shown to be not obligatory, i.e. they can be easily disrupted by a large number of 323 experimental manipulations (e.g. expectation, stimulus repetition, presentation of non-painful but iso-324 salient stimuli, congenital insensitivity to pain; (60-62)). Second, almost all nociceptive-evoked neural 325 responses fail to track pain sensitivity across subjects, in both human and animal studies (42). 326 It is important to note that the relationship between entrainment and pain sensitivity across subjects were 327 not driven by the rating task, and were also present when examining conditions entailing different 328 intensities of nociceptive stimulation. Most notably, the relationship between oscillatory amplitude and 329 ratings persisted following the end of rhythmic stimulation (Fig 6). continuously rate their perceived stimulus intensity on a visual analogue scale (VAS) ranging from 0 to 10 365 using a custom-built vertical slider controlled with their right hand. The slider was connected to a 366 potentiometer to record their ratings. For nociceptive stimuli, the lower and upper ends of the slider were 367 defined as "no pain" and "the maximum pain tolerable" respectively. For the auditory stimuli, they were 368 defined as "no sound" and "the loudest sound tolerable". Ratings were recorded from the onset of the 369 periodic stimulation, and lasted for 30 s in auditory trials and 45 s in pain trials. Rating data was digitized at 370 1024 Hz (USB-1408FS, Measurement Computing Corporation, Norton, MA, USA) and synchronized with 371 stimulation triggers and EEG recordings. In all trials, participants were instructed to focus their attention on 372 the stimuli and keep their gaze on a fixation cross placed centrally in front of them, at a distance of 373 approximately 60 cm and 30° below eye level. The order of stimulus presentation and rating task is detailed 374 in Table S2. The experimenters started trials manually after ensuring that the participant was ready and 375 instructed whether to rate the sensation. As a result, the time between the end of a trial and the beginning 376 of the following trial ranged between 10.0 and 78.2 s (average 21.2 s). Participants were allowed to rest for 377 approximately 2 minutes after every 15 trials. 378

Sensory stimuli 379
Nociceptive tonic stimuli were generated by a temperature-controlled CO2 laser stimulator (Laser 380 Stimulation Device, SIFEC, Ferrières, Belgium) with a wavelength of 10.6 μm and a beam diameter of 6 mm. 381 The output power of the laser was continuously regulated by a feedback control loop based on an online 382 monitoring of skin temperature at the target site (67). The laser stimulation target was changed after each 383 stimulus to avoid nociceptor fatigue, sensitization, and skin damage. Laser power modulation resulted in a 384 0.1-Hz sinusoidal modulation of skin temperature, starting with an initial phase of π (i.e., trough), and lasting 385 for 30 s (Fig 1A). The temperature difference between peaks and troughs was 3°C. Each participant received 386 laser stimuli of two intensity levels (high-pain and low-pain stimuli), individually adjusted as described below. 387 To avoid saliency-related brain responses during the periodic stimulation, the skin temperature was first 388 brought to the desired trough level in a 1-s heating ramp and then maintained at this level for 5 s before the 389 onset of the periodic stimulation. After the periodic stimulation, the skin temperature was maintained at the 390 trough level for 10 additional seconds. 391 Before the experiment, stimulus temperatures were determined individually as follows. First, pain detection 392 threshold and pain tolerance were estimated. To measure the pain detection threshold, participants 393 received linearly increasing stimuli at 1°C/s (with a cutoff at 54°C) on the dorsum of their left hand, and were 394 instructed to press a button with their right hand as soon as they felt a painful sensation. The button-press 395 immediately terminated the stimulation. To measure the pain tolerance, participants were instructed to 396 press the button as soon as they felt the painful sensation become intolerable. The laser target was changed 397 after each stimulus. Both the pain detection threshold and pain tolerance were measured three times, in 398 consecutive trials, and their corresponding temperatures were estimated as the mean of the three 399 consecutive measurements. The trough temperature of the low-pain stimuli was set to 1°C above the pain 400 detection threshold, and the trough temperature of the high-pain stimuli was set to 1°C above that of the 401 low-pain stimuli. The peak temperature of the stimulus was always below the pain tolerance. The resulting 402 peak temperature in the high-pain stimuli was 48.3 ± 1.9°C (mean ± SD across participants). 403 Auditory stimuli were generated by MATLAB (MathWorks, Natick, MA, USA) and presented through 404 headphones binaurally. Auditory stimuli consisted of a pure tone (frequency of 280 Hz) whose amplitude 405 was sinusoidally modulated at 0.1 Hz for 30 s (Fig 1B). As for the nociceptive stimulation, the sinusoidal 406 modulation started with an initial phase of π (i.e., trough). The sound intensity was individually adjusted to 407 ensure that perceived intensity was similar to the high-pain condition. Auditory stimuli were presented using 408 the Psychophysics Toolbox (68). 409

Analysis of subjective sensations 410
Single-trial rating time series were down-sampled to 512 Hz and smoothed using a moving mean filter with a 411 1-s window. Filtered data was averaged across trials, for each subject and condition. The peak latency and 412 amplitude of each rating cycle were measured from the average waveforms. The rating peak latencies were 413 measured with respect to the corresponding peak latencies in the stimulus time series (i.e., 5, 15, and 25 s; 414 Fig 1) simultaneously recorded using two surface Ag/AgCl electrodes, one placed below the lower eyelid and one 423 laterally to the outer canthus of the right eye. Electrode impedances were kept below 10 kΩ. 424 EEG preprocessing was conducted using the EEGLAB toolbox (69) and custom-written MATLAB scripts. 425 Continuous EEG data was notch filtered at 50 Hz and harmonics to remove power line noise, and then 426 segmented into epochs ranging from -10 s to 45 s relative to the onset of the sinusoidal stimulation. 427 Differences in EEG baseline across trials were removed by demeaning each trial. The EEG data was re-428 referenced to the common average. Signals contaminated by eye blinks, eye movements, or muscle 429 activities were corrected using independent component analysis (70). Trials containing excessive signal 430 fluctuations in at least one electrode (amplitude exceeding ± 500 μV) were excluded from further analyses. 431 These trials constituted 4.6% of the total number of trials. The corresponding trials in the rating data were 432 also excluded. 433 Both frequency and phase measures of entrained neural oscillations can be confounded by transient 434 "evoked-type" responses that repeat at the stimulus frequency (6, 7, 41, 44-47). We did observe a transient 435 response (lasting ~0.5 to 2 s) locked to the increase of auditory stimulation intensity (Fig S4). Although 436 interesting, this response could contaminate our measures of ultra-low frequency neural oscillations at the 437 stimulus frequency. We therefore applied a cascade of filters at specifically-defined scales in the time 438 domain to both pain and auditory trials, to minimize the potential confounding effects of such regularly 439 occurring transient responses. We first denoised single-trial data after the above processing steps (s0) using 440 a moving mean filter with a 0.5-s window (s1). We then applied a 2-s median filter and a Gaussian filter with 441 full-width at half-maximum (FWHM) of 1 s to s1, yielding a signal (s2) in which the transient responses were 442 removed while the long-period signals were kept. Finally, we reconstructed the EEG signal (s', without the 443 transient responses) as s' = s0 -s1 + s2. This algorithm was effective in removing the transient responses 444 while leaving other features of the signal largely intact (Fig S4). Thus, the power increase and phase locking 445 of the EEG responses revealed by the following analyses were most likely due to a true 0.1-Hz oscillation 446 rather than transient responses that repeated at this frequency. 447

Power analysis of EEG data 448
A fast Fourier transform was applied to single-trial signals ranging from 0 to 30 s after the onset of the 449 sinusoidal stimulation, yielding power spectra with a frequency resolution of 0.0333 Hz (a frequency 450 resolution of 0.01 Hz, i.e., spectral interpolation, was achieved by zero padding in the time domain and was 451 used for illustrative purposes). Power estimates were log-transformed and averaged across trials for each 452 subject and condition. To reveal the frequency of power increases, for each subject, condition, electrode 453 and frequency point, the contribution of background activities (e.g., spontaneous brain activities or slow eye 454 movements) was removed by subtracting the average power at surrounding frequencies (-0.0333 Hz and 455 +0.0333 Hz) (20, 36, 37). Scalp topographies of this background-subtracted power (BSP) were computed by 456 spline interpolation. 457 To identify the frequencies at which power increase occurred, a one-sample one-tailed t-test was performed 458 at each frequency point to test whether the BSP was consistently greater than zero across subjects (FDR 459 corrected for multiple comparisons across frequencies). This analysis was first performed on signals from a 460 central electrode cluster (Cz and its closest neighbours FCC1h, FCC2h, CCP1h, and CCP2h) and then extended 461 to all electrodes (FDR corrected for multiple comparisons across frequencies and electrodes). An additional 462 one-sample one-tailed t-test was performed separately for each subject and condition, to examine whether 463 the BSP at 0.1 Hz in the central electrode cluster was consistently greater than zero across single trials. 464 To test the effects of modality and rating task on the power increase detected at 0.1 Hz, we conducted, for 465 each electrode, a two-way repeated-measures ANOVA with factors Modality (two levels: high pain and 466 sound) and Rating (two levels: rating and no rating) on 0.1-Hz BSP (FDR corrected for multiple comparisons 467 across electrodes). Post hoc paired-sample two-tailed t-tests were performed when a significant main effect 468 or interaction was found. 469

Phase analysis of EEG data 470
We examined the phase locking of the EEG signal (0-30 s) across trials by calculating the inter-trial phase 471 coherence (ITPC) (71) for each subject, condition, and electrode. Briefly, given the Fourier phase φ n for trial 472 n, we define the mean vector of phase angles across trials as m = −1 ∑ =1 , where N is the number of 473 trials. The ITPC value is given by the modulus of m, i.e., ITPC = |m|. To determine at which frequencies 474 phase locking occurred, we evaluated the ITPC as a function of frequency by calculating the ITPC values in 475 steps of 0.0333 Hz (a frequency resolution of 0.01 Hz was achieved as described in the previous paragraph). 476 ITPC scalp topographies were computed by spline interpolation. For each subject, the significance of ITPC 477 was determined using the Rayleigh's test for circular uniformity (P<0.05) (39). The percentage of subjects 478 with significant ITPC was calculated for each frequency point, electrode, and condition. 479 To further assess the significance of the ITPC at the group level (i.e., whether ITPC was greater than what 480 one would expect by chance), the mean ITPC across subjects and the percentage of subjects with significant 481 ITPC were compared to randomized data. Specifically, we added random phase values drawn from circular 482 uniform distribution to the single trial phases, recalculated the ITPC for each subject, and determined its 483 significance using the Rayleigh's test. We then computed the mean ITPC across subjects and the percentage 484 of subjects with significant ITPC. This process was repeated 1,000 times, yielding null distributions of the 485 mean ITPC and of the percentage of subjects with significant ITPC. P values of the actual data were 486 determined by comparing the mean ITPC and percentage of subjects to the respective null distributions. This 487 analysis was also first conducted for the 0.1-Hz oscillation in the central electrode cluster, and then 488 extended to other frequencies and electrodes (FDR corrected for multiple comparisons across frequencies  489 and electrodes). 490 To test the effects of modality and rating on the ITPC at 0.1 Hz, for each electrode, we performed a two-way 491 repeated-measures ANOVA with factors Modality (two levels: high pain and sound) and Rating (two levels: 492 rating and no rating) (FDR corrected for multiple comparisons across electrodes). Post hoc paired-sample 493 two-tailed t-tests were performed when a significant main effect or interaction was found. 494 To examine the dependence of the degree of phase locking on the phase of oscillations occurring before 495 stimulus onset, we applied a causal, linear-phase bandpass FIR filter with cutoff frequencies at 0.05 and 0.15 496 Hz to single-trial EEG signals from the central electrode cluster, and then extracted the instantaneous phase 497 of the filtered signals using the Hilbert transform. The causal filter was used to avoid the influence of signals 498 after stimulus onset on the prestimulus phase. We sorted the single trials into six bins from 0 to 2π 499 according to the instantaneous phase at time 0 s in the filtered signal. Within each subject, we pooled rating 500 and no-rating trials from the same modality, to increase the number of trials in each bin. We then calculated 501 0. introduced a delay in the filtered signal, we shifted the phase bins accordingly (by ½ cycle of a 0.1-Hz 507 oscillation) to estimate the relationship between the ITPCz and the instantaneous phase at the onset of the 508 rhythmic stimulation (importantly, this was not a shift of the filtered signal and did not introduce non-509 causality). We performed a two-way repeated-measures ANOVA on the ITPCz with factors Bin (six levels: 510 equal-sized bins from 0 to 2π) and Modality (two levels: pain and sound). Post hoc paired-sample two-tailed 511 t-tests were performed when a significant main effect or interaction was found. 512

Analysis of relationship between the strength of neural entrainment and perceived pain intensity across 513
individuals 514 For each condition entailing a rating task, we calculated across-subject Pearson correlation coefficients 515 between the intensity rating at each time point in the rating time series and the 0.1-Hz BSP as well as the 516 0.1-Hz ITPC in the central electrode cluster, yielding time series of the correlation coefficient r and the P 517 value (FDR corrected for multiple comparisons across time points). To test whether the correlations were 518 consequent to the rating task, we performed the same correlation analyses between the intensity ratings 519 and the BSP/ITPC measured in the conditions without rating task. Finally, the same correlation analyses 520 were also performed between conditions entailing different intensities of nociceptive stimulation. 521 For each subject and condition, we calculated the phase of the entrained 0.1-Hz oscillations as the 522 orientation of the above-defined mean vector m (i.e., arg(m)). To evaluate the across-subject relationship 523 between the intensity rating and the phase of entrained oscillations, we fitted a single-cycle cosine to the 524 subject-mean peak rating (i.e., the peak rating averaged across the three cycles) as a function of the phase 525 of 0.1-Hz oscillations in the central electrode cluster. Significance of the cosine fit was estimated with 526 permutation testing: we randomly permuted the phase values across subjects, fitted a cosine function, and 527 calculated the coefficient of determination R 2 as a measure of the goodness of fit. This procedure was 528 repeated 1,000 times, yielding null distribution of the R 2 . The P value of R 2 obtained from the actual data was 529 determined by comparing it to the null distribution. This analysis was performed between intensity rating 530 and 0.1-Hz phase within each condition entailing the rating task, between intensity rating and 0.1-Hz phase 531 in the conditions without rating task, and between conditions entailing different intensities of nociceptive 532 stimulation. 533 To ensure that the results from the above analyses were not due to individual variability in stimulus 534 temperature, we performed similar analyses but using stimulus temperature instead of ratings, and also 535 analyzed the correlation between stimulus temperature and pain ratings. Thus, we analyzed across-subject 536 relationships between the laser stimulation temperature and (i) mean peak pain rating (i.e., the peak rating 537 averaged across the three cycles), (ii) the BSP, (iii) the ITPC, and (iv) the phase of 0.1-Hz oscillations in the 538 central electrode cluster. 539

Analysis of neural oscillations outlasting the stimulus 540
The EEG time series from each subject was smoothed by a moving mean filter with a 2-s window, and 541 linearly detrended. No zero padding was used at the signal edges, to ensure that any oscillation after 542 stimulus offset was not due to the temporal smoothing. A one-sample two-tailed t-test of EEG amplitude 543 against zero was performed at each point of the time series (FDR corrected for multiple comparisons across 544 time points). Scalp topographies of the t value were computed over a 1-s window centered around the peak 545 and trough of each cycle. Finally, Pearson correlation coefficients were calculated across subjects between 546 the mean peak rating and the peak-to-trough amplitude of the cycle occurring after the sinusoidal 547 stimulation in the central electrode cluster. 548  was smoothed with a moving mean filter with a 2-s window, linearly detrended, and finally averaged across 783 subjects (N=30) in the conditions entailing high-pain stimulation without rating task (B), high-pain 784 stimulation with rating task (C), and low-pain stimulation with rating task (D). Shaded regions indicate the 785 time windows in which the EEG amplitude is significantly different from 0 (P<0.05, point-by-point one-786 sample t-test against 0, FDR corrected across time points). Scalp maps show the t value topographies within 787 1-s window around the peak and trough of each cycle. The similarity of scalp topographies at the peak of the 788 cycle after the end of rhythmic stimulation with those of the previous cycles further indicates self-sustained 789 activity of the same underlying neural process. Plots on the right show that the amplitude difference 790 between the peak and trough of the cycle after the end of rhythmic stimulation was correlated to the mean 791 of the peak pain ratings across subjects. 792 793  Table S1. Across-subjects relationship between stimulation temperature and pain ratings, as well as between stimulation temperature and different features of neural entrainment at 0.1 Hz