Introduction

Given the significant number of people experiencing sleep issues in modern society, there is a growing need for a better understanding of human sleep and its function (Lancet, 2022). Sleep is characterized by relative stationary states, each believed to serve specific functions. To characterize these states, human sleep research has historically classified sleep into a set of stages using polysomnography (PSG) (Berry et al., 2020; Rechtschaffen and Kales, 1968), which combines electroencephalography (EEG) measures of brain activity with several physiological measures. These sleep stages include the progressively deeper sleep stages of N1, N2, and N3 non-rapid eye movement (NREM), as well as REM. Stages are characterized by different patterns of cortical excitability, as a result of varying levels of modulatory neurotransmitters (Jones, 2020). Across a full night of sleep, these stages cyclically alternate, with REM sleep typically occurring 90 minutes after falling asleep and becoming longer as the night progresses. This cycling is thought to be related to homeostasis and memory consolidation (Diekelmann and Born, 2010; Strauss et al., 2022). Neuroimaging studies using techniques such as Positron Emission Tomography (PET) and functional MRI (fMRI) have identified unique activity patterns for each PSG stage, contributing to our understanding of sleep’s functional role (Braun et al., 1997; Damaraju et al., 2020; Picchioni et al., 2013; Rué-Queralt et al., 2021; Tagliazucchi and Laufs, 2014; Tagliazucchi and van Someren, 2017; Zhou et al., 2019).

While these PSG-guided neuroimaging studies provided new information about sleep function, our understanding of brain dynamics is limited by the low temporal resolution of PSG-based sleep scoring rules (i.e., 30-second epochs), low spatial resolution (i.e., limited EEG channels on the scalp), and the subjective visual inspection rules (Decat et al., 2022; Himanen and Hasan, 2000; Lambert and Peter-Derex, 2023). Alternative to PSG-based sleep staging, applying an unsupervised learning method, the Hidden Markov Model (HMM) (Stevner et al., 2019; Vidaurre et al., 2017), to sleep fMRI data can objectively model the time series of sleep and infer sleep brain states that recur at different points during sleep. A recent study demonstrated promising results in capturing NREM sleep transitions by applying HMM to relatively short bouts of sleep (<1 hour) fMRI data (Stevner et al., 2019). However, because the REM stage typically occurs 90 minutes after falling asleep and lasts progressively longer over time, capturing brain dynamics associated with sleep cycling requires whole-night data.

In addition, given that studies on sleep stage transitions have shown promising results in diagnoses of various sleep disorders, including narcolepsy (Christensen et al., 2015), chronic fatigue syndrome (Kishi et al., 2011), and insomnia (Wei et al., 2017), it is of great interest to establish an objective and reliable measurement of brain states transitions within and between PSG sleep stages. To achieve this goal, we applied HMM to a unique and extensive dataset of EEG-fMRI concurrent recordings acquired over 8 hours of sleep each night for two consecutive nights (Moehlman et al., 2019). This analysis revealed 21 unique brain states, surpassing the number of PSG-defined sleep stages. For potential application in clinical settings, we tested whether our HMM model trained using night 2 data can predict night 1 data. As it turned out, the identified brain states were highly consistent between night 1 and night 2. Furthermore, analyzing the transition probabilities between HMM states revealed a significant subdivision within N2 and within REM sleep stages. This data-driven, PSG-blind analysis of fMRI data provides reproducible brain states and their transition probabilities, potentially serving as a biomarker of sleep transitions in both normal and clinical settings.

Results

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.

Figure 1.

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.

Figure 2.

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.

Figure 3.

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 state, we calculated the spatial fMRI activation map and FC pattern of each HMM states 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.

Figure 4:

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.

Discussion

By applying an unsupervised learning method to Night 2 of two-night fMRI sleep recordings, we deduced 21 HMM states and their transition probabilities, independently of PSG-defined sleep stages. The identified HMM states showed excellent reproducibility to Night 1 data in a semi-supervised manner, a feat not previously demonstrated. Moreover, through modular analysis focused solely on transition probabilities, a duality within REM-related and N2-related HMM states was found. These findings offer unique new information about brain sleep states and their transitions that extend beyond previous PSG-based research, as well as fMRI research without whole-night recordings.

Our work addressed well-known shortcomings of PSG-based sleep staging (Abeysuriya and Robinson, 2016; Decat et al., 2022) by integrating insights from whole-brain fMRI recordings. First and foremost, while traditional PSG-based sleep staging is based on 30-second epochs, HMM analysis allows for a state-specific duration as short as the fMRI temporal resolution (here 3 seconds). On average, the duration of HMM states is 12 seconds (see Supplementary Figure 1), suggesting a more detailed characterization of brain states compared to PSG-based sleep stage analysis. Secondly, in terms of spatial resolution, the functional atlas used in our study encompassed 300 ROIs, offering a more detailed view of activation patterns across the entire brain, including subcortical and cerebellar regions that are ignored in PSG-based sleep staging. Third, our approach is mostly automated and objective, eliminating concerns related to inter-rater reliability issues and human error (Lambert and Peter-Derex, 2023; Lee et al., 2022; Rosenberg and Van, 2013). Lastly, identifying transitions between sleep stages can pose challenges when relying solely on PSG data. In contrast, the HMM is explicitly designed to model these transitions between states, providing a better understanding of the dynamic shifts that occur throughout the sleep cycle, especially when the sleep stages transition is not linear from Wake to NREM to REM in the second half of night.

Previous research suggests that the analysis of sleep data at a finer temporal resolution than PSG-based sleep staging may be valuable. For example, distinct and recurring states of waking brain activity may be as brief as 100 ms during wake (Baker et al., 2014; Koenig et al., 2005). In mice, rapid (seconds-scale) fluctuations in brain-wide neuronal spiking activity has been reported during states of low alertness, attributed to fluctuation in adrenergic and cholinergic neuromodulation from basal forebrain and locus coeruleus (Aston-Jones and Bloom, 1981; Collins et al., 2023; Kjaerby et al., 2022; Osorio-Forero et al., 2023). Capturing second-scale changes of brain states with the analysis approach employed in the current study may therefore allow a more comprehensive investigation of the functional roles of sleep and shed light on the mechanisms by which these roles are accomplished.

The modular analysis, which solely relied on transition probabilities between states, uncovered a significant discovery. This analysis clustered HMM states into modules closely associated with PSG-defined sleep stages. This suggests that transition probabilities contain essential information about sleep states and PSG-based sleep stages. For instance, a module predominantly linked to the N3 stage consisted of two N3-related states and one undefined state. Importantly, all three states also exhibited the highest slow-wave density among all states (see Supplementary Figure 2A).

Within the Wake module, there were four HMM states, each representing pre-sleep wake (states 16 and 20), post-sleep wake (state 13), and N1-2 (state 11). The absence of a dedicated module/state representing the N1 stage is unsurprising, considering that N1 does not distinctly manifest as a well-defined sleep stage (Carskadon and Dement, 2011) and it has the lowest inter-rater reliability (0.24 vs 0.76 overall) among all the PSG-defined sleep stages (Lee et al., 2022).

Two modules were associated with the N2 stage. One of these termed the “Deep-N2 module,” exhibited a low likelihood of transitioning to REM and Wake while showing a slightly higher probability of transitioning to N3-related states when compared to the other module, referred to as the “Light-N2 module”. This finding aligns with previous studies (Brandenberger et al., 2005; Decat et al., 2022), which separated the N2 stage into a quiet type (before the transition into N3 stages, which resembles the Deep-N2 module in the current study) and an active type (preceding the transition to REM, related to the Light-N2 module).

The two REM-related states (6 and 19) within the REM module were notably different in several aspects. First, state 19 displayed a higher propensity for transitioning to the Wake module in contrast to state 6. Second, state 6 tended to occur towards the end of sleep and also late within the fMRI run (see Supplementary Figure 2 E&F). Third, the FC patterns between sensory networks (Vis, lSMN, AUD) and other cortical networks, i.e., DMN and MTL, were different in those two states. For example, in state 19, the MTL was positively connected to lSMN/AUD, while exhibiting negative connections with those sensory networks in state 6 (see Supplementary Figure 5). These differences suggest an alignment of the HMM REM states along the previously defined microstates of REM, i.e., “phasic” and “tonic” episodes (Simor et al., 2020). Tonic REM is thought to be an intermediate state between wakefulness and phasic REM and is associated with a higher environmental awareness. Phasic REM occurs more often at the end of the night and is associated with a higher level of brain activity (Simor et al., 2020). Taken together, HMM state 19 might represent tonic REM given the high transition probability to Wake-related HMM states, while HMM state 6 might be related to phasic REM with higher FC and occurring later in the night.

In terms of both BOLD activation and FC patterns, a notable divergence between N3-related states and Wake-related states is observed in the interaction between DMN and its ACNs (SAL/CON/FPN, etc.). It is plausible that the degree of correlation or anticorrelation between DMN and its ACNs is a pivotal factor influencing the transitions from wakefulness to light sleep and, subsequently, to deep sleep. The SAL is considered crucial for cognitive control, as it handles the perception and response to homeostatic demands (Menon, 2011; Peters et al., 2016; Seeley, 2019). It further acts as a mediator for dynamic interactions among other prominent large-scale brain networks engaged in externally focused attention (FPN) and internally directed self-referential cognitive processes (DMN). It is plausible that during sleep, the mediating function of the SAL is temporarily suspended to allow for its restoration. Recent findings have indicated that disruptions in SAL connections were observed following one night of sleep deprivation (Fang et al., 2015) or in individuals with insomnia disorder (Cheng et al., 2022; Li et al., 2022; Wei et al., 2020).

There are a few limitations worth mentioning. Firstly, we made an arbitrary selection of 13 principal components for PCA, accounting for 40.7% of the total variance. While this percentage of explained variance may seem low, it was a necessary step to stabilize the fitting of the Hidden Markov Model (HMM) in the current study. Notably, the trained HMM demonstrated generalization to Night 1 data, validating the chosen principal components as they encompass sufficient information about the fMRI signals. Secondly, while our study involved a relatively small number of participants (12), it included a large amount of fMRI data (∼16 hours scan) per participant. While the HMM trained on data from 12 participants was robust, future research may benefit from utilizing a larger number of participants to test the out-of-sample generalization of the current results.

To summarize, we demonstrated how a data-driven analysis of an extensive sleep fMRI dataset can reproducibly characterize the full pattern of arousal state changes that recur during a whole night’s sleep. The findings underscore the advantages of the whole-night fMRI data, over the traditional PSG sleep staging and previous fMRI sleep studies, in achieving a fine-grained characterization of brain sleep states and their transitions. The successful generalization of our approach trained on Night 2 to Night 1 data shows its robustness, reliability, and objectivity across multiple nights. Our exploration of transitions between HMM states unveiled modules closely linked to distinct sleep stages, revealing a duality within N2-related modules that further dissects N2 stages into ‘light’ and ‘deep’ N2 modules. We identified a duality with REM-related HMM states, which resembles the ‘phasic’ versus ‘tonic’ REM. Additionally, we separated pre-sleep from post-sleep Wake states. Analysis of brain activation and FC patterns of HMM states indicated that the connections between DMN and ACNs, especially SAL, may play a critical role in the transition from wake to light sleep and subsequently to deep sleep. Collectively, this enriched comprehension of brain dynamics during nocturnal sleep holds the potential for identifying promising biomarkers associated with sleep disorders that significantly impact sleep-stage transitions.

Methods

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). 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 on 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, which were then concatenated across all participants. 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. There two examples showed that how HMM state timecourse (top panel) contained fine-grained information compared to the tranditiaonal PSG-based sleep stages (second panel) and also associated with PPG (third panel) and Respiratory signals (last panel).

Data Availability

The datasets analyzed for this study are available from the last author upon reasonable request and in accordance with institutional policies. The codes are available at https://github.com/nilsyang/Codes.

Acknowledgements

This research was supported by the Intramural Research Program of the NIH, NINDS. Sussan Guttman, Steve Newman, and Tina T. Liu are acknowledged for advice and assistance.

Author Contributions

F.N.Y. conceptualized the study, analyzed the data, generated figures, and wrote the original draft. D.P. collected the data and edited the manuscript. J.A.D., Y.W., and P.V. participated in experiment setups and edited the manuscript. J.H.D. supervised the study and edited the manuscript.

Financial Disclosure

The authors declare no competing financial interests.

Non-financial Disclosure

The authors declare no competing non-financial interests.

Data and code availability statement

The datasets analyzed for this study are available from the last author upon reasonable request and in accordance with institutional policies. The codes are available at https://github.com/nilsyang/Codes.