Individual Topological Analysis of Synchronization-Based Brain Connectivity

: Functional connectivity analysis aims at assessing the strength of functional coupling between the signal responses in distinct brain areas. Usually, functional magnetic resonance imaging (fMRI) time series connections are estimated through zero-lag correlation metrics that quantify the statistical similarity between pairs of regions or spectral measures that assess synchronization at a frequency band of interest. Here, we explored the application of a new metric to assess the functional synchronization in phase space between fMRI time series in a resting state. We applied a complete topological analysis to the resulting connectivity matrix to uncover both the macro-scale organization of the brain and detect the most important nodes. The synchronization metric is also compared with Pearson’s correlation coefﬁcient and spectral coherence to highlight similarities and differences between the topologies of the three functional networks. We found that the individual topological organization of the resulting synchronization-based connectivity networks shows a ﬁner modular organization than that identiﬁed with the other two metrics and a low overlap with the modular partitions of the other two networks suggesting that the derived topological information is not redundant and could be potentially integrated to provide a multi-scale description of functional connectivity.


Introduction
A wide range of biological phenomena are characterized by synchronization of systems whose complex cooperation and integration explain the emergence of vital functions [1]. Several synchronization metrics have been introduced to quantify the coupling behavior of interacting dynamic systems [2,3]. In particular, there has been a growing body of literature addressing the application of complex network analysis for the characterization of dynamical systems based on time series. These combined approaches have clarified fundamental issues about the organization of nonlinear dynamics in different fields. As an example, the complex and nonlinear nature of brain dynamics is widely discussed in the literature [4][5][6]. The human brain is modeled as a complex network composed of anatomical and functional sub-systems whose interactions allow the performance of both high-and low-level cognitive functions. Such interactions have been extensively modeled through the mathematical framework of the complex networks, which allows the treating of each subsystem as a node of a network and to identify specific topological properties of nodes and links in the network [7,8]. Topological data analysis have been largely applied to investigate the architecture of brain networks and connectivity differences between individuals at several levels of the data hierarchy [9][10][11]. These techniques can be applied both to functional data, which describe the dynamic brain activity through functional magnetic resonance imaging (fMRI) or electrophysiological acquisitions [12,13]. In particular, over the past decade, there has been an increasing interest in inferring connectivity properties from fMRI data [14,15]. Functional connectivity (FC) analysis aims at assessing the strength of functional coupling between the signal responses in distinct brain areas [16]. According to the complex network framework, the anatomical regions of interest are modeled as the nodes of the network, linked by edges resulting from the selected synchronization metric. Pairwise Blood oxygenation level dependent (BOLD) time series connections can be estimated through different metrics including cross-correlation, cross-coherence and causal modeling [17]. Each FC metric quantifies different forms of dependency between the BOLD signals: zero-lag correlation metrics quantify the statistical similarity between pairs of region, other metrics assess synchronization at frequency band of interest [18,19], others capture direct coupling between time series [20,21].
Recently, new metrics to assess the coupling behavior of BOLD time series have been introduced to capture both linear and nonlinear interactions in human brain. In particular, a generalized synchronization index (SYNC) has been proposed to assess the functional connectivity with fMRI data [22,23]. Such metric has been proven effective in detecting active sub-networks during a high-level cognitive process, describing age-related functional patterns [24] and identifying and quantifying pathological states [25].
The analyses of functional brain connectivity in the state of rest have revealed different modular configurations of the brain activity, which reflect specific functions and varied spatial topologies [26]. Moreover, the role of specific regions acting as network hubs for the efficient functioning of the human brain has been extensively documented [27,28]. Brain hubs facilitate the integration of functionally specialized neural systems by supporting long-range connections [29,30] and by routing a large amount of neural information [31]. The dysfunction of some brain hubs has been also related to pathological conditions and psychiatric disorders [32,33].
In this work, we performed a comparison of the synchronization connectivity pattern derived from the SYNC metric, statistical Pearson correlation and spectral coherence functional connectivity. In particular, we exploited the resting state fMRI data of a single subject to construct three functional connectivity networks and explore the topological organization of their macro-scale architectures by means of a modularity analysis. We investigated the topological properties of the most important nodes of the networks to detect stable hubs. Here, the SYNC hub nodes represent regions with a high level of synchronization over the entire time interval with different functional sub-systems. Such nodes differ from the hubs of correlation and coherence networks as statistical correlation is related to the average synchronization level over the observation period, while spectral coherence can identify the frequencies that drive the linear association between the two time courses. One of the principal objectives of this work is to verify whether an overlap of these nodes exists with respect to those identified in the Pearson and coherence networks in order to better elucidate the synchronization index method in functional connectivity analysis. We stress that the main focus of this exploratory analysis is a detailed comparison of the three network topologies concerning their modular organization and the roles of the regions acting as hubs in such modular configurations. Recent precision functional mapping has shown that individual brain organization is qualitatively different from group average estimates and that individuals can exhibit distinct brain network topologies [34,35]. For this reason, we have chosen a single subject-based analysis. Indeed, the individual analysis address a targeted comparison of the regions included in each module and the hub roles between the three functional networks obtained respectively with SYNC metric, Pearson statistical correlation and spectral coherence. We also tested the reliability of the results by analyzing ten scans of the same subjects.
The rest of this paper is structured as follows. In Section 2 a description of the data is provided. In Section 3.1 the analysis of synchronization of dynamic processes in the phase space is described with an illustrative example of the proposed metric. Section 3.2 is devoted to both topological and modularity analysis of fMRI connectivity matrices. In Section 4 the results of the analysis for the SYNC, the Pearson and coherence networks are presented with a discussion on the modular organization and the most important nodes of the three networks and a comparison between them. Finally, the conclusion of the work is presented in Section 5.

fMRI Data
We included ten resting-state baseline fMRI sessions of two subjects (first subject 0025427, sex: male, age: 23 years; second subject 0025432, sex: male, age: 21 years) from the Hangzhou Normal University (HNU) cohort of the Consortium for Reliability and Reproducibility (CoRR) [36]. During functional scanning, subjects of the sample were presented with a fixation cross and were instructed to keep their eyes open, relax and move as little as possible while observing the fixation cross. Subjects were also instructed not to engage in breath counting or meditation. More details about recruitment, inclusion criteria and assessment protocols are available on the project website (http://fcon_1000.projects.nitrc.org/indi/CoRR/html/hnu_1.html).
The data sessions were preprocessed using the configurable pipeline for the analysis of connectomes (CPAC). CPAC preprocessing includes slice time correction, motion correction, skull strip, global mean intensity normalization and nuisance signal regression with 24-motion parameters, signals from white matter (WM) and cerebrospinal fluid (CSF), and linear and quadratic low-frequency drifts. After band-pass filtering (0.01-0.1 Hz), the brain volume of each session was divided in 116 non-overlapping anatomical regions of interest according to the AAL Atlas [37] and the average fMRI time series of each ROI was collected, obtaining 116 time series for each session.

Methods
Functional synchronization connectivity between all pairwise combinations of ROI time series was assessed by estimating their cross-recurrence plots (CRP) and then by calculating the SYNC index as detailed in Section 3.1 resulting in a 116 × 116 weighted connectivity matrix for each fMRI session. We also computed the correlation between each couple of time series by means of the Pearson correlation coefficient and the coherence metric, resulting in a 116 × 116 Pearson connectivity matrix and a 116 × 116 coherence connectivity matrix. We performed a modularity analysis to determine the most representative community structure of each functional network across sections and compare them as explained in Section 3.2.1. Finally, we investigated the topological roles of the most important hubs in each connectivity networks to verify the overlap across the three networks as reported in Section 3.2.2. We repeated the same processing steps and analysis on another subject of the same dataset whose results are included in Supplementary Information file.

Synchronization in Phase Space
A system could be represented in the phase space, i.e., a topological representation of its behavior under different initial conditions. This method assumes that each signal represents a projection of a higher-dimensional dynamical system evolving in time, whose trajectories are embedded into a manifold, i.e., a region of its phase space. The phase space of a system could be reconstructed by means of the Takens's Theorem [38]. For a generic single temporal observation of a system u(t), the trajectory x i is expressed as: where m is the embedding dimension and τ is the time delay. Both parameters must be properly selected to avoid redundancy in the phase space. The dimension m of the reconstructed phase space should be large enough to preserve the properties of the dynamical system: m ≥ 2D + 1, where D is the correlation dimension of the original phase space. The correct time delay τ should be chosen by determining when the samples of the time series are independent enough to be useful as coordinates of the time delayed vectors. For the estimation of the embedded parameters m and τ several techniques have been proposed. As an example, the first local minimum of average mutual information algorithm [39] can be used to select the proper time delay. The minimum embedding dimension is usually estimated through the false nearest neighbors (FNN) algorithm [40].
After the reconstruction of the dynamic trajectories of the two systems in the phase space, it is possible to quantify their interacting behavior by projecting the phase space into the bidimensional cross-recurrence plots (CRPs) [41]: where Θ is the Heaviside function, is a threshold for closeness, N is the number of considered states for each system and || · || the maximum norm function. A CRP represents an effective method to reduce the dimensionality of the phase space and compare the trajectories of the interacting systems as it is a matrix whose entries include information on the degree of closeness of each state of the first system with each state of the second system. A CRP exhibits characteristic patterns that show local time relationships of the segments of the trajectories of the two interacting systems. As an example, diagonal lines occur when the evolution of the states is similar at different times and their lengths are related to the periods during which the two systems move in similar ways remaining close to each other [42]. The main diagonal of a CRP represents a line of synchronization (LOS), since it implies the identity of the states of the two systems in the same time intervals. The LOS structure can be analyzed to extract information about the synchronization of the two time series [43]. In particular, the presence of LOS suggests that the two time series are fully synchronized, while discontinuities appear when the two signals are not fully synchronized. We defined the synchronization time metric (SYNC) to quantify the mean period during which the two systems are synchronized in order to reflect the dynamical synchronization behavior of the series throughout the observation period [22,25]. SYNC is proportional to the ratio of the sum of the lengths of the subsegments l j along the LOS to the total number of samples N: The coupling behavior in phase space can be illustrated by referring to a pair of simple sine waves. We compare the SYNC approach to Pearson correlation coefficient computed as: are the two time series andx andȳ denote their sample means. We also compared the SYNC index to spectral coherence C xy ( f ) across the low-frequency band (0.01-0.1 Hz): where S x ( f ) and S y ( f ) are the power spectrums of the time series estimated using Welch's method [44] and S xy ( f ) is their cross spectrum. In Figure 1, the CRPs and the extracted LOS are represented for three different cases: 1. a sine wave synchronized with itself throughout the observation period 0-T; 2. a sine wave and a second sine wave fully synchronized during half of the observation period 0-T/2; 3.
a sine wave and a second intermittently synchronized sine wave during a first observation interval 0-T/4 and a second observation interval T/2-3T/4.
For the first case, R = 1, C xy = 1 and SYNC = 1, since the LOS is complete. For the second case R = 0.7, C xy = 0.16 and SYNC = 0.5, while for the third case R = 0.7, C xy = 0.15 and SYNC = 0.25. The Person correlation measures the overall statistical similarity between the two time series, the coherence metric measures their linear relationship in the frequency domain while the SYNC metric takes into account also the intermittency.

Identification and Analysis of Functional Modules
Louvain algorithm [45] has been used to find communities of ROIs in the three functional networks obtaining three partitions for each fMRI session. This method aims at maximizing a modularity metric that evaluates the quality of a partition by comparing the density of connections within a community to that expected in a random network. This algorithm optimizes the modularity function defined as: where A i,j is the edge weight between i and j, k i is the sum of the weights of the links attached to node i , c i is the community assigned to the node i , m is the sum of all of the links of the networks and δ(α, β) = 1 if α = β and δ(α, β) = 0 if α = β.
An overview of the analysis is shown in Figure 2. To select the best resolution value, the algorithm was applied for γ = (0, 0.1, ..., 1.7). We performed 1000 iterations of the Louvain algorithm for each value of γ. Finally, we selected the γ value corresponding to the maximum Q value averaged over the 1000 partitions. This procedure returns a robust resolution index value, but also a set of different partitions corresponding to the optimal γ value. To identify the most representative community partition, we applied a consensus algorithm. This approach removes the statistical fluctuations of the modularity algorithm and provides a stable output, as well as identifying only one partition for each network [46]. We compared the community partitions resulting from the three networks within each session by using the normalized mutual information (NMI) [47]. For two networks with partitions, respectively, A and B, it is defined as: where I(A, B) is the mutual information between the two partitions and H(A) and H(B) are the entropy of A and B. This metric ranges between zero (if and are completely independent) and one (if and are identical).

Identification of Hub Nodes
Two graph metrics have been used to characterize the topological organization of the synchronization-based functional network, Pearson correlation network and coherence network. More specifically, the participation coefficient and within-module degree z score were computed to assess the regional importance of the brain nodes. For a community partition, • the participation coefficient PC i is defined as: where K i is the sum of the weighted links attached to i, K is is the sum of weighted links from i to community s and N M is the total number of modules. Thus, the participation coefficient quantifies how the links of a node are distributed across modules. A node's participation coefficient is maximum if it has an equal sum of edge weights to each module in the network. A node's participation coefficient is 0 if all its links are within a single community. • the within-module degree z score for each node is computed as: where k i is the number of links of node i to other nodes in its community s i ,k s i is the average of k over all of the nodes in s i and σ k s i is the standard deviation of k in s i . Thus, the within-module degree z score measures how well-connected node i is to other nodes in the community relative to other nodes in the community.
These graph metrics are often used to capture the centrality of nodes in brain networks and to map regions that act as hubs [48]. Here, we used a cartographic system mapping each node in a plane through the ordered pair of values (PC i , z i ) to identify the connective hubs, i.e., those nodes with higher connection strength and high centrality values [49][50][51]. This cartographic system allows identification through a statistical threshold level, different roles of the nodes in a network. Hub nodes are central nodes, so they are characterized by high z values. In addition, hubs can be divided into three different categories: provincial hubs, i.e., hub nodes with the vast majority of links within their module (lower values of PC i ); connector hubs with many links to most of the other modules and kinless hubs with links homogeneously distributed among all modules (highest PC values). We used the 25th and 75th percentile of both metrics to divide the plane (PC, z) into nine regions to identify the hubs of the networks. To compare the topologies of the three functional networks across all the sections of the subject, we included into the final three sets only the hubs that appear with a specific role in at least six sessions of each functional connectivity network.

Network Structure
The three functional connectivity networks are shown in Figure 3 with the empirical probability distributions of the weights in semi-logarithmic scale. The three matrices are very different form each other: the SYNC network is sparser, with only a small number of entries with high synchronization values. This aspect is also emphasized by the probability distribution which shows a right long-tailed trend, with few outliers between 0.2 and 1. This finding highlights a fist important difference between the connectivity matrices, i.e., most regions exhibit a very intermittent synchronization behavior, although on average the pairs of ROIs are moderately statistically correlated in time and frequency (median value of Pearson matrix ∼0.3 and median value of coherence matrix ∼0.25). The same results were found for all the sessions of the subjects as reported in Figures SI-1 and SI-2. Additionally, the Wilcoxon rank sum test was used to compare the distributions of each couple of metrics, i.e., SYNC-Pearson, SYNC-coherence, Pearson-coherence within each session. For each comparison we found p < 0.0001 confirming that the connectivity matrices are statistically different at 5% significance level.

Modularity Analysis
To explore the topological functional organization of the brain regions, the Louvain algorithm was computed 1000 for each network. Figure 4 shows the error bars of modularity index Q over the 1000 iterations of the algorithm for the optimum resolution value γ = 1. It is worth noting that the modularity index is related to the difference between the within-module interactions and the between-module interactions of a network, so it statistically quantifies the goodness of a network partition into modules [52]. Here, the SYNC network, exhibits a significantly higher value of Q with respect the Pearson and coherence networks for all sessions, pointing out a better hard partition into modules. This result confirms the findings of our previous work in which a modularity analysis was applied to networks defined with both the SYNC index and Pearson's correlation coefficient of task-related fMRI data [22]. In particular, both a greater Q value and a higher number of functional communities activated during the working memory task were better identified in SYNC networks, while in Pearson's correlation networks modularity patterns reflecting the working memory load were not detected during the task.
The NMI values reported in Figure 5 highlight that the three partitions are quite different from each other for all the sessions. However, the NMI values between the Pearson and coherence networks are slightly higher, showing a higher degree of agreement between the partitions of these two networks. Figure 6 shows the most representative network partition across the fMRI session of the subject. Seven modules were detected in the SYNC matrix; the first module comprises areas almost becoming from the cingulate cortex, e.g., cincgulate gyrus, anterior and mid cingulate. Several specialized control areas of frontal regions such as medial orbital superior frontal gyrus are included in the second module. The third module resembles both peripheral and central visual networks as it includes all the occipital lobe, the calcarine fissure and surrounding cortex. The fourth module resulted the largest as it includes all the other regions from the frontoparietal network. The fifth module overlaps with most regions of the cerebellar network; the sixth module includes: superior frontal gyrus, middle temporal gyrus, inferior temporal gyrus and cerebellar crus areas. These modules are equally distributed across the default and control networks previously described in the resting-state literature [53,54]. The seventh is mainly composed by areas from the limbic system such as amygdala and para-hippocampal gyrus.  In contrast, the Pearson and coherence network are composed by only four and three modules, respectively. The first module includes some central structures such as caudate nucleus, putamen and thalamus, the inferior frontal and orbital gyrus and several regions from cerebellar networks. The second module resulted the most specialized as it covers all the occipital lobe and several regions from cerebellar network. Most of the regions from parietal, frontal and temporal lobe are equally distributed across the first and third module. In particular, the third module comprises the cerebellar crus 1 and 2, all the areas of the superior and middle frontal lobe and orbital gyrus. The fourth module resulted the smallest one including most of the regions from cerebellum and vermis.
The first module of the coherence network covers all the parietal lobe and middle frontal areas. The second module includes the areas of the occipital lobe, superior frontal, orbital and temporal regions and the cerebellar crus 1 and 2. Finally, the third module is mainly composed of the cerebellar network, inferior frontal gyrus and sub-cortical regions.
The three network configurations resemble different functional organization of the brain during the resting state. The SYNC network shows very specialized co-synchronized patterns, also including very small modules clearly related to salient control functions of the brain; on the other hand, both Pearson and coherence connectivity capture functional organization at macro-scale highlighting patterns of different degree of average statistical correlation between regions. These findings suggest that the modular organization of the three functional connectivity could be integrated to describe the hierarchical behavior of the multi-scale synchronization between the interacting regions during different cognitive states or conditions.

Hub Identification
A cartographic system with the statistical thresholds of 25th and 75th percentiles for both distributions of PC and z was used to identify the hub regions. As an example, in Figure 7 nine sectors of the resulting partitions are shown for the SYNC network related to the first session of the subject together with the distributions of the two topological metrics. The hub regions are confined in the upper part of the plane. In details, kinless hubs, associated with more modules, are localized in the upper right, peripheral hubs with most of the connections within their own module are localized in the upper left side, while the other connector hubs are identified as between these two extreme categories. In Table 1 is reported the complete list of all the regions identifies as hubs with their categories for the three functional networks. It is interesting to note that beyond the specific role designated to each hub in the three network configurations, in the Pearson's network many connector hubs are located in the temporal lobe, which instead show just two hubs in the SYNC network, while most of the hubs of the coherence are located in the frontal lobe. Several nodes we detected as connector hubs were also identified as integrative nodes in frontal, inferior parietal areas, cingulate and some regions of both default-mode and executive control networks in the adult brain [55][56][57]. In addition, in both SYNC and Pearson networks, the cerebellar crus regions are detected as important hubs, confirming the complex associations among the cerebellar regions and the cerebral cortex captured during the resting state by means of functional connectivity MRI [58].
Another interesting aspect concerns the presence of kinless hubs: such kind of nodes belong to any modules, hence they particularly facilitate the integration of information among the brain areas. These nodes are identified in the upper right side of the plane. However, in many works, advanced methodologies are proposed to identify functional connectivity hubs [59]. In this work, we explored an intuitive approach to depict central and influential nodes in functional connectivity networks. In general, hubs are identified by using consensus scores across several centrality metrics: the nodes with the highest scores on a set centrality measures are usually designated as hubs [60]. However, the level of influence of such nodes depends on the distribution of the metrics used to define hubness. As an example, if the distribution is long-tailed, the network contains a subset of outlier nodes that are considerably more central than others; therefore, they can be considered influential hubs [57]. Here, the range limits for the participation coefficient PC are defined by simple statistical thresholds, hence these ranges correspond to the extreme values assumed by the coefficient in both networks. Actually, in the literature, a strict definition of the kinless hub role corresponds to PC > 0.7. According to many studies, these nodes are almost completely absent in the cortex of primates [60]. As shown in Figure 8, we found lower values of the upper limit of PC in both Pearson and coherence networks in all the sessions, while the SYNC networks exhibit values of PC corresponding to the formal definition, effectively detecting the presence of highly versatile hub nodes. These results are confirmed in the connectivity analysis of another subject reported in the Supplementary Information file.   In summary the three networks showed significant differences among them. First, the modularity analysis showed higher values of modularity index Q in SYNC networks than the other two networks thus indicating a clearer division into communities. Second, the analysis of NMI values highlighted that the three networks showed modular partitions with a low overlap between them (below 50%). The different topological configurations in the three networks were confirmed by the analysis of the hubs, which outlined the presence of different central roles of functional areas. The results suggest that statistical correlation and coherence networks show more evenly distributed synchronization patterns of comparable size in the brain, while SYNC networks exhibit more granular partitions, highlighting more varied synchronization patterns.

Conclusions
In this work, we analyzed the topological organization of SYNC-based, Pearson correlation and spectral coherence networks of a single subject. We showed that even if the metrics are mutually linked, the SYNC metric emphasizes the intermittent behavior of the interacting systems showing a greater matrix sparsity. The resulting network topologies were examined to identify the macro-scale modular organization of the nodes and compared to each other. Our analysis revealed some interesting aspects of the structure of SYNC network during resting rest such as the presence of multi-domains highly synchronized brain regions. In addition, the hub nodes of the networks were detected showing only a partial overlap across the three networks. It is important to note that the three measures provide different descriptions of the interaction behaviors between BOLD signals. Indeed the results suggest that the topological information derived from each matrix is not redundant and could be integrated into a multi-view logic to retrieve different information on the interaction and integration of functional sub-systems by using a multilayer connectivity matrix. Future work will explore more in detail the structure of the SYNC connectivity and different features for several tasks and cognitive states in larger populations in order to further elucidate the association between the proposed metric and different synchronization mechanisms in the brain.