Reproducible, data-driven characterization of sleep based on brain dynamics and transitions from whole-night fMRI

HMM brain states

To study brain activity representative of the entire Wake-NREM-REM-Wake sleep cycle, we analyzed data from concurrent whole-brain EEG-fMRI measurements on healthy, non-sleep-deprived participants (n=12, age 24 ± 3.5, 8 female) over two successive, entire nights of sleep (Moehlman et al., 2019). This data was acquired for an independent project that included eight randomly timed acoustical arousals to gauge sleep depth (Moehlman et al., 2019). PSG-based sleep staging was conducted by a sleep technologist, utilizing data from EEG, EMG, ECG, and EOG, following the criteria outlined by the AASM (Berry et al., 2020).

Following data preprocessing (see Methods section for details), the fMRI time courses from voxels were spatially averaged within each of the 300 regions of interest (ROIs). These ROIs encompassed cortical, subcortical, and cerebellar areas from the Seitzman 300-ROI atlas (Seitzman et al., 2019). To ensure consistency and comparability, the ROI time courses were demeaned and variance-normalized for each participant and then concatenated along the temporal dimension. Of note, all 12 participants exhibited at least one complete sleep cycle, encompassing all four sleep stages (N1-3 and REM), during both Night 1 and Night 2 (Moehlman et al., 2019). This uniquely comprehensive dataset provided a robust foundation for our analyses.

The HMM estimated from the Night 2 data encompassed a collection of whole-brain states. Each of these states was characterized as a multivariate Gaussian distribution, incorporating two key components: (i) a mean activation distribution, signifying the average activity levels within each ROI when a state was active, and (ii) a functional connectivity (FC) matrix, representing the temporal co-variations among ROIs while in that state.

Furthermore, the HMM featured a transition probability matrix that detailed the likelihood of transitioning between every pair of states. Each state was accompanied by a state timecourse, delineating the specific time points (defined by the fMRI temporal resolution of 3 seconds) when the state was active. Notably, the HMM was constructed with 21 distinct states and was devoid of any prior knowledge regarding PSG staging during its estimation. For a comprehensive visual representation of the analytical process, please refer to Figure 1 (see the Methods section for a detailed explanation). Also, there is no HMM state that was participant-specific. That is, all 21 HMM states can be found in each participant’s fMRI timecourse.

• Download figure
• Open in new tab

Figure 1.

Whole-brain activity dynamic identified from fMRI sleep recording using a Hidden Markov Model. A. Participants slept inside a scanner from ∼23:00 to ∼07:00 for two consecutive nights, with concurrent EEG-fMRI recording. During each night, the fMRI experiments were intermittently disrupted by either acoustical arousals (8 random arousals) or spontaneous awakenings. Sleep stages and slow wave density were derived from EEG signals alone. B. HMM was trained on the principal components of fMRI signals of Night 2. Then the identified HMM states were generalized to Night 1 fMRI signals. Finally, we studied the state-related variations in fMRI activation, FC patterns, and EEG measures. Notes: EEG, electroencephalographic; TR: repetition time; FC, functional connectivity; ROI, region of interest; PCA, principal component analysis.

HMM states show PSG stage specificity

The 21 brain states (see Figure 2B), identified solely from fMRI, exhibited a mixture of six PSG-based sleep stages: N1, N2, N3, REM, Wake, and an “Undefined” stage for epochs that could not be confidently assigned to one of the four following sleep stages: N1-3 and REM.

• Download figure
• Open in new tab

Figure 2.

PSG-based sleep stages and HMM states for each night. A. Distribution of sleep stages for all 12 participants during Night 2. B. Distribution of sleep stages for 21 HMM states during Night 2. C. Distribution of sleep stages for all 12 participants during Night 1. D. Distribution of sleep stages for 21 HMM states during Night 1. The correlation coefficient between sleep stage distributions of HMM states during Night 2 and those during Night 1 is 0.94, p < 0.0001.

To investigate the relationship between HMM states and PSG-based sleep stages, we adopted a “winner-takes-all” approach that assigned HMM states to the sleep stage where they most frequently occurred. Thirteen of the 21 brain states were most frequently associated with N2 sleep stages. HMM states 8 and 10 predominantly occurred during N3 sleep, while HMM states 6 and 19 were prevalent during REM sleep. HMM state 4 corresponded to the undefined sleep stage, and HMM states 13, 16, and 20 were primarily observed during Wake. Intriguingly, none of the HMM states were predominantly linked to N1 sleep. However, HMM state 11 was active for a comparable duration during both N1 and N2 sleep stages. See Figure 2B.

The temporal characterization of these brain states enabled us to investigate the subtle details of brain dynamics within the traditional PSG-based sleep stages. The average duration, referred to as “Lifetime,” of the HMM states varied from 8.7 to 36 seconds. Specifically, the mean Lifetime in states associated with N2 stages tended to be shorter compared to those linked to N3, REM, and Wake (with exceptions of state 13), as illustrated in Supplementary Figure 1.

Sleep states as modules of HMM state transitions

The use of a data-driven approach empowered us to explore the temporal dynamics of HMM states, and enabled us to investigate whether the fMRI-driven HMM states reveal novel dimensions of the Wake-NREM-REM-Wake sleep cycle that are hidden from traditional PSG analyses. We examined the transition probabilities among HMM states, identifying modules of HMM states that exhibited more frequent transitions between each other than to other states (Stevner et al., 2019; Vidaurre et al., 2017).

The transition probabilities of HMM brain states were organized into a 21 × 21 transition matrix. To explore the potential clustering of states with prevalent mutual transitions, a modularity analysis was performed on this matrix based solely on transition probabilities (see Methods section for details). As illustrated in Figure 3, this analysis identified five distinct transition modules, encompassing N3-, REM-, Wake-, and two different N2-modules. Importantly, this modularity analysis was conducted independently of PSG-based sleep stages. Interestingly, it revealed a natural clustering of states associated with the same sleep stages. For instance, two HMM states, 6 and 19, both linked to REM sleep, were grouped within the same module.

• Download figure
• Open in new tab

Figure 3.

Results of the modular analysis based solely on transition probability between HMM states. Each row represents the transition probability of the current HMM state (y-axis) to other states (x-axis). Twenty-one HMM states were categorized into five modules (black boxes): from left to right, light-N2 module (states 5, 7, 9, 14, 15, 18, 21), N3 module (states 4, 8, 10), deep-N2 module (states 1, 2, 3, 12, 17), REM module (states 6 and 19), and Wake module (states 11, 13, 16, 20). The pie chart under each state represents the sleep stage distribution for the state.

Twelve N2-related HMM states were divided into two separate modules. The first module is characterized as the light-N2 module, with higher transition probabilities to REM and Wake modules compared to the other module. The second module exhibited low transition probabilities to both the REM and Wake modules and is referred to as the deep-N2 module.

A similar duality was evident within the REM module. HMM state 19 displayed a notably higher transition probability to states in the Wake module compared to HMM state 6.

Within the Wake module, four HMM states were observed. State 11 was found to be linked to both N1 and N2 sleep stages, while the other three states (13, 16, and 20) were associated with the Wake stages. Further investigation revealed that state 13 typically occurred later in the night and later within an MRI run (see Supplementary Figure 2 E&F), suggesting it represents post-sleep wakefulness, whereas states 16 and 20 were pre-sleep wakefulness. State 13 also showed higher PPG variation, respiratory variation, and heart rates than states 16 and 20 (see Supplementary Figure 2 B-D). This observation was confirmed by the transition probability matrix, that only HMM state 13 has a lower chance of transition into N2- or N3-related states, especially for the states within the light-N2 module, compared to HMM states 16 and 20.

HMM states generalize to Night 1 fMRI data

Next, to test the robustness of our HMM approach, we employed a semi-supervised learning approach (https://github.com/OHBA-analysis/HMM-MAR/wiki/Theory#semi) to predict Night 1 data based on the model trained on Night 2 data. Specifically, we maintained state assignments from Night 2 and applied the model to Night 1. The resulting model indicated that despite having fewer REM and N3 stages during Night 1 (See Figure 2A&C), there was a significant correlation between the sleep stage proportions of the HMM states for Night 1 and those for Night 2 (r = 0.94, p < 0.0001, see Figure 2 B&D). Moreover, the physiological variables displayed similar patterns between Night 1 and Night 2 (see Supplementary Figures 2&3).

fMRI activation and FC patterns of HMM states

To investigate brain activity patterns specific to individual HMM states, we calculated the spatial fMRI activation map and FC pattern of each HMM state relative to the averages over all HMM states. Supplementary Figure 4 showcases the mean fMRI activation for each state, while the associated FC patterns are depicted in Supplementary Figure 5.

For mean fMRI activation, Wake-related HMM state 20 demonstrated the classic opposite activation pattern between the default-mode network (DMN) and its anti-correlated networks (ACNs), see Figure 4 and Supplementary Figure 4. In contrast, during sleep-related HMM states, e.g., states 8 and 10, DMN and FPN showed the same activation direction.

• Download figure
• Open in new tab

Figure 4:

Mean fMRI activation in ROIs within DMN and FPN for each HMM state. Bottom right panel: illustration of DMN (purple) and FPN (green) nodes. Note: DMN: Default Mode Network; FPN: Frontoparietal Network.

For FC patterns, similar anti-correlated patterns were found (see Supplementary Figure 5). In wake-related HMM states 16 and 20, the FCs between DMN and Salience Network (SAL)/ Control Network (CON) were negative, while during N3-related HMM states 8 and 10, these FCs were positive.

As expected, the FC patterns between the Visual Network (VIS) and other sensory networks (Auditory Network, AUD, and lateral/dorsal Somatomotor Network, lSMN/dSMN) were positive during wake-related HMM states but were negative during sleep-related HMM states. One notable exception was HMM state 6 (REM-related), in which VIS had a positive correlation with lSMN and AUD, mirroring those in wake-related states. During REM-related HMM states 6 and 19, the Basal Ganglia (BG) and Thalamus (THAL) had a strong positive correlation with lSMN and AUD.

When we correlated the FC patterns of each state to those of another state, the FC patterns of states that belong to the same module or are related to the same PSG-based sleep stages were highly correlated (similar to the modular results in Figure 3), see Supplementary Figure 6.

Data acquisition and processing

All the data used in this study followed approved human subjects research protocols approved by the National Institutes of Health Combined Neuroscience Institutional Review Board (USA, Protocol Number 16-N-0031), and informed consent was obtained from the participants. Data acquisition was conducted as part of a previously described sleep experiment (Moehlman et al., 2019), encompassing two consecutive nights of concurrent fMRI-EEG data collection while participants slept inside a 3 T Siemens Prisma MRI scanner. To ensure a consistent sleep schedule, participants were instructed to adhere to regular sleep patterns for two weeks before the experiments, and compliance was verified with wearable devices. No sleep deprivation protocols were implemented during the course of the study.

The fMRI data encompassed whole-brain scans consisting of 50 axial slices, captured at a spatial resolution of 2.5 mm (2.5 × 2.5 mm² in-plane), with a 2.0 mm slice thickness and a 0.5 mm slice gap. The data was acquired at a temporal resolution of 3 seconds, employing a 90° flip angle and an echo time of 36 ms. Data acquisition utilized a multi-slice echoplanar imaging approach in an interleaved manner. Simultaneously, EEG data was recorded at a digitization rate of 5 kHz, employing 64 channels to comprehensively cover the scalp. The MR-compatible EEG system used was from Brain Products (Gilching, Germany).

Additionally, concurrent peripheral physiological measures were acquired, including a chest belt to monitor respiratory chest excursion and finger skin photoplethysmography (PPG) to monitor cardiac rate and peripheral vascular volume. These physiological parameters were collected using a Biopac acquisition system with TSD200-MRI and TSD221-MRI transducers, combined with an MP 150 digitizer sampling at 1 kHz, sourced from Biopac in Goleta, CA, USA. To ensure accurate synchronization, data collection for EEG was timed using the 10 MHz clock from the MR instrument. The Biopac device also recorded volume triggers from the MRI scanner to facilitate synchronization of peripheral physiology recordings.

A total of 12 subjects (aged 18–35 years, including 8 females), out of 16 attempts, completed both nights of scanning (from 23:00 to 07:00 roughly). Throughout each night, the fMRI experiments were intermittently disrupted by either acoustically stimulated or spontaneous awakenings. As a consequence, a series of experimental runs was generated, with durations ranging from 5 minutes to 3 hours. Detailed fMRI, EEG, and peripheral physiological measures preprocessing steps can be found elsewhere (Moehlman et al., 2019; Picchioni et al., 2022). Briefly, a tailored version of the ‘afni_proc’ script in AFNI software was used (Cox, 1996), including outlier removal, detrend, RETRIOCOR (Glover et al., 2000), slice timing correction, motion correction, normalization, registration, global signal removal, and censoring (when a motion parameter exceeded 0.3 mm or 0.3 degrees). Previous analysis of the same data indicated that motion during extended sleep scans is comparable to the motion observed in shorter resting-state scans (Moehlman et al., 2019). We also found that motion is lower during deep sleep compared to wake, see Supplemental Results. The EEG signal underwent correction for MRI gradient and cardio-ballistic artifacts and was subsequently down-sampled to a rate of 250 Hz using the Analyzer software (Brain Vision, Morrisville, USA). The process of sleep scoring was carried out using a central electrode in 30-second epochs, in accordance with established criteria with standard filters, and channel references (Berry et al., 2020). ICA cleaning and slow wave auto-detection script were applied to EEG signals (Betta et al., 2021; Mensen et al., 2016; Riedner et al., 2007). Sleep score, slow wave density, and peripheral physiological measures were resampled into a 3-second resolution aligned with the BOLD signal.

HMM overview

In pursuit of a data-driven approach to understanding the brain dynamics in the fMRI signals, we employed an HMM(Vidaurre et al., 2018, 2017) to analyze timecourses extracted from 300 ROIs based on the Seitzman 300-ROI atlas (Seitzman et al., 2019). To prepare the data for analysis, we first standardized the participant-specific sets of 300 ROI timecourses (scaled to a mean of 0, and a standard deviation of 1), which were then concatenated across all participants. This standardization was performed separately for each night. This resulted in a data matrix with dimensions of 300 × (12 × ∼5500) for each night, with approximately 5500 repetition time (TR), excluding breaks between runs and censored TR), accounting for 8 hours of scan time based on a 3-second TR.

The HMM inference process sought to find a sequence of recurring discrete states, each characterized by a distinct statistical arrangement of data. We employed a Gaussian HMM using the Matlab toolbox HMM-MAR (https://github.com/OHBA-analysis/HMM-MAR), where each state was modeled as a multivariate normal distribution encompassing both first-order statistics (mean activity) and second-order statistics (covariance matrix). These state parameters were determined collectively at the group level, while the state timecourses were individually defined for each subject. As a result, the HMM identified periods of quasi-stationary activity, during which the 300 ROI timecourses displayed specific configurations of mean activity and FC.

Given the high spatial dimensionality of fMRI data, we employed Principal Component Analysis (PCA) to reduce the number of parameters in the decomposition process as a common practice. This not only improves the signal-to-noise ratio but also enhances the overall robustness of HMM results (Stevner et al., 2019; Vidaurre et al., 2018, 2017). By selecting the top 13 principal components, we retained 40.7% of the signal variance, resulting in a data matrix with dimensions of 13 × (12 × ∼5500). This matrix was then input into the HMM. For a more detailed overview of the analytical workflow, please refer to Figure 1.

Choice of the number of HMM states

Our analysis involved running the HMM across a range of model orders, specifically spanning from 4 to 25. The assessment of each solution encompassed various summary statistics, with the most pertinent findings illustrated in Supplementary Figure 7.

Supplementary Figure 7A displays the minimum free energy plotted against the HMM model order. This free energy, functioning as a statistical metric, undergoes minimization in the Bayesian optimization process, approximating the model evidence. It encapsulates two crucial factors: the model’s alignment with the data and its complexity, assessed by its deviation from the prior distribution. A lower value of free energy indicates a better model. The first negative peak is observed at K = 21.

To provide insights into the temporal aspects, we defined fractional occupancy as the proportion of time in which an HMM state was active. In Supplementary Figure 7B (and C), we present the evolution of maximum (median) fractional occupancy across HMM states as a function of the model order. We observe a rapid decline in this curve for low values of K, suggesting that, as anticipated, the contribution of each HMM state to the total recording time decreased with an increasing number of states. However, this trend stabilizes at approximately K = 21. This phenomenon is also mirrored in the development of the mean Lifetime of HMM state, which exhibits a similar stabilization pattern at around K = 21, as indicated in Supplementary Figure 7E.

To assess the relationship between the fMRI-based HMM states and PSG-based sleep scoring, we conducted a multivariate analysis of variance (MANOVA). The MATLAB function manova1 was employed to compute Wilk’s Λ, which provides insights into how effectively the K HMM state timecourses can be categorized according to sleep scoring (the lower the better), as depicted in Supplementary Figure 7D. There is a local minimum at K = 21.

Taken together, we chose the model order K = 21 as the number of HMM states. It should be noted that free energy is weighted most among those five model evaluation statistics.

Analysis and visualization of HMM transitions

The transition probability matrix, a fundamental element explicitly modeled by HMM, exhibited a discernible structure characterized by subnetworks of HMM states that displayed more frequent transitions among themselves than to states external to their respective subnetworks. Essentially, this transition matrix could be viewed as a directed graph marked by a modular organization. This characteristic was effectively demonstrated by applying the transition matrix (depicted in Figure 3) to a modularity analysis. This modular analysis was performed using MATLAB functions sourced from the Brain Connectivity Toolbox (https://sites.google.com/site/bctnet/Home) (Rubinov and Sporns, 2010), which relies on Newman’s spectral community detection method (Leicht and Newman, 2008).

Visualizing mean fMRI activation maps and FC patterns of HMM states

The mean distributions and covariance matrices specific to each state were subsequently projected back onto the MNI space utilizing the mixing matrix derived from the PCA. We generated mean fMRI activation maps and FC patterns for every HMM state relative to the baseline averaged over all HMM states. For FC patterns, within- or between-network connectivities were calculated as the average Fisher-transformed functional connectivity between each pair of ROIs within or between networks. For visualization purposes, we grouped 300 ROIs into 14 networks based on the Seitzman Atlas (Seitzman et al., 2019). In addition, we assigned subcortical and cerebellar regions to the additional four Networks: Posterior Hippocampus (pHIP, anterior hippocampus is included in MTL network), Basal Ganglia (BG), Thalamus (THAL), and Cerebellum (CB). Hence, a total of 18 networks were used.

Visualizing state timecourse of HMM states and its associations with PSG stages, PPG amplitude, and Respiratory signals

Two example runs have been shown in Supplementary Figures 8 and 9. These two examples showed how HMM state timecourse (top panel) contained fine-grained information compared to the traditional PSG-based sleep stages (second panel) and also associated with PPG (third panel) and Respiratory signals (last panel).

See Supplementary Results for EEG spectral features across HMM states.