Large-Scale Intrinsic Functional Brain Organization Emerges from Three Canonical Spatiotemporal Patterns

The past decade of functional neuroimaging research has seen the application of increasingly sophisticated advanced methods to characterize intrinsic functional brain organization. Accompanying these techniques are a patchwork of empirical ndings highlighting novel properties of intrinsic functional brain organization. To date, there has been little attempt to understand whether there is an underlying unity across this patchwork of empirical ndings. Our study conducted a systematic survey of popular analytic techniques and their output on a large sample of resting-state fMRI data. We found that the apparent complexity of intrinsic functional brain organization can be seamlessly reduced to three fundamental low-frequency spatiotemporal patterns. Our study demonstrates that a long list of previously observed phenomena, including functional connectivity gradients, the task-positive/task-negative pattern, the global signal, time-lag propagation patterns, the quasiperiodic pattern and the network structure of the functional connectome are simply manifestations of these three spatiotemporal patterns. An in-depth characterization of these three spatiotemporal patterns using a novel time-varying complex pattern analysis revealed that these three patterns may arise from a single hemodynamic mechanism. present conducted a systematic survey of a host of popular analytic techniques and their output on a large-sample of resting-state data. Surprisingly, we found that all analytic approaches we largely consistent in their results, converging on a set of three, low-frequency spatiotemporal patterns. Furthermore, with these three spatiotemporal patterns, we can account for a wide variety of previously observed properties of intrinsic functional brain organization. time-lag


Introduction
Many of us have heard the Indian parable of the blind men and the elephant. This anecdote teaches the perils of missing the 'bigger picture' due to our own limited observations. A group of blind men encounter an elephant for the rst time and try to acquire an overall description of the elephant by discussing which part of the elephant each blind man is touching (trunk, tusk, foot, etc.). Each blind man knew a piece of the truth, but failed to grasp how their individual observations came together in a uni ed whole. This parable may be a suitable metaphor for the current state of our theoretical understanding of intrinsic functional brain organization.
Since the discovery of intrinsic low-frequency blood-oxygenation-level dependent (BOLD) uctuations in the late 1990s, increasingly complex analytic techniques have been applied to understand the spatial and temporal structure of these signals. These applications have generated a patchwork of novel empirical ndings into intrinsic functional brain organization. However, there has been little attempt to assess any underlying unity across this patchwork of ndings. Further, the eld is lacking a unifying conceptual framework for translating empirical ndings across analytic techniques. Our present study conducted a systematic survey of a host of popular analytic techniques and their output on a large-sample of restingstate fMRI data. Surprisingly, we found that all analytic approaches we surveyed are largely consistent in their results, converging on a set of three, low-frequency spatiotemporal patterns. Furthermore, with these three spatiotemporal patterns, we can account for a wide variety of previously observed properties of intrinsic functional brain organization.
We operationalize 'intrinsic functional brain organization' as the spatial correlation structure between cortical BOLD signals in the low-frequency range (~0.01 -0.1 Hz). Spatial dependence is by far the most studied property of low-frequency intrinsic BOLD signals. Importantly, the pattern of dependencies between neural activity across the brain are thought to re ect functional systems for cognition, perception and action (Biswal et al., 1995;Damoiseaux et al., 2006;Fox et al., 2005;Smith et al., 2009).
For the purpose of the present study, we distinguish between three different descriptions of intrinsic functional brain organization: 1) low-dimensional (latent) representations of zero-lag synchrony between brain regions -'functional connectivity topographies', 2) low-dimensional representations of time-lagged synchrony between brain regions -'spatiotemporal patterns', and 3) the network structure of zero-lag synchrony between brain regions over the entire brain -the 'functional connectome'.
Zero-lag synchrony is de ned as in-phase or anti-phase statistical dependence between intrinsic BOLD signals -e.g., the correlation between two intrinsic BOLD signals with no time-lag. Following the standard terminology of the functional magnetic resonance imaging (fMRI) literature, we refer to zero-lag synchrony among intrinsic BOLD uctuations as 'functional connectivity' (FC) (Reid et al., 2019). The pattern of zero-lag synchrony between brain regions across the brain is often summarized using lowdimensional latent representations, what we refer to as 'FC topographies'. A myriad of empirical ndings have been generated from studies of FC topographies, including the task-positive/task-negative pattern (Fox et al., 2005), the primary or principal FC gradient (Margulies et al., 2016) and global signal topography (Li et al., 2019).
Zero-lag synchrony can also be studied from a network-centric or graph-based approach, where the unit of analysis is pairwise relationships between brain regions (Bassett et al., 2011;van den Heuvel and Hulshoff Pol, 2010). Rather than reducing pairwise relationships to low-dimensional representations, the network-centric approach analyzes the structures of these relationships in terms of a graph with brain regions as nodes and edges to represent their pair-wise relationships. We refer to the pattern of pair-wise relationships between brain regions across the entire brain as the 'functional connectome'. Study of the functional connectome has revealed important ndings regarding patterns of inter-communication between brain regions (Power et al., 2011), 'hub' regions of the human brain (van den Heuvel and Sporns, 2013), and the community structure of the connectome (Meunier et al., 2010).
More recently, relationships between BOLD time courses at temporal lags (time-lag synchrony) have been investigated. Recent research has observed replicable patterns of time-lagged relationships between brain regions, including wave-like propagation patterns (Mitra et al., 2014(Mitra et al., , 2015 and large-scale pseudoperiodic patterns of activity (Gu et al., 2020;Majeed et al., 2011). We refer to low-dimensional latent representations of time-lag synchrony between brain regions as 'spatiotemporal patterns'. Studies of spatiotemporal patterns have revealed 1) prominent wave-like, propagation patterns in cortical neural activity at short-time scales (~ 0 -2s) (Mitra et al., 2014(Mitra et al., , 2015, and 2) recurring spatiotemporal patterns of BOLD activity changes at larger time scales (~ 20s) (Majeed et al., 2011). At short time scales, Mitra and colleagues (2014) have observed replicable, time-lagged propagation of intrinsic BOLD activity that are uncoupled from hemodynamic delays. At longer time scales, Majeed et al. (2011) has described a pseudo-periodic spatiotemporal pattern involving an alteration in BOLD amplitudes between the frontoparietal (FPN) and default mode networks (DMN), known as the 'quasi-periodic pattern' (QPP).
To date, investigations of these three properties of intrinsic functional brain organization have been largely independent efforts. Further, there has been little attempt to reconcile empirical ndings across different analytic techniques. For example, what are the interrelationships between different approaches for constructing FC topographies, such as principal component analysis, non-linear manifold learning (Margulies et al., 2016), temporal independent component analysis (ICA) (Smith et al., 2012), spatial ICA (Calhoun et al., 2001), hidden markov models (Vidaurre et al., 2017), seed-based correlation analysis (Fox et al., 2005), and discrete co-activation patterns (Liu and Duyn, 2013b)?
Our study conducted a systematic survey of these three properties of intrinsic functional brain organization using a host of popular analytic techniques applied to a cohort of participants with restingstate fMRI recordings from the Human Connectome Project (HCP, n = 50). Our primary goal was to reduce the complex diversity of empirical ndings on intrinsic functional brain organization to a simpler set of fundamental or canonical patterns. What we found was surprising: all three properties of intrinsic functional brain organization -large-scale FC patterns ('FC topographies'), time-lagged propagation and pseudo-periodic patterns ('spatiotemporal patterns'), and the network structure of pairwise FC relationships ('functional connectome') -can be accounted for by three fundamental or canonical spatiotemporal patterns. More speci cally, we found that all zero-lag FC analyses we surveyed consistently produced one of three spatially overlapping FC topographies. These three FC topographies account for several previously observed properties of functional brain organization, such as the global signal topography, the task-positive/task-negative pattern and primary FC gradient. Using a novel application of a time-varying PCA algorithm, we nd that these three FC topographies correspond to 'static' or time-invariant representations of three fundamental spatiotemporal patterns. We further show that these three spatiotemporal patterns account for two other spatiotemporal phenomena of functional brain organization: cortical lag projections (Mitra et al., 2014) and the QPP (Majeed et al., 2011). Finally, we show that the network structure of the functional connectome can be explained by the shared dynamics of these three spatiotemporal patterns.

Three Dominant Functional Connectivity Topographies in Intrinsic BOLD Fluctuations
We conducted a quantitative survey of widely-used zero-lag functional connectivity (FC) analyses with the primary aim of comparing the spatial overlap between the FC topographies produced by each analysis. These analyses were applied to a random sample (n=50) of human subject resting-state scans (~15min each; n = 1200 time points) from the Human Connectome Project (HCP). For input to all quantitative analyses, subject resting-state scans were temporally concatenated and reshaped into a 2D matrix of time points (rows) by cortical vertices (columns). All analyses in our investigation successfully replicated in an independent sample (n=50) of HCP subjects.
Our survey included several latent dimension-reduction methods, as well as seed-based correlation and co-activation methods. Latent dimension-reduction methods included principal component analysis (PCA), PCA with simple structure rotation (varimax) (Andersen et al., 1999;Thomas et al., 2002), Laplacian Eigenmaps (LE) (Vos de Wael et al., 2020), the commonly used spatial independent component analysis (SICA) (Calhoun et al., 2001), and the more recent temporal independent component analysis (TICA) (Smith et al., 2012). Hidden Markov models (HMM) are a commonly used latent state space model for estimating brain states (Vidaurre et al., 2017), and were also included in our study. Seedbased analysis methods included the traditional seed-based correlation/regression analysis (Fox et al., 2005) and co-activation pattern (CAP) analysis (Liu and Duyn, 2013b) with k-means clustering of suprathreshold time points into two clusters. Seed-based methods were run for three key seed locations corresponding to major hubs in the somatomotor network, default mode, and frontoparietal network -the somatosensory cortex, precuneus, and dorsolateral prefrontal cortex (Yeo et al., 2011). Results were found to be identical for alternatively placed seed regions within these three networks (Supplementary

Results B).
To determine a useful number of dimensions in our latent dimension-reduction analyses (PCA, varimax PCA, LE, SICA, TICA and HMM) (Bzdok et al., 2016;Eickhoff et al., 2015), we examined the explained variance of the principal component solution at a range of dimension numbers (see 'Methods and Materials'). As we were interested in large-scale cortical FC topographies, our survey focuses on the lower end for the number of estimated latent dimensions. To estimate the number of dimensions, we examined the drop-off in explained variance (i.e. eigenvalues) associated with each successive principal component, a procedure known as Catell's scree plot test (Carlson et al., 2011;Cattell, 1966;Ecker et al., 2007;Stetter et al., 2000). This widely used component number selection criterion indicated a clear dropoff in explained variance after three principal components ( Figure 1C). Based on these assessments, we committed to a granularity of three latent neural activity dimensions for the dimension-reduction algorithms for comparability of subsequent analysis steps.
Each zero-lag analysis produced one or more FC topographies with weights for each cortical vertex, representing the degree to which that topography is expressed at that vertex. To compare the spatial similarity between two FC topographies, we used the spatial correlation (Pearson's correlation) between the cortical vertex weight values of each topography. To summarize the similarities among the FC topographies, we compared each topography to the rst three principal component maps computed from PCA. The rst three principal components represent the top three dimensions of variability across cortical BOLD time series. The rst principal component represents the most dominant/leading axis of variance across cortical BOLD time series. The rst principal component explains over 20% (R 2 = 20.4%) of the variance in BOLD time series, greater than three times the variance explained by the second (R 2 = 6.8%) or third principal component (R 2 = 6.1%). As can be observed from Figure 1A, each FC topography exhibits strong similarities with one or more of the rst three principal components (r > 0.6). In other words, all FC topographies in our survey resembled one of the top three dimensions of variability in cortical BOLD time series.
The three principal components can be differentiated most clearly with reference to three cortical brain networks: the default mode network (DMN), the frontoparietal or 'executive control' network (FPN) and the sensorimotor and medial/lateral visual cortices (SMLV) ( Figure 1C). Note that the spatial extent of these three networks changes between principal components, and the reference to these networks is for descriptive simplicity. The rst principal component is distinguished by

Three Dominant Spatiotemporal Patterns in Intrinsic BOLD Fluctuations
All FC topographies in our survey were found to resemble the rst three principal components of cortical BOLD time series. FC topographies are produced from measures of zero-lag synchrony between BOLD time series (e.g. Pearson's correlation coe cient). Thus, FC topographies are unable to represent correlations between cortical BOLD time series at any time-lag. Time-lag relationships may re ect cortical propagation patterns or rapidly-changing sequences of cortical activity patterns. We refer to these timelag structures as 'spatiotemporal patterns.' Examples of such spatiotemporal patterns include short timescale lag projections (Mitra et al., 2014(Mitra et al., , 2015 and the quasi-periodic pattern (QPP) (Majeed et al., 2011). In the following section, we demonstrate that the three dominant FC topographies discovered from our survey of zero-lag FC analyses correspond to 'static' or zero-lag descriptions of three spatiotemporal patterns. We utilized a simple modi cation of PCA for detection of time-lag relationships between cortical BOLD time series. Speci cally, we apply PCA to complex BOLD signals obtained by the Hilbert transform of the original BOLD signals (see 'Methods and Materials'). We refer to this analysis as complex PCA (cPCA).
We applied cPCA to the same resting state fMRI dataset. The scree plot criterion again motivated the choice of three complex principal components (Supplementary Figure C). The relative explained variance between the rst three complex principal components from cPCA was similar to that observed between the rst three principal components from PCA: component 1 (R 2 = 21.4%), component 2 (R 2 = 6.8%) and component 3 (R 2 = 5.7%). Associated with each complex principal component is a time-lag delay map, re ecting the time-delay (in seconds) between cortical vertices (see 'Methods and Materials' for construction of the time-lag delay maps). To examine the temporal progression of each complex principal component, we sampled the reconstructed BOLD time courses from each complex principal component at multiple, equally-spaced phases of its cycle (n=30; see 'Methods and Materials'). Movies of the reconstructed BOLD time courses are displayed in Movie 1.
The time-lag delay maps and reconstructed time courses of each complex principal component are displayed in Figure 2. The rst complex principal component describes a spatiotemporal pattern that begins with negative BOLD amplitudes in the SMLV complex ( Figure 2A). This is followed by a propagation of negative BOLD amplitudes towards cortical regions overlapping primarily with the FPN, but also with the DMN and primary visual cortex. This is followed by a mirrored propagation of positive BOLD amplitudes with the same dynamics. Given this pattern of propagation, we refer to the rst complex principal component as the 'SMLV-to-FPN' spatiotemporal pattern. Because the explained variance of the rst complex principal component is three times greater than the subsequent complex principal components (Supplementary Figure C), we also refer to the SMLV-to-FPN as the 'dominant spatiotemporal pattern' in intrinsic BOLD uctuations. The second principal component describes a spatiotemporal pattern that begins with positive BOLD amplitudes in the DMN and primary visual cortex, and negative BOLD amplitudes in the FPN. This is followed by the onset of positive BOLD amplitudes in the SMLV that quickly propagates towards the FPN, with a simultaneous propagation of negative BOLD amplitudes from the FPN to the DMN. We refer to the second complex principal component as the 'FPNto-DMN' spatiotemporal pattern. The third principal component describes a spatiotemporal pattern that begins with positive BOLD amplitudes in regions of the FPN and negative BOLD amplitudes in regions of the DMN and SMLV. This is followed by a fast propagation of positive BOLD amplitudes from the FPN to the SMLV, with a simultaneous propagation of negative amplitudes from the DMN to the FPN. Note, activation of the SMLV complex occurs slightly before the regions of the DMN. We refer to the third complex principal component as the 'FPN-to-SMLV' spatiotemporal pattern.
Examination of the reconstructed time courses reveals that the pattern of spatial weights in the three dominant FC topographies (i.e. principal component maps. Figure 1B) resembles the pattern of BOLD activity at individual time points of the three spatiotemporal patterns. The pattern of weights of the rst principal component (PC1) occurs within the rst spatiotemporal pattern (r = 0.998, t = 11.9s), PC2 occurs within the second spatiotemporal pattern (r = 0.986, t = 12.6s), and PC3 occurs within the third spatiotemporal pattern (r = 0.972, t = 3.7s). Further, the time courses of the three spatiotemporal patterns components closely tracks the time courses of the rst three principal component time courses, respectively: PC1 (r = 0.98), PC2 (r = 0.95), and PC3 (r = -0.83, at a temporal lag of ~3 TRs). This nding suggests that the three dominant FC topographies are 'static' or 'stationary' representations of three temporally-extended, dynamic patterns of BOLD activity.

Steady States and Propagation Events in Spatiotemporal Patterns
To visualize the temporal dynamics of the three spatiotemporal patterns, we projected the reconstructed time points (see above) into the 3-dimensional space formed by the rst three principal components, corresponding to the three dominant FC topographies ( Figure 1B). This projection allowed for a simple visualization of the temporal progression of BOLD activity within each spatiotemporal pattern. The structure of the resulting projection is interpreted as follows: reconstructed time points with a spatial pattern of BOLD activity resembling the spatial weights of one of the three principal components have higher scores on the axis of that principal component. The bene t of this representation of spatiotemporal patterns is two-fold: 1) examination of the movement of consecutive time points provides information regarding the 'speed' of change in BOLD activity between time points (e.g., steady states vs. rapid propagation or transition events), and 2) time points between two spatiotemporal patterns that are close together in this space indicate common patterns of BOLD activity between those spatiotemporal patterns (see below). The reconstructed time points of each spatiotemporal pattern is also displayed in Movie 1.
The temporal cycle of each spatiotemporal pattern forms an oval in the three-dimensional principal component space ( Figure 3B), corresponding to a full cycle of the spatiotemporal pattern. For all three spatiotemporal patterns, most consecutive time points cluster closely together, indicating a 'steady state' of BOLD activity with relative stability of BOLD activity over that period. The steady states of the SMLVto-FPN, FPN-to-DMN and FPN-to-SMLV vary strongest along the rst, second and third principal component axes, respectively. This is apparent from the location of steady state time points (i.e. consecutive time points clustered closely together in space) in the 2-dimensional plots formed by two principal component axes in the three-dimensional space. For example, the steady state time points of the SMLV-to-FPN (time points displayed in blue) exhibit the highest positive and negative scores on the rst principal component axis with smaller scores on the rst and third principal component axes. This simply re ects the fact that the spatial pattern of BOLD activity during the steady states of the SMLV-to-FPN is strongly correlated with the pattern of spatial weights of the rst principal component. These steady state periods are interrupted by large movement between consecutive time points that correspond to rapid propagation of BOLD activity towards another steady state. All three spatiotemporal patterns spend most of their cycle in a period of steady synchronous activity that is interrupted by rapid propagation events.
In the case of the SMLV-to-FPN, the speed of propagation is relatively slower compared with the more abrupt propagation of the FPN-to-DMN and FPN-to-SMLV. This becomes apparent by tracing consecutive propagation time points (i.e. time points with long distances from their previous time point) between the steady states of each spatiotemporal pattern. The propagation events of the FPN-to-DMN and FPN-to-SMLV travel the same distance in approximately three time points as the SMLV-to-FPN in approximately ve time points. For each spatiotemporal pattern, a full-cycle contains two mirrored steady-states and two mirrored propagation events. Mirrored steady-states and propagation events are the same spatial pattern of BOLD activity with a sign-ip -i.e. ipped positive and negative values. Another notable observation in this representation is that propagation time points of the SMLV-to-FPN and FPN-to-DMN vary most strongly along the third principal component axis. This indicates that the pattern of BOLD activity during propagation events of the SMLV-to-FPN and FPN-to-DMN resemble the pattern of spatial weights of the third principal component. The opposite is true of the FPN-to-SMLV. In this spatiotemporal pattern, steady-state time points vary strongest along the third principal component axis, and the propagation time points vary strongest along the second principal component axis.

Recurring Spatial Topographies in Spatiotemporal Patterns
There is overlap in cortex-wide BOLD activity at certain time points between the three spatiotemporal patterns. For example, the pattern of BOLD activity at 8.6 seconds into the SMLV-to-FPN corresponds closely (r = 0.909) to the pattern of activity at 12.2s into the FPN-to-SMLV. This suggests that the same spatial topography of BOLD activity may appear across more than one spatiotemporal pattern. To examine repeating spatial topographies across the three spatiotemporal patterns, we applied a clustering algorithm to the reconstructed time points from all three spatiotemporal patterns. We concatenated the reconstructed time points from each spatiotemporal pattern (n=90, 30 time points per spatiotemporal pattern) and clustered the time points from the concatenated matrix using a k-means clustering algorithm. To avoid scaling differences in the distance calculations between time points, the BOLD values within each time point were z-score normalized. We chose a six-cluster solution from the k-means algorithm, as this was found to capture the three mirrored pairs of steady states from each spatiotemporal pattern ( Figure 3E). Examination of the cluster assignments of each time point across spatiotemporal patterns ( Figure 3C and D) yields several important insights. First, the six clusters correspond to the three pairs of mirrored or sign-ipped steady-states of the three spatiotemporal patterns. The rst two clusters correspond to the steady states of the SMLV-to-FPN. Clusters three and four correspond to the steady states of the FPN-to-DMN, and clusters ve and six correspond to the steady states of the FPN-to-SMLV. Second, the pattern of BOLD activity in the steady-states of one spatiotemporal pattern occurs within propagation events of the other two spatiotemporal patterns. For example, the pattern of BOLD activity in the steady-states of the FPN-to-SMLV (cluster three and four) occur within propagation events of the SMLV-to-FPN and the FPN-to-DMN. Further, the pattern of BOLD activity in the steady-states of the FPN-to-DMN (cluster ve and six) occur within propagation events of the FPN-to-SMLV. In other words, the same pattern of BOLD activity occurs as a steady state or a propagation event depending on the spatiotemporal pattern. Third, the FPN-to-DMN and FPN-to-SMLV spatiotemporal patterns are mirror images of each other. The pattern of BOLD activity in the FPN-to-DMN steady-states occurs as propagation events in the SMLV-to-DMN, and vice versa. In fact, the FPN-to-DMN can be converted to the FPN-to-SMLV by swapping the steady-states (clusters ve and six) and the propagation events (clusters three and four), and vice versa.

Relationships With Previously Observed Phenomena in Intrinsic BOLD Fluctuations
A further aim of this study was to understand the relationship between these three spatiotemporal patterns and previously observed phenomena in intrinsic BOLD signals. First, we consider spatiotemporal patterns discovered by previous approaches. Lag projections (Mitra et al., 2014) and the QPP (Majeed et al., 2011) correspond to time-lagged phenomena at shorter (~2s) and longer (~20s) time scales, respectively. Lag projections are computed as the column average of the pairwise time-lag matrix. The time-lag between a pair of BOLD time courses is the difference in time at which the correlation between the BOLD time courses is maximal. The column average of the pairwise time-lag matrix, or lag projection, provides the average 'ordering' in time of cortical BOLD time courses. We hypothesized that the average time-delay, represented by the lag projection, would match the dominant spatiotemporal pattern, the SMLV-to-FPN. To test this hypothesis, we computed the lag projection of all cortical BOLD time courses from the same 50 subject sample of resting-state scans. We found that the spatial correlation between the time-lag delay map of the SMLV-to-FPN and the lag projection map is strong (r = 0.81), and both exhibit the same direction of propagation ( Figure 4). Thus, the time-lag dynamics described by lag projections map closely onto that described by the SMLV-to-FPN. However, there is a discrepancy between the estimated duration between the two approaches -the estimated duration of the SMLV-to-FPN from the cPCA is ~22s, and the full duration of the lag projection is ~2.5s. This may suggest that the time-lag structure of the SMLV-to-FPN exists at shorter time scales. Note, the lag projection we computed only partially resembles the average lag projection in Mitra et al. (2014) -the differences are due to preprocessing differences, which we discuss in Supplementary Results F.
While lag projections describe short-time scale propagation structures, the QPP is a much longer temporally-extended pattern (>20s). Visual comparison of the spatiotemporal pattern of the QPP (Majeed et al., 2011) with the SMLV-to-FPN revealed a super cial similarity. Thus, we hypothesized that both the QPP and the SMLV-to-FPN describe the same spatiotemporal dynamics. We derived the QPP from a repeated-template-averaging procedure with a commonly used window size (~21s; 30TRs; Youse et al., 2018) on the resting-state fMRI data. We then computed the correlation between the time course of the QPP and the time course of the SMLV-to-FPN. We found that the time courses of the SMLV-to-FPN and QPP were strongly correlated (r = 0.72) at a time-lag of 7TRs (~5s). The similarity in spatiotemporal dynamics between the QPP and SMLV-to-FPN can also be illustrated visually. We visualized the spatiotemporal template of the QPP from the repeated template-averaging procedure, and compared it to the time point reconstruction of the SMLV-to-FPN described above (Movie 2). As can be seen from the visualization, the temporal dynamics of the SMLV-to-FPN overlap signi cantly with the dynamics of the QPP.
During the steady states of SMLV-to-FPN, the distribution of weights is roughly unipolar, meaning it is either all positive or all negative ( Figure 2A). This suggests that the SMLV-to-FPN may track the global mean time course of cortical vertices, otherwise known as the 'global signal'. We found that this is indeed the case -the time course of the SMLV-to-FPN and the global mean time course are statistically indistinguishable (r = 0.97). This would also suggest that the temporal dynamics surrounding the time points before and after the peak of the global signal correspond to the dynamics of the SMLV-to-FPN. We constructed a dynamic visualization of the global signal through a peak-averaging procedure. Speci cally, all BOLD time courses within a xed window (15TRs on each side) were averaged around randomly-sampled peaks (N=200, > 1 standard deviation above the mean) of the global mean time course. Visually comparing the spatiotemporal visualization of the global signal to the SMLV-to-FPN, we found that the temporal dynamics surrounding peaks of the global signal precisely match those of the SMLV-to-FPN (Movie 2).
The temporal dynamics of the FPN-to-DMN largely represents an anti-correlated pattern between the FPN and DMN -i.e. when regions of the DMN exhibit negative BOLD activity, the regions of the FPN exhibit positive BOLD activity (and vice versa). This resembles the "task-positive" (i.e. FPN) vs. "task-negative" (i.e. DMN) anti-correlation pattern originally observed by Fox et al. (2005) and Fransson (2005). We reproduced these results by correlating each cortical vertices' BOLD time course with a seed time course from the left and right precuneus, a key node of the DMN. Note that the same results were observed with a seed placed in the left and right inferior parietal cortex. As expected, an anti-correlated pattern was observed between the FPN and DMN (Figure 4). We compared the precuneus-seed correlation map to the time points of the FPN-to-DMN using spatial correlation. We found that the pattern of correlations in the precuneus-seed map precisely corresponds to the pattern of BOLD activity in the beginning phase of the FPN-to-DMN (r = 0.92, t = 1.8s). Thus, this would seem to suggest that the task-positive vs. task-negative pattern arises from the anti-correlated dynamics between the FPN (task-positive) and DMN (tasknegative) represented by the FPN-to-DMN spatiotemporal pattern.
A similar anti-correlation pattern to the task-positive/task-negative pattern has been observed in the FC gradient literature (Margulies et al., 2016;Vos de Wael et al., 2020), known as the 'principal' or 'primary' FC gradient (PG). In our zero-lag FC topography survey (Figure 1), we computed the PG as the rst component derived from the Laplacian Eigenmaps (LE) algorithm, consistent with Vos de Wael et al. (2020). As opposed to the task-positive/task-negative pattern, the PG exhibits an anti-correlated pattern of spatial weights between the SMLV complex and the DMN (Figure 4). Further, the PG has been referred to as the principal direction of variance in cortical functional connectivity (Margulies et al., 2016). However, the results from both PCA ( Figure 1) and cPCA ( Figure 2) identify the SMLV-to-FPN as the principal direction of variance in cortical functional connectivity, which does not exhibit the anti-correlated pattern between the SMLV complex and the DMN observed in the PG. In fact, none of the three spatiotemporal patterns exhibit an anti-correlated dynamic between the SMLV complex and the DMN.
With no clear correspondence between the PG and the three spatiotemporal patterns, we sought to identify the factors that might explain the uniqueness of the PG. We discovered that the spatial topography of the PG is due to the con uence of two factors: 1) global signal regression and/or time point normalization (i.e., z-scoring or de-meaning without unit-variance scaling), and 2) thresholding of FC matrices. First, as has been previously observed by Liu et al. (2017), regression of the global mean time course, and de-meaning of cortex-wide BOLD values within a time point (i.e. time-point centering) have very similar effects on cortical time series. Implicit in the computation of LE for functional connectivity gradients, as well as other manifold-learning techniques, is a time-point centering step (Ham et al., 2004 see Supplementary Discussion E). This is relevant because the global mean time course precisely tracks the time course of the SMLV-to-FPN (r = 0.96). This would suggest that removal of the global mean time course through global signal regression and/or time-point centering effectively removes the SMLV-to-FPN from BOLD time courses. What is left over is the second most dimension of variance in FC, the FPN-to-DMN. In fact, this would explain the appearance of the task-positive vs. task-negative pattern after global signal signal regression in seed-based correlation analysis (Fox et al., 2005). We tested this possibility by comparing the output of PCA and cPCA with and without a time-point centering preprocessing step ( Figure 4B). Consistent with our hypothesis, PCA of time-point centered BOLD time courses produces a pattern of spatial weights for the rst principal component that resembles the second principal component from PCA of non-time-point centered BOLD time courses (PC2: r = 0.94). Further, the rst complex principal component of cPCA on time-point centered BOLD time courses exhibits a timedelay map that resembles the second complex principal component time-delay map on non-time-point centered data (cPC2:= 0.49 vs cPC1: = 0.10). Note, the correlation between the time-lag maps was computed using a circular correlation coe cient due to the circular properties of the spatiotemporal patterns (e.g. 0 and 2 are identical angles). Thus, at least one effect of time-point centering and/or global signal regression of BOLD time courses is the removal of the rst principal component and/or the SMLVto-FPN from BOLD time courses.
It is the dual effect of time-point centering and FC matrix thresholding that resolves the discrepancy between the DMN-to-FPN and the PG observed in our study. The FC matrix represents Pearson's correlation of BOLD time courses between all pairs of cortical vertices (i.e. correlation matrix). It is standard practice in computation of the PG that a threshold is performed row or column-wise on the FC matrix (e.g. 90th percentile of correlation values within that row) before the creation of an a nity matrix to input to the manifold learning algorithm (Margulies et al., 2016;Vos de Wael et al., 2020). This preprocessing step is intended to remove noisy or arti cial correlation values from the FC matrix. In our survey of zero-lag FC topographies (Figure 1), we applied a 90th percentile threshold across rows of the FC matrix prior to LE. We found that this preprocessing step obscures the relationship between the PG and DMN-to-FPN. In fact, if no thresholding of the FC matrix is performed, the rst eigenmap produced from LE precisely resembles the FPN-to-DMN contrast observed in the second principal component of non-time-point centered BOLD time courses (r = 0.83; Figure 4A and 4C). As one raises the percentile threshold applied to the FC matrix, the spatial weights of vertices within the FPN, DMN and SMLV complex become more uniform, and the spatial weights of the vertices within the FPN fall to zero ( Figure  4C). At the higher end of percentile thresholds (e.g. 90th percentile) a contrast between the SMLV and DMN begins to appear that is almost equally similar to the unipolar contrast of the rst principal component (r = 0.83) and the anti-correlation contrast of the second principal component (r = 0.82).
6. Network-Based Representations of Functional Connectivity FC topographies are low-dimensional representations of zero-lag synchronous relationships among BOLD time courses. In the network or graph-based approach to FC analysis, the unit of analysis is pairwise relationships between BOLD time courses. Rather than reducing pairwise relationships to lowerdimensional representations, the network-based approach analyzes the structure of these relationships.
We sought to identify the degree to which the structure of pairwise zero-lag synchronous relationships between BOLD time courses arises from the shared dynamics of the three distinct spatiotemporal patterns. A network representation of FC was constructed by computing the correlations between all pairs of cortical BOLD time courses to create a correlation or FC matrix ( Figure 5). We compared this FC matrix to a FC matrix that was reconstructed from the three spatiotemporal patterns. Reconstructed cortical BOLD time courses were created from the spatiotemporal patterns by projecting the time courses of each pattern back into the cortical vertex space. A 'reconstructed' FC matrix was computed from these reconstructed time courses in the same manner as the original BOLD time courses. We estimated the similarity between the two FC matrices by computing the correlation coe cient between the lower triangles of each matrix. Despite a larger mean correlation value in the reconstructed FC matrix, we found that the patterns of correlations between the FC matrices were highly similar (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 rst thresholded the FC matrices by setting those correlation values below the top 80% of correlation values to zero. We estimated network communities from the original FC matrix using the Louvain modularity-maximization algorithm with a resolution parameter value of 1. To assess the degree of community structure in the original FC matrix, we computed the modularity value of the partition of vertices into communities from the Louvain algorithm. The modularity value varies from -1 to 1 and represents the ratio of the summed intracommunity correlation coe cients to that expected at random, such that higher values indicate a 'higher quality' partition. The modularity value of the original FC matrix was Q = 0.34. We then examined whether the same community structure of the original FC matrix was present in the FC matrix reconstructed from the spatiotemporal patterns. We assigned the vertices of the reconstructed FC matrix to the community assignments derived from the original FC matrix and re-calculated the modularity value. We found that the modularity value of the community assignments applied to the reconstructed FC matrix was almost as strong as the original FC matrix (Q = 0.29). In other words, the community structure of the original FC matrix is present in the FC matrix constructed from the shared dynamics of the three spatiotemporal patterns.

Discussion
Over the past decade, intrinsic functional brain organization has been characterized by a myriad of analytic methods (Bijsterbosch et al., 2020). This study attempted to synthesize this complex landscape into a set of fundamental patterns and organizational principles. We found three canonical spatiotemporal patterns of BOLD activity that are consistently observed across analytic methods, what we refer to as the 'SMLV-to-FPN', 'FPN-to-DMN' and 'FPN-to-SMLV' spatiotemporal patterns. These spatiotemporal patterns are best represented by time-lag analyses capable of characterizing temporallyextended patterns of BOLD activity. However, we also found that many zero-lag FC analyses are capable of representing parts of this spatiotemporal pattern, albeit in the form of 'static' or 'stationary' snapshots.
Further, we found that a signi cant proportion of cortex-wide FC network structure is explained by the shared dynamics of these three spatiotemporal patterns.
The three spatiotemporal patterns account for a wide-variety of previous ndings in resting-state fMRI.
Previous research has demonstrated that travelling waves or propagatory patterns are ubiquitous features of intrinsic cortical BOLD activity (Gu et al., 2020;Majeed et al., 2011;Mitra et al., 2015). We found that the 'SMLV-to-FPN', representing the dominant axis of variance in cortical BOLD time courses, captures the same pattern of propagatory activity in lag projections originally discovered by Mitra et al. (2014). Further, the 'SMLV-to-FPN' precisely matches the spatiotemporal dynamics of the QPP originally discovered by (Majeed et al., 2011) (though recent ndings suggest there may be more than one QPP, Youse and Keilholz, 2021). Consistent with ndings that the time course of the QPP is closely correlated with the global BOLD signal (Youse et al., 2018), we nd that the time course of the SMLV-to-FPN is statistically indistinguishable from that of the global BOLD signal. Further, the spatiotemporal pattern of BOLD activity around peaks of the global signal precisely matches the spatiotemporal pattern of the SMLV-to-FPN.
Another ubiquitous feature of intrinsic cortical BOLD activity is anti-correlated FC between the DMN and FPN, sometimes referred to as the 'task-positive' and 'task-negative' pattern (Fox et al., 2005;Fransson, 2005). Some have argued that the anti-correlated pattern of FC between the DMN and FPN is arti cially introduced by regression of the global BOLD time course (Murphy et al., 2009). Others have argued that features of the anti-correlated pattern are independent of global signal regression (Fox et al., 2009) Despite their independence, the FPN-to-DMN shares at least two features in common with the SMLV-to-FPN. First, the all-positive contrast between SMLV (high amplitude) and DMN (low amplitude) in the beginning phase of the SMLV-to-FPN has a strong degree of spatial similarity with the anti-correlated contrast in the later phase of the FPN-to-DMN (r = 0.7). The discrepancy between the two patterns is the lack of high amplitude values in the SMLV in the beginning phase of the FPN-to-DMN. However, the SMLV does appear at a later phase of the FPN-to-DMN, during the propagation of the SMLV to the FPN. In the SMLV-to-FPN, SMLV activation is simultaneous with activation of the FPN. In the FPN-to-DMN, SMLV activation occurs slightly after that of the FPN. This suggests that the primary difference between the two spatiotemporal patterns is a differential onset time of BOLD activity within the SMLV. Second, the direction of propagation is consistent between the two spatiotemporal patterns. Both patterns exhibit propagation from the SMLV towards the FPN. However, as noted above, the propagation of BOLD activity from the SMLV to FPN occurs later in the FPN-to-DMN than the SMLV-to-FPN. The similarity between the two spatiotemporal patterns suggests that they may share common physiological or neuronal mechanisms.
The temporal dynamics of the three spatiotemporal patterns can be divided into steady-states and propagation events. Steady-states are periods of relative stability of BOLD activity across consecutive time points of the spatiotemporal pattern. These steady-states are interrupted by propagation events exhibiting rapid changes in BOLD activity between consecutive time points. The three spatiotemporal patterns can be distinguished by their ratio of time points classi ed into propagation events versus steady-states. The FPN-to-DMN and FPN-to-SMLV spend most of their time in steady-states with very few time points exhibiting propagation. The SMLV-to-FPN exhibits a smooth, spatially continuous propagation of BOLD activity and contains a higher ratio of propagation time points versus steady-state time points than the FPN-to-DMN and FPN-to-SMLV.
Analysis of the temporal evolution of the three spatiotemporal patterns revealed a surprising pattern: the same spatial topography of BOLD activity can appear across more than one spatiotemporal pattern. For example, the spatial topography of BOLD activity that occurs during the propagation events of the SMLVto-FPN and FPN-to-DMN appears in the steady-state time points of the FPN-to-SMLV. Comparison of the temporal evolution between the FPN-to-DMN and SMLV-to-FPN revealed another striking observation -the dynamics of the FPN-to-DMN and SLMV-to-FPN appear to be mirror images of one another. Speci cally, the steady-state time points of the FPN-DMN resemble the propagation time points of the FPN-to-SMLV, and vice versa. However, these spatiotemporal patterns exhibit opposite directions of propagation: propagation of BOLD activity from the FPN to DMN and SMLV to FPN in the FPN-to-DMN, and propagation of BOLD activity from the DMN to FPN and FPN to SMLV in the FPN-to-SMLV.
Our ndings suggest that intrinsic functional brain representations are markedly consistent across analytic methods. This fact is even more surprising considering the wide array of mathematical and statistical assumptions of the analytic methods surveyed in this study. This observation does not imply that all analytic methods we surveyed produce the same insights. Each method affords a unique perspective on the spatial and temporal properties of the three spatiotemporal patterns discovered in this survey. Further, the level of functional brain organization explored in this study is an important quali cation of our ndings. This study has shown consistency at the level of widely-distributed cortical BOLD activity patterns among analytic approaches. Most of the analysis approaches we surveyed can reveal ner-grained spatial insights at higher component or cluster numbers. Thus, we do not expect the same consistency in analytic approaches at ner-grained levels of analysis (e.g. an ICA solution of 50 components vs. a PCA solution of 50 components). Despite this quali cation, the consistency of widelydistributed representations of cortical functional brain organization, in the form of three common spatiotemporal patterns, attests to the dominance of these structures at higher spatial levels of analysis.
Another limitation of this study was the exclusion of subcortical areas from our analyses. This study was primarily interested in cortical functional brain organization. However, previous research suggests that the large-scale spatiotemporal patterns observed in this study have signi cant subcortical contributions as well (Youse and Keilholz, 2021;Youse et al., 2018).
The appearance of these three large-scale spatiotemporal patterns across a wide-variety of analytic methods leads to the question of what these patterns may represent in terms of neuronal or hemodynamic processes. The similar time-scales and temporal dynamics between these three spatiotemporal patterns suggest they may emerge from similar mechanisms. Momentary uctuations in arousal and/or vigilance are known to be related to the global BOLD signal (Liu et al., 2015(Liu et al., , 2018Schölvinck et al., 2010), BOLD propagation dynamics (Gu et al., 2020), and the anti-correlated dynamic between the FPN and DMN (Kucyi et al., 2020). All three of these features of intrinsic BOLD activity were found to be related to one or more of three spatiotemporal patterns. This may suggest that the three spatiotemporal patterns emerge from neurometabolic processes associated with momentary uctuations in arousal. Future research may be directed towards a more complete understanding of the common or distinct neuronal or physiological mechanisms that give rise to these spatiotemporal patterns.

Resting-State fMRI Data
Our study utilized resting-state fMRI scans from the Human Connectome Project (HCP) S1200 release (Van Essen et al., 2013). Participants were non-related, healthy young adults (ages 22-37). Resting-state fMRI data was collected over two consecutive days for each subject and two sessions, each consisting of two 15 minute runs, amounting to four resting-state scans per subject. 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 min scan from a random sample of participants (n=50; 21 males) on the rst 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 results were not biased by acquisition from any given phase encoding direction. We chose a single 15 min scan per participant to ensure that the phase encoding/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 that had been previously preprocessed with the HCP's ICA-based artefact removal process (Smith et al., 2013) to minimize effects of spatially structured noise in our analysis. All brain-imaging data were acquired on a customized Siemens 3 T 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 s temporal resolution. Further details of the data collection and preprocessing pipelines of the HCP can be found elsewhere (Smith et al., 2013;Van Essen et al., 2013). Informed consent was obtained from all subjects. All methods were carried out in accordance with relevant guidelines and the University of Miami Institutional Review Board approved the study.

Resting-State fMRI Preprocessing
Resting-state fMRI scans were spatially smoothed with a 5mm FWHM kernel using the surface-based smoothing algorithm in Connectome Workbench Version 1.4.2. Resting-state fMRI signals from each vertex were then temporally ltered to the conventional low-frequency range of resting-state fMRI studies using a Butterworth bandpass zero-phase lter (0.01-0.1Hz). 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 (Dale et al., 1999). This step down-sampled the total number of vertices in the left and right cortex to 4800 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 subject resting-state scans.

Description of Zero-lag FC Analyses
Our study distinguished between two different descriptions of intrinsic functional brain organizationzero-lag synchrony between brain regions, and time-lag synchrony between brain regions. Zero-lag synchrony is de ned as in-phase or anti-phase statistical dependence between intrinsic BOLD signalse.g. the correlation between two intrinsic BOLD signals with no time-lag. Following the standard terminology of the functional magnetic resonance imaging (fMRI) literature, we refer to zero-lag synchrony among intrinsic BOLD uctuations as 'functional connectivity' (FC) (Reid et al., 2019). FC between cortical brain regions organize into global, cortex-wide patterns, referred to as 'FC topographies'.
All analyses were conducted so as to be consistent as possible with previous studies. For some of these analyses, results were compared with and without global signal regression. Global signal regression was Varimax rotation of principal components: consists of an orthogonal rotation of the principal component spatial weights, such that the simple structure of the spatial weights are maximized.
Simple structure is de ned such that each vertex loads most strongly one component, and weakly on all others. We used the implementation of varimax rotation in the FactorAnalyzer Python package (https://github.com/EducationalTestingService/factor_analyzer).
Laplacian Eigenmaps (spectral embedding): is a nonlinear manifold learning algorithm popular in the FC gradient literature (Vos de Wael et al., 2020). The input to the Laplacian eigenmaps algorithm was the vertex-by-vertex cosine similarity matrix (Margulies et al., 2016), representing the similarity in BOLD time series between all cortical vertices. Of note, cosine similarity is equivalent to Pearson correlation in mean-centered and unit normalized time series (i.e. z-score normalization), as was the case with our data. Laplacian Eigenmaps performs an eigendecomposition of the transformed similarity matrix, known as the normalized Laplacian matrix. We also computed Laplacian Eigenmaps 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 scikitlearn (V0.23) Python package, and details can be found at (https://scikitlearn.org/stable/modules/generated/sklearn.manifold.SpectralEmbedding.html).
Spatial and temporal independent component analysis (ICA): estimates linearly mixed, statistically independent sources from a set of input variables. In the case of spatial ICA, principal component axes derived from PCA of the time point-by-time point covariance matrix are rotated to enforce statistical independence in the spatial domain. In the case of temporal ICA, principal component 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 principal component solution for both temporal and spatial ICA. We used the FastICA algorithm implemented in the scikit-learn (V0.23) Python package. Details can be found at (https://scikitlearn.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 dominant spatiotemporal patterns -SMLV, FPN and DMN. We chose seeds in the somatosensory cortex (SMLV), precuneus (DMN), and supramarginal gyrus (FPN) (Supplementary Figure B2). The spatial outline of the SMLV, DMN and FPN for guiding the selection of seed regions were determined through a k-means clustering analysis of the temporally-concatenated group time series with cortical vertices as observations and BOLD values at each time points as input variables (i.e. 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 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.
Co-activation pattern (CAP) analysis: Three CAP analyses were performed for the same three seed regions used in the seed-based regression analysis. CAP analyses rst identify time points with the highest activation for a seed time course. Consistent with previous studies (Liu and Duyn, 2013b), we chose the top 15% of time points from the seed time course. The BOLD values for all cortical vertices in the top 15% time points 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.
Hidden Markov modeling (HMM): is 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/deactivations, somewhat similar to CAPs. To avoid over tting and to reduce noise in the high-dimensional input data, we conducted a PCA of the cortical BOLD time series. The rst 100 principal component 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 principal component (N = 100), and a covariance matrix between the 100 principal component 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 time point conditional on a hidden brain state (i.e. 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 (V0.2.5) (https://github.com/hmmlearn/hmmlearn).

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 Laplacian Eigenmaps, 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 ne-grained distinctions of the data description -i.e. 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 (Bzdok et al., 2016;Schaefer et al., 2018).
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 rst T dominant axes of variation (i.e. principal components) of the correlation matrix formed between all pairs of cortical BOLD time series. Speci cally, we examined the drop-off in explained variance (i.e. eigenvalues) associated with neighboring principal components, a procedure known as Catell's scree plot test (Cattell, 1966). 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 principal component solution of one and three ( Figure 1A). We chose the higher-dimensional solution of three components. Note, the elbow in explained variance after three components was independent of the functional resolution (i.e. 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 preprocessing pipeline). Thus, three latent dimensions were estimated for all dimension-reduction algorithms, and three hidden states were estimated for the HMM.

Time-lagged Analyses of Resting-State fMRI Data
Quasiperiodic Pattern and Lag Projections Time-lag analyses capture relationships between two or more time series at past or future lags of the time series. We refer to representations of time-lag relationships between cortical time series as 'spatiotemporal patterns'. There are two widely-used algorithms for the study of spatiotemporal patterns in BOLD signals: 1) interpolated cross-covariance functions (Mitra et al., 2014(Mitra et al., , 2015 for the detection of lag/latency projections (~0-2s) and 2) a repeated-template-averaging algorithm of similar spatiotemporal segments (Majeed et al., 2011) for detection of the QPP (~20s).
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 de ned as the peak of their lagged cross-covariance function. Lag projections are derived as the column average of the pairwise time-lag matrix between all cortical vertex time courses.
To estimate the QPP, the template-autoregressive matching algorithm of Majeed et al. (2011) was used.
The algorithm operates in the following manner: start with a random window of BOLD TRs, compute a sliding window correlation of the window across the temporally concatenated group data at each time point, and then average this segment with similar segments of BOLD TRs (de ned 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 nal sliding window correlation time series. The nal sliding window time series is the same length as the original subject 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 modi ed from the C-PAC toolbox (https://fcp-indi.github.io/). Consistent with previous studies (Majeed et al., 2011;Youse et al., 2018), 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 10 times. The template with the highest average sliding window correlation time series across the 10 runs was chosen as the nal result.

Complex Principal Component Analysis
To tie the three dominant FC topographies to spatiotemporal processes in the cortex, we used a simple modi cation of PCA for detection of spatiotemporal patterns. Speci cally, we apply PCA to complex BOLD signals obtained by the Hilbert transform of the original BOLD signals. We refer to this analysis as complex PCA (cPCA). This technique has been referred to as complex Hilbert empirical orthogonal functions in the Atmospheric and Climate sciences literature (Horel, 1984).
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=(i.e. 90 degree), 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 courses -i.e. the time-lag at which the correlation is maximum. In the same manner, cPCA produces complex spatial weights for each principal component that can give information regarding the time-lags between BOLD time courses. In the same manner that a complex signal is composed of real and imaginary signal components, the complex principal component's spatial weights are composed of real and imaginary components. and the phase delay of the spatial weights.  To examine the temporal progression of each complex principal component, we sampled the reconstructed BOLD time courses from each complex principal component at multiple, equally-spaced phases of its cycle (N=30; Figure 2). For each complex principal component, the reconstruction procedure was as follows : 1)  Form   information regarding the temporal dynamics of the spatiotemporal pattern. C) The same twodimensional slices of each spatiotemporal pattern in the 3-dimensional principal component space colored according to their cluster assignment by a k-means clustering algorithm. K-means clustering was used to identify recurring spatial patterns of BOLD activity across time points of the three spatiotemporal patterns. Six clusters were estimated. D) The cluster assignments (color) by time (y-axis) of each spatiotemporal pattern (x-axis). Note, that the same cluster assignment can occur across more than one spatiotemporal pattern. E) The cluster centroids from the k-means clustering algorithm, corresponding to the average spatial pattern of BOLD activity for the time points that belong to that cluster. Note, the cluster centroids of the rst two clusters are mean-centered versions of the original unimodal (all-positive or all-negative) steady-state of the SMLV-to-FPN, as z-score normalization of the time-points across vertices was performed beforehand.

Figure 4
Similar Time-lag Dynamics between SMLV-to-FPN and Lag Projection. Comparison between the time-lag delay map of the SMLV-to-FPN (left) and the average lag projection (right). As in Figure 2, the SMLV-to-FPN time-lag delay map represents the phase delay (in radians) of cortical BOLD time series. The lag projection map represents the average time-lag delay (in seconds) between each vertex of the cortex. The spatial correlation between the SMLV-to-FPN time-lag delay map and lag projection is r = 0.81, indicating a strong similarity in time-lag dynamics.

Figure 5
The Task  Modules with < 100 vertices were considered 'junk' modules and are placed in the upper right hand corner of each sorted correlation matrix. 'Junk' modules were not represented in the module assignments (bottom). Beside the title in parentheses of each correlation matrix is the modularity value. The modularity value represents the degree of community structure for the modular partition of that correlation matrix (i.e. degree of intra vs. inter module correlation weights). The modular structure estimated from the original correlation matrix (Q = 0.34) retains a high degree of modularity when used to partition the reconstructed correlation matrix (Q =0.29), indicating the three spatiotemporal patterns explain a signi cant amount of the modular structure in the original correlation matrix. Further, despite a higher mean value of correlations in the reconstructed correlation matrix, the pattern of correlations between the two correlation matrices is highly similar (r = 0.77).