# A Novel Synchronization-Based Approach for Functional Connectivity Analysis

Authors: Angela Lombardi; Sabina Tangaro; Roberto Bellotti; Alessandro Bertolino; Giuseppe Blasi; Giulio Pergola; Paolo Taurisano; Cataldo Guaragnella
Journal: Complexity
 (2017)

---

## Abstract

Complex network analysis has become a gold standard to investigate functional connectivity in the human brain. Popular approaches for quantifying functional coupling between fMRI time series are linear zero-lag correlation methods; however, they might reveal only partial aspects of the functional links between brain areas. In this work, we propose a novel approach for assessing functional coupling between fMRI time series and constructing functional brain networks. A phase space framework is used to map couples of signals exploiting their cross recurrence plots (CRPs) to compare the trajectories of the interacting systems. A synchronization metric is extracted from the CRP to assess the coupling behavior of the time series. Since the functional communities of a healthy population are expected to be highly consistent for the same task, we defined functional networks of task-related fMRI data of a cohort of healthy subjects and applied a modularity algorithm in order to determine the community structures of the networks. The within-group similarity of communities is evaluated to verify whether such new metric is robust enough against noise. The synchronization metric is also compared with Pearson’s correlation coefficient and the detected communities seem to better reflect the functional brain organization during the specific task.


---

## Body

## 1. Introduction

The human brain, as many biological systems, can be seen as a complex network of interacting components whose integration leads to a hierarchical architecture of highly specialized modules [1]. A network formulation simplifies the analysis of a complex system by providing mathematical tools able to capture different aspects of its organization in a compact and straightforward manner. Graph theoretical methods have been extensively applied to many neuroimaging datasets in order to describe the topological properties of both functional and structural networks [2, 3].In particular, over the past few years, there has been an increasing interest in inferring connectivity properties from fMRI data. Functional connectivity analysis aims at assessing the strength of functional coupling between the signal responses in distinct brain areas [4]. According to the complex network framework, the anatomical regions of interest are the nodes of the network, connected by edges resulting from the adopted interregional interaction metrics. Pairwise fMRI time series connections are usually estimated through zero-lag correlation metrics, leading to a weighted network whose links quantify the statistical similarity between pairs of regions. Different preprocessing techniques and strategies are also applied in order to extract only relevant information from the functional network, for example, by considering only a range of weights or by applying several thresholds to filter out weak connections [3]. Functional connectivity studies have revealed interesting insights on normal functional brain organization such as property of small-worldness [5], modularity and presence of hub nodes [6], and the existence of critical alterations of low-frequency neural activity patterns in pathological conditions [7]. Among the proposed strategies, some techniques are more established than others, even if there is still no agreement on which ones are the most effective or appropriate.A number of important questions regarding the identification of networks have to be addressed before considering any analysis technique. Recent studies have demonstrated that different edge definitions could affect the topological properties of brain networks obtaining variable findings [8, 9]. Thus, properties like time resolution of the physiological time series under investigation, the effect of the observational noise, and the presence of nonlinear effects should been taken into account for selecting measures for edge definition. The low temporal resolution of fMRI data limits the number of methods that can be used to assess the statistical interactions between the time series. Linear correlation metrics, including Pearson’s correlation and partial correlation, have been used in simulation environment and resting state studies, showing good performances in estimating functional connections in both cases [9, 10]. On the other hand, nonlinear phenomena in the human brain have been explored at various scales, revealing complex coupling mechanisms in both resting state and task-based neural activity [11, 12]. Most of the functional connectivity studies are focused on configurations of intrinsic connectivity networks (ICNs) and therefore did not assess complex connectivity patterns that can arise in the presence of a cognitive task. Indeed, even if a steady intrinsic network architecture has been found at rest and across a large number of tasks and conditions, task-evoked changes of functional connectivity have been also documented, proving the existence of task-specific network configurations [13]. Exploring topological changes in functional networks when the neural activity is modulated by a cognitive task could improve the understanding of some important mechanisms of human cognition, for example, the dynamic balancing of specialization and integration of brain regions for supporting different cognitive loads [4] and the trade-off between connection cost and topological efficiency in information processing [14]. Assessing functional interactions during external tasks should require metrics that (i) are sensitive to nonlinear coupling between time series and (ii) are more robust with respect to noise.In this work, we propose a novel approach for quantifying functional coupling between fMRI time series and constructing functional brain networks. We use a phase space framework to map pairs of signals in their reconstructed phase space, that is, a topological representation of their behavior under all possible initial conditions [15]. 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, that is, a region of its phase space. Cross recurrence plots (CRPs) [16] are then employed to reduce the dimensionality of the phase space and compare the trajectories of the interacting systems. A synchronization metric is finally extracted from the CRP to assess the coupling behavior of the time series.The proposed metric and Pearson’s correlation coefficient are applied to the fMRI data of a cohort of healthy subjects acquired during performing a working memory task to construct weighted networks.At macroscopic level, functional related brain regions exhibit similar BOLD responses. These groups of regions form dense communities that reflect the functional organization of the brain and whose properties can be linked to the topological features of the task-evoked network configuration [17, 18]. The analysis carried out in this work aims at investigating some properties of the modular structure of task-evoked functional networks obtained with Pearson’s correlation metric and the proposed synchronization index in order to understand which index can better highlight the functional organization of distinct subsystems involved in the specific working memory task. Therefore, a modularity algorithm is used to determine the community structure of each functional network. The within-group similarity of communities is evaluated and exploited to verify whether the metrics are sufficiently robust against noise and effective in revealing correlation even in presence of external stimuli. The rationale underlying this choice is that community structure of a group of healthy subject is expected to be highly consistent in presence of the same task.


## 2. Materials

### 2.1. Subjects

We studied50 healthy subjects (age: mean=25, standard deviation SD=6; 24 females) in the analysis. All of them were evaluated using the Non-Patient Structured Clinical Interview for DSM-IV [19] to exclude any psychiatric condition. Other exclusion criteria were a significant history of drug or alcohol abuse, active drug abuse in the previous year, experience of a head trauma with loss of consciousness, and any other significant medical condition. Socioeconomic status (Hollingshead Four Factor Index, [20]), handedness (Edinburgh Inventory) [21], and total IQ (WAIS-R [22]) were also measured (see Table 1). The present study was approved by the local ethics committee (Comitato Etico Locale Indipendente Azienda Ospedaliera Ospedale Policlinico Consorziale Bari). Written informed consent was obtained by all participants after a complete description of the procedures, in accordance with the Helsinki Declaration.Table 1
Demographic data of the healthy cohort (mean ± standard deviation).
Demographic Data
Age (years)
25
±
6
Gender (M/F)
26
/
24
Handedness
0.60
±
0.55
Socioeconomic status
39
±
17
IQ
114
±
4


### 2.2. fMRI Task

Participants performed the N-Back working memory task, in which a sequence of stimuli is presented and the subject has to remember the stimulus from “N” steps earlier. The stimuli consisted of numbers (1–4) presented in random sequence and displayed at the points of a diamond-shaped box. The control condition (0-back) simply required the subjects to identify the current stimulus. In the working memory condition, the task required the collection of a stimulus seen two stimuli earlier (2-back). The task was organized in a block design, consisting of eight alternating 0-back and 2-back conditions, each lasting 30 seconds. Each 30 sec. block includes 14n-back trials with an interstimuli interval of 2000 ms. Each run lasted 4 minutes and 8 seconds, from which dummy scans were acquired and discarded, obtaining 120 volumes.


### 2.3. fMRI Data

Echo planar imaging blood oxygenation level dependent fMRI data were acquired on a GE Signa 3T scanner (GE Healthcare) equipped with a standard quadrature head coil. A gradient-echo planar imaging sequence (repetition time,2000 ms; echo time, 30 ms; thickness, 4 mm; gap, 1 mm; flip angle, 90°; field of view, 24 cm; and matrix, 64×64) was used to acquire images while the subjects performed the tasks.Images were preprocessed using Statistical Parametric Mapping 8 software (SPM8;http://www.fil.ion.ucl.ac.uk/spm). Images were realigned to the first volume in the time series to correct for head motion (<2 mm translation, <1° rotation), resampled to a 2 mm isotropic voxel size, spatially normalized into a standard stereotactic space (Montreal Neurological Institute template) using a 12-parameter nonlinear warping, and smoothed to minimize noise and residual differences in gyral anatomy with a Gaussian filter, set at 6 mm full-width at half-maximum.


## 2.1. Subjects

We studied50 healthy subjects (age: mean=25, standard deviation SD=6; 24 females) in the analysis. All of them were evaluated using the Non-Patient Structured Clinical Interview for DSM-IV [19] to exclude any psychiatric condition. Other exclusion criteria were a significant history of drug or alcohol abuse, active drug abuse in the previous year, experience of a head trauma with loss of consciousness, and any other significant medical condition. Socioeconomic status (Hollingshead Four Factor Index, [20]), handedness (Edinburgh Inventory) [21], and total IQ (WAIS-R [22]) were also measured (see Table 1). The present study was approved by the local ethics committee (Comitato Etico Locale Indipendente Azienda Ospedaliera Ospedale Policlinico Consorziale Bari). Written informed consent was obtained by all participants after a complete description of the procedures, in accordance with the Helsinki Declaration.Table 1
Demographic data of the healthy cohort (mean ± standard deviation).
Demographic Data
Age (years)
25
±
6
Gender (M/F)
26
/
24
Handedness
0.60
±
0.55
Socioeconomic status
39
±
17
IQ
114
±
4


## 2.2. fMRI Task

Participants performed the N-Back working memory task, in which a sequence of stimuli is presented and the subject has to remember the stimulus from “N” steps earlier. The stimuli consisted of numbers (1–4) presented in random sequence and displayed at the points of a diamond-shaped box. The control condition (0-back) simply required the subjects to identify the current stimulus. In the working memory condition, the task required the collection of a stimulus seen two stimuli earlier (2-back). The task was organized in a block design, consisting of eight alternating 0-back and 2-back conditions, each lasting 30 seconds. Each 30 sec. block includes 14n-back trials with an interstimuli interval of 2000 ms. Each run lasted 4 minutes and 8 seconds, from which dummy scans were acquired and discarded, obtaining 120 volumes.


## 2.3. fMRI Data

Echo planar imaging blood oxygenation level dependent fMRI data were acquired on a GE Signa 3T scanner (GE Healthcare) equipped with a standard quadrature head coil. A gradient-echo planar imaging sequence (repetition time,2000 ms; echo time, 30 ms; thickness, 4 mm; gap, 1 mm; flip angle, 90°; field of view, 24 cm; and matrix, 64×64) was used to acquire images while the subjects performed the tasks.Images were preprocessed using Statistical Parametric Mapping 8 software (SPM8;http://www.fil.ion.ucl.ac.uk/spm). Images were realigned to the first volume in the time series to correct for head motion (<2 mm translation, <1° rotation), resampled to a 2 mm isotropic voxel size, spatially normalized into a standard stereotactic space (Montreal Neurological Institute template) using a 12-parameter nonlinear warping, and smoothed to minimize noise and residual differences in gyral anatomy with a Gaussian filter, set at 6 mm full-width at half-maximum.


## 3. Methods

### 3.1. Network Construction

The brain volume of each subject was divided into 246 nonoverlapping anatomical regions of interest (ROIs) according to the Brainnetome Atlas [23]. Thirty regions from the most ventral part of the brain not acquired during scans were discarded and are not included in the following analysis. For each of the 216 remaining ROIs, a single time series was extracted by averaging the fMRI time series over all the voxels within the ROI. The time series were high-pass filtered (cutoff frequency 1/128 s). For each subject, functional connectivity between all pairwise combinations of ROI time series was assessed:(i)
by calculating their Pearson’s correlation coefficient;(ii)
by computing their CRP and then by calculating their synchronization (SYNC) index as described in the following subsection.Finally, for each subject, we identified two undirected weighted networks, whose edges resulted from(1)
the signed pairwise Pearson’s correlation coefficients;(2)
the SYNC indexes.


### 3.2. Synchronization Index

A state of a system is defined by the values of the variables that describe it at a given time. When such system evolves in time, the sequence of all its states forms a trajectory in the phase space, that is, a multidimensional space whose dimension depends on the number of the variables of the system. Starting from different initial conditions, a real physical dissipative system tends to evolve in similar ways, such that its trajectories converge in a region of the phase space called attractor which represents the steady-state behavior of the system [15].In experimental contexts, where the time seriesuii=1N obtained from the sampling of a single observable variable is available, it is possible to reconstruct the phase space of the system under investigation by means of Takens’s Theorem [24]. Accordingly, a state in the reconstructed phase space is given by a m-dimensional time delay embedded vector obtained from time delayed versions of the output signals as(1)x→i=ui,ui+τ,…,ui+m-1τ,where m is the embedding dimension and τ is the time delay.Both parameters have to be properly selected to avoid redundancy in the phase space. The dimensionm 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 [25] can be used to select the proper time delay. The minimum embedding dimension is usually estimated through the false nearest-neighbors (FNN) algorithm [26].The trajectories of two distinct systems with the same embedded parameters can be compared in a CRP [16], 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. In detail, for two systems with trajectories, respectively, x→ii=1,…,N and y→j(j=1,…,N), the CRP is defined as(2)CRi,jϵ=Θϵ-x→i-y→j,where Θ is the Heaviside function, ϵ is a threshold for closeness, N is the number of considered states for each system, and · is a norm function. A generic entry CRi,j in the resulting N×N array is set to one if the distance between the points x→i and y→j is smaller than the threshold ϵ or to zero elsewhere.The value of the parameterϵ must be estimated carefully, as it influences the creation of structures in the plot. The selection of an appropriate value for the threshold ϵ can be made by taking into account the influence of the observational noise that could affect the experimental measures and the minimum distance between the trajectories of the two systems. In general, choosing ϵ equal to few percent of the maximum phase space diameter could ensure a sufficient number of structures in the cross recurrence plot [27], while the appearance of artifacts could be avoided by considering the signal-to-noise ratio for the underlying physical systems [28].A CRP exhibits characteristic patterns that show local time relationships of the segments of the trajectories of the two interacting systems. Typical structures include single dots, diagonal lines, and vertical and horizontal lines. 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 [29]. A CRP can also exhibit the main diagonal known as line of synchronization (LOS). The presence of LOS implies the identity of the states of the two systems in the same time intervals, that is, the (i,i) states, so its structure can be analyzed to extract information about the synchronization of the two time series [30]. In particular, the presence of LOS suggests that the two time series are fully synchronized, while discontinuities appear when the two signals do not have the same frequency and the same phase. Hence, thesynchronization time (SYNC) has been defined as a metric 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.SYNC is proportional to the ratio of the sum of the lengths of the subsegments lj along the LOS to the total number of samples N:(3)SYNC=1Nd∑j=1NdljN,where Nd is the total number of subsegments.For a visual reference, see Figure1. In Figure 1(a), two fMRI unsynchronized time series are compared and in Figure 1(b) are shown two fully synchronized fMRI time series. It is worth noting that in the first case there are discontinuities of LOS, while in the second case the LOS is continuous. Their SYNC values are, respectively, 0.05 and 1.Figure 1
Pairs of fMRI time series and their CRPs for (a) occipital inferior L and frontal medial orbital L (SYNC = 0.05); (b) occipital superior L and occipital superior R (SYNC = 1).
(a)
(b)


### 3.3. Modularity Detection

Several community detection methods have been proposed to find an optimum partition of the nodes into nonoverlapped communities, that is, clusters of nodes that are more densely connected to each other than to other nodes in the network [31–33]. All these methods aim 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. Here, the Louvain algorithm [33] has been used to find communities of ROIs in the two functional networks obtaining two partitions for each subject. The Louvain method is divided into two phases that are repeated iteratively. The first step favors local optimizations of modularity, while during the second step the communities found in the first step define a new coarse-grained network to be evaluated. This algorithm was chosen because it is fast and seems to be less affected by the resolution limit problem (i.e., the capability to detect modules smaller than a certain size) thanks to its multilevel nature. This method optimizes the modularity function defined as(4)Q=12m∑i,jAi,j-kikj2mδci,cj,where Ai,j is the link between nodes i and j, ki is the sum of the weights of the links attached to node i, ci is the community assigned to the node i, m is the sum of all of the links of the networks, and δ is the δ function.


### 3.4. Statistical Analysis of Modularity

A statistical framework was adopted in order to compare the partitions of all the subjects for each functional network [34].The normalized mutual information (NMI) [35] was used to assess the similarity between a couple of community partitions. For two networks with partitions, respectively, A and B, it is defined as(5)NMIA,B=2IA,BHA+HB,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 A and B are completely independent) and one (if A and B are identical).The statistical relevance of the within-group community structure similarity was evaluated through a permutation test. First, a randomly rewired version of each functional network was generated preserving weights, density, and degree sequence, resulting in two groups of networks: the actual and its randomized matching network. Then, the NMI was calculated between all the possible pairs of network partitions within each group. A null distribution was generated by randomizing group labels10000 times and by calculating the permuted within-group mean NMI at each permutation. Finally, a p value was assigned as the number of times that the permuted within-group mean-similarity was greater than the actual within-group similarity, divided by the number of permutations.In order to inspect the consistency of node assignments to specific functional communities, we carried out further analyses on the networks. Since the labels of modules are arbitrarily assigned by the community detection algorithm at each iteration, it is necessary to match the partition values across the subjects for visualizing the group level community structure. This problem can be overcome by finding a template partition as a reference and by reassigning the labels of communities to match the template, while preserving the distinctions between different modules in each partition [34]. In this work, the partitions of each network for both metrics were matched to the most representative network partition of the group, that is, the median determined by pairwise NMI. Once the labels of partitions are reassigned, it is possible to assess the within-group consistency of each ROI in community membership by counting the number of occurrences with which a ROI appears with a particular label.


## 3.1. Network Construction

The brain volume of each subject was divided into 246 nonoverlapping anatomical regions of interest (ROIs) according to the Brainnetome Atlas [23]. Thirty regions from the most ventral part of the brain not acquired during scans were discarded and are not included in the following analysis. For each of the 216 remaining ROIs, a single time series was extracted by averaging the fMRI time series over all the voxels within the ROI. The time series were high-pass filtered (cutoff frequency 1/128 s). For each subject, functional connectivity between all pairwise combinations of ROI time series was assessed:(i)
by calculating their Pearson’s correlation coefficient;(ii)
by computing their CRP and then by calculating their synchronization (SYNC) index as described in the following subsection.Finally, for each subject, we identified two undirected weighted networks, whose edges resulted from(1)
the signed pairwise Pearson’s correlation coefficients;(2)
the SYNC indexes.


## 3.2. Synchronization Index

A state of a system is defined by the values of the variables that describe it at a given time. When such system evolves in time, the sequence of all its states forms a trajectory in the phase space, that is, a multidimensional space whose dimension depends on the number of the variables of the system. Starting from different initial conditions, a real physical dissipative system tends to evolve in similar ways, such that its trajectories converge in a region of the phase space called attractor which represents the steady-state behavior of the system [15].In experimental contexts, where the time seriesuii=1N obtained from the sampling of a single observable variable is available, it is possible to reconstruct the phase space of the system under investigation by means of Takens’s Theorem [24]. Accordingly, a state in the reconstructed phase space is given by a m-dimensional time delay embedded vector obtained from time delayed versions of the output signals as(1)x→i=ui,ui+τ,…,ui+m-1τ,where m is the embedding dimension and τ is the time delay.Both parameters have to be properly selected to avoid redundancy in the phase space. The dimensionm 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 [25] can be used to select the proper time delay. The minimum embedding dimension is usually estimated through the false nearest-neighbors (FNN) algorithm [26].The trajectories of two distinct systems with the same embedded parameters can be compared in a CRP [16], 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. In detail, for two systems with trajectories, respectively, x→ii=1,…,N and y→j(j=1,…,N), the CRP is defined as(2)CRi,jϵ=Θϵ-x→i-y→j,where Θ is the Heaviside function, ϵ is a threshold for closeness, N is the number of considered states for each system, and · is a norm function. A generic entry CRi,j in the resulting N×N array is set to one if the distance between the points x→i and y→j is smaller than the threshold ϵ or to zero elsewhere.The value of the parameterϵ must be estimated carefully, as it influences the creation of structures in the plot. The selection of an appropriate value for the threshold ϵ can be made by taking into account the influence of the observational noise that could affect the experimental measures and the minimum distance between the trajectories of the two systems. In general, choosing ϵ equal to few percent of the maximum phase space diameter could ensure a sufficient number of structures in the cross recurrence plot [27], while the appearance of artifacts could be avoided by considering the signal-to-noise ratio for the underlying physical systems [28].A CRP exhibits characteristic patterns that show local time relationships of the segments of the trajectories of the two interacting systems. Typical structures include single dots, diagonal lines, and vertical and horizontal lines. 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 [29]. A CRP can also exhibit the main diagonal known as line of synchronization (LOS). The presence of LOS implies the identity of the states of the two systems in the same time intervals, that is, the (i,i) states, so its structure can be analyzed to extract information about the synchronization of the two time series [30]. In particular, the presence of LOS suggests that the two time series are fully synchronized, while discontinuities appear when the two signals do not have the same frequency and the same phase. Hence, thesynchronization time (SYNC) has been defined as a metric 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.SYNC is proportional to the ratio of the sum of the lengths of the subsegments lj along the LOS to the total number of samples N:(3)SYNC=1Nd∑j=1NdljN,where Nd is the total number of subsegments.For a visual reference, see Figure1. In Figure 1(a), two fMRI unsynchronized time series are compared and in Figure 1(b) are shown two fully synchronized fMRI time series. It is worth noting that in the first case there are discontinuities of LOS, while in the second case the LOS is continuous. Their SYNC values are, respectively, 0.05 and 1.Figure 1
Pairs of fMRI time series and their CRPs for (a) occipital inferior L and frontal medial orbital L (SYNC = 0.05); (b) occipital superior L and occipital superior R (SYNC = 1).
(a)
(b)


## 3.3. Modularity Detection

Several community detection methods have been proposed to find an optimum partition of the nodes into nonoverlapped communities, that is, clusters of nodes that are more densely connected to each other than to other nodes in the network [31–33]. All these methods aim 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. Here, the Louvain algorithm [33] has been used to find communities of ROIs in the two functional networks obtaining two partitions for each subject. The Louvain method is divided into two phases that are repeated iteratively. The first step favors local optimizations of modularity, while during the second step the communities found in the first step define a new coarse-grained network to be evaluated. This algorithm was chosen because it is fast and seems to be less affected by the resolution limit problem (i.e., the capability to detect modules smaller than a certain size) thanks to its multilevel nature. This method optimizes the modularity function defined as(4)Q=12m∑i,jAi,j-kikj2mδci,cj,where Ai,j is the link between nodes i and j, ki is the sum of the weights of the links attached to node i, ci is the community assigned to the node i, m is the sum of all of the links of the networks, and δ is the δ function.


## 3.4. Statistical Analysis of Modularity

A statistical framework was adopted in order to compare the partitions of all the subjects for each functional network [34].The normalized mutual information (NMI) [35] was used to assess the similarity between a couple of community partitions. For two networks with partitions, respectively, A and B, it is defined as(5)NMIA,B=2IA,BHA+HB,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 A and B are completely independent) and one (if A and B are identical).The statistical relevance of the within-group community structure similarity was evaluated through a permutation test. First, a randomly rewired version of each functional network was generated preserving weights, density, and degree sequence, resulting in two groups of networks: the actual and its randomized matching network. Then, the NMI was calculated between all the possible pairs of network partitions within each group. A null distribution was generated by randomizing group labels10000 times and by calculating the permuted within-group mean NMI at each permutation. Finally, a p value was assigned as the number of times that the permuted within-group mean-similarity was greater than the actual within-group similarity, divided by the number of permutations.In order to inspect the consistency of node assignments to specific functional communities, we carried out further analyses on the networks. Since the labels of modules are arbitrarily assigned by the community detection algorithm at each iteration, it is necessary to match the partition values across the subjects for visualizing the group level community structure. This problem can be overcome by finding a template partition as a reference and by reassigning the labels of communities to match the template, while preserving the distinctions between different modules in each partition [34]. In this work, the partitions of each network for both metrics were matched to the most representative network partition of the group, that is, the median determined by pairwise NMI. Once the labels of partitions are reassigned, it is possible to assess the within-group consistency of each ROI in community membership by counting the number of occurrences with which a ROI appears with a particular label.


## 4. Results

### 4.1. CRP Parameters

We randomly selected a subset of5000 BOLD time series from the whole dataset and applied the FNN algorithm for estimating the embedded dimension and the first local minimum of the averaged mutual information for selecting the proper time delay. We obtainedm = 5.2 ± 0.75 and τ=1.4±0.66, so the embedded parameters were set to m=6 and τ=1. Following the criteria reported in [27, 28], we identified the range [1.2–1.8] for the threshold ϵ. The analysis was carried out with the average value of the range, setting ϵ=1.5.


### 4.2. Statistical Analysis of Modularity

Permutation tests reveal significant differences of modularity structures between all the functional networks and their randomly rewired versions (p=0 for both the couples), indicating different modular decompositions compared to the null models. However, as shown in Figure 2(a), Pearson’s networks exhibit within-group NMI values much lower than those obtained by means of the SYNC metric (see Table 2 for mean, median, and interquartile range quantities). The nonparametric Wilcoxon rank sum test confirmed significant differences between NMI values of SYNC and those of Pearson’s metric (p<0.0001, α=0.05). The ranges of NMI values of SYNC networks are also comparable to those found among control healthy subjects in resting studies at different threshold values of network density [34]. These results suggest that the functional networks constructed with the SYNC metric share more modularity structures than Pearson’s networks and exhibit also a higher signal-to-noise ratio.Table 2
Mean and median (interquartile range) quantities of NMI andQ distributions for the metrics synchronization and Pearson’s correlation.
Distributions
NMI
Q
Synchronization
0.24
;
0.23 (0.11)
0.17
;
0.17 (0.04)
Pearson
0.15
;
0.14 (0.15)
0.11
;
0.11 (0.07)Figure 2
Violin plots of (a) within-group NMI distributions and (b)Q distributions for synchronization and Pearson’s correlation connectivity metrics with their randomized matching networks (in gray).
(a)
(b)In addition, we evaluated the modularity indexQ. This index ranges between 0 and 1 and measures the density of links inside communities as compared to links between communities. As shown in Figure 2(b), the two distributions are significantly different from their random versions (permutation tests: p=0 for both pairings) and the modularity index of the networks obtained with the SYNC metric is higher than that of Pearson’s networks (p value resulting from the nonparametric Wilcoxon rank sum test p=3.52·10-8, α=0.05).


### 4.3. Comparison of Modular Partitions in SYNC and Pearson’s Networks

Since we computed connectivity measures on a time series derived from a working memory task, we expected to find modules related to working memory performance involving the frontoparietal network [36] to motor activity related to the 0-back task [37] and to the default mode network, which is deactivated when performing the task [38]. Figure 3 shows the five modules detected by the Louvain algorithm at group level. The first module includes areas critical for visuospatial memory and closely resembles the classical frontoparietal network. In contrast, the second module includes more medial regions, with nodes belonging to both the anterior and the posterior default mode networks [39, 40]. The third module overlaps widely with the sensory-motor network, including pre- and postcentral nodes, but also areas of the temporal lobe involved in auditory perception. Interestingly, the fourth and fifth module map almost exclusively to subcortical regions, including the dorsal basal ganglia and the thalamus with the ventral striatum, respectively. These regions are involved in working memory performance [41, 42], but it is intriguing to notice that the technique here employed parsed the connectivity of cortical and subcortical regions based on the time series of activation, yielding anatomic information just based on functional activity patterns. Figure 4 shows the two communities identified at group level for Pearson’s networks. The first module comprises most of the ROIs mapped in the first community of the SYNC networks, while the rest of the ROIs are included in the second module.Figure 3
The five group level functional communities detected in SYNC networks. In each row, a single community is shown in four brain views (left side, right side, top side, and bottom side).Figure 4
The two group level functional communities detected in Pearson’s networks. In each row, a single community is shown in four brain views (left side, right side, top side, and bottom side).The consistency of the assignment of brain regions to functional modules for the SYNC networks is shown in Figure5. As can be seen, all the ROIs within the frontoparietal network are the most consistent among the subjects; in contrast, some nodes from the medial temporal lobes, insular gyrus, and globus pallidus are assigned less uniformly to the same community across the subjects. These findings are in line with the crucial involvement of the frontal parietal network in working memory processing [43]. As this map resembles closely an activity group map, these findings highlight that the connectivity assessment we developed is sensitive to the functional role of the modules identified. Overall, the network parsing obtained by the novel technique we reported reveals a pattern of coupling between brain regions consistent with known models of activation and deactivation during task performance. In Figure 6 is shown the within-group consistency of each ROI in community membership for Pearson’s networks. Although the overall consistency seems generally higher due to the lower number of communities (two versus five), a direct comparison with the SYNC networks is possible only for the first module. The one-sided hypothesis Wilcoxon rank sum test confirmed a greater consistency of the ROIs within the first module for the SYNC matrices (median values of consistency: msync=85.5, mpearson=82, p=0.0096, α=0.0) proving a better identification of the frontoparietal network across the subjects in such matrices.Figure 5
Consistency of the assignment of brain regions to modules measured as the frequency of occurrence of the node with a specific label (in percent) for SYNC networks.Figure 6
Consistency of the assignment of brain regions to modules measured as the frequency of occurrence of the node with a specific label (in percent) for Pearson’s networks.


## 4.1. CRP Parameters

We randomly selected a subset of5000 BOLD time series from the whole dataset and applied the FNN algorithm for estimating the embedded dimension and the first local minimum of the averaged mutual information for selecting the proper time delay. We obtainedm = 5.2 ± 0.75 and τ=1.4±0.66, so the embedded parameters were set to m=6 and τ=1. Following the criteria reported in [27, 28], we identified the range [1.2–1.8] for the threshold ϵ. The analysis was carried out with the average value of the range, setting ϵ=1.5.


## 4.2. Statistical Analysis of Modularity

Permutation tests reveal significant differences of modularity structures between all the functional networks and their randomly rewired versions (p=0 for both the couples), indicating different modular decompositions compared to the null models. However, as shown in Figure 2(a), Pearson’s networks exhibit within-group NMI values much lower than those obtained by means of the SYNC metric (see Table 2 for mean, median, and interquartile range quantities). The nonparametric Wilcoxon rank sum test confirmed significant differences between NMI values of SYNC and those of Pearson’s metric (p<0.0001, α=0.05). The ranges of NMI values of SYNC networks are also comparable to those found among control healthy subjects in resting studies at different threshold values of network density [34]. These results suggest that the functional networks constructed with the SYNC metric share more modularity structures than Pearson’s networks and exhibit also a higher signal-to-noise ratio.Table 2
Mean and median (interquartile range) quantities of NMI andQ distributions for the metrics synchronization and Pearson’s correlation.
Distributions
NMI
Q
Synchronization
0.24
;
0.23 (0.11)
0.17
;
0.17 (0.04)
Pearson
0.15
;
0.14 (0.15)
0.11
;
0.11 (0.07)Figure 2
Violin plots of (a) within-group NMI distributions and (b)Q distributions for synchronization and Pearson’s correlation connectivity metrics with their randomized matching networks (in gray).
(a)
(b)In addition, we evaluated the modularity indexQ. This index ranges between 0 and 1 and measures the density of links inside communities as compared to links between communities. As shown in Figure 2(b), the two distributions are significantly different from their random versions (permutation tests: p=0 for both pairings) and the modularity index of the networks obtained with the SYNC metric is higher than that of Pearson’s networks (p value resulting from the nonparametric Wilcoxon rank sum test p=3.52·10-8, α=0.05).


## 4.3. Comparison of Modular Partitions in SYNC and Pearson’s Networks

Since we computed connectivity measures on a time series derived from a working memory task, we expected to find modules related to working memory performance involving the frontoparietal network [36] to motor activity related to the 0-back task [37] and to the default mode network, which is deactivated when performing the task [38]. Figure 3 shows the five modules detected by the Louvain algorithm at group level. The first module includes areas critical for visuospatial memory and closely resembles the classical frontoparietal network. In contrast, the second module includes more medial regions, with nodes belonging to both the anterior and the posterior default mode networks [39, 40]. The third module overlaps widely with the sensory-motor network, including pre- and postcentral nodes, but also areas of the temporal lobe involved in auditory perception. Interestingly, the fourth and fifth module map almost exclusively to subcortical regions, including the dorsal basal ganglia and the thalamus with the ventral striatum, respectively. These regions are involved in working memory performance [41, 42], but it is intriguing to notice that the technique here employed parsed the connectivity of cortical and subcortical regions based on the time series of activation, yielding anatomic information just based on functional activity patterns. Figure 4 shows the two communities identified at group level for Pearson’s networks. The first module comprises most of the ROIs mapped in the first community of the SYNC networks, while the rest of the ROIs are included in the second module.Figure 3
The five group level functional communities detected in SYNC networks. In each row, a single community is shown in four brain views (left side, right side, top side, and bottom side).Figure 4
The two group level functional communities detected in Pearson’s networks. In each row, a single community is shown in four brain views (left side, right side, top side, and bottom side).The consistency of the assignment of brain regions to functional modules for the SYNC networks is shown in Figure5. As can be seen, all the ROIs within the frontoparietal network are the most consistent among the subjects; in contrast, some nodes from the medial temporal lobes, insular gyrus, and globus pallidus are assigned less uniformly to the same community across the subjects. These findings are in line with the crucial involvement of the frontal parietal network in working memory processing [43]. As this map resembles closely an activity group map, these findings highlight that the connectivity assessment we developed is sensitive to the functional role of the modules identified. Overall, the network parsing obtained by the novel technique we reported reveals a pattern of coupling between brain regions consistent with known models of activation and deactivation during task performance. In Figure 6 is shown the within-group consistency of each ROI in community membership for Pearson’s networks. Although the overall consistency seems generally higher due to the lower number of communities (two versus five), a direct comparison with the SYNC networks is possible only for the first module. The one-sided hypothesis Wilcoxon rank sum test confirmed a greater consistency of the ROIs within the first module for the SYNC matrices (median values of consistency: msync=85.5, mpearson=82, p=0.0096, α=0.0) proving a better identification of the frontoparietal network across the subjects in such matrices.Figure 5
Consistency of the assignment of brain regions to modules measured as the frequency of occurrence of the node with a specific label (in percent) for SYNC networks.Figure 6
Consistency of the assignment of brain regions to modules measured as the frequency of occurrence of the node with a specific label (in percent) for Pearson’s networks.


## 5. Discussion

In the current study, a modularity analysis is applied to networks defined with both the proposed SYNC index and Pearson’s correlation coefficient in order to investigate the task-related functional organization of the brain. Modularity is implicitly related to significant self-regulating mechanisms of the human brain: efficient dense within-module processing and sparse fast integration among subsystems reduce noise propagation and latency [44]. Thus, this feature is strictly connected to critical functional organization between brain systems that are specialized to carry out different tasks: modularity is expected to be greater for optimal system organizations, while decreased modularity implies that there are less intramodular edges than intermodular edges [45]. A low level of modularity would not be compatible with a fast adaptation of the human brain in response to external stimuli. Indeed, lack of highly specialized modules may not allow a rapid execution of complex cognitive task [45, 46]. Consistently, a decreased modularity has been associated with brain disorders characterized by abnormal cognitive processing and has been found as a marker of abnormal brain network development [47–49]. Moreover, there is evidence that while the adaptation speed of the functional organization of the brain is not critical among healthy individuals that perform a specific task, modularity is stable across time, suggesting the existence of latent specific task-related modular configurations [17, 50]. The statistical analysis of modularity reveals that a greater structure homogeneity and a higher number of functional communities activated during the working memory task seem to be better identified in SYNC networks, while Pearson’s correlation does not reflect such features expected in a healthy population. In detail, the SYNC networks showed both higher NMI and Q values thus indicating that the extracted modular partitions are more similar to each other across the population and exhibit a clearer division into communities. Indeed the modularity index Q statistically quantifies the goodness of a hard partition as its value is related to the difference between the within-module interactions and the between-module interactions [51]. Furthermore, the consistency analysis in which the partition of each subject is compared with the median partition of the population points out two results: (i) both networks show at the group level a similar first community that resembles the frontoparietal network, but in the SYNC networks other modules that map to systems engaged during working memory performance are detected; (ii) the statistical comparison of the ROIs within the first module highlights a greater consistency of such task-related regions in SYNC networks. These findings suggest that a problem of community resolution is evident in Pearson’s networks, whereas all the regions not included in the first module are identified in a single community without distinction among sensory-motor network, default mode network, and subcortical areas and even the frontoparietal network is identified more weakly across the population.In our framework, the same community detection algorithm was applied to both kinds of networks. Since the algorithm generates a node partition of a connectivity matrix, some properties of the index used to identify the network such as sensitivity to noise and to complex interaction mechanisms occurring among the brain regions could affect the degree of partition of the network into communities. Several brain connectivity metrics have been proposed as alternatives to Pearson’s correlation coefficient. Coherence and partial coherence analysis were applied to fMRI data to extend linear metrics of zero-lag correlation. These spectral measures estimate the linear time-invariant relationship between time series by using phase and magnitude information for all the time lags [8]. Both coherence and partial coherence were proved effective in overcoming an important limitation of the zero-lag correlation, that is, its sensitivity to the shape of the regional hemodynamic response function that could result in spurious correlations of the underlying neural activity. In the last two decades, there has been a growing interest in developing new connectivity metrics sensitive to both linear and nonlinear interactions in human brain. In fact, the spatiotemporal nonlinearity was shown to be an important feature of the BOLD signal that should be considered to properly characterize the complex interactions between brain regions. In [52] a phase space multivariate approach was adopted to investigate the nonlinear properties of resting state fMRI data. The dynamics of the signals were reconstructed by using the time delay embedding of some principal components of the fMRI data and the correlation dimension and the spatiotemporal Lyapunov exponents were calculated to assess the nonlinear fractal property and the chaotic dynamic behavior of the signals. A surrogate data test confirmed an inherent deterministic nonlinear behavior in fMRI fluctuations. Other methods for exploring the dynamic behavior of physiological signals have been proposed. Recurrence plots and recurrence quantitative analysis of the structures therein contained were used to examine the recurrence properties of dynamic systems [29]. As an example, in [53] RQA was employed as a univariate data-driven technique to quantify recurrent patterns in fMRI data. This technique involves the projection of each time series in the phase space from which a recurrence plot is obtained. Several numerical descriptors are then used to quantify recurrent patterns in each time series. This method has been developed as an alternative to general linear model and probabilistic independent component analysis in activation studies. The underlying idea is that single-voxel signals become more regular in response to a stimulus, so RQA can detect the most active voxels without any model assumption. Recurrence plots and RQA were proved successful in analyzing very noisy and nonstationary signals. These methods afford a set of metrics able to capture comprehensively the dynamic behavior of a system in the phase space. Some studies confirmed their effectiveness also for the analysis of EEG and MEG data, particularly for detecting functional anomalies in several diseases [54–57]. Cross recurrence plots are bivariate extensions of the recurrence plots that consist in two-dimensional matrices showing the interactions of pairs of signals in the phase space. The proposed index, extracted from the CRP, represents an intuitively interpretable generalized dynamic synchronization metric that could be used to extend the set of known RQA measures.


## 6. Methodological Limitations

These results are promising with respect to the value of the novel technique we are proposing, even though they are not free of limitations. For example, since we used task-dependent time series, we do not know yet whether these results extend to resting state data, and this will be the object of future studies. We chose to examine task-driven functional connectivity as done in several other studies [17, 18, 51], by analyzing the modular architecture during working memory. In particular, we used both the presence of known task-related functional modules and their high consistence across a healthy cohort of subjects to evaluate the proposed synchronization metric. An advantage of the block-designed task we considered is that BOLD activity presents cyclostationary properties due to the ON-OFF periods of the task. Instead, spontaneous BOLD fluctuations are intrinsically dynamic over time and thus nonstationary [58]. For this reason, studying modularity with resting state data will require a modified dynamical framework to correctly identify stable ICNs for the declared purposes.Another relevant issue concerns the modularity properties used to perform the comparison between the SYNC metric and Pearson’s correlation index. Indeed, in our analysis we found both higher modularity and higher consistency of task-related communities in the SYNC matrices. These features are related to a greater homogeneity of the functional organization across the subjects in response to the same task and although they are compatible with behaviors expected in a healthy cohort, a more rigorous assessment of the sensitivity of the proposed synchronization metric should require further analysis. Future studies could employ alternative topological properties of SYNC networks and their correlation with task performance or behavioral data to uncover additional insights into the suitability of the SYNC index as a functional connectivity metric for fMRI time series.Finally, our study has focused on an alternative method to define functional connectivity between pairs of BOLD time series. Generally, functional connectivity refers to a larger spectrum of neuroimaging techniques including EEG, MEG, and NIRS. As discussed above, recurrence plots have been used to explore dynamical properties of EEG and MEG, providing interesting features on complex phenomena in human brain. Although the SYNC metric is extracted from cross recurrence plots, a separate and accurate analysis may be needed to assert the validity of the index in a broader context and extend its use to more functional imaging techniques.


## 7. Conclusion

In this work, a new synchronization-based metric is proposed to assess functional connectivity in human brain. The metric is a generalized synchronization measure that takes into account both the amplitude and phase coupling between pairs of fMRI series. This method differs from the correlation measures used in the literature, as it is more sensitive to nonlinear coupling phenomena between time series and it is more robust against the physiological noise. In order to probe these latter two aspects, we performed a modularity analysis of task-related fMRI networks of a cohort of healthy subjects built with the new proposed metric. The aim was to verify whether the new metric was able to return networks whose functional modules were coherent with the actual organization of the brain regions during the task-based activity. We considered unthresholded complete connectivity matrices to test the effectiveness of the synchronization against noise and spurious correlations. Indeed unthresholded networks have lower signal-to-noise ratio as the most important links do not stand out among all the weights. By comparing the networks constructed by means of the proposed metric with those obtained through Pearson’s coefficient, it seems that the synchronization metric better reflects the task-related network structure for number of detected communities, for the functional organization of the ROIs, and for greater consistency of communities across the subjects.



---

*Source: 7190758-2017-10-30.xml*