Altered Topological Properties of Functional Network Connectivity in Schizophrenia during Resting State: A Small-World Brain Network Study

Aberrant topological properties of small-world human brain networks in patients with schizophrenia (SZ) have been documented in previous neuroimaging studies. Aberrant functional network connectivity (FNC, temporal relationships among independent component time courses) has also been found in SZ by a previous resting state functional magnetic resonance imaging (fMRI) study. However, no study has yet determined if topological properties of FNC are also altered in SZ. In this study, small-world network metrics of FNC during the resting state were examined in both healthy controls (HCs) and SZ subjects. FMRI data were obtained from 19 HCs and 19 SZ. Brain images were decomposed into independent components (ICs) by group independent component analysis (ICA). FNC maps were constructed via a partial correlation analysis of ICA time courses. A set of undirected graphs were built by thresholding the FNC maps and the small-world network metrics of these maps were evaluated. Our results demonstrated significantly altered topological properties of FNC in SZ relative to controls. In addition, topological measures of many ICs involving frontal, parietal, occipital and cerebellar areas were altered in SZ relative to controls. Specifically, topological measures of whole network and specific components in SZ were correlated with scores on the negative symptom scale of the Positive and Negative Symptom Scale (PANSS). These findings suggest that aberrant architecture of small-world brain topology in SZ consists of ICA temporally coherent brain networks.

ICA, which was developed to solve problems similar to the ''cocktail party'' scenario in which individual voices must be resolved from microphone recordings of many people speaking simultaneously [19,24], has typically been applied to fMRI data by determining a set of maximally spatially independent brain networks each with associated time courses [25][26][27]. This approach is useful to examine brain activity during resting state in both healthy controls (HCs) and various patient groups [28][29][30][31][32]. Jafri et al [33] evaluated functional temporal connectivity among ICA component time courses of resting state fMRI data. The interrelationship among multiple brain networks (components) was defined as functional network connectivity (FNC). They examined differences in FNC between HCs and SZ and found greater occurrences of higher correlation among networks in SZ. However the topological properties of these FNC relationships have not yet been studied.
Studies implementing graph theory to neuroimaging data, especially fMRI, are growing rapidly [31,[34][35][36][37][38]. Small-world properties are consistently revealed in human brain networks, which may suggest that the brain generates and integrates information with high efficiency [39][40][41][42][43][44]. Several studies document aberrant small-world network metrics in SZ [45][46][47]. For example, Yu et al [48] discovered altered topological measures in auditory oddball task-related small-world brain networks in SZ. Bloch et al [49] evaluated disrupted modularity in SZ. Several other resting state fMRI studies revealed less hierarchy, less smallworld properties, less clustering and less efficiency in SZ [34,50,51]. However, all these results were obtained from networks constructed based on brain voxels-or regions-of-interest [52]. The topological properties of brain networks consisting of spatial components in SZ versus HC are yet not known.
The aim of the present study was to identify such differences of small-world network measures in FNC between fMRI data acquired in HCs and SZ during the resting state. We hypothesized that SZ would show abnormal topological properties in this kind of small-world brain networks based on previous studies [33,[53][54][55].

Ethics Statement
This study has been approved by Hartford Hospital and Yale ethics committee. All subjects gave their written informed consent.

Participants
Subjects consisted of 19 ( [56] administered by a research nurse and by review of the medical records. All patients had chronic schizophrenia (positive and negative syndrome scale (PANSS) [57]: positive score 1665; negative score 1565) and all were taking medication (including the atypical antipsychotic medications aripiprazole, clozapine, risperidone, quetiapine and olanzapine, first-generation antipsychotics including fluphenazine, and miscellaneous moodstabilizing, hypnotic and anti-cholinergic medications including zolpidem, zaleplon, lorazepam, benztropine, divalproex, trazodone, clonazepam). All participants except 1 healthy control and 2 patients were right-handed. Exclusion criteria included auditory or visual impairment, mental retardation (full scale IQ,70), traumatic brain injury with loss of consciousness greater than 15 min, and presence or history of any central nervous system (CNS) neurological illness. Participants were also excluded if they met criteria for alcohol or drug dependence within the past 6 months or showed a positive urine toxicology screen (screening was for cocaine, opioids including methadone, cannabis, amphetamine, barbiturates, PCP, propoxyphene, and benzodiazepines) on the day of scanning. Healthy participants were free of any DSM-IV TR Axis I disorder or psychotropic medication and had no family histories of Axis I disorders.

Image acquisition
One 5-min resting state run for each subject was acquired at the Olin Neuropsychiatry Research Center at the Institute of Living/ Hartford Hospital on a Siemens Allegra 3T dedicated head scanner equipped with 40 mT/m gradients and a standard quadrature head coil. The functional scans were acquired transaxially using gradientecho echo-planar-imaging with the following parameters: repeat time (TR) 1.50 s, echo time (TE) 27 ms, field of view 24 cm, acquisition matrix 64664, flip angle 70u, voxel size 3.756 3.7564 mm 3 , slice thickness 4 mm, gap 1 mm, 29 slices, ascending acquisition. Six ''dummy'' scans were acquired at the beginning to allow for longitudinal equilibrium, after which the paradigm was automatically triggered to start by the scanner.

Preprocessing
FMRI Data were preprocessed using the SPM5 (http://www.fil. ion.ucl.ac.uk/spm/software/spm5/) software package. Data were motion corrected using INRIalign-a motion correction algorithm unbiased by local signal changes [58], spatially normalized into the standard Montreal Neurological Institute (MNI) space, and spatially smoothed with a 10610610 mm 3 full width at halfmaximum Gaussian kernel. Following spatial normalization, the data (originally acquired at 3.7563.7564 mm 3 ) were resliced to 36363 mm 3 , resulting in 53663646 voxels.

Group ICA
Group spatial ICA [55] was conducted for all the data using the infomax algorithm [24]. Data of all 38 participants were decomposed into 75 components using the GIFT software (http://icatb.sourceforge.net/). Single subject time courses and spatial maps were then back-reconstructed [59,60]. We chose the relatively high model order ICA as previous studies demonstrated that such models yield refined components which correspond to known anatomical and functional segmentations [61][62][63][64]. The Infomax ICA algorithm was repeated 10 times in ICASSO (http://www.cis.hut.fi/projects/ica/icasso) and resulting components were clustered to estimate the reliability of the decomposition. Fifty-seven components that did not contain large edge effects or ventricles by visual inspection were selected for further analysis. Temporal band-pass filtering (0.01,f,0.10 Hz) [65,66] was performed for component time courses of each subject before computing partial correlations.

Estimation of inter-component partial correlations
Partial correlation are useful as a measure of connectivity between a given pair of brain regions because they attenuate the contribution of other sources of covariance [67,68]. Given a set of N random variables, the partial correlation matrix is symmetric, where each off-diagonal element is the correlation coefficient between a pair of variables after filtering out the contributions of all other variables included in the dataset [50]. In this study, we evaluated temporal connectivity between each pair of ICs using partial correlation of ICA time courses to reduce the effects of the other 55 brain networks [69], and built undirected graphs respectively for each subject.
The first step was to estimate the sample covariance matrix S from the data matrix Y = (x i ), i = 1,…,57, of observations for each individual. Here x i was the time courses of each i th component. If we introduce X = (x j , x k ) to denote the observations in the j th and k th components, Z = Y\X denotes the other 55 time courses matrices. Each component of S contains the sample covariance value between two components (say j and k). If the covariance matrix of [X, Z] was in which S 11 was the covariance matrix of X, S 12 was the covariance matrix of X and Z and S 22 was the covariance matrix of Z, then the partial correlation matrix of X, controlling for Z, could be defined formally as a normalized version of the covariance matrix [50], Finally, a Fisher's r-to-z transformation [50,70,71] was applied to the partial correlation matrix in order to induce normality on the partial correlation coefficients. Connectivity strength [45,50] of network was computed by absolute z values. Two sample two-tailed t-test was performed to test for group difference in the strength of functional connectivity.

Graph analysis
An N6N (N = 57 in the present study) binary graph brain network, G, consisting of nodes (brain components) and undirected edges (connectivity) between nodes, could be constructed by applying a correlation threshold T (Fisher's r-to-z) to the partial correlation coefficients: That is, if the absolute z(i, j) (Fisher r-to-z of the partial correlation coefficient) of a pair of components, i and j, exceeds a given threshold T, an edge is said to exist; otherwise it does not exist. We defined the sub-graph G i as the set of nodes that were the direct neighbors of the i th node, i.e. directly connected to the i th node with an edge [50]. The degree of each node, K i, i = 1,2,…,57, was defined as the number of nodes in the sub-graph G i . The degree of connectivity, K net , of a graph is the average of the degrees of all the nodes in the graph: which is a measure to evaluate the degree of sparsity of a network. The total number of edges in a graph, divided by the maximum possible number of edges N(N-1)/2: is called the cost (connection density) of the network.
The clustering coefficient of a node was the ratio of the number of existing connections to the number of all possible connections in the subgraph G i is: where E i is the number of edges in the sub-graph G i [72,73]. The clustering coefficient of a network is the average of the clustering coefficients of all nodes: where C net is a measure of the extent of the local density or cliquishness of the network. The mean shortest path length of a node was: in which min{L i,j } is the shortest path length between the i th node and the j th node, and the path length was the number of edges included in the path connecting two nodes. The characteristic path length of a network is the average of the shortest path lengths between the nodes: L net is a measure of the extent of average connectivity of the network. E global , a measure of the global efficiency of parallel information transfer in the network, is defined as the inverse of the harmonic mean of the minimum path length between each pair of nodes [43,74,75]: The local efficiency of the i th node can be calculated as follows: In fact, since the i th node is not an element of the subgraph G i , the local efficiency can also be understood as a measure of the fault tolerance of the network, indicating how well each subgraph exchanges information when the index node was eliminated [43].
In addition, based on its definition, it is a measure of the global efficiency of the subgraph G i . The mean local efficiency of a graph, is the mean of all the local efficiencies of the nodes in the graph. For more information about the uses and interpretations of the complex brain network measures to see Rubinov and Sporns [35].

Small-Worldness
Compared with random networks, small-world networks have similar path lengths but higher clustering coefficients, that is c = C net /C random .1, l = L net /L random <1 [73]. These two conditions can also be summarized into a scalar quantitative measurement, small-worldness, s = c/l, which is typically .1 for small-world networks [41,76,77]. To examine the small-world properties, the values of C net and L net of the functional brain network need to be compared with those of random networks. The theoretical values of these two measures for random networks are C random = K/N, and L random <ln(N)/ln(K) [42,76,78]. However, as suggested by Stam et al [78], statistical comparisons should generally be performed between networks that have equal (or at least similar) degree sequences; whereas theoretical random networks have Gaussian degree distributions that may differ from the degree distribution of the brain networks that we discovered in this study. To obtain a better control for the functional brain networks, we generated 100 random networks for each degree K of each individual network by the same Markov-chain algorithm [44,79,80] used in previous studies [50,81]. And averaged across all 100 generated random networks to obtain a mean C random and a mean L random for each cost K cost .

Small-World Regime
Since the topological indices of the FNC are computed as a function of threshold, the results can be influenced by differences in the number of edges between the two groups [50,78]. In this study, we chose thresholds to keep the same number of edges in networks of all participants and investigated the topological properties of the FNC as a function of cost, K cost_min #K cost #K cost_max . The range was determined following Liu et al [50] and Liao et al [50,81]: (1) the minimum cost assures that each network was fully connected with N = 57 nodes. This allowed us to compare the topological properties between the two groups in a way that was relatively independent of the size of the network. (2) the maximum cost was selected to ensure that the brain networks have a lower global efficiency and a larger local efficiency compared to random networks with relatively the same distribution of the degree of connectivity [43]. We selected the small-world regime as 0.351#K cost #0.417 and repeated the full analysis for each value of K cost with increments of 0.004 (the corresponding number of edges was 560-665 with steps of 7 and the corresponding degree was 19.649-23.333 with steps of 0.246). The range of K cost was higher than prior studies might indicate standard deviation of partial correlation among brain ICs is higher than among brain regions. The corresponding z-score threshold value range for all of the subjects was 0.244-0.513.

Statistical Analysis
Following Liu et al [50], two sample two-tailed t-tests were performed to test for group differences in C net , L net , E global , E local , c, l and s at each of the 16 selected K cost values. A statistical significance level of P,0.05 was used. 1660.05 = 0.80,1 which means there was less than one false-positive result at this [45]. If a significant difference was found, Pearson linear correlation coefficients were used to evaluate the relationship between topological indices and PANSS for each selected cost point in the SZ group. In addition, as each node of the network had its own value of degree (K node ), cluster coefficient (C node ), shortest path length (L node ), global efficiency (E node_global ), and local efficiency (E node_local ), the distribution of components which showed significant differences in these measures were investigated by post-hoc two sample t-tests (P,0.05). If significant difference was found for any brain component, Pearson linear correlation coefficients were used to examine the relationship between topological properties of that brain component and PANSS scores in the SZ group.

Effect of number of components on small-world results
To explore the effect of ICs' number on results of small-world networks, C net , L net , E local , and E global were computed in HCs and SZ after building networks by changing the number of components as follow. Three more times of group ICA were performed by splitting the fMRI data into 65, 80, 85 ICs. 52, 56, 59 components of interest were selected to build small-world networks respectively. Small-world regimes were 0.324#K cost #0.392 with steps of 0.004; 0.364#K cost #0.432 with steps of 0.004; and 0.317#K cost #0.396 with steps of 0.005. Corresponding number of edges in the networks were 430-520 with steps of 6; 560-665 with steps of 7; and 542-677 with steps of 9.

Small-world network metrics
For topological indices as a function of cost, as cost increases, the clustering coefficient, local efficiency and global efficiency also increase whereas the characteristic path length decreases because more and more edges are added into the network. As shown in Figure 2, two sample t-tests indicated that clustering coefficients were higher (P,0.05, uncorrected) in SZ at most of the cost points; local efficiencies were higher in SZ only at higher cost values (7 highest values of the selected cost); whereas the character path length and global efficiency showed significant group difference at lower cost values (4 lowest values of the selected cost), SZ had higher characteristic path length and lower global efficiency. After changing the number of components (52,56,59 ICs from total of 65,80,85), SZ were showing the same trend for all the network properties. And the group differences of C net and E local at some K cost points were significant (P,0.05, uncorrected) or marginally significant (P,0.1, uncorrected).
The small-world attribute was evident in the FNC of both groups. c was significantly greater than 1 while l was near to 1 over the whole range of K cost . Because C random and L random show similar group differences as Cnet and Lnet (see Figure 2), no statistical significant differences in the values of c, l or s between the two groups were found when the same cost was used.

Relationships between network measures and PANSS scores in SZ
Significant (P,0.05, uncorrected) correlations between topological metrics and PANSS values in SZ were only found between characteristic path length (L net ), global efficiency (E global ) and negative scale score by Pearson correlation analysis. L net was positively correlated with negative scale of PANSS whereas E global was negatively correlated with the PANSS negative scale. Figure 3 shows the patterns at a typical cost point (K cost = 0.390).

Distribution of components in which topological metrics are altered in SZ
When using two sample t-tests to explore statistical differences in the topological properties including degree, clustering coefficient, local efficiency, shortest path length and global efficiency between SZ and HC for each of the 57 components at each selected cost, we found significant differences (P,0.05, uncorrected) in eleven components (IC6, IC8, IC14, IC25, IC28, IC29, IC32, IC40, IC45, IC46 and IC52) which involve frontal, parietal, occipital and cerebellar brain areas ( Figure 4 shows the activation of these components). Although the results do not pass false discovery rate (FDR) correction (P,0.05) for multiple comparisons of 57 nodes, these tests are implemented after finding significant group differences in the metrics of the whole network. Thus the node specific comparisons are useful to report as post-hoc tests.

Relationships between network properties of individual components and PANSS scores in SZ
Based on Pearson correlation coefficients, shortest path length and local efficiency of IC14, shortest path length of IC32 were significantly (P,0.05, uncorrected) positively correlated with negative scale of PANSS. Patterns at a typical cost point (K cost = 0.351) were shown in Figure 7.

Discussion
In the present study, resting state fMRI data for both HC and SZ were decomposed into spatially independent components by group ICA. Undirected graphs were built based on partial correlation matrices which were computed using ICA time courses. The brain networks showed ''small-world'' patterns in the selected range of cost in both groups that was consistent with previous studies that found small-world features in complex brain networks consisting of brain regions [36,44]. The finding that topological properties were altered in SZ and that several metrics were significantly correlated with negative PANSS scores in SZ provides further evidence for aberrant FNC associated with this disease [33]. One novel aspect of this study is that all results were obtained from brain topology consisting of ICA temporally coherent brain networks.

Altered topological metrics in patients
Short characteristic path lengths and high global efficiencies have been demonstrated to promote effective interactions between and across different cortical regions [42,43,50]. In our study, characteristic path length showed significant higher value and global efficiency showed lower value only at lower cost points (see Figure 2), possibly indicating that information interactions between temporally interconnected brain components were slower and less efficient in SZ when the whole network of FNC was sparse. Clustering coefficients were significantly increased at most of the selected cost points and local efficiencies were increased only at the higher cost points in SZ (see Figure 2) in the networks constructed based on brain components in this study. This finding implied relatively dense local connectedness and robust local information processing of the brain networks in SZ [50,74,82] suggesting abnormal FNC in schizophrenia. The findings that clustering coefficients and local efficiencies were higher in SZ are consistent with previous studies which involved patients with other brain disorders. For example, Wang et al [83] found increased local efficiency in children with ADHD (attention-deficit/hyperactivity disorder); He et al [84] revealed increased clustering coefficient in Alzheimer's disease. However, as far as we know, all prior studies involving schizophrenia in which networks were constructed based on brain regions found less clustering coefficient and less local efficiency in SZ [34,45,50,51]. It is notable that this is the first work which finds increased clustering coefficient and local efficiency in schizophrenia. Moreover, that HCs and SZ showed similar small-worldness values distinguished our results from prior brain region-based graph studies [45,50]. Reasons for getting the distinct results are not clear, but the different graph building method (ICs-based graphs) of this study may play a role.

Distribution of altered brain components
Consistent with previous findings [45,50], only some nodes were altered in SZ. The network metrics of some components involving frontal, parietal, occipital and cerebellar brain areas were significantly altered in SZ (see Figure 4,5). For example, the degree of connectivity for IC8 (occipital region) was smaller in SZ indicating a lower connectivity between this component and other brain components. Our findings that SZ show aberrant connectivity in vision (IC6, IC8, IC14), motor (IC28, IC32, IC46), attention (IC25) and default (IC40, IC45, IC52) networks are consistent with previous studies which reported abnormal activation [85,86] or disturbed connectivity in those brain regions in SZ [33,[87][88][89].

Relationship between topological measurements and PANSS scores in SZ
Interestingly, we found characteristic path length and global efficiency of the whole brain network were correlated with PANSS negative scale values in SZ. Higher negative scale scores were associated with longer character path length and lower global efficiency ( Figure 3). These might indicate the more severe these symptoms, the lower information interactions among brain components. In addition, clustering coefficient, local efficiency of IC14 (occipital region) and shortest path length of IC32 (parietal region) were correlated with negative PANSS scores in SZ (see Figure 7). Higher negative scale scores were associated with a higher clustering coefficient, higher local efficiency of IC14 (occipital region) and longer shortest path length of IC32 (parietal region). These findings are in line with studies which found psychopathology is associated with aberrant intrinsic organization of functional brain networks in schizophrenia [21] and provide further evidence for this illness as a disconnection syndrome [50,[90][91][92][93].

Methodological limitations
The main limitation of this study is our use of a statistical significance level of P,0.05 (uncorrected). We did not use any stringent type I error control such as Bonferroni or false discovery rate correction especially for the results of nodal level analysis. Further studies could increase the statistical power by increasing the sample of subjects. The second possible confound is medication use in patients. It is not clear if medication in SZ alters graph parameters. In addition, this study examined the For nodal definition, compared with previous studies which used a predefined anatomical template, ICA-based nodes are datadriven. Each IC consisted of intrinsic connected brain voxels which share the same time course. By using this kind of graphs, altered topological metrics among intrinsic connected brain networks (but not predefined brain regions which are not intrinsic connected) were revealed. However, there are also some disadvantages by using ICs as graph nodes. For example, the potential for artificial splitting of a network into spurious subnetworks when the ICA dimensionality is highly chosen (here is 75 components based on 5 minutes of resting state data); the potential for lack of coverage of certain cortical regions or overrepresentation (spatial overlap) of certain other cortical regions; 57 non-artifactual components were selected arbitrarily. The effect of specified number of components on small-world results should be further studied by recruiting larger sample of subjects in future.

Conclusions
To our knowledge, this is the first study using ICA and graph theory methods to explore abnormal small-world brain network properties in SZ during the resting state. Small-world topological metrics including clustering coefficient, local efficiency, characteristic path length and global efficiency of the graph built based on FNC were altered in SZ. Characteristic path length and global efficiency measures were correlated with negative scale scores on the PANSS. The network parameters of some individual brain components involving frontal, parietal and occipital areas were disturbed in SZ. These findings provide further evidence for aberrant FNC in SZ and expand our understanding of brain disconnection in schizophrenia.