Introduction

During sleep, the human brain engages in highly synchronized activity driven by distinct neural rhythms. In non-rapid eye movement (NREM) sleep, the emergence and coordination of slow oscillations (SOs) and sleep spindles are key features (Hahn et al., 2020; Helfrich et al., 2018; Ngo, Fell, & Staresina, 2020; Schreiner et al., 2022; Staresina et al., 2015; Staresina et al., 2023). However, the brain activation patterns associated with these rhythms, and especially their coordination, remain largely unexplored, and their functional implications are even less understood.

SOs (∼1 Hz) are characterized by fluctuations in neuronal membrane potential between periods of silence (’hyperpolarization’, or ‘DOWN’ state) and excitation (’depolarization’, or ‘UP’ state) (Amzica & Steriade, 2002; Steriade, McCormick, & Sejnowski, 1993). These oscillations originate in neocortical regions, particularly the medial prefrontal cortex (mPFC), during NREM sleep (Dang-Vu et al., 2008; Massimini et al., 2004; Nir et al., 2011). The depolarization phase of SOs is believed to trigger sleep spindles in the thalamus, producing 11–16 Hz oscillations (Fernandez & Lüthi, 2020; Mak-McCully et al., 2017). Thalamic spindles tend to align with the excitable UP-states of SOs and can synchronize with hippocampal ripples (Buzsáki, 2015; Helfrich et al., 2019; Joo & Frank, 2018; Ngo, Fell, & Staresina, 2020). Rodent optogenetics studies have shown that only thalamic spindles phase-locked to the UP-state of cortical SOs enhance memory consolidation, while out-of-phase spindles do not (Latchoumane et al., 2017). This finding is consistent with intracranial evidence from human epilepsy patients, suggesting the thalamus acts as a key relay between the hippocampus and cortex, coordinating their interaction (Coulon, Budde, & Pape, 2012; Ferraris et al., 2021; Schreiner et al., 2022; Staresina et al., 2015).

The triple coupling of SOs, spindles, and ripples is believed to facilitate memory consolidation by synchronizing neuronal activity across brain regions and replaying memory traces from the hippocampus to the neocortex (Diekelmann & Born, 2010; Latchoumane et al., 2017; Liu et al., 2022; Singh, Norman, & Schapiro, 2022; Staresina et al., 2015; Staresina et al., 2023), possibly through relay from thalamic spindles. However, the exact mechanisms of this inter-regional communication during coordinated sleep rhythms remain unclear. Understanding the brain-wide activation associated with these rhythms is crucial, as it reveals how information flows among brain regions, especially during the critical time windows when these rhythms align.

Investigating these sleep-related neural processes in humans is challenging because it requires tracking transient sleep rhythms while simultaneously assessing their widespread brain activation. Recent advances in simultaneous EEG-fMRI techniques provide a unique opportunity to explore these processes. EEG allows for precise event-based detection of neural signal, while fMRI provides insight into the broader spatial patterns of brain activation and functional connectivity (Horovitz et al., 2008; Huang et al., 2024; Laufs, 2008; Laufs, Walker, & Lund, 2007; Schabus et al., 2007; Spoormaker et al., 2010). Previous EEG-fMRI studies on sleep have focused on classifying sleep stages or examining the neural correlates of specific waves (Bergmann et al., 2012; Caporro et al., 2012; Czisch et al., 2009; Fogel et al., 2017; Hale et al., 2016; Ilhan-Bayrakcı et al., 2022; Moehlman et al., 2019; Picchioni et al., 2011). These studies have generally reported that slow oscillations are associated with widespread cortical and subcortical BOLD changes, whereas spindles elicit activation in the thalamus, as well as in several cortical and paralimbic regions. Although these findings provide valuable insights into the BOLD correlates of sleep rhythms, they often do not employ sophisticated temporal modeling (Huang et al., 2024), to capture the dynamic interactions between different oscillatory events, e.g., the coupling between SOs and spindles.

Understanding sleep rhythm-related brain activation is also essential for uncovering its potential cognitive functions. Since performing cognitive tasks during sleep is limited, most neuroimaging studies on human sleep rely on the targeted memory reactivation (TMR) paradigm, where participants are re-exposed to sensory cues during sleep to trigger associated memory reactivation (Cousins et al., 2016; Hu et al., 2020; Oudiette & Paller, 2013; Siefert et al., 2024). For instance, fMRI studies have shown significant hippocampal activation in response to odour re-exposure during sleep (Rasch et al., 2007), while EEG studies suggest that this reactivation is more likely during the UP-state of SOs (Cousins et al., 2016; Göldi et al., 2019; Xia et al., 2023). However, while TMR provides insight into externally induced memory processing during sleep, it offers limited understanding of the endogenous neural mechanisms underlying natural sleep rhythms.

To address this gap, our study investigates brain-wide activation and functional connectivity patterns associated with SO-spindle coupling, and employs a cognitive state decoding approach (Margulies et al., 2016; Yarkoni et al., 2011)—albeit indirectly—to infer potential cognitive functions. In the current study, we used simultaneous EEG-fMRI recordings during nocturnal naps (detailed sleep staging results are provided in the Methods and Table S1) in 107 participants. Although directly detecting hippocampal ripples using scalp EEG or fMRI is challenging, we expected that hippocampal activation in fMRI would coincide with SO-spindle coupling detected by EEG, given that SOs, spindles, and ripples frequently co-occur during NREM sleep. We also anticipated a critical role of the thalamus, particularly thalamic spindles, in coordinating hippocampal-cortical communication.

We found significant coupling between SOs and spindles during NREM sleep (N2/3), with spindle peaks occurring slightly before the SO peak. This coupling was associated with increased activation in both the thalamus and hippocampus, with functional connectivity patterns suggesting thalamic coordination of hippocampal-cortical communication. These findings highlight the key role of the thalamus in coordinating hippocampal-cortical interactions during human sleep and provide new insights into the neural mechanisms underlying sleep-dependent brain communication. A deeper understanding of these mechanisms may contribute to future neuromodulation approaches aimed at enhancing sleep-dependent cognitive function and treating sleep-related disorders.

Results

Simultaneous EEG-fMRI data during nocturnal sleep were collected from 138 healthy subjects (age, 22.56 ± 0.30; 81 females). After excluding subjects due to excessive head motion (mean framewise displacement > 0.5) or insufficient scanning sleep duration (<10 min), 107 participants (age, 22.40 ± 0.34; 63 females) were retained for subsequent analysis, with an average scanning time of sleep of 3.01 ± 0.14 hours (see Methods and Table S1). Since NREM sleep predominantly occurs in early sleep (Carskadon & Dement, 2005), we conducted our scans during the first half of the night to ensure sufficient coverage of NREM sleep, this is also the period when sleep rhythms are most consistently observed and studied in relation to memory consolidation (Maingret et al., 2016; Molle et al., 2011; Schreiner et al., 2021; Staresina et al., 2015; Staresina et al., 2023).

Sleep stages and sleep rhythms

The primary objective of the EEG analysis was to examine the intrinsic characteristics of sleep rhythms and to identify the timing of these SO and spindle events for subsequent fMRI analyses.

We first removed MRI gradient noise from the EEG signal (detailed in Methods, see Fig. S1) and then applied an automated sleep staging algorithm (Vallat & Walker, 2021). The staging results were manually reviewed by two sleep experts to ensure the accurate classification of each sleep stages (Fig. 1b). Upon validation, we confirmed that N2/3 stages predominated in the dataset, with a significantly higher proportion of N2/3 sleep (75.77% ± 1.50%) compared to N1 (14.57% ± 1.06%, t₍₁₀₆₎ = 25.04, p < 1e-4, paired-samples t test) and REM sleep (9.66% ± 0.90%, t₍₁₀₆₎ = 29.60, p < 1e-4, Fig. 1d, Table S1).

Fig. 1

Each sleep stage is characterized by distinct spectral properties and rhythmic waveforms, serving as physiological markers (Fig. 1c). As a sanity check, we compared SO and spindle amplitudes and densities across stages. Both were significantly higher during N2/3 sleep (SO: 25.59 ± 1.49 μV, spindle: 7.39 ± 0.27 μV) compared to N1 (SO: 18.87 ± 1.08 μV, spindle: 6.72 ± 0.25 μV) and REM (SO: 16.67 ± 1.14 μV, spindle: 6.23 ± 0.23 μV, all p < 1e-4, Fig. 1e, Fig. S2). This aligns with the established understanding that N2/3 stages feature more pronounced slow-wave activity and robust spindles (Gaillard & Blois, 1981). Furthermore, we found that the amplitude of the SO DOWN-state (13.62 ± 0.67 μV) was significantly higher than the UP-state (12.98 ± 0.66 μV, t₍₁₀₆₎ = 4.82, p < 1e-4), consistent with previous work (Cash et al., 2009; Dang-Vu et al., 2008). Based on this, we used the DOWN-state to lock the timing of SOs for subsequent EEG-informed fMRI analysis. Descriptive information on sleep rhythms during these stages for all 107 subjects can be found in Table S2-S4.

SO-spindle coupling

SO-spindle coupling is considered important for sleep-dependent memory consolidation. In current study, we found that this coupling primarily occurs during N2/3 sleep stages (2.46 ± 0.06 times/min), significantly more than in N1 (2.23 ± 0.09 times/min, t₍₁₀₆₎ = 2.81, p = 0.0059) and REM stages (2.32 ± 0.09 times/min, t₍₁₀₆₎ = 1.99, p = 0.0492, Fig. 2b, Table S2-S4), consistent with previous findings (Ngo et al., 2013; Staresina et al., 2015). After extracting all EEG epochs where SO-spindle coupling occurred, we analysed their spectral and phase characteristics. The spindles were most likely to occur slightly before the UP state peak of SOs (Fig. 2a, e), aligning with results from both animal studies (Maingret et al., 2016) and human research (Staresina et al., 2015). In our data, this pattern was consistent across subjects (Fig. 2d, Rayleigh test: z = 9.51, p < 1e-4), with the peak of the spindle aligned at an SO phase of −41.61 ± 0.86° (the SO UP state peak is 0°).

Fig. 2

Interestingly, spindle events coupled with SOs had higher amplitudes (7.69 ± 0.26 μV) compared to those not coupled (7.33 ± 0.25 μV, t₍₁₀₆₎ = 4.65, p < 1e-4, Fig. 2c, left). Similarly, SO events coupled with spindles exhibited higher amplitudes (26.28 ± 1.34 μV) than those not coupled (24.56 ± 1.33 μV, t₍₁₀₆₎ = 6.86, p < 1e-4, Fig. 2c, right). These findings suggest coordinated neural processing during SO-spindle coupling. Next, we explore the brain-wide activation associated with these sleep rhythms.

Sleep rhythms associated brain activation

To investigate brain-wide activity during sleep rhythms, we modeled the fMRI blood oxygen level-dependent (BOLD) signal based on EEG-identified timing of SOs and spindle events during N2/3 sleep stages. The design matrix included the main effects of SOs, spindles, and their coupling (Fig. 3a, detailed in Methods). This method is validated in our previous work to effectively capture replay-aligned brain activation during rest (Huang et al., 2024). We found that SOs were associated with positive activation in the thalamus (ROI analysis, t₍₁₀₆₎ = 2.38, p = 0.0096, one-sample t test) and significant deactivation in the neocortex (Fig. 3b, see also Fig. S5), particularly in the default mode network (DMN), which includes the mPFC and posterior cingulate cortex (PCC), consistent with previous research (Dang-Vu et al., 2008; Picchioni et al., 2011; Tagliazucchi et al., 2013). This deactivation of the neocortex, particularly the DMN, at the trough of SOs is intriguing, given that the DMN is thought to encode the internal model of the world or cognitive map (Baldassano, Hasson, & Norman, 2018; Constantinescu, O’Reilly, & Behrens, 2016; Park, Miller, & Boorman, 2021).

Fig. 3

Spindles were correlated with positive activation in the thalamus (ROI analysis, t₍₁₀₆₎ = 15.39, p < 1e-4), the anterior cingulate cortex (ACC), and the putamen, alongside deactivation in the DMN (Fig. 3c). Notably, SO-spindle coupling was linked to significant activation in both the thalamus (ROI analysis, t₍₁₀₆₎ = 3.38, p = 0.0005) and the hippocampus (ROI analysis, t₍₁₀₆₎ = 2.50, p = 0.0070, Fig. 3d). However, no decrease in DMN activity was found during SO-spindle coupling, and DMN activity during SO was significantly lower than during coupling (ROI analysis, t₍₁₀₆₎ = -4.17, p < 1e-4). For more detailed activation patterns, see Table S5-S7. We also varied the threshold used to detect SO events to assess its effect on hippocampal activation during SO-spindle coupling and observed that hippocampal activation remained significant when the percentile thresholds for SO detection ranged between 71% and 80% (see Fig. S6).

To understand the functional implications of these sleep rhythm-related brain activation, we performed open-ended cognitive state decoding using the NeuroSynth database (Yarkoni et al., 2011), examining associations between topic terms and regions of interest derived from the activation patterns of the sleep rhythms. Topic terms were ranked by their z-scores across the rhythms to assess overall functional relevance (see Methods for details). For positive activation patterns (Fig. S7a), SOs were associated with declarative (z = 6.31, p < 1e-4) and episodic memory (z = 5.08, p < 1e-4), while spindles aligned most with working memory (z = 4.08, p < 1e-4). SO-spindle coupling showed the strongest similarity to episodic (z = 3.66, p = 0.0003) and declarative memory (z = 2.21, p = 0.0271). For negative (deactivation) patterns (Fig. S7b), SOs were linked to declarative memory (z = 7.61, p < 1e-4), while spindles resembled task-related activations (z = 7.05, p < 1e-4). SO-spindle coupling had the highest similarity to task-related (z = 9.40, p < 1e-4) and working memory (z = 9.15, p < 1e-4). These findings suggest that internal memory-related regions are more likely to reactivate during SO-spindle coupling, whereas task-related activations, like working memory, are less prominent.

Functional connectivity during SO-spindle coupling

Previous studies have indicated that the hippocampal-cortical dialogue during sleep memory consolidation relies on the thalamus (Coulon, Budde, & Pape, 2012; Ferraris et al., 2021; Latchoumane et al., 2017). To explore this in the human brain, we examined how functional connectivity (FC) changes when spindles couple (versus not) with SOs during sleep, using psychophysiological interaction (PPI) analysis (Friston et al., 1997). Building on the previous GLM (Fig. 3a), this model incorporated two additional regressors: the BOLD signal of the ROI seed and its interaction with SO-spindle events (Fig. 4a, detailed in Methods).

Fig. 4

We first examined whole-brain FC (physiological term) in three ROIs: the hippocampus (bilateral, anatomically defined), thalamus (bilateral, anatomically defined), and mPFC (bilateral, functionally defined). The hippocampus-based FC was mainly connected with the DMN, particularly the mPFC (Fig. 4b). The thalamus-based FC was primarily linked to the basal ganglia and ACC (Fig. 4c), and the mPFC-based FC showed strong connections with the hippocampus and PCC (Fig. 4d).

Crucially, the interaction between FC and SO-spindle coupling revealed that only the functional connectivity of hippocampus -> thalamus (ROI analysis, t₍₁₀₆₎ = 1.86, p = 0.0328) and thalamus -> mPFC (ROI analysis, t₍₁₀₆₎ = 1.98, p = 0.0251) significantly increased during SO-spindle coupling, with no significant changes in all other pathways (Fig. 4e). We also conducted PPI analyses for the other two events (SOs and spindles), and neither yielded significant connectivity changes in the three ROIs, as all failed to survive whole-brain FWE correction at the cluster level (p < 0.05). Together, these findings suggest that the thalamus, likely via spindles, coordinates hippocampal-cortical communication selectively during SO-spindle coupling, but not isolated SOs or spindle events alone.

Discussion

Using simultaneous EEG and fMRI during nocturnal naps in 107 subjects, we identified robust coupling between spindles and the UP-state of SOs in deep NREM sleep (N2/3 stages). This coupling was linked to distinct neural activation patterns, including significant thalamic activation during spindle occurrences, reduced activity in the DMN during SOs trough, and selective hippocampal activation when spindles and SOs coupled. Furthermore, functional connectivity analyses revealed increased connectivity from the hippocampus to the thalamus and from the thalamus to the mPFC, only during SO-spindle coupling, but not isolated SOs or spindle events alone. These results shed light on how the thalamus may coordinate hippocampal-cortical interactions during NREM sleep.

Our data provide insights into the neurobiological underpinnings of these sleep rhythms. SOs, originating mainly in neocortical areas such as the mPFC, alternate between DOWN- and UP-states. The thalamus generates sleep spindles, which in turn couple with SOs. Our finding that spindle peaks consistently occurred slightly before the UP-state peak of SOs (in 83 out of 107 participants), concurs with prior studies, including Schreiner et al. (2021). Yet it differs from some results suggesting spindles might peak right at the SO UP-state (Hahn et al., 2020; Helfrich et al., 2018). Such discrepancies could arise from differences in detection algorithms, participant age (Helfrich et al., 2018), or subtle variations in cortical-thalamic timing. Nonetheless, these results underscore the importance of coordinated SO-spindle interplay in supporting sleep-dependent processes.

When modeling the timing of these sleep rhythms in the fMRI, we observed hippocampal activation selectively during SO-spindle events. This suggests the possibility of triple coupling (SOs–spindles– ripples), even though our scalp EEG was not sufficiently sensitive to detect hippocampal ripples—key markers of memory replay (Buzsáki, 2015). Recent iEEG evidence indicates that ripples often co-occur with both spindles (Ngo, Fell, & Staresina, 2020) and SOs (Staresina et al., 2015; Staresina et al., 2023). Therefore, the hippocampal involvement during SO-spindle events in our study may reflect memory replay from the hippocampus, propagated via thalamic spindles to distributed cortical regions.

The thalamus, known to generate spindles (Halassa et al., 2011), plays a key role in producing and coordinating sleep rhythms (Coulon, Budde, & Pape, 2012; Crunelli et al., 2018), while the hippocampus is found essential for memory consolidation (Buzsáki, 2015; Diba & Buzsáki, 2007; Singh, Norman, & Schapiro, 2022). The increased hippocampal and thalamic activity, along with strengthened connectivity between these regions and the mPFC during SO-spindle events, underscores a hippocampal-thalamic-neocortical information flow. This aligns with recent findings suggesting the thalamus orchestrates neocortical oscillations during sleep (Schreiner et al., 2022). The thalamus and hippocampus thus appear central to memory consolidation during sleep, guiding information transfer to the neocortex, e.g., mPFC.

An intriguing aspect of our findings is the reduced DMN activity during SOs when modeled at the SO trough (DOWN-state). This reduced DMN activity may reflect large-scale neural inhibition characteristic of the SO trough. The DMN is typically active during internally oriented cognition (e.g., self-referential processing or mind-wandering) and is suppressed during external stimuli processing (Yeshurun, Nguyen, & Hasson, 2021). It is unlikely, however, that this suppression of DMN during SO events is related to a shift from internal cognition to external responses given it is during deep sleep time. Instead, it could be driven by the inherent rhythmic pattern of SOs, which makes it difficult to separate UP-from DOWN-states (the two temporal regressors were highly correlated, and similar brain activation during SOs events was obtained if modelled at the SO peak instead, Fig. S5). Since the amplitude at the SO trough is consistently larger than that at the SO peak, the neural activation we detected may primarily capture the large-scale inhibition from DOWN-state. Interestingly, no such DMN reduction was found during SO-spindle coupling, implying that coupling may involve distinct neural dynamics that partially re-engage DMN-related processes, possibly reflecting memory-related reactivation. Future research using high-temporal-resolution techniques like iEEG could clarify these possibilities.

To explore functional relevance, we employed an open-ended cognitive state decoding approach using meta-analytic data (NeuroSynth: Yarkoni et al. (2011)). Although this method usefully generates hypotheses about potential cognitive processes, particularly in the absence of a pre- and post-sleep memory task, it is inherently indirect. Many cognitive terms showed significant associations (16 of 50), such as “episodic memory,” “declarative memory,” and “working memory.” We focused on episodic/declarative memory given the known link with hippocampal reactivation (Diekelmann & Born, 2010; Staresina et al., 2015; Staresina et al., 2023). Nonetheless, these inferences regarding memory reactivation should be interpreted cautiously without direct behavioral measures. Future research incorporating explicit tasks before and after sleep would more rigorously validate these potential functional claims.

Despite providing new insights, our study has several limitations. First, our scalp EEG did not directly capture hippocampal ripples, preventing us from conclusively demonstrating triple coupling. Second, the combination of EEG-fMRI and the lack of a memory task limit our ability to parse fine-grained BOLD responses at the DOWN- vs. UP-states of SOs and link observed activations to behavioral outcomes. Third, the use of large anatomical ROIs may mask subregional contributions of specific thalamic nuclei or hippocampal subfields. Finally, without a memory task, we cannot establish a direct behavioral link between sleep-rhythm-locked activation and memory consolidation. Future studies combining techniques such as ultra-high-field fMRI or iEEG with cognitive tasks may refine our understanding of subregional network dynamics and functional significance during sleep.

In summary, we found that SO-spindle coupling during deep NREM sleep is associated with distinct activation in both the thalamus and hippocampus, accompanied by increased connectivity from the hippocampus to the thalamus and from the thalamus to the mPFC. These results emphasize the central role of the thalamus in orchestrating hippocampal-cortical interactions and highlight potential neural substrates underlying the complex interplay of SOs, spindles, and, possibly, ripples. Our open-ended decoding analyses suggest that these coordinated sleep events may be relevant for memory-related processes; however, future experiments incorporating explicit behavioral measures are necessary to confirm their functional significance. By revealing how these rhythms modulate widespread brain networks, our findings open avenues for developing targeted neuromodulatory strategies to enhance sleep-dependent cognitive functions and potentially address sleep-related disorders.

Methods

Participants & Protocol

A total of 138 healthy adults (age: 22.56 ± 0.30; 81 females) were recruited for the study. All participants had either normal vision or vision corrected to normal standards. None had any past psychiatric or neurological disorders. Prior to participation, they were screened for eligibility to undergo magnetic resonance imaging (MRI). The experiment received approval from Peking University ethics committee (reference number: #2015-26, #2018-11-05). Every participant provided written informed consent.

At the start of the experiment, subjects arrived early at the laboratory for registration. They signed informed consent forms, washed their hair, changed their clothes, and put on the EEG cap in preparation for the experiment. After these preliminary steps, participants entered the MRI scanning room by around 20:30, settled down for sleep, and data collection began. Given concern for potential data drift and decline in the signal-to-noise ratio due to extended and continuous fMRI scanning, the data acquisition was segmented into several sessions for data collection. A “sleep” session was concluded either when a participant was fully awake (unable to fall asleep again or needed to use the restroom) or after 8192 seconds (the maximum scanning time of the scan, equating to 4096 volumes with 2 seconds TR).

After removing subjects with excessive head movements (mean FD > 0.5) or those whose sleep durations were under 10 minutes, the final participant count for further analyses was 107 (age: 22.40 ± 0.34; 63 females) with a mean scanning duration of sleep of 3.01 ± 0.14 hours.

EEG data acquisition

EEG was recorded simultaneously with fMRI data using an MR-compatible EEG amplifier system (BrainAmps MR-Plus, Brain Products, Germany), along with a specialized electrode cap. The recording was done using 64 channels in the international 10/20 system, with the reference channel positioned at FCz. In order to adhere to polysomnography (PSG) recording standards, six electrodes were removed from the EEG cap: one for electrocardiogram (ECG) recording, two for electrooculogram (EOG) recording, and three for electromyogram (EMG) recording. EEG data was recorded at a sample rate of 5000 Hz, the resistance of the reference and ground channels was kept below 10 kΩ, and the resistance of the other channels was kept below 20 kΩ. To synchronize the EEG and fMRI recordings, the BrainVision recording software (BrainProducts, Germany) was utilized to capture triggers from the MRI scanner. The EEG sampling was synchronized to the MR scanner’s 10 MHz gradient system clock, ensuring a stable gradient artifact shape over time and enabling accurate artifact removal. This was achieved via the standard clock synchronization interface of the EEG amplifier, minimizing timing jitter and drift.

MRI data acquisition

All MRI data were acquired using a 20-channel head coil on a research-dedicated 3-Tesla Siemens Magnetom Prisma MRI scanner. Earplugs and cushions were provided for noise protection and head motion restriction. We chose the 20-channel head coil because it was compatible with our EEG system and ensured sufficient signal-to-noise ratio (SNR) for both EEG and fMRI acquisition. The EEG electrode cables were guided through the lateral and posterior openings of the head coil, secured with foam padding to reduce motion and minimize MR-related artifacts. Moreover, given the extended nature of nocturnal sleep recordings, the 20-channel coil helped maintain participant comfort while still achieving high-quality simultaneous EEG-fMRI data.

For registration purposes, high-resolution T1-weighted anatomical images were acquired using a three-dimensional magnetization prepared rapid acquisition gradient-echo sequence (sequence specification: number of slices = 192, repetition time (TR) = 2,530 ms, echo time (TE) = 2.98 ms, inversion time = 1,100 ms, voxel size = 0.5 × 0.5 × 1 mm³, flip angle (FA) = 7°). The participants were asked to lie quietly in the scanner during data acquisition.

Then, the “sleep” session began after the participants were instructed to try and fall asleep. For the functional scans, whole-brain images were acquired using k-space and steady-state T2*-weighted gradient echo-planar imaging (EPI) sequence that is sensitive to the BOLD contrast. This measures local magnetic changes caused by changes in blood oxygenation that accompany neural activity (sequence specification: 33 slices in interleaved ascending order, TR = 2000 ms, TE = 30 ms, voxel size = 3.5 × 3.5 × 4.2 mm³, FA = 90°, matrix = 64 × 64, gap = 0.7 mm). A relatively long TR and larger voxel size were chosen to optimize SNR and reduce susceptibility-related distortions, which are critical in EEG-fMRI sleep studies where head motion and physiological noise can be substantial. The longer TR allowed whole-brain coverage with sufficient temporal resolution, while the larger voxel size helped enhance BOLD sensitivity and minimize partial volume effects in deep brain structures (e.g., the thalamus and hippocampus), which are key regions of interest in this study.

EEG data preprocessing and sleep stage scoring

EEG data collected inside MRI scanner was contaminated by imaging, ballistocardiographic (BCG) and ocular artifacts. To deal with these artifacts, the preprocessing processes was divided into two parts: automated processing and manual processing.

The automated processing was carried out using custom Python scripts and the MNE toolkit (Gramfort et al., 2013). This encompassed the following procedures: I. The raw EEG data were downsampled from 5000 Hz to 2048 Hz. This downsampled rate accelerates subsequent analyses while retaining the necessary precision to distinguish MR noise. II. Imaging artifacts were removed using the Average Artefact Subtraction method (Allen, Josephs, & Turner, 2000). III. We then applied 0.1 - 30 Hz band-pass finite impulse response (FIR) filters to all channels. The EEG and EOG electrodes were re-referenced by subtracting the value from the M1 or M2 electrode (whichever was from the opposite hemisphere). The outer EMG electrodes were re-referenced to the middle one. IV. The signal space projection algorithm was applied to extract all epochs based on TRs and remove residual MR noise, as well as automatically extracting all heartbeat epochs using ECG electrode data for the removal of BCG noise. V. The data were further downsampled to 100 Hz and subjected to Independent Component Analysis (ICA).

Following this, the manual processing is carried out with labelling and removing the physiological artifacts ICs belonging to eye, muscle and residual BCG through the use of custom MATLAB scripts and the EEGLAB toolkit (Delorme & Makeig, 2004). Sleep stages were then scored by YASA toolkit (Vallat & Walker, 2021) after preprocessing, necessitating data from electrodes C3, left EOG, and left EMG. The classification of N2/3 stages were further verified by experts for validity.

MRI data preprocessing

Results included in this manuscript come from preprocessing performed using fMRIPrep 21.0.4 (Esteban et al., 2020; Esteban et al., 2019), which is based on Nipype 1.6.1 (Gorgolewski et al., 2011). Many internal operations of fMRIPrep use Nilearn 0.8.1 (Abraham et al., 2014), mostly within the functional processing workflow. For more details of the pipeline, see https://fmriprep.readthedocs.io/en/latest/workflows.html.

Conversion of data to the brain imaging data structure standard

To facilitate further analysis and sharing of data, all study data were arranged according to the Brain Imaging Data Structure (BIDS) specification using dcm2bids tool, which is freely available from https://unfmontreal.github.io/Dcm2Bids/.

Anatomical data preprocessing

The T1-weighted (T1w) image was corrected for intensity non-uniformity (INU) with N4BiasFieldCorrection (Tustison et al., 2010), distributed with ANTs 2.3.3 (Avants et al., 2008), and used as T1w-reference throughout the workflow. The T1w-reference was then skull-stripped with a Nipype implementation of the antsBrainExtraction.sh workflow (from ANTs), using OASIS30ANTs as target template. Brain tissue segmentation of cerebrospinal fluid (CSF), white-matter (WM) and gray-matter (GM) was performed on the brain-extracted T1w using fast (FSL 6.0.5.1:57b01774, RRID:SCR_002823) (Zhang, Brady, & Smith, 2001). Volume-based spatial normalization to one standard space (MNI152NLin2009cAsym) was performed through nonlinear registration with antsRegistration (ANTs 2.3.3), using brain-extracted versions of both T1w reference and the T1w template. The following template was selected for spatial normalization: ICBM 152 Nonlinear Asymmetrical template version 2009c (Fonov, 2009) [RRID:SCR_008796; TemplateFlow ID: MNI152NLin2009cAsym].

Functional data preprocessing

First, a reference volume and its skull-stripped version were generated using a custom methodology of fMRIPrep. Head-motion parameters with respect to the BOLD reference (transformation matrices, and six corresponding rotation and translation parameters) are estimated before any spatiotemporal filtering using mcflirt (Jenkinson et al., 2002) [FSL 6.0.5.1:57b01774]. BOLD runs were slice-time corrected to 0.96s (0.5 of slice acquisition range 0s-1.92s) using 3dTshift from AFNI (Cox & Hyde, 1997). The BOLD time-series (including slice-timing correction when applied) were resampled onto their original, native space by applying the transforms to correct for head-motion. These resampled BOLD time-series will be referred to as preprocessed BOLD in original space, or just preprocessed BOLD. The BOLD reference was then co-registered to the T1w reference using mri_coreg (FreeSurfer) followed by flirt (Jenkinson & Smith, 2001) [FSL 6.0.5.1:57b01774] with the boundary-based registration (Greve & Fischl, 2009) cost-function. Co-registration was configured with six degrees of freedom. Several confounding time-series were calculated based on the preprocessed BOLD: framewise displacement (FD), DVARS and three region-wise global signals. FD was computed using two formulations following Power (absolute sum of relative motions) (Power et al., 2014) and Jenkinson (relative root mean square displacement between affines, Jenkinson et al. (2002)). FD and DVARS are calculated for each functional run, both using their implementations in Nipype (following the definitions by Power et al. (2014)). The three global signals are extracted within the CSF, the WM, and the whole-brain masks. Additionally, a set of physiological regressors were extracted to allow for component-based noise correction (CompCor) (Behzadi et al., 2007). Principal components are estimated after high-pass filtering the preprocessed BOLD time-series (using a discrete cosine filter with 128s cut-off) for the two CompCor variants: temporal (tCompCor) and anatomical (aCompCor). tCompCor components are then calculated from the top 2% variable voxels within the brain mask. For aCompCor, three probabilistic masks (CSF, WM and combined CSF+WM) are generated in anatomical space. The implementation differs from that of Behzadi et al. (2007) in that instead of eroding the masks by 2 pixels on BOLD space, the aCompCor masks are subtracted a mask of pixels that likely contain a volume fraction of GM. This mask is obtained by thresholding the corresponding partial volume map at 0.05, and it ensures components are not extracted from voxels containing a minimal fraction of GM. Finally, these masks are resampled into BOLD space and binarized by thresholding at 0.99 (as in the original implementation). Components are also calculated separately within the WM and CSF masks. For each CompCor decomposition, the k components with the largest singular values are retained, such that the retained components’ time series are sufficient to explain 50 percent of variance across the nuisance mask (CSF, WM, combined, or temporal). The remaining components are dropped from consideration. Frames that exceeded a threshold of 0.5 mm FD or 1.5 standardised DVARS were annotated as motion outliers.

The BOLD time-series were resampled into standard space, generating a preprocessed BOLD run in MNI152NLin2009cAsym space. First, a reference volume and its skull-stripped version were generated using a custom methodology of fMRIPrep. All resamplings can be performed with a single interpolation step by composing all the pertinent transformations (i.e., head-motion transform matrices, susceptibility distortion correction when available, and co-registrations to anatomical and output spaces). Gridded (volumetric) resamplings were performed using antsApplyTransforms (ANTs), configured with Lanczos interpolation to minimize the smoothing effects of other kernels (Lanczos, 1964). Non-gridded (surface) resamplings were performed using mri_vol2surf (FreeSurfer).

EEG sleep rhythms analysis

In line with other studies examining the coordination between SOs and spindles, particularly their spatial distribution (Massimini et al., 2004; Molle et al., 2011), EEG data from the F3 electrode were used to quantify the occurrence of SOs and spindles. Using the established detection algorithms (detailed below), we identified SOs and sleep spindles for each subject.

Detection of SOs

Data were first bandpass-filtered between 0.16 and 1.25 Hz (Butterworth filter, order 3, bidirectional filtering for zero phase). After identifying all positive-to-negative zero crossings, potential SOs were defined based on the interval between consecutive zero crossings (ranging from 0.8 s to 3 s). For each potential SO, we calculated the amplitude range (peak minus trough). Only those exceeding the 75th percentile in amplitude were labeled as SOs, following previous work (Schreiner et al., 2021).

Detection of sleep spindles

Data were bandpass-filtered between 12 and 16 Hz (Butterworth filter, order 3, bidirectional filtering for zero phase). The root mean square (RMS) of the filtered signal was computed with a 200 ms sliding time window. The amplitude threshold for spindles was set at the 75th percentile of RMS values. Any RMS values that surpassed this threshold for 0.5 s to 3 s were identified as spindles (Staresina et al., 2015).

Detection of SO-spindle couplings

From the detected SOs and spindles, we identified the peak time of each spindle. Within each SO interval, we checked for a spindle peak; if found, it was recorded as an SO-spindle coupling event. For every SO-spindle coupling event, an epoch was created time-locked to the SO trough (as the central reference), following Schreiner et al. (2021). We extracted data in a [-4 s to 4 s] window around this point, forming the epoch for each coupling event.

We note that the above methods for detecting SOs, spindles, and their couplings were originally developed for N2 and N3 sleep data, based on the specific characteristics of these stages. These methods are widely recognized in sleep research (Hahn et al., 2020; Helfrich et al., 2019; Helfrich et al., 2018; Ngo, Fell, & Staresina, 2020; Schreiner et al., 2022; Schreiner et al., 2021; Staresina et al., 2015; Staresina et al., 2023). However, because this percentile-based detection approach will inherently identify a certain number of events if applied to other stages (e.g., N1 and REM), the nature of these events in those stages remains unclear compared to N2/N3. We nevertheless identified and reported the detailed descriptive statistics of these sleep rhythms in all sleep stages, under the same operational definitions, both for completeness and as a sanity check. Within the same subject, there should be more SOs, spindles, and their couplings in N2/N3 than in N1 or REM (see also Figure S2-S4, Table S1-S4).

Event-related potentials (ERP) analysis

After completing the detection of each sleep rhythm event, we performed ERP analyses for SOs, spindles, and coupling events in different sleep stages. Specifically, for SO events, we took the trough of the DOWN-state of each SO as the zero-time point, then extracted data in a [-2 s to 2 s] window from the broadband (0.1–30 Hz) EEG and used [-2 s to -0.5 s] for baseline correction; the results were then averaged across 107 subjects (see Fig. S2a). For spindle events, we used the peak of each spindle as the zero-time point and applied the same data extraction window and baseline correction before averaging across 107 subjects (see Fig. S2b). Finally, for SO-spindle coupling events, we followed the same procedure used for SO events (see Fig. 2a, Figs. S3–S4).

Time-frequency analysis

We performed this analysis on the SO-spindle coupling epochs using the FieldTrip toolbox (Oostenveld et al., 2011). Fourier analysis was executed on the data, utilizing a sliding time window with a forward shift of 50 ms. The window length was defined as five cycles of the present frequency (frequency range: 1-30 Hz, step size: 1 Hz). Before Fourier analysis, windowed data segments underwent multiplication by a Hanning taper. Subsequently, power values were normalized to z-scores within the time range of [-4 s to 4 s], selecting a longer duration to prevent edge artifacts from low-frequency activity in the target time window of [-1.5 s to 1.5 s].

Preferred phase analysis

All SO-spindle coupling epochs were filtered to the SO frequency band (0.16-1.25 Hz, 3rd-order Butterworth filter), and the instantaneous SO phase was determined through the Hilbert transform. Using the saved spindle peak times for each SO-spindle coupling, the preferred phase was deduced. The distribution of these phases across participants was then tested for uniformity using the Rayleigh test.

fMRI GLM analysis

All fMRI GLM analyses were conducted using only data from the N2/3 stages of NREM sleep. Segments of sleep data from these stages were concatenated for each participant. After pre-processing, the fMRI data was smoothed with an 8 mm FWHM kernel. All images underwent a temporal high-pass filter (width of 100 s), and the autocorrelation of the hemodynamic response was modelled via an AR(1) model. We then conducted General Linear Model (GLM) analyses to identify whole-brain activations associated with sleep rhythms and their coupling. We incorporated nuisance regressors estimated during the preprocessing with fMRIprep: six rigid-body motion-correction parameters identified during realignment (comprising three translation and rotation parameters), the mean White Matter, and the mean Cerebral Spinal Fluid signals. The cosine drift model was employed to detect and eliminate data fluctuations, thereby enhancing the accuracy of actual experimental effects. The GLMs were designed to evaluate the influence of experimental conditions on BOLD responses. All whole-brain analyses were thresholded using whole-brain FWE corrected at the cluster level (p < 0.05), with a cluster-forming voxel threshold of p_unc. < 0.001, unless otherwise specified.

The fMRI GLM analysis was designed to look for brain activations pertinent to SO, spindle, and SO-spindle coupling. Drawing from the EEG results concerning the timing of SO and spindle, we identified the trough of the SO as the onset for SO events and the peak of the spindle as the onset for spindle events. We chose the trough of the SO and the peak of the spindle as the onset points, primarily due to the larger amplitudes observed at these moments (thereby more robust to noise). This was also the common choice in previous EEG studies on SO and spindles (Ngo, Fell, & Staresina, 2020; Staresina et al., 2015; Staresina et al., 2023). To model SO-spindle coupling in fMRI, we followed Schreiner et al. (2021), and defining SO-spindle complex as the peak of spindle appeared within 1.5 seconds after SO trough. When such a peak was found, that particular SO was labelled as an SO-spindle event, with its SO trough labelled as the event onset. After convolving these events with the canonical Hemodynamic Response Function (HRF) of SPM and downsampling them to align with the fMRI time resolution, they were integrated as regressors in the design matrix for GLM analysis, allowing us to determine the effects of SO and spindle, as well as the SO-spindle coupling (Fig. 3b-d).

We have also conducted a control analysis aimed to distinguish brain activations associated with DOWN- vs. UP state of SOs. Drawing again from EEG data on SO and spindle, we randomly selected half of the SO events using the peak as the onset and the other half using the trough, to reduce collinearity between the two regressors. As with the main analysis, we identified the spindle peak as the onset for spindle events (Fig. S5).

ROI analysis

The objective of our ROI analysis was to test the activation of specific brain regions during distinct sleep rhythms. Within each ROI, the beta values at the subject-level were averaged across all voxels to facilitate subsequent statistical inferences (one-sample t test). We defined the mDMN and mPFC ROIs functionally based on the meta-analysis results in Neurosynth (Yarkoni et al., 2011), by the keyword “medial default mode network” and “ventromedial prefrontal cortex” respectively. In addition to mDMN and mPFC, we defined other ROIs anatomically, using the high-resolution Harvard-Oxford Atlas probabilistic atlas (Desikan et al., 2006).

Open-ended Cognitive State Functional Decoding

For functional decoding, we followed the approach in Margulies et al. (2016). We used the NeuroSynth (Yarkoni et al., 2011) meta-analytic database (www.neurosynth.org) to assess topic terms associated with the sleep rhythms activation patterns. ROI masks were created from the whole-brain activation patterns associated with the three sleep rhythms in Fig. 3b-d. These activation patterns were binarized using a threshold of p < 0.001, forming a brain mask for each sleep rhythm.

For each ROI mask, the analysis output was a z-statistic associated with a feature term, and the terms were ranked based on the weighted sum of their z-scores. The feature terms were derived from the most general set of 50 topic terms available in the NeuroSynth v3 dataset (https://github.com/NeuroanatomyAndConnectivity/gradient_analysis/blob/master/gradient_data/neurosynth/v3-topics-50-keys.txt).

From the 50 topics, 16 exceeded the threshold of z > 1.96 (p < 0.05). Two terms were removed as they did not represent coherent cognitive functions, leaving 14 relevant topic terms (Fig. S7). These included: episodic memory, declarative memory, working memory, task representation, language, learning, faces, visuospatial processing, category recognition, cognitive control, reading, cued attention, inhibition, and action.

PPI analysis

PPI (psychophysiological interaction) analysis is designed to quantify context-dependent alterations in effective connectivity between brain regions. For instance, it can determine whether hippocampus-based functional connectivity with the thalamus varies based on the occurrence of an SO-spindle coupling event. We conducted nine whole-brain PPI analyses with nilearn to investigate interactions based in the hippocampus, thalamus or mPFC during SO-spindle coupling, isolated SO and isolated spindle, utilizing the same ROIs. Both the main effects of physiological and physiological condition, as well as their interaction were modelled in the design matrix for PPI analysis. The main effect of physiological term produced functional connectivity maps of the entire brain in relation to the specified ROI (Fig. 4b-d). While the interaction term asked for the change of functional connectivity along with the SO-spindle coupling (Fig. 4e).

Data and materials availability

The data that support the findings can be provided by the corresponding author Y.L. pending scientific review and a completed material transfer agreement. Requests for the data should be submitted to: https://doi.org/10.5061/dryad.2fqz612x0. The analysis code will be publicly available on https://gitlab.com/liu_lab/eeg-fmri-sleep upon publication.

Supplementary figures and tables

Fig. S1

Fig. S2

Fig. S3

Fig. S4

Fig. S5

Fig. S6

Fig. S7

Table S1

Table S2

Table S3

Table S4

Table S5

Table S6

Table S7

Acknowledgements

Author contributions: Conceptualization, Y.L., H.W., J.G.; Investigation, H.W., Q.Z., J.Z., Y.L., J.G.; Writing – Original Draft, H.W., Y.L.; Writing – Review & Editing, H.W., Y.L., J.Z., Q.Z., J.G. Funding: This study is supported by the National Science and Technology Innovation 2030 Major Program (2022ZD0205500), the National Natural Science Foundation of China (32271093), the Beijing Municipal Science and Technology Commission (Z230010, L222033), Beijing United Imaging Research Institute of Intelligent Imaging Foundation (CRIBJZD202101), the Open Research Fund of the National Center for Protein Sciences at Peking University, and the Fundamental Research Funds for the Central Universities.