Bursty properties revealed in large-scale brain networks with a point-based method for dynamic functional connectivity

The brain is organized into large scale spatial networks that can be detected during periods of rest using fMRI. The brain is also a dynamic organ with activity that changes over time. We developed a method and investigated properties where the connections as a function of time are derived and quantified. The point based method (PBM) presented here derives covariance matrices after clustering individual time points based upon their global spatial pattern. This method achieved increased temporal sensitivity, together with temporal network theory, allowed us to study functional integration between resting-state networks. Our results show that functional integrations between two resting-state networks predominately occurs in bursts of activity. This is followed by varying intermittent periods of less connectivity. The described point-based method of dynamic resting-state functional connectivity allows for a detailed and expanded view on the temporal dynamics of resting-state connectivity that provides novel insights into how neuronal information processing is integrated in the human brain at the level of large-scale networks.


Supplementary Methods
Justifying the number of chosen clusters (k) using an independent dataset.
It first needs to be said that in any type of clustering analysis performed on fMRI data when the number of clusters (k) has to be empirically determined, there are two problems which might occur. (1) The value of k is set too low, resulting in a too small number of clusters being defined which would potentially miss out on some aspects of the brain dynamics.
(2) The value of k is set too high which results in too many clusters being defined, leading to an "over-splitting" of the fMRI data. The case of over-splitting could prove to be highly detrimental, as an additional split of a cluster could possible result in that the magnitude of correlation between the BOLD signal intensity values in a cluster switch from being positive to negative (or vice versa), inducing spurious connectivity. Thus, providing as solid experimental support as possible for the chosen value of k is of paramount interest. As an alternate route, the analysis can be extended to include a validation of the results for a range of values of k to show a robustness of an effect that is largely independent on the size of k.
In this work we have both developed a strategy that includes a separate, independent data to validate our choice of k as well as validating our results for different choices of k. With respect to the first approach, we performed additional analyses on the second dataset from the human connectome dataset (i.e. the so called "LR" dataset from the first resting scan in which the direction of the phase encoding gradient (left-right versus right-left) was flipped with respect to the "RL dataset" used in the main analysis, see Ugurbil et al., 2013 andvan Essen et al., 2012 for further details) obtained from the same subjects. Identical pre-processing and post-processing steps were applied. A schematic diagram of the work-flow used here is presented in Supplementary figure 3A. Prior to clustering, the complete dataset was evenly split into three blocks of data. The time-points (or equivalently image volumes) 1 to 400 for all subjects were placed in the first block (here named A), 401 to 800 were placed in the second block (B), and time-points 801 to 1200 was allocated to the third block (C).
For all three data blocks, a clustering using k-means algorithm was computed for k =2 to 20. In the next step, for each k, we paired the clusters between each combination of block dataset (i.e. possible parings in this case are AB, AC and BC) based on the average BOLD activity for each node in each cluster (see Supplementary Figure 3B that shows an example for k = 8). The paring was achieved by computing distance matrices between each data block's clusters using the averaged taxicab distance (see equation (1)) for the average BOLD activity of each node. This procedure yielded a k x k distance matrix for each unique paring of block dataset. Subsequently, the so called Hungarian Algorithm (Kuhn 1955;Munkres 1957) was used to find the optimal cluster pairing for each distance matrix. The Hungarian Algorithm was applied to the clustering paring distance matrices and the output is a k x 2 solution matrix for each of the block dataset pairings. Since there are three possible block dataset pairings, it is conceivable that the optimal solutions for the different block dataset pairings are not transitive and thereby resulting in incompatible clustering solutions. This possibility was accounted for by checking if the pairing of clusters in the block dataset pairing AC was identical to the pairing of clusters that could be derived for AC by using AB and BC. After checking for transitivity among possible pairing of data blocks, nine of the tested k values had transitive solutions (k = 2, 3, 4, 5, 6, 7, 8, 11 and 12). This suggests that any choice of k that has a transitive solution is to some extent reproducible in terms of k-means clustering with similar spatial patterns independently in different parts of the second data set. Hence, by the validation of our approach using an independent data set, we believe that our choice of setting k=8 is justified, since it presented us with similar global spatial patterns for the three different parts of an independent dataset. Additionally, our key results were replicated using other choices of k (for results setting k to 5 and 12, see results shown in Supplementary Figure S7). It should be noted that our justification of choice of k that includes the splitting of an independent dataset into three parts relies on the assumption that the number of derived states in the same subjects will identical. .   Fig. S1. An illustration of the concept of estimating the co-variance between ROI signal intensity time-courses based on non-neighboring time-points in different clusters and the creation of functional connectivity time-series. In this example, the BOLD signal intensity time courses from two seed regions located in the ACC (red line) and insula (blue line) are first extracted (A). For simplicity, we here show for a single subject and, in contrast to Figure 1, how the PBM method effects a single connection between two brain regions, rather than the entire graph. The signal intensity time-series are subsequently divided into different clusters using the k-means clustering algorithm (see Figure 1 and Methods section) which are here shown for the case of k=8 using a coloring scheme displayed above the signal intensity time-series in panel (A). In a subsequent step, the co-variance is computed for all time-points that are allocated for each cluster (B). The final step is shown in panel (C) where a functional connectivity time-series between the ACC and the insula is created based on the co-variance computed in (B) for the eight different clusters. An additional descriptive illustration that includes the clustering step is given in Figure 1 and a detailed account is given in the methods section.    The eight states/s-graphlets derived from the k-means clustering are here shown using the spring-embedded Kamada-Kawai representation. Only the top 5 percent of all edges, similar to Figure 2D, were included in the plots. Fig. S6 No state corresponds exclusively to movement. Distribution of timepoints assigned to each state at time-points for which FD>0.5, indicating high micro-movements. There is no state which is exclusively found to be associated with estimated micro-movement artifacts.  Replication of the burstiness of brain connectivity for different choices of k (k = 5 and 12 in the k-means algorithm) used to derive s-graphlets. Similar to the results shown in Figs. 6C-F, the percentage of edges that were significantly bursty or tonic/periodic (p ≤ 0.05, two tailed) is shown using two different edge thresholds (withholding either the top 5 or 10 percent of all edges).