Introduction

During sleep, the human brain engages in highly synchronized activity driven by distinct rhythms. In non-rapid eye movement (NREM) sleep, the emergence and coordination of slow oscillations (SOs) and sleep spindles are key features (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 their coordination remain largely unexplored, and their functional implications are even less understood.

SOs (∼1 Hz) represent shifts 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 spontaneously 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 oscillations in the 11–16 Hz range (Fernandez & L üthi, 2020; Mak-McCully et al., 2017). These spindles align with the excitable UP-states of SOs and synchronize with hippocampal ripples (Buzsáki, 2015; Helfrich et al., 2019; Joo & Frank, 2018; Ngo, Fell, & Staresina, 2020). Rodent studies using optogenetics have shown that only thalamic spindles phase-locked with the UP-state of cortical SOs enhance memory consolidation, while those out of phase 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 (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). 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.

Studying these processes in the human brain is challenging because it requires tracking both the transient nature of sleep rhythms and their widespread brain activation simultaneously. Recent advances in simultaneous EEG-fMRI techniques, which allow researchers to model temporal events captured by EEG and explain the associated brain activity observed through fMRI, offer promising opportunities to capture these brain-wide activations (Huang et al., 2024). In the current study, while directly capturing hippocampal ripples with scalp EEG or fMRI is difficult, we expect to observe hippocampal activation in fMRI whenever SOs-spindles coupling is detected by EEG, if SOs-spindles-ripples triple coupling occurs during human NREM sleep.

Another advantage of identifying sleep rhythm-related brain activation is the potential to infer their potential cognitive functions through open-ended cognitive state decoding (Yarkoni et al., 2011). 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, this approach provides limited insight into the neural mechanisms underlying natural sleep rhythms.

By using cognitive state decoding, we can infer the functional implications of natural sleep rhythms by comparing their brain activation patterns with known task-related patterns (Margulies et al., 2016; Yarkoni et al., 2011). This method offers a promising way to bridge the gap between observed brain activity during sleep and its functional significance, providing new insights into the relationship between sleep and cognition.

In this study, we investigated brain-wide activation during sleep using simultaneous EEG and fMRI in 107 participants, combined with cognitive state decoding. Our findings revealed significant coupling between SOs and spindles during NREM sleep, particularly at the transition to UP-state of SOs. This coupling was associated with increased activation in the thalamus and hippocampus, and a brain-wide activation pattern resembling episodic memory processing, yet distinct from task-related activations. Furthermore, SO-spindle coupling was linked to increased functional connectivity from the hippocampus to the thalamus and from the thalamus to the mPFC.

These findings suggest that the thalamus plays a crucial role in coordinating hippocampal-cortical dialogue during human sleep, highlighting the importance of intrinsic brain-wide activation when sleep rhythms align. Understanding these mechanisms paves the way for developing neuromodulation interventions to enhance cognitive function during sleep and treat 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). 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).

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

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), 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 during the transition to the UP state 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. S2), 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). 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). For more detailed activation patterns, see Table S1-S3.

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. 3e, right), SOs were linked to patterns associated with declarative memory (z = 6.31, p < 1e-4) and episodic memory (z = 5.08, p < 1e-4), while spindles were most similar to activation patterns related to working memory (z = 4.08, p < 1e-4). SO-spindle coupling aligned most closely with activation patterns related to episodic memory (z = 3.66, p = 0.0003) and declarative memory (z = 2.21, p = 0.0271).

For negative (deactivation) patterns (Fig. 3e, left), SOs were linked to patterns related to declarative memory (z = 7.61, p < 1e-4), while spindles were most similar to task-related activation patterns (z = 7.05, p < 1e-4). SO-spindle coupling showed the strongest similarity to patterns related to task-related activation (z = 9.40, p < 1e-4) and working memory (z = 9.15, p < 1e-4). These findings suggest that brain regions involved in internal memory processing are the most likely to reactivate during SO-spindle coupling events, while task-related activations, such as working memory, are the least likely to occur.

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 hippocampus -> thalamus FC (ROI analysis, t₍₁₀₆₎ = 1.86, p = 0.0328) and thalamus -> mPFC FC (ROI analysis, t₍₁₀₆₎ = 1.98, p = 0.0251) significantly increased during SO-spindle coupling, with no significant changes in all other FC pathways (Fig. 4e). These findings suggest that the thalamus, likely through thalamic spindles, coordinate hippocampal-cortical communication during human NREM sleep.

Discussion

Using simultaneous EEG and fMRI during nocturnal sleep in 107 subjects, we identified robust coupling between sleep spindles and the (transition to) UP-state of SOs during N2 and N3 stages of NREM sleep. This coupling was linked to distinct neural activation patterns, including significant thalamic activation during spindle occurrences, reduced activity in the DMN during SOs, and elevated hippocampal activity during SO-spindle events. Moreover, we observed increased functional connectivity from the hippocampus to the thalamus and from the thalamus to the mPFC during these SO-spindle events. Through cognitive state decoding, we found that brain-wide activation during SO-spindle events most closely resembles that of internal memory processing, while deactivation patterns align with those of working memory and task-related representation. These findings offer insights into the brain mechanisms and associated functional implications during human NREM sleep.

We found that spindle amplitude peaks were nested within the transition to the UP-state of SOs in most subjects (83 out of 107), a coupling likely crucial for sleep-dependent memory consolidation. This is consistent with Schreiner et al. (2021), showing that memory reactivation is stronger when spindles occur closer to the SO UP-state. When modeling the timing of these sleep rhythms, we observed hippocampal activation during SO-spindle events, which is intriguing as it may indicate the presence of triple coupling, including hippocampal ripples. While our scalp EEG was not sensitive enough to detect hippocampal ripples—key markers of memory replay (Buzsáki, 2015)—recent intracranial EEG (iEEG) studies have observed ripples during human sleep, showing their coupling 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.

Our observation of reduced DMN activity during SOs is intriguing. SOs, generated in cortical networks, reach subcortical structures like the hippocampus and alternate between DOWN- and UP-states of synchronized neural activity (Nir et al., 2011). Due to the low temporal resolution of fMRI and the rhythmic pattern of SOs, it is difficult to distinguish the activation patterns of the DOWN- vs. UP-states (see also Fig. S2). It is likely that our findings of diminished DMN activity reflect brain activity during the SO DOWN-state, as this state consistently shows higher amplitude compared to the UP-state within subjects, which is why we modelled the SO trough as its onset in the fMRI analysis. The DMN is typically active during wakeful rest and is associated with self-referential processes like mind-wandering, daydreaming, and task representation (Yeshurun, Nguyen, & Hasson, 2021). Its reduced activity during SOs may signal a shift towards endogenous processes such as memory consolidation. Future research using techniques like iEEG, with both temporal and spatial precision, could clarify this possibility.

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 essential for memory consolidation during sleep, facilitating the transfer of memories to the neocortex for long-term storage (Buzsáki, 2015; Diba & Buzsáki, 2007; Singh, Norman, & Schapiro, 2022). The increased hippocampal and thalamic activity, along with strengthened functional connectivity between the hippocampus, thalamus, and mPFC during SO-spindle events, highlights the hippocampal-thalamic-neocortical information flow during human NREM sleep. This aligns with recent findings that suggest the thalamus orchestrates neocortical oscillations during sleep (Schreiner et al., 2022). These results emphasize the central role of thalamus in sleep-related memory consolidation, guiding the routing of information from the hippocampus to the neocortex (e.g., mPFC), a process critical to understanding memory consolidation during sleep.

Our open-ended cognitive state decoding analysis revealed that brain activation patterns linked to different sleep rhythms were associated with distinct cognitive functions. SOs and SO-spindle coupling were strongly connected to declarative and episodic memory, particularly during coupling events. In contrast, brain deactivation patterns were associated with task representations and working memory, reflecting the reduction in DMN activity (Baldassano, Hasson, & Norman, 2018; Smallwood et al., 2021). It is intriguing to speculate that during SOs, especially when aligned with spindle events, spontaneous replay of past episodic memories may occur, facilitating the consolidation of hippocampal memories into long-term declarative memory in the neocortex. This internal processing contrasts with the brain patterns associated with external tasks, such as working memory. This hypothesis warrants further studies to decode the content of sleep activity, capturing the replay of specific memories during sleep.

Together, our findings provide insights into the neural mechanisms of sleep rhythms, emphasizing their coordination and the roles of the hippocampus and thalamus. They also highlight the importance of the role of thalamus in generating and coordinating sleep rhythms, the implications of reduced DMN activity during SOs, and the temporal dynamics of hippocampal-neocortical communication. Uncovering these brain-wide activation patterns is crucial for understanding the functional implications of sleep activity. These insights will lead to more targeted interventions that improve sleep quality and cognitive performance.

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

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

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 et al., 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, especially the spatial distribution of these sleep rhythms (Massimini et al., 2004; Molle et al., 2011), EEG data from the F3 electrode was 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 underwent a 0.16-1.25 Hz bandpass filtering (Butterworth filter, order 3, bidirectional filtering for zero-phase). After identifying all positive-to-negative zero crossings, potential SOs were determined based on the duration between consecutive zero crossings (ranging from 0.8 s to 3 s). We computed the amplitude range (peak minus trough) for each potential SO. Only those surpassing the 75th percentile in amplitude were marked as SOs (Schreiner et al., 2021). The onsets and durations of these SOs were saved.

Detection of sleep spindles

Data underwent a 12-16 Hz bandpass filtering (Butterworth filter, order 3, bidirectional filtering for zero-phase). The root mean square (RMS) of the filtered signal was calculated with a sliding time window of 200 ms. The amplitude threshold for spindles was set at the 75th percentile of RMS values. RMS values exceeding this threshold and lasting between 0.5 s and 3 s were identified as spindles (Staresina et al., 2015). We saved all their onsets and durations.

Detection of SO-spindle couplings

From the detected SOs and spindles, we identified the peak time of each spindle. Within each SO range, we checked for a spindle peak. If detected, 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). Data were extracted in a window of [-4 s to 4 s] around this point, forming the epoch for each coupling event. All SO-spindle epochs and their associated spindle peak times were saved.

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

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. 3e). 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 three whole-brain PPI analyses with nilearn to investigate interactions based in the hippocampus, thalamus or mPFC during SO-spindle coupling, 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, c, 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://datadryad.org/stash/dataset/doi:10.5061/dryad.2fqz612×0. The analysis code will be publicly available on https://gitlab.com/liu_lab/eeg-fmri-sleep upon publication.

Supplementary materials

Fig. S1

Fig. S2

Table S1

Table S2

Table S3

Additional information

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.