1. 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. Seed-based 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 drop-off 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 first three principal component maps computed from PCA. The first three principal components represent the top three dimensions of variability across cortical BOLD time series. The first principal component represents the most dominant/leading axis of variance across cortical BOLD time series. The first principal component explains over 20% (R2 = 20.4%) of the variance in BOLD time series, greater than three times the variance explained by the second (R2 = 6.8%) or third principal component (R2 = 6.1%). As can be observed from Figure 1A, each FC topography exhibits strong similarities with one or more of the first 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 first principal component is distinguished by its globally positive topography, with positive values (or negative, due to the sign ambiguity of PCA) in the SMLV and less positive values in the DMN. The second principal component is distinguished by negative values in the DMN and positive values in the FPN. The third principal component is distinguished by positive values in the SMLV and DMN, and negative values in the FPN. Detailed description of the output of dimension reduction analyses, as well as seed-based regression and CAP analyses are provided in Supplementary Results A and Supplementary Results B, respectively.
2. Three Dominant Spatiotemporal Patterns in Intrinsic BOLD Fluctuations
All FC topographies in our survey were found to resemble the first 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 coefficient). Thus, FC topographies are unable to represent correlations between cortical BOLD time series at any time-lag. Time-lag relationships may reflect cortical propagation patterns or rapidly-changing sequences of cortical activity patterns. We refer to these time-lag structures as ‘spatiotemporal patterns.’ Examples of such spatiotemporal patterns include short time-scale lag projections (Mitra et al., 2014, 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 modification of PCA for detection of time-lag relationships between cortical BOLD time series. Specifically, 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 first three complex principal components from cPCA was similar to that observed between the first three principal components from PCA: component 1 (R2 = 21.4%), component 2 (R2 = 6.8%) and component 3 (R2 = 5.7%). Associated with each complex principal component is a time-lag delay map, reflecting 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 first 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 first complex principal component as the ‘SMLV-to-FPN’ spatiotemporal pattern. Because the explained variance of the first 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 fluctuations. 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 ‘FPN-to-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 first principal component (PC1) occurs within the first 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 first 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 finding suggests that the three dominant FC topographies are ‘static’ or ‘stationary’ representations of three temporally-extended, dynamic patterns of BOLD activity.
3. 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 first 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 benefit 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 SMLV-to-FPN, FPN-to-DMN and FPN-to-SMLV vary strongest along the first, 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 first principal component axis with smaller scores on the first and third principal component axes. This simply reflects 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 first 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 five 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-flip - i.e. flipped 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.
4. 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-flipped steady-states of the three spatiotemporal patterns. The first 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 five 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 five 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 five and six) and the propagation events (clusters three and four), and vice versa.
5. 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 superficial 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; Yousefi 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 significantly 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. Specifically, all BOLD time courses within a fixed 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 (task-negative) 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 first 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 confluence 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 first principal component that resembles the second principal component from PCA of non-time-point centered BOLD time courses (PC2: r = 0.94). Further, the first complex principal component of cPCA on time-point centered BOLD time courses exhibits a time-delay 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 coefficient 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 first principal component and/or the SMLV-to-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 affinity 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 artificial 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 first 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 first 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 lower-dimensional 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 coefficient 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 first 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 intra-community correlation coefficients 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.