Transient brain activity disentangles fMRI resting-state dynamics in terms of spatially and temporally overlapping networks

Dynamics of resting-state functional magnetic resonance imaging (fMRI) provide a new window onto the organizational principles of brain function. Using state-of-the-art signal processing techniques, we extract innovation-driven co-activation patterns (iCAPs) from resting-state fMRI. The iCAPs' maps are spatially overlapping and their sustained-activity signals temporally overlapping. Decomposing resting-state fMRI using iCAPs reveals the rich spatiotemporal structure of functional components that dynamically assemble known resting-state networks. The temporal overlap between iCAPs is substantial; typically, three to four iCAPs occur simultaneously in combinations that are consistent with their behaviour profiles. In contrast to conventional connectivity analysis, which suggests a negative correlation between fluctuations in the default-mode network (DMN) and task-positive networks, we instead find evidence for two DMN-related iCAPs consisting the posterior cingulate cortex that differentially interact with the attention network. These findings demonstrate how the fMRI resting state can be functionally decomposed into spatially and temporally overlapping building blocks using iCAPs.


Supplementary Figure 7. iCAPs' sustained-activity time courses
Figure S7: iCAPs' sustained-activity time courses for each subject.  0 20 50  100  150  200  250  0  0 20  200  400  600  0  0 20  200  400  600  800  1  0 20  200  400   There is a range of reasonable values for the number of clusters (indicated by the arrow). We opted for 20 clusters, but limited the detailed analysis to 13. 1 There is a range of reasonable values for the number of clusters (indicated by the arrow). We opted for 20 clusters, but limited the detailed analysis to 13.   Table 1: We compute the average z-score and the total number of voxels occupied in brain areas defined with Talairach Client.

Supplementary Method
Participants. Fourteen healthy volunteers participated in the study. Data acquisitions were obtained with a Siemens 3T Trio TIM scanner, using a 32-channel head coil. The structural images were acquired using a high resolution three- FMRI Data Processing. FMRI data is preprocessed using in-house MATLAB code combined with SPM8 (FIL,UCL,UK) and IBASPM toolboxes 1 . First, fMRI volumes were realigned to the first scan and spatially smoothed with Gaussian filter (FWHM=5mm). We used further motion correction to mark the time points with high frame-wise displacement 2 , and performed cubic spline interpolation around these time points. We did not remove those frames since TA algorithm exploits the continuity of fMRI time courses to deconvolve the effect of the hemodynamic response. Finally, we excluded two subjects with high motion, therefore, the results were obtained from twelve healthy controls in total.
The anatomical images are coregistered onto the functional mean image and segmented (NewSegment, SPM8) for the six different MNI templates. The anatomical automatic labeling (AAL) atlas 3 , composed of 90 regions without the cerebellum, was mapped onto each subject's coregistered anatomical image and further downsampled to match the native space of the functional images. The anatomical atlas is only used to guide the spatial regularization of the Total Activation framework.
Total Activation. The TA framework 4 allows obtaining denoised and well-behaving reconstructions of the activityrelated, activity-inducing, and innovation signals from noisy fMRI measurements by using state-of-the-art regularization that takes the L 1 -norm of the activity-related signal after applying a well-chosen differential operator. The differential operator Δ L includes the inverse operator of the hemodynamic system together with a first-order derivative. The inverse hemodynamic operator is adapted from the formulation in 5 based on the first-order Volterra-series 2 2 2 approximation of non-linear Balloon model 6,7 . TA also combines the temporal with spatial regularization; i.e., we impose mixed-norm constraint to promote (but not enforce) coherent activations inside anatomically defined regions whereas sparse activations across regions.
In essence, TA reverts to convex optimization that combines least-square data fitting with the two regularization terms: where where Δ Lap is the spatial Laplacian operator, λ 1 , 2 are the regularization parameters, ||x|| F is the Frobenius norm and Δ Lap is the Laplacian filter. We use generalized forward-backward splitting 8 , for denoising case also known as parallel Dykstra-like proximal algorithm 9 , to solve the optimization problem in (1). The joint solution is obtained by incorporating the proximal maps of both spatial and temporal regularizations 4,10 .
TA is applied to every subject in his own native functional space, and the results are subsequently normalized (using SPM8) to a common MNI space as to be able to determine group results.

Surrogate data analysis.
For each subject, we generated a surrogate dataset by phase randomizing every voxel time course, and then apply TA with exactly the same settings. In Supplementary Fig. S1, we depict the histograms of the surrogate data and the original data. The histograms of phase-randomised BOLD signal and the original BOLD signal are very similar whereas theys change considerably after TA regularization, which is indicative of the non-random nature of block-type activity in the real data. We plot the log-plot of innovations in order to highlight the difference between the two histogram profiles to accommodate for the large number of data points.

Temporal clustering of transients.
There is no baseline in the innovation signals, that is, every activation is relative with respect to the previous time point. We first separate the innovation signals between activation and de-activation, which are represented by the positive and negative peaks, respectively. Then, we flipped the sign of negatives and concatenated these time points to perform a temporal k-means clustering. In order to determine the time points of interest, we followed a two-step procedure (for space and time) based on the surrogate data that accounts for the subject effect. Specifically, for each subject, we generated a surrogate data by phase randomization of the original BOLD signals of each voxel and run TA on the surrogate data. Then, again for each subject we selected two thresholds as 99% and 1% from the histogram of innovations signals of each subjects surrogate data. Similar to 11 , after thresholding, we kept the voxels that are connected with 26 neighborhood and 6 voxels. Finally, we sum the 'active' time points at each time point and only selected the time points where there were at least 500 active voxels in space (that threshold corresponds to around 4% of the whole volume) 12 . In Supplementary Fig. S2  We then performed k-means clustering using cosine distance as the similarity criteria. The cosine distance between two vectors can be regarded as the non-normalised version of correlation where the mean is kept, as Here, the reasons of using cosine distance are twofolds: 1) it is well adapted for non-negative signals, and 2) it does not modify the baseline. Finally, we empirically set the number of clusters to 20 by using a leave-one-subject-out crossvalidation scheme (see iCAP activation maps in Supplementary Fig. S3-S4 and the cost function in Supplementary   Fig. S9). Inspection of the iCAPs' spatial and temporal patterns, we deduced that the first 13 clusters were particularly stable across subjects and the last 7 clusters had less consistency.
The iCAPs' maps were computed by combining and averaging the clusters of positive and negative (sign i flipped) transients. Since the distribution of the maps are not symmetric, in order to be able to compute z-scores, we subtracted the mode (the maximum value of the histogram) instead of the mean and divide by the standard deviation. The time course of each cluster was computed by back-projecting the iCAPs onto the sustained activity-inducing signals. The back-projection was computed separately for positive and negative weights in order to minimize the effect of spatial linear dependency.
Amount of spatial overlap through Jaccard Distance. The similarity measure was computed via Jaccard distance, which is the total intersection between two binary patterns normalised by each individual pattern; i.e., ratio of intersection to union.
Binarized maps are obtained by thresholding according to z-score ≥ 1.5. Significant intersections are determined as follows: we generate surrogate data by spatially permuting the binary maps, where instead of voxel-wise permutation we exploited a block-permutation to keep the spatial structure (blocks of size 6 × 10 × 5 voxels). Then, we performed a non-parametric test using the maximum statistics and computed the significance p≤0.01. Supplementary Fig. S5 shows the spatial overlap matrix; stars indicate statistically significant connections.

Spatial correlates of traditional DMN
We computed the conventional DMN using seed-region based correlation with the seed in PCC; MNI coordinates (0,-53,26), averaged over a 7 × 7 × 11mm 3 neighborhood. The conventional DMN was derived for each subject and then averaged to obtain the group-level conventional DMN. We computed the spatial similarity between the iCAPs and group-level conventional DMN using cosine distance. The subject-level similarities are shown in Supplementary Fig. 8.

Temporal overlap.
Total and average durations of iCAPs were computed from the normalized iCAPs time courses (absolute |z| ≥ 1, see Supplementary Fig. S8 for the histograms). We counted the number of 'active' iCAPs at each time point regardless of the sign of the weight, and observed the distribution over number of iCAPs versus total time. We then specified which iCAP combinations occur mostly for different number of overlapping iCAPs (Fig. 5 and Supplementary Fig. 9 for weight dependent combinations and histogram at level of overlap). For overlapping iCAPs (i.e., from two to five), the iCAP distribution and the most frequent 20 iCAP combinations are illustrated with their respective activated (red) and deactivated (blue) state. For instance, the most frequent combination at level two is the DMN (positive) and anterior salience network (negative), followed by DMN (negative) and secondary visual (positive).