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.

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],

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

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

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

in which 2π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.

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:

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

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:

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.

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):

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

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.

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.

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.

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.

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.

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.

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.

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).

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.

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.