FINE-SCALE PATTERNS DRIVING DYNAMIC FUNCTIONAL CONNECTIVITY PROVIDE MEANINGFUL BRAIN PARCELLATIONS

Dynamic functional connectivity (dFC) derived from resting-state functional magnetic resonance imaging (fMRI) allows identifying large-scale functional brain networks based on spontaneous activity and their temporal reconﬁgurations. Due to limited memory and computational resources, these pairwise measures are typically computed for a set of brain regions from a pre-deﬁned brain atlas, which choice is non-trivial and might inﬂuence results. Here, we ﬁrst leverage the availability of dynamic information and new computational methods to retrieve dFC at the ﬁnest voxel level in terms of dominant patterns of ﬂuctuations, and, second, we demonstrate that this new representation is informative to derive meaningful brain parcellations that capture both long-range interactions and ﬁne-scale local organization. We analyzed resting-state fMRI of 54 healthy participants from the Human Connectome Project. For each position of a temporal window, we determined eigenvector centrality of the windowed fMRI data at the voxel level. These were then concatenated across time and subjects and clustered into the most representative dominant patterns (RDPs). Each voxel was then labeled according to a binary code indicating positive or negative contribution to each of the RDPs. We obtained a 36-label parcellation displaying long-range interactions with remarkable hemispherical symmetry. By separating contiguous regions, a ﬁner-scale parcellation of 448 areas was also retrieved, showing consistency with known connectivity of cortical/subcortical structures including thalamus. Our contribution bridges the gap between voxel-based approaches and graph theoretical analysis.


INTRODUCTION
Functional connectivity (FC) based on resting-state (RS) functional magnetic resonance imaging (fMRI) allows to investigate networks underlying fluctuations of spontaneous brain activity [1].Typically, measures of FC such as Pearson correlation are computed between different regions of the brain over a whole resting-state scan of several minutes, yielding a scan-averaged, stationary estimate of connectivity.Since the number of pairwise FC measures grows quadratically with the number of spatial locations, it is common to rely upon brain atlases; in fact, FC matrices considering all gray-matter (GM) voxels could end up containing 10 10 entries.In addition, recent evidence has pointed at the intriguing temporal fluctuations of FC and led to the development of dynamic functional connectivity (dFC) methods, where a sliding-window approach is typically adopted to explore networks reconfigurations over time [2].Dynamic analysis further exacerbates dimensionality issues, as 10 2 -10 3 temporal windows need to be considered for each scan.
Whole-brain parcellations have been proposed based on different criteria, such as anatomical [3], functional [4][5][6], or multimodal [7,8].However, the use of different parcellation schemes could lead to substantially different network structure and statistics [6,9], and there is no consensus on the ideal parcellation resolution, which is a compromise between spatial variability of the functional signals and anatomical interpretability [6,10].Despite this, due to computational limitations only few studies have explored FC at the voxel resolution with direct assessment of the FC matrix [11][12][13][14][15][16][17].Others [18,19] circumvented the problem by investigating eigenvector centrality of functional networks, with a fast computation method introduced by [20] that avoids the explicit calculation of the FC matrix.
In this work, we build upon the concept of eigenvector centrality [15,20] to integrate it into a dynamic analysis, and yield a data-driven voxel-wise estimation of dFC.The retrieved dominant patterns of dynamic FC offer a good approximation of the voxel-level connectivity matrix and their changes over time.Then, we investigate to what extent FC fluctuations would be able to drive themselves a meaningful parcellation, and whether they would be informative both on long-range and fine-scale organization.To do so, dFC dominant patterns are concatenated across participants and clustered to obtain the most representative dominant patterns (RDPs).Each voxel is then labeled according to its negative/positive contribution to each RDP, yielding a data-driven parcellation guided by fluctuations of FC.

Preprocessing
Minimally preprocessed HCP images were used as starting point [21].Additionally, spatial smoothing (Gaussian kernel with FWHM=5mm) was performed with SPM8 (FIL,UCL,UK).The first 10 volumes were discarded so that the fMRI signal achieves steady-state magnetization, resulting in T = 1190 time points considered.Voxel timecourses were detrended and nuisance variables were regressed out (6 head motion parameters, average cerebrospinal fluid and white matter signal computed in standard masks mapped to the subject's fMRI space and masked with individual segmentation maps).Then, the preprocessed voxel time courses were band-pass filtered ([0.0167 − 0.15Hz]) and the volumes where masked using a standard GM parcellation in MNI coordinates (IIT GM Destrieux mask, http://www.iit.edu/),previously resliced to fMRI resolution and masked with the population mean functional volume (N v = 109 ′ 783 voxels considered).Finally, the global signal at every time point was removed to discard a global component present in the dFC-driven clusters retrieved in the following.

Voxel-level dFC and Representative Dominant Patterns
Fig. 1 illustrates the full method pipeline.First, we retrieved subject-specific dominant patterns of dFC with the technique detailed hereafter, that we previously tested on a preliminary analysis conducted on a restricted sample [22].For each subject s, s = 1, . . ., N S , X (s) denotes the N V × T fMRI data matrix , containing in its rows the preprocessed timecourses of GM voxels.A sliding window approach with window length N T = 83 TR (equivalent to 59.76 sec) and step ∆ = 2 TR (1.44 sec) yielded N W windows to be analyzed (Fig. 1,  A).After normalization of the timecourses using z-scoring and division by the square root of the window length, the fast eigenvector centrality method -initially proposed by [20] for stationary FC-was applied for each windowed data matrix i (Fig. 1, A), which is optimal in terms of explained variance: i ) T , and then implementing products from right to left: first (X , which is then premultiplied by X (s) .Therefore, the matrix C (s) i is never explicitly computed nor stored, as we express all operations on C (s) i as equivalent operations using the implicit representation provided by u (s) i .Further, dFC values are usually temporally centered for each connection through a row wise demeaning of the connections' timecourses, to focus exclusively on the contribution of dynamics [23].We achieved this centering by subtracting from each windowed connectivity matrix C (s) i the subject-wise stationary connectivity matrix C (s) .As also C (s) is too large to be computed/stored, we first approximated it by its M largest eigenvectors (we choose M = 50): where (s) .The centralities after temporal centering can then be computed considering the difference C (s) i − C (s) in the eigenvalue equation, leading to: which can be incorporated in the eigenvector centrality computation and where operations are implemented from right to left and only involve matrix-vector products.The ARPACK software library that is available in Matlab was used to compute eigenvectors.
The eigenvector centralities u (s) i were then aggregated across windows and subjects into the N V × N W N S matrix U (Fig. 1, B), and k-means clustering with 10-fold crossvalidation based on cosine similarity was applied to retrieve the most representative dominant patterns (RDPs) (Fig. 1, C).

dFC-driven parcellation
A label was assigned to every cortical voxel corresponding to its positive/negative contribution to each of the RDPs; i.e. voxels with the same dFC behaviour were aggregated together (Fig. 1, D).These labels identify dFC patterns reflecting longrange interactions, since they are made out of clusters distributed across the brain.Each of the dFC patterns was then further split into its contiguous regions.Regions with less than 20 voxels (0.16 cm 3 ) were removed.To verify consistency with previous knowledge, a specific analysis of the thalamus was performed.An existent thalamus mask [3] was used to mask the atlas and the atlas regions of interests belonging to this structure were analyzed in terms of subdivisions and connectivity.

RESULTS AND DISCUSSION
The clustering cross-validation identified K = 6 as optimal number of clusters.The resulting six RDPs (Fig. 1,  C), display long-range patterns of connectivity including the main known RS networks (e.g., default mode (DMN), frontoparietal, salience, attention) and were expressed uniformly across subjects.Interestingly, these networks appear sometimes fully, sometimes only partially recruited in the different RDPs; e.g., the DMN is included in its entirety in RDP5, opposed to primary sensory / sensorimotor regions, while it is segregated in its ventral and dorsal parts in RDP4 and 6.This goes well with the view seeing the DMN as a heterogeneous and dynamical network, which can also partially co-occurr with attention network in specific moments [24,25].Further, the selective recruitment of different DMN subportions could relate to the involvement of this network in different types of internal processes during mind-wandering [24].
Voxel labeling based on their contribution to the six RDPs led to a parcellation including 36 unique labels.These identify patterns of long-range interactions, such as the RDPs from which they originate, and characterized by striking inter-hemispherical symmetry.Splitting these 36 patterns into contiguous regions led to a finer-scale parcellation including 448 areas, which can be interpreted as the finest scale at which dFC remained meaningful.Avoiding the need of a-priori choosing regions of interest, these parcellations have the advantage of being whole brain, completely data-driven and taking into account the dynamic nature of FC.So far, only few previous reports consider dFC to determine a parcellation [26,27], but no previous work has performed the dFC analysis in the whole brain at the voxel level.It should be emphasized that the parcellation is obtained from temporal clustering of the dFC dominant patterns, and thus no spatial constraint (e.g., voxel neighborhood) was taken in account.
The thalamus was subdivided by our parcellation into 10 subregions.A noteworthy inter-hemispherical symmetry, despite the fully data-driven procedure, can be observed for this structure as well (Fig. 3) and is encouraging about the parcellation significance.The obtained parcels do not always precisely resemble the known thalamic anatomy, suggesting that the dynamic behavior of the thalamus can lead to finer subdivisions or aggregate known anatomically distinct regions.Only one previous study attempted a parcellation of the thalamus based on dFC [27], by clustering (spatially and temporally) patterns of connectivity between thalamic voxels and five main cortical regions.In that case, however, the segmentation was dependent on the a-priori chosen regions, therefore a comparison is not straightforward.The extracted thalamus dFC network (Fig. 3, A) shows a strong interplay between thalamus and cerebellum (7 of the 10 atlas labels involving the thalamus also include the cerebellum), as expected from previous knowledge describing afferent connections from the cerebellar deep nuclei passing through the thalamus and proceeding to the cortical areas.In addition, the distributed network of cortical regions involved in thalamic connectivity, including basal ganglia (caudate, putamen), frontal (superior medial, middle frontal, prefrontal areas) and parietal regions, hippocampus and amygdala (Fig. 3, B), well fit with the vision of the thalamus as a functional hub integrating multimodal information across diverse cortical functional networks [28].Subcortical structures like basal ganglia [29] were also already known in literature to be strongly interconneted with the thalamus, as well as the hyppocampus and amygdala [30] which were found here.
As methodological considerations, several controversies have been generated in the past by the application of global signal regression (GSR) to resting-state fMRI data.Here, we repeated the analysis without GSR [31], which showed similar results apart from a strong global component emerging as RDP, previously also found [23] and reflecting global fluctuations of dFC.The effect of GSR was then only discarding this component to facilitate the extraction of dynamical networks present on top of it.Finally, possible future improvements could see the use of higher-rank approximations for the connectivity at every timepoint, as well as the use of soft instead of hard clustering to provide a parcellation with overlapping regions.

CONCLUSION
We introduced a novel method considering for the first time voxelwise dominant patterns of dFC.The results show how large-scale functional networks interact dynamically during a resting-state scan.The dominant patterns of dFC can then be used to build a meaningful data-driven whole-brain parcellation.The initial 36 labels reflect long-range patterns with striking hemispherical symmetry, and are reminiscent for known resting-state networks.Further imposing contiguity led to 448 regions at a strikingly fine scale, which highlights the richness of spontaneous brain activity.

Fig. 1 .
Fig. 1.Data processing pipeline.(A) The sliding-window connectivity is approximated at each window by eigenvector centrality.(B) Centralities are concatenated across time and subjects and (C) k-means clustering is applied to retrieve dFC representative dominant patterns (RDPs).(D) Voxels are labeled according to their positive/negative contribution to each RDP, yielding a dFC-driven parcellation.

Fig. 3 .
Fig. 3. Thalamus subdivisions and their connectivity.(A) The regions of the atlas which cover the thalamus are shown for the whole brain.(B) Close up on the four largest thalamic subregions, specifying the cortical/subcortical regions to which they connect.Abbreviations: DMN: default mode network, FPN: fronto-parietal network, MCC: middle cingulate cortex, ACC: anterior cingulate cortex.