Discovery of key whole-brain transitions and dynamics during human wakefulness and non-REM sleep.

The modern understanding of sleep is based on the classification of sleep into stages defined by their electroencephalography (EEG) signatures, but the underlying brain dynamics remain unclear. Here we aimed to move significantly beyond the current state-of-the-art description of sleep, and in particular to characterise the spatiotemporal complexity of whole-brain networks and state transitions during sleep. In order to obtain the most unbiased estimate of how whole-brain network states evolve through the human sleep cycle, we used a Markovian data-driven analysis of continuous neuroimaging data from 57 healthy participants falling asleep during simultaneous functional magnetic resonance imaging (fMRI) and EEG. This Hidden Markov Model (HMM) facilitated discovery of the dynamic choreography between different whole-brain networks across the wake-non-REM sleep cycle. Notably, our results reveal key trajectories to switch within and between EEG-based sleep stages, while highlighting the heterogeneities of stage N1 sleep and wakefulness before and after sleep.

This led to an increase in the correlation and specificity values of the HMM states compared to when the raw non-convolved occurrences of spindles and K-complexes were used, but overall did not change the differences on these scores between the states. HMM states 3 and 6 still accounted for the majority of the HRF-convolved sleep spindles and K-complexes. Within HMM states 3 and 6 we did not find marked differences in their relationship to the graphoelements. Both with and without the HRF convolution HMM state 3 tended to show higher correlation and specificity values than HMM state 6 for both spindles and K-complexes, although these differences did not survive correction for multiple comparisons, as shown in Supplementary Figures 20 and 21. The theoretical scenario that the two HMM states (3 and 6), accounting for the majority of N2 sleep, might represent direct reflections of different sleep graphoelements does thus not seem to be supported.
Rather, both of these HMM states included each of their share of both spindles and K-complexes. In light of previous event-related demonstrations of robust effects of sleep graphoelements on the BOLD signal 18,19,20,21,22 , as well as the similarity between these event-related patterns and the mean activation maps of HMM states 3 and 6, it seems unlikely that spindles and K-complexes did not influence the HMM state description, however this did not result in any HMM state coding exclusively for one or the other graphoelement. There is of course the possibility that HMM states 3 and 6 are suggesting a categorisation of spindles and K-complexes beyond what the scalp EEG is able to resolve, or at least beyond the classical interpretation of these graphoelements. Future mapping of the spectral and spatial properties of these graphoelements at higher resolutions than the AASM criteria, that were used here, could bring more insights in this regard, as could a combination or comparison with intracortical evidence 23,24,25 .
There is growing evidence that neuroimaging timecourses contain long-range temporal dependencies 26,27,28 , i.e. they are non-Markovian 29 . The HMM used here follows the Markovian assumption in the sense that the probability of a state transition at a given time point depends only on the state that is active at the preceding time point, and hence it does not parametrically model long-range temporal dependencies. Importantly, however, it does not preclude them either. This means that the HMM state timecourses can in fact exhibit non-Markovian dynamics and long-term dependencies; see e.g. 30 . Notably, our finding of HMM states grouping into modules of transitions represents an analysis that goes beyond Markovianity, and demonstrates non-Markovian dynamics (i.e. long-term dependencies) at the system level of the HMM states. In light of this, our finding that N3 sleep was modelled almost exclusively by a single HMM state, while several states grouped into modules during wakefulness, is in line with the study by Tagliazucchi and colleagues, showing that long-range temporal dependencies in fMRI signals decreases from wakefulness to N3 sleep 16 .

Methodological considerations
For the HMM analysis we chose to make use of the full dataset of 57 participants, when inferring the states. Subsequently we analysed the part of the HMM solution that corresponded to the 18 participants that included all PSG stages. We chose to include as much data as possible for the initial inference in order to maximise the signal-to-noise ratio and amount of evidence for the HMM parameter estimation. Yet, as it is clear from Supplementary Table 1, the full dataset included rather uneven distributions of PSG stages. While PSG stages are more evenly distributed in the subset of 18 participants, and we diligently made sure to normalise the relevant summary measures by number of samples of PSG stages within participants, there still exists a possibility that the HMM could be biased by having more data available from certain PSG stages than others, even if the HMM remained uninformed of the PSG staging. For instance, one could imagine that more data from a certain PSG stage would lead to more HMM states per time being assigned to data from that PSG stage. Such an effect is not immediately present in our results, however, since for instance we found both switching and range of HMM states to be higher in N1 sleep, even though the original full dataset included more than twice as much data from wakefulness (see Supplementary Table 1 and Figure 3d-e).
The number of HMM states was set to 19 in our analyses based on an evaluation of a range of HMM solutions with varying numbers of states (see Methods). It is important to note, however, that it is difficult to determine a 'correct' number of states, when decomposing continuous recordings of brain activity. The recent study by Haimovici and colleagues aimed to identify individual FC states for each PSG stage, and thus chose 4 states for their sliding-window analysis 12 . Our aim was to extract as much temporal resolution as possible from the BOLD ROI timecourses. Ideally, a higher number of states should provide more temporal detail, however increasing the number of states above 19 was associated with a higher occurrence of 'sporadic' HMM states, modelling very specific subparts of the data and not generalising across participants (see Methods and Supplementary Note 3). Based on this it is important to emphasise that we do not suggest the number 19 as definitive, but simply a tool to resolve as much temporal information as possible from the current dataset.
It should be noted that the HMM framework was chosen over other methods for extracting dynamic states from multivariate neuroimaging datasets, such as sliding-window clustering 12, 31, 32 , point-process analysis 33 , and co-activation pattern analysis 34,35 (for reviews, see 36,37 ). The employed HMM framework has been successfully applied to resting state data of wakefulness in both MEG 38,39,40 and fMRI 30,41 , and was particularly suitable for our purpose by virtue of its explicit modelling of temporal dynamics, resulting in states that repeat in a predictable way. Although the HMM is not a mechanistic model of brain activity (a limitation shared with the alternative approaches mentioned above) we have shown how the explicitly modelled HMM transition matrix was fundamental to suggest new partitions of dynamic whole-brain states, which future mechanistic frameworks of NREM sleep and wakefulness should take into account 42,43,44,45 .
Sleep is of course a process associated with profound physiological changes, not merely those reflected in brain activity 46 . Despite our use of the RETROICOR method to reduce the effects on the fMRI data of cardiac and respiratory signals (see Methods), it is currently not possible to completely isolate the neural effects of sleep in fMRI. Our results hence share the limitation with other neuroimaging studies of sleep of potentially being influenced by physiological changes not directly linked to brain activity. On the other hand, certain sleep-dependent peripheral changes such as those of the autonomic nervous system will also induce genuine activities in the brain, which in future studies would be important to investigate and with the proper recordings could potentially be evaluated within an HMM framework.
In addition, it should be noted that there could be potentially confounding effects of spatially smoothing the fMRI data, which can create artificial dependencies between regions of interest. In the context of the HMM, which focuses on the aspects of the data that represent more variance, this confound is however likely to be minor.
Another potential caveat of our analyses pertains to the initialisation of the HMM, which is not deterministic. In Supplementary Note 2, we provide a summary analysis showing that the HMM infers consistent states across independent initialisations and splits of data (see Supplementary

Perspectives
Features identified by the HMM could prove to be essential supplements to PSG and other conventional methods when trying to understand phenomena like the subjective perception of sleep 47,48 , mental content during sleep 49,50 , such as the hypnagogic or even hallucinogenic character of sleep onset 51,52 , sleep inertia of the awakening process 53 , sleep-dependent processes related to memory and learning 54 , and disordered sleep, like insomnia 55 . Such studies should explore the theoretical potential of applying the current HMM, parameterised on the present sleep fMRI data, to identify the presence of the same dynamical whole-brain network states and transition modules in data from different cohorts, potentially even at the individual level. This new data could then be linked to behaviour and cognition through sophisticated measures of arousal, such as eyelid-closure 56 , sleep mentation 49 , post-sleep memory-and learning performance 57 , and careful clinical examination of sleep disorders 58,59,60 .
The wake-NREM sleep cycle merely represents a sub-part of a continuum of activities that the brain supports. Other important brain processes should be sought integrated with the presented transition map (Figure 7), most obviously including REM sleep, but also other altered states of consciousness, such as anaesthesia 61,62,63 , the psychedelic experience 64 , and even different contents of consciousness during wakefulness 65 .
Finally, there is scope for an even more detailed examination of sleep within the HMM framework, given that BOLD data is not the most temporally sensitive modality available. Recently developed methods combining the HMM framework with source-reconstructed MEG data could prove capable of providing an even more fine-grained picture of sleep's evolution in whole-brain networks, and allow for an examination of microstructural EEG elements of sleep, such as spindles and K-complexes 38,39,40,66 , as well as EEG-markers of vigilance fluctuations during wakefulness 67 .

Robustness across different parcellations
We chose the AAL over other possible parcellations because it is the most frequently used in previous fMRI studies of FC during NREM sleep 15,16,17,18,19,20 . Alternative parcellations, such as those derived from FC configurations in the data, could be problematic, since FC has been shown to robustly vary across the sleep cycle 9, 10, 11 . Being anatomically defined, the AAL is essentially agnostic to potentially changing FC configurations within the data. In order to make sure that the use of the HMM generalises to different levels of spatial granularity and that the interpretation following from our results were not specific to the use of the AAL atlas, we re-ran the HMM with a different parcellation. While the field of proposed parcellations for large-scale neuroimaging is rapidly expanding 68 , we opted for the Brainnetome atlas, originally published by Fan and colleagues 69 . Unlike many of the most popular parcellation schemes, the Brainnetome is not solely derived from fMRI FC, but also depends on structural connectivity information for its partitioning of the brain volume. As mentioned above, the use of an FC-derived atlas could bias results, since FC has a well-established dependence on vigilance. Another advantage of the Brainnetome is that it, like the AAL, includes sub-cortical regions, which, as shown in the Results and Discussion of the main text, undergo important changes in activity across NREM sleep. Finally, the 246 regions of the Brainnetome atlas compared to the 90 regions of the AAL provides a good test for the robustness of the HMM across different levels of spatial granularity.
We followed exactly the same steps as explained in the Methods section, but extracted ROI timecourses from the Brainnetome atlas instead of the AAL. It became clear that the increase in spatial detail, going from the AAL to the Brainnetome, had an impact on the ability of the HMM to track the sleep scoring. As such, when using 90% of the variance from the PCA on the Brainnetome ROI timecourses (Figure 1), the performance of the HMM, as quantified through MANOVA between the resulting HMM state timecourses and the sleep scoring, was inferior to the original results using the AAL (see Supplementary Figure 22). However, a slightly stronger regularisation of the ROI timecourses, using only 85% of the variance from the PCA, made the results from the Brainnetome highly comparable to the original results using the AAL. At 85% of the variance the HMM on the Brainnetome data performed in a very similar fashion to the HMM on the 90% of the AAL data, in terms of MANOVA and the development of median fractional occupancy across number of HMM states (see Supplementary Figure 22b and c). The difference between using 90% and 85% of the variance was importantly also evident in the number of HMM states that were consistent across participants for a given HMM solution (K = 19, see Supplementary Figure 22d).
For 90% of the variance, only 6 HMM states occurred in more than 25% of the participants, whereas this number increased to 12, when 85% of the variance was used. In Supplementary Figure   23 we have re-constructed Figure 3 of the main text, but with the results using 19 HMM states on 85% of the variance of the Brainnetome data. The results are highly consistent, with individual HMM states showing sensitivity and specificity to different sleep stages, and in terms of the differences in dynamics found between sleep stages. Regarding the spatial configuration of the HMM states resulting from the Brainnetome data, these were also highly consistent with the original HMM states using the AAL. This is illustrated in Supplementary Figures 24 and 25, where we have matched HMM states from the Brainnetome to the original HMM states, based on their specificity profiles to sleep stages and spatial patterns.
Overall, the above analysis shows that increasing the spatial granularity by introducing a different parcellation comes at the cost at decreasing the signal-to-noise ratio on the HMM estimation. However when this is controlled through PCA, results can be brought to convergence.

Robustness across different HMM initialisations
The initialisation of the HMM includes a stochastic element. To make sure that the states inferred by the HMM were not contingent on the initialisation, we ran the HMM with 19 states an additional four times on the full dataset (N = 57), and five times on each of the two half-splits of the data (N = 29 and N = 28). The 19 resulting states of each HMM repetition were matched to the states of the original HMM. Each state of a repetition was thus paired to an original HMM state, based on the similarity between their respective Gaussian distributions. The similarity was estimated using the Bhattacharyya distance 70 , and the matching of states across repetitions were carried out using the Munkres algorithm 71 .
Following the pairing of states, all resulting states were compared in an all-to-all manner, again using the Bhattacharyya distance as a measure of similarity. The resulting matrix [(n dataset + n states + n repetitions ) × (n dataset + n states + n repetitions )] is shown in Supplementary Figure 15a. The common pattern in the data-set-specific sub-matrices indicates that consistent HMM-state distributions were inferred across initialisation repetitions and data-splits. Another main result of this study is presented in the transition map of the HMM states (see Whereas these overall configurations of the HMM transitions were found robust to the chosen number of states, the more fine-grained details of the transition map appeared more variable. The gateway-like quality of a DMN-like configuration of brain activity was thus particularly clear for the originally chosen 19 states.

Relationship between HMM states and sleep graphoelements
In order to determine the effect of micro-structural features in the sleep EEG, (sleep graphoelements) on the HMM states we used information on the occurrence of sleep spindles and K-complexes during the fMRI recordings. The procedure for obtaining this information from the EEG for the present data has previously been described in Jahnke et al. 21 . Briefly, sleep graphoelements were manually identified according to the criteria set out in the AASM guidelines 72 . This included the use of an EEG montage with frontal, central, and occipital electrodes rereferenced to the contra-lateral mastoid electrodes (TP9, TP10). The resulting temporal markings of sleep spindles and K-complexes were re-sampled to the sampling frequency of the fMRI acquisition (TR = 2.08 seconds) and collected in the variables SS-timecourse and KC-timecourse (for illustration of the SS-and KC-timecourse in an example participant, see Supplementary Figure 19).
To account for the delay in the BOLD response, we also created versions of the SS-and KCtimecourses convoluted with the canonical hemodynamic response function (HRF). We used the HRF included in the SPM12 function spm_hrf.m (http://www.fil.ion.ucl.ac.uk/spm/). Specifically, SS-timecourse and KC-timecourse were binary and of the same length as the fMRI data, with ones representing the fMRI samples during which the respective graphoelement occurred, while the HRF-convolved versions were scaled between 0 and 1 with the canonical delays and undershoots (an example of the HRF-convolved timecourses is provided in Supplementary Figure 19c). We evaluated in turn the temporal association of each HMM activity timecourse to both sleep spindles and K-complexes. Three summary measures of association were used: i) Pearson's correlation was computed between each of the HMM state timecourses and the SS-and KC-timecourses within the set of participants that included the given graphoelement (see Supplementary Table 3 Figure 19 for an illustration of the permutation principle).

FC maps of whole-brain network states
The HMM characterised each state by a vector of mean activity and a covariance matrix. These two data statistics may be understood as a dyadic hierarchy of FC information. The mean distribution of a given state gets estimated from the demeaned and standardised timecourses, and hence describes a change away from the grand-average activity level. ROIs that change their activity in the same direction (positive or negative) are thus functionally connected. The covariance matrix of the state then describes the pairwise ROI-to-ROI co-fluctuations within the time periods where the state is active after subtracting the mean activity of the state (i.e. it is the covariance matrix of the residual).
Therefore, the state-wise covariance matrices reflect FC within states, above and beyond the largest global FC trends as accounted for by the mean parameter of the Gaussian distributions. The mean activation maps have been overlaid on brain surfaces in Figure 5  For the WASO-related HMM states 5 and 17 the frontal increases in mean activation seen in The whole-brain network states of the white N1-related module are represented in Supplementary Figure 5a. These states generally exhibited opposite polarities in their mean activation between subcortical areas (thalamus and parts of the basal ganglia) and primary sensory cortical areas (see Figure 6a). In terms of differential FC this subcortico-cortical decoupling was mainly evident for HMM state 1, while increases in connections towards the medial prefrontal cortex and the anterior cingulate was common for all three states. N2 sleep was dominated by HMM states 3 and 6, and the differential FC maps of these wholebrain network states are shown in Supplementary Figure 5b. Relative increases in FC between superior temporal/inferior parietal areas and medial frontal areas were common, and a few increases in thalamo-cortical connections could be related to association of these HMM states to sleep spindles and the thalamo-cortical mechanism behind the generation of these 73,74 (although see Supplementary Note 4 for a closer examination of the relationship between the HMM states and sleep graphoelements).
HMM state 16 accounted for the majority of time spent in N3 sleep, and its differential FC is shown in Supplementary Figure 5c. Interestingly, the decreases in mean activation, exhibited broadly throughout the frontal cortices, were complemented by relative increases in frontal and temporo-frontal connectivity. In the same way that the frontal decreases in mean activation could be explained by slow-wave activity (see main text), the relative increases in frontal FC are also in line with the localisation found in EEG-based source modelling of slow waves 75 .
To make sure that the differential FC maps were not biased by the effect of the HMM states having different baseline mean activation patterns (as modelled by the mean vector of the Gaussian distribution), we produced equivalent maps using the cosine similarity instead of the covariance matrices outputted by the HMM. Unlike Pearson's correlation, the cosine similarity does not demean the time series, and, therefore, the differences in baselines are accounted for. These maps are shown in Supplementary Figures 17 and 18. This analysis yielded 19 matrices of 90 × 90 cosine similarity values. By taking each cosine similarity matrix and subtracting from it the average of the remaining 18 cosine similarity matrices, we obtained maps equivalent to the differential FC maps, which were based on the covariance information modelled directly by the HMM. As may be seen,

Supplementary Table 3. Summary statistics of sleep graphoelements in dataset
The EEG acquired simultaneously with the fMRI was used to identify sleep graphoelements. This scoring information was resampled to the fMRI, such that volumes during which either a sleep spindle or a K-complex occurred were marked. The first column of the table shows the mean number of occurrences of each sleep graphoelement after this re-sampling. The mean value is calculated within the participants that included at least one of the given graphoelement. The number of participants including at least on sleep spindle or K-complex is shown in the first and second row, respectively, of the second column of the table.

Supplementary Figures
Supplementary Figure 1. Summary measures of HMM solutions across a range of model orders.  Error bars indicate the standard errors across participants. Note how the curve stagnates around K ~ 18. Figure 3).a and b corresponds to a and b of Figure 3 respectively, but for the 31 participants that woke up after having reached N2 sleep. Wake after sleep onset (WASO) was defined as polysomnographically estimated wakefulness following N2 sleep, and is represented by the black bars, while wakefulness (W) was defined as periods of PSG-estimated wakefulness prior to N2 sleep. The coloured bars and error bars show the average and standard error, respectively, across participants. To highlight the differences between W and WASO, we have excluded the remaining sleep stages from this plot. Horizontal lines show significant differences with p-values < 0.01, as evaluated with paired t-tests and permutation testing. Note how HMM states forming the black module in the transition map of Figure 4, showed higher specificity for WASO compared to W. To show that this was a feature of the ROI timecourses and not an effect of the HMM, we have included histograms of the super-diagonal elements of the static FC matrices computed directly from the fMRI data within the PSG stages. c The ROIs of the FC matrices in a have been reordered such that the first and third quadrants include cross-hemispheric connections, while inter-hemispheric connections are shown in the second and fourth quadrants. Figure  5). Differential FC maps were computed by taking an FC matrix of a given HMM state and subtracting the average of the FC matrices of the remaining states. Maps show the 1% most negative and 1% most positive weights of each state a Differential FC maps of wakefulness prior to sleep. Note the relative increases in connections between occipital and temporal areas, as well as the decreases between supramarginal gyrus and the posterior cingulate cortex. b Differential FC maps of wakefulness after sleep onset (WASO). These states were in general characterised by relative increases in frontal connections. Figure 6). a Differential FC maps of HMM states related to N1 sleep. Note the decreased FC between subcortical and temporal/parietal areas of HMM state 1. b Differential FC maps of N2-related HMM states. A common trait of these three states was increased FC in temporal/inferior parietal areas and medial frontal. c N3-related differential FC maps. Note the relative increases in frontal FC. Maps show the 1% most negative and 1% most positive weights of each state

Supplementary Figure 6. Summary of inconsistent HMM states across participants. a
The graph shows the proportion of HMM states that were present in all participants as a function of model order. Included here are the 18 participants that reached all four PSG stages during their recording session. As the model order increased the probability of finding all HMM states in all participants dropped. Error bars show the standard error across participants. b The solution presented in the main text (K = 19) is summarised in terms of the presence of the individual HMM states. The bars show the proportion of participant in which the HMM states were not present. As evident, six HMM states (7, 9, 11, 12, 14, and 19) were present in less than 60 % of the participants. We call these 'sporadic' HMM states, as they did not model data traits consistent across participants. Figure 3). a Fractional occupancies of each of the 15 HMM states computed within the four PSG stages corresponded to the PSG-sensitivity of the whole-brain network states. In d and e error bars represent standard error across participants and significant differences between PSG stages are denoted by stars: * p < 0.05, ** p < 0.01, *** p < 0.001.

Supplementary Figure 11. K = 15 Investigating transitions between whole-brain network states.
(supplement to Figure 4). a The figure shows the 15 × 15 transition probability matrix of the HMM states calculated for the 18 participants that included all four PSG stages in their respective scanning session. This quantifies the likelihood of transitioning from any given state to any other state, giving each matrix entry: probability of (departure state, destination state). b A few HMM states were 'sporadic' and did not occur consistently across participants. HMM states not occurring in more than 25% of the participants were excluded. c The strongest transitions of the consistent HMM states were partitioned through a modularity analysis, and reorganised in a matrix according to the four resulting modules. d The transitions shown in c are presented as a transition map with each state depicted as a pie plot expressing its specificity for each of the four PSG stages. Arrows show the direction of the transitions with thickness proportional to the transition probability. Note the similar overall structure of the transition map to the one presented in Figure 4 of the main text. Unlike the other numbers of states tested, K = 15 only yielded 3 modules, suggesting that the white and the blue module of Figure 4 have merged together in one. e Pie chart showing the total proportion PSG stages within the 18 participants.

Supplementary Figure 12. K = 17 Investigating transitions between whole-brain network states.
(supplement to Figure 4). a The figure shows the 17 × 17 transition probability matrix of the HMM states calculated for the 18 participants that included all four PSG stages in their respective scanning session. This quantifies the likelihood of transitioning from any given state to any other state, giving each matrix entry: probability of (departure state, destination state). b A few HMM states were 'sporadic' and did not occur consistently across participants. HMM states not occurring in more than 25% of the participants were excluded. c The strongest transitions of the consistent HMM states were partitioned through a modularity analysis, and reorganised in a matrix according to the four resulting modules. d The transitions shown in c are presented as a transition map with each state depicted as a pie plot expressing its specificity for each of the four PSG stages. Arrows show the direction of the transitions with thickness proportional to the transition probability. Note the similar overall structure of the transition map to the one presented in Figure 4 of the main text, including 4 separated transition modules. e Pie chart showing the total proportion PSG stages within the 18 participants.

Supplementary Figure 13. K = 21 Investigating transitions between whole-brain network states.
(supplement to Figure 4). a The figure shows the 21 × 21 transition probability matrix of the HMM states calculated for the 18 participants that included all four PSG stages in their respective scanning session. This quantifies the likelihood of transitioning from any given state to any other state, giving each matrix entry: probability of (departure state, destination state). b A few HMM states were 'sporadic' and did not occur consistently across participants. HMM states not occurring in more than 25% of the participants were excluded. c The strongest transitions of the consistent HMM states were partitioned through a modularity analysis, and reorganised in a matrix according to the four resulting modules. d The transitions shown in c are presented as a transition map with each state depicted as a pie plot expressing its specificity for each of the four PSG stages. Arrows show the direction of the transitions with thickness proportional to the transition probability. Note the similar overall structure of the transition map to the one presented in Figure 4 of the main text, including 4 separated transition modules. e Pie chart showing the total proportion PSG stages within the 18 participants.

Supplementary Figure 14. K = 23 Investigating transitions between whole-brain network states.
(supplement to Figure 4). a The figure shows the 23 × 23 transition probability matrix of the HMM states calculated for the 18 participants that included all four PSG stages in their respective scanning session. This quantifies the likelihood of transitioning from any given state to any other state, giving each matrix entry: probability of (departure state, destination state). b A few HMM states were 'sporadic' and did not occur consistently across participants. HMM states not occurring in more than 25% of the participants were excluded. c The strongest transitions of the consistent HMM states were partitioned through a modularity analysis, and reorganised in a matrix according to the four resulting modules. d The transitions shown in c are presented as a transition map with each state depicted as a pie plot expressing its specificity for each of the four PSG stages. Arrows show the direction of the transitions with thickness proportional to the transition probability. Note the similar overall structure of the transition map to the one presented in Figure 4 of the main text, including 4 separated transition modules. e Pie chart showing the total proportion PSG stages within the 18 participants.
Supplementary Figure 15. Robustness across different random initialisation of the HMM. a To make sure that the states inferred by the HMM were not contingent on the initialisation, we ran the HMM with 19 states an additional four times on the full dataset ('Full data set', N = 57), and five times on each of the two half-splits of the data ('Half-split 1', N = 29, and 'Half-split 2' N = 28). The 19 HMM states from a given repetition were matched to the HMM states of the original solution, using the Munkres algorithm, based on the intra-solution state-distances (measured using the Bhattacharyya distance between the states' Gaussian distributions). The matrix shows the pairwise Bhattacharyya distances between the Gaussian distributions following the matching of the HMM states from the different repetition runs. As indicated by the labelling above the matrix, the white borders demarcate states from the various data splits. The smaller black squares surround 5 repetitions of each of the 19 HMM states. The consistent appearance within each of the 9 white squares is a sign that the HMM inferred states with consistent Gaussian distributions. b For a pair of HMM states (one original and one from a repetition run) the temporal correspondence was quantified as the ratio between time points of overlap (simultaneous activity or inactivity) and time points of misses. The bar plot shows mean values and error bars show the standard deviations within data-splits. Note how that temporal overlaps outweighed misses for all runs of the HMM.
Supplementary Figure 16. Robustness of HMM results without the use of temporal filter (supplement to Figure 3). All plots were computed in the same way as for Figure 3, using an HMM with 19 states. The only difference was that the BOLD data was not temporally filtered. Note in a and b how, similarly to the results of the main text, select HMM states were sensitive and specific for certain PSG stages (although not for N1 sleep). In c the overall mean life time of the HMM states is decreased compared to the results using a lowpass temporal filter (compare with Figure 3). As a consequence, d shows increased switching frequencies within all of the four PSG stages. Importantly, the relative differences between the PSG stages were effectively unchanged (even if the significant difference between 'Wake' and 'N2' was no longer evident). Interestingly, in e, the number of unique HMM states visited per unit time were numerically quite stable with or without the use of temporal filter.  significantly higher with spindles than most of the other HMM states, while no significant difference were found between the two. This was true regardless of HRF convolution of the spindles. HMM states 3 and 6 also accounted for the majority of spindle occurrences, as quantified through their sensitivity. Given the generally low specificity values, it is also clear that no HMM state occurred exclusively during spindles.

Supplementary Figure 21. Relationship between the 19 HMM states and the presence of K-complexes. a
The grey violin plot shows the distribution of Pearson's correlation values computed between the timecourse of each HMM state and the raw timecourse of K-complexes for each of the 57 participants that included Kcomplexes (see Supplementary Table 3 Figure 20. K-complexes correlated higher with the HMM states with high specificity for N2 sleep. HMM states 3 and 6 were thus found to correlate significantly higher with K-complexes than most of other HMM states, while no significant difference were found between the two. HMM states 3 and 6 also accounted for the majority of K-complex occurrences, as quantified through their sensitivity. Given the generally low specificity values, it is also clear that no HMM state occurred exclusively during K-complexes.  Figure 1b for an equivalent plot for the AAL data). Going from 90% of the variance (red) to 85% of the variance (blue) had a significant effect on how well the HMM states related to the PSG scoring. Notice how the red line rarely goes below the zone representing the permuted, random cases, whereas the blue line emulates the original analysis on the AAL data (see curve in Supplementary Figure 1b). c Tracking of the median fractional occupancy of the HMM states across model orders, within the 18 participants that included all four PSG stages. In blue is shown the curve for 40 PC's ~ 85%, while the result using 70 PC's ~ 90% is shown in red (An equivalent plot for the AAL data may be found in Supplementary Figure 1c). The fact that the red line is consistently lower than the blue suggests that using the higher percentage of variance implied a high occurrence of 'sporadic' HMM states that accounted for only small portions of the data. This is also evident in d where the HMM solution using 19 states are shown in more detail, for 40 PC's on the left and for 70 PC's on the right. These plots are equivalent to that of Supplementary Figure 6b, which pertains to the original HMM on the AAL data, and show the percentage of participants that did not include each of the 19 HMM states. Using 40 PC's ~ 85 % of variance produced a result more similar to the original HMM on the AAL, with 12 HMM states being included in more than 25% of the participants, whereas including 70 PC's ~ 90% meant that only 6 HMM states were represented in more than 25% of the participants. Error bars in b and c represent standard error across HMM states within a model order.

Supplementary
Supplementary Figure 23. Robustness of HMM results when using an alternative parcellation (Brainnetome) (Equivalent to Figure 3 but for the HMM run on 40 PC's ~ 85% of the variance of the Brainnetome ROI timecourses with 19 states). a Select HMM states account for the majority of different PSG stages as quantified through their fractional occupancies. b In the same way as for the original analysis on the AAL data there is an overlap between the HMM states with high sensitivity for a given PSG stage and the HMM states with high specificity for the same PSG stage. c HMM states with high specificity for N3 sleep expressed higher mean life times, as was the case for the original analysis. d The relative as well as the absolute values of switching were very similar to those of the original analysis on the AAL data. e Similarly to the switching dynamics in E, the ranges of unique HMM states visited within each PSG stage were similar to the original analysis. Please note that all plots were calculated from the 18 participants that included all PSG stages, and that significant differences between HMM states or PSG stages were calculated in the same way as for Figure 3.
Supplementary Figure 24. Spatial correspondence between HMM states from AAL and Brainnetome in wakefulness-related HMM states. (Brain plots from the Brainnetome data are extracted from the HMM solution with 19 states on the 40 PC's ~ 85% of the variance). To demonstrate the correspondence between the original HMM solution on the AAL data and the HMM solution on the Brainnetome data, the original brain plots of mean activation distributions (from Figure 5) are shown together with brain plots from the HMM on Brainnetome. These have been matched based on visual similarity between the spatial maps and their specificity profiles for PSG stages, represented in pie plots. a The original wake-related HMM states from the AAL together with HMM states from the Brainnetome. Notice the high correspondence not only in spatial distribution but also in the PSG-specificity. The original HMM state 8 appeared to show similarity to two HMM states from the Brainnetome analysis (HMM states 16 and 3). b The original three WASO-related HMM states appeared to have two equivalents in the Brainnetome solution.
sleep) A) The original N1-related HMM states found two equivalents from the Brainnetome HMM states. B) Each of the three original N2-related HMM states had equivalents from the Brainnetome HMM, which was also the case for the original N3-related HMM state shown in C.