Abstract
Both passive tactile stimulation and motor actions result in dynamic changes in beta band (15–30 Hz Hz) oscillations over somatosensory cortex. Similar to alpha band (8–12 Hz) power decrease in the visual system, beta band power also decreases following stimulation of the somatosensory system. This relative suppression of α and β oscillations is generally interpreted as an increase in cortical excitability. Here, next to traditional single-pulse stimuli, we used a random intensity continuous right index finger tactile stimulation (white noise), which enabled us to uncover an impulse response function of the somatosensory system. Contrary to previous findings, we demonstrate a burst-like initial increase rather than decrease of beta activity following white noise stimulation (human participants, N = 18, 8 female). These β bursts, on average, lasted for 3 cycles, and their frequency was correlated with resonant frequency of somatosensory cortex, as measured by a multifrequency steady-state somatosensory evoked potential paradigm. Furthermore, beta band bursts shared spectro-temporal characteristics with evoked and resting-state β oscillations. Together, our findings not only reveal a novel oscillatory signature of somatosensory processing that mimics the previously reported visual impulse response functions, but also point to a common oscillatory generator underlying spontaneous β bursts in the absence of tactile stimulation and phase-locked β bursts following stimulation, the frequency of which is determined by the resonance properties of the somatosensory system.
SIGNIFICANCE STATEMENT The investigation of the transient nature of oscillations has gained great popularity in recent years. The findings of bursting activity, rather than sustained oscillations in the beta band, have provided important insights into its role in movement planning, working memory, inhibition, and reactivation of neural ensembles. In this study, we show that also in response to tactile stimulation the somatosensory system responds with ∼3 cycle oscillatory beta band bursts, whose spectro-temporal characteristics are shared with evoked and resting-state beta band oscillatory signatures of the somatosensory system. As similar bursts have been observed in the visual domain, these oscillatory signatures might reflect an important supramodal mechanism in sensory processing.
Introduction
Oscillatory dynamics are assumed to play a critical role in the processing of perceptual information (Buzsaki, 2006), with different sensory modalities operating preferentially in different frequency bands (Klimesch et al., 2007; Spitzer and Haegens, 2017; Kayser, 2019). In response to tactile stimulation, as well as during movement preparation and execution, the power of β oscillations (15–30 Hz) over somatosensory cortices decreases (Salmelin and Hari, 1994; Hari et al., 1997; Salenius et al., 1997; Neuper and Pfurtscheller, 2001; Parkkonen et al., 2015). This reduction in beta band power has generally been interpreted as an increase in excitability, allowing the somatosensory system to process external and interoceptive information during perception, motor coordination, and feedback (Cassim et al., 2000, 2001; Gaetz and Cheyne, 2006). Hence, β oscillations might have a similar role to that of α oscillations in the visual system, namely, cortical inhibition (Linkenkaer-Hansen et al., 2004; Klimesch et al., 2007; Engel and Fries, 2010; Jensen and Mazaheri, 2010; Jones et al., 2010; Haegens et al., 2011; Händel et al., 2011; Bonnefond and Jensen, 2012; Shin et al., 2017).
The interpretation of α and β as inhibitory brain rhythms is based on observations of sustained (∼500 ms duration) power decreases in trial-average waveforms following the stimulus (Hari et al., 1997; Neuper and Pfurtscheller, 2001; Parkkonen et al., 2015). A more nuanced spatiotemporal dynamic and functional role of low-frequency rhythms has been revealed by focusing on their transient nature. Single-trial analyses, in case of β oscillations in the somatosensory system, revealed a short-lived (∼150 ms, 3 cycles) burst-like behavior of β oscillations (Shin et al., 2017; Lundqvist et al., 2018; Little et al., 2019). Similar oscillatory transients have been uncovered in the visual system using reverse correlation between visual input and output (∼10 Hz, lasting up to 1 s) (VanRullen and Macdonald, 2012). Both of these approaches revealed a transient increase in oscillatory power hidden in global sustained decreases in power, which suggests a more active than traditionally assumed and temporally precise role of these rhythms (e.g., in sequence learning and predictive processing) (Huang et al., 2018; Alamia and VanRullen, 2019).
In this study, we aimed to investigate the transient nature and frequency tuning of β oscillations in the somatosensory system using several complementary approaches: (1) tactile random white noise (WN) stimulation (all frequencies in a 1–75 Hz range), which allowed us to compute participant-specific impulse response function (IRF) of the somatosensory system; (2) frequency-by-frequency rhythmic tactile stimulation in a 12–39 Hz frequency range (steady-state somatosensory evoked potentials [SSSEPs]), which allowed us to estimate the resonance, or maximal response amplitude, frequency of somatosensory cortex for each participant; (3) pulsed tactile stimulation, a typical event-related potential protocol; and (4) spontaneous prestimulus beta band activity. This broad approach allowed us to compare multiple β oscillatory signatures within somatosensory cortex and thus provide comprehensive evidence for the link between all of them.
We found that somatosensory IRFs, just like their α-rhythmic visual counterparts, contained β-frequency components and were unique to each participant in frequency (mean = 24.78 Hz; SD = 6.18 Hz). Tactile stimulation elicited β bursts that lasted on average for ∼3 cycles (mean = 2.81 cycles; SD = 0.35 cycles). This is in line with reports of transient 3-cycle-long movement-related β bursts (Lundqvist et al., 2016; Little et al., 2019). The peak frequency of somatosensory IRFs was strongly correlated with the somatosensory resonance frequency determined using an SSSEP approach, as well as with the peak frequency of stimulus-evoked and spontaneous β oscillations.
Our findings point to a common oscillatory generator underlying spontaneous β bursts in the absence of tactile stimulation and phase-locked β bursts following stimulation, the frequency of which is determined by the resonance properties of the somatosensory system. We propose an active role of beta band bursts in processing tactile information.
Materials and Methods
Participants
Eighteen participants (aged 22-41 yr, 8 females) with normal or corrected-to-normal vision enrolled in the study. The number of participants was chosen based on previous literature investigating somatosensory resonance (Snyder, 1992: N = 17; Tobimatsu et al., 1999: N = 10; Müller et al., 2001: N = 10; Moungou et al., 2016: N = 12). Informed consent forms were signed before the experiment. The experiment was conducted in accordance with the protocol approved by the Center National de la Recherche Scientifique ethical committee and followed the Code of Ethics of the World Medical Association (Declaration of Helsinki). Subjects were compensated with 10 Euro/h.
Experimental design
Tactile stimulation was delivered via vibrotactile electromagnetic solenoid-type actuator (Dancer Design) driven by an audio amplifier connected to a microcontroller board (Fig. 1A). We used amplitude-modulated stimulation using 150 Hz carrier frequency. Such approach is commonly used to study human touch perception, as time-varying stimulation strongly drives several types of mechanoreceptors (Moungou et al., 2016).
Stimulation protocol and task. A, Tactile stimulation was applied to participants' right index finger with a vibrotactile actuator. B, Three protocols were used: (1) WN stimulation was created by random amplitude modulation of the 150 Hz carrier sine-wave. The spectral content of WN sequences was flattened between 0 and 75 Hz. (2) SSSEP stimulation was created by rhythmic amplitude modulation of the 150 Hz carrier sine-wave at one specific frequency on each trial (12-39 Hz). Perceptual targets were embedded in both WN and SSSEP sequences by decreasing the frequency of the carrier sine-wave from 150 to 127.5 Hz. (3) Single-pulse stimulation consisted of 200 ms fixed amplitude increases of a carrier 150 Hz wave. C, Experimental procedure and stimulation protocols.
We used three vibrotactile stimulation protocols: (1) random WN stimulation, (2) steady-state rhythmic stimulation, and (3) short pulses (Fig. 1B). Tactile stimulation was applied to the participants' right index finger. WN sequences were created by randomly modulating the amplitude of a 150 Hz carrier sine-wave every cycle (i.e., every 1000/150 = 6.7 ms). The amplitude modulated signal was created by generating a random number series and normalizing the amplitude of its Fourier components before applying an inverse Fourier transform. Steady-state stimuli were created by amplitude-modulating a 150 Hz carrier frequency with a constant frequency sine-wave signal at one of the frequencies. In total, we used 18 stimulation frequencies covering the beta band frequency range and a few frequencies below and above: 12–20 Hz in steps of 2 Hz, 21–30 Hz in steps of 1 Hz, and 33–39 Hz in steps of 3 Hz. Each frequency was presented on 12 trials. The order of frequencies was randomized across trials.
Before the experiment, participants performed a practice session to become acquainted with a sensation of vibrotactile stimulation and target detection task. The task involved detecting a 1-s-long 15% decrease in carrier frequency (from 150 to 127.5 Hz). At the end of each trial, participants had to press with their left hand the appropriate keyboard button indicating whether the trial contained a target or not. The responses were not speeded and only 20% of trials contained targets.
Trial length for both WN and sine-wave stimulation was 6.7 s (6.6625 s). The practice session consisted of 5 sine-wave stimulation trials all containing targets. After the practice session, participants were first presented with 200 pulse stimuli (200 ms duration, 150 Hz vibration) separated by a random ITI of 1-1.5 s. Thereafter, participants completed 616 experimental trials: 400 trials of WN sequences and 216 trials of sine-wave stimulation. Experimental trials were divided into 8 blocks, with 80 trials each (7 blocks of 80 trials and 56 trials in the last block). The experimental part ended with an additional 200 single-pulse stimuli. To make sure that participants could not perform the target detection task by using auditory information (vibrotactile electromagnetic solenoid-type stimulators make a faint yet audible sound), they were using earplugs.
Data acquisition and preprocessing
We recorded participants' EEG using a 64-channel ActiveTwo Biosemi system. Data analysis was performed in MATLAB using the Fieldtrip toolbox and custom-written MATLAB scripts. Before all preprocessing steps, we visually identified and removed noisy channels. The EEG data were then rereferenced to the channel average, bandpass filtered between 1 and 200 Hz and line noise was removed using a DFT filter (50, 100, and 150 Hz). We performed independent component analysis to identify and remove eye movement-related artifacts. Thereafter, the data were epoched for SSSEP (−1000 to 8000 ms relative to stimulus onset), WN sequences (−1000 to 6625 ms relative to stimulus onset), and pulse stimulation (−400 to 600 ms relative to stimulus onset) trials separately. Baseline correction was applied by subtracting average values calculated from −1000 to 0 ms (WN), −200 to 0 ms (SSSEP), or −400 to 0 ms (ERP) from individual trials. For SSSEP analysis, EEG data were high-pass filtered at 0.5 Hz and downsampled to 512 Hz. Based on visual inspection, trials containing large muscle and head-movement related artifacts were removed.
Statistical analysis
ERP analysis
ERPs were calculated by averaging epochs from 400 ms before to 600 ms after stimulus presentation and applying absolute baseline correction (baseline window −200 to 0 ms).
Impulse Response Functions.
IRFs were calculated by first cross-correlating the z-scored single-trial EEG signal with the WN sequence presented at that particular trial and then averaging the result over trials. EEG trials were downsampled to 150 Hz before cross-correlation to match the rate of vibrotactile stimulation amplitude changes (Fig. 2A). To evaluate statistical significance of IRFs, we created a null hypothesis distribution by cross-correlating EEG signal with WN sequences from random trials (surrogate impulse response function). This randomization procedure was repeated 10,000 times, and the resulting trial-averaged IRFs served as the distribution of expected correlation values under the null hypothesis (H0: EEG and WN are uncorrelated). All stimulus time points were entered in the cross-correlation, except the first 500 ms and the last 1000 ms to avoid influences of the onset and offset ERP to the estimate of tactile IRF. The cross-correlation was calculated for lags between −300 and 500 ms as follows:
A, Cross correlation procedure. The cross-correlation between WN sequences and EEG data from corresponding (same) trials results in the impulse response function: the brain response to a single impulse (right, blue line). Cross-correlation of WN sequences with randomly shuffled (different) trials (surrogate IRF) serves as a null hypothesis (EEG and WN are uncorrelated) and produces an almost flat line (right, purple line). C, IRF duration was estimated by extracting beta band power envelope from TFR of IRF signal. B, IRF duration was defined as FWHM between beta power peak at fmax and intersects x1 and x2 between the envelope and 95% quantile of the surrogate IRF power distribution at the same frequency.
SSSEPs
Frequency-specific SSSEP responses were isolated using a spatiotemporal source separation method (Cohen and Gulbinaite, 2017; Cohen, 2022), which is based on generalized eigenvalue decomposition and allows to maximize signal-to-noise ratio of steady-state responses by exploiting information present in interchannel covariance matrices. Importantly, such approach also allowed us to account for individual differences in somatotopic organization and to separate narrow-band tactile stimulation-related activity from temporally co-occurring broadband movement artifacts in beta band. Thus, instead of analyzing SSSEPs from a subset of electrodes with maximum power at the stimulation frequency, we analyzed a linearly weighted combination of signal from all electrodes.
For each participant and stimulation frequency, a separate spatial filter was constructed by temporally narrow-bandpass filtering (Gaussian filter) the raw data (X) around the stimulation frequency f (FWHmean = 0.5 Hz) and at the two neighboring frequencies (f ± 1 Hz; FWHmean = 2 Hz). Temporally filtered data (500-6500 ms relative to stimulus onset) were then used to compute covariance matrices: one “signal” matrix (S covariance matrix) and two “reference” matrices that were averaged (R covariance matrix). The first 500 ms contain evoked potentials that have different sources than SSSEPs (Nangini et al., 2006), and thus were excluded to not compromise the quality of the spatial filter (Cohen and Gulbinaite, 2017). Generalized eigenvalue decomposition (MATLAB function eig) performed on “signal” and “reference” covariance matrices returned matrices of eigenvalues and eigenvectors. To increase the robustness of the spatial filters, we applied a 1% shrinkage regularization to the average “reference” covariance matrix. Shrinkage regularization involves adding a percent of the average eigenvalues onto the diagonals of the average “reference” covariance matrix (Cohen, 2022). This reduces the influence of noise on the resulting eigen decomposition. The eigenvectors (column vectors with values representing electrode weights [w]) were used to obtain component time series (eigenvector multiplied by the original unfiltered single-trial time series [wT X]). The component with the highest signal-to-noise ratio in the power spectra at the stimulation frequency was selected for further analysis (of total 324 components = 18 subjects × 18 stimulation frequencies, in 264 of cases the first component had the highest SNR, in 48 cases the second, in 11 cases the third, and only 1 case the fourth component). Topographical representation of each component was obtained by left-multiplying the eigenvector by the signal covariance matrix (wTS). The obtained topographical maps were normalized, and the sign of eigenvector was flipped for 5 subjects that showed spatial peaks opposite to that of the group average. The sign of the components affects only the representation of the topographical maps and has no effect on component time series (Cohen, 2022).
Power at each stimulation frequency was computed using FFT on single-trial component time series in the 500-6500 ms time window (relative to the stimulus onset) and zero-padded to obtain power exactly at the stimulus frequency with 0.1 Hz resolution. The absolute value of FFT coefficients was squared and averaged across trials. To facilitate comparison across SSSEPs elicited by different stimulation frequencies, SSSEP power values were expressed in SNR units:
For each participant, we plotted SSSEP power and eigenvalues as a function of frequency (frequency tuning curves) and selected the frequency with the highest power as resonance frequency for that participant. In cases when frequency tuning curves contained double peaks, the frequency with the highest eigenvalue was selected. The eigenvalue indicates how well the associated eigenvector separates the “signal” and “reference” matrices and thus can be used as additional selection criterion when comparing components. The statistical significance of SSSEPs was evaluated by repeating the above procedure (constructing spatial filters, obtaining component time series, computing power spectra) using trials in which neither the stimulation frequency nor its harmonics were present. Eigenvector associated with the highest eigenvalue was selected to construct spatial filters expected when no SSSEPs are present. This procedure allowed to determine the SNR values that can be expected when no stimulation at a particular frequency was delivered and to test whether the stimulation protocol induced any significant SSSEP response.
IRF and ERP time-frequency analysis
Time-frequency decomposition for IRFs and evoked potentials in response to single-pulse stimulation were performed via continuous wavelet transformation. The total power of the ERP was calculated by multiplying the power spectrum of single-trial EEG data with the power spectrum of complex Morlet wavelets (
To extract the non–phase-locked part of the EEG response (Fig. 4G), we first averaged the time-domain signal and subtracted it from all trials before time-frequency decomposition. The phase-locked part of the signal (Fig. 4D) was calculated by subtracting non–phase-locked power from the total power. Decibel baseline correction was applied to ERP TFRs according to the following formula
Given some variability for maximum amplitude IRF channel across subjects (individual differences in anatomy and thus projection on the scalp), subject-specific channels of interest were defined as follows: The TFR of IRFs was averaged over all channels, and the time-frequency window of interest (0-200 ms, 15-40 Hz), showing increased beta band power, was defined based on visual inspection (Fig. 3B,E). The channel with the largest power value within this window was selected as the channel of interest for this subject (Fig. 3C,F). The frequency at which the largest power was observed was taken as the individual's IRF peak frequency.
Duration of IRF at the channel of interest was quantified in the following steps. First, we extracted the power envelope of IRF at individuals' IRF peak frequency (fmax) determined from the TFR of IRF (Fig. 2C). Then, the IRF duration was defined as the FWHM between the maximum IRF power and intercepts x1 and x2, which were defined as the time points where IRF power envelope exceeded the 95% quantile of IRF surrogate power values (null hypothesis distribution).
TFRs of IRFs were statistically assessed by comparing each individual observed time-frequency power value to the corresponding distribution of power values from the surrogate distribution (Fig. 4A). If the observed power exceeded 95% quantile of the surrogate distribution, it was considered significant. TFRs of phase-locked and non–phase-locked EEG responses were assessed by comparing all participants time-frequency power values (N = 18) to the baseline power values at the corresponding frequency using t tests. Correction for multiple comparisons was done using false discovery rate (FDR). Before comparison, baseline power values were averaged between −400 and −100 ms for each individual frequency. We only considered the largest cluster of significant time-frequency points for further analysis.
Spontaneous β burst analysis
Spontaneous β bursts were extracted from the −500 to 0 ms baseline activity from all 3 paradigms (WN stimulation, SSSEP, and single-pulse stimulation). Single-trial TFRs were calculated at subject-specific IRF channels of interest as described in section IRF and ERP time-frequency analysis. Power thresholds for β burst detection were set to be 3 times the SD of all beta band (15-40 Hz) power values calculated over all trials. Before calculating the SD, we removed 1% of the largest power values to reduce noise stemming from high-power artifacts. We used MATLAB's imregionalmax function to detect peaks (between 15 and 40 Hz) in the TFRs that exceeded the power threshold. From the bursts that matched our selection criteria, we calculated individual frequency histograms. Individual peak frequencies were defined as the highest peak in the frequency histograms.
Results
The power of ongoing (non–phase-locked) β oscillations is generally decreased in response to tactile stimulation and motor actions (Salmelin and Hari, 1994 ; Hari et al., 1997; Salenius et al., 1997; Neuper and Pfurtscheller, 2001; Parkkonen et al., 2015). More recent work capitalizing on single-trial analysis has demonstrated that spontaneous β bursts occurring in close temporal proximity to tactile targets deteriorate stimulus detection (Shin et al., 2017; Lundqvist et al., 2018; Little et al., 2019). These and other findings have led to the hypothesis that somatosensory β oscillations serve as an inhibitory brain rhythm that is generally absent during periods of increased processing demands for incoming information (Karvat et al., 2021). Given the aforementioned reports on burst-like nature of β oscillations and previous findings of the impulse response function in the visual domain (VanRullen and Macdonald, 2012), we set out to test whether the somatosensory system also responds to WN vibrotactile stimulation with a rhythmic IRF, a transient increase in beta band power. We first characterized the frequency and duration of the IRF in response to WN, and then we tested whether the properties of the IRF are related to: (1) the resonance properties of the somatosensory system measured in a SSSEP paradigm; (2) dynamic changes in beta power in response to single-pulse stimulation; and (3) spontaneous β bursts occurring during the intertrial interval (a proxy for the resting state beta band activity).
Random noise stimulation reveals oscillatory impulse response function
Cross-correlation of WN tactile stimulation sequences with concomitantly recorded EEG signal revealed an oscillatory response present for lags up to ∼200 ms, the IRF of the somatosensory system (Fig. 3). Topographically, the IRF amplitude was maximal over the left parietal and frontal channels. The TFR of IRF showed that these fluctuations in response to tactile impulses were constrained to the beta band. Spatial and frequency characteristics of the somatosensory IRF were similar to previous reports on vibrotactile stimulation, with the highest beta power observed contralaterally to the stimulated hand (Nangini et al., 2006; Bardouille and Ross, 2008; Bardouille et al., 2010; Spitzer et al., 2010; Langdon et al., 2011; Colon et al., 2012; Vlaar et al., 2015), thus indicating a genuine somatosensory response (Fig. 4B).
IRFs and their corresponding TFRs from 2 representative subjects. A, The observed IRFs (blue lines) were calculated by cross-correlating tactile WN sequences with the recorded EEG signal. Red lines indicate cross-correlation result after randomly shuffling trials (null hypothesis: EEG and WN are uncorrelated). B, TFR of the single subject IRF at the individual channel of interest (CP3). C, Topography of average IRF beta power (15-40 Hz, 0-200 ms). D-F, Same as in A-C, but for a different subject (channel CP5).
A, Grand average TFR of IRF (averaged over individual channels of interest). Black contours represent contiguous time-frequency points at which the observed IRF was significantly different from surrogate IRFs at p < 0.05 (corrected for multiple comparisons across time-frequency points using FDR). B, Grand-average power spectrum of the IRF calculated in the 0-200 ms time window. Individual subjects showed spectral peaks distributed around a mean of 25.1 Hz. C, Topography of IRF beta band power (15-40 Hz, 0-200 ms). Black dots indicate the distribution of electrodes that were identified as a subject-specific channel of interest. D, TFR (averaged over individual channels of interest) and β topography of the phase-locked EEG response (ERP). E, Grand-average power-spectra calculated in the 0-400 ms time window over individual channels of interest for the phase-locked EEG response. G, Grand average TFR of non–phase-locked EEG response. Black contours represent regions in which contiguous time-frequency points were significantly different from the prestimulus baseline at p < 0.05 (corrected for multiple comparisons using FDR). Purple contours represent a region, in which the real IRF was significant (the same as black contour in A). H, Same as in E, but for non–phase-locked EEG response. Individual subjects showed local minima at frequencies distributed around a mean of 23.1 Hz. F, I, Topography of beta band power (15-40 Hz, 0-400 ms) for the phase-locked and non–phase-locked part of the EEG response.
To test whether the β oscillatory response to WN tactile stimulation was merely a result of a general increase in beta power in response to tactile stimulation, we repeated the same analysis steps; however, this time we cross-correlated WN sequences with random EEG trials 10,000 times and reperformed spectral analyses. We statistically compared the observed average power values at every time-frequency point to the distribution of power values from the surrogate distribution. Real IRF power values were significantly higher for lags between 0 and ∼200 ms and frequencies between 10 and 60 Hz. The power spectra calculated by averaging over all time points between 0 and 200 ms showed a distinct peak at 20.1 Hz, which was not present in the surrogate distribution (Fig. 4B). Together, this initial assessment suggests that the somatosensory system responds to tactile stimulation with an oscillatory β burst that is: (1) phase-locked to tactile stimulation and (2) specific to the exact WN sequence that was presented on that particular trial.
Next, we quantified individual participants' IRF peak frequency to estimate the length of individual IRFs and to compare its spectral properties to other β signatures. Analyzing the maxima in participants TFRs at individual channels of interest allowed us to identify individual peak frequencies ranging from 18 to 39 Hz (mean = 25.1, SD = 6.12). These frequencies are well within the beta band, although IRF peak frequency of 2 subjects (39 Hz) neared the upper edge of classically defined beta band.
Duration of IRFs
Visual IRFs can be observed for up to 1 s, and last on average 4 or 5 α cycles (Brüers, 2017). We were curious if the short (∼200 ms) duration of tactile IRFs compared with the visual ones was a result of oscillatory frequency difference between tactile and visual IRFs: Albeit temporally shorter, tactile IRFs might oscillate for the same number of cycles as their visual counterpart. We therefore quantified the length of each subject's tactile IRFs (at the channel of interest) based on the power envelope at individual peak frequency (Fig. 2B). The average duration (FWHM relative to 95% quantile surrogate distribution) of individual tactile IRFs was 116.58 ms (SD = 19.13). We used the following formula to calculate the number of cycles for a given individual peak frequency
Non–phase-locked beta power and IRF share spectral properties
To directly compare the spectral properties of IRFs with beta band activity in response to traditional sparse stimulation protocols, we performed time-frequency analyses of ERP data. We separately analyzed non–phase-locked and phase-locked activity (for details, see Materials and Methods). Similar to previous reports (Salmelin and Hari, 1994; Hari et al., 1997; Salenius et al., 1997; Neuper and Pfurtscheller, 2001; Parkkonen et al., 2015), non–phase-locked beta power decreased significantly as a result of tactile stimulation (Fig. 4G). This decrease was most pronounced over the left parietal/frontal channels (Fig. 4I). For each participant, we identified the peak frequency at which the strongest decrease in beta power occurred by detecting power troughs following the stimulus onset (Fig. 4H). Individual minima occurred at frequencies ranging from 19 to 33 Hz (mean = 23.1 Hz, SD = 3.34 Hz). To test whether the non–phase-locked β signatures and IRFs might share a common neural generator, we correlated individual β trough frequencies. Spearman correlation between IRF and decreases in non–phase-locked beta power revealed a significant relationship (r = 0.55, p = 0.0216; see Fig. 7), indicating that both β signatures might rely on a common oscillatory generator.
Phase-locked beta power and IRF share spectrotemporal properties
The phase-locked part of the EEG response showed a broad band increase in power over somatosensory areas (Fig. 4D,F). Statistical analysis of the TFRs revealed a cluster of significantly increased power across a large frequency band (15-60 Hz) between ∼0 and ∼200 ms that was followed by a more frequency-specific increase in the 15-30 Hz range (∼200 to 400 ms; Fig. 4D).
Intriguingly, the TFRs of phased locked ERPs closely resembled the TFRs of the IRFs (Fig. 4A,D). Related, we found beta band topographies to be highly similar (Fig. 4C,F). To test whether IRFs constitute the initial part of the commonly investigated ERPs, we sought to identify spectral peaks in the TFR of phase-locked ERPs and compare them with the IRF peak frequency. We could identify spectral peaks in all subjects (mean = 16.28 Hz, SD = 5.04); however, for 4 of these subjects, the peaks did not fall within the significant cluster (Fig. 4D) but were significantly lower (6–12 Hz). Furthermore, the peak frequencies did not correlate with IRF peak frequencies (r = −0.21, p = 0.4, Spearman). This null finding reflects the fact that TFR of ERP is dominated by the low frequencies (Figs. 4D, 5E), and frequency-based analyses are suboptimal to uncover short-lived high-frequency spectral components. We thus performed time-frequency decomposition of ERPs and directly compared TFRs of IRFs and ERPs. Our rationale was as follows: If phase-locked ERPs contain IRFs and if IRFs vary between participants, then within-participant ERP and IRF correlation should be higher compared with between-participant ERP and IRF correlation.
A, Within- and between-subject correlation of IRF and phase-locked EEG response (ERP) TFRs. Pixel values within the black outlined area were correlated using Pearson correlation. B, Resulting correlation matrix. C, Individual within-subject correlation values (light gray circles) and between-subject correlations (dark gray bars; error bars indicate 95% CIs). D, t test between within-subject and average between-subject correlation values (error bars indicate 95% CIs). E, Characterization or ERP in time and frequency domain (group average, channel CP5).
All time-frequency values of an individual IRF cluster (Figs. 4A, 5A) were correlated (Pearson's r) with the corresponding time-frequency values of the ERP (Figs. 4D, 5A) within and between subjects yielding a correlation matrix, in which diagonal compared with off-diagonal correlation values should be significantly higher (Fig. 5B). We averaged all row and column vectors corresponding to individual diagonal entries and compared the resulting values to the diagonal using a t test (Fig. 5C,D). TFRs of individuals IRF were significantly higher correlated with the same individuals ERP compared with the ERP of all other individuals (Fig. 5D, t(17) = 2.28, p = 0.035, effect size d = 0.437). This indicates that some spectral components of the phase-locked ERP might indeed be oscillatory IRFs, embedded in the ERPs low-frequency components.
Somatosensory resonance frequency matches the frequency of somatosensory IRFs
From the linear dynamical systems perspective, the frequency response of a system, the system's response to a sinusoidal input with variable frequency, should be similar to the frequency transform of its impulse response function. To test this hypothesis and determine subject-specific resonance frequencies of the somatosensory cortex, we used sine-wave stimulation in 12–39 Hz frequency range (SSSEP paradigm). Steady-state evoked potentials were extracted using a spatiotemporal source separation method that allows to linearly combine information from all the channels instead of analyzing SSSEPs from channels with maximum power at the stimulation frequency (Cohen and Gulbinaite, 2017).
As expected, vibrotactile right index finger stimulation elicited maximal response over the contralateral somatosensory cortex (Fig. 6A). Although for all participants SSSEPs were strongest over the frontocentral left-hemifield channels, substantial across-subject variability was evident. To accommodate for individual differences in SSSEP scalp projections, which are likely caused by anatomic differences, we linearly combined information from all the channels by computing spatial filters (channel weights) that optimally isolate frequency-specific SSSEPs (for details, see Materials and Methods). This was done separately for each participant and each stimulation frequency.
Steady-state somatosensory potential analysis. A, Electrode-level analysis results. Each topographical map represents single-subject power (expressed in SNR units) averaged across all stimulation frequencies (12-39 Hz, 18 frequencies). Note the individual variability in scalp projections of the SSSEP responses. B, Topographical representation of single-subject spatial filter corresponding to each subject's resonance SSSEP frequency. Note the high similarity across subjects and dipolar appearance of the components. C, Subject-average normalized spatial filter weights depicted in B (left) and the associated putative anatomic generators of SSSEPs. D, Subject-average SSSEP power expressed in SNR units plotted as a function of vibrotactile stimulation frequency (left): Input frequencies are depicted on the y axis and output frequencies on the x axis. Brighter colors represent stronger response, with maximal responses observed, as expected, at the stimulus frequency. Right, SSSEP power at stimulation frequencies (diagonal from the matrix of power spectra on the left), with the resonance peak at 26 Hz. Black line “H0 hypothesis” indicates SSSEP amplitude at each of the tested frequencies on trials when stimulation at that frequency was not delivered and can be expected as a result of overfitting. Black line on the x axis indicates that SSSEPs at all frequencies were statistically significant at p < 0.001 (corrected for multiple comparisons across frequencies using cluster-based permutation testing).
Figure 6B depicts normalized spatial filter weights for the resonance stimulation frequencies averaged across participants. The putative anatomic estimate of SSSEPs determined from the subject-average spatial filter was localized, as expected, over the left somatosensory cortex (Fig. 6B).
The subject-average response power (expressed in SNR units) plotted as a function of stimulus frequency showed the strongest responses at the stimulation frequency (diagonal). Harmonic responses that are sometimes reported for rhythmic vibrotactile stimulation paradigms are not visible here because spatial filters were designed to isolate responses to narrow-band rhythmic stimulus and suppress temporally co-occurring activity at other frequencies. All stimulation frequencies elicited statistically significant SSSEPs (for all stimulation frequencies, p < 0.001, corrected for multiple comparisons across frequencies using cluster-based permutation testing), as indicated by the small SNR at each of the 18 tested frequencies on trials when stimulation frequency was not present (Fig. 6D; for details, see Materials and Methods). The highest amplitude SSSEPs, or resonance frequency, on average was ∼26 Hz (22–39 Hz range, SD = 5.32 Hz). These findings are consistent with previous reports that used a “best” electrode(s) approach and reported a resonance of human somatosensory cortex in the 20-26 Hz range (Snyder, 1992; Tobimatsu et al., 1999; Müller et al., 2001).
We extracted individual peak resonance frequencies from all subjects and compared them with the frequency of IRFs and non–phase-locked beta power decreases following stimulation. Resonant frequencies were correlated significantly with peak frequencies of IRF and stimulus-related beta power decreases (p = 0.0003 and 0.0195, respectively, Spearman correlation; Fig. 7). This result suggests that the frequency of increasing and decreasing beta power signatures during stimulus processing might be determined by the fundamental resonance frequency of the somatosensory system.
Peak frequency correlation analysis. Spearman correlations were calculated between individual peak frequencies of IRFs, non–phase-locked EEG response, and spontaneous bursts measured during baseline. Correction for multiple comparisons was performed using FDR.
Spontaneous β bursts and IRFs share spectral properties
Recent studies focused on beta band oscillations at a single-trial level and revealed that instead of sustained β oscillations present in trial-average TFRs, β oscillations are transient burst-like events (Shin et al., 2017; Lundqvist et al., 2018; Little et al., 2019). The duration (3 cycles), frequency range (15-30 Hz), and location (somatosensory cortex) of the reported transients is highly similar to the spectro-temporal properties of IRFs that we report here (mean = 2.8, frequency range 18-39 Hz). To investigate whether the two share neural generators, we quantified average β burst frequency in our dataset and compared them with individual IRF peak frequencies. Beta burst analysis was performed in the prestimulus window (−500 to 0 ms) from all trials (SSSEP, WN, and single-pulse ERP trials). Bursts were detected from single-trial TFRs by using a threshold of 3 SDs (for details, see Materials and Methods). We could identify peaks in all subjects in the 17–30 Hz range (mean = 19.44, SD = 3.06). Individual burst peak frequencies were significantly correlated with peak frequencies of SSSEP resonance (p = 0.038), IRF (p = 0.014), and non–phase-locked beta power decreases (p = 0.0003) (Spearman correlation, Fig. 7). Our analysis provides evidence that several well-investigated β phenomena are all strongly correlated, and potentially reflect activity of the same neural generators, linked to the fundamental resonant properties of the somatosensory system.
Discussion
The visual system's response to an impulse of light is oscillatory (“perceptual echo”) and strongly correlates with the endogenous alpha rhythm (VanRullen and Macdonald, 2012), a feature revealed by cross-correlating stimulus WN sequences with the corresponding EEG trials. These findings suggested that response to visual stimulation is phase-locked to the stimulus and may play an active role in processing stimulus-specific information. Here we extend these findings to the somatosensory domain by showing that the somatosensory cortex responds to tactile stimulation with a phase-locked 3-cycle-long reverberation in the beta band akin to a burst. These stimulus-locked β bursts were significantly correlated in peak frequency with the non–phase-locked decrease in beta power following tactile stimulation, the individual participants' tactile resonance frequency (determined in SSSEP paradigm), and the frequency at which spontaneous β bursts occur during rest. Our findings suggest a common oscillatory generator underlying all β signatures characterized in this study, and reveal a supramodal mechanism in sensory processing.
Phase-locked β bursts in response to tactile stimulation
The cross-correlation between tactile WN sequences and the corresponding EEG signal revealed oscillatory signatures that showed strong power in the beta band (∼22 Hz) over somatosensory cortex. Theoretically, the IRF revealed using cross-correlation should resemble a classical time-domain averaged ERP (Polge and Mitchell, 1970). Our findings, however, revealed significant differences. The spectral content of the phase-locked part of the ERP was dominated by frequencies in the 1–15 Hz range (Fig. 4D). The IRF, on the other hand, showed very little low-frequency components and contained a burst-like increase in beta power, that lasted for ∼3 cycles.
The most surprising finding, in the context of previous research, is the phase-locked nature of somatosensory echoes. In contrast, previous reports have highlighted spontaneous occurrence of beta band bursts in the empty intervals before target detection (Shin et al., 2017), working memory delay intervals (Lundqvist et al., 2018), or movement preparation and response errors (Little et al., 2019). This precise time-locking strongly suggests that somatosensory echoes are directly involved in the processing of tactile stimuli and do not likely reflect secondary processes related to movement planning or working memory maintenance. Which functional role could such an oscillatory signature play? Communication through coherence assumes precise phase alignment to be the key factor in information transfer (Bastos et al., 2015; Fries, 2015). While high-frequency oscillations, such as γ, usually increase in power following stimulation, low-frequency rhythms, such as α or β, decrease in power in these critical moments of information transfer (Fig. 3G). Our finding of phase-aligned low-frequency β oscillations during this critical early stage of processing indicates that β might be fundamentally involved in bottom-up transfer of tactile information.
The relationship between ERPs and IRFs
Tactile stimulation of the skin or electrical stimulation of the median nerve is generally accompanied by two EEG signatures in the somatosensory system: One is a broad-band, phase-locked increase in power, and the other is a non–phase-locked initial suppression of μ and beta rhythm followed by an increase (so-called “β rebound”) (Hari et al., 1997; Neuper and Pfurtscheller, 2001; Parkkonen et al., 2015). It has been a matter of debate whether the initial phase-locked broad-band EEG response (ERP) is generated, at least partially, by phase reset of ongoing oscillations (Sauseng et al., 2007). Somatosensory evoked potentials (SEPs) commonly show two early components (N20 and N60), the interpeak difference between which is ∼40 ms and variable across individuals (Niedermeyer and da Silva, 2005). The interpeak timing could potentially correspond to a 25 Hz oscillation and is similar to average somatosensory IRF frequency, as well as the peak of somatosensory frequency tuning curves reported here. High correspondence in spectro-temporal characteristics between these different electrophysiological measures of somatosensory system (SEPs, somatosensory IRFs, peak of somatosensory tunning curves) indicates that they all capture resonance characteristics of the underlying neural network.
While ERPs could also be considered IRFs, there are some important differences between the two: First, ERPs are responses to sparse tactile stimuli whereas IRFs reflect responses to continuous increments and decrements of tactile stimulation intensity. Second, the WN stimulation technique is much more efficient: While in the computation of ERPs, one sample per ERP-time point is collected per trial, the cross-correlation considers all samples to calculate similarity between the signals. Third, short and sparse presentation of stimuli used in experiments strongly differs from natural tactile stimulation and perception. Natural perception is much more continuous and incoming information fluctuates dynamically at multiple frequencies. Thus, IRFs might be a more accurate characterization of neural EEG correlates of sensory processing. Notably, calculating the IRF via cross-correlation assumes a linear relationship between input and output which is likely violated in neural processing (e.g., because of the involvement of multiple mechanoreceptors that bandpass filter the input frequencies) (Handler and Ginty, 2021). The IRF will hence only capture a linear part of brain response and ERP and IRF might therefore reflect different aspects of the same underlying physiological process. Furthermore, although arguably nonlinear, the visual IRF does, to a certain degree, capture functionally relevant components that can predict perceptual outcomes (Brüers and VanRullen, 2018).
Characterizing the resonant properties of the somatosensory system
Experiments using rhythmic stimulation have shown that the somatosensory system responds more strongly to beta band frequency stimulation (Snyder, 1992; Tobimatsu et al., 1999; Müller et al., 2001). Using an SSSEP paradigm in combination with a spatiotemporal source separation method, we investigated the resonance properties of the somatosensory system and compared it with the spectral characteristics of tactile IRFs, and spontaneous β bursts occurring during intertrial intervals. In line with previous reports, we found the mean resonance frequency (highest SSSEP power) to be ∼25 Hz (Snyder, 1992: 26 Hz; Müller et al., 2001: 27 Hz; Tobimatsu et al., 1999: 21 Hz). Statistical analysis revealed a significant positive correlation between individual SSSEP peak frequencies, individual IRF peak frequencies, and spontaneous β burst frequencies. We interpret these consistent correlations as evidence for a common neural generator that is primarily determined by the resonance properties of the somatosensory system.
A remaining question is to which extent peripheral mechanoreceptors could give rise to the observed β oscillations. The time-varying stimuli we used (150 Hz carrier frequency and 12–39 Hz amplitude modulated sine-wave stimulation) likely activate two types of mechanoreceptors responsive to vibrations of the skin: Pacinian corpuscles are most sensitive to high frequency 100–400 Hz vibrations and Meissner's corpuscles, which respond optimally to 40–60 Hz vibrations (Handler and Ginty, 2021). The here reported rhythmic responses (IRF: 25 Hz, resonance: 26 Hz) fall below the optimal (resonance) frequencies of these. The amplitude modulation contained frequencies between 0 and 75 Hz and could therefore potentially drive Meissner's corpuscles. These receptors, however, have phasic responses and only fire on initial contact with an object to the vibrations caused by the movement between object and skin (Abraira and Ginty, 2013). As we discarded the initial 500 ms and the last 1000 ms of each individual epoch, we find it unlikely that the responses of these receptors could have given rise to the specific beta band oscillations observed in the EEG. Last, our stimulation could have activated Merkel's discs as they respond to sustained point pressure. The firing rate of these cells lies in the 60 Hz range and is hence also unlikely to lead to consistent 20 Hz oscillatory activity in the scalp EEG signal (Iggo and Muir, 1969). We conclude that the observed rhythms most likely hallmark active processing and originate from cortical or subcortical regions.
IRFs as multimodal signatures of cortical processing
Beta oscillations have been associated with a multitude of processes, such as motor stiffening (Baker et al., 1997), attention (Buschman and Miller, 2007; Lee et al., 2013), status quo maintenance (Engel and Fries, 2010), sensorimotor integration (Gilbertson et al., 2005; Androulidakis et al., 2006, 2007; Baker, 2007), working memory (Haegens et al., 2017; Spitzer and Haegens, 2017; Lundqvist et al., 2018; Schmidt et al., 2019), time perception (Baumgarten et al., 2015; Wiener et al., 2018), and top-down communication (Buschman and Miller, 2007; Lee et al., 2013; Cannon et al., 2014; Bressler and Richter, 2015). Importantly, these theories propose an active role of β oscillations in temporally organizing incoming information. The stimulus specific nature and strict temporal alignment of beta band somatosensory IRFs, as we demonstrate here, thus could provide a rhythmic temporal structure for these processes.
Two recent studies have implicated the IRFs in regularity learning and predictive coding processes. Repeated presentation of a single visual WN sequence leads to a gradual increase in oscillatory alpha power in the corresponding IRF (Chang et al., 2017). This increase is persistent even when the repetition is interleaved with a novel WN sequence. The authors hypothesized that, in line with the idea that the visual system dynamically encodes visual sequences, the IRF is a signature of this regularity learning mechanism. The authors discussed the possibility that IRFs might reflect a rhythmic updating of predictions about expected stimulation patterns, thus giving rise to 10 Hz oscillations. This line of arguments was later supported by a simple computational predictive coding model capable of generating physiologically plausible IRFs, which propagated as traveling waves along the visual hierarchy (Alamia and VanRullen, 2019). Their results were verified in human EEG studies that showed similar traveling wave dynamics that altered their direction depending on whether visual input was present or not (Pang et al., 2020). Together, these results suggest that visual IRFs might reflect predictive coding processes that continuously generate and update predictions about the incoming information at an alpha rhythm (Friston, 2005; Linares et al., 2009; Clark, 2013; Seth, 2014). Their locus of generation might therefore lie in the interactions between different columns in the visual hierarchy hypothesized to implement predictive coding computations, potentially through canonical microcircuits (Bastos et al., 2012).
The fact that similar oscillatory IRFs are found in two distinct modalities (tactile and visual) might reflect a general process that is fundamental to cortical processing. Previous work has attempted to quantify auditory IRFs and found ERP like responses (Lalor et al., 2009; Lalor and Foxe, 2010). Notably, however, oscillatory IRFs have not been found in the auditory domain, potentially because of architectural differences between visual and auditory cortical hierarchies (İlhan and VanRullen, 2012). It was also suggested that the auditory system could rely on fundamentally different strategies to process sensory inputs as auditory information, in contrast to visual and tactile information, is defined in terms of temporal fluctuations across multiple frequency bands (VanRullen et al., 2014). It might well be that auditory IRFs could be revealed using more naturalistic stimuli, such as speech, as speech contains multifrequency “perceptual units” (in the words of Giraud and Poeppel, 2012) that range from low frequencies (syllabic rates; 1–8 Hz) to high frequencies (phonetic features; 30–40 Hz) (Dikker et al., 2020).
In conclusion, we find that, similarly to the visual domain, the impulse response function of the somatosensory system is oscillatory, with a maximum in the beta band. This means that the response to a single tactile impulse is reverberated in the beta band for ∼3 cycles. These reverberations significantly correlate in frequency with the resonance frequency of the somatosensory system and spontaneously occurring β bursts during rest, pointing to a common oscillatory generator. Our findings provide evidence that the beta rhythm is actively involved in tactile stimulus processing. Furthermore, the existence of rhythmic IRFs in the somatosensory domain as well as the visual domain supports the idea that they reflect a more basic, modality-agnostic signature of sensory processing.
Footnotes
This work was funded by an ERC consolidator grant P-Cycles (ref. 614244) awarded to Rufin VanRullen. Rasa Gulbinaite has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 843379. We thank Barkın İlhan for the MATLAB code used as a basis to generate stimulus sequences and microcontroller-based device used in the İlhan and VanRullen (2012) study, which was connected to vibrotactile electromagnetic solenoid-type actuator and used to deliver stimuli in this study.
The authors declare no competing financial interests.
- Correspondence should be addressed to Samson Chota at samson.chota{at}googlemail.com