Replay bursts in humans coincide with activation of the default mode and parietal alpha networks

1A key mechanism by which the brain forms and stores new knowledge is neural replay, where the patterns of neural activity associated with specific items are spontaneously reinstated in structured sweeps. These sweeps project to widespread regions of the cortex, with physiological signatures known as sharp wave ripples that have been described as the most synchronous events in the mammalian brain. Such spatially dispersed patterns of activity are often initiated during specific states, such as slow-wave sleep, presumably to preclude interference with ongoing wakeful processes. Replay was originally discovered during sleep but is now known to also occur during wakefulness, particularly within periods of immobility and rest. An unanswered question in neuroscience is how the brain orchestrates these structured events in a manner that minimizes interference with ongoing cognition.

2Neuroimaging studies have long highlighted that, during wakeful rest, the brain displays an intrinsic activity structure, cycling through a series of canonical resting-state networks (RSNs). One network that has drawn particular attention is the default mode network (DMN), a disparate set of brain regions that coactivate during “offline” periods, such as between trials in the absence of specific tasks as well as during sleep. The DMN has since been identified as correlating with a number of introspective cognitive states, such as episodic memory and future-oriented thought, suggesting a functional role during wakeful rest in mediating internally generated cognition and inhibiting bottom-up sensory processing. Magnetoencephalography (MEG) resting-state studies have likewise detected DMN activation (along with that of many other networks) with millisecond temporal resolution, demonstrating that these networks activate transiently within specific spectrally defined modes. An ability to measure spontaneous replay noninvasively in humans has now been demonstrated, and this allowed us to investigate a potential link between replay and resting brain network activity.

3Replay during slow-wave sleep is associated with specific electrophysiological patterns; low-frequency oscillations synchronize widespread regions of the cortex, whereas high-frequency sharp wave ripples propagate between the hippocampus and cortical areas. These widespread patterns appear to be integral for effective function of replay in consolidating memories. In contrast, the wakeful brain displays a markedly different profile, with transient periods of synchronous activity interspersed with widespread desynchronization associated with distributed processing of parallel cognitive tasks. Here we wanted to find out whether replay events as detected by are linked to specific changes in whole-brain neural activity, changes that might explain how wakeful replay could reinstate distributed cortical patterns from memory without interference from competing cognitive demands.

4We investigated whether the human replay events discovered by, each representing rapid serial reactivation of learned stimulus representations (Figure 1), coincided with specific macro patterns of resting-brain network activity studied widely in the literature. The focus of our analysis was the same MEG scan data of, collected during resting periods of their experiment.

5We first repeated the analysis of to identify specific moments when replay occurred within these resting-state data. Briefly, we trained multivariate classifiers to recognize each experimental stimulus, applied these classifiers to the resting data, and found that the times when classifiers detected stimulus representations played out in rapid sequences in an order defined by the task (STAR methods; ).

6Next, in the same data, we determined which of a set of canonical RSN states were active at each point in time (STAR methods; Figure 1). We inferred 12 RSN states in a data-driven way using an established hidden Markov model framework (STAR methods; ). These RSN states were labeled according to multidimensional scaling of their distances from each other (STAR methods; Figure S1); thus, RSN states 1 and 12 represented opposing extremes of a single major axis of differentiation between networks. We then conducted an evoked response analysis to find out whether activation of the RSN states was modulated around replay events.

7As shown in Figure 2A, a strong relationship emerged between replay onset and two RSN states in particular: RSN states 1 and 2. This relationship peaked at around t = 0, the exact time of replay onset, but exhibited a decay at either side of this time, with a statistically significant association up to 0.5 s before and after each estimated replay onset time (non-parametric cluster significance test, p < 2e−4 for both RSN states). In addition, a weaker relationship was also evident between RSN states 3 and 4 and the observed replay times (p < 2e−4 for both RSN states). The remaining RSN states display a pattern of decreasing activation that appears approximately in ordinal sequence. Given that the RSN state labels were based on the transition matrix (STAR methods; Figure S1), we then correlated the strength of each RSN state’s replay association with its probability of activation immediately following RSN state 2 (the state most strongly associated with replay) and found a very strong relationship (Pearson’s Rho = 0.77, p = 5.1e−3). Thus, the probability of any RSN state activating at the time of replay is determined by that state’s probability of proceeding fromRSN state 2.

8Each RSN state can be described by its distinct spectral power and phase-locking profile; these profiles are summarized by averaging over frequency bands in Figure 2B for RSN states 1–4. This highlights that RSN state 1 was associated with activity over the parietal cortex; the equivalent network in simultaneous electroencephalogram (EEG)-fMRI studies has been shown to be anticorrelated with the dorsal attention network (DAN); thus, activation of RSN state 1 corresponds to the DAN switching off. RSN state 2 combines high-power signals in frontal and temporal regions with coherent oscillations in the lateral parietal cortex, regions that comprise the DMN. RSN state 3 can be interpreted as activation of the visual cortex, and RSN state 4 can be interpreted as activation of frontal cortical regions. The profiles of the remaining RSN states are shown in Figure S1.

9The alignment of replay with specific RSN states would not be surprising if some of these networks simply reflected the patterns of activity present during the original stimulus encoding in the functional localizer data. Because our methods are based on linear classifiers, we expect the instantaneous patterns between stimulus encoding and replay to reflect activity in roughly the same visual areas (and show that this is indeed the case in Figure S2E); but it is not clear whether whole-brain network patterns of activity would be the same. Figure 2C shows that they are not; if we take the same canonical RSN states, hold them fixed, and fit them to the original functional localizer data on which the classifiers were trained, then we identify a markedly different relationship with the RSN states compared with replay. No significant increase was observed to RSN state 1 (p = 0.88 at t = 0.2 s), and RSN state 2 showed a mild increase that was not significant after Bonferroni correction (p = 0.006 at t = 0.2 s). Directly comparing the distribution of RSN states evoked by replay and by stimulus identified a clear differential distribution with highly significant increases in RSN states 1 and 2 (Figure 2D; one-sided t test, p = 8.9e–7 for RSN state 1 and p = 1.9e−4 for RSN state 2). Although RSN states 3 and 4 do show a significant association with replay, this is only in the sense that they replicate patterns of activity associated with the original stimulus-encoding data. Thus, we conclude that brain-wide patterns of resting RSN state activity that are associated with replay encompass both common and unique patterns of activity with respect to the original stimulus-encoding data. Patterns unique to spontaneous replay are better characterized by activation of the default mode and parietal alpha RSN states.

10To ensure replicability of these results, we conducted the same analysis in a second independent study that examined replay data in MEG using a very similar but slightly amended paradigm. We replicated the exact findings, showing that RSN states 1 and 2 had a strong association with replay (Figure 2A, inset; Figure S2; p < 2e−4 for RSN states 1–4; inflexible direct cluster comparison, p < 0.05 for RSN state 1 and RSN state 2). Neither RSN state 1 nor RSN state 2 displayed a significant increase in relation to the functional localizer data (p = 0.67 for RSN state 1 and p = 0.14 for RSN state 2). A paired t test also confirmed that these RSN states were more strongly associated with replay than with the original training data (p = 5.4e−3 for RSN state 1 and p = 9.1e−3 for RSN state 2).

11The correlation between replay and the RSN states shown in Figure 2B was maintained for over 0.5 s either side of each replay event. This is difficult to immediately reconcile with the highly transient nature of MEG RSN states, which typically activate for less than 100 ms but also display dispersive long-term temporal statistics (Figure 3B). Replicating previous findings, we found that RSN states 1 and 2 displayed longer average lifetimes than any other RSN states and also quiesced for longer periods than any other RSN states (Figure 3B). A qualitative assessment suggested that neither replay nor RSN state activations were distributed evenly over long timescales, with clusters of heightened DMN and parietal alpha network activity in particular coinciding with bursts of replay events (Figure 3A).

12To test whether replay events were concentrated into transient bursts, as suggested by, we first computed the Fano factor over the time course of replay events; Fano factors equal to 1 correspond to a regular, non-bursting process, whereas Fano factors greater than 1 correspond to increasingly irregular bursting. As shown in Figure 3D, the observed replay Fano factor increased as a function of window size and exceeded 1 for all window sizes tested, showing that the occurrence of replay events was increasingly irregular over longer time periods. The bursting nature of replay events was further supported by rejection of a broader null hypothesis that intervals between replay events were independent and identically distributed (p < 1e−3; STAR methods). For further characterization of the time course of replay and spectral properties of the classifier outputs, see Figures S3A–S3C.

13We conducted the same analysis on the time courses of RSN state occurrences (Figure 3E). As with replay events, this showed that visits to different RSN states were also increasingly irregular over longer time periods but displayed a degree of irregularity that was not uniform over the different RSN states (one-way ANOVA, p < 2e−8 for all window sizes tested; this was significant after multiple comparison correction for RSN states and number of windows). This was particularly pronounced for RSN states 1 and 2 (two-sample t test, p < 6.4e−9 for RSN state 1 and p < 0.01 for RSN state 2), the DMN and parietal alpha RSN states that were most strongly correlated with replay. Thus, the RSN states that were most inclined to cluster together into periods of increased intensity over long timescales were similarly the most correlated with replay.

14We have shown evidence that the replay events and RSN state visits display bursty behavior, resulting in clusters of intense activity interspersed with long periods of quiescence. In addition, we have evidence that replay events and RSN state occurrence temporally coincide (cf. Figure 2A). However, this alone does not necessarily mean that the bursting itself temporally coincides. To test whether this was the case, we computed the inter-replay interval time conditioned upon the active RSN state at that time. We found that the interval to the next replay event was significantly determined by the currently active RSN state (one-way ANOVA, p = 2e−9) so that, when RSN states 1 and 2 were active, there were shorter intervals between replay events (p = 2.6e−3 and p = 5e−4, respectively). This suggests that replay events are packaged into bursts that occur selectively during clusters of intense DMN and parietal alpha RSN state activity (for further evidence of RSN state clustering around replay events, see Figure S3G).

15To assess reproducibility, we again replicated all results reported here on a second dataset of 22 subjects (; STAR methods). As shown in Figure S3, replay Fano factors again exceeded 1, increased with window size, and exceeded non-parametric permutation test thresholds (p < 1e−3). Visits to RSN states displayed a similar bursty profile, the degree of which was RSN state dependent (p < 1.5e−8, one-way ANOVA) with the parietal alpha and DMN RSN states displaying the highest amount (two-sample t test, p < 7e−6 for RSN state 1, p < 1.1e−4 for RSN state 2). Inter-replay intervals were again significantly determined by active RSN state (p = 0.01, one-way ANOVA). RSN state 1 was significantly associated with shorter intervals (p = 6e−3), and RSN state 2 was trending in the same direction but not statistically significant (p = 0.07).

16Having established a strong temporal association between replay and specific RSN state activation, we next sought to characterize the nature of the brain-wide activity in the replay-associated RSN states. For each RSN state, we calculated the spatial patterns of oscillatory power over all brain regions and the degree of synchronization (coherence) between all brain regions (according to; STAR methods). Figure 4A plots a wide-band summary measure of power and coherence per region of interest (ROI) and RSN state, showing that the replay-associated RSN states (RSN states 1 and 2 and, to a lesser extent, RSN states 3 and 4) were associated with widespread increases in the power and coherence of oscillatory activity compared with the other RSN states (one way ANOVA for group-wise variation, p < 1e−50 for power and coherence; two-sample t tests, p < 1e−50 for RSN states 1–4 for power and coherence).

17This wide-band summary measure aggregates over frequencies and may therefore be comprised of multiple narrow-band modes in different spatial regions. To better characterize the spatial distribution of activity in distinct frequency bands, we decomposed the spectral activity patterns using non-negative matrix factorization. This identified three prominent frequency modes reflecting activity in canonical frequency bands of delta/theta, alpha, and beta bands, respectively (Figure 4B). RSN state 2, the DMN, showed a prominent elevation in network coherence across all three frequency bands compared with all other RSN states, whereas RSN state 1, the parietal alpha network, showed increased coherence in the alpha band compared with the other frequency bands. The two other replay-associated RSN states were associated with activity concentrated in the alpha band (RSN state 3) and delta/theta band (RSN state 4).

18We also used the RSN state description of the resting state data to calculate the brain-wide patterns of oscillatory power and synchronization occurring specifically around replay events. Figure 5 shows this first as time-frequency plots of power (Figure 5A) and coherence (Figure 5B) averaged over all brain regions, again highlighting a strong increase in power and coherence associated with each replay event. Figure 5C shows the brain-wide patterns of oscillatory power (left) and synchronization (right) that occurred at the time of onset of the replay events. This highlights activity across the frontal and parietal regions of the DMN, with activity in each brain region dissociated into two distinct frequency bands. Frontal nodes of the DMN, including the medial prefrontal cortex and temporal poles, were associated with coherent oscillations in a low delta/theta frequency band; parietal nodes, taken to include the posterior cingulate and lateral parietal cortex, were associated with coherent oscillations in the alpha band.

19Again, in the interests of replicability, all of these results were replicated on a second dataset of 22 subjects (; STAR methods), where all activation maps and spectro-spatial profiles were remarkably consistent (Figure S4).

20found an association between the onset of replay and increases in high frequency (>100 Hz) power, consistent with a model of sharp wave ripple activity that coincides with the detected replay events. Given the strong association between replay events and activity in DMN and parietal alpha networks in our work, we wanted to find out whether power in higher-frequency bands might correlate more generally with activity in these networks. Crucially, our approach for estimating the RSN states applied a low-pass filter to the data with a 45-Hz cutoff frequency; thus, any correlation with power spectra above this cutoff can be interpreted as entirely independent of the original RSN state estimation.

21Figure 6A plots the average high-frequency power spectra over all periods when a given RSN state is active. This demonstrates a significantly elevated power spectral density across all frequencies between 50 and 150 Hz associated exclusively with the DMN RSN state (one-way ANOVA, p < 2.5e−6 for each frequency band between 52 and 148 Hz; p > 0.05 when ANOVA excludes the DMN). This relationship mimics that observed in the power spectral density averaged over 30-ms windows around replay events (Figure 6B) while accentuating the power increase in much higher frequencies. Plotting the spatial distribution of this high-frequency power in the DMN RSN state (Figure 6C) identifies a concentration of power in temporal regions, areas that encompass the hippocampus and neocortical regions known to couple to ripple activity.

22Notwithstanding the broadband nature of this signal and the limitations of MEG in imaging deep sources, our results suggest a relationship between RSN state activity and the high-frequency bursts associated with sharp wave ripples, with such bursts occurring exclusively during activation of the DMN RSN state.

23Figure S5 demonstrates replication of these results on a second dataset of 22 subjects (; STAR methods) where the DMN was found to be significantly associated with elevated high-frequency power (one-way ANOVA; p < 2.4e−6 for all frequencies between 52 Hz and 148 Hz; p > 0.05 when the DMN state was omitted), with a spatial distribution of power concentrated over temporal cortices.

24These results bridge two quite separate fields of enquiry in neuroscience. The study of replay has been predominantly characterized at the level of cellular connections and Local Field Potential (LFP) oscillations in animals, whereas the study of resting brain networks has largely been the preserve of human neuroimaging. Therefore, the link we now establish between these has the potential to extend not only our understanding of replay but also our understanding of the functional role played by human resting brain networks.

25The DMN and parietal alpha activity have parallel histories in the scientific literature, initially interpreted as reflecting idling or default patterns of activity and only subsequently understood to have functional roles supporting attention and cognition. The DMN in particular has since been linked to a role broadly defined as internally oriented cognition, encompassing functions such as episodic memory and future-oriented thought. But in the same way that the brain uses sleep to replay past experience and consolidate memories, our results suggest that healthy waking brain activity may undergo periods of heightened DMN and alpha activity to perform the same function alongside ongoing cognition. Given our more refined mechanistic understanding of the role of replay, our new findings could extend our interpretation of the functional relevance of the DMN. Replay itself is fundamentally understood as a mechanism for memory consolidation but has also been proposed to have more expansive roles in building cognitive maps, preparing neural structures for learning (preplay), and transferring knowledge from the hippocampus to the cortex. This suggests a more expansive role of the DMN in executive control, with the regular transient activations of the DMN associated with building and maintaining stable representations of recently acquired information.

26Replay occurs at highest intensity during slow-wave sleep, when large-scale synchronized oscillations provide an environment conducive to large-scale propagation of sharp wave ripples. Notably, the RSN states that correlate with replay in our study are characterized by large increases in oscillatory coherence, itself hypothesized to support integration of signals from disparate regions of the cortex. The DMN itself appears to be preserved at least into light sleep stages, with further evidence showing that the DMN correlates directly with sharp wave ripples under light anesthesia. Our findings appear to reflect both of these phenomena, low-frequency coherent oscillations coupled to high-frequency power bursts, coinciding directly with estimated replay events during an awake state. Such high-frequency power bursts are taken to reflect increased aggregate rates of action potential discharge, which may or may not coincide with sharp wave ripples. Importantly, the high-frequency band power increases we see alongside DMN activation appear to be much more broadband than what has been established by studies of human sharp wave ripple activity using invasive methods. We do not have a good explanation for this at present, but it is notable that comparable broadband spectra have been detected by studies analyzing ripple band activity in neocortical areas locked to hippocampal ripple events. Furthermore, the replay-evoked spectrum in higher frequencies appears to be more narrowband when baseline correction is applied, indicating that this narrowband/broadband distinction may be sensitive to such methodological choices (Figure S5).

27The role of the parietal alpha network is perhaps more readily understood through the unique requirements of awake-state replay. Our results are based on spontaneous replay of visual stimulus representations, raising the problem of how the brain could reinstate these items without interference from ongoing perception. Alpha oscillations are widely interpreted as an inhibitory signal that acts to gate irrelevant stimuli from active processing. One possibility is that strong alpha activity may combine with DMN activation to inhibit bottom-up sensory perception during inward-oriented attention, supporting replay of items from memory. This could provide a crucial mechanism for how replay plays out without interference from competing sensory inputs during the awake state.

28We have focused on sequential reactivations consistent with task sequence, but similar analyses work when we simply consider all reactivations (Figure S6). It has been established in these datasets that reactivations are more likely to occur in relevant sequences than in control sequences. However, because reactivations are bursty, we cannot make selective claims about sequential reactivations here. Increases in the density of sequential reactivations coincide with increases in the density of all reactivations, and these times align with periods of DMN and parietal alpha activity.

29In our work, we have drawn a clear distinction between spontaneous replay and cued reactivations but acknowledge that it is not clear how important this distinction might be in general. Notably, a number of MEG and EEG studies that have defined replay differently in the context of cued memory recall paradigms or event boundary analyses have reported findings that appear to overlap with our own; in particular, finding strong dependencies on neural dynamics in the theta/alpha frequency range. In the broader context of our findings, this might suggest common neural mechanisms for cued memory recall and spontaneous replay, mechanisms that may be a core function of the DMN.

30The temporal profile of replay activation that we have characterized may also explain replay-related signals across different recording modalities. In particular, fMRI studies have shown reliable behavioral correlations between the reinstatement of blood-oxygen-level-dependent (BOLD) traces associated with experimental stimuli and subsequent task performance. This has supported an interpretation of the BOLD signal as a reflection of cellular replay despite the disparity of timescales between cellular replay events (known to be temporally compressed on the order of milliseconds) and the hemodynamic response (assumed to reflect sustained activations on the order of seconds). However, our results characterize replay as occurring in transient bursts of high intensity interspersed with long periods of quiescence. Such a temporal profile may bridge this disparity of timescales and explain how a hemodynamic signal could arise from clusters of replay bursts in quick succession.

31Furthermore, our results can help to bridge the understanding of RSNs studied in electrophysiology and fMRI. A long-standing challenge in unifying findings across modalities has been to understand how the BOLD response relates to activity in canonical frequency bands. Activity in the gamma band (>30 Hz) has consistently shown a strong correlation with a subsequent BOLD signal; however, a more complex relationship emerges between the BOLD signal and activity in lower-frequency bands, in which electrophysiological RSNs are largely defined. It has now been shown that different RSN states show markedly different hemodynamic profiles; in particular, the DMN and DAN evoke BOLD signals that are opposed in polarity and distinct in their temporal decay profile. If the DMN state visits cluster together in time while being linked to high-frequency power increases, as our results indicate, then this may explain these distinct profiles and provide a key bridge between the understanding of RSNs recorded across distinct modalities.

32Finally, our results suggest measures that could potentially serve as non-invasive proxy measures of replay, potentially opening the door to a broader set of replay experimental paradigms. Until very recently, replay has been predominantly studied in animal models using spatial navigation paradigms because of the necessity of highly invasive electrophysiology to detect replay and the sophisticated understanding of the entorhinal-hippocampal spatial navigation systems. The development of methods to detect spontaneous reactivation and replay in humans noninvasively has enabled experiments that test how these theories generalize to non-spatial domains and other abstract cognitive tasks that are unique to human neuroscience. However, these methods require demanding experimental designs and cognitive paradigms. Our results further broaden the range of tools, providing additional non-invasive measures that could provide an estimate of aggregate replay activity under simple experimental conditions (rest) that can be potentially be studied in large populations or patient groups.

33Overall, our results highlight an important link between two influential domains of research in modern neuroscience, the study of replay and the study of RSNs, and provide a potential connection between noninvasive human imaging studies and invasive cellular physiology.

34Further information and requests for resources should be directed to Cameron.Higgins@ohba.ox.ac.uk.

35This study did not generate any new reagents.

36The code used to run the analyses in this paper is publicly available at the following GitHub repository: https://github.com/OHBA-analysis/Higgins2020_Neuron

37This paper analyses data from three separate datasets, referred to as Dataset A, B and C (see STAR methods). Dataset A and B will be freely available upon request (subject to participant consent) to yunzhe.liu.16@ucl.ac.uk. Dataset C was collected and held by the MEG UK Partnership, and access to this data can be requested at https://meguk.ac.uk/contact

38We make reference throughout this publication to three separate datasets; we will refer to these as Dataset A, the replay data collected by that forms the main focus of analysis in this paper; Dataset B, a secondary replay dataset that was used to prove replication of the analysis conducted on dataset A, which is referred to throughout the paper and is the focus of Figures S2–S5; and Dataset C, an established, larger dataset comprising resting state MEG data that we used to train the canonical resting state networks.

39Dataset A, the primary replay dataset, was acquired from 25 participants (aged 19-34, mean age 24.89; eleven male, two left-handed); four subjects were then excluded due to large motion artifacts or missing trigger information, leaving 21 subjects for the analysis conducted. All participants signed written consent in advance; ethical approval for the experiment was obtained from the Research Ethics Committee at University College London under ethics number 9929/002.

40Dataset B, the replication dataset, was acquired from 26 participants (aged 19-34, mean age 25.48; ten male, two left-handed); four participants were later excluded due to motion artifacts or failure to complete the task, leaving 22 subjects for the analysis conducted. All participants signed written consent in advance; ethical approval for the experiment was obtained from the Research Ethics Committee at University College London under ethics number 9929/002.

41Dataset C, the RSN-state training dataset, comprised resting state MEG data from a larger group of subjects that has previously been used to characterize MEG resting state network dynamics. This study acquired resting state MEG scans and structural MRI scans for 55 participants (mean age 26.5 years, maximum age 48 years, minimum age 18 years; 35 males). All participants gave written informed consent and ethical approval was granted by the University of Nottingham Medical School Research Ethics Committee.

42For full details of the experimental protocol, readers are directed to Liu et. al. 2019; key details only are summarized here. In Dataset A, each participant attended two days, the first for learning the structure of the task and the second for completing the task while undergoing MEG scanning. The task was based around 8 visual stimuli, and participant’s objective was to correctly unshuffle these into the two correct four item-long sequences (eg, A- > B- > C- > D and A’ - > B’ - > C’ - > D’), using a set of unshuffling rules learned on the first day of the experiment. Once in the scanner, participants observed multiple presentations of the novel visual stimuli in a randomized order to act as training data for the multivariate classifiers (the Functional Localizer data. They were then presented with the visual stimuli in the shuffled order from which they could infer the correct sequence. Participants then underwent the first 5 minute resting state scan; one of the two four-item sequences was paired with a reward; participants underwent a second resting state scan, and were tested on their correct recall of the item sequences. All analysis in this paper focuses only on the data from the two resting state scans within this overall experiment. Dataset B used a very similar but slightly amended paradigm; the shuffling rule that participants learned was different, as was the ordering of individual blocks (see ).

43This paper builds on the results of, using their computed replay onset timings. For clarity we outline their methods to detect replay here, but direct readers to the original publication for further details. For each of the eight visual stimuli, sparse logistic regression classifiers were trained using L1 regularization on the functional localizer data to identify the neural patterns associated with each visual stimulus. To select an appropriate time point on which to train the classifiers, the authors used cross validation over the functional localizer data to plot the overall classification accuracy; this identified the time point 200msec following stimulus presentation to be the time point corresponding to the highest classification accuracy averaged over all stimuli, trials and subjects. The trained classifiers associated with this time point were then fit to the resting state data, producing eight time series that represented the probability of reinstatement of the activity patterns associated with each visual stimulus:Where is the recorded resting state data at time; is the sparse classifier coefficients associated with the th visual stimulus; is the probability of reinstatement of the patterns associated with stimulus at time in the resting state scan; and denotes the logistic sigmoid transform which maps from the real number plane to a probability on the interval [0,1].

44The authors analyzed the temporal cross-correlation of these scores and found significant evidence for the reinstatement of stimulus patterns in specific sequence ordering that matched the task structure; that is, scores followed in the patterns A- > B- > C- > D and A’- > B’- > C’- > D’ (see Figure 1B). This patterning was strongest for a time lag of msec, indicating very rapid serial reinstatement of visual stimulus representations. Finally, the authors estimated a single overall time course of replay using the following logical operation:Where denotes the logical AND operation, and denotes the logical OR operation. Consequently, the replay time course,, represents the probability of any task-relevant two item sequence occurring with the specified time lag msec. This probabilistic output was thresholded at the 99th percentile to provide the estimated replay event times used throughout this paper.

45The RSN-states referred to in Figure 1 and throughout this paper use an established Hidden Markov Model with Time Delay Embedding (HMM-TDE) approach outlined in detail in. We outline the approach here and direct readers to the original publication for further details.

46This model makes use of the Hidden Markov Modeling framework, a generative modeling approach that describes the data at each time point as being generated from a latent state variable. That is, the latent state at any point in time is an integer between and, where K is a parameter controlling the total number of states. Therefore, the DMN state being active at time would be denoted by, with the DMN corresponding to RSN-state 2 in this case (Figure 2B).

47To complete the model specification, we must define the observation model – that is, the probabilistic model relating the activation of a particular latent state to the underlying data. As in, we make use of a temporal embedding with a Gaussian observation model:In this notation, the operator performs the temporal embedding; that is, we take a matrix of points centered on the time point, where P is the number of data dimensions and is the length of the temporal embedding - and stack it into a vector with the operator. We then model this vector of datapoints with a Gaussian distribution, with zero mean and a covariance matrix determined by the active state. Crucially, each RSN-state is now parameterized by a unique covariance matrix which reflects the autocovariance on each channel, as well as the cross covariance across channels. As outlined in, this is an efficient parameterization of the power spectrum and the cross power spectrum respectively. Consequently, each distinct RSN-state corresponds to a distinct distribution of power on each channel and coherence between channels.

48In both Dataset A and B, MEG scans were acquired at 600 samples/second on a 275 channel CTF MEG system. For Dataset C, MEG data was acquired at 1200Hz using a 275 channel CTF MEG system operating at third order synthetic gradiometry configuration; MRI data were acquired using a Phillips Achieva 7T system. MRI data, used only for the purpose of MEG coregistration, were acquired using a Philips Achieva 7T scanner.

49Across the three datasets and the two analysis pipelines of Figure 1B we sought to minimize differences in preprocessing however minor deviations were necessary. The replay identification conducted by (Figure 1) filtered sensor space MEG data to a pass band of 0.5 to 50 Hz; data were downsampled to 100 samples / second. All analysis reported by in this paper was conducted in sensor space; for further details see. The Resting State Network analysis pipeline filtered sensor space data to a pass band of 1 to 45 Hz and downsampled data to 250 samples / second. This slightly amended filter passband was to ensure RSN-state dynamics were not driven by low frequency sensor drift effects or mainline power noise effects, which the HMM modeling approach is more sensitive to compared to the replay identification methods introduced by. Similarly, the higher sampling rate was to ensure sufficient resolution of the time embedding to ensure good estimation of spectral content for each RSN-state definition. The only exception to this was for the final analysis of spectral content in higher frequencies (Figure 6), for which we returned to the raw data to utilize the highest frequency information available; this data had been acquired at 600 samples / second and low pass filtered with cutoff frequency 160Hz.

50One of the motivations for using Dataset C to train the resting state networks was to ensure the highest possible confidence and replicability of the anatomical distributions of activity unique to each RSN-state. Given the limited spatial resolution of MEG, confidence can be increased through the use of larger datasets which have been coregistered using MRI acquired structural scans. MRI scans were not acquired for either dataset A or dataset B, motivating our focus on dataset C to determine high confidence spatial topographies of each RSN-state.

51Dataset C was coregistered to MRI structural information using a multiple local sphere forward model. Datasets A and B, which did not have associated MRI structural information, were coregistered using fiducial markers. All data then followed a common pipeline of analysis thereafter, with source space reconstruction performed using an LCM-V beamformer projecting to an 8mm MNI grid. This grid was then parcellated into 38 anatomically defined Regions of Interest (ROIs) derived from an independent component analysis of fMRI resting state data from the Human Connectome Project. Source leakage was then corrected for by orthogonalization as outlined in.

52The full model outlined above is amenable to variational Bayesian methods, which support fast scalable inference for large datasets. The datasets we have used are moderately large – the full model was trained on dataset C, comprising 5 minutes of resting state scans from 55 subjects at 250samples/second, totaling 3.9 million unique data samples. This defined our set of canonical RSN-states, which we then fixed and subsequently fit to datasets A and B, each comprising a total of 10 minutes of resting state scans from 21 and 22 subjects respectively, totaling 3.2 million unique samples and 3.7 million unique samples, respectively. To enable scaleable inference over such datasets, we utilized stochastic gradient variational Bayes methods that iteratively learn a full model through batch training. Furthermore, we use PCA dimensionality reduction steps as outlined in when training the model; to avoid potential misrepresentation of the data, the dimensionality reduced data is only used for model training, all spectral features plotted in figures are based on a multitaper fit back to the original, full dimensional data. Finally, hierarchical models fit with variational methods can be sensitive to local minima. To ensure good convergence, the full model inference was run five times; only the model with the lowest final free energy was kept, the others were discarded.

53As a guiding principle we have maintained the same parameter choice as that previously used in. The parameter controls the total number of RSN-states and therefore the granularity of the solution; we set this to, consistent with previous publications; the same results can be broadly reproduced with other choices of K.

54The HMM-TDE model infers a set of mutually exclusive latent states, which we number from 1 to 12 for ease of reference throughout. The RSN-state numbers themselves are arbitrary, however for ease of interpretation we assign labels using a data-driven approach that minimizes a distance metric between ordinal states. This can be visualized as in Figure S1, and has the result that RSN-states that are nearer in ordinal label are considered nearer in the corresponding state space; for example RSN-state 2 can be considered nearer to RSN-states 1 or 3 than to RSN-state 12.

55We derive this distance metric from the transition matrix directly. Let be the Markov state transition matrix, such that each entry; that is, each entry contains the probability of observing a direct transition from state into state. First, we exclude self-transitions, which are uninformative for the purposes of distances between states; this is done by setting the diagonal entries to zero, then normalizing so each row sums to one, such that our new matrix has entries. Noting that high transition probabilities should correspond to small distances, we then define the distance between two states as the probability of not observing that state transition:Finally, to convert this to a symmetrical matrix of distances (as currently ), we simply average corresponding off-diagonal entries.

56Given the 12 × 12 matrix of distances, we then use multidimensional scaling to identify the single axis that accounts for the most variance in the distance matrix. RSN-states are then labeled 1-12 in the order of their appearance on this axis – see visualization in Figure S1.

57Each RSN-state is defined by a distinct spatial distribution of spectral power and coherence. These characteristics can be interpreted directly from the observation model parameters - however this may emphasize certain traits of the model fitting procedure, in particular the effects of the PCA dimensionality reduction methods applied, rather than the underlying spectra that is observable in each RSN-state. For this reason, we extract the spectral information by fitting a multitaper to the raw data, conditioned on the active RSN-state as introduced by which provides an empirical assessment of the power and coherence for each ROI as a function of frequency. We fit the multitaper using the following parameters: a taper window length of 2 s; frequency resolution of 0.5Hz; frequency range of 1 to 45 Hz; and applying 7 Slepian tapers. As outlined above, due to the increased confidence of spatial RSN-state parameters derived from dataset C, a larger dataset for which MRI structural information was available, the spatial and spectral distributions per RSN-state referred to in this paper were derived from this dataset, supporting their interpretation as a set of canonical RSN-states.

58This information is still very high dimensional so we summarize it through a spectral mode decomposition, obtaining spatial power and coherence maps for a data-driven set of frequency band modes. This decomposition is implemented by non-negative matrix factorization as detailed in and. For summary purposes in Figure 2B and Figure S1, we fit this with two modes, to separate a single wideband mode from higher frequency noise. To explore further the spatial breakdown of different frequency modes, we fit the same non-negative matrix factorization procedure with four modes as in Figure 5 (use of four modes was found to produce broadly stable decompositions and mirrors findings in previous work). These can be derived for each RSN-state; alternatively, to more accurately reflect the full power and coherence distribution at a particular event time, we can use the evoked RSN-state distribution to weight the power and coherence information, producing HMM regularized time frequency plots in Figure 5 and Figure S4.

59Figure 2A analyses the average RSN-state distribution evoked by individual replay events. At the single subject level, we took the probabilistic RSN-state time course; that is, a timeseries with K distinct values between 0 and 1 at each timestep reflecting the probability of each RSN-state’s activation. We then baseline corrected by subtracting the average over all time points for that subject, such that values greater than zero reflect probabilities greater than the session average, and values less than zero reflect activation probabilities below the session average for that subject. We then epoched this time series using the replay times identified above, extracting a window extending from 0.5 s prior to each event to 0.5 s after each event. For each subject, we then computed the evoked RSN-state distribution over this time window by averaging the epoched RSN-state distributions over all replay events. Let us denote by the evoked distribution for subject and RSN-state at timelag from the estimated replay events. Figure 2A plots the mean of these values over subjects, reflecting the expected increase or decrease in RSN-state probability evoked by individual replay events.

60To test the significance of these results two types of statistical test were utilized; to support claims of a significant result at a specific point in time we used a one-sided t test, whereas to support claims of a broad peak of significant points across time we used a non-parametric sign flipping cluster permutation test. The cluster permutation test determined the probability of finding by chance a cluster of consecutive significant points that was the same size or larger than that observed in the data. We defined significance in this context as observing a group level t-statistic exceeding 3; we then generated a null distribution of cluster sizes by randomly permuting 5000 times the sign of each subject’s evoked response parameters before computing group t-statistics, and determining the number of consecutive points exceeding the significance threshold. We could thus compute the chance probability of generating a cluster of the same size observed in the data.

61When testing whether this result replicated on a second study, we applied an inflexible analysis to determine cluster significance. Specifically, we took the clusters in time found in the analysis of Dataset A. We then extracted the corresponding data for that time cluster in Dataset B and defined our cluster threshold as the minimum t-statistic within this temporal cluster. In a final step we then computed (using the same sign-flipping permutation tests) the probability of a cluster of points the same or greater in size than this threshold occurring by chance in Dataset B.

62Figure 2C plots the results of the same analysis on the functional localizer data; epochs were defined by the time of visual stimulus presentation, corresponding to t = 0 in the plot. The thick red line indicates the time at which the replay classifiers were trained.

63Figure 2D then directly compares the evoked RSN-state distribution between the replay and functional localizer data sessions; for ease of presentation, we only focused on the distribution at the exact time of replay onset (t = 0) and the exact time point the replay classifiers were trained on (t = 200msec); for each of the 12 RSN-states we conducted a paired t test to check whether the evoked RSN-state distributions were significantly different across the two conditions.

64RSN-states 1 and 2, which correlate with replay, also have the highest variance in the data (as in Figure 4, these RSN states have the highest broadband power and therefore the highest variance). This leads to increased variance in the classifier outputs when these states are on, which could lead to spurious “reactivations” occurring solely as a result of the increased variance in the classifier outputs.

65To control for this, we note that such spurious effects must be unbiased; that is they should not selectively bias a classifier output toward a particular value. It therefore follows that, when classifier outputs are plotted by percentile, any spurious effects attributed to increased variance should be symmetric across the median. Thus, we repeated the replay evoked RSN analysis while varying the threshold value from 1% to 99% (see Figure S6). This allowed us to characterize the degree to which the relationship was dependent upon the choice of threshold value, and also whether the effect was symmetric over thresholds (suggesting it is driven by the variance of classifier outputs) or strongly weighted toward higher threshold values (suggesting it is driven by the classifier values). These are shown in Figure S6.

66Our main result focuses on the replay evoked distributions; in Figures S6F and S6G we compare this with the reactivation evoked distribution to relate our findings to the extensive literature on reactivations. This analysis was performed by repeating the precise steps outlined above for Replay Evoked RSN Analysis, but replacing the replay time courses with the reactivation time courses – that is, the actual stimulus classifier outputs denoted by … in the Replay Detection methods outline.

67These two measures are however strongly correlated (Pearson’s Rho = 0.45), as we have defined the replay time course itself as a linear superposition of the reactivation time courses. Furthermore the autocorrelation structure of the reactivation time courses and the bursty nature of the replay time course characterized below impose fundamental methodological constraints this specific analysis, such that Figures S6E and S6F should be interpreted accordingly.

68A growing field of literature focuses on the dynamics of classifier outputs using spectral analysis methods. Such studies may, for example, use some cued paradigm, training multiple classifiers on each successive time point after stimulus onset and then analyze in detail the Fourier spectrum of the cross validated accuracy time course over trials. Readers familiar with such analyses may be interested in how they would apply to our paradigm.

69Our methods are different in a number of important ways that limit the potential of these common analysis techniques when applied to our data. First, these analyses rely upon the use of aligned trial data, allowing independent classifiers to be applied to the data at each consecutive timestep from the stimulus presentation time. In contrast, we take a single classifier for each stimulus (taken from the point of maximum decoding accuracy, 200msec following stimulus onset in the functional localizer data) and apply that classifier to the entire timeseries of resting state data. Consequently, whereas the approach as applied by typically identifies complex spectral characteristics emerging in the decoded signal, in our analysis the spectra of the classifier outputs is by definition a static linear function of the cross spectrum (because the decoding weights being used are fixed over time).

70Nonetheless, this can still be informative, and when we apply this analysis to our data in Figure S3E we find the spectrum closely mirrors that of the RSN-states that align to replay; specifically, it reproduces two prominent peaks in the alpha and delta/theta range. These peaks are accentuated when we time-lock this analysis to epochs surrounding replay events; extracting the same modes in this data as used in the main manuscript text, we find a prominent increase in the alpha range (paired t test, p = 0.002) and a weaker but significant increase in the delta/theta range (paired t test, p = 0.045). Due to the methodological differences with the analyses of, as highlighted above, the interpretation of these peaks is more limited. These oscillatory patterns are stronger at replay times, but we do not know whether these reactivations are themselves a function of these oscillations or merely coincident with them.

71Some readers may be interested in what such an analysis would conclude when applied to the replay time courses. However, such an analysis would be inappropriate in this case. Spectral analysis applies well to data that has a roughly Gaussian distribution; it is typically applied to the ‘classifier evidence’, namely the classifier output before a sigmoid nonlinearity is applied. However, as the replay time course is determined by a logic operation (see section on Replay Detection), there is no equivalent to the ‘classifier evidence’ time course upon which these analyses are conducted; the replay time course itself is highly nonlinear as shown in Figures S3B and S3C. The shape of this time course led us to conclude that it was more appropriately modeled as a Poisson point process, whose dynamics are better represented by statistics such as the Fano Factor, than by a spectral analysis such as pursued by, and by others in similar papers.

72To analyze the temporal statistics of replay onset, we modeled the estimated replay event times for each subject and as a Poisson point process. Wherever an interval between events encompassed MEG data samples that had been marked as bad, the whole interval was excluded. We then partitioned the replay time course into non-overlapping windows, for a given window length W, and computed the total sum of replay events observed in each window. The Fano Factor, for a given window of length W, is then given by:Where and are the mean and variance of event counts in each window. The Fano Factor was computed individually for each subject, Figure 3D plots the mean ± ste over subjects for each window length. We then computed parametric t tests on the Fano factor to test the null hypothesis that it was equal to one.

73This alone does not prove conclusively that replay occurrences are transient burst events, but rather that they are not homogeneous Poisson processes. Hence, we furthermore tested a much broader null hypothesis that the intervals between replay events were independent and identically distributed, following the approach of. We generated surrogate data by shuffling the interval times; that is, for each subject, we took the replay event time course, computed the intervals between consecutive events, and randomly permuted the intervals to generate a new event time course. This maintained the exact distribution of replay intervals while removing any dependence between consecutive intervals. We then computed the subject and group Fano Factors exactly as above. With this surrogate data, after 1000 permutations the maximum Fano Factor computed (plotted as the dotted line in Figure 3D) remained well below that of the true time course, confirming that the longer term dispersion arises from correlation over consecutive intervals (p < 1e-3) and allowing us to conclude that replay is characterized by irregular transient bursting behavior.

74We then repeated the analysis for RSN-state activations (Figure 3E). Treating each RSN-state visit as an event, we removed consecutive activations and just used the interval between events to create the equivalent RSN-state visit event time course. As above, we computed the Fano Factor versus window length for each subject.

75To compute the significance of these results, we used a one-way ANOVA to test the null hypothesis that there was no significant difference between the Fano factors of different RSN-states; this was conducted for each window length resulting in 100 multiple ANOVA tests. We only report the highest p value obtained, which remained significant with Bonferroni correction, allowing us to reject the null hypothesis for all window lengths observed. We could then assess a two-sample t test, taking each RSN-state one by one and testing the hypothesis that the observed Fano Factors for that RSN-state were no different from the entire population of remaining RSN-states. For example, we would take the 21 subject level Fano Factors for RSN-state 1, and compare that with the population of size 21 by 11 of Fano Factors for each subject for RSN-states 2-12, then repeat this for RSN-state 2. Again, this test was computed for each value of the window size; only the highest obtained p value was reported, all others were below this and therefore more significant.

76Finally, Figure 3F considers the interval between replay events given the active RSN-state; that is, if a replay event occurs at time t we compute the length of time to the successive replay event, and also the most likely RSN-state to be active at time t. In this analysis, we omit any observations where the interval to the next replay event overlaps with segments marked as bad. For each subject, we thus compute a mean replay interval given the active RSN-state, generating a collection of 21 observations of replay intervals given the active RSN-state. Certain subjects had very few observed replay events in particular RSN-states, resulting in high variance estimates; we therefore omitted from further analysis any means that were derived from fewer than ten observations. We then could then test for significant variation by active RSN-state with a one-way ANOVA followed by two sample t tests, in which we tested whether the replay intervals given RSN-state k was active were significantly different from the replay intervals of all other RSN-states combined.

77Figure S4 is another way to characterize the spectral information unique to each RSN-state. For these plots, we take the RSN-state spectral information as outlined by the multitaper approach outlined above.

78To summarize the wideband response over all frequencies, we take the NNMF decomposition with two modes – that is, with one mode that we interpret as higher frequency noise, and the remaining wideband mode interpreted as a denoised average of the overall frequency spectrum – and plot the RSN-state specific power and coherence for each ROI, given the active RSN-state. The power is given directly by the multitaper approach; coherence is defined per pair of ROIs, so the values plotted on the y axis in these plots is the sum of coherences between a given ROI and all others.

79We can also follow the same procedure for each of the three spectral decompositions corresponding to individual spectral modes; that is, we derive the power and coherence for each ROI and for each given RSN-state, as computed for each of the delta / theta, alpha and beta spectral modes.

80To accurately convey the spectral profile that is linked to each replay event, we can combine the replay evoked RSN-state distribution plotted in Figure 2A with each RSN-state’s spectral information. As outlined above, we have values for the spectral power and coherence of each RSN-state from dataset C. We can initially visualize these by averaging over all recorded ROIs; the left-hand side of Figures 5A and 5B plot the state specific power and coherence values, respectively. To compute the power and coherence patterns that most closely reflect replay times, we take the subject specific replay evoked RSN-state distributions, referred to as above and plotted in the lower panel of Figures 5A and 5B. For consistency of interpretation with Figure 2 these are plotted after baseline correction, however to compute the replay evoked power and spectral density plots we revert to the non-baseline corrected values, such that each is a probability between 0 and 1, with the sum of values over all different states equal to 1 for each time point and subject. We can then compute the expected power and coherence as a function of both time and frequency around each replay event by weighting the state specific PSD and coherence values by the replay evoked state distribution. Note that this analysis was done at the subject level; that is, we computed subject-specific power and coherence estimates as a function of both time and frequency using each subject’s specific replay evoked state distribution. This produced maps of dimension [time x frequency x subjects x channels] in the case of PSD and [time x frequency x subjects x channels x channels] in the case of coherence. Given the high dimensionality of this data we separately summarize the time-frequency information (Figures 5A and 5B) and the spatial frequency distribution (Figure 5C). For visualization in the heatmaps of Figures 5A and 5B we baseline corrected at the subject level, such that zero on the color scale denotes an average PSD and coherence. Figure 5C then expresses how this information is distributed over frequency and space. This plot uses only the distribution defined at t = 0, the exact time of replay estimated by. We decompose the spectral information using the previously introduced non-negative matrix factorization, a data-driven approach that identifies the main spectral modes explaining the data, producing a vector of power values over channels and a matrix of coherence values between channels for each spectral mode and each subject. For the purposes of visualization, PSD values are thresholded by the highest 90% of values. The higher dimensionality of the coherence matrix required a more stringent thresholding approach; for this we applied a Gaussian mixture model threshold test, only plotting connections if they fall into a distinct mixture separated from the global distribution. All these spectral estimation and thresholding steps follow those previously introduced.

81To analyze the higher frequency correlates of specific RSN-state activations, we first took the timeseries of state activations that had been fit to the resting state data as outlined above. Importantly, this analysis had been conducted on data that was downsampled to 250 samples/second and band pass filtered between 1 and 45Hz. For this analysis we therefore returned to the raw data which had been acquired at 600samples/second with a low pass filter applied with cutoff 160Hz. This data did contain power line noise with a baseline frequency at 50Hz and some harmonic noise at integer multiples of 50Hz; we did not attempt to filter this out to avoid introducing any undesirable filter effects. Using the exact state activation times and multitaper approach outlined above, but instead using data with a higher sampling rate that included high frequency content, we estimated the power spectral density associated with each state over all frequencies from 1 to 160Hz. Figure 6A plots the mean ± standard error over subjects of the PSD across this broader frequency band for all states. To compare this to the PSD estimated at replay times, we reused the same multitaper approach but replaced the state activation times with estimated replay times. Specifically, we defined a window of width 30msec around the time of each estimated replay event, and computed the PSD across all frequencies for the data in these windows using the same multitaper approach. To estimate the baseline non-replay high frequency power, we followed the same approach but now using randomly sampled replay times that were shuffled over subjects – this ensured consistency in epoch lengths and temporal profiles while extracting data that was not aligned to replay events in any consistent way. These PSD values are plotted in Figure 6B. Finally, given the exclusive correlation between high frequency power increases and RSN-state 2, we explored the spatial distribution of power in higher frequencies in this specific RSN-state. Figure 6C averages the PSD values over all subjects between 102-148Hz (this band was defined to exclude mainline power noise effects) for each ROI, and plots the spatial topography over all ROIs.

82The replay evoked spectrum in higher frequencies in Figure 6B may appear more broadband than previously reported in. This is in fact due to the baseline correction measures used, with a more narrowband signal emerging when the signal is plotted relative to baseline than when the absolute signal is plotted.

83To compare our findings directly, we recreated Figure 6B using a wavelet decomposition applied to the entire timeseries. This allowed us to compare our replay evoked distribution of Figure 6B – that averages over a time window around replay events – with the full epoch analysis reported by which plots the replay evoked spectrum as a function of both time and frequency around replay events. To produce Figure S5D we took the same 30msec window around replay events which replicates the broadband high frequency spectrum of Figure 6B. We then plotted the difference between this signal and baseline as in Figure S5E. Finally, we plotted this spectrum as a function of both time and frequency in Figure S5F, replicating the figures reported by.