Alterations of Functional and Structural Networks in Schizophrenia Patients with Auditory Verbal Hallucinations

Background: There have been many attempts at explaining the underlying neuropathological mechanisms of auditory verbal hallucinations (AVH) in schizophrenia on the basis of regional brain changes, with the most consistent findings being that AVH are associated with functional and structural impairments in auditory and speech-related regions. However, the human brain is a complex network and the global topological alterations specific to AVH in schizophrenia remain unclear. Methods: Thirty-five schizophrenia patients with AVH, 41 patients without AVH, and 50 healthy controls underwent resting-state functional magnetic resonance imaging (fMRI) and diffusion tensor imaging (DTI). The whole-brain functional and structural networks were constructed and analyzed using graph theoretical approaches. Inter-group differences in global network metrics (including small-world properties and network efficiency) were investigated. Results: We found that three groups had a typical small-world topology in both functional and structural networks. More importantly, schizophrenia patients with and without AVH exhibited common disruptions of functional networks, characterized by decreased clustering coefficient, global efficiency and local efficiency, and increased characteristic path length; structural networks of only schizophrenia patients with AVH showed increased characteristic path length compared with those of healthy controls. Conclusion: Our findings suggest that less “small-worldization” and lower network efficiency of functional networks may be an independent trait characteristic of schizophrenia, and regularization of structural networks may be the underlying pathological process engaged in schizophrenic AVH symptom expression.


INTRODUCTION
Hallucinations, which are defined as any perceptual experience in the absence of corresponding stimuli in the external world (David, 2004), are one of the most debilitating symptoms of schizophrenia. Although they may involve any of the senses, auditory verbal hallucinations (AVH) or "hearing voices" are the most prevalent type of hallucinations in schizophrenia (Andreasen and Flaum, 1991), and affect about 60-80% of schizophrenic patients. In about 25% of these patients, AVH are chronic and resistant to medication (Shergill et al., 1998), causing an impaired quality of life (Gaite et al., 2002).
In the past two decades, the advent of neuroimaging techniques has allowed researchers to explore the neurobiological basis of schizophrenic AVH in vivo (Allen et al., 2008(Allen et al., , 2012. Research into brain structural and functional alterations [including gray matter volume (Allen et al., 2008(Allen et al., , 2012Palaniyappan et al., 2012;Modinos et al., 2013), cortical thickness (van Swam et al., 2012), white matter integrity (Geoffroy et al., 2014;Wigand et al., 2015), state activation (Jardri et al., 2011;Kompus et al., 2011;Allen et al., 2012;Kuhn and Gallinat, 2012), and resting-state functional connectivity (Hoffman and Hampson, 2011)] related to AVH has been accumulating in recent years. Some key brain regions involving auditory processing, speech generation and speech perception have been identified. These AVH-related brain regions include the superior temporal gyrus, Wernicke's area, Broca's area, and arcuate fasciculus bundle. However, previous studies only focus on regional changes of the brain, and it is largely unknown whether there are global topological alterations specific to AVH in schizophrenia.
Recently, an increasing number of brain network studies based on graph theory have emerged using modern brain mapping techniques, such as functional magnetic resonance imaging (fMRI) and diffusion tensor imaging (DTI; Bullmore and Sporns, 2009). The application of graph theoretical approaches has provided a powerful tool to characterize topological properties of complex brain networks on a whole brain scale by defining a graph as a set of nodes (brain regions) and edges (functional or structural connections; Bullmore and Sporns, 2009;He and Evans, 2010;Bullmore and Bassett, 2011). Through these approaches, both structural and functional networks of human brain have been found to have a optimum small-world topology (Watts and Strogatz, 1998), characterized by a high local specialization and a high global integration between brain regions (Salvador et al., 2005;Achard et al., 2006;Hagmann et al., 2007;He et al., 2007). Global measures of brain networks are often represented in multiple ways, such as small-world properties and network efficiency, which variously quantify the capability of integration and segregation of information processing, test the global efficiency of parallel information transfer, and characterize the fault tolerance of the network (Rubinov and Sporns, 2010). There is an expanding body of evidence that patients with schizophrenia exhibited disrupted functional and structural brain networks, predominantly characterized by altered global topological properties (Fornito et al., 2012;Micheloyannis, 2012;van den Heuvel and Fornito, 2014).
Because AVH are one of the most characteristic and frequent symptoms with definite brain damage in schizophrenia patients, we hypothesized that there would be abnormality specific to AVH in global topological properties of brain networks. To test our hypothesis, we collected resting-state fMRI and DTI data from 35 schizophrenia patients with AVH, 41 patients without AVH, and 50 healthy controls. Functional and structural networks were constructed and analyzed using graph theoretical approaches.
Inter-group differences in global metrics of both networks were investigated.

Subjects
A total of 126 individuals were enrolled in the present study, including 76 patients with schizophrenia and 50 healthy controls. The diagnoses of schizophrenia were determined by the consensus of two experienced clinical psychiatrists using the Structured Clinical Interview for the DSM-IV Axis I Disorder, Patient Edition (SCID-P). All healthy controls were screened using the non-patient edition of the SCID (SCID-NP) to confirm an absence of psychiatric illnesses. The inclusion criteria for all participants were age (20-60 years) and right-handedness. The exclusion criteria for all participants were a history of head trauma with consciousness disturbances lasting more than 5 min, a history of drug or alcohol abuse, pregnancy, and any physical illness such as cardiovascular disease or neurological disorders, as diagnosed by an interview and medical records review. Additional exclusion criteria for healthy controls included a history of psychiatric disease and first-degree relatives who had a history of psychotic episodes. According to the experience of AVH, schizophrenia patients were subdivided into two groups: AVH group (n = 35) included patients who experienced AVH at least once daily, and nAVH group (n = 41) included patients who had never experienced AVH or had not experienced AVH within 12 months prior to MRI scanning. Clinical symptoms of psychosis were quantified using the Positive and Negative Syndrome Scale (PANSS; Kay et al., 1987). The Auditory Hallucination Rating Scale (AHRS; Hoffman et al., 2003) was used to assess AVH on seven characteristics: frequency, reality, loudness, number of voices, length, attention dedicated to the hallucinations, and hallucination-induced arousal. The study was conducted in accordance with the Declaration of Helsinki and was approved by the Medical Research Ethics Committee of Tianjin Medical University General Hospital. After a complete description of the study, written informed consent was obtained from each subject.

fMRI Data Preprocessing
Resting-state BOLD data were preprocessed using Statistical Parametric Mapping 8 (SPM8 1 ). The first 10 volumes for each participant were discarded to allow the signal to reach equilibrium and the participants to adapt to the scanning noise. The remaining volumes were corrected for the acquisition time delay between slices. Then, realignment was performed to correct the motion between time points. All participants' BOLD data were within the defined motion thresholds (i.e., translational or rotational motion parameters less than 2 mm or 2 • ). We also calculated frame-wise displacement (FD), which indexes the volume-to-volume changes in head position. There were no significant differences in mean FD (one-way ANOVA, F = 0.421, P = 0.657) among the AVH (0.106 ± 0.046), nAVH (0.112 ± 0.057), and control (0.102 ± 0.052) groups. Several nuisance covariates (six motion parameters, their first time derivations, the global brain signal, the white matter signal, and the cerebrospinal fluid signal) were regressed out from the data. Recent studies have reported that the signal spike caused by head motion significantly contaminated the final resting-state fMRI results even after regressing out the linear motion parameters (Power et al., 2012). Therefore, we further regressed out spike volumes when the FD of the specific volume exceeded 0.5. The datasets were then band-pass filtered in a frequency range of 0.01-0.08 Hz. In the normalization step, individual structural images were firstly co-registered with the mean functional image; then the transformed structural images were segmented and normalized to the Montreal Neurological Institute (MNI) space using a high-level non-linear warping algorithm, that is, the diffeomorphic anatomical registration through the exponentiated Lie algebra (DARTEL) technique (Ashburner, 2007). Finally, each filtered functional volume was spatially normalized to MNI space using the deformation parameters estimated during the above two steps and resampled into a 3-mm cubic voxel.

Functional Network Construction
GRETNA software 2 was used to construct the whole-brain networks (Wang et al., 2015b). Nodes and edges are two basic elements of a brain network. To define the nodes of brain networks, we used the automated anatomical labeling (AAL) template (Tzourio-Mazoyer et al., 2002) to parcellate the whole brain into 90 (45 for each hemisphere) cortical and subcortical regions of interest (ROIs). For each subject, the representative mean time series of each ROI was obtained by averaging the BOLD time series over all voxels within that region. To define the edges of functional brain networks, we computed the Pearson correlation coefficients and their significance levels (i.e., P-values) between the regional mean time series of all possible pairs of nodes, resulting in a 90 × 90 correlation matrix for each subject (Figure 1). Because of the ambiguous interpretation of negative correlations (Fox et al., 2009;Murphy et al., 2009), we restricted our analysis to positive correlations. To further de-noise spurious interregional correlations, only those correlations whose corresponding significance levels survived a statistical threshold of P < 0.05 (Bonferroni corrected) were retained. This significance level-defined correlation threshold method has been used in previous brain network studies Wang et al., 2013Wang et al., , 2015a. Finally, the functional weighted networks were constructed with considering functional connectivity as the weight of each edge. To test the effects of the choices of correlation thresholds on our network analyses, we also constructed the functional networks using other two different correlation thresholds, that is, a looser one whose corresponding significance level survived a statistical threshold of P < 0.001 (uncorrected) and a stricter one whose corresponding significance level survived a statistical threshold of P < 0.01 (Bonferroni corrected).

DTI Data Preprocessing
The software packages FMRIB Software Library (FSL 3 ) and Diffusion Toolkit (DTK 4 ) were used for the DTI preprocessing steps. The diffusion-weighted images were first registered to a reference volume (i.e., the first b0 image) by using affine transformations to minimize distortions caused by the eddy currents and head motions. Then, the three-dimensional maps of the diffusion tensor and the fractional anisotropy (FA) were calculated by using DTK. The whole-brain deterministic fiber tracking was performed via Fiber Assignment by Continuous Tracking (FACT) algorithm (Mori et al., 1999) with the FA threshold of 0.2 and the angle threshold of 50 • .

Structural Network Construction
Because structural networks were constructed in the diffusion image native space, the AAL template in MNI space was needed to be transformed to individual native space. Briefly, individual structural images were firstly co-registered to their first b0 images. Then the transformed structural images were segmented and normalized to the MNI space using a non-linear warping algorithm, i.e., the DARTEL technique. The derived deformation parameters were inverted and used to transform the AAL template from MNI space to diffusion image native space. To define the edges of structural brain networks, we calculated the number of fibers (with end-points in both nodes during the fiber FIGURE 1 | Functional connectivity (correlation coefficient) and structural connectivity (fiber number) matrices for the construction of functional and structural weighted networks for AVH (A), nAVH (B), and HC (C). Both axes represent the 90 AAL brain regions used in the analysis. Color of each entry represents the level of connectivity. AVH, schizophrenia patients with auditory verbal hallucinations; nAVH, schizophrenia patients without auditory verbal hallucinations; HC, healthy controls.
tracking) between any pairs of nodes, resulting in a 90 × 90 fiber number (FN) matrix for each subject (Figure 1). Considering the limited resolution of DTI and the capacity of the FACT algorithm, there is a risk that some pseudoconnections will be included. To reduce influence of possible pseudoconnections, a threshold of three fibers was applied to all FN matrices, that is, two nodes were considered connected with an edge if at least three fibers existed Shu et al., 2009;Lo et al., 2010). Finally, the structural weighted networks were constructed with considering FN as the weight of each edge Bai et al., 2012;Baggio et al., 2015;Zhao et al., 2015).

Network Analysis
For both functional and structural weighted brain networks, global network metrics were calculated. Small-world properties of the brain networks are the most frequently used metrics, including clustering coefficient Cp, characteristic path length Lp, normalized clustering coefficient γ, normalized characteristic path length λ, and small-worldness σ (Watts and Strogatz, 1998). C p is a measure of the extent of the local density or cliquishness of the network, which characterizes network segregation. L p is a measure of the extent of average connectivity or overall routing efficiency of the network, which characterizes network integration. They were scaled against the mean clustering coefficient and characteristic path length obtained from 1000 matched random networks which had the same number of nodes, edges, and degree distribution as the real networks (Maslov and Sneppen, 2002), resulting in normalized clustering coefficient γ and normalized characteristic path length λ. These two measurements can also be summarized into a simple quantitative metric, small-worldness, σ = γ/λ. The brain networks have been found to have efficient smallworld properties in support of efficient transfer of parallel information at a relatively low cost (Achard and Bullmore, 2007). Thus, network efficiency is a biological relevant metric to describe brain networks from the perspective of information transfer (Latora and Marchiori, 2001;Achard and Bullmore, 2007). We calculated the network efficiency at global and local levels. The global efficiency Eg represents the capability of parallel information transfer over the network. The local efficiency Eloc reflects the fault tolerance of the network, which indicates how well the information is communicated within the neighbors of a given node when this node is eliminated. They were also scaled against the mean Eg and Eloc obtained from 1000 matched random networks, resulting in normalized Eg and normalized Eloc.

Statistical Analysis
To determine the inter-group differences in the global topological properties of the functional and structural networks, a one-way analysis of variance (ANOVA) were separately performed on each network metric (small-world properties and network efficiency) using non-parametric permutation tests. The non-parametric permutation tests were conducted using MATLAB language and the number of permutations was set to 10000. To examine the relationship between global topological properties and severity of AVH, a Spearmen's correlation analysis between each network metric and AHRS total score was performed in the AVH group. Because these analyses were exploratory in nature, we used a statistical significance level of P < 0.05. To validate our main findings, the same ANOVA was applied to network metrics of the functional networks constructed at the other two correlation thresholds (P < 0.001, uncorrected and P < 0.01, Bonferroni corrected), respectively.

Intergroup Differences in Global Metrics of Functional Networks
At the significance level-defined correlation threshold (i.e., P < 0.05, Bonferroni corrected), functional brain networks of the AVH, nAVH, and control groups exhibited higher clustering coefficients (γ AVH = 1.51 ± 0.21, γ nAVH = 1.55 ± 0.20, γ control = 1.52 ± 0.16) but almost identical characteristic path lengths (λ AVH = 1.10 ± 0.01, λ nAVH = 1.10 ± 0.02, λ control = 1.10 ± 0.01) relative to comparable random networks, which suggests that the three groups had typical small-world properties in functional networks. Despite common smallworld architecture of the functional networks, statistical analyses revealed significant differences in global metrics across the AVH, nAVH, and control groups (Figure 2). Compared with the control group, functional networks of both AVH and nAVH groups showed significantly decreased Cp, Eg, and Eloc, and increased Lp. However, no significant differences were observed in Cp, Lp, Eg, and Eloc between AVH and nAVH groups. In addition, we found that γ, λ, σ, normalized Eg and normalized Eloc had no significant differences across the three groups. Moreover, we found that our main results were reproducible after considering the effects of different correlation thresholds (Supplementary Figures S1 and S2). Inter-group differences in basic network properties of functional networks including density, connectivity strength and size of largest component are shown in Supplementary Material (SI Text) and Supplementary Figure S3.

Intergroup Differences in Global Metrics of Structural Networks
At the threshold of three fibers, structural brain networks of the AVH, nAVH, and control groups had higher clustering coefficients (γ AVH = 3.88 ± 0.32, γ nAVH = 3.91 ± 0.34, γ control = 3.86 ± 0.30) but almost identical characteristic path lengths (λ AVH = 1.12 ± 0.03, λ nAVH = 1.13 ± 0.03, λ control = 1.13 ± 0.03) relative to matched random networks, which indicates that the three groups exhibited a typical smallworld topology in structural networks. Differences in global metrics of structural networks across the AVH, nAVH, and control groups are shown in Figure 3. Compared with the control group, structural networks of AVH group exhibited significantly increased Lp. However, there was no significant difference in Lp between nAVH and control groups. Additionally, we found that Cp, γ, λ, σ, Eg, Eloc, normalized Eg and normalized Eloc had no significant differences across the three groups.
Inter-group differences in basic network properties of structural networks including density, connectivity strength and size of largest component are shown in Supplementary Material (SI Text) and Supplementary Figure S4.

Associations between Global Network Metrics and AVH Symptom Severity
For both functional and structural networks, we did not find any statistical correlations between global network metrics and ARHS total score in the AVH group.

DISCUSSION
To our knowledge, the present study is the first to investigate global topological alterations of functional and structural networks related to AVH in schizophrenia. We found that (1) schizophrenia patients with and without AVH exhibited common disruptions of functional networks, characterized by decreased clustering coefficient, global efficiency and local efficiency, and increased characteristic path length; (2) structural networks of only schizophrenia patients with AVH showed increased characteristic path length compared with those of healthy controls. These findings suggest that less "small-worldization" and lower network efficiency of functional networks may be an independent trait characteristic of schizophrenia, whereas regularization of structural networks may be the underlying pathological process engaged in schizophrenic AVH symptom expression. The human brain is structurally and functionally organized into a complex network, the main feature of which is a smallworld topology facilitating the segregation and integration of information processing . Complex network analysis based on graph theory has been applied to quantify brain networks with a rich array of neurobiologically meaningful and easily computable metrics (Rubinov and Sporns, 2010), such as the most frequently used small-world properties and network efficiency. According to graph theory, the brain can be depicted as a graph composed of a number of nodes interconnected by a set of edges. The nodes represent brain regions or even voxels. The edges represent functional connectivity measured by resting-state fMRI or fiber bundles tracked by diffusion tensor tractography. Several functional (Salvador et al., 2005;Achard et al., 2006) and diffusion (Hagmann et al., 2007;Gong et al., 2009) MRI studies have consistently demonstrated that brain networks in healthy subjects have small-world topologies and high efficiency at a low wiring cost. A combined analysis of fMRI-based functional and DTIbased structural networks would provide integrated information on the underlying neuropathological mechanisms of brain disorders.
In the current study, we observed that functional networks of both schizophrenia patients with and without AVH had less "small-worldization" (characterized by lower local specialization and global integration) and lower network efficiency compared with the healthy controls, indicating that aberrant functional networks may be an trait characteristic of schizophrenia regardless of the presence or absence of AVH. Our findings are in line with several previous functional network studies. For example, Micheloyannis et al. (2006) found that patients with schizophrenia showed smaller absolute clustering coefficients but relatively longer absolute path lengths in the alpha, beta and gamma bands during resting-state using electroencephalogram (EEG). In a resting-state fMRI study, Liu et al. (2008) reported that functional networks of schizophrenic patients exhibited altered small-world attributes characterized by decreased clustering coefficient, increased characteristic path length and decreased small-worldness, as well as decreased global and local efficiency (Liu et al., 2008). Brain networks with small-world properties confer resilience against pathological attack and impairments of these properties may lead to the development of schizophrenia. These findings are also consistent with the increasing evidence that schizophrenia is a disorder of abnormal brain connectivity, that is, dysfunctional integration among multiple, spatially distributed brain regions (Fitzsimmons et al., 2013;van den Heuvel and Fornito, 2014;Fornito and Bullmore, 2015).
Compared with healthy controls, structural networks of schizophrenia patients showed no significant changes in the FIGURE 2 | The differences in global metrics of the functional weighted networks across the AVH, nAVH, and HC groups. Error bars represent standard errors. Cp, clustering coefficient; Lp, characteristic path length; γ, normalized clustering coefficient; λ, normalized characteristic path length; σ, small-worldness; Eg, global efficiency; Eloc, local efficiency; AVH, schizophrenia patients with auditory verbal hallucinations; nAVH, schizophrenia patients without auditory verbal hallucinations; HC, healthy controls.
global network metrics, except that schizophrenia patients with AVH exhibited a lower global integration suggestive of a shift toward regular networks. A previous DTI study demonstrated that although a trend toward altered clustering, path length and small-worldness was found in schizophrenia patients, these differences were not significant (Zalesky et al., 2011), which is nearly consistent with our results. Different clinical profiles in patients (mixed vs. homogeneous) may account for the subtle differences between that study and ours. Our findings suggest that the regularization of structural networks might be the underlying neuropathology specific to AVH in schizophrenia. However, we did not find any correlations between global network metrics and AVH symptom severity, indicating that the network regularization process may be a stable pattern of schizophrenic AVH independent of its severity. However, we cannot exclude the possibility that there is a complex relationship between global network metrics and AVH symptom severity beyond a simple linear correlation.
The relationship between brain function and structure is traditionally illustrated by the hypothesis that neural activity is shaped, but not limited, by the underlying anatomy. In this study, however, converging evidence points out that functional and structural networks exhibit completely segregated changes related to AVH in schizophrenia. Functional networks are shown to be comprehensively disrupted and these impairments are specific to schizophrenia, indicating that functional network metrics are more sensitive biomarkers which can be used to detect whether schizophrenia exists. Structural networks present a slight regularization and this alteration is specific to AVH, indicating that structural network metrics are more specific biomarkers which can be used to predict and monitor AVH in schizophrenia.
There are several limitations to be considered in this study. First, most of the schizophrenia patients were receiving antipsychotic drug treatment. Although there was no difference in antipsychotic dosage between the schizophrenia patients with and without AVH, we cannot absolutely rule out the effect of antipsychotic medication on our results. Second, we did not assess the psychotic symptoms of patients during the MRI scan, which means that whether or not patients were experiencing AVH was unknown. Instead, AVH and other psychotic symptoms of all the schizophrenia patients were assessed before the MRI procedure. Third, the schizophrenic patients who had not experienced AVH within 12 months prior to MRI scanning were included in this study. These patients may have experienced AVH earlier in the psychotic illness, and they may still have trait characteristics of AVH. Fourth, the angular resolution of DTI was relatively low (25 directions) in this study, which may affect the robustness of the results of the tensor-based deterministic tracking method. Finally, we used AAL template to segment the whole brain into 90 regions for brain network construction. Recent studies have pointed out that the node definition by different parcellation strategies may lead to considerable variations in the graph theoretical metrics (Wang et al., 2009;Fornito et al., 2010;Hayasaka and Laurienti, 2010;Zalesky et al., 2010). In future studies, different parcellation schemes should be used to test the reproducibility of our findings.

CONCLUSION
The present study used a combined analysis of functional and structural networks to investigate global topological alterations specific to AVH in schizophrenia. We found that functional networks of schizophrenia patients had less "small-worldization" and lower network efficiency regardless of the presence or absence of AVH. More importantly, only schizophrenia patients with AVH exhibited regularization of structural networks, which may be the underlying pathological mechanisms of AVH in schizophrenia. These findings suggest that functional network measures might be used to determine whether schizophrenia exists, and structural network measures can be used to predict and monitor AVH in schizophrenia.

AUTHOR CONTRIBUTIONS
JZ, CW, and CZ designed the study. JZ and JL acquired the data, which JZ, FL, and WQ analyzed. JZ, CW, and CZ wrote the article, which all authors reviewed and approved for publication.