Abstract
Nonrapid eye movement (NREM) sleep is characterized by the interaction of multiple oscillations essential for memory consolidation, alongside a dynamic arrhythmic 1/f scale-free background that may also contribute to its functions. Recent spectral parametrization methods, such as fitting oscillation and one-and-over-F and irregular resampling auto-spectral analysis, enable the dissociation of rhythmic and arrhythmic components in the spectral domain; however, they do not resolve these processes in the time domain. Instantaneous measures of frequency, amplitude, and phase–amplitude coupling (PAC) are thus still confounded by fluctuations in arrhythmic activity. This limitation represents a pitfall for studies of NREM sleep relying on instantaneous estimates to investigate oscillatory coupling. To address this limitation, we introduce “Rhythms and Background” (RnB), a novel wavelet-based methodology designed to dynamically denoise time series data of arrhythmic interference. This enables the extraction of purely rhythmic time series, suitable for enhanced time-domain analyses of sleep rhythms. We first validate RnB through simulations, demonstrating it accurately estimates the spectral profiles of individual and multiple oscillations across a range of arrhythmic conditions. We then apply RnB to publicly available intracranial electroencephalogram sleep recordings, showing that it provides an improved spectral and time-domain representation of hallmark NREM rhythms. Finally, we demonstrate that RnB significantly enhances the assessment of PAC between cardinal NREM oscillations, outperforming traditional methods that conflate rhythmic and arrhythmic components. This methodological advance offers a substantial improvement in the analysis of sleep oscillations, providing greater precision in the study of rhythmic activity critical to NREM sleep functions.
Significance Statement
The Rhythms and Background (RnB) algorithm introduces a novel approach to signal processing in electrophysiology by disentangling rhythmic activity from the arrhythmic background at the time series level. RnB denoise brain rhythms from arrhythmic interference in both the time and spectral domains, providing clearer insights into cerebral oscillatory processes. This breakthrough has direct applications in studying brain connectivity and oscillatory dynamics during sleep. Additionally, its application in clinical populations where pathological changes in arrhythmic activity are common, such as neurodevelopmental and neurodegenerative disorders, will help to better understand abnormal oscillatory processes. By improving the accuracy of rhythmic signal analysis, RnB opens new avenues for understanding brain function and dysfunction in research and clinical settings.
Introduction
Brain rhythms from electroencephalograms (EEG) and other electrophysiological techniques reveal the synchronized rhythmic activity from neural ensembles. These rhythms are often identified using the spectral peaks of the power spectrum (Buzsáki, 2004). A noteworthy case is nonrapid eye movement (NREM) sleep, which is characterized by cardinal rhythms such as sleep slow waves (SSWs; 0.5–4 Hz, delta), theta bursts (TBs; 6–10 Hz), sleep spindles (8–16 Hz, sigma), and sharp-wave ripples (100–200 Hz). NREM sleep rhythms are organized in complex wave sequences, with organized phase–amplitude coupling (PAC) across frequencies occurring locally and between remote brain regions (Klinzing et al., 2016; Helfrich et al., 2019). For instance, sleep spindles are organized by the phase of SSWs and occur preferentially during the transition period of SSWs toward their upstate. PAC between delta and sigma rhythms during NREM sleep is increasingly recognized as a crucial component of overnight memory consolidation, and the loss of this phase coupling in aging may contribute to impaired memory consolidation (Helfrich et al., 2018).
New evidence suggests that “arrhythmic” brain activity recorded from the EEG plays complementary roles to brain rhythms for NREM sleep functions (Helfrich et al., 2021). Arrhythmic activity recruits a broad range of frequencies and is usually expressed as a
At the signal level, arrhythmicity manifests as ubiquitous “desynchronized” scale-free fluctuations in electrophysiological recordings visible “in between” transient rhythmic oscillations and on which they are also superimposed (Helfrich et al., 2021). These scale-free brain fluctuations introduce noise into multiple signal processing methods by reducing the reliability of phase and amplitude narrowband estimates associated with brain rhythms (Jurkiewicz et al., 2021; Donoghue et al., 2022). As scale-free power dominates in lower frequencies, slow fluctuations create transient events that are erroneously detected as slow waves (Samaha and Cohen, 2022). Although methods such as FOOOF and IRASA are increasingly used to distinguish rhythmic and arrhythmic activity at the level of the power spectrum (Pani et al., 2022), no method currently addresses the disentanglement of rhythmic processes from the arrhythmic background at the signal level.
Here, we introduce Rhythms and Background (RnB), a wavelet-based method to extract a rhythmic time series in which the arrhythmic contribution is removed. Our method attenuates the
Materials and Methods
From the spectral to the wavelet paradigm
The objective of this work is to extract as time series the rhythmic activity embedded in electrophysiological signals, in which the ubiquitous
Spectral whitening approaches such as IRASA (Wen and Liu, 2016) proposed to estimate a data-driven spectral filter
RnB translates the model above into a wavelet framework, using the “timescale” decomposition of the signal. The upcoming section will give a detailed mathematical description of the algorithm; here, we provide a didactic and integrative overview.
Wavelets are transient “little waves” of a specific scale localized in time (Fig. 1A). Like a “zoom” function, wavelets at small temporal scales capture fast transient fluctuations, whereas wavelets at larger scales capture slower oscillations. The wavelet coefficients are thus defined by correlating a signal with wavelets at different scales and time position. This results in a “timescale” decomposition of a signal, with wavelets coefficients expanded across scales and time (Fig. 1B,C). Both the scale-free (arrhythmic) background and the scale-dependent (rhythmic) component contribute to each wavelet coefficient in the timescale decomposition. The original contribution of RnB is to remove the scale-free contribution from each coefficient, preserving rhythmic activity in the coefficients.
Timescale wavelet framework. A, Example of a fractional spline wavelet
This is achieved in two successive steps: first, the analysis of the original time series, and second, the synthesis of the rhythmic time series. During the analysis of the “original” time series, the scaling exponent (β) will be estimated from the wavelet power across all scales within each epoch (see Fig. 1C,D for NREM sleep). Because of their sensitivity to scaling properties, wavelets are uniquely suited to characterize the scale invariance inherent to electrophysiological recordings (Abry et al., 2019). This scale-free contribution will then be attenuated by removing the influence of β from all wavelet coefficients in the timescale expansion, preserving only scale-dependent rhythmic activity in the “filtered” representation (Fig. 1E,F). This procedure yields a “rhythmic time series,” providing a substrate for further analyses of oscillatory activity (Fig. 1G).
Mathematically, the scale-free contribution is removed from the coefficients via a fractional differentiation operation grounded in the properties of the spline wavelet basis (Unser and Blu, 1999), as detailed in the next section.
“Mother wavelet,” wavelet family and timescale expansion of a signal
A “mother wavelet”
Localization of coefficients. A wavelet coefficient at scale
Low-frequency residue. The expansion also includes a residue [
Thus, wavelet coefficients provide a timescale representation of the signal, whereas the residual captures an unavoidable low-frequency approximation in finite-length recordings.
The RnB procedure
To extract a rhythmic time series, RnB will proceed in three steps using the analysis of an original time series: (1) estimate β per epoch from the log–log slope of wavelet power across scales; (2) attenuate the
Estimation of the scaling exponent β in the wavelet domain
Rather than estimating
β in the Fourier domain, we use the scaling of wavelet coefficient power given by Abry et al. (2019). Let
Formally, this “wavelet-based” scaling exponent is defined as follows:
Attenuation of 1/f activity and denoising of wavelet coefficients
In analogy with the IRASA approach introduced in Equation 3, RnB adjusts wavelet coefficients in the “timescale” domain to account for ubiquitous, scale-free activity. To this end, we will use the fractional spline wavelets introduced by Blu and Unser (2000), whose regularity parameter
α controls the smoothness of the basis and sets the expected scale law of the coefficients
For a signal with aperiodic exponent
To remove the arrhythmic contribution, RnB will rescale the coefficients as
In addition to the previous
The “rhythmic” wavelet coefficients are therefore given by the following expression:
Synthesis of the rhythmic time series
Since the previous processing removed the
The regularity parameter
We set
The resulting timescale representation and rhythmic time series (Fig. 1E,G) shows a marked attenuation of the arrhythmic component and highlight scale-specific rhythmic activity across multiple NREM3 epochs (Fig. 1F). Note that the rhythmic time series
Spectroscopy of NREM sleep recordings
Spectroscopy refers to the systematic analysis of the intrinsic oscillatory content of a signal or observation, providing a basis for unambiguous interpretation. In RnB, the rhythmic signal is obtained from the synthesis of the denoised wavelet coefficients, in which the scale-free component is attenuated. However, because wavelet power is organized by dyadic scales (i.e., discrete frequency intervals), it lacks a fine quasicontinuous frequency resolution when considered as a power spectrum (as illustrated on Fig. 1F). However, the Fourier power spectra from the rhythmic signals will provide this finer description when averaged across epochs (akin to the IRASA method described in Eq. 3; Wen and Liu, 2016), yielding a reliable representation of dominant oscillatory activity.
In context of this paper, we will define the NREM sleep spectroscopy as the joint analysis of (1) the distribution of scale-free exponents
The “NREM sleep spectroscopy.” A, The probability density distribution of the β from processing 300 20 s NREM3 sleep epochs (same epochs as in D and F). B, The Fourier spectra from the rhythmic signal of the same 300 20 s NREM3 sleep epochs after RnB processing. Light gray lines show individual epochs, while bold black is the average across all epochs.
The RnB algorithm: code accessibility
In the following sections, we apply the RnB algorithm to simulations and intracranial NREM sleep recordings, demonstrating its validity through simulations and application using NREM sleep analysis. A guide how to use this algorithm on user-provided data can be found in a publicly available GitHub (https://github.com/michaelcfoti/RnB-Wavelet). The analyzed dataset for the figures presented in the following sections is available in the Extended Data. The RnB algorithm on the GitHub is available as a MATLAB package licensed under an Apache 2.0 license. The module supports MATLAB 2021a or later, with dependencies upon the “Signal Processing” toolbox. The package is openly developed and maintained. The project’s repository includes the codebase, a demo which allows for illustration use of the algorithm on publicly available and user-provided data and the documentation materials. This pipeline allows user-provided time series to extract a rhythmic time series and perform basic spectroscopy. Rhythmic time series are extracted as a matrix that can be used by the user to perform time-resolved analyses in other pipelines.
Two parameters need to be set to use the algorithm on user-provided data (see the wRnBMain script).
The j1 and j2 parameters specify the range in which the scaling exponent
The parameter J refers to the largest scale (lower frequency) upon which the wavelet analysis/synthesis operations are performed to obtain the rhythmic signal.
These parameters can be changed according to the input data properties, and the documentation on GitHub allows the user to understand the impact of such changes. As default parameters used in this work, the scale ran from 1 to 9 to determine the arrhythmic
Simulation analyses of RnB
Simulations were conducted to assess RnB performance in characterizing background and rhythms. A set of eight second epochs were constructed, combining various scale-free backgrounds with transient 4 s oscillations. Considering fractional Brownian motion signals, the scale-free backgrounds were obtained by filtering white-noise time series to obtain a power law
Four-second alpha oscillations were added to these scale-free backgrounds using a classical phenomenological neural mass model (Jansen and Rit, 1995; Pinotsis et al., 2012). This model simulates the dynamics of coupled inhibitory and excitatory neuronal populations interacting with a pyramidal cell. Examples of the temporal and spectral responses of this pure “rhythmic” model are presented in Figure 3. These neural mass responses were superimposed onto pure scale-free backgrounds to construct realistic time series (see Fig. 4A for an example).
Methods: neural mass simulations. A, A 4 s alpha oscillation produced by a standard neural mass model (with the original neurophysiological parameters onto the Figure). B, The Fourier power spectrum of the model response (A). Notice that the alpha peak is accompanied with a lower-frequency content that corresponds to the first seconds of the simulated oscillation, which does not correspond to a “pure” harmonic oscillation.
Methods: RnB estimations of neural masses simulations. A, Example of a simulated time series of 8 s containing a mixture of a scale-free activity (
Two sets of simulations were conducted to assess RnB's accuracy in estimating the maximal rhythmic amplitude of alpha and delta rhythms, as well as the scaling exponents:
The first set of simulations tested the robustness of RnB estimation for an alpha oscillation (10.5 Hz) simulated at 15 arbitrary amplitude levels within three arrhythmic background conditions (β = 1.7, 2.1, and 2.5).
A second set of simulations assessed the ability of RnB to perform accurate estimations in the context of two concurrent oscillations (delta at 3 Hz and alpha at 13 Hz). The delta oscillation varied across 15 incremental arbitrary amplitude levels, and the alpha amplitude remained constant.
For each condition (oscillation amplitude × background β), we first generated 300 pure scale-free epochs with β sampled around the target value (SD = 0.1). Oscillatory activity was then superimposed on 200 randomly chosen epochs per condition, leaving the 100 remaining ones as arrhythmic controls. To mimic a realistic scenario where rhythms may or may not be present, 60 epochs per condition were randomly selected among the total 300 for statistical analyses. RnB estimates of maximal rhythmic amplitude and scaling exponents were then extracted from these analyzed epochs.
Intracranial NREM sleep analyses using RnB
From NREM sleep recordings, we investigated whether rhythmic time series exhibit cardinal NREM sleep rhythms, notably their spectral signature and interfrequency (delta–sigma) coupling. To do so, we first performed spectral analyses of the NREM sleep rhythmic time series produced by RnB. Second, we detected SSWs in the rhythmic time series and characterized how these events locally group and synchronize other rhythms in the theta and sigma frequency bands. Finally, we assessed whether delta–sigma PAC during NREM sleep differed between rhythmic and original time series.
Intracranial recordings
The intracranial recordings used in this work are extracted from a publicly available repository of intracranial EEG (iEEG) recordings in patients with epilepsy (MNI open iEEG atlas: https://mni-open-ieegatlas.research.mcgill.ca). From this multicentric database providing physiological activity in nonlesional tissues of the sleeping brain, we used the recordings of 16 patients. Analyses were performed on NREM2 and NREM3 sleep recordings, including all electrodes (sampling rate, 200 Hz) localized in the precuneus and the anterior cingulate, two key hubs of cortical synchronization for slow waves and sleep spindles (Murphy et al., 2009; Zerouali et al., 2014). All recordings were segmented into 20 s epochs. The duration of recordings in each region and sleep stage is described in Table 1.
Dataset for the intracranial sleep recordings
Spectroscopy of NREM sleep recordings
Each 20 s NREM sleep epoch was processed with the RnB algorithm to estimate the scaling exponent and extract a rhythmic time series. From RnB outputs, we estimated the distribution of scaling exponents and the average power spectral density separately in NREM2 and NREM3 sleep, what will be referred to as “NREM sleep spectroscopy” (Fig. 2). We compared these results to model fits from FOOOF, a standard spectral parameterization method (Donoghue et al., 2020).
For FOOOF analyses, we used the following parameters: an extraction of three peaks in the 0–45 Hz range, with fixed knee aperiodic modes (single slope). Aperiodic knees were initially parametrized in a subset of models but did not yield stable results and were therefore excluded from subsequent analyses.
Slow-wave detection in the NREM sleep rhythmic signals
We adopted a proof-of-concept approach to detect SSWs in the rhythmic time series from NREM3 sleep, to assess how they impacted the delta–sigma coupling during NREM sleep.
Rhythmic SSWs (rSSW) were identified in the rhythmic time series using a custom-made detector using region-adaptive frequency filtering. This filtering served to use the locally expressed oscillatory content in each region, in line with the recent approaches (Donoghue et al., 2022). This detector operates in three steps.
Step 1: Canonical slow-wave candidates and polarity in iEEG
We first used a standard intracranial slow-wave detector on the rhythmic times series, as defined by Ellenrieder et al. (2016). Briefly, candidates slow-wave half-waves were identified in the rhythmic time series after delta filtering (0.5–4 Hz) by detecting zero crossings. Only half-waves with a duration falling within physiologically plausible limits (125–1,000 ms) were retained. Because intracranial bipolar polarity varies according to neuronal geometry and therefore cannot be used to determine the phase of SSWs, we inferred rSSWs phase by comparing γ-band (30–80 Hz) power between successive half-waves; the half-wave exhibiting lower gamma power was labeled as the rSSWs downstate.
Step 2: region-adaptive rhythmic filtering
From the NREM sleep spectroscopy in each region (rhythmic time series), we defined a stage-specific delta bandpass filter centered on the main rhythmic delta peak which stood out locally during NREM3 (see Fig. 6 for relevant results). For the anterior cingulate, this range was 0.5–2.5 Hz; for the precuneus, it was 1–3 Hz. We then reapplied the same detection criteria as in Step 1 (zero-crossing identification, duration limits, and γ-based polarity determination) within this region-adaptive band to isolate events with maximal rhythmic content for each region and stage.
Step 3: conjunction of detected SW
The final rSSW inventory was defined as the conjunction of events detected in both passes. This procedure acted as a region-adaptative filter to lower the number of rSSWs retained from standard Step 1, in accordance with the frequencies observed in the NREM sleep spectroscopy of each region.
Original-signal SSW inventory
We also produced a second SSW inventory from the original NREM sleep recordings, using the standard intracranial detection criteria described in Step 1 above (von Ellenrieder et al., 2016). This set of SSWs detected in the original intracranial signal (without RnB processing) served to compare NREM sleep findings between the rhythmic time series provided by RnB and the original signal.
Windowing and classification of SSWs
For all analyses (spectroscopy, time–frequency, and PAC), we extracted 4 s windows centered on the maximal hyperpolarization peak of each SSW or rSSW.
As in Bouchard et al. (2021), we characterized the transition frequency of each SSWs and rSSWs by measuring the time delay
Spectroscopy analyses of rSSWs
In the first set of analyses, the spectroscopy of rSSWs was compared with an equivalent number of baseline periods, using 4 s NREM3 sleep time windows without rSSWs. Second, we directly compared the spectroscopy of rSlowS with that of rFastS, as SSWs with a slower transition frequency are known to better synchronize NREM sigma (≈13 Hz) rhythms (Chylinski et al., 2022).
Time–frequency analyses of rSSWs
To examine the timing of other NREM rhythmic activity around rSSWs, we used the rhythmic time series to perform time–frequency analyses around two distinct characteristic markers of rSSWs: (1) for standard rSSWs analyses, the peak of hyperpolarization was used as the central landmark consistent with previous studies (Cox et al., 2014); and (2) when specifically comparing rSlowS and rFastS, the window was centered on the zero-crossing point between hyperpolarization and depolarization, representing the midpoint of the transition phase for both types. This alignment ensures equivalence between the two types of switchers, as depolarization occurs more rapidly in fast switchers.
To assess the impact on rSSWs on other rhythms, we computed two types of scalograms, following Tallon-Baudry et al. (1999): (1) the induced power scalogram and (2) the phase-locked power scalogram. Both provide amplitude estimates in a time–frequency plane but differ in their computation. Induced power is obtained by averaging the amplitude of wavelet coefficients computed on a trial-by-trial basis. In contrast, phase-locked power is computed by first averaging the coregistered signal across rSSWs and then extracting the amplitude of the resulting wavelet coefficients (Tallon-Baudry et al., 1999). Phase-locked power thus isolates activity that is systematically aligned in phase across events, whereas induced power reflects the total event-related power regardless of phase consistency.
For all rSSWs analyses, we calculated the time–frequency wavelet transform within a 1 s window around each central landmark (indicated as t = 0 in Fig. 9). The scalograms were calculated using the continuous wavelet transform with the Morse analytic wavelet (default parameters,
PAC analyses of rSSWs
To investigate how the rhythmic and original time series differ in their sensitivity to the PAC of delta and sigma rhythms during NREM sleep, we computed delta–sigma PAC from 4 s epochs centered on SSW epochs in the original time series and around rSSWs epochs in the rhythmic time series (rhythmic PAC, rPAC).
For SSWs, we used the phase
For SSWs or rSSWs, the PAC and rPAC are defined as follows:
Detection of TBs associated with rSSWs in the rhythmic time series
The impact of TBs on delta–sigma PAC was investigated specifically in the precuneus (Jiang et al., 2019) because we identified a theta peak in the precuneus spectroscopy (Fig. 7). This was done because TBs were previously identified as precursors and potential modulators of sleep spindles (sigma) activity (Jiang et al., 2019). To do so, we adapted an existing NREM TB detection algorithm for the rhythmic time series. We used the specific θ-band from the NREM sleep spectroscopy of the precuneus (theta mode
Statistics
We used statistical tests to determine the accuracy of RnB in our simulation and to evaluate the impact of the SSW switcher types on rhythmic power and rPAC. Data distribution was inspected with Shapiro–Wilk tests. Parametric tests (two-tailed) were used when normality held; otherwise, nonparametric equivalents were chosen. All statistical analyses were performed using SPSS 26. All statistical tests, including confidence intervals, are reported in the statistical table (Table 2).
Statistical table
Linear regressions were used to assess the relationship between simulation parameters and RnB estimates (estimation of maximal rhythmic amplitude and arrhythmic scaling exponents). Simulation parameters (arrhythmic beta, rhythmic amplitude) and their interactions were entered as predictors of RnB outcomes in the SPSS Statistics 26 software.
For intracranial data of NREM sleep, bilateral t tests (with 1,000 bootstraps to obtain bias-corrected confidence intervals) or nonparametric Mann–Whitney U tests investigated the impact of the switcher type (rSlowS vs rFastS) on the NREM sleep spectroscopy (arrhythmic beta and rhythmic power). For rhythmic power, we integrated power in adaptive frequency bands for identified peaks in the NREM sleep spectroscopy of the anterior cingulate and precuneus (delta, 0.5–4 Hz in both regions; sigma, 8–14 Hz in the anterior cingulate, 11.5–14.5 Hz in the precuneus; theta, 5–9 Hz in the precuneus; Fig. 6).
Mann–Whitney U tests were used to test the difference in delta–sigma PAC between the original and rhythmic time series separately for both types of switchers in the precuneus and in the anterior cingulate. PAC and rPAC difference between switchers were also investigated separately in each region. Finally, the impact of TB on PAC differences was investigated in the preceding classification (rSSWs + TBs vs rSSWs − TBs) in the precuneus.
Results
We used realistic simulations to evaluate RnB performance in estimating rhythmic (oscillation amplitude) and scaling (β) parameters. As a first illustrative example of RnB application, we considered a simulated alpha (10.5 Hz) oscillation embedded in an arrhythmic background (Fig. 4). The initial signal (Fig. 4A) exhibited a clear arrhythmic drift, which was evident in its Fourier spectrum (Fig. 4B). After RnB processing, the rhythmic signal shows a clear transient oscillation, and the arrhythmic drift is removed (Fig. 4C). The spectral analysis of the rhythmic signal reveals a clear alpha peak in the rhythmic spectral density (Fig. 4D).
Regression analyses then examined how varying the amplitude of simulated alpha oscillations (15 amplitude levels × 3 arrhythmic background conditions; Fig. 5A,B) predicted RnB alpha amplitude estimates. RnB alpha amplitude scaled significantly and positively with the simulated alpha amplitude (b = 0.99; p < 0.0001), accounting for nearly all variance in RnB outputs (R2 = 0.98; Fig. 5A). A small amplitude × background interaction was observed (t = 8.58; p < 0.0001; ΔR2 < 0.01), indicating a weaker relationship in the background condition with higher β values (Δb = −0.05; p < 0.0001). Importantly, the association between simulated and measured alpha amplitude remained strong and significant in all three background conditions (slopes ranging from b = 0.94 to 0.99; all p < 0.0001).
Simulation results: RnB estimation of rhythms and background. A1–A3, Examples of simulated “pure” arrhythmic time series with three different scale-free β exponents (β in each condition, red, 1.69; blue, 2.11; yellow, 2.49); 300 simulations were generated for each arrhythmic condition. Simulated alpha oscillations (10.5 Hz) at 15 fixed amplitude levels (x-axes on the B and C subplots) were added to a subset of simulations in each of the arrhythmic conditions (see Materials and Methods for details). B, Median and kernel density of RnB estimated alpha amplitudes (Y-axis) for the 15 fixed simulated alpha amplitude levels (X-axis), in each arrhythmic condition (red, β = 1.69; blue, β = 2.11; yellow, β = 2.49). C, Median and kernel density of RnB estimated β exponents (Y-axis) for the 15 fixed simulated alpha amplitude levels (X-axis), in each arrhythmic condition (red, β = 1.69; blue, β = 2.11; yellow, β = 2.49) D, Example of two concurrent simulated oscillations (delta at 3 Hz and alpha at 13 Hz). Simulated delta oscillations were varied at 15 incremental amplitude levels (x-axes in F1 and F2 subplots), whereas the delta oscillation was fixed. E, Power spectral density (after RnB processing) of the dual-oscillation simulation model showed in D. F1, Median and kernel density of RnB estimated delta amplitude (y-axis) for the 15 fixed delta amplitude levels (x-axis), concurrent with an alpha oscillation. F2, Median and kernel density of RnB estimated alpha amplitude (y-axis) for the 15 fixed amplitude levels of the simulated delta oscillations (x-axis), concurrent with alpha oscillation. G, Median and kernel density of RnB estimations of the β exponents (y-axis) for the 15 fixed amplitude levels of the simulated delta oscillations (x-axis), concurrent with an alpha oscillation.
We next tested whether varying the amplitude of simulated alpha oscillations predicted RnB scale-free background (β) estimates across the same 15 amplitude levels × 3 background conditions (Fig. 5C). In contrast to the strong effects on alpha amplitude estimates, RnB β estimates scaled only modestly with simulated alpha amplitude (b = 0.10; t = 91.1; p < 0.001), which accounted for <1% of the variance in β estimates (R2 < 0.01; Fig. 5C). A small amplitude × background interaction was again observed (t = 9.47; p < 0.001; ΔR2 < 0.01), indicating a stronger positive association in the high-β background condition (Δb = 0.01; p < 0.001). Across all models, background condition remained by far the dominant predictor of RnB β estimates (t = 663.3; p < 0.0001; R2 = 0.99).
Finally, we evaluated RnB estimates in simulations with two concurrent oscillations (delta at 3 Hz and alpha at 13 Hz; Fig. 5D–G). Simulated delta amplitude varied across 15 incremental levels, while alpha amplitude remained constant. RnB estimated delta amplitude scaled strongly and positively with simulated delta amplitudes (b = 0.87; p < 0.0001; R2 = 0.85; Fig. 5F1). In contrast, RnB estimated alpha amplitude was unaffected by simulated delta amplitude (b = −0.01; p = 0.27; Fig. 5F2). By comparison, β estimates increased significantly with higher simulated delta amplitude, explaining ∼31% of the variance of β estimates (b = 0.03; p < 0.001; R2 = 0.31; Fig. 5G).
Overall, RnB provided reliable estimates of both rhythmic amplitudes and scaling exponents across single and concurrent oscillations (alpha and delta). Two minor caveats should be noted: (1) RnB scaling exponent (β) estimates are slightly overestimated in the presence of high amplitude oscillation (especially delta), and (2) rhythmic amplitude estimates (alpha) are slightly attenuated in high-β arrhythmic conditions. Together, these patterns indicate that while rhythmic and arrhythmic components are disentangled, minimal residual cross-influence persists.
The spectroscopy of sleep recordings
We now demonstrate the potential of RnB to characterize NREM sleep rhythms using intracranial recordings of NREM sleep. Our analyses focused on the anterior cingulate and the precuneus, two key regions involved in NREM delta (slow waves) and sigma (sleep spindles) rhythms (Murphy et al., 2009; Zerouali et al., 2014).
For illustrative purposes, we first show how the NREM sleep power spectrum can be parameterized using FOOOF (Donoghue et al., 2020). As expected, results from FOOOF highlight that intracranial recordings of NREM2 and NREM3 sleep spectrum are characterized by different scaling exponents, with steeper average arrhythmic slope in NREM3 compared with NREM2 in both anterior cingulate and precuneus (β shown on Fig. 6B1,B2). FOOOF also reveals sigma peaks in the average NREM spectrum in both stages and regions, with an additional theta peak in the precuneus during NREM2 and NREM3 sleep (Fig. 6B1,B2).
Results: RnB application to NREM sleep recordings. Results from spectral factorization using FOOOF (top row) and spectroscopy analyses using RnB (bottom row) for NREM sleep recordings in the anterior cingulate (two left columns) and precuneus (two right columns). A1, A2, Standard spectral density per region in NREM2 (blue) and NREM3 (pink) sleep, highlighting a significant arrhythmic 1/f background alongside spectral peaks. B1, B2, Residual periodic power is estimated by FOOOF in NREM2 and NREM3 sleep. Dashed vertical lines indicate spectral peaks. C1, C2, Probability density of scaling exponents (per epoch) as assessed by RnB in NREM2 (red) and NREM3 (blue) sleep. D1, D2, Spectral analysis from the RnB rhythmic time series, highlighting distinctive peaks in each region. Thick horizontal lines in the delta band indicate the frequency bands to be used for the detection of SSWs in the rhythmic time series.
Next, we assessed whether the rhythmic time series extracted by RnB provided a reliable spectral representation of NREM sleep. Spectroscopy analyses were performed on rhythmic time series from NREM2 and NREM3 sleep, considering power spectral density and the distribution of scaling exponents. Like FOOOF, the distributions of scaling exponents β assessed by RnB showed higher scaling exponents in NREM3 as compared with those in NREM2 sleep, for both anterior cingulate and precuneus (β distributions shown respectively on Fig. 6C1,C2). Furthermore, higher rhythmic delta power was observed in NREM3 compared with that in NREM2 sleep in the anterior cingulate and precuneus, although the specific frequency bands varied between the two regions (Fig. 6D1,D2). In the anterior cingulate, a single low-frequency delta peak (1.17 Hz) predominated in NREM3 as compared with NREM2 sleep, whereas two delta peaks (1.95 and 3.71 Hz) were found in the precuneus during NREM3 and NREM2 sleep. A rhythmic theta peak (≈7 Hz) was also found specifically in the precuneus during NREM2 and NREM3 sleep. Additionally, rhythmic sigma peaks were found in both regions, with slightly faster frequencies in the anterior cingulate (12.5 Hz) as compared with the precuneus (13.5 Hz) and lower maximal power in NREM3 compared with NREM2. The spectroscopy of NREM recordings is thus concordant with expected standard NREM delta and sigma rhythms, in addition to theta rhythms in the precuneus.
rSSWs modulate rhythmic and arrhythmic activity
Given that rhythmic time series showed higher delta activity during NREM3 compared with lighter (NREM2) sleep in both recorded regions, we detected SSWs in the rhythmic time series from NREM3 sleep using spectroscopy-adapted frequency criteria (see Materials and Method), which we named “rhythmic sleep slow waves” (rSSWs). We wanted to assess how they modulate other rhythms during NREM3 sleep using time–frequency analyses (time-resolved level) and spectroscopy (spectral level).
Detection rates of SSWs and rSSWs were lower in the rhythmic time series than in the original signal during NREM3 sleep (≈5 rSSW/min vs ≈14 SSW/min), which is expected since arrhythmic fluctuations inflate SSW counts in raw data. rSSW average waveform retained the expected biphasic shape of SSWs, with fast-frequency suppression during the downstate (Fig. 7A1,A2). Interestingly, time–frequency analyses of the rSSW in NREM3 indicate that they were associated with sigma activity in their depolarization period, predominantly in the anterior cingulate (Fig. 7B1,B2).
Detection of rSSWs from the anterior cingulate and precuneus during NREM3 sleep. A1, A2, Average waveforms of rSSWs during NREM3 sleep in the anterior cingulate (top) and precuneus (bottom). B1, B2, Time–frequency analyses (induced power) of the rhythmic time series associated with rSSWs coregistered with respect to the hyperpolarization phase. Red/yellow indicates higher normalized power than blue. C1, C2, Distribution of scaling exponents across epochs with rSSW (continuous lines) and epochs without rSSW (dashed lines). D1, D2, Rhythmic power density with rSSWs (continuous lines) and without rSSWs (dashed lines).
The spectroscopy of rSSWs was compared with baseline NREM3 sleep epochs (without SSWs). The probability density distributions showed lower scaling exponents for epochs with rSSWs than for baseline epochs, particularly in the anterior cingulate (Fig. 7C1,C2). Rhythmic power density showed higher delta and lower sigma rhythmic power for epochs with rSSWs compared with epochs without rSSWs in both regions, although the anterior cingulate showed stronger stage-related difference in rhythmic delta (Fig. 7D1,D2). Rhythmic sigma peak was wider in the anterior cingulate as compared with the precuneus. The theta peak previously observed in the precuneus also showed lower power during rSSWs as compared with the baseline, suggesting that theta activity in the precuneus may not be specifically linked with rSSWs.
The spectroscopy of SSWs differs according to transition frequency
SSWs can be categorized as slow or fast switchers, depending on the duration of their transition frequency, the inverse of their transition period between their hyperpolarization and depolarization phase (Bouchard et al., 2021). Sigma rhythms are preferentially nested in the transition of SSWs with longer transition periods (slow switchers) as compared with fast switchers during NREM sleep (Chylinski et al., 2022). We thus evaluated in the rhythmic time series how slow switchers group sigma activity during NREM3 sleep.
We categorized each rSSW detected in the rhythmic time series from NREM3 sleep as a rhythmic slow switcher (rSlowS) or a rhythmic fast switcher (rFastS; see Materials and Methods). Depending on their transition frequency, we clustered rSSWs between rSlowS and rFastS and found a probability distribution skewed toward rFastS in both regions (Fig. 8A1,A2). The transition frequency separating rSlowS and rFastS was ∼1.32 Hz in the anterior cingulate and 1.67 Hz in the precuneus, indicating a slower transition from hyperpolarization to depolarization in the anterior cingulate (see Fig. 8B1,B2 for examples).
The spectroscopy of rSSWs in function of their transition frequency during NREM3. rSSWs can be classified according to the duration of their transition period between hyperpolarization and depolarization (i.e., transition frequency). Slow waves with a shorter transition period between their hyperpolarization and depolarization peaks are considered “fast switchers” (rFastS in the rhythmic time series), whereas rSSWs with a longer transition periods are considered “slow switchers” (rSlowS in the rhythmic time series). A1, A2, The probability density distributions of the transition frequencies of rSSWs in the anterior cingulate (A1) and precuneus (A2), highlighting a bimodal distribution of transition frequency in both regions. The thresholds to cluster rFastS and rSlowS are shown as dashed vertical lines in both regions. B1, B2, Average waveforms of rSlowS and rFastS epochs in the rhythmic time series. The dashed lines show the transition midpoint. C1, C2, Probability density distributions of the scale-free exponents for rSlowS and rFastS in each region. D1, D2, Spectral power density of the rhythmic time series associated with rFastS and rSlowS in NREM3 for both regions. The horizontal shaded bars indicate local frequency bands (delta, theta, sigma) used in statistical analyses (t tests). Bilateral significance for each band is shown separately using stars (*p < 0.05; **p < 0.01). Significant differences in rhythmic power between rSlowS and rFastS were found in the anterior cingulate only.
The spectroscopy of rSlowS differed from rFastS particularly in the anterior cingulate. In this region, the scaling exponents of rSlowS showed a distribution with lower median values compared with rFastS (Fig. 8C1; standardized Z, −4.5; p < 0.001 on Mann–Whitney standardized U test). In the precuneus, no difference between the distributions of scaling exponents was found between slow and rFastS (Fig. 7C2; p > 0.05). In the anterior cingulate (Fig. 8D1), rSlowS showed significantly higher rhythmic delta power (t test, p < 0.001) and higher rhythmic sigma power (t test, p < 0.05) as compared with rFastS. In the precuneus, rhythmic power spectral density did not significantly differ between rFastS and rSlowS in either considered frequencies (all t tests, p > 0.05; Fig. 8D2).
RSlowS and rFastS regroup differently transient sigma and theta rhythms
To evaluate how other rhythms were regrouped in time by the transition of rSSWs, we examined induced and phase-locked power in time–frequency analyses from the rhythmic time series. Induced power reflects the overall amplitude from the average of the wavelet coefficients across all events. In contrast, phase-locked power is the wavelet amplitude computed on the average signal coregistered to a common landmark. Phase-locked power thus provides an indication of how consistently a given event or rhythm occurs in phase at a specific point across coregistered event (here, relative to the transition midpoint of rSSWs), whereas induced power captures increase or decreases in oscillatory energy regardless of the precise phase alignment of the signal. Previous work shows that sigma rhythms tend to occur specifically in phase during the transition from hyperpolarization to depolarization of SSWs and especially in waves with a slower transition such as rSlowS (Chylinski et al., 2022). We thus hypothesized that rSlowS would exhibit stronger phase-locked sigma power during their transition when compared with rFastS.
Induced sigma power was observed around the transition of both rSlowS and rFastS, in the anterior cingulate (Fig. 9A1–C1) and in the precuneus (Fig. 9A2–C2). As expected, rSlowS in the anterior cingulate further showed noticeable phase-locked sigma power during their transition, reaching maximal power around the depolarizing peak (Fig. 9B1). rFastS in the anterior cingulate also exhibited some phase-locked power around the depolarizing peak (Fig. 9D1), although less pronounced than for rSlowS. In the precuneus, phase-locked sigma power was not strongly modulated around the maximal depolarization peak, for either rSlowS or rFastS. However, some phase-locked sigma emerged ∼400 ms later after the depolarization peak of rFastS in the precuneus (Fig. 9D2).
Time–frequency analyses of rSSWs in function of their transition frequency during NREM3 sleep. Time–frequency maps (scalograms) of rSSWs in the anterior cingulate (top) and precuneus (bottom), coregistered to the transition midpoint. Slow switchers (longer transitions; left panels A, B) are contrasted with fast switchers (shorter transitions; right panels C, D). Induced power is shown in A and C, phase-locked power is shown in B and D. Vertical dashed lines mark the transition (bold) and subsequent depolarization peak (light dash). White arrows highlight transition-locked sigma in the anterior cingulate (B1, D1) and precuneus (D2) and transition-induced theta in the precuneus (A2, C2). Each panel is individually normalized for visualization (red, high power; blue, low power).
In the precuneus specifically, induced theta power was observed with the transition of rSlowS and rFastS (Fig. 9A2,C2). However, theta power was not sufficiently synchronized across events to yield a robust phase-locked power signature around the transition of either rSlowS (Fig. 9B2) or rFastS (Fig. 9D2).
The rhythmic time series highlights coupled delta and sigma rhythms in NREM3 sleep
We compared delta–sigma coupling from the rhythmic time series versus the original time series. PAC was assessed from SSWs (original iEEG signal) and from rSSWs (rhythmic time series), referred to as rhythmic PAC (rPAC). We hypothesized that rSlowS and rFastS events would show higher PAC than in the original signal, due to minimized arrhythmic interferences in the rhythmic time series.
As shown in Figure 10, A and B, rPAC values were significantly higher than PAC values for both types of switchers in the two brain regions (Mann–Whitney U tests, p < 0.001 for all events). Moreover, rSlowS showed higher rPAC as compared with rFastS in the anterior cingulate during NREM3 (Fig. 10A; U = 2.82; p = 0.005). This difference between switchers was not significant with standard PAC in the anterior cingulate (U = 1.85; p = 0.06). In the precuneus, no significant differences in rPAC were found between switchers (Fig. 10B; U = 0.18; p > 0.1).
Delta–sigma PAC in the original and rhythmic time series. Average values and standard error of the mean for rPAC (rhythmic) and PAC (raw), for slow and fast switchers during NREM3 sleep in the anterior cingulate (A) and in the precuneus (B). In the precuneus, rPAC values are shown in C for all rSSWs without TBs (rSSW − TBs) and for rSSWs codetected with a TBs (rSSW + TBs). For each event type, arrows indicate significant differences between PAC and rPACs values. Stars indicate significance value (*: p ≤ 0.01).
Since induced theta power was observed around the transition of rSlowS and rFastS in the precuneus (Fig. 9A2,C2), we also explored the impact of TBs events on delta–sigma rPAC during rSSWs in the precuneus. This was prompted by previous literature suggesting that NREM TBs precede slow-wave-spindle complexes and facilitate delta–sigma PAC (Gonzalez et al., 2018). We identified rSSWs co-occurring with TBs (<500 ms; rSSWs + TBs events) and rSSWs that did not (rSSWs − TBs). We hypothesized higher delta–sigma rPAC for rSSWs + TBs events as compared with rSSWs − TB in the precuneus.
As shown in Figure 10C, delta–sigma rPAC was higher for rSSWs + TBs compared with that for rSSWs − TBs (U = 2.46; p = 0.01). Importantly, rSSW + TB events represented ∼6% of all rSSWs in NREM3 (116/1,815 rSSWs). Although rSSWs + TBs are rare, they constitute a significant modulator of rPAC in NREM3.
In summary, delta–sigma PAC was consistently higher in the rhythmic compared with the original time series. Importantly, RnB revealed that delta–sigma PAC is differentially modulated by specific sleep events: rSSWs increased rPAC in the anterior cingulate, whereas TBs enhanced delta–sigma PAC in the precuneus.
Discussion
The RnB algorithm introduces a novel approach to signal processing in electrophysiology, by disentangling rhythmic activity from the arrhythmic background at the time series level. Spectral methods such as FOOOF and IRASA can separate rhythmic peaks from the scale-free background in the power spectrum (Pani et al., 2022), but they remain limited for time-resolved analyses and signal processing: instantaneous phase, amplitude, and PAC estimates remains biased when the background deviates from a simple power law. To address this limitation, we leveraged the wavelet time–scale paradigm to attenuate scale-free activity at the level of wavelet coefficients, producing rhythmic time series with preserved oscillatory components. Simulations and intracranial NREM sleep recordings confirmed that RnB reliably recovers transient oscillations across multiple frequency bands. Importantly, rhythmic time series extracted by RnB revealed stronger delta–sigma PAC during NREM3 sleep, underscoring its value for studying cross-frequency interactions.
Dynamic of the scale invariance and NREM sleep EEG spectroscopy
Growing evidence indicates that the scale-free background is neither constant nor unitary: it changes across sleep stages and often exhibits frequency-specific “knees,” where the spectrum bends away from a single power law (Gao, 2016; Gao et al., 2020; Brake et al., 2024). In NREM sleep, these knees overlap with the delta range, “hiding” slow oscillatory peaks (Bódizs et al., 2021; Gerster et al., 2022; Schneider et al., 2022). The wavelet-based RnB framework is particularly robust to this problem for several reasons.
First, wavelet-based estimates provide more reliable characterization of the scale-free regime and of low power than Fourier-based approaches (Abry et al., 2019). This advantage arises from key wavelet properties: (1) their vanishing moments, which cancel slow signal drifts or polynomial trends according to their regularity parameter
Simulations confirm that RnB preserves delta and sigma amplitudes across a range of 1/f exponents (Fig. 5). While RnB provides robust estimates of both oscillatory and arrhythmic activity, our simulations revealed residual cross-component influences. These were modest (in general, R2 < 0.01) but became more substantial in the delta range: strong delta oscillations increased β estimates, explaining up to ∼30% variance in RnB β estimates. Given the importance of delta activity in sleep research, this coupling must be acknowledged as a potential confound. At the same time, the systematic nature of the effect makes it tractable, suggesting that future refinements (e.g., calibration, iterative separation) could mitigate this bias and enhance the reliability of RnB in high-delta contexts.
A case study: detection of rhythmic SW and switchers during NREM sleep
Our intracranial data further showed that the rhythmic time series provided by RnB recovered the frontal dominance of rhythmic delta in NREM3, consistent with established slow-wave topography (Steriade, 2006). Using spectroscopy-adapted criteria (Muehlroth and Werkle-Bergner, 2020; Donoghue et al., 2022), RnB enabled detection of SSWs in the rhythmic time series (rSSWs). Importantly, rSSWs in anterior regions modulated sigma power more strongly than those in the precuneus, corroborating results from EEG studies indicating that frontal waves serve as more potent triggers for spindles (Cox et al., 2014) and supporting the validity of our results.
Notably, we further observed lower scaling exponents during rSSWs in the anterior cingulate compared with NREM sleep epochs without rSSWs, which is aligned with earlier spectral parametrization findings (Lendner et al., 2020). This result also converges with studies linking lower scaling exponents to higher neuronal excitability (Gao et al., 2017), which is also observed during the depolarization phase of SSWs (Steriade, 2006). Interestingly, although simulations indicated inflation of β estimates under higher delta amplitudes, RnB remains sensitive to transient reductions in β during NREM3 SSWs and rSlowS, reinforcing the validity of our estimates in relevant physiological ranges.
We also reproduced the fast–slow switchers classification of SSWs in the rhythmic signal (Bouchard et al., 2021). In the anterior cingulate, the transition of slow switchers (rSlowS) was associated with higher phase-locked sigma power, and rSlowS showed higher delta–sigma rPAC when compared with rFastS. These results support the idea that SSWs with prolonged depolarization periods, such as rSlowS, promote stronger spindle generation (Botella-Soler et al., 2012; Cox et al., 2018, 2020; Chylinski et al., 2022). In the precuneus, rSSW–sigma rPAC depended more on co-occurring TBs, consistent with their proposed role as an “early” downstate that favors subsequent spindle emergence (Gonzalez et al., 2018; Jiang et al., 2019). In our data, both event types were relatively rare: most of rSSWs were classified as rFastS (Fig. 8), and TBs occurred in ∼6% of all rSSWs (see Results, The rhythmic time series highlights coupled delta and sigma rhythms in NREM3 sleep). The scarcity of these specific events likely contributes to the high variability of delta–sigma PAC across studies (Cox et al., 2019). Collectively, our findings extend MEG observations suggesting that locally synchronized transient sigma rhythms associated with rSlowS and TBs may act as specific drivers of broader spindle synchronization (Zerouali et al., 2014).
Enhancing PAC measures with RnB
Finally, delta–sigma PAC was higher for SSWs detected in RnB's rhythmic signal compared with those detected in the raw signal. This enhancement likely reflects that the rhythmic time series isolate genuine oscillatory activity, thereby increasing the specificity of rSSW detection and PAC estimation, albeit at the cost of sensitivity, resulting in a lower number of detected events in rhythmic compared with original time series. Our trade-off thus favors precision over exhaustivity. Our finding aligns with the notion that scale-free signal fluctuations can introduce spurious SSW detections and lower PAC estimates, thereby increasing PAC variability (Aru et al., 2015; Donoghue et al., 2022). Such variability has repeatedly been noted as a source of small and heterogeneous PAC effect sizes in sleep research (Mander et al., 2017; Jurkiewicz et al., 2021). By filtering out these arrhythmic confounds, RnB enhances the isolation of rhythmic events, improving the robustness of both slow-wave detection and PAC measures. This is particularly valuable for studying rhythm-dependent interactions across hippocampo–thalamo–cortical networks, where arrhythmic activity may differ between regions (Staresina et al., 2023).
Future directions
RnB offers a practical shift in EEG signal processing, offering a wavelet-based framework to extract “true” oscillations from complex
Conclusion
By leveraging the wavelet-based multiresolution “timescale” paradigm, RnB disentangles NREM sleep rhythms from their arrhythmic background, enabling extraction of rhythms from scale-free dominated signals. Our findings show that this approach preserves physiologically meaningful oscillatory dynamics to study interfrequency interactions of neural rhythms. RnB allows to shift the spotlight: in the raw signal, scale-free dynamics dominate the stage, whereas RnB brings the oscillatory actors to the forefront. In the long run, RnB may support the development of a refined detection and classification system of rhythmic sleep events across the lifespan and in clinical populations, where altered scale-free dynamics may obscure oscillatory markers.
Data 1
Download Data 1, ZIP file.
Footnotes
The authors declare no competing financial interests.
This work was supported by Natural Sciences and Engineering Research Council of Canada Discovery Grant (RGPIN-2022-05351, J.M.L.).
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.
References
In this issue
