One of the most common techniques for examining time-varying rsfMRI is a sliding window (SWC) correlation analysis (Chang and Glover, 2010; Keilholz et al., 2013; Preti et al., 2017). Based on widely-used methods for calculating average functional connectivity (FC) as the correlation between brain areas over the course of the entire scan, in SWC, the fMRI timeseries is segmented using a time window of length w (typically 30–60s) and the correlation between pairs of brain regions is calculated for each window at every timestep.

Phase synchrony (PS) is similar to the sliding window approach, except the window length is reduced to a single timeframe by using Hilbert (Glerean et al., 2012; Cabral et al., 2017) or wavelet (Chang and Glover, 2010; Yaesoubi et al., 2015) transformations to obtain the instantaneous phase of the timeseries at each timepoint, then the PS between pairs of brain regions can be calculated. As with the sliding window approach, k-means clustering (Yaesoubi et al., 2015; Cabral et al., 2017) can then be applied to identify temporal dynamics.

Another approach to the analysis of single timeframe rsfMRI is to examine co-activation patterns (CAP) (Liu et al., 2013; Chen et al., 2015). In contrast to the phase synchrony approach, CAP analysis identifies simultaneous occurrence of BOLD signal peaks or troughs in different brain regions independent of signal phase. The relationship between the BOLD signal and neural activity has been shown to arise from the temporally sparse events which form the basis of CAP analysis (Zhang et al., 2020). CAP analysis typically utilizes k-means clustering to identify different CAPs and this allows temporal dynamics to be identified and potentially compared to other time-varying rsfMRI methods including SWC and PS.

Finally, one alternative to temporal clustering of single timeframes is to identify spatiotemporal patterns which occur repeatedly over sequences of frames (Majeed et al., 2011). Quasi-periodic patterns (QPPs) observed through this type of analysis represent a known dynamic feature of rsfMRI (Yousefi et al., 2018; Bolt et al., 2021) and have been shown to correlate with neural activity (Thompson et al., 2014). The strongest QPP typically consists of a sequence displaying a transition between a period of strong activation of default-mode network (DMN) and deactivation of sensory and attention networks to a period of DMN deactivation coupled with activation of sensory and attention networks (Abbas et al., 2019a; Yousefi and Keilholz, 2021).

We examine the relationship between brain states formed using SWC, PS, or CAP analysis, and investigate how each one classifies QPPs observed in the dataset. Our findings highlight commonalities and disparities across the methods and provide some insight into the sensitivity of various approaches.

The minimally preprocessed grayordinate and FIX de-noised rsfMRI of the HCP S900 release (Van Essen et al., 2012) consisting of 817 individuals with four complete rsfMRI scans (1,200-timepoints, TR = 0.72s) were downloaded and preprocessed for previous studies (Yousefi et al., 2018; Keilholz et al., 2020). In brief, each timeseries was demeaned and bandpass filtered (0.01–0.1 Hz), and global (white matter, cerebrospinal fluid (CSF), and gray matter) signals were regressed. The spatial dimension was reduced to 360 cortical parcels (Glasser et al., 2016) and each parcel’s timeseries was z-standardized. The first scan of the first day for each subject was used in all calculations of time-varying FC.

Square 60s windows (83 timepoints), overlapping for all but a single timepoint, were utilized for the sliding-window analysis (Shakil et al., 2016). The Pearson correlation for each parcel-pair was calculated and Fisher transformed to normalize the variance during each 60s time-window from the first timepoint and at each ensuing time-step (1TR = 0.72s) with windows overlapping except for a single timepoint. A symmetrical 360 × 360 matrix of pairwise values for all brain regions was obtained for each window (1,117 per scan = 1,200 total timepoints–window length) and a vector-by-time matrix containing a single value for each of the 64,620 parcel-pairs (equivalent to the lower triangle of the 360 × 360 matrix) at each timepoint served as the input for the clustering algorithm to define brain states.

The instantaneous phase was calculated at every timepoint by taking the Hilbert transform for each parcel. The phase synchrony between each parcel-pair was then computed as the cosine of the difference in phase angles at each timepoint. A vector-by-time matrix of 64,620 pairwise phase synchrony values was obtained for each timepoint and used as the input for the clustering algorithm.

The original CAP analysis applied an activation threshold (often top 15% or z-score > 1) to a seed region and averaged BOLD frames across all suprathreshold timepoints (Liu and Duyn, 2013). Further work extended CAP analysis to a more data driven, whole-brain approach by performing k-means clustering directly on the BOLD timeseries data (Liu et al., 2013). For the present study we performed k-means clustering of the 360 parcels at each timepoint. In order to display cluster centers in a comparable format to the other methods, after clustering, the 360 elements of each centroid were multiplied with themselves generating a 360 × 360 outer-product matrix for each cluster center.

Clustering was performed in MATLAB using the k-means function with k = 5 as the cluster number, correlation distance as the distance metric, and 30 replications per run for each of the SWC, PS, and CAP methods. Each of these methods generated a 2D vector-by-time matrix of values for all 360 parcels (CAP) or 64,620 parcel-pairs (SWC/PS) at every timepoint (1,200/scan for PS/CAP) or every time window (1,117/scan for SWC). For each dynamic analysis method, the vector by time matrices were initially clustered for each individual scan. The cluster centroids (means) from all individual scans were then clustered by k-means using the same parameters to generate group-level cluster centroids. The group-level centroids were the used to initialize clustering of all timepoints across all individuals. This iterative clustering process has been utilized in previous studies of brain dynamics (Allen et al., 2014, 2018) and greatly reduces the computational time and resources required for large datasets.

Selecting an appropriate number of clusters is a critical part of this type of dynamic analysis. Larger numbers of clusters can potentially provide more information but increasing the cluster number greatly increases computation time and can make interpretation more challenging by increasing complexity and reducing the distinction between clusters. For this study, preliminary results were generated using a range of 3–10 clusters on a subset of 20 subject scans (but otherwise the same methods described above). These preliminary results are displayed in Supplementary Materials for k = 3, 7, and 10 clusters. Based on these initial observations, k = 5 clusters were selected as a sufficient number to provide an effective comparison between the different dynamic analysis methods. While the ideal cluster number may vary for the different approaches, it was important to use a consistent number for this comparative analysis.

The resulting clusters, termed brain states, for each dynamic analysis method were visualized by plotting of the cluster centroids as 360 × 360 pairwise matrices with the parcels grouped into seven cortical functional networks as defined by Yeo et al. (2014). Weighted cluster averages were computed for each method as a comparison to time-averaged functional connectivity by taking the mean of the five cluster centroids weighted by the state occurrence rates. The clusters from each method were also compared against each other using several dynamic metrics, including the occurrence rate of each state, the probability of a transition from one state to each of the other states [after concatenating across subjects and ignoring within-state transitions; method described by Cornblath et al. (2020)], and the mean dwell time spent within one state before a transition to another state. The state composition was also compared for each method to determine whether timepoints identified as belonging to a specific brain state by one method were likely to belong to a specific state identified by one of the other methods.

Quasi-periodic patterns analysis identifies repeating spatiotemporal sequences of a specified duration using a pattern finding algorithm (Majeed et al., 2011). The QPP can be thought of as a sequence of specific brain states that repeats over time. For this analysis, a duration of 33 timepoints (24s) was used based on previous work showing prominent QPPs with a duration of 20-24s occurring in resting-state data (Abbas et al., 2019a). Over the full dataset a total of 6,717 occurrences (defined as 33-timepoint sequences with a correlation of r ≥ 0.3 to the QPP template) of the strongest QPP were identified. The QPP was used to compare across dynamic analysis methods by determining the probability of each of the 33 timepoints in the QPP sequence being identified as belonging to a specific brain state using each method.

Several comparisons between analysis methods were tested for statistical significance using MATLAB functions. Permutation tests were performed to quantify the significance of sequential state transitions differing from chance in a method similar to the one described by Cornblath et al. (2020). A total of 100,000 permutations were used to calculate the probability distribution of transitioning to the next (different) state from the current state given the overall occurrence rates of each state with Bonferroni correction for multiple comparisons. The Kruskal-Wallis test was performed to quantify differences in the distribution of subject-mean dwell times across states within and between each of the three analysis methods. The Mann-Whitney test was performed for further pairwise comparisons of dwell times using Bonferroni correction for multiple comparisons. Chi-square goodness of fit tests were used to test for overlap in the composition of states identified by the different analysis methods. The occurrence rates of each state over all 980,400 timepoints were input as the expected frequency and compared against observed frequencies to determine overlap in state composition. The same Chi-square procedure was used to quantify the likelihood that timepoints belonging to a QPP would be clustered into a specific state.

Sliding window connectivity may still be advantageous relative to time-averaged functional connectivity for comparisons between groups (Damaraju et al., 2014; Hindriks et al., 2016; Menon and Krishnamurthy, 2019) and window lengths shorter than the 60s used for this study may produce results closer to the PS method. However, PS and CAP methods do provide noteworthy advantages for dynamic analysis. In additional to capturing information at shorter time scales, each can potentially be performed without an arbitrary window length selection. While the CAP analysis for this study did not employ an activation threshold, performing clustering directly on the BOLD data, seed-based CAP analysis commonly utilizes an arbitrarily selected threshold that could influence the results (Liu et al., 2013; Chen et al., 2015; Gutierrez-Barragan et al., 2019; Zhang et al., 2020). Likewise, PS does not require the investigator to set any specific parameters, but the data is commonly bandpass filtered (Cabral et al., 2017) into frequency bands (typically 0.01–0.1 hz or lower) which is likely to reduce sensitivity to dynamics at higher frequencies. CAP analysis also offers the advantage of reduced dimensionality and computation time given the clustering is performed on data from N parcels rather than (N-1)*N/2 parcel pairs. However, Cabral et al. (2017) have shown that PS (and potentially SWC) analysis can be effectively performed on the leading eigenvector of the parcel pairs, reducing the dimensionality to the same number N and producing similar results.

Investigators may also wish to consider performing QPP analysis for examining brain dynamics as this method has demonstrated clear neural correlates (Thompson et al., 2014), reproducibility across individuals (Yousefi et al., 2018) and across species (Belloy et al., 2018a), and the potential to differentiate patient groups from controls (Abbas et al., 2019b). Further, as PS and CAP analysis have now been shown to be sensitive to QPPs, it is possible that regression of one or more QPPs could benefit a PS or CAP analysis by uncovering additional features that might be obscured by the widespread spatial changes of the QPP (Belloy et al., 2018b; Yousefi and Keilholz, 2021). Meanwhile, for studies explicitly seeking to identify dynamic changes occurring on slower time scales, the SWC method could be advantageous specifically because it is less sensitive to QPPs that could potentially confound the hypothesized differences. Indeed, there is evidence suggesting that SWC methods featuring window lengths ≥30s perform better than shorter windows or single timepoint methods at segmenting different cognitive tasks (Xie et al., 2019).

Finally, there is recent evidence that the correlation between BOLD signal and underlying neural activity arises specifically from the type of peak amplitude events utilized for CAP analysis (Zhang et al., 2020), potentially making CAPs the most sensitive method for studying neural brain dynamics. However, more research is certainly needed to elucidate the relationship between BOLD signal dynamics and neural activity.

One significant limitation of this study is that only k = 5 clusters/brain states were fully tested. Choosing an appropriate number of states for a particular study is important and can pose a challenge as there is not a well-established number of brain dynamics features and no perfect method for determining the optimal number to select. Some studies (Allen et al., 2014, 2018; Damaraju et al., 2014; Gutierrez-Barragan et al., 2019) have employed the “elbow method” in which the ratios of within-cluster distances to between-cluster distances are plotted for a series of k cluster numbers (see Supplementary Figure 11) and a number k is selected which protrudes from the plot such that a greater contribution to optimizing the clustering is apparent for k than for k + 1. However, the number of clusters featuring optimized distances may not always be the optimal number for investigating brain dynamics and ideally studies will investigate a range of k clusters as shown here in the Supplementary Materials for a subset of the data (N = 20). For this study, the subset results do not indicate an optimal cluster number for any method using the elbow criterion but do appear consistent across a range of 3–10 clusters with each analysis method displaying common features and clustering QPP sequences in a consistent manner.

An important methodological consideration for this type of analysis is the distance metric employed for k-means clustering. While the correlation distance was used for the present study, some have suggested that L1 distance may offer advantages for this type of analysis (Allen et al., 2014). However, in performing this study the L1 distance was found to require greater computation time than the L2 or correlation distance metrics, and in limited testing the correlation distance appeared to offer the most qualitatively consistent results.