Whole-brain network states identified by HMM

In order to extract the large-scale networks inherent to whole-brain recordings of the wakefulness-NREM sleep cycle, we estimated an HMM on fMRI data from 57 healthy participants (age 23.5 ± 3.3 years, 39 females). Participants were instructed to lie still in the scanner with their eyes closed. Each recording had a duration of 52 min, and was accompanied by acquisitions of EEG, EMG, ECG and EOG, based on which PSG staging was performed by an expert, according to the AASM criteria2 (see Supplementary Table 1). Following preprocessing, the voxel-wise BOLD timecourses were temporally averaged over 90 region-of-interest (ROI) timecourses, using the cortical and subcortical regions of the automated anatomical labelling (AAL) atlas26. ROI timecourses were demeaned and variance-normalised for each participant, and subsequently concatenated across participants along the temporal dimension.

The estimated HMM contained a set of whole-brain network states, each defined as a multivariate Gaussian distribution, including: (i) a mean activation distribution, representing the mean level of activity in each ROI when a state is active; and (ii) an FC matrix, summarising the pairwise temporal co-variations occurring between the ROIs during that state. The HMM also contained a transition probability matrix with the probabilities of transitioning between each pair of states. Each state also had an associated state timecourse describing the points in time (defined by the fMRI sampling, TR = 2.08 s) where the state was active24,27. The HMM was endowed with 19 states, and, crucially, was given no information about the PSG staging for its estimation. An illustration of the analysis workflow is given in Fig. 1 (see Methods for details).

Fig. 1 Dynamic whole-brain networks from fMRI sleep recordings using a Hidden Markov Model. a ROI timecourses were extracted by averaging BOLD signals across voxels within each of the 90 cortical and subcortical AAL areas for each participant. Each ROI timecourse was demeaned and normalised by its standard deviation. b The data were concatenated across participants, and the dimensionality was reduced using PCA (principal component analysis), such that ~90% of the variance of the ROI timecourses was retained. The HMM was run on the PCA timecourses, resulting in K number of states with associated timecourses, each describing the points in time each state is active and inactive. c Each HMM state was characterised by a multivariate Gaussian distribution comprising a covariance matrix, Σ K , and a mean distribution, μ K. The state-specific mean distributions and covariance matrices were back-projected to the MNI space of the AAL by using the mixing matrix, MT from the PCA decomposition, yielding a mean activation map and an FC matrix for each HMM state Full size image

To allow for unbiased within-participant testing, when comparing the HMM output to the PSG scoring we considered the subset of the HMM output that corresponded to the data from the 18 participants that reached all four PSG stages (wakefulness, N1, N2 and N3, see Supplementary Table 1).

Whole-brain network states underlie PSG sleep stages

The 19 whole-brain network states, inferred solely from the fMRI by the HMM, contained most of the temporal information given by the PSG stages that were scored from the EEG independently of the fMRI. The HMM state timecourses and the PSG scoring are plotted together in Fig. 2a for the 18 participants that included all four PSG stages, illustrating how the activity of the different HMM states co-varied with specific PSG stages. We have highlighted the HMM state timecourses of two participants to ease visual inspection, however, the temporal relationships between HMM states and PSG stages were consistent across the group. It can be observed, for example, that HMM state 8 occurred most often during wakefulness, HMM state 3 occurred during N2 sleep and HMM state 16 occurred during N3 sleep.

Fig. 2 State timecourses of whole-brain network states and their association to polysomnography. a The figure shows the 19 HMM state timecourses describing each state’s probability of being active at each sample point of the fMRI sessions in the 18 participants that reached all four PSG stages. Below the HMM state timecourses are shown the independently obtained PSG sleep scoring (based on the simultaneously acquired EEG). The coloured overlay shows periods scored as wakefulness (red), N1 (white), N2 (blue) and N3 (green). The two dashed boxes highlight the HMM state timecourses and PSG scoring of two representative participants. Note, how the majority of HMM timecourses varied with the PSG stages, in highly consistent ways across participants. A few ‘sporadic’ HMM states, occurring mainly in a few participants, are also visible (e.g., states 11 and 12). b Quantifying the multivariate relationship between the HMM states and the PSG scoring, through the use of MANOVA, revealed a hierarchical grouping of the HMM states, in which wakefulness and N1 sleep were separated from N2 sleep, which in turn was separated from N3 Full size image

We quantified the temporal association between the PSG stages and the HMM state timecourses using multivariate analysis of variance (MANOVA). This allowed us to ask if the 19 HMM states were significantly grouped in time by the four PSG stages (for the 18 participants that included all four PSG stages). Through non-parametric testing (see Methods) we confirmed this temporal relationship (p < 0.05, permutation testing, see Supplementary Figure 1b). The MANOVA placed the PSG stages in the space of the HMM state timecourses, resulting in the clustering dendrogram of Fig. 2b, with wakefulness and N1 sleep significantly separated from N2 sleep, which in turn was further separated from N3 sleep.

Whole-brain network states track different PSG stages

Next, we examined the contribution of the individual whole-brain network states to the multivariate relationship, established above, between the HMM and the PSG scoring. We quantified the temporal sensitivity and specificity of the HMM states for each of the PSG stages. For each of the 18 participants that included all four PSG stages, we defined the sensitivity of each HMM state as the proportion of total time spent in a PSG stage, in which this HMM state was active. Specificity was defined as the likelihood of finding each HMM state active during a given PSG stage. We compared the sensitivity and specificity for each PSG stage within each of the HMM states, using paired t tests and a randomisation scheme of the PSG scoring (see Methods). The results are presented in Fig. 3a, b. HMM state 8 occupied a large proportion of PSG-scored wakefulness, i.e., it exhibited high sensitivity for wakefulness (see Fig. 3a). Since this whole-brain network state was significantly more sensitive to wakefulness than to any of the other PSG stages, i.e., it rarely occurred outside of wakefulness, its specificity for wakefulness was also high (see Fig. 3b). This combined sensitivity and specificity for wakefulness was also found for HMM states 10 and 18.

Fig. 3 Sensitivity and specificity of HMM states and dynamics within polysomnography stages. a Fractional occupancies of each of the 19 HMM states computed within the four PSG stages corresponded to the PSG-sensitivity of the whole-brain network states. The coloured bars and error bars show the average and standard error, respectively, across the 18 participants that included all four PSG stages. b PSG-specificity of the HMM states for each of the four PSG stages. Specificity corresponds to the probability of an HMM state occurring within a PSG stage. The bars represent the group average and the error bars the standard error (n = 18). In a and b horizontal lines show significant differences within HMM states, with p values < 0.01 as evaluated through paired t tests and permutation testing. c The mean life times of the 19 HMM states are shown by the bars, representing values averaged across the 18 participants. Error-bars represent the standard error across participants. Each HMM state is coloured according to the probability of finding it within each of the four PSG stages, i.e., their PSG specificity. Note how HMM states with high specificity for N3—and to a lesser extent N2—exhibit longer mean life times. d The dynamics of the HMM transitions were calculated within each of the four PSG stages, in terms of switching frequency (‘Switching’), and e the number of different HMM states visited per time (‘Range of HMM states’). These measures significantly separate the four PSG stages suggestive of a higher dynamical repertoire during wakefulness and N1. In d and e error bars represent standard error across participants and significant differences between PSG stages are denoted by stars: one star: p < 0.05, two stars: p < 0.01 and three stars p < 0.001; all evaluated using paired t tests and permutations. W: wakefulness, N1: N1 sleep, N2: N2 sleep, N3: N3 sleep Full size image

Select whole-brain network states displayed similarly exclusive sensitivity and specificity profiles for N2 (HMM states 3 and 6) and N3 sleep (HMM states 16). Notably, this was not the case for N1 sleep. The whole-brain network states occupying most of N1 sleep, such as HMM states 1, 4 and 15, were not found specific for this PSG stage. Instead these states would also occur with considerable likelihood outside of N1 sleep, although rarely during N3 sleep.

In summary, wakefulness was found to correspond to a collection of whole-brain networks states, while N2 and N3 were characterised by less state-diversity, and dominated by two and one whole-brain states, respectively. In contrast, no single whole-brain states were found specific for N1 sleep, which instead was modelled by a collection of HMM states with mixed PSG profiles.

Changes in whole-brain network dynamics between PSG stages

Having the whole-brain network states temporally defined allowed us to investigate the large-scale brain dynamics of the traditionally defined PSG stages in the 18 participants that reached all PSG stages during their recordings.

In Fig. 3c, the HMM states are represented by a bar plot showing their mean lifetimes, i.e., the average duration of the state visits. The bars have been overlaid with colours depicting the PSG specificity averaged across the corresponding HMM states. HMM states with high specificity for N2 and N3 (HMM states 3, 6 and 16) generally expressed longer mean lifetimes than those related to wakefulness and N1. The mean lifetimes of the HMM states ranged from seconds to tens of seconds.

Figure 3d, e shows two summary measures for the dynamics of the whole-brain network states during the individual PSG stages: (i) the amount of switching defined as the average number of transitions between HMM states during a given PSG stage divided by the total time a participant spent in this PSG stage and (ii) the range of HMM states defined as the number of unique states visited during the given PSG stage divided by the total time a participant spent in this PSG stage. Both measures were estimated for each PSG stage, within each of the 18 participants that included all four PSG stages, and normalised by time. Wakefulness and N1 sleep expressed significantly higher values than N2 and N3. Interestingly, the amount of switching was particularly low for N3 sleep.

In summary, unique state visits per time were few and of long durations during N2 and N3 relative to wakefulness and N1 sleep. Consequently, the switching between and range of HMM states were significantly higher in wakefulness and N1.

Sleep stages as modules of whole-brain network transitions

So far, we have used the traditional PSG stages to organise and evaluate the temporally resolved whole-brain network states. Yet, the data-driven nature of the HMM also allowed us to perform reverse inference, and consider the temporal progression of HMM states, taking this—rather than the PSG staging—as a starting point. This way, we were able to ask if the high-resolution, fMRI-based, HMM suggests new aspects of the wake-NREM sleep cycle, hidden from the EEG-based PSG. For this purpose, we examined the transition probabilities of the HMM states, extracting modules of HMM states that transitioned more often between each other than to other states—as recently identified for the waking resting state in ref. 24.

The whole-brain network states organised into a transition map as presented in Fig. 4, where the 19 × 19 transition probability matrix (Fig. 4a) was submitted to a modularity analysis (see Methods). By considering the most frequent transitions between the HMM states that were consistent across participants (see Fig. 4b), the thresholded transition matrix organised into four partitions or transition modules (see Fig. 4c, and Methods), suggestive of a lower time scale (see ref. 24 and Supplementary Discussion 1). When these most consistent transitions are presented as a transition map, and each whole-brain network state is represented by a circle plot indicating its specificity for each of the four PSG stages, it can be seen that the HMM states exhibit a strong temporal structure (Fig. 4d). In line with the MANOVA results above, this transition map describes an overall progression from states with high specificity for PSG-defined wakefulness (red module) through states with more activity during, albeit not significant specificity for, N1. From here transitions lead towards states specific to N2 sleep and finally to a single whole-brain network state modelling N3 sleep. The N2- and N3-related HMM states thus grouped together in the blue module.

Fig. 4 Investigating modules of transitions between whole-brain network states. a The figure shows the 19 × 19 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, yielding each matrix entry: transition probability from departure state to 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. The transitions describe a passage from HMM states with high activity during wakefulness in the top, further down through HMM states including more N1, and down to HMM states specific to N2 and finally N3. Interestingly, wakefulness appears to be represented by two modules (red and black). Even though no individual HMM state showed clear specificity for N1 sleep, a module (white) is evident between HMM states with specificity for wakefulness and HMM states with specificity for N2 and N3 sleep. e Pie chart showing the total proportion PSG stages within the 18 participants. Transition prob.: transition probability Full size image

Interestingly a collection of HMM states with mixed PSG-specificity formed a transition module of their own. This white module was intercalated between the red module of wakefulness in the top and the blue module of N2/N3 sleep. Even if the included HMM states were not specific for PSG-defined N1 sleep, the white module appears in the location of the transition map, where one would expect to find N1 or rather sleep onset.

The transition map suggested two sub-divisions of HMM states with high specificity for wakefulness. The red module in close proximity to the white module of N1-related states, and the black module sending transitions to the blue module of consolidated NREM sleep. This apparent separation of wakefulness and the asymmetric relationship to the sleep-related HMM states led us to the hypothesis that one of these could represent wakefulness after sleep onset (WASO). Given the poor correspondence between the HMM states and the general uncertainty associated with the staging of PSG-defined N1 sleep (see Discussion), we chose to define WASO as PSG-staged wakefulness, which followed after visits to N2 sleep28. By computing the sensitivity and specificity of the whole-brain network states in the subset of the data corresponding to the 31 participants who woke up after having reached N2 sleep (see Supplementary Table 2), we were able to confirm this hypothesis. As shown in Supplementary Figure 2, HMM states 5, 17 and 18 were all more sensitive and specific to WASO compared to wakefulness prior to N2 sleep. Whereas periods of wakefulness prior to and after sleep are scored equally in PSG, the whole-brain network states separated these into two different transition modules.

Although this transition map suggests multiple pathways from wakefulness (red module) to the white module of NREM sleep, it is interesting to note that HMM state 8 has direct access to HMM state 15, which in turn guards the transition to the blue module of N2/N3 sleep. Similarly, waking-up relates to a transition from HMM state 4 to HMM state 10, which in turn connects with HMM state 18 of the black WASO module. Further it is worth noticing the strong triangular transition structure within the blue module between the N2-specific whole-brain network states (HMM states 3 and 6) and the N3-modelling HMM state 16.

In summary, while agreeing with the overall sequence of PSG stages, the organisation of the transition modules also points to aspects of sleep-related brain activity that the PSG scoring cannot access, including the data-driven suggestions of N1 sleep, WASO-related whole-brain network states, and multiple transition pathways between wakefulness and sleep.

Spatial activation and FC maps of whole-brain network states

We present the spatial maps of the whole-brain network states in the order suggested by the transition modules of Fig. 4d. Figure 5 and Fig. 6 show the mean activation maps, while the corresponding FC information is presented in Supplementary Figures 3–5 and 17–18 (see also Supplementary Note 5).

Fig. 5 Mean activation distributions of wakefulness-related HMM states. a Three HMM states identified as sensitive and specific to wakefulness prior to sleep. Note HMM state 8 with DMN-like configuration increases, and concomitant decreases in DAN/CEN-like areas. b Three HMM states associated to wakefulness after sleep. There are marked frontal increases in mean activation in HMM states 5 and 17. All maps were thresholded above the 50% strongest positive and negative changes, respectively Full size image

Fig. 6 Mean activation distributions of sleep-related HMM states. a HMM states associated with N1 sleep showed opposite signs in mean activation in subcortical areas and primary sensory areas of the cortex. b Three HMM states related to N2 sleep. HMM states 3 and 6 in particular showed peak increases and decreases, respectively, in areas previously identified as fMRI-correlates of sleep spindles. c HMM state 16 is dominating slow wave sleep (N3). Interestingly, there are marked decreases in mean activation in frontal areas and insula, and very localised increases in the supplementary motor area and paracentral lobule. All maps were thresholded above the 50% strongest positive and negative changes, respectively Full size image

In Fig. 5a, which shows the red module of wakefulness, the mean activation maps of HMM states 2 and 8 resemble resting-state network (RSN) configurations9,10. The main increases of HMM state 8 were thus seen in key areas of the default-mode network (DMN)29, including the bilateral posterior cingulate cortex, bilateral angular cortex, bilateral middle temporal cortex, and bilateral medial prefrontal cortex. These DMN-like increases in HMM state 8 were accompanied by decreases in the so-called anti-correlated network (ACN), involving the supramarginal gyrus and the dorsolateral part of the frontal cortex30. In contrast, HMM state 2 was characterised by increases in many of these ACN-areas, including the bilateral supramarginal gyrus, middle cingulate cortex and dorsolateral part of the frontal cortex. These results suggest an inverse relationship between the activity of the DMN and the ACN, which is an established trait of these RSNs30. Since the discovery of these RSN patterns they have been hypothesised to reflect complex cognitive processes. The DMN has been linked to inwardly directed mentation, such as autobiographical memory and mind wandering31,32. The ACN overlaps with areas also referred to as the dorsal attention network33 or the central executive network (CEN)34, and has been proposed to be involved in more externally directed processes, including attention35. In agreement with this, we found these high-order RSNs to be relatively exclusive for wakefulness. However, previous investigations have suggested a rather ubiquitous presence of both the DMN and the ACN, not just in wakefulness but, in all stages of NREM sleep15,17,36 (see Supplementary Discussion 1).

Figure 5b shows the mean activation maps of whole-brain network states with higher sensitivity and specificity for WASO (black module). The mean activation map of HMM state 18 expressed a distribution similar to that of HMM state 8 (see Fig. 5a), but with opposite signs. Hence, HMM state 18 showed decreases in DMN-related areas, and increases in regions overlapping the ACN. HMM states 5 and 17 were both characterised by mean activation increases in the frontal cortices. Interestingly, findings from high-density EEG studies of participants waking up from sleep show that the posterior parts of the cortex are particularly ‘slowʼ at returning to levels of activity seen prior to sleep37. Consistent with this, converging evidence from PET and fMRI have indicated frontal cortical activity to be increased relative to that of posterior areas upon awakening28,38,39.

The whole-brain network states of the white N1-related module are represented in Fig. 6a. A general observation for these spatial maps is the inverse relationship between mean activation in subcortical areas (thalamus and parts of the basal ganglia) and primary sensory cortical areas. Increases in subcortical activity were accompanied by decreases in primary sensory areas of the cortex and vice versa. This was true for HMM states 4 and 15 (and HMM state 1 although its decreases were not confined to subcortical areas, but supplemented by decreases in the anterior and middle cingulate cortex). This is consistent with intracortical studies of the sleep onset process in rats40 and in humans41 showing that thalamic changes in dynamics precede those of cortical areas near the onset of NREM sleep. Previous fMRI studies of NREM sleep have suggested decreased connectivity between the thalamus and cortical regions as perhaps the most consistent trait of FC during N1 sleep16,18,19,23.

N2 sleep was dominated by HMM states 3 and 6, and the mean activation maps of these whole-brain network states are shown in Fig. 6b. The supplementary motor area was involved in both of these states; in HMM state 3 as increases in concert with the bilateral precuneus and primary motor cortices; and in HMM state 6 as decreases together with the bilateral thalamus, middle cingulate, supramarginal cortex, and the rolandic operculum. Interestingly, these configurations overlap considerably with those previously reported in studies mapping fMRI-correlates of sleep spindles42,43, which represent a defining EEG-feature of N2 sleep. However, no HMM state appeared to be driven solely by either sleep spindles or K-complexes. By identifying sleep spindles and K-complexes from the EEG data, we assessed the temporal relationships between these graphoelements and the HMM states. In summary the HMM states that were dominant during N2 sleep showed comparable sensitivity and specificity to both types of graphoelements, and hence the HMM did not appear to have assigned individual states for either spindles or K-complexes (for further information please see the Supplementary Discussion 1, Supplementary Note 4, Supplementary Table 3, and Supplementary Figures 19–21).

HMM state 16 accounted for the majority of time spent in N3 sleep. The corresponding mean activation map is shown in Fig. 6c. Apart from some very localised increases in the paracentral lobule and adjacent supplementary motor area the mean activation was characterised mainly by decreases, particularly in the bilateral middle and superior temporal pole, the orbital part and the operculum of the inferior frontal cortex, bilateral insula as well as medial temporal areas. These frontal decreases are consistent with previous PET findings of decreased metabolism in these areas during N3 sleep, which in turn are believed to reflect the high, localised concentration of slow-wave activity44.