A parsimonious description of global functional brain organization in three spatiotemporal patterns

Resting-state functional magnetic resonance imaging (MRI) has yielded seemingly disparate insights into large-scale organization of the human brain. The brain’s large-scale organization can be divided into two broad categories: zero-lag representations of functional connectivity structure and time-lag representations of traveling wave or propagation structure. In this study, we sought to unify observed phenomena across these two categories in the form of three low-frequency spatiotemporal patterns composed of a mixture of standing and traveling wave dynamics. We showed that a range of empirical phenomena, including functional connectivity gradients, the task-positive/task-negative anti-correlation pattern, the global signal, time-lag propagation patterns, the quasiperiodic pattern and the functional connectome network structure, are manifestations of these three spatiotemporal patterns. These patterns account for much of the global spatial structure that underlies functional connectivity analyses and unifies phenomena in resting-state functional MRI previously thought distinct. The whole-brain organization of functional MRI signals has been studied in myriad ways. An in-depth study of these signals suggests a parsimonious description with a small number of spatiotemporal patterns.

S ince the discovery of spontaneous low-frequency blood oxygenation level dependent (BOLD) fluctuations in the 1990s 1 , increasingly complex analytic techniques have been applied to understand the spatial and temporal structure of these brain signals. A notable feature of these signals is their organization into global patterns that span across functional systems for cognition, perception and action [2][3][4] .
Here we distinguish between two characterizations of this global structure: zero-lag synchrony and time-lag synchrony between brain regions. Zero-lag synchrony is defined as instantaneous statistical dependence between two time courses or the correlation between two BOLD signals with no time-lag. The zero-lag analysis approach has identified several global patterns spanning across functional networks that have generated sustained research interest: the global signal 5 , the task-positive/task-negative pattern 6 and the principal functional connectivity (FC) gradient 2 .
Time-lag synchrony is defined as the statistical dependence between two time courses, where one time course is delayed in time. Two prominent global patterns with coherent time-lag structure have emerged. Short spontaneous global propagation BOLD fluctuations that extend across cortical and subcortical brain regions (~0-2 seconds) are known as lag projections 7 . This propagation pattern varies according to experimental manipulations of task demands and sensory inputs, suggesting that at least some of this structure is uncoupled from hemodynamic delays 7 . A pseudo-periodic spatiotemporal pattern at a longer time scale (~20 seconds), known as the 'quasi-periodic pattern' (QPP), involves an alteration in BOLD amplitudes between the task-positive (TP) and default mode networks (DMNs). The shift in BOLD amplitudes between TP and DMN regions is marked by a large-scale propagation of BOLD activity between the two networks.
There may be an underlying unity to these representations that has heretofore remained overlooked. We hypothesized that the vast majority of widely used zero-lag and time-lag representations of intrinsic functional brain organization capture different aspects of a small number of spatiotemporal patterns that exhibit both zero-lag and time-lag structure. Our specific hypotheses were that (1) global patterns of zero-lag and time-lag synchrony are describing different facets of the same underlying spatiotemporal patterns and (2) a small set of spatiotemporal patterns can explain a large number of previous findings in the literature describing spontaneous BOLD signal fluctuations.
Three lines of evidence support the first hypothesis. First, time-lag representations have a spatial distribution that precisely maps to the spatial weights of the principal FC gradient [8][9][10] . Second, the cortical global signal spatial topography is not entirely constituted by zero-lag spatial structure but has time-lag structure 10,11 . Third, removal of time-lag synchronous patterns, such as the QPP, from spontaneous BOLD fluctuations substantially alters patterns of zero-lag synchrony in FC network representations 12 . These findings suggest that there may be a common pattern of global BOLD activity that unifies these zero-lag and time-lag representations.
To develop intuitions on the proposed relationship between zero-lag and time-lag synchrony patterns, we use the concepts of 'standing' and 'traveling' waves 13,14 . Standing waves refer to stationary oscillations exhibiting no time-lagged statistical dependencies across space (of the kind captured by FC analyses). Traveling waves refer to oscillations in a spatial field with non-zero time-lag statistical dependence across space (of the kind captured by the QPP and lag projection algorithm). We suggest that global BOLD spatiotemporal patterns consist of a mixture of standing and traveling wave A parsimonious description of global functional brain organization in three spatiotemporal patterns Taylor  Resting-state functional magnetic resonance imaging (MRI) has yielded seemingly disparate insights into large-scale organization of the human brain. The brain's large-scale organization can be divided into two broad categories: zero-lag representations of functional connectivity structure and time-lag representations of traveling wave or propagation structure. In this study, we sought to unify observed phenomena across these two categories in the form of three low-frequency spatiotemporal patterns composed of a mixture of standing and traveling wave dynamics. We showed that a range of empirical phenomena, including functional connectivity gradients, the task-positive/task-negative anti-correlation pattern, the global signal, time-lag propagation patterns, the quasiperiodic pattern and the functional connectome network structure, are manifestations of these three spatiotemporal patterns. These patterns account for much of the global spatial structure that underlies functional connectivity analyses and unifies phenomena in resting-state functional MRI previously thought distinct. spatial structure. Zero-lag analyses capture the standing wave structure of these patterns, whereas time-lag analyses capture the traveling wave structure of these patterns. To capture these patterns in a single latent representation, we use a complex-valued extension of a popular dimension reduction technique: complex principal component analysis (CPCA).
In support of the second hypothesis, we begin with the observation that the resting-state functional MRI (fMRI) literature reveals very similar patterns of global BOLD activity across analytic approaches, including FC gradients 2 , co-activation patterns (CAPs) 15 , independent component analysis (ICA) 16 and other latent brain state methods 17 as well as seed-based correlation analyses 6 . We propose that these similar patterns across analysis methods are descriptions of the same underlying spatiotemporal patterns proposed in our first hypothesis. To test the second hypothesis, we compared a systematic survey of zero-lag and time-lag analyses to a set of spatiotemporal patterns derived from CPCA.
Our analyses revealed that three spatiotemporal patterns constitute the dominant large-scale spatial structure in spontaneous low-frequency BOLD fluctuations. With these three patterns, we can unify a range of previous findings in resting-state fMRI, including lag projections 7 , the QPP 18 , the topography of the global signal 19 , the task-positive/task-negative pattern 6 , the principal FC gradient 2 and FC network structure. We demonstrate that all of these previous observations are manifestations of three spatiotemporal patterns captured within a unifying framework capable of modeling standing and traveling oscillatory BOLD phenomena. This novel framework allows for a parsimonious description of global functional brain organization that can inspire new hypotheses about the mechanisms underlying coordination of spontaneous activity across the brain.

results
Standing and traveling wave simulation. We posit that the relative mixture of standing and traveling waves in cortical BOLD signals explains the spatial similarity between outputs of zero-lag and time-lag analyses. To test this claim, we first conducted a simulation study of varying degrees of traveling and standing spatiotemporal wave patterns.
Standing and traveling wave simulations (Supplementary Modeling Note 1) consisted of a back-and-forth sinusoidal oscillation of Gaussian curves on a two-dimensional grid (Fig. 1). This approach allowed us to systematically vary the degree of traveling wave behavior in each oscillation by adjusting the distance between Gaussian curve peak locations, from a distance of zero ('pure standing' motion) to a large distance ('pure traveling' motion). Zero-lag dimension-reduction techniques were applied to the time series of this simulated oscillation.
Zero-lag dimension-reduction analyses applied to pure standing waves effectively captured the oscillation in one latent factor ( Fig. 1a; zero peak distance). Using Catell's scree plot test [20][21][22] , we identified the point where the successive extraction of latent components exhibits a flattening in explained variance as the optimal number of components to extract. For pure standing motion, only the first eigenvalue was non-zero, confirming a single latent factor. Seed-based correlation analysis from a seed at the center of one Gaussian curve returned the same spatial pattern.
At non-zero distance between Gaussian curve peak locations, zero-lag analyses separated the non-synchronous motion into two components. The second eigenvalue of the correlation matrix increased with larger distance, indicating the presence of two latent components. The spatial patterns of the two components were largely consistent across methods. However, methods favoring spatial weight sparsity separated the traveling motion into isolated Gaussian curves, whereas non-sparse methods extracted separate phases or 'snapshots' of the overall oscillation.
Dimension-reduction techniques were extended into the complex-valued domain by augmenting the real-time courses of the grid into a complex signal via the Hilbert transform. To demonstrate the utility of complex-valued dimension-reduction methods in extracting traveling wave oscillations, we applied CPCA to the same empirical data simulations of traveling motion (Supplementary Modeling Note 1). Analogous to our diagnostics of the zero-lag correlation matrix, we constructed a scree plot consisting of the ordered eigenvalues from the complex-valued correlation matrix constructed from the complex grid time courses.
For pure traveling motion, the scree plot from the complex-valued correlation matrix revealed that only the first eigenvalue was non-zero, indicating the existence of a single latent factor ( Fig. 1a; large peak distance). The grid amplitudes of the first complex principal component (PC) reflected the spatial distribution of the two Gaussian curves, indicating that the coherent fluctuations of the two Gaussian curves were captured by a single latent component. Notably, the phase-lag values of the first component precisely reflected the back-and-forth motion of the Gaussian curves.
In simulations with zero distance between peak locations, CPCA and zero-lag techniques converged on a similar solution. Furthermore, at smaller distances, the scree plots of the zero-lag and complex-valued correlation matrix both yielded a onefactor solution.
By construction, this simulation had access to the 'ground truth' mixture of traveling and standing wave components. No such ground truth is available in BOLD signals recorded from the brain. Thus, we invoked a 'traveling index' 14 that measures the presence of traveling waves in CPCA components, varying from 0 (pure standing waves) to 1 (pure traveling waves). Such a metric provides a quantitative estimate of the mixture of traveling and standing waves in the absence of ground truth simulated data.
From the observation that zero-lag methods tend to split traveling Gaussian waves into separate latent components, we used the percentage of explained variance of the first eigenvector as a quality metric. The greater the explained variance, the more effectively zero-lag methods can capture traveling wave motion in a single latent component. We found a linear decrease in explained Gaussian curves at three peak distances (large, moderate and zero distance). For each simulation scenario, the bottom Gaussian curve moves upward (bottom to top) toward the top Gaussian curve. In the top-left panel of each simulation are four arbitrarily sampled timepoints. Note, the travel distance between the two Gaussian curves grows smaller at smaller peak distances (moving from the top to the bottom panel). The amplitude and phase-lag maps of the first complex PC from CPCA are shown in the bottom-left panel of each simulation. For all simulations, the amplitude and phase maps of the first complex component accurately describe the spatial distribution and phase-lag between the two Gaussian curves. In addition, the scree plot of the zero-lag (blue) and complex-valued (green) correlation matrix is displayed. Results of various zero-lag FC analyses are displayed in the right panel. For the dimension-reduction techniques, two components were estimated at non-zero peak distances, motivated by the non-zero second eigenvalue of the correlation matrix. For zero-peak distance, only a single component was estimated. At non-zero peak distances, zero-lag dimension-reduction techniques tend to split the traveling wave oscillation into either (1) separate Gaussian curves (varimax) or (2) separate phases of the oscillation. b, From left to right, plots of the traveling index by peak distance, the variance explained by the first eigenvalue (zero-lag correlation matrix) by peak distance and the variance explained by the first eigenvalue by traveling index. variance moving from values beyond 0.2 to larger values of the traveling index (Fig. 1b). Overall, we found that, for moderate values of the traveling index (<0.5), the explained variance of the first eigenvector is greater than 80%. This suggests that zero-lag methods effectively capture a large majority of the variance of a spatiotemporal pattern with moderate traveling wave behavior. For systems a b    consisting predominantly of standing waves or a moderate degree of traveling waves, we expect solutions of time-lag and zero-lag FC methods to converge.
Three prominent spatiotemporal patterns. To understand the standing and traveling wave components that underlie empirical spontaneous BOLD fluctuations, we applied CPCA to a random  The first component ('pattern one'), representing the leading axis of variance, explains 21.4% of the variance in intrinsic BOLD time series, more than three times the variance explained by the second (6.8%) or third (5.7%) component. The traveling index of the first component was 0.25, indicating a largely standing spatiotemporal pattern, with some traveling wave behavior. The dynamics can be separated into two phases. In the first phase, negative cortex-wide BOLD amplitudes are observed that are strongest in the sensorimotor (SM) cortex, superior parietal (SP) lobule, lateral visual (LV) cortex and superior temporal (ST) gyrus ( Fig. 2a and Movie 1). We refer to these brain regions as the somato-motor-visual (SMLV) complex, noting that it also includes some regions outside sensorimotor cortices (for example, SP and ST). The second phase exhibited a propagation of strong negative BOLD amplitudes in the SMLV toward cortical regions overlapping primarily with the frontoparietal network (FPN) but also with the DMN and V1. This entire spatiotemporal sequence of negative BOLD amplitudes was followed by a spatiotemporal sequence with positive BOLD amplitudes with the same dynamics.
The second component ('pattern two') was the most stationary of the first three components, with a traveling index of 0.14. The overall spatiotemporal pattern can be described as an anti-correlated oscillation between SMLV regions and the DMN. Interestingly, this pattern of anti-correlation or bipolar contrast appeared to delineate the unipolar (all-negative) contrast in the first phase of the first complex PC. Visual inspection revealed a minor traveling wave pattern propagating out from the SM region to the premotor cortex (anterior direction) and SP (posterior direction) in between peak amplitudes of the anti-correlated oscillation.
The third component ('pattern three') was similar to the first component with a traveling index of 0.27. The dynamics can be split into two phases. In the first phase, strong negative amplitudes are observed in the SM, ST and LV, with weak negative amplitudes in the DMN; strong positive amplitudes are observed in the inferior parietal (IP) lobule, inferior temporal (IT) gyrus, premotor cortex, dorsolateral prefrontal cortex (DLPFC) and V1. The second phase was marked by propagation from the IP lobule (anterior direction) and premotor (posterior direction) toward the SM cortex and the IT gyrus (posterior direction) toward the LV cortex.
To further understand what properties of the cortical BOLD time courses give rise to the time-lag structure of the three spatiotemporal patterns, we conducted two null model exercises. First, BOLD time courses were randomly shuffled to demonstrate that the time-lag structure of the spatiotemporal patterns depends on properties of the time courses beyond their zero-lag correlations. Random shuffling of the time courses preserves the zero-lag correlation structure while eliminating time-lag correlation structure. As expected, CPCA of timepoint shuffled data removes the intermediate phase delays (between 0 and π), leaving behind a pure standing wave structure ( Supplementary Fig. 6).
Second, we simulated time courses from a first-order multivariate autoregression (VAR(1)) model fit to the cortical BOLD time series. This model leaves the time-lag (that is, autoregressive) structure of the cortical BOLD time courses intact while assuming Gaussianity, linearity and stationarity of the time courses 23 . CPCA of the simulated time courses generated from the VAR(1) model replicated both the zero-lag and time-lag structure of the three spatiotemporal patterns ( Supplementary Fig. 6). These findings suggest that the three spatiotemporal patterns arise from properties of the BOLD time courses that are mostly Gaussian, linear and stationary.

FC topographies reflect the three spatiotemporal patterns.
The three spatiotemporal patterns are predominantly composed of standing waves. This finding, and the results of the simulation study ( Fig. 2), suggests that zero-lag FC methods should capture most of the variance in these spatiotemporal patterns. We selected a large number of commonly used zero-lag FC methods and applied them to the original time courses from the same resting-state fMRI data. Dimension-reduction methods included principal component analysis (PCA), PCA with Varimax rotation (varimax) 24 , Laplacian eigenmaps (LEs) 25 , spatial independent component analysis (SICA) 26 and temporal independent component analysis (TICA) 16 . Hidden Markov models (HMMs), which are commonly between each vertex of the cortex. The circular average map (middle) represents the average phase delay (in radians) between each vertex of the cortex. As in Fig. 2, the pattern one phase delay map (right) represents the phase delay (in radians) between the vertices within the dynamics of pattern one. The two methods for computing average latency structure exhibiting strong agreement in their propagation patterns (r = 0.83). The strong similarity between the average latency structure and pattern one phase delay map (r = 0.98) indicates that the average latency structure is primarily driven by the spatiotemporal dynamics of pattern one.
used latent state space models for estimating brain states 17 , were also included. Seed-based FC analysis methods included seed-based correlation analysis 6 and CAP analysis 15 with k-means clustering of suprathreshold timepoints into two clusters. Seed-based methods were run for three key seed locations corresponding to major hubs in the SMLV, DMN and FPN-the somatosensory cortex, precuneus and dorsolateral prefrontal cortex 4 . Results were identical for alternatively placed seed regions within these three networks ( Supplementary Fig. 7).
To determine a meaningful number of dimensions in our latent dimension-reduction analyses (PCA, varimax PCA, LE, SICA, TICA and HMM), we again used the scree plot criterion. The scree plot of the zero-lag correlation matrix indicated a clear flattening in explained variance after three PCs (Fig. 2c).
To compare FC topography spatial similarity, we used the spatial correlation between the cortical vertex weight values of each topography. To summarize FC topography similarities, we compared each topography to the first three PC maps computed from PCA. Similarly to the first complex PC of CPCA, the first PC explains 20.4% of the variance in BOLD time series, more than three times the variance explained by the second (6.8%) or third (6.1%) PC. Each FC topography exhibits strong similarities with one or more of the first three PCs (r > 0.6) (Fig. 3a). Statistical significance of the spatial correlation between each FC topography and its most strongly correlated PC was computed using spin permutation tests 27 (n samples = 1,000). All correlation pairs were statistically significant (P = 0.001, lower limit of permutation test). However, there was notable variability in the strength of correlations between each FC topography and one or more of the PCs. Furthermore, strong similarities between FC topographies are seen for more than one PC. Detailed results of seed-based regression/CAP and dimension-reduction analyses are provided in Supplementary Figs. 8 and 9. Overall, our survey revealed a considerable consistency in FC topographies across methods.
Comparison of the spatial weights of the three PCs ( Fig. 3b) with the spatiotemporal patterns (patterns one, two and three) of the three complex PCs (Fig. 2) illustrates that both methods capture similar spatial patterns. The correspondence was one to one-the first PC matching pattern one, the second PC matching pattern two and the third PC matching pattern three. The static spatial weights of each PC appeared as a 'snapshot' within their corresponding spatiotemporal pattern: PC1-pattern one (r = 0.998, t = 11.9 seconds), PC2-pattern two (r = 0.986, t = 12.6 seconds) and PC3-pattern three (r = 0.972, t = 3.7 seconds). Furthermore, the three spatiotemporal pattern time courses closely tracked the first three PC time courses: PC1-pattern one (r = 0.98), PC2-pattern two (r = 0.95) and PC3-pattern three (r = −0.83; temporal lag of ~3 TRs). Given the similarity, we refer to the patterns produced by standard PCA and CPCA as patterns one, two and three, interchangeably.
Comparison with previously observed resting-state phenomena. A further aim of this study was to understand the relationship between these three spatiotemporal patterns and previously observed phenomena in intrinsic BOLD signals. Lag projections 7 and the QPP 18 correspond to time-lagged phenomena at shorter (~2-second) and longer (~20-second) time scales, respectively. Lag projections are computed from the average pair-wise time delays between BOLD time courses and represent the average 'ordering' in time of BOLD amplitude peaks across the brain. An analogous time delay representation can be computed from the average of the . At zero to low thresholding of the FC matrix, the first LE resembles pattern two (PC2). As the threshold is raised, the spatial weights of vertices within the DMN and SMLV become more uniform. At higher thresholds, this results in an eigenmap that resembles pattern one.
pair-wise phase delays of the complex correlation matrix used by CPCA. The spatial correlation between the lag projection and the CPCA-derived lag projection was very similar (r = 0.83), and both exhibit the same direction of propagation (Fig. 4). Interestingly, both average latency structures exhibited a strong spatial similarity with the phase delay map of pattern one (lag projection: r = 0.81; circular mean: r = 0.98) (Fig. 4). This suggests that the average latency structure of spontaneous BOLD fluctuations is largely driven by the first spatiotemporal pattern. Note, due to pre-processing differences, the lag projection that we computed only partially resembles the average lag projection in Mitra et al. 7 , as discussed in Supplementary Fig. 10. The QPP is conventionally derived from BOLD time courses after the application of global signal regression. As shown by Yousefi et al. 28 , two types of QPP can be observed across individuals: one QPP pattern exhibits a global pattern of activity with positive correlation across brain regions ('global' QPP), and the other QPP exhibits a global pattern of anti-correlation between the TP and DMN ('anti-correlated' QPP). Global signal regression was found to consistently produce the anti-correlated QPP across all individuals.
We hypothesized that the global QPP and anti-correlated QPP would correspond to patterns one and two, respectively. To derive the global QPP, we applied the QPP algorithm without global signal regression. To derive the anti-correlated QPP, we applied the QPP algorithm after global signal regression was applied to all BOLD time courses. As hypothesized, the time course of the global QPP derived from non-global signal regressed time courses was strongly correlated with the time course of pattern one (r = 0.71) at a time-lag of 7 TRs (~5 seconds). Furthermore, the time course of the anti-correlated QPP was strongly correlated with pattern two (r = 0.80) at a time-lag of 16 TRs (~12 seconds). We visualized the spatiotemporal template of the global and anti-correlated QPP from the repeated template-averaging procedure and compared it to the timepoint reconstruction (see above) of pattern one and pattern two, respectively. The temporal dynamics of pattern one and the global QPP overlapped significantly (Supplementary Movie 2) as did the temporal dynamics of pattern two and the anti-correlated QPP (Supplementary Movie 3).
At peak amplitudes of pattern one, the distribution of weights is roughly all positive or all negative (Fig. 2a). This suggests that pattern one may track the global mean time course of cortical vertices, or the 'global signal' . Indeed, the time course of pattern one and the global mean time course are statistically indistinguishable (r = 0.97). This also suggests that the temporal dynamics surrounding the timepoints before and after the peak of the global signal correspond to the dynamics of pattern one. We constructed a dynamic visualization of the global signal through a peak-averaging procedure. Specifically, all BOLD time courses within a fixed window (15 TRs on each side) were averaged around randomly sampled peaks (n = 200, >1 standard deviation above the mean) of the global mean time course. This revealed that the temporal dynamics surrounding peaks of the global signal precisely match those of pattern one (Supplementary Movie 2).
The temporal dynamics of pattern two largely represent an anti-correlated pattern between the SMLV and the DMN. This resembles the task-positive/task-negative anti-correlation pattern originally observed by Fox et al. 6 and Fransson 29 . We reproduced these results by creating whole-brain vertex-wide correlation maps using a seed time course from the left and right precuneus, a key DMN node. As expected, an anti-correlated pattern was observed between the SMLV and DMN (Fig. 5a). Additionally, the precuneus seed whole-brain correlation spatial map precisely corresponds to the pattern of BOLD activity at peak amplitudes of pattern two (CPC1: r = 0.92, t = 1.8 seconds). This suggests that the task-positive versus task-negative pattern arises from the anti-correlated dynamics between the SMLV and DMN represented by pattern two.
A similar anti-correlation pattern to the task-positive/ task-negative pattern has been observed in the FC gradient literature 2 and is known as the principal FC gradient (PG). In our zero-lag FC survey (Fig. 3), we computed the PG as the first component derived from the LE algorithm, consistent with Vos de Wael et al. 25 . The PG resembles the anti-correlated spatial structure of pattern two, but this similarity depends on the level of thresholding of the FC matrix input to the LE algorithm (Fig. 5b). The FC matrix represents Pearson's correlation of BOLD time courses between all pairs of cortical vertices (that is, correlation matrix). Generally, thresholding is performed row-wise or column-wise on the FC matrix (for example, 90th percentile of correlation values within that row) 2 . If no FC matrix thresholding is performed, the first eigenmap produced from LE precisely resembles the BOLD activity in pattern two (PC2: r = 0.83; Fig. 5a). As the percentile threshold applied to the FC matrix increases, the spatial weights of vertices within the DMN and SMLV become more uniform (Fig. 4b). At higher percentile thresholds (for example, 90th percentile), a contrast between the SMLV and DMN appears that is similar to the unipolar contrast of pattern one (PC1: r = 0.83) and the anti-correlation contrast of pattern two (PC2: r = 0.82). The reason that the first component of LE returns pattern two, over the more prominent pattern one, is due to timepoint centering ( Supplementary Fig. 11). Network-based representations of FC. The network or graph-based approach to FC analysis models the structure of pair-wise relationships between BOLD time courses. We investigated the degree to which the structure of pair-wise relationships between BOLD time courses arises from the dynamics of the three spatiotemporal patterns. An FC matrix was constructed by computing correlations between all pairs of cortical BOLD time courses (Fig. 6). We compared this matrix to the reconstructed FC matrix from the three spatiotemporal patterns. We estimated the similarity between the two FC matrices by computing the correlation coefficient between the lower triangles of each matrix. Despite a larger mean correlation in the reconstructed FC matrix, we found that the FC matrices were strongly correlated (r = 0.77).
We also sought to determine whether the community structure of cortical BOLD time courses can arise from the shared dynamics of the three spatiotemporal patterns. We estimated network communities from the original and reconstructed FC matrix using the Louvain modularity maximization algorithm. A measure of similarity between community structures, normalized mutual information (NMI), showed strong community structure similarity between the original and reconstructed FC matrix (NMI = 0.73).

Discussion
The goal of this study was to provide a parsimonious taxonomy of prominent global spatiotemporal patterns in spontaneous BOLD activity to enable insight into the functional architecture of the human brain. Using a complex-valued extension of PCA, we identified three prominent spatiotemporal patterns consisting primarily of standing waves with some traveling wave behavior. Consistent with our simulations, the relative predominance of standing over traveling wave structure ensured that zero-lag FC methods effectively captured much of their spatiotemporal structure. This finding explains previous observations that traveling waves of BOLD activity resemble patterns of large-scale FC topographies 8,9 .
A notable finding from this study was the ubiquitous presence of three spatiotemporal patterns across a wide variety of zero-lag FC techniques. Several zero-lag FC methods (ICA, varimax and HMM) involved an initial PCA dimension-reduction step. Thus, the common patterns resulting from zero-lag FC methods may be caused by a common reliance on PCA. However, an ample number of studies demonstrate that the output similarity across these methods is not guaranteed but depends on the area of application 13,30 . Furthermore, the same patterns emerged from zero-lag FC methods that do not rely on PCA (for example, CAP and seed-based correlation) and also from time-lag analyses with no relation to PCA (QPP and lag projection).
Another core finding of our study is the identification of the brain's average latency structure (Fig. 4) with pattern one. This latency structure, originally discovered by averaging the pair-wise time delays between brain regions 7 , was found to map directly to the phase delay pattern of pattern one. In other words, the average latency structure of spontaneous BOLD fluctuations reflects the traveling wave structure of the dominant spatiotemporal pattern (pattern one). Furthermore, a repeated peak-averaging procedure demonstrated that this time-lag structure corresponds to local propagation patterns surrounding peaks of the global mean time course. This is consistent with previous findings of a global propagating wave event surrounding peaks of the global mean time course 10 .
Although time-varying or dynamic FC was not examined in this study, fluctuations between the spatiotemporal patterns identified here may explain the variability in zero-lag FC structure over time. For example, a new time-varying FC analysis, known as edge time series analysis 31 , has found that patterns of co-fluctuations in zero-lag FC resemble the spatial structure of the three spatiotemporal patterns 32 . Although non-stationarity of zero-lag FC is a controversial topic 33 , modeling work has shown that time-varying FC is consistent with stationary time series models 23 . Our null modeling results (Supplementary Fig. 6) suggest that a stationary, linear autoregressive process sufficiently describes the zero-lag and time-lag correlation structure of the three spatiotemporal patterns.
As the primary aim of our study was descriptive, we have avoided any explanatory or causal explanation of the possible neuronal or non-neuronal origins of these three spatiotemporal patterns. However, the identification of pattern one with the global mean time course suggests potential candidate causal mechanisms. One promising candidate is a systemic circulatory origin 34 . A considerable portion of low-frequency BOLD signal variance is correlated with systemic oxygenation fluctuations in the periphery 11 . This low-frequency peripheral oscillation tracks the global mean fMRI time course 35 and exhibits traveling wave behavior induced by differential blood transit time in the cerebral vasculature 36 . Potential sources of this systemic circulatory effect include vasomotion, fluctuations in arterial CO 2 and/or Mayer waves 34 . However, other studies support a neural origin for the global signal 37,38 . We recently demonstrated that individual differences in global signal topography in the human brain are related to life outcomes and psychological function 19 . This work follows studies using pharmacological intervention 38 and electrophysiological recordings in macaques 37 providing evidence for neural origins of the global signal. Momentary fluctuations in arousal and/or vigilance are also related to the global BOLD signal 37 . Some physiological processes strongly co-vary with cortical excitability 39 , making it even more difficult to disentangle the neural versus non-neural sources of the BOLD fluctuations. Taken together, these findings suggest that the global signal (and, by extension, pattern one), although influenced by vascular and other systemic physiological sources, is at least partially driven by neuronal sources. The current study does not resolve these questions about the sources of the BOLD fluctuations and underscores the importance of understanding vascular components of the global signal to effectively denoise fMRI data while preserving neuronal signals 40 .
Candidate origins of patterns two and three are more difficult to identify. Pattern two closely resembles the task-negative/ task-positive pattern originally discovered in resting-state fMRI 6 . A similar anti-correlated pattern is consistently observed in response to task stimuli 41 . Pattern two strongly resembles the QPP obtained after global signal regression, which has been linked to infraslow electrical activity 42,43 . Patterns of propagating activity have also been observed in optical fluorescence imaging 44 and electrophysiological recordings 45 , which demonstrates that time-lagged relationships can arise from neural as well as hemodynamic processes. Future research may be directed toward a more complete understanding of the common or distinct neuronal or physiological mechanisms that give rise to these spatiotemporal patterns.
Although these three spatiotemporal patterns explain a large portion of variance in spontaneous BOLD fluctuations (~32%), a considerable portion of variance is left unexplained. Notably, these three components did account for much of the observed whole-brain functional connectivity. Thus, the remaining, larger fraction of variance may contribute less to the global spatial structure of BOLD signals. The sources of this remaining signal component remain unclear. Although some portion is likely to arise from non-neural effects, there is an intriguing possibility that it also contains activity related to ongoing cognition 46 . We speculate that it may reflect more spatially focal variance in spontaneous BOLD fluctuations, possibly more closely related to neurovascular coupling than the global patterns observed in this study.
It is important to qualify the assumptions on which our comparisons were based. First, the primary metric that determined the number of latent dimensions in our analysis was explained variance. The number three was chosen based on the observation of diminishing explained variance with increasingly higher-dimensional solutions. Although alternative methods for the choice of dimensionality have been proposed, we found that the criterion of explained variance achieved high robustness. Second, most of the analysis approaches that we surveyed can reveal finer-grained spatial insights at higher component or cluster numbers (for example, ICA or varimax PCA). However, we do not expect the same consistency in analytic approaches at finer-grained levels of analysis. Empirical examination of consistency in zero-lag FC analyses at solutions greater than three components confirms this: higher-dimensional solutions yield less consistent FC topographies (Supplementary Fig. 12). Notably, the consistency of analysis approaches at our low-dimensional level of analysis suggests that these three spatiotemporal patterns are robust phenomena. Finally, these three spatiotemporal patterns were modeled at the population level; inter-participant variability in the spatial structure of these patterns was not considered in our analyses. However, an examination of the first three PCs across a sample of participants reveals notable inter-participant variability in their spatial structure (Supplementary Fig. 13).
In summary, we identified three spatiotemporal patterns that parsimoniously recapitulate major findings from a wide range of analytical techniques. Furthermore, these patterns account for much of the large-scale structure that underlies all FC analysis and, therefore, affect how everything from functional networks to graphs are interpreted. As the study of the brain as a complex system advances, these three spatiotemporal patterns have potential applications in better constraining generative models of brain activity, predicting variability in response to external stimuli, informing targeted modulation of brain activity to achieve a particular state and understanding interactions between global brain states and localized activity in neuronal circuits. The concise characterization of systems-level brain activity in three spatiotemporal patterns will facilitate the cross-scale research needed to link fundamental neuroscience studies and human behavior.

online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41593-022-01118-1. Resting-state fMRI data. Our study used resting-state fMRI scans from the HCP S1200 release 47 . Participants were unrelated, healthy young adults (ages 22-37 years). Resting-state fMRI data were collected over two consecutive days for each participant and two sessions, each consisting of two 15-minute runs, amounting to four resting-state scans per participant. Within a session, the two runs were acquired with opposite phase-encoding directions: L/R encoding and R/L encoding. We selected a single 15-minute scan from a random sample of participants (n = 50, 21 males) on the first day of scanning. We balanced the number of L/R and R/L phase-encoding scans across our participants (n = 25 for each encoding direction) to ensure that results were not biased by acquisition from any given phase-encoding direction. We chose a single 15-minute scan per participant to ensure that the phase-encoding/phase-decoding parameter and the imaging session (two resting-state scans per imaging session) did not differ within the same participant. A second independent random sample of participants (n = 50, 22 males) was used as a validation sample. We selected surface-based CIFTI resting-state fMRI scans (MSMall registered) that had been previously pre-processed with the HCP's ICA-based artifact removal process 48 to minimize effects of spatially structured noise in our analysis. All brain imaging data were acquired on a customized Siemens 3T Skyra at Washington University in St. Louis using a multi-band sequence. The structural images were 0.7-mm isotropic. The resting-state fMRI data were at 2-mm isotropic spatial resolution and with TR = 0.72 seconds. Further details of the data collection and pre-processing pipelines of the HCP can be found elsewhere 47,48 . Informed consent was obtained from all participants, and all methods were carried out in accordance with relevant guidelines.
Resting-state fMRI pre-processing. Resting-state fMRI scans were spatially smoothed with a 5-mm full width at half maximum kernel using the surface-based smoothing algorithm in Connectome Workbench version 1.4.2. Resting-state fMRI signals from each vertex were then temporally filtered to the conventional low-frequency range of resting-state fMRI studies using a Butterworth band-pass zero-phase filter (0.01-0.1 Hz). Due to (1) the computational complexity of our analytic pipeline, owing to the large number of analyses studied, and (2) our interest in global, spatially distributed patterns, resting-state fMRI scans were then resampled to the fs4 average space from FreeSurfer 49 . This step downsampled the total number of vertices in the left and right cortex to 4,800 vertices. In group analyses, we Z-scored (to zero mean and unit variance) the BOLD time series from all vertices before temporal concatenation of individual scans. All analyses were applied to group-level data formed by temporal concatenation of participant resting-state scans.

CPCA.
To extract traveling wave patterns, we applied PCA to complex BOLD signals obtained by the Hilbert transform of the original BOLD signals. We refer to this analysis as CPCA. This technique has been referred to as complex Hilbert empirical orthogonal functions in the atmospheric and climate sciences literature 50 or complex orthogonal decomposition in the engineering and physics literature 14 .
CPCA allows the representation of time-lag relationships between BOLD signals through the introduction of complex correlations between the Hilbert transformed BOLD signals. The original time courses and their Hilbert transforms are complex vectors with real and imaginary components, corresponding to the non-zero-lagged time course (t = 0) and the time course phase shifted by t = pi/2 radians (that is, 90°), respectively. The correlation between two complex signals is itself a complex number (composed of a real and imaginary part) and allows one to derive the phase offset (and magnitude) between the original time coursesthat is, the time-lag at which the correlation is maximum. CPCA applied to the complex-valued correlation matrix produces complex spatial weights for each PC that can give information regarding the time-lags between BOLD time courses. In the same manner that a complex signal can be represented in terms of amplitude and phase components (via Euler's transform), the real and imaginary components of the complex PC can be represented in terms of amplitude and phase spatial weights. Of interest in this study is the phase delay spatial map that represents the time-lag between pairs of BOLD time courses-that is, those cortical vertices with a low phase value activate earlier than cortical vertices with a high phase value. Notably, the PCs from the CPCA retain the same interpretive relevance as the original PCA-the first n PCs represent the top n dimensions of variance in the Hilbert transformed BOLD signals. CPCA was implemented with singular value decomposition of the group-wise temporally concatenated complex-valued time series using the fast randomized singular value decomposition (SVD) algorithm developed by Facebook (https://github.com/facebookarchive/fbpca).
Estimating the time scale of complex PCs. For simplicity, the phase spatial maps of each complex PC are displayed in seconds (Fig. 2) as opposed to radians. However, the conversion of phase values (in radians) to time units (seconds) requires an estimation of the time scale of each complex PC. The phase spatial maps of the complex PCs have no characteristic time scale other than that imposed by our band-pass filtering operation (0.01-0.1 Hz-that is, 100 seconds to 10 seconds) in the pre-processing stage. To approximate a unique time scale within this frequency range for each component, we calculated the average duration for a full oscillation of each complex PC using the temporal phase of the complex component time series. This was calculated by fitting a linear curve to the unwrapped temporal phases of the complex PC time series. The slope of the curve was then used as an estimate of the average duration in radians of a TR (0.72 seconds) or timepoints. To estimate the average duration in TRs of a full oscillation, we divided a full oscillation (2 radians) by the duration in radians of a TR. For example, for a TR duration of 0.5 radians, the duration of a full oscillation (2 pi radians) would be approximately 12.6 TRs. Using this procedure, we found that the average duration of the first three complex PCs were ~28 seconds (38.7 TRs), ~27 seconds (37.4 TRs) and ~28 seconds (39.1 TRs), respectively. Using this duration as an estimate of the characteristic time scale of each complex PC allows us to provide an estimate of the time delay in seconds of the spatial phase map. For example, for the first complex PC, a 360° (2 pi radians) phase difference between two cortical BOLD time series would correspond to a ~28-second time-lag between the time series. A smaller phase difference between two cortical BOLD time series, such as 1 radian, would correspond to a ~14-second time-lag between the time series and so forth.
Temporal reconstruction of complex PCs. To examine the temporal progression of each complex PC, we sampled the reconstructed BOLD time courses from each complex PC at multiple, equally spaced phases of its cycle (n = 30; Fig. 2). For each complex PC, the reconstruction procedure was as follows: (1) the complex PC time series was projected back into the original vertex-by-time space to produce time courses of the complex PC at each vertex; (2) the temporal phase of the complex PC time course was segmented into equal-width phase bins (n = 30) spanning a full oscillation of the spatiotemporal pattern (0-2 pi radians); and (3) the vertex values within each bin were averaged to produce a 'snapshot' of BOLD activity at each phase bin (n = 30) of the spatiotemporal pattern. The end result is a spatiotemporal representation of each complex PC in terms of time-varying BOLD activity at equally spaced phases of its cycle.
Traveling index of complex PCs. The real and imaginary parts of a complex PC correspond to the spatial weights of the component at zero and pi/2 (90°) phase shift of the original time courses. In a sense, they encode the temporal evolution of the complex PC from one configuration (the real part) to a subsequent configuration (the imaginary part). By definition, a pure standing wave would exhibit the same spatial configuration from zero to pi/2 (90°) phase shifts of its cycle. A pure traveling wave would exhibit a different spatial configuration from zero to pi/2 (90°) phase shifts of its cycle. This observation suggests a means to quantify the degree of 'traveling' wave behavior of a complex PC using the statistical dependence between its real and imaginary parts. A coherent traveling wave (that is, propagation) of BOLD amplitudes across the cortex would exhibit one spatial configuration at one point in time and a different spatial configuration at another point in time. Thus, a complex PC that encodes this traveling wave behavior would exhibit differing spatial configurations in its real and imaginary spatial weights. Using a metric developed by Feeny 14 , we define the 'traveling' index of a complex PC as the reciprocal of the condition number of the matrix whose two columns are the real and imaginary parts of the complex PC. This metric simply encodes the statistical dependence between the real and imaginary parts of the complex PC. Pure traveling waves would exhibit completely orthogonal real and imaginary parts and a traveling index of one. Pure standing waves would exhibit completely dependent real and imaginary parts and a traveling index of zero.
Analysis of null models. We performed two null model exercises to determine the statistical properties of the cortical BOLD time courses that are necessary to capture the time-lag structure of the complex PCs. Null models allow the selective removal of features of the time courses, such as removing autocorrelation structure and preserving zero-lag correlation structure. First, to illustrate that the spatiotemporal patterns depend on properties of the time courses beyond zero-lag correlation, we randomly shuffled the timepoints of the time courses; this procedure selectively preserves the zero-lag correlation structure of spontaneous BOLD fluctuations and removes the time-lag (autocorrelation) structure. We randomly shuffled the group-concatenated BOLD time courses and estimated three complex PCs from cPCA. Second, we simulated time courses from a first-order multivariate autoregression model-that is, VAR(1)-that was fit to the cortical BOLD time series. A VAR(1) model predicts each time course from preceding timepoints (lag of one timepoint) of itself and all other time courses (that is, all cortical vertices-4,800 time courses). This model leaves the time-lag (that is, autoregressive) structure of the cortical BOLD time courses intact while assuming Gaussianity, linearity and stationarity of the time courses 23 . Due to the large number (4,800 cortical vertices) of highly collinear time courses, we extracted 200 dimension-reduced time courses from the original time courses using PCA. The VAR(1) model was fitted on the PC time courses, and simulated time courses (the same length as the original time courses) were generated from the fitted model. The simulated PC time courses were then projected back into the original cortical vertex space. CPCA was then applied to these simulated time courses.
Zero-lag FC analyses. Description of zero-lag FC analyses. Following the standard terminology of the fMRI literature, we refer to zero-lag synchrony among intrinsic BOLD fluctuations as functional connectivity (FC) 51 . FC between cortical brain regions organize into global, cortex-wide patterns, referred to as 'FC topographies' . All analyses were conducted so as to be as consistent as possible with previous studies. For some of these analyses, results were compared with and without global signal regression. Global signal regression was performed by regression of the global mean time series (averaged across all cortical vertices) on all cortical time series. Residual time series from each regression were then used for subsequent analysis. All analyses were conducted using custom Python scripts and are publicly available at https://github.com/tsb46/BOLD_WAVES. The following zero-lag FC analyses were conducted: • PCA: consists of eigendecomposition of the empirical covariance matrix of the vertices' time series or, alternatively, SVD of the mean-centered group data matrix (time series along rows, vertices as columns). The first T PCs represent the top T dimensions of variance among cortical BOLD time courses. By construction, the first PC is the latent direction of variation with the largest explained variance across all input variables, followed by the second most explanatory component and so forth. The PC spatial weights on each vertex were used to interpret the spatial patterns of each PC. PC scores were obtained from the projection of the temporally concatenated group time series onto the PC space and represent the time course of each PC. • Varimax rotation of PCs: consists of an orthogonal rotation of the PC spatial weights, such that the simple structure of the spatial weights are maximized. Simple structure is defined such that each vertex loads most strongly one component and weakly on all others. We adapted code from the Python package xmca (https://github.com/nicrie/xmca) for implementation of varimax rotation. • LEs (spectral embedding): a non-linear manifold learning algorithm popular in the FC gradient literature 25 . The input to the LE algorithm was the vertex-by-vertex cosine similarity matrix 2 , representing the similarity in the BOLD time series among all cortical vertices. Of note, cosine similarity is equivalent to Pearson correlation in mean-centered and unit-normalized time series (that is, Z-score normalization), as was the case with our data. LE performs an eigendecomposition of the transformed similarity matrix, known as the normalized Laplacian matrix. We also computed LEs with a Gaussian radial basis function (gamma = 1), and the results were virtually identical to the cosine similarity metric. We used the spectral embedding algorithm implemented in the scikit-learn Python package, and details can be found at https://scikit-learn.org/stable/modules/generated/sklearn.manifold.Spec-tralEmbedding.html. • Spatial and temporal ICA: estimates linearly mixed, statistically independent sources from a set of input variables. In the case of spatial ICA, PC axes derived from PCA of the timepoint-by-timepoint covariance matrix are rotated to enforce statistical independence in the spatial domain. In the case of temporal ICA, PC axes derived from PCA of the vertex-by-vertex covariance are rotated to enforce statistical independence in the temporal domain. As with varimax rotation, we input a three-PC solution for both temporal and spatial ICA. We used the FastICA algorithm implemented in the scikit-learn Python package. Details can be found at https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.FastICA.html. • Seed-based correlation analysis: consists of correlations between a seed brain region time course and time courses of all cortical vertices. Seed-based correlation analysis was performed for three seed locations. There are various methods for determining the location of seed regions. In our analysis, we chose seed regions within the three most prominent networks in the three prominent spatiotemporal patterns: SMLV, FPN and DMN. We chose seeds in the somatosensory cortex (SMLV), precuneus (DMN) and supramarginal gyrus (FPN) (Supplementary Fig. 7). The spatial outline of the SMLV, DMN and FPN for guiding the selection of seed regions was determined through a k-means clustering analysis of the temporally concatenated group time series with cortical vertices as observations and BOLD values at each timepoint as input variables (that is, features). We found that a three-cluster k-means clustering solution precisely delineated the spatial outline of the three networks. This spatial outline was used to ensure that the seeds were placed within their appropriate location of each network. In addition, we also tested the robustness of our results for different seed locations in the three networks-medial insula (SMLV), inferior parietal cortex (DMN) and dorsolateral prefrontal cortex (FPN)-and found that the results were identical. • CAP analysis: Three CAP analyses were performed for the same three seed regions used in the seed-based regression analysis. CAP analyses first identify timepoints with the highest activation for a seed time course. Consistent with previous studies 15 , we chose the top 15% of timepoints from the seed time course. The BOLD values for all cortical vertices in the top 15% timepoints are then input to a k-means clustering algorithm to identify recurring CAPs of BOLD activity. We chose a two-cluster solution for all CAP analyses. For each seed, the two cluster centroids from the k-means clustering analysis represent two CAPs associated with the seed time course. • HMM: a probabilistic generative model used to infer the sequence and form of discrete hidden states as well as their transition probabilities from an unobserved sequence of latent states. HMM construes the data-generating process based on multivariate Gaussian distributions conditioned on unknown latent 'brain states' that are assumed to generate the observed cortical BOLD time series. Each brain state represents a recurring pattern of BOLD co-activations/ de-activations, somewhat similar to CAPs. To avoid overfitting and to reduce noise in the high-dimensional input data, we conducted a PCA of the cortical BOLD time series. The first 100 PC projections of the time series served as input to the HMM algorithm. Associated with each brain state is a mean amplitude vector with a value for each PC (n = 100) and a covariance matrix between the 100 PC time courses. The mean amplitude vector represents the pattern of BOLD activity amplitudes associated with that brain state. For interpretation, the mean amplitude vector is projected back into cortical verex space for interpretation. A variety of potential 'observation models' are frequently used in HMM models. As cortical time series are measured on a continuous scale (as opposed to discrete measurements), the probability of a timepoint conditional on a hidden brain state (that is, emission probabilities) is modeled as a mixture of Gaussian distributions. We used the HMM algorithm with Gaussian mixture emission probabilities implemented in the Python package hmmlearn (https://github.com/hmmlearn/hmmlearn). • Modularity analysis: To determine whether the community structure of cortical BOLD time courses can be explained by the three spatiotemporal patterns, we applied the Louvain modularity maximization algorithm to the original FC matrix and an FC matrix reconstructed from the three spatiotemporal patterns. FC matrices were calculated by computing the Pearson correlation coefficient between all pairs of cortical vertex time courses. The reconstructed FC matrix was created from projecting the time courses of the three complex PCs back on to vertex space and computing an FC matrix. The Louvain modularity maximization algorithm was applied with a resolution parameter value of 1 (with asymmetric treatment of negative weights). To compare the degree of similarity between the community structure of the original and reconstructed FC matrix, we computed the NMI between their community assignments from the Louvain algorithm. The NMI varies from 0 (completely independent communities) to 1 (completely identical communities). This analysis was performed using the Python package bctpy (https://github.com/aestrivex/bctpy).

Model selection: choice of number of dimensions in dimension-reduction algorithms.
The dimension-reduction algorithms used in this study, including PCA, PCA with varimax rotation, spatial and temporal ICA and LE, as well as HMMs, require a choice of the number of latent dimensions/hidden states to estimate. For PCA with varimax rotation, spatial and temporal ICA and HMM, this controls the degree of richness and/or fine-grained distinctions of the data description-that is, how many separate unobserved hidden phenomena are assumed and quantitatively modeled to underlie each given data point or observation. We did not assume or try to derive a single 'best' number of latent dimensions to represent intrinsic functional brain organization 52,53 . As we were interested in large-scale cortical patterns of FC, our survey focuses on low-dimensional latent solutions. As an initial estimate of the number of latent dimensions for all choices of dimension-reduction algorithms, we examined the first T dominant axes of variation (that is, PCs) of the correlation matrix formed between all pairs of cortical BOLD time series. Specifically, we examined the flattening (or diminishing return) in explained variance (that is, eigenvalues) associated with neighboring PCs, a procedure known as Catell's scree plot test 20 . According to this test, the number of components to extract is indicated by an 'elbow' in the plot, representing a 'diminishing return' in extracting more components. Clear elbows in the scree plot were observed after a PC solution of one and three (Fig. 3c). We chose the higher-dimensional solution of three components. Note, the elbow in explained variance after three components was independent of the functional resolution (that is, vertex size) of the cortex. We found the same elbow after three components in a scree plot constructed from high-resolution functional scans (~60,000 vertices without downsampling to 4,800 vertices as described above in our pre-processing pipeline). Thus, three latent dimensions were estimated for all dimension-reduction algorithms, and three hidden states were estimated for the HMM.
QPP and lag projections. There are two widely used algorithms for the study of spatiotemporal patterns in BOLD signals: (1) interpolated cross-covariance functions (Mitra et al. 7,54 ) for the detection of lag/latency projections (~0-2 seconds) and (2) a repeated-template-averaging algorithm of similar spatiotemporal segments 18 for detection of the QPP (~20 seconds). Lag projections represent the average time-lag between a brain region's time course and all other brain regions. It provides an estimate of the average temporal 'ordering' of brain region time courses, such that a brain region with a greater average time-lag occurs after a brain region with a smaller average time-lag. For our study, we applied the lag projection algorithm to all cortical vertex time courses. The time-lag between a pair of cortical vertex time courses is defined as the peak of their lagged cross-covariance function. Lag projections are derived as the column average of the pair-wise time-lag matrix between all cortical vertex time courses. An analogous column-wise averaging operation can be applied to the complex correlation matrix used by CPCA. Specifically, we computed the column-wise circular mean of the pair-wise phase delay values of the complex correlation matrix. The circular mean was computed due to the circular nature of the phase values of a complex correlation (that is, −pi and pi are identical phase differences).
To estimate the QPP, the template autoregressive matching algorithm of Majeed et al. 18 was used. The algorithm operates in the following manner. Start with a random window of BOLD timepoints, compute a sliding window correlation of the window across the temporally concatenated group data at each timepoint and then average this segment with similar segments of BOLD timepoints (defined using a correlation threshold). This process is repeated iteratively until a level of convergence is reached. The result is a spatiotemporal averaged template of BOLD dynamics (that could be displayed in a movie, for example), along with the final sliding window correlation time series. The final sliding window time series is the same length as the original participant or group concatenated time series and provides a time index of the appearance of the QPP in BOLD data. Python code for this analysis was modified from the C-PAC toolbox (https://fcp-indi. github.io/). Consistent with previous studies 18, 28 , the following parameters were chosen for the template matching algorithm: the window length was 30 TRs; the maximum correlation threshold for identifying similar segments was r > 0.2; and the algorithm was repeated ten times. The template with the highest average sliding window correlation time series across the ten runs was chosen as the final result.

Statistics and reproducibility.
No calculations were used to predetermine sample size. It is important to note that this analysis was conducted at the within-participant level, with 1,200 sampled timepoints per participant. Thus, the true 'sample size' was substantially larger than 50 data points. Our analyses suggest that n = 50 is sufficient to replicate the three spatiotemporal patterns across independent samples. No experimental conditions were randomized. No data were excluded from analysis. The investigators were not blinded to allocation during experiments and outcome assessment.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Code availability
All code for pre-processing and analysis is provided at https://github.com/tsb46/ BOLD_WAVES. Last updated by author(s): Jun 5, 2022 Reporting Summary Nature Portfolio wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Portfolio policies, see our Editorial Policies and the Editorial Policy Checklist.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Portfolio guidelines for submitting code & software for further information.

Data
Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A description of any restrictions on data availability -For clinical datasets or third party data, please ensure that the statement adheres to our policy Data from the Human Connectome Project (HCP) is publicly available at http://www.humanconnectomeproject.org/data/. Instructions for accessing HCP data can be found at https://www.humanconnectome.org/. All metadata is provided at https://github.com/tsb46/BOLD_WAVES.