Multiscale intracranial EEG dynamics across sleep–wake states: toward memory-related processing

Sleep is known to support memory consolidation through a complex interplay of neural dynamics across multiple timescales. Using intracranial EEG (iEEG) recordings from patients undergoing clinical monitoring, we characterize spectral activity, neuronal avalanche dynamics, and temporal correlations across sleep-wake states, with a focus on their spatial distribution and potential functional relevance. We observe increased low-frequency power, larger avalanches, and enhanced long-range temporal correlations—quantified via Detrended Fluctuation Analysis—during N2 and N3 sleep. In contrast, REM sleep and wakefulness show reduced temporal persistence and fewer large-scale cascades, suggesting a shift toward more fragmented and flexible dynamics. These signatures vary across cortical regions, with distinctive patterns emerging in medial temporal and frontal areas—regions implicated in memory processing. Rather than providing direct evidence of consolidation, our results point to a functional neural landscape that may favor both stabilization and reconfiguration of internal representations during sleep. Overall, our findings highlight the utility of iEEG in revealing the multiscale spatio-temporal structure of sleep-related brain dynamics, offering insights into the physiological conditions that support memory-related processing.

1 Introduction

Understanding the neural dynamics that characterize different states of consciousness, such as wakefulness, non-rapid eye movement (NREM) sleep, and rapid eye movement (REM) sleep, is crucial for elucidating how the brain establishes conditions conducive to memory consolidation. Previous work has established hallmark electrophysiological signatures across vigilance states—posterior dominant alpha during restful wakefulness, slow waves and sleep spindles in NREM (especially N2/N3), and REM-specific patterns such as sawtooth waves and regionally localized delta activity (Berry et al., 2012; Troester et al., 2023; Steriade et al., 1993; De Gennaro and Ferrara, 2003; Fernandez and Lüthi, 2020; Frauscher et al., 2020). Yet most studies isolate a single domain or scale; a unified framework that jointly quantifies spectral content, scale-dependent temporal structure, and spatial organization within the same intracranial dataset remains under development.

This work investigates how neural signals across wakefulness and sleep exhibit state-specific dynamics that may create a functional environment favorable for memory-related processes, employing a multiscale and multimodal approach. Using intracranial EEG (iEEG) recordings from the Multicenter iEEG Sleep Atlas (von Ellenrieder et al., 2020), we apply complementary methods—including spectral analysis, neuronal avalanche statistics, and detrended fluctuation analysis (DFA)—to examine how brain activity patterns, as operationally defined by sleep stages, vary across time scales and anatomical regions.

Studies have demonstrated that sleep stages exhibit distinct neural features. For instance, slow oscillations in the delta band during NREM sleep have been associated with synaptic downscaling, potentially facilitating the stabilization of neural circuits involved in declarative memory (Diekelmann and Born, 2010; Rasch and Born, 2013). Sleep spindles, characteristic of NREM sleep, have been implicated in the facilitation of synaptic plasticity, thereby supporting the consolidation of motor skills and procedural memories (Fogel and Smith, 2011). Conversely, REM sleep has been associated with the reactivation and reorganization of memory traces—processes that may contribute to the integration of emotional content and the refinement of neural representations (Walker, 2009).

Beyond traditional spectral methods, neuronal avalanche analysis enables the characterization of spontaneous cascades of activity across time and space. Empirical observations show that avalanche size and duration distributions display approximate scale-invariant behavior across wake, NREM, and REM sleep (Ribeiro et al., 2010; Wilting and Priesemann, 2019). These patterns have been interpreted in various theoretical frameworks, but here we focus on characterizing how avalanche properties vary empirically across states.

In addition, Detrended Fluctuation Analysis (DFA) is introduced here as one of three principal analytical techniques to quantify long-range temporal structure in iEEG signals, providing complementary insights into the scale-dependent dynamics of neural activity.

By integrating frequency-specific spectral signatures, neuronal avalanche statistics, and DFA-derived temporal structure, this study offers a multiscale characterization of brain dynamics across sleep–wake states, with potential implications for sleep-dependent memory processing and neural organization.

2 Materials and methods

2.1 The atlas dataset

In this study, we utilized data from the Multicenter iEEG Sleep Atlas (ATLAS). This international collaboration has compiled a comprehensive open-access dataset of intracranial EEG (iEEG) recordings from patients with drug-resistant focal epilepsy. Importantly, these recordings encompass both wakefulness and various sleep stages, curated to reflect physiologically normal activity (von Ellenrieder et al., 2020).

The iEEG Sleep ATLAS provides an opportunity to examine cortical dynamics during NREM (N2, N3), REM, and waking states across 38 distinct brain regions, including frontal, temporal, parietal, occipital, and insular cortices. Through precise anatomical localization and expert sleep scoring, this dataset enables the exploration of how different cortical areas contribute to state-dependent neural activity patterns associated with sleep processes.

The recordings consist of artifact-free, 60-second epochs selected during resting wakefulness (eyes closed) and each annotated sleep stage. Signals were resampled at 200 Hz and filtered (0.5–80 Hz). Note that these preprocessing steps are part of the ATLAS pipeline; infra-slow activity (< 0.5 Hz) is therefore unavailable. Analyses are reported in 0.7–75 Hz to ensure robust estimates within the 0.5–80 Hz acquisition bandwidth (fs = 200 Hz). The standardized ATLAS format facilitates reproducibility and integration with computational pipelines, supporting analyses of local and distributed activity patterns. Representative raw iEEG examples across wake, N2, N3, and REM are available from the ATLAS reference used here (e.g., Figure 1A in Kalamangalam et al., 2021).

The richness and granularity of this dataset are particularly suited for investigating neural dynamics at multiple spatiotemporal scales. Specifically, it allows for comparisons of regional activity by means of measures such as power spectra, neuronal avalanches, and temporal correlations across sleep–wake states, particularly those implicated in memory consolidation.

Table 1 summarizes the number of channels per brain region and lobe across sleep and waking conditions. While electrode coverage varies across subjects—particularly during sleep—this variability is addressed by conducting the analysis at the population level. This approach allows us to identify region-specific trends and explore how different cortical areas may contribute to the neural dynamics associated with memory processes.

Table 1

Table 1. The 38 regions classified by brain lobes, and the number of channels per region for the non-REM N2, N3, and REM sleep states, and the waking state.

2.2 Spectral analysis

In this section, we provide a brief overview of the spectral analysis technique applied in this work. The primary tool utilized in this section is the power spectral density (PSD), which is defined as the Fourier transform of the autocorrelation function (Semmlow and Griffel, 2014). For discrete signals, such as those found in the ATLAS, denoted as x[n],

$\begin{array}{ll}\text{PSD}\left[m\right]={\displaystyle \sum _{n=1}^N}{r}_{xx}\left[n\right]exp(-\frac{i2\pi mn}{N}),m=0,1,\dots ,N/2,& (1)\end{array}$

in which N is the total number of samples and r_xx[n] is the autocorrelation function for a discrete-time signal. This function indicates the similarity of a signal with a delayed copy of itself as a function of the delay time k. It is defined by

$\begin{array}{ll}{r}_{xx}\left[k\right]=\frac{1}{N}{\displaystyle \sum _{n=1}^N}x\left[n\right]x\left[n-k\right].& (2)\end{array}$

Alternatively, the spectral power PSD[m] can be calculated using the Parseval relation

$\begin{array}{ll}\text{PSD}\left[m\right]=|X(m){|}^2,& (3)\end{array}$

where X[m] is the discrete Fourier transform (Semmlow and Griffel, 2014):

$\begin{array}{ll}X\left[m\right]=\frac{1}{N}{\displaystyle \sum _{n=0}^{N-1}}x\left[n\right]exp(-\frac{i2\pi mn}{N}),m=0,1,\dots ,N-1,& (4)\end{array}$

in which $\frac{2\pi }{N}$ is the spacing between samples.

In general, when dealing with biosignals, we typically have only a sample of a longer signal; therefore, spectral analysis becomes an estimation process. A commonly employed approach is to compute the Power Spectral Density (PSD) over multiple segments of the sample. These segments can then be averaged to generate a spectrum with improved overall characteristics. When the PSD is computed directly using the Fourier Transform (FT) followed by averaging, it is referred to as a periodogram. One of the most widely used methods for calculating the average periodogram is the Welch technique (as seen, for example, in Semmlow and Griffel, 2014). In this approach, overlapping segments are employed, and a shaping window (i.e., a non-rectangular window) is applied to each segment. Traditional periodograms typically average the spectra of overlapping segments, usually with a 50% overlap. Averaging introduces a trade-off between spectral resolution, which decreases with averaging, and statistical reliability. We report results both with and without frequency-wise normalization. In the normalized case, the PSD at each frequency are rescaled by the corresponding across-epoch mean, which emphasizes relative shifts and crossover points across conditions. The unnormalized PSD, in contrast, preserve the absolute power distribution.

To estimate band-specific features, we followed the procedure of Kalamangalam et al. (2020). For each 60 s epoch, the power spectrum was smoothed with a Savitzky–Golay filter to reduce high-frequency noise while preserving local spectral structure (Virtanen et al., 2020). The smoothed spectrum was then modeled as the sum of five Gaussians (delta, theta, alpha, beta, gamma), fitting amplitude, center frequency, and bandwidth (σ) for each component. Fits were accepted only if R² ≥ 0.95. All spectral estimates were restricted to 0.7–75 Hz; bands were defined as delta (0.7–4 Hz), theta (4–8 Hz), alpha (8–12 Hz), beta (12–30 Hz), and gamma (30–75 Hz). This procedure yields robust, reproducible estimates while minimizing contamination by noise or non-oscillatory features.

2.3 Neuronal avalanches

In terms of iEEG activity, a neuronal avalanche is defined as a consecutive series of time slots with a width of Δt that contains at least one active electrode (Beggs and Plenz, 2003; Plenz and Thiagarajan, 2007). Each avalanche is preceded and followed by at least one time slot without activity. The size of an avalanche is determined by the sum of the absolute amplitudes of the iEEG signal at the active electrodes or simply by the number of active electrodes (s). An electrode is considered active when the signal exceeds a certain threshold. In this study, we identify avalanches in the iEEG signals of the Atlas using the number of active electrodes, separated by periods of inactivity, as a measure of avalanche size.

Specifically, for each electrode, the signal was normalized by subtracting its mean and dividing by its standard deviation, resulting in a z-scored signal with zero mean and unit variance. Discrete events were defined as individual time points where the normalized signal crossed below a threshold of −2.0 standard deviations. This thresholding strategy is widely used to detect statistically significant deflections from baseline activity, while minimizing the influence of background fluctuations and noise (Plenz and Thiagarajan, 2007; Beggs and Plenz, 2003).

Binary event trains were constructed from these threshold crossings, and the activity was subsequently binned into consecutive, non-overlapping 10 ms time windows—a commonly used window size in the literature for capturing scale-invariant neuronal dynamics (Beggs and Plenz, 2003). A bin was considered active if at least one electrode registered an event within that window. A neuronal avalanche was then defined as a sequence of contiguous active bins, bounded before and after by at least one empty bin. The avalanche size was computed as the total number of events (across all electrodes) occurring within the sequence.

The distribution of avalanche sizes is a key indicator of scale-invariant behavior in complex systems. Empirical distribution of avalanche sizes allows the characterization of scaling behavior in brain activity. Previous empirical work shows that avalanche sizes often follow approximate power-law distributions across species and human data (Ribeiro et al., 2010; Priesemann et al., 2013). In this work we focus on the operational definition and empirical estimation of avalanche statistics.

We analyzed the statistical properties of avalanche sizes and durations by fitting their probability density functions (PDFs) to power-law models. Specifically, the PDF of avalanche sizes, P(s), was fitted to:

$\begin{array}{ll}P(s)\propto s^{-\tau },& (5)\end{array}$

where τ is the exponent governing the decay pattern of the size distribution. Similarly, the PDF of avalanche durations, P(T), was fitted to:

$\begin{array}{ll}P(T)\propto T^{-\alpha },& (6)\end{array}$

where α is the exponent shaping the duration distribution.

In addition, we examined the scaling relationship between the average avalanche size, 〈s〉(T), and its corresponding duration T:

$\begin{array}{ll}\langle s\rangle (T)\propto T^{\gamma },& (7)\end{array}$

where γ is the scaling exponent that relates size and duration (Friedman et al., 2012). The value of γ was obtained via linear regression on logarithmically transformed data.

To estimate this exponent, we fitted the avalanche size distributions using the discrete maximum likelihood method implemented in the powerlaw Python package (Alstott et al., 2014), which provides robust estimation for power-law exponents in empirical datasets. The value of x_min, defining the lower cutoff above which power-law behavior is assumed, was selected as the one that provides the best compromise between the distributions of avalanche exponents across all sleep states, aiming to minimize systematic biases in the comparison between them.

A fixed temporal binning was used to discretize the neural activity time series, independent of the overall activity level. This approach ensures consistent time resolution across different sleep stages, avoiding potential biases introduced by adaptive binning methods that may conflate changes in activity rate with changes in the temporal structure of avalanches. In addition, we quantified the mean activity per active bin, defined as the total number of detected events divided by the number of bins in which at least one channel registered an event, thereby capturing the average event rate conditional on periods of activity.

Note that we use τ for the avalanche size exponent and α for the avalanche duration exponent. The symbol α is also used in Sections 2.4–3.3, 4.3 for the DFA scaling exponent; meanings are disambiguated by context (avalanche duration vs. DFA) and by explicit references in captions and text.

Finally, we use the compact format x = a(b) to denote x = a±b in the last digit(s). For example, α = 0.84(2) means 0.84 ± 0.02. Unless noted otherwise, b is one standard error of the estimate. The same notation is employed for DFA scaling exponent estimations.

2.4 Detrended fluctuation analysis

In this study, we use Detrended Fluctuation Analysis (DFA) to assess temporal correlations in iEEG signals from the ATLAS, particularly focusing on scaling properties.

The steps for calculating the DFA in this work are as follows (Hardstone et al., 2012):

- Calculate the cumulative sum (signal profile): The cumulative sum {X_i} is computed by successively summing the values of the iEEG signal {x_j} and subtracting the mean value 〈x〉:

$\begin{array}{ll}{X}_i={\displaystyle \sum _{j=1}^i}(x_j-\langle x\rangle ).& (8)\end{array}$

In Equation 8, x_j denotes the j-th sample of the demeaned discrete iEEG signal within a 60-s epoch; X_i is the cumulative-sum (profile) up to index i. Here i ∈ {1, …, N} is the running index of the profile, j is the summation index, and N is the number of samples in the epoch.

- Define windows: A set of window sizes, denoted as Δt, is defined. These sizes are evenly distributed on a logarithmic scale, ranging from a minimum of four samples to the length of the signal.

- Divide the cumulative sum into N overlapping windows: A window of length Δt is chosen, and a series of consecutive windows {w_n(Δt)} of the same size is created from the signal profile, with a 50% overlap.

- Remove linear trend in each window: For each window w_n(Δt), fit a straight line (trend) using a least-squares fit and subtract it from the window to create the detrended window w_nD(Δt):

$\begin{array}{ll}{w}_{nD}(\Delta t)=w_n(\Delta t)-{\text{Linear-Trend}}_n.& (9)\end{array}$

- Calculate the standard deviation of detrended windows: For each detrended window w_nD(Δt), calculate the standard deviation (Fluctuation Function F(Δt)):

$\begin{array}{ll}F(\Delta t)=\sqrt{\frac{1}{N}{\displaystyle \sum _{n=1}^N}{\left[w_{nD}(\Delta t)-\langle w_{nD}(\Delta t)\rangle \right]}^2}.& (10)\end{array}$

- Repeat for different window sizes: Iterate the above steps for all window sizes Δt.

- Plot log(F(Δt)) vs. log(Δt): The slope of this log-log plot (estimated via linear regression) yields the DFA exponent α.

The scaling exponent α characterizes the scaling properties of mean fluctuations across time scales in the iEEG signal. As described in Hardstone et al. (2012):

α < 0.5: Anti-correlated time series (smaller fluctuations at larger scales).

α = 0.5: Uncorrelated (random) process.

0.5 < α < 1: Positively correlated time series.

1 < α < 2: Non-stationary process.

In this work, we use the term long-range temporal correlations (LRTC) strictly for the stationary regime 0.5 < α < 1 (Hardstone et al., 2012). Values α ≈ 1 correspond to 1/f-like (pink) scaling, whereas α > 1 indicate non-stationary, integrated dynamics; in the latter case, long-range dependence typically characterizes the increments of the process rather than the raw signal itself (Poil et al., 2012; Dalla Porta and Copelli, 2019). This convention is adopted to avoid conflating stationary LRTC with non-stationary persistence.

3 Results

3.1 Spectral analysis of wake and sleep stages

To assess how different stages of the sleep-wake cycle modulate brain dynamics we performed spectral analysis on iEEG data using power spectral density (PSD) estimates computed via Welch's method (600-sample windows, 50% overlap). Both absolute and normalized PSDs were analyzed globally and by cortical region. Details of the preprocessing and spectral estimation pipelines are provided in Section 2.2.

3.1.1 Global spectral profiles

Figure 1 summarizes the global spectral patterns across vigilance states. Figure 1A displays the grand-average PSDs (log–log scale), obtained by averaging across all 60-s epochs in the dataset. Figure 1B shows the same spectra after normalizing each curve by the across-epoch mean at each frequency, emphasizing relative shifts and crossover points.

Figure 1

Figure 1. (A) Grand-average power spectral density (PSD) for each vigilance state (log–log scale), obtained by averaging across all 60-s epochs in the dataset; the green line shows the across-epoch mean at each frequency (pooled across states). (B) Same spectra as the left panel, after normalizing each curve by the across-epoch mean at each frequency, highlighting relative shifts and crossover points. During sleep, there is a significant increase in delta band power around 1 Hz, indicating non-REM sleep. Theta band activity increases during wakefulness, particularly in its higher range. Additionally, REM sleep shows a reduction in spectral density between 8 and 20 Hz. Furthermore, non-REM N2 and N3 stages exhibit a distinct peak around 12 Hz in the spectral profile, followed by a drop in the beta band, distinguishing them from wakefulness and REM sleep.

As expected, non-REM sleep stages—particularly N3—exhibit a marked increase in delta power around 1 Hz, with differences exceeding an order of magnitude compared to wakefulness. Theta activity is most pronounced in wakefulness near 7–8 Hz. A distinct peak around 12 Hz, corresponding to sleep spindles, emerges during N2 and N3. REM sleep shows a broadband reduction in power between alpha and low beta bands (8–20 Hz). In the gamma range, all sleep stages display lower power than wakefulness, though with smaller effect sizes.

3.1.2 Regional spectral features

To assess spectral differences across brain states at the local level, we conducted independent (unpaired) Mann–Whitney (U) tests (Virtanen et al., 2020) on Gaussian-fitted amplitudes for canonical frequency bands (delta, theta, alpha, beta, and gamma) across predefined lobes (see Table 1). We compared (i) Wake vs REM and (ii) N2 vs. N3. To control for multiple comparisons over band–lobe combinations, we applied the Benjamini–Hochberg false discovery rate (FDR) procedure (Benjamini and Hochberg, 1995) at α = 0.05.

The contrast between Wake and REM revealed focal differences mainly in alpha amplitude in temporal, parietal and occipital cortices, and delta differences in frontal regions. Conversely, the N2 vs. N3 comparison exhibited more widespread and robust differences, especially in beta and gamma bands, with temporal, frontal, and parietal regions showing strong modulation. All reported differences correspond to statistically significant contrasts, with FDR-corrected p-values below 0.05.

To further explore spatial heterogeneity in spectral organization, we examined region-specific PSDs. Several representative cases are highlighted below.

Figure 2 illustrates Gaussian-fitted alpha-band amplitudes in the calcarine cortex (primary visual area). Black points denote single-epoch (60 s) Gaussian fits, and violins summarize their distribution across epochs for each state. A clear separation is observed, with higher alpha during wakefulness and reduced amplitude during REM. The same plotting convention is used in Figures 3–5.

Figure 2

Figure 2. Gaussian-fitted alpha-band amplitudes in the calcarine cortex (primary visual cortex) during wakefulness and REM sleep. Higher alpha power in wakefulness reflects the characteristic occipital alpha rhythm, attenuated in REM. Black points indicate individual single-epoch (60 s) Gaussian fits; violin plots summarize their distribution for each state. Group differences were assessed with a two-sided Mann–Whitney U test; multiple comparisons were controlled via FDR. The same plotting convention is used in Figures 3–5.

Figure 3

Figure 3. Gaussian-fitted delta-band amplitudes in the medial frontal cortex across sleep stages N2 and N3. Delta amplitude increases significantly in N3 relative to N2 (two-sided Mann–Whitney U; FDR-controlled). This pattern reflects enhanced cortical slow-wave synchrony during deep sleep in frontal-medial regions.

On the other hand, Figure 3 shows that Gaussian-fitted delta-band amplitudes in the medial frontal cortex increase in N3 relative to N2. This elevation reflects enhanced slow-wave synchrony in frontal-medial regions during deep sleep.

Interestingly, the amygdala exhibits a different spectral pattern. As shown in Figure 4, delta amplitude is lower during REM than during wakefulness, diverging from the cortical trend. This result aligns with previous reports of limbic-specific modulation during REM, potentially related to emotional processing and memory reorganization (Kalamangalam et al., 2021).

Figure 4

Figure 4. Gaussian-fitted delta-band amplitudes in the amygdala across wakefulness, REM, N2, and N3. A relative decrease in delta power is observed during REM compared to wakefulness, suggesting distinct modulatory dynamics in this limbic structure.

A similar divergence is found in the hippocampus (Figure 5), where delta amplitude is higher during REM than during wakefulness. This contrasts with most cortical and subcortical areas and supports the view that the hippocampus exhibits unique REM-related dynamics. Prior studies have linked such patterns to memory consolidation mechanisms and the replay of waking experiences (Kalamangalam et al., 2021; von Ellenrieder et al., 2020).

Figure 5

Figure 5. Gaussian-fitted delta-band amplitudes in the hippocampus across wakefulness, REM, N2, and N3. Unlike most cortical regions, delta amplitude is higher during REM than wakefulness.

We present a lobe-level overview with selected, illustrative region-level PSDs, emphasizing the most challenging contrasts (Wake vs. REM; N2 vs. N3). An exhaustive region-by-region survey is beyond the scope of this work.

3.2 Neuronal avalanches across wakefulness and sleep states

Here, we investigate neuronal avalanches—spontaneous cascades of neural activity—across wakefulness and sleep stages using intracranial EEG data from the iEEG Atlas. Our analysis focuses on quantifying how avalanche size distributions vary systematically with vigilance state, without assuming a specific underlying dynamical regime.

3.2.1 Neuronal avalanche results

Neuronal avalanches were identified using a discrete event-based approach, following the established methodology of binning thresholded neural activity in short time windows (Section 2.3).

Given the variability in the number of implanted electrodes across patients, we observed that the power-law scaling of avalanche size distributions was often disrupted in subjects with lower electrode coverage. To mitigate the inevitable bias introduced by spatial sub-sampling, only patients with at least 16 electrodes were included in the analysis. For each selected patient, a uniform subset of 16 electrodes was considered; in cases where more than 16 electrodes were available, a random selection was performed to ensure comparability across subjects. This approach reduces variability due to differences in electrode coverage and allows for a more consistent estimation of avalanche statistics across patients.

Avalanche size distributions were fitted using the discrete maximum likelihood method implemented in the powerlaw package, which estimates the power-law exponent τ and the exponential cutoff parameter λ. To ensure consistency across all conditions and avoid statistical instabilities associated with very small avalanches, the lower cutoff x_min was fixed at 3. This choice excludes the smallest events, which are often dominated by noise or are undersampled due to the temporal discretization, and focuses the analysis on more reliable segments of the avalanche size distribution. Similar fixed thresholds have been employed in previous studies to enhance the robustness of power-law fitting and cross-condition comparisons.

The analysis used Patient 20 as an example of the statistical procedure, showing distinct avalanche dynamics across sleep stages in Figure 6. The probability density function (PDF) of avalanche sizes (s) followed power-law distributions with exponents τ as follows: Wake [τ = 2.30(7), KS = 0.1], REM [τ = 2.12(7), KS = 0.06], N2 [τ = 2.00(7), KS = 0.03], and N3 [τ = 1.78(6), KS = 0.08], with Kolmogorov-Smirnov (KS) distances indicating fit quality. The PDF of avalanche durations (T) also showed power-law behavior with exponents α: Wake [α = 2.70(9), KS = 0.06], REM [α = 2.40(9), KS = 0.03], N2 [α = 2.32(9), KS = 0.03], and N3 [α = 1.95(7), KS = 0.06], supported by KS distances. The scaling relation between mean avalanche size (〈s〉) and duration (T) in log-log scale was linear, with slopes γ: Wake [γ = 1.15(1), R² = 0.95], REM [γ = 1.17(1), R² = 0.95], N2 [γ = 1.18(1), R² = 0.95], and N3 [γ = 1.16(1), R² = 0.97], and R² values confirming a consistent trend across sleep states. These results suggest varying neural avalanche properties depending on sleep stage.

Figure 6

Figure 6. Representative avalanche statistics across sleep stages for Patient 20. (Left) Probability density function (PDF) of avalanche sizes (s), fitted to power-law distributions with exponents τ and Kolmogorov-Smirnov (KS) distances quantifying the goodness-of-fit for each sleep stage: Wake (W), REM (R), N2, and N3. (Middle) PDF of avalanche durations (T), similarly fitted to power laws with exponents α and corresponding KS distances. (Right) Scaling relation between mean avalanche size (〈s〉) and duration (T) in log-log scale, with linear regression fits yielding slopes γ (indicating the scaling exponent) and R² values assessing the fit quality.

The scaling exponent τ—estimated as the slope of the fit in log–log space—indicates a higher frequency of large avalanches during non-REM sleep relative to REM and wakefulness.

Extending this analysis to the entire cohort, Figure 7 shows histograms of the scaling size exponent τ (upper left panel) and duration exponent α (middle left panel) across all patients and vigilance states. Avalanche size exponent (τ) values decreased progressively from wakefulness W (M = 2.21, SD = 0.11, n = 44) and REM (M = 2.12, SD = 0.09, n = 29), to N2 (M = 1.94, SD = 0.09, n = 41) and N3 (M = 1.78, SD = 0.084, n = 41). A similar pattern appeared for the duration exponent (α): W (M = 2.58, SD = 0.15), REM (M = 2.39, SD = 0.14), N2 (M = 2.13, SD = 0.13), and N3 (M = 1.90, SD = 0.11).

Figure 7

Figure 7. State-dependent modulation of avalanche scaling exponents across sleep stages. (a) Probability density functions (PDFs) of avalanche size exponents (τ) for vigilance states Wake (W, red), N3 (yellow), N2 (cyan), and REM (purple). Group means (M), standard deviations (SD), and sample sizes (n) are shown in the inset, together with one-way ANOVA results (F = 175.34, p < 0.0001), indicating significant differences across states. (b) PDFs of avalanche duration exponents (α) for the same states as in (a). Both means and ANOVA results (F = 200.08, p < 0.0001) show a systematic decrease in α with increasing sleep depth. (c) Relationship between α and τ for all subjects, color-coded by vigilance state. Theoretical predictions α = 1+(τ−1)γ for γ = 1.10 (red dashed) and γ = 1.40 (blue dashed) are plotted, with the shaded region covering 1.10 < γ < 1.40. The exponent γ increases from ~1.10 in N3 and N2, through intermediate values in REM, to ~1.40 in Wake, reflecting a state-dependent shift in scaling. (d) Distributions of Kolmogorov–Smirnov (KS) distances for avalanche size fits across states. Boxplots display medians, interquartile ranges, whiskers, and individual data points. Fits are best (lowest KS) during N2, intermediate in W and R, and poorest in N3, with deeper sleep stages showing stronger deviations from power-law scaling. (e–h) Dependence of avalanche measures on analysis parameters. (e) The size exponent τ shows modest variation with bin width (10–20 ms), with changes below 5% across all states. (f) The duration exponent α is more affected, particularly in N3, where variations approach 9%. (g) Mean activity per active bin displays state-dependent differences, being higher in W and R and lower in N2 and N3. (h) Changes in detection threshold (−2.2, −2.0, −1.8 SD) shift the absolute values of τ, but the relative ordering of states is preserved, with W and R consistently above N2 and N3.

Prior to applying one-way ANOVA, assumptions were confirmed: Shapiro–Wilk tests within groups (n ≈ 29–44) yielded p > 0.05, indicating no departures from normality; Levene's test confirmed homogeneity of variances (p ≥ 0.14), conditions under which ANOVA is robust for balanced designs (Zhou et al., 2023; Agbangba et al., 2024; Shaw and Mitchell-Olds, 1993). One-way ANOVA results were highly significant: size exponent (F = 175.34, p = 0.0000) and duration exponent (F = 200.08, p = 0.0000), indicating a systematic shift in avalanche dynamics with increasing sleep depth.

To assess the scaling relationship between avalanche size (τ) and duration (α), we calculated the scaling exponent γ from the slope relating τ and α (Figures 7a, b). Across vigilance states, γ values ranged from approximately 1.10 in N3 sleep, increasing through N2 and REM, and reaching about 1.40 during wakefulness. The points followed a clear positive association between τ and α, with all states falling within the reference band (1.10–1.40; Figure 7c).

The KS distance (Figure 7d), quantifying deviations from power-law dynamics in neuronal avalanches, varied significantly across sleep/wake states (Figure 7). N3 exhibited the strongest deviation (median KS = 0.104, IQR: 0.072–0.116), suggesting the poorest fit to power-law scaling. In contrast, N2 showed the closest alignment (median KS = 0.063, IQR: 0.053–0.075), followed by R (0.076) and W (0.079). The broader distribution in N3 (SD = 0.032, max = 0.167) implies state-dependent variability in neuronal dynamics, with deeper sleep stages (N3) departing most markedly from scale-free dynamics.

The bottom row of Figures 7e–h shows the dependence of avalanche measures on temporal binning within the restricted range of 10–20 ms. For the size exponent τ (Figure 7e), mean values changed only slightly across bin widths, with relative variations below 5% in all states (maximum ≈4.8% in N3). The duration exponent α (Figure 7f) showed larger shifts in mean values, with the strongest variation in N3 (≈8.8%) compared to smaller changes in W, R, and N2 (< 6%).

Mean activity per active bin (Figure 7g) exhibited clear state differences in average levels: N2 (1.514-1.585) showed the highest mean activity, followed by N3 (1.457-1.512), whereas W (1.340-1.490) and R (1.326-1.456) were lower (N2 > N3 > W ≈ R).

For the threshold dependence of τ (Figure 7h), relative variations were modest in N2 (≈5.6%) and N3 (≈5.9%), but substantially larger in R (≈16.1%) and W (≈19.2%). Despite this increased sensitivity, the relative ordering of states was preserved across thresholds, indicating that the state-dependent hierarchy in τ remains robust to moderate changes in detection criteria.

3.3 Detrended fluctuation analysis of sleep and wake states

We applied Detrended Fluctuation Analysis (DFA) to intracranial EEG (iEEG) recordings from the Multicenter iEEG Sleep ATLAS to investigate long-range temporal correlations (LRTCs) across wakefulness and sleep stages. DFA quantifies the presence of LRTCs by estimating a scaling exponent (α), where α = 0.5 indicates uncorrelated noise, 0.5 < α < 1 denotes persistent correlations, and α > 1 suggests non-stationary behavior. Detailed methodological parameters are provided in Section 2.4.

We focused on the intermediate scale range (0.04s–2.5s), which shows the highest variability (as will be shown below) and physiological relevance, in line with previous studies (Hardstone et al., 2012; Meisel et al., 2017; Zhang et al., 2021). Other regions, while potentially interesting, are not the focus here: besides the uniformity of behavior observed both below and above this range, for scales above 2.5s the behavior not only becomes more uniform but also yields poorer statistics due to the limited number of available window sizes inherent to the method. This intermediate range therefore provides sufficient richness to explore meaningful state-dependent differences.

Each of the 5,720 available iEEG channels (60 s epochs sampled at 200 Hz) was analyzed individually. Figure 8 shows the fluctuation function F(Δt) for a representative 60 s epoch. While short time windows (< 0.04s) and very long ones (>2.5s) display relatively uniform patterns across vigilance states, the intermediate range clearly reveals the strongest and most distinctive state-dependent differences, making it the central focus of our analysis.

Figure 8

Figure 8. DFA results for wakefulness and sleep stages from channel 1 of a representative patient (9). Short windows (< 0.04 s) and long windows (>2.5 s) show relatively uniform slopes across states. In contrast, the intermediate range (0.04-2.5s) displays the largest variability and the most distinctive state-dependent differences, making it the focus of our analysis. The legend reports the fitted α values for each state within this range.

3.3.1 DFA global features

We computed grand-average fluctuation functions by pooling all available 60 s epochs and channels within each state, and then estimated global α exponents by fitting a line in the intermediate window range (0.04s–2.5s). Results are displayed in Figure 9. Wakefulness and REM sleep exhibited α values between 0.8 and 1.0, consistent with stationary LRTC. In contrast, N2 and N3 showed α > 1, indicating non-stationary, integrated dynamics with strong temporal persistence, in line with the predominance of slow oscillations in NREM. To avoid terminological ambiguity, we reserve “LRTC” for the stationary regime 0.5 < α < 1; α > 1 reflects non-stationary persistence (with long-range dependence in the increments), rather than stationary LRTC (Hardstone et al., 2012).

Figure 9

Figure 9. Global DFA curves for each vigilance state (grand-averaged across all channels and all available 60-s epochs per state). The intermediate range (0.04-2.5s) shows the greatest variability across states, with α highest for N3, followed by N2, REM, and wakefulness. The legend reports the fitted α values for each state within this range.

3.3.2 DFA variability across regions

To better understand how temporal dynamics distribute across brain regions and sleep stages, we first inspect the spatial patterns of the DFA exponent (α). From Figure 10, it is evident how the DFA exponent reflects the transition from wakefulness and REM to non-REM sleep states, where correlations are positive, moving toward non-stationary states. Only in the wakefulness state can values near 0.5 be observed, indicating a lack of correlations and a tendency toward randomness in the signal. The electrode contact locations (aligned in a straight line) are also discernible, particularly in cases where data are obtained from a single patient.

Figure 10

Figure 10. Brain mapping of the DFA exponent (α) depicted in a transverse plane of the brain. The color bar indicates the value of the exponent α. (a) Waking state, (b) REM, (c) non-REM N2, and (d) non-REM N3 sleep states. The exponent values across the brain, reveal transitions from wakefulness and REM to non-REM sleep states. Positive correlations characterize wakefulness, while non-REM states show non-stationary behavior. Values near 0.5 represent randomness, with distinct electrode contact locations visible, especially in single-patient data.

To quantify these observations, we then calculated, for each lobe and state, key dispersion metrics of the DFA exponent (α), specifically: mean standard deviation (SD), and interquartile range (IQR). SD highlights regions showing greater variability across epochs or subjects, while IQR serves as a robust spread measure less sensitive to outliers. Particularly during deep sleep (N3), regions with elevated SD or IQR may indicate functional heterogeneity or state-dependent modulation.

Figure 11 illustrates the distribution of α across brain lobes and sleep states (Wakefulness, REM, non-REM N2, and non-REM N3). As sleep deepens, α increases monotonically across all lobes, reflecting enhanced temporal persistence and integration.

Figure 11

Figure 11. Boxplots showing the distribution of the Detrended Fluctuation Analysis (DFA) exponent α across six brain lobes (Temporal, Frontal Medial, Frontal Lateral, Insular, Occipital, Parietal) and four sleep states (Wakefulness, REM, non-REM N2, non-REM N3). Each box represents the interquartile range (IQR), with the central line showing the median and whiskers extending to 1.5×IQR; outliers are plotted individually. Colors denote sleep states as indicated in the legend. Across most lobes, α increases progressively from wakefulness to deep non-REM sleep (N3). An exception is the insular lobe, which shows no significant difference between Wake and REM.

To evaluate the statistical significance of these inter-state differences, paired Mann–Whitney U tests (Virtanen et al., 2020) were conducted for each lobe comparing Wake vs. REM and N2 vs. N3, with False Discovery Rate (FDR) correction using the Benjamini–Hochberg method (Benjamini and Hochberg, 1995). Table 2 presents the resulting −log₁₀(p)-values.

Table 2

Table 2. FDR-corrected significance levels (−log₁₀p) for DFA exponent α comparisons between sleep states by brain lobe.

These results indicate that differences between Wake and REM are generally less significant compared to those between N2 and N3. Notably, the insular lobe does not show a significant difference between Wake and REM (denoted “n.s.”), whereas the N2 vs. N3 comparison reveals highly significant differences across most lobes—especially in temporal and frontal lateral regions. This pattern suggests that while deeper sleep stages involve substantial reorganization of cortical temporal dynamics, particularly in medial and temporal areas, the differences between Wake and REM are more subtle and regionally stabilized.

Finally, although a full regional analysis is not the main objective here, to obtain a more complete picture of the DFA exponents α variability across different brain regions, the data from all analyzed channels were categorized by region and plotted in Figure 12. Each colored data point represents the exponent α in a specific brain region. The mean value (depicted by a white-centered circle) and the standard deviation (represented by the bounds of the shaded area) were calculated for each sleep and wakefulness state.

Figure 12

Figure 12. DFA exponent α values across various brain regions are depicted in the figure. Each colored data point corresponds to the α exponent in a specific state: wakefulness (red), non-REM N2 (light blue), non-REM N3 (yellow), and REM (dark blue). The mean value is represented by a white-centered circle, and the standard deviation is indicated by the bounds of the shaded area, calculated for each sleep and wakefulness state. The data is arranged from left to right in ascending order based on the mean standard deviation of the DFA alpha across the four states, computed per region from all 60-s epochs and channels. Within the figure, data related to each region is organized into vertical blocks. Each block, arranged from left to right, includes the region's name along with the states Wakefulness, non-REM N2, non-REM N3, and REM. Below each block, there is a corresponding square in a distinct color, signifying the lobular region to which it belongs. Beneath each square, there is a number that identifies the region according to Table 1. The association of each color with its respective lobular region is depicted in the horizontally extended bar beneath the figure. The most significant variations are observed in the precentral and postcentral gyri, linked to motor and sensory functions. The cuneus region, associated with visual processing, also exhibit notable differences. The hippocampus displays the least variation, with REM sleep showing higher α values than non-REM states. Wakefulness α values are close to 1, while non-REM N2 and N3 states surpass 1.

In Figure 12, the data are displayed from left to right in ascending order according to the mean standard deviation of the DFA alpha across the four states, calculated per region using all 60-s epochs and channels. Each region's data is presented in vertical blocks, including the region's name and states (Wakefulness, non-REM N2, N3, REM). Corresponding colored squares beneath each block indicate the lobular region. The color-region association is shown in the extended bar below the figure.

Upon inspecting Figure 12, notable variability in DFA exponent values is prominent in the precentral and postcentral gyri regions, associated with the motor and sensory cortices, respectively. Subsequently, is the cuneus region, which is situated in the posterior part of the occipital lobe, and stands out for its involvement in fundamental visual processing. At the opposite extreme of Figure 12 lies the hippocampus within the temporal lobe, which plays a pivotal role in behavior inhibition, memory, and spatial perception.

4 Discussion

4.1 Spectral dynamics across vigilance states

Consistent with established sleep physiology, we observed a marked increase in low-frequency (especially delta-band) power during NREM sleep compared to wakefulness. These slow oscillations (below 5 Hz) are hallmark features of deep NREM sleep and are most prominent in frontal–medial cortex, consistent with widespread synchrony observed in intracranial recordings (Brancaccio et al., 2020).

In parallel, sleep spindles (12–15 Hz), were particularly elevated in N2 (and attenuated in N3). Sleep spindles are characteristic of stage-N2 thalamocortical architecture (Staresina et al., 2023). Although we did not assess memory behaviorally, the delta–sigma distribution we observed matches canonical stage-specific spectral profiles.

Wakefulness, in contrast, was dominated by faster theta and alpha oscillations, with reduced delta power. This spectral profile reflects an externally oriented brain state, engaged in sensory processing and encoding of new information (Klimesch, 1999). REM sleep displayed a hybrid profile. Higher-order cortical areas showed fast activity, while primary sensory regions remained in a delta-rich state, consistent with the idea of sensory disconnection during REM (Hobson et al., 2000). This REM profile—fast activity in higher-order cortices alongside persistent limbic delta—is consistent with accounts that REM may protect internally generated reactivation from external interference while permitting the integration of recently acquired information (Hobson et al., 2000; Diekelmann and Born, 2010). In our data, visual regions such as the calcarine cortex showed reduced α-band power during REM stage, compared to wake, suggesting diminished visual input, while limbic areas behaved differently. The amygdala showed decreased delta power in REM (relative to wake), whereas the hippocampus maintained elevated delta power, hinting at continued internal replay or integration (von Ellenrieder et al., 2020).

These differences between wakefulness and REM sleep are particularly relevant—not only because they reflect distinct functional roles, but also because they may serve as robust, physiologically grounded spectro-anatomical markers for distinguishing the two states. Our findings—highlighting, for example, elevated alpha power in the primary visual cortex and increased delta activity in the amygdala during wake, as well as enhanced delta power in the hippocampus during REM—suggest that anatomically specific spectral features (e.g., occipital alpha versus limbic delta) could improve physiologically informed state classification in iEEG datasets.

Importantly, the observed elevation of medial-frontal delta in N3 relative to N2 aligns with findings by Bernardi et al. (2019), who documented source-localized slow-wave amplitude increases in frontal regions during deep sleep. This suggests regional slow-wave amplification as a marker of sleep depth rather than tissue-specific dysfunction. Similarly, the persistence of high delta power in the hippocampus during REM—contrasting with cortical desynchronization—parallels intracranial observations (Moroni et al., 2007; Ferrara et al., 2012), emphasizing a unique hippocampal signature of REM characterized by tonic delta synchronization.

Moreover, the decrease in amygdala delta during REM mirrors findings from emotion-related intracranial studies, where REM selectively attenuates amygdala reactivity overnight (e.g. Wassing et al., 2019). Such limbic divergence exemplifies functionally heterogeneous spectral modulation across subcortical structures.

Finally, the band-specific differences observed in beta and gamma between N2 and N3—particularly their reduction in N3 over fronto-parietal cortices—are consistent with Brancaccio et al. (2020). Taken together, these findings underscore that transitions from N2 to N3 are accompanied not only by slow-wave augmentation but also by broadband desynchronization of higher-frequency activity, possibly reflecting global network downscaling.

4.2 Neuronal avalanches

We analyzed neuronal avalanche size distributions across wakefulness, REM, and NREM sleep stages. A key finding was that the fitted power-law exponent τ varied systematically with vigilance state. Avalanche scaling analysis was performed using a fixed binning approach with a consistent detection threshold of −2.0 standard deviations. The size exponent (τ) ranged from 1.6–2.0 during N2/N3 sleep to 2.0–2.4 in wakefulness and REM states. This indicates a relative increase in the frequency of larger avalanches during wakefulness and REM compared to NREM sleep. Duration exponents (α) showed a clear progression, from 1.7–2.1 in NREM sleep to 2.3–2.9 in REM and wakefulness, reflecting longer avalanche lifetimes in more activated brain states. Spatial scaling (γ) exhibited state-dependent variations, with values of 1.1 in N3 sleep and 1.4 during wakefulness, consistent with increased spatial propagation and integration of avalanches during wakefulness. Furthermore, the increased Kolmogorov–Smirnov (KS) distance observed in power-law fitting during N3 indicates that in this state, avalanche size distributions deviate more significantly from an ideal power-law, suggesting reduced adherence to scale-free dynamics. These quantitative changes in scaling exponents reflect systematic modulations in the spatiotemporal organization of neuronal avalanches across sleep–wake states.

Our results align with prior intracranial and EEG studies showing that vigilance state modulates avalanche statistics (Priesemann et al., 2013).

Limitations in spatial sampling inherently bias the estimation of scaling exponents toward higher values, as previously reported studies (Priesemann et al., 2009). Despite this systematic upward shift, the relative differences in scaling exponents between vigilance states remain preserved and are more reliably distinguished due to the consistent sampling methodology applied across conditions. Consequently, spatial undersampling, while affecting absolute exponent values, does not compromise the comparative analysis of state-dependent avalanche dynamics.

That said, cortical networks during deep sleep are strongly shaped by slow oscillations and synchrony—features that can yield similar scaling without fine tuning. Empirical and modeling work on driven or adaptive networks shows that scale-invariant avalanche patterns can emerge without invoking a phase transition (Priesemann et al., 2014; Martinello et al., 2017; Lombardi et al., 2023). Accordingly, we keep our interpretation descriptive—as state-dependent variation in scaling—without committing to a specific underlying mechanism or operating regime.

Wakefulness and REM sleep showed higher τ values and fewer large avalanches, consistent with more fragmented scaling. This matches observations from MEG and EEG studies (Shriki et al., 2013) and may reflect a functional preference for adaptability and responsiveness over large-scale synchronization.

Our findings are also consistent with longitudinal rodent studies, which have reported a re-expansion of avalanche size distributions during sleep following wake-related shifts in scaling (Xu et al., 2024). These patterns are consistent with sleep-dependent network recalibration, but do not, on their own, demonstrate fine-tuned operation in any specific regime.

Neuronal avalanches have been proposed as signatures of scale-free activity in cortical networks, motivating extensive theoretical and empirical work (Chialvo, 2010; Beggs and Plenz, 2003). Our aim here is not to adjudicate between competing frameworks, but to characterize how avalanche size distributions vary systematically across sleep stages and how these shifts relate to state-dependent processing.

More recent models show that scale-invariant avalanche statistics can emerge in driven, neutral, or adaptive networks without fine tuning (Martinello et al., 2017; Matin et al., 2021; Lombardi et al., 2023). In light of this, we keep our interpretation descriptive—reporting empirical scaling relationships—without assigning them to a particular dynamical context.

4.3 Temporal correlations and information integration windows

DFA revealed state-dependent organization of neural dynamics. At intermediate timescales (0.04-2.5s), wakefulness and REM typically showed 0.6 < α < 1.0, consistent with stationary long-range temporal correlations (LRTC). In contrast, N2 and N3 exhibited α > 1, reflecting non-stationary, integrated persistence consistent with the predominance of slow oscillatory activity in NREM. These scale-dependent patterns align with previous EEG studies (Lee et al., 2002, 2004; Goshvarpour et al., 2013; Ma et al., 2018), and the iEEG dataset used in this study confirms their robustness across cortical and subcortical territories. We keep our analysis descriptive, without interpreting α within any particular theoretical framework.

These temporal regimes reflect a shift in the intrinsic correlation structure of neural dynamics across vigilance states. In wakefulness and REM, moderate long-range correlations suggest persistent but stationary activity patterns. The similarity between REM and wake in short-scale DFA profiles highlights shared dynamic features, despite their different sources of activation.

In NREM, α > 1 denotes non-stationary, integrated persistence, consistent with stronger slow-timescale autocorrelation during deep sleep (Zhang et al., 2021; Meisel et al., 2017).

Our region-specific DFA analysis further revealed heterogeneous modulations of temporal dynamics. While sensory and prefrontal cortices showed marked shifts in α between wakefulness and NREM, the hippocampal α remained comparatively stable across states (Tagliazucchi et al., 2013), a pattern compatible with uninterrupted internal reactivation in two-stage accounts, although we did not assess replay or behavioral memory outcomes.

These findings align with intracranial studies showing asynchronous NREM to REM transitions, starting in occipital regions before reaching frontal areas (Peter-Derex et al., 2023). They also reflect known increases in DFA variability during deep sleep (Goshvarpour et al., 2013). The current results provide a refined anatomical granularity, identifying specific regions with functionally meaningful dynamics.

Overall, our DFA findings demonstrate that sleep influences not only the spectral and spatial organization of neural activity, but also the intrinsic temporal architecture of cortical and subcortical dynamics. The scaling exponent α provides a compact summary of state-dependent changes in temporal structure, highlighting enhanced persistence during deep sleep and dynamic variability across brain regions.

4.4 Integration of multiscale dynamics and implications

By combining oscillatory analysis with avalanche statistics and scaling measures, this work highlights how multiple levels of neural organization collectively contribute to cognitive function across vigilance states. The convergent picture that emerges is that sleep is not a quiescent period, but rather an active neurophysiological process of reorganization—one that potentially lays the groundwork for memory consolidation and brain-wide plasticity.

Slow-wave sleep, in particular, appears to create the necessary conditions for stabilization and integration of new memory representations: strong delta oscillations could synchronize widespread networks, neuronal avalanches skew toward large-scale coactivations, and non-stationary, integrated persistence. These features are those expected to support the repeated reactivation of recently learned neuronal assemblies and the strengthening of synaptic connections according to Hebbian principles (Tononi and Cirelli, 2003; Diekelmann and Born, 2010).

REM sleep, on the other hand, appears to share many dynamical traits with wakefulness, yet it operates in an isolated environment shielded from external perturbations. This unique combination in REM may allow the brain to autonomously explore and refine neural representations—potentially integrating new memories with older networks, resolving interference, and modulating emotional tone—while maintaining a high level of network plasticity (Fagerholm et al., 2015).

Such a role is consonant with hypotheses that REM sleep supports memory consolidation in a complementary fashion to NREM (Diekelmann and Born, 2010).

4.5 Limitations

It is important to recognize the limitations of our interpretations. The links we draw between the observed neural signatures and memory functions are based on consistency with prior empirical and theoretical work (e.g., Tononi and Cirelli, 2003; Diekelmann and Born, 2010; Shew et al., 2009; Tagliazucchi et al., 2013). We thus refrain from making definitive claims that “memory consolidation occurred” in our participants; rather, we argue that the brain states we characterized are compatible with the neurophysiological requirements of consolidation as established in the literature. Notably, extensive research in animals has demonstrated causal connections between specific sleep phenomena and memory—for example, the work of Girardeau et al. (2009) showing that disrupting hippocampal ripples in NREM sleep impairs spatial memory, and evidence that enhancing slow oscillation–spindle coupling in humans can improve recall (Staresina et al., 2015).

These studies bolster the interpretation that the patterns we observe (e.g., abundant NREM slow oscillations with spindles, and an active hippocampus during REM) are indeed meaningful for memory processing. However, future studies should explicitly link our multiscale metrics to behavioral memory performance in order to confirm their functional significance.

Second, our data were obtained from individuals with drug-resistant epilepsy who had intracranial electrodes implanted for clinical monitoring. Although epileptic spikes and seizures were carefully excluded and focused on periods of normal sleep, the generalizability of the findings to healthy brains must be made with caution. Prior work suggests that many of the phenomena we report (such as spectral signatures of sleep stages and avalanche dynamics) have also been observed with noninvasive recordings in healthy populations (Tagliazucchi et al., 2013).

Nonetheless, the potential influence of pathological hyperexcitability or medication effects cannot be completely ruled out. We acknowledge this caveat and view our findings as a first step that should be replicated in further datasets, including those from healthy sleepers if intracranial recordings become available (e.g., during presurgical mapping in non-epileptic tumor cases or with emerging minimally invasive techniques).

Third, our conclusions on temporal correlations and avalanche distributions are based on empirical results validated by statistical significance tests. Linking these patterns to specific theoretical models is beyond the scope of this work.

In all cases, we emphasize a cautious approach: the multiscale descriptors we report are robust and novel characterizations of sleep–wake states, but their exact roles in memory processing remain to be fully determined.

4.6 Conclusions and future directions

This work demonstrates the power of a multiscale analytic approach to sleep neurophysiology, combining spectral, temporal, and network-level assessments to reveal how the sleeping brain balances the demands of information storage and responsiveness.

Beyond physiological interpretation, our findings also highlight opportunities for feature-engineered sleep staging. Region- and band-specific spectral signatures (e.g., posterior alpha in primary visual cortex; delta modulation in limbic and frontal areas) provide physiologically grounded inputs that can aid automated discrimination—particularly for the challenging pairs Wake vs. REM and N2 vs. N3. Complementing these with local temporal descriptors from DFA—specifically the intermediate-scale exponent α (0.04–2.5 s) computed per region, which tends to be ~0.6-1.0 in Wake/REM and >1 in N2/N3—yields a compact, interpretable feature set. Integrating (i) Gaussian-fitted band amplitudes and (ii) region-wise DFA α (optionally alongside spatial co-activation metrics) could improve classifier performance while keeping models grounded in neurophysiology and tractable on 60-s epochs.

On the other hand, future studies should aim to link these multiscale neural signatures with behavioral measures of memory performance, and extend these analyses to larger and more diverse populations, including healthy individuals. Furthermore, collaborative efforts between clinical researchers, computational neuroscientists, and sleep specialists will be essential for uncovering the full cognitive relevance of these dynamic sleep states. The integration of iEEG with interventional techniques such as stimulation or pharmacological modulation may help establish causal links between the observed neural patterns and memory outcomes.

Looking ahead, future efforts could focus on patient-specific modeling and predictive analysis, particularly in clinical contexts where disrupted sleep impacts cognitive function. Establishing direct links between iEEG dynamics and behavioral or memory task performance will be key to clarifying their functional relevance. Additionally, incorporating spatially-aware machine learning methods—such as graph neural networks—may enhance our ability to decode inter-regional interactions. Extending these approaches to longitudinal and developmental datasets could provide new insights into how sleep-related neural representations evolve across the lifespan or in pathological conditions.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: https://ieegatlas.loris.ca/.

Ethics statement

Ethical approval was granted at the Montreal Neurological Institute (MNI) as lead ethics organization (https://ieegatlas.loris.ca/). The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required from the participants or the participants' legal guardians/next of kin in accordance with the national legislation and institutional requirements.

Author contributions

JT: Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – review & editing. MP: Data curation, Investigation, Methodology, Software, Validation, Visualization, Writing – review & editing. MB: Conceptualization, Formal analysis, Funding acquisition, Methodology, Supervision, Validation, Writing – review & editing. OR: Conceptualization, Formal analysis, Funding acquisition, Supervision, Writing – review & editing. FM: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Validation, Writing – review & editing. MA: Conceptualization, Funding acquisition, Investigation, Project administration, Resources, Supervision, Validation, Writing – original draft.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This work was partially supported by CONICET, Argentina (grants PIP 2332 and PIP 0511), UNLP (grants X11/987, X11/895, and X11/978), and IFLP-CONICET (PUE 22920170100066CO).

Acknowledgments

We thank UnCaFiQT-INIFTA (SNCAD) for computational resources.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Generative AI statement

The author(s) declare that Gen AI was used in the creation of this manuscript. Generative AI tools were used to support the writing and editing process of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Agbangba, C. E., Aide, E. S., Honfo, H., and Kakai, R. G. (2024). On the use of post-hoc tests in environmental and biological sciences: a critical review. Heliyon 10:e25131. doi: 10.1016/j.heliyon.2024.e25131

PubMed Abstract | Crossref Full Text | Google Scholar

Alstott, J., Bullmore, E., and Plenz, D. (2014). powerlaw: a python package for analysis of heavy-tailed distributions. PLoS ONE 9:e85777. doi: 10.1371/journal.pone.0085777

PubMed Abstract | Crossref Full Text | Google Scholar

Beggs, J. M., and Plenz, D. (2003). Neuronal avalanches in neocortical circuits. J. Neurosci. 23, 11167–11177. doi: 10.1523/JNEUROSCI.23-35-11167.2003

PubMed Abstract | Crossref Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Statist. Soc. 57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

Crossref Full Text | Google Scholar

Bernardi, G., Betta, M., Ricciardi, E., Pietrini, P., Tononi, G., and Siclari, F. (2019). Regional delta waves in human rapid eye movement sleep. J. Neurosci. 39, 2686–2697. doi: 10.1523/JNEUROSCI.2298-18.2019

PubMed Abstract | Crossref Full Text | Google Scholar

Berry, R. B., Brooks, R., Gamaldo, C. E., Harding, S. M., Marcus, C., Vaughn, B. V., et al. (2012). The aasm manual for the scoring of sleep and associated events. Rules, Terminology and Technical Specifications, Darien, Illinois, American Academy of Sleep Medicine 176, 7.

Google Scholar

Brancaccio, A., Tabarelli, D., Bigica, M., and Baldauf, D. (2020). Cortical source localization of sleep-stage specific oscillatory activity. Sci. Rep. 10:6976. doi: 10.1038/s41598-020-63933-5

PubMed Abstract | Crossref Full Text | Google Scholar

Chialvo, D. R. (2010). Emergent complex neural dynamics. Nat. Phys. 6, 744–750. doi: 10.1038/nphys1803

Crossref Full Text | Google Scholar

Dalla Porta, L., and Copelli, M. (2019). Modeling neuronal avalanches and long-range temporal correlations at the emergence of collective oscillations: continuously varying exponents mimic m/eeg results. PLoS Comput. Biol. 15:e1006924. doi: 10.1371/journal.pcbi.1006924

PubMed Abstract | Crossref Full Text | Google Scholar

De Gennaro, L., and Ferrara, M. (2003). Sleep spindles: an overview. Sleep Med. Rev. 7, 423–440. doi: 10.1053/smrv.2002.0252

PubMed Abstract | Crossref Full Text | Google Scholar

Diekelmann, S., and Born, J. (2010). The memory function of sleep. Nat. Rev. Neurosci. 11, 114–126. doi: 10.1038/nrn2762

PubMed Abstract | Crossref Full Text | Google Scholar

Fagerholm, E. B., Lorenz, R., Scott, G., Hellyer, P. J., Shanahan, M., Sharp, D. J., et al. (2015). Cascades and cognitive state: focused attention incurs subcritical dynamics. J. Neurosci. 35, 4626–4634. doi: 10.1523/JNEUROSCI.3694-14.2015

PubMed Abstract | Crossref Full Text | Google Scholar

Fernandez, L. M., and Lüthi, A. (2020). Sleep spindles: mechanisms and functions. Physiol. Rev. 100, 805–868. doi: 10.1152/physrev.00042.2018

PubMed Abstract | Crossref Full Text | Google Scholar

Ferrara, M., Iaria, G., and De Gennaro, L. (2012). Hippocampal sleep features: relations to human memory performance. Front. Neurol. 3:57. doi: 10.3389/fneur.2012.00057

PubMed Abstract | Crossref Full Text | Google Scholar

Fogel, S. M., and Smith, C. T. (2011). The function of the sleep spindle: a physiological index of intelligence and a mechanism for sleep-dependent memory consolidation. Neurosci. Biobehav. Rev. 35, 1154–1165. doi: 10.1016/j.neubiorev.2010.12.003

PubMed Abstract | Crossref Full Text | Google Scholar

Frauscher, B., von Ellenrieder, N., Dolezalova, I., Bouhadoun, S., Gotman, J., and Peter-Derex, L. (2020). Rapid eye movement sleep sawtooth waves are associated with widespread cortical activations. J. Neurosci. 40, 8900–8912. doi: 10.1523/JNEUROSCI.1586-20.2020

PubMed Abstract | Crossref Full Text | Google Scholar

Friedman, N., Ito, S., Brinkman, B. A., Shimono, M., DeVille, R. L., Dahmen, K. A., et al. (2012). Universal critical dynamics in high resolution neuronal avalanche data. Phys. Rev. Lett. 108:208102. doi: 10.1103/PhysRevLett.108.208102

PubMed Abstract | Crossref Full Text | Google Scholar

Girardeau, G., Benchenane, K., Wiener, S. I., Buzsáki, G., and Zugaro, M. B. (2009). Selective suppression of hippocampal ripples impairs spatial memory. Nat. Neurosci. 12, 1222–1223. doi: 10.1038/nn.2384

PubMed Abstract | Crossref Full Text | Google Scholar

Goshvarpour, A., Abbasi, A., and Goshvarpour, A. (2013). Analysis of electroencephalogram signals in different sleep stages using detrended fluctuation analysis. Int. J. Image, Graph. Signal Proc. 5:49. doi: 10.5815/ijigsp.2013.12.07

Crossref Full Text | Google Scholar

Hardstone, R., Poil, S.-S., Schiavone, G., Jansen, R. M., Nikulin, V. V., Mansvelder, H. D., et al. (2012). Detrended fluctuation analysis: a scale-free view on neuronal oscillations. Front. Physiol. 3:450. doi: 10.3389/fphys.2012.00450

PubMed Abstract | Crossref Full Text | Google Scholar

Hobson, J. A., Pace-Schott, E. F., and Stickgold, R. (2000). Dreaming and the brain: toward a cognitive neuroscience of conscious states. Behav. Brain Sci. 23, 793–842. doi: 10.1017/S0140525X00003976

PubMed Abstract | Crossref Full Text | Google Scholar

Kalamangalam, G. P., Long, S., and Chelaru, M. I. (2020). A neurophysiological brain map: spectral parameterization of the human intracranial electroencephalogram. Clin. Neurophysiol. 131, 665–675. doi: 10.1016/j.clinph.2019.11.061

PubMed Abstract | Crossref Full Text | Google Scholar

Kalamangalam, G. P., Long, S., and Chelaru, M. I. (2021). Neurophysiological brain mapping of human sleep-wake states. Clin. Neurophysiol. 132, 1550–1563. doi: 10.1016/j.clinph.2021.03.014

PubMed Abstract | Crossref Full Text | Google Scholar

Klimesch, W. (1999). Eeg alpha and theta oscillations reflect cognitive and memory performance: a review and analysis. Brain Res. Rev. 29, 169–195. doi: 10.1016/S0165-0173(98)00056-3

PubMed Abstract | Crossref Full Text | Google Scholar

Lee, J.-M., Kim, D.-J., Kim, I.-Y., Park, K. S., and Kim, S. I. (2004). Nonlinear-analysis of human sleep eeg using detrended fluctuation analysis. Med. Eng. Phys. 26, 773–776. doi: 10.1016/j.medengphy.2004.07.002

PubMed Abstract | Crossref Full Text | Google Scholar

Lee, J. M., Kim, D. J., Kim, I. Y., Park, K. S., and Kim, S. I. (2002). Detrended fluctuation analysis of eeg in sleep apnea using mit/bih polysomnography data. Comput. Biol. Med. 32, 37–47. doi: 10.1016/S0010-4825(01)00031-2

PubMed Abstract | Crossref Full Text | Google Scholar

Lombardi, F., Pepić, S., Shriki, O., Tkačik, G., and De Martino, D. (2023). Statistical modeling of adaptive neural networks explains co-existence of avalanches and oscillations in resting human brain. Nat. Comput. Sci. 3, 254–263. doi: 10.1038/s43588-023-00410-9

PubMed Abstract | Crossref Full Text | Google Scholar

Ma, Y., Shi, W., Peng, C.-K., and Yang, A. C. (2018). Nonlinear dynamical analysis of sleep electroencephalography using fractal and entropy approaches. Sleep Med. Rev. 37, 85–93. doi: 10.1016/j.smrv.2017.01.003

PubMed Abstract | Crossref Full Text | Google Scholar

Martinello, M., Hidalgo, J., Maritan, A., Di Santo, S., Plenz, D., and Muñoz, M. A. (2017). Neutral theory and scale-free neural dynamics. Phys. Rev. X 7:041071. doi: 10.1103/PhysRevX.7.041071

Crossref Full Text | Google Scholar

Matin, S., Tenzin, T., and Klein, W. (2021). Scaling of causal neural avalanches in a neutral model. Phys. Rev. Res. 3:013107. doi: 10.1103/PhysRevResearch.3.013107

Crossref Full Text | Google Scholar

Meisel, C., Olbrich, E., Shriki, O., and Achermann, P. (2017). Critical slowing down governs the transition to neuron spiking. PLoS Comput. Biol. 13:e1005287.

PubMed Abstract | Google Scholar

Moroni, F., Nobili, L., Curcio, G., De Carli, F., Fratello, F., Marzano, C., et al. (2007). Sleep in the human hippocampus: a stereo-EEG study. PLoS ONE 2:e867. doi: 10.1371/journal.pone.0000867

PubMed Abstract | Crossref Full Text | Google Scholar

Peter-Derex, L., von Ellenrieder, N., van Rosmalen, F., Hall, J., Dubeau, F., Gotman, J., et al. (2023). Regional variability in intracerebral properties of nrem to rem sleep transitions in humans. Proc. Nat. Acad. Sci. 120:e2300387120. doi: 10.1073/pnas.2300387120

PubMed Abstract | Crossref Full Text | Google Scholar

Plenz, D., and Thiagarajan, T. C. (2007). The organizing principles of neuronal avalanches: cell assemblies in the cortex? Trends Neurosci. 30, 101–110. doi: 10.1016/j.tins.2007.01.005

PubMed Abstract | Crossref Full Text | Google Scholar

Poil, S.-S., Hardstone, R., Mansvelder, H. D., and Linkenkaer-Hansen, K. (2012). Critical-state dynamics of avalanches and oscillations jointly emerge from balanced excitation/inhibition in neuronal networks. J. Neurosci. 32, 9817–9823. doi: 10.1523/JNEUROSCI.5990-11.2012

PubMed Abstract | Crossref Full Text | Google Scholar

Priesemann, V., Munk, M. H., and Wibral, M. (2009). Subsampling effects in neuronal avalanche distributions recorded in vivo. BMC Neurosci. 10:40. doi: 10.1186/1471-2202-10-40

PubMed Abstract | Crossref Full Text | Google Scholar

Priesemann, V., Munk, M. H., and Wibral, M. (2013). Neuronal avalanches differ from wakefulness to deep sleep-evidence from intracranial depth recordings in humans. PLoS Comput. Biol. 9:e1002985. doi: 10.1371/journal.pcbi.1002985

PubMed Abstract | Crossref Full Text | Google Scholar

Priesemann, V., Valderrama, M., Wibral, M., and Le Van Quyen, M. (2014). Spike avalanches in vivo suggest a driven, slightly subcritical brain state. Front. Syst. Neurosci. 8:108. doi: 10.3389/fnsys.2014.00108

PubMed Abstract | Crossref Full Text | Google Scholar

Rasch, B., and Born, J. (2013). About sleep's role in memory. Physiol. Rev. 93, 681–766. doi: 10.1152/physrev.00032.2012

PubMed Abstract | Crossref Full Text | Google Scholar

Ribeiro, T. L., Copelli, M., Caixeta, F., Belchior, H., Chialvo, D. R., Nicolelis, M. A., et al. (2010). Spike avalanches exhibit universal dynamics across the sleep-wake cycle. PLoS ONE 5:e14129. doi: 10.1371/journal.pone.0014129

PubMed Abstract | Crossref Full Text | Google Scholar

Semmlow, J. L., and Griffel, B. (2014). Biosignal and Medical Image Processing. New York: CRC Press.

Google Scholar

Shaw, R. G., and Mitchell-Olds, T. (1993). Anova for unbalanced data: an overview. Ecology 74, 1638–1645. doi: 10.2307/1939922

Crossref Full Text | Google Scholar

Shew, W. L., Yang, H., Petermann, T., Roy, R., and Plenz, D. (2009). Neuronal avalanches imply maximum dynamic range in cortical networks at criticality. J. Neurosci. 29, 15595–15600. doi: 10.1523/JNEUROSCI.3864-09.2009

PubMed Abstract | Crossref Full Text | Google Scholar

Shriki, O., Alstott, J., Carver, F., Holroyd, T., Henson, R. N., Smith, M. L., et al. (2013). Neuronal avalanches in the resting meg of the human brain. J. Neurosci. 33, 7079–7090. doi: 10.1523/JNEUROSCI.4286-12.2013

PubMed Abstract | Crossref Full Text | Google Scholar

Staresina, B. P., Bergmann, T. O., Bonnefond, M., van der Meij, R., Jensen, O., Deuker, L., et al. (2015). Hierarchical nesting of slow oscillations, spindles and ripples in the human hippocampus during sleep. Nat. Neurosci. 18, 1679–1686. doi: 10.1038/nn.4119

PubMed Abstract | Crossref Full Text | Google Scholar

Staresina, B. P., Niediek, J., Borger, V., Surges, R., and Mormann, F. (2023). How coupled slow oscillations, spindles and ripples coordinate neuronal processing and communication during human sleep. Nat. Neurosci. 26, 1429–1437. doi: 10.1038/s41593-023-01381-w

PubMed Abstract | Crossref Full Text | Google Scholar

Steriade, M., Nunez, A., and Amzica, F. (1993). A novel slow (< 1 hz) oscillation of neocortical neurons in vivo: depolarizing and hyperpolarizing components. J. Neurosci. 13, 3252–3265. doi: 10.1523/JNEUROSCI.13-08-03252.1993

PubMed Abstract | Crossref Full Text | Google Scholar

Tagliazucchi, E., von Wegner, F., Morzelewski, A., Brodbeck, V., Jahnke, K., and Laufs, H. (2013). Breakdown of long-range temporal dependence in default mode and attention networks during deep sleep. Proc. Nat. Acad. Sci. 110, 15419–15424. doi: 10.1073/pnas.1312848110

PubMed Abstract | Crossref Full Text | Google Scholar

Tononi, G., and Cirelli, C. (2003). Sleep and synaptic homeostasis: a hypothesis. Brain Res. Bull. 62, 143–150. doi: 10.1016/j.brainresbull.2003.09.004

PubMed Abstract | Crossref Full Text | Google Scholar

Troester, M., Quan, S., of Sleep Medicine, A. A., and Berry, R. (2023). The AASM Manual for the Scoring of Sleep and Associated Events, Version 3. New York: American Academy of Sleep Medicine.

Google Scholar

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., et al. (2020). Scipy 1.0: fundamental algorithms for scientific computing in python. Nat. Methods 17, 261–272. doi: 10.1038/s41592-019-0686-2

PubMed Abstract | Crossref Full Text | Google Scholar

von Ellenrieder, N., Gotman, J., Zelmann, R., Rogers, C., Nguyen, D. K., Kahane, P., et al. (2020). How the human brain sleeps: direct cortical recordings of normal brain activity. Ann. Neurol. 87, 289–301. doi: 10.1002/ana.25651

PubMed Abstract | Crossref Full Text | Google Scholar

Walker, M. P. (2009). The role of sleep in cognition and emotion. Ann. N. Y. Acad. Sci. 1156, 168–197. doi: 10.1111/j.1749-6632.2009.04416.x

PubMed Abstract | Crossref Full Text | Google Scholar

Wassing, R., Lakbila-Kamal, O., Ramautar, J. R., Stoffers, D., Schalkwijk, F., and Van Someren, E. J. (2019). Restless rem sleep impedes overnight amygdala adaptation. Curr. Biol. 29, 2351–2358. doi: 10.1016/j.cub.2019.06.034

PubMed Abstract | Crossref Full Text | Google Scholar

Wilting, J., and Priesemann, V. (2019). 25 years of criticality in neuroscience—established results, open controversies, novel concepts. Curr. Opin. Neurobiol. 58, 105–111. doi: 10.1016/j.conb.2019.08.002

PubMed Abstract | Crossref Full Text | Google Scholar

Xu, Y., Schneider, A. M., Wessel, R., and Hengen, K. B. (2024). Sleep restores an optimal computational regime in cortical networks. Nat. Neurosci. 27, 328–338. doi: 10.1038/s41593-023-01536-9

PubMed Abstract | Crossref Full Text | Google Scholar

Zhang, J., Chang, W.-C., Yu, H.-F., and Dhillon, I. S. (2021). “Fast multi-resolution transformer fine-tuning for extreme multi-label text classification,” in Advances in Neural Information Processing Systems (NeurIPS), 26676–26689.

Google Scholar

Zhou, Y., Zhu, Y., and Wong, W. K. (2023). Statistical tests for homogeneity of variance for clinical trials and recommendations. Contemp. Clin. Trials Commun. 33:101119. doi: 10.1016/j.conctc.2023.101119

PubMed Abstract | Crossref Full Text | Google Scholar