Voxel-based, brain-wide association study of aberrant functional connectivity in schizophrenia implicates thalamocortical circuitry

Background: Wernicke’s concept of ‘sejunction’ or aberrant associations among specialized brain regions is one of the earliest hypotheses attempting to explain the myriad of symptoms in psychotic disorders. Unbiased data mining of all possible brain-wide connections in large data sets is an essential first step in localizing these aberrant circuits. Methods: We analyzed functional connectivity using the largest resting-state neuroimaging data set reported to date in the schizophrenia literature (415 patients vs. 405 controls from UK, USA, Taiwan, and China). An exhaustive brain-wide association study at both regional and voxel-based levels enabled a continuous data-driven discovery of the key aberrant circuits in schizophrenia. Results: Results identify the thalamus as the key hub for altered functional networks in patients. Increased thalamus–primary somatosensory cortex connectivity was the most significant aberration in schizophrenia (P=10−18). Overall, a number of thalamic links with motor and sensory cortical regions showed increased connectivity in schizophrenia, whereas thalamo–frontal connectivity was weakened. Network changes were correlated with symptom severity and illness duration, and support vector machine analysis revealed discrimination accuracies of 73.53–80.92%. Conclusions: Widespread alterations in resting-state thalamocortical functional connectivity is likely to be a core feature of schizophrenia that contributes to the extensive sensory, motor, cognitive, and emotional impairments in this disorder. Changes in this schizophrenia-associated network could be a reliable mechanistic index to discriminate patients from healthy controls.


INTRODUCTION
One of the earliest mechanistic notions proposed to account for the myriad of symptoms seen in individuals with psychotic disorders is the concept of 'sejunction' (fragmentation of localizable association links) put forward by Wernicke in 1894. 1 Wernicke believed that a disjunction between distinct functional modules that involve both sensorimotor and association areas of the brain generate symptoms of psychosis. 2 The notion of sejunction has given rise to the current concept of dysconnectivity in schizophrenia, 3,4 resulting in numerous attempts to locate the hypothesized aberrant networks. Extant neuroimaging studies have highlighted a number of abnormal regional connections in patients contributing to the general acceptance of a vaguely defined 'widespread dysconnectivity'. However, the putative 'sejunction circuitry'-the most consistent and characteristic aberration in functional connections that define the syndrome of schizophrenia-is still elusive as no constant patterns that reliably explain schizophrenia's complex and heterogeneous symptoms have emerged to date. This poses a challenge to the reliability of schizophrenia as a diagnostic entity.
Notwithstanding its narrow diagnostic definition, schizophrenia is remarkably heterogeneous in outcome, and widespread variable patterns of changes in neural networks can contribute to its varied clinical expressions. When searching for final common pathways in any multifactorial complex disease, it is imperative that the noise induced by heterogeneity is tackled by increasing sample size. 5 Unfortunately, most studies investigating dysconnectivity in schizophrenia have modest sample sizes, often insufficient for observations to survive even lenient corrections for multiple comparisons. 4 To some extent, validation of results from small samples can be achieved through a meta-analytic approach. For example, consistent structural changes in the anterior insula and anterior cingulate cortex have emerged when morphometric data from more than 2,500 patients were synthesized by several groups. [6][7][8] Unfortunately, the two common methods used for studying functional connectivity-seed region approach and independent component analysis-do not allow unbiased pooling of results reported in individual studies. Seedbased analysis with a priori specification of brain regions is not a 'discovery approach', but allows for testing hypothesis that are cultivated by prior observations. The independent component analysis approach, though a data-driven approach well suited for novel discoveries of aberrant connectivity, is only partially independent of prior assumptions, [9][10][11] as the components are assumed to arise from statistically independent sources and are often selected on the basis of prior expectations of plausible signals, e.g., large-scale networks such as the default mode network. As a result, across studies, even similarly named components have a diverse anatomical distribution that again precludes pooled synthesis of individual studies.
The research into the pathophysiology of schizophrenia in general, and functional neuroimaging in particular, is plagued by partial replications and lack of convincing refutations. A summary of systematically selected resting-state connectivity studies on schizophrenia provided by Patterson-Yeo et al. 4 highlight this issue (an update of this summary is provided in Supplementary  Table S1). There is an urgent need to use methods that will allow large-scale pooling of data to both reduce the impact of heterogeneity and allow the study of stratified subgroups. To provide both greater confidence and accuracy in identifying which specific regions and their functional connections contribute most to schizophrenia, it is important to be able to use a voxelbased brain-wide analysis strategy, in addition to region-based. In this paper, we therefore report the first meta-study using both a region and voxel-based unbiased brain-wide association study (BWAS) approach on resting-state functional magnetic resonance imaging (MRI) data. Three resting-state schizophrenia data sets from China, USA, and the UK comprising of 415 patients and 405 matched controls were collated to produce the largest restingstate study in schizophrenia. The BWAS is modeled along the lines of the genome-wide associate study (GWAS) where large genetic data sets are pooled to identify significant genetic variations in specific disorders. The data-driven GWAS approach has contributed to firmly refuting a major role for certain genetic loci derived from small studies. The BWAS approach used here also uses a higher-order community discovery algorithm to account for the inherent modular nature of brain networks.
In addition, we investigate clinical associations among the so-called 'sejunction circuitry' (i.e., the most consistent abnormal pattern of connectivity discovered from the whole-brain datadriven search), symptom severity, and illness duration, and the extent to which it can reliably discriminate between patients and controls using a pattern classification approach.

MATERIALS AND METHODS Participants
Participants' demographic and clinical characteristics are summarized in Table 1. More details of subjects and the collection of these data are provided in the Supplementary Information section.

Image acquisition and preprocessing
All the Taiwan, COBRE, and Xiangya imaging data were acquired using a 3-T Siemens Trio Tim MRI scanner with an eight or a twelve channel phased array head coil, whereas that from Huaxi data were acquired using a 3-T General Electric MRI scanner. The Nottingham data set was acquired using a 3-T Philips Achieva MRI scanner. During data acquisition, the subjects were instructed to keep their eyes closed but not fall asleep (Huaxi, Taiwan, and Xiangya sites) or to remain with their eyes open staring at a fixation cross (COBRE and Nottingham sites). Image preprocessing steps included slice timing, within-subject realignment, spatial normalization to the stereotactic space of the Montreal Neurological Institute with voxel size of 3 × 3 × 3 mm 3 , linear detrending, band-pass filtering (0·01~0·08 Hz), and scrubbing. In all cases the data were smoothed spatially (FWHM 8 mm) and nuisance signal were regressed (including six motion parameters, the global, white matter, and cerebrospinal fluid signals). Any data affected by head motion (maximal motion between volumes in each direction, and rotation about each axis) of 43 mm or rotation of 43°was excluded. The Supplementary Information provides additional details of the imaging acquisition and data preprocessing.
Voxel-wise and atlas-based brain-wide association studies Step 1: analysis within each imaging center.
In the present study, each resting-state fMRI image included 47,636 voxels, which is based on the automated anatomical labeling atlas. For each pair of voxels, the time series were extracted and their correlation was calculated for each subject followed by z-transformation. Two-tailed, twosample t-tests were performed on the 1,134,570,430 (47,636 × 47,635/2) Fisher's z-transformed correlation coefficients to identify significantly altered functional links in schizophrenia patients compared with controls within each imaging center. The effect of age, gender ratios, antipsychotic dose, and head motion were regressed within each data set in this step.
Step 2: combination of results from all imaging centers. The Liptak-Stouffer z-score method, 12 which is a well-validated method for multisite data sets and has previously been used widely in multisite gene [13][14][15][16] and MRI data analysis 7,17 was then used to combine the results from the individual data sets. Specifically, the P value of each functional connectivity result from two-sample t-test in step 1 was converted to its corresponding z-score. This was calculated first as in equation: where Φ is the standard normal cumulative distribution function and i represent the i site. Next, a combined z-score for a functional connectivity was calculated using the Liptak-Stouffer formula: which follows a standard normal distribution under the null hypothesis; where w i ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi samplesize p is the weight of the i data set. Finally, The Z is transformed into its corresponding P value and a Bonferroni procedure was used to correct for multiple comparisons.
Step 3: calculating a measure for the association. A measure for the association (MA) between a voxel i and the brain disorder was then defined as: MA = N α , where N α is the number of links between voxel i and every other voxel in the brain that have a P value of less than α (in the present study α = 0.05/(47,363 × 47,635/2)) in t-tests. Voxel-level analysis of functional connectivity changes in schizophrenia W Cheng et al A larger value of MA implies a more significant alteration in functional connectivity.
To verify the voxel-wise results, we also parcellated the whole brain into 90 regions of interest according to the automated anatomical labeling atlas. 18 The names of the regions of interest and their corresponding abbreviations and anatomical classification 19,20 are listed in Table 2. The functional connectivity was evaluated between each pair of regions using a Pearson correlation coefficient.
Consistency of effect-size across studies was examined using the Cochran's Q test. 21 Between-study (effect) heterogeneity was indicated by a Q statistic P o0.05. The links with significant heterogeneity were additionally assessed using random-effects model. 22

Clinical correlates
We investigated whether the emerging pattern of dysconnectivity correlated with clinical variables (positive and negative syndrome scale, PANSS and illness duration) using partial correlation. The data without PANSS score (Nottingham) was excluded in our analysis. Specifically, we calculated the partial correlation between the strength of each altered functional connection and PANSS score and illness duration after removing the effect of sites. To further test the clinical correlates of the 'sejunction' concept in classification of schizophrenia, we applied a support vector machine 23 approach using the altered pairwise associations from BWAS to discriminate patients from controls. We used a leave-one-out crossvalidation strategy to estimate the generalizability of this classifier and estimated accuracy, sensitivity, and specificity. Permutation tests were used to estimate the statistical significance of the observed classification accuracy. Figure 1b, the region with the highest number of abnormal functional connections in schizophrenia patients is the thalamus, followed by the postcentral gyrus (PoCG). The smallest P value is around 10 − 18 , with a total of 9,044 altered links being significant after Bonferroni correction (Po 0.05, the significance level uncorrected had to be P o5.4 × 10 − 7 ). The most significantly altered cluster was in the thalamus (  Supplementary  Table S2.

Voxel-based BWAS As shown in
Atlas-based BWAS Using the automated anatomical labeling template, the most significantly changed region was again the thalamus: with a total of 33 changed links out of 74 links (Figure 2 and Supplementary  Table S3, Bonferroni correction P o0.05). The region-wise results are therefore consistent with the voxel-based meta-analysis. In the atlas-based 'sejunction map' plotted in Figure 2, one main hub is involved: the thalamus. Another notable feature is that the strength of functional links between the thalamus and the frontal cortical regions are mostly reduced in patients, other than those involving the motor and temporal cortical regions that are increased ( Figure 2). Moreover, the altered patterns of functional connections are very consistent across the data sets from the three different centers, despite the differences in the clinical populations: the mean functional links all increase or decrease in the same direction when comparing schizophrenia patients with healthy controls. More than 70% of altered links showed a consistent pattern; among the 29.7% links showing across-sites variation, the effect was due to a single center. Figure 2b show the altered patterns of top 11 functional connections that were most consistent across the data sets from the five sites. Another notable feature is that most changes involve links between the major functional divisions/lobes of the brain rather than within them, e.g., aberrant connectivity between subcortical regions and prefrontal or motor cortices, rather than among subcortical regions themselves. Clinical correlates of the altered circuit Studying the correlation between thalamocortical connections and clinical variables (PANSS scores and illness duration) at a false discovery rate-corrected P o 0.05, we noted several significant relationships. In general, lower thalamofrontal connectivity predicted higher burden of positive symptoms; higher thalamosensory (PoCG, STG) connectivity predicted higher burden of negative symptoms (Table 3). Interestingly, none of the altered nonthalamic links predicted any of the clinical variables tested (Supplementary  Table S4) at the corrected statistical threshold, indicating the primacy of thalamic aberrations in schizophrenia. Furthermore, we also performed correlation analyses between altered functional connectivities and all PANSS subscores and the result is summarized in Supplementary Figure S1. Interestingly, the thalamus-PoCG and thalamus-PCG functional connectivities were most strongly associated with difficulty in abstract thinking and thalamussuperior frontal gyrus connections with delusions. The support vector machine analysis ( Figure 3) indicated a high degree of accuracy in discriminating patients and controls on the basis of the aberrations in connectivity identified from the atlas-based BWAS (overall accuracy from all data sets pooled = 75.81%). The accuracy rates for distinguishing patients from healthy controls varied within a small range (73.5 to 80.9%, i.e.,7.4%) across the data sets. Permutation tests revealed that within and across data sets, discrimination accuracies were highly significant (Po0.01 in all cases). Support vector machine results are summarized in Supplementary Table S5.

DISCUSSION
We have presented a large systematic investigation of altered brain-wide functional connectivity in schizophrenia. To our knowledge, this is the first time that both a regional and voxelbased analysis that covers the entire brain with no a priori selections have been conducted in such a large multisite sample. The results have revealed highly replicable aberrations in connectivity involving thalamocortical connections. Importantly, the putative 'sejunction' that is mostly anchored upon thalamus varies with symptom severity as well as illness duration and provides a good discrimination accuracy between patients and controls (73.5-80.9%), despite variations in the data acquisition  Table 2. Results shown are obtained using a meta-analytical approach correcting for inter-site variation in effect-size, age, gender ratios as well as antipsychotic dose. AAL, automated anatomical labeling; BWAS, brain-wide association study.
Voxel-level analysis of functional connectivity changes in schizophrenia W Cheng et al Voxel-level analysis of functional connectivity changes in schizophrenia W Cheng et al 5 across several independent sites. This emphasizes the reliability of resting fMRI as an investigative tool to study the pathophysiology of heterogeneous syndromes such as schizophrenia.
Many previous studies investigating brain changes associated with schizophrenia have mainly been based on seed-based or independent component analysis approaches. Although these studies have suggested the importance of anomalies in the default mode network 24,25 and other intrinsic networks, 26-28 both resting-state and task-related studies have reported varying findings that has not converged on a single pattern. Nevertheless, extant literature on functional connectivity in schizophrenia points towards increased default mode network and decreased prefrontal connectivity, 29 but significantly biased by the excessive focus on default mode network for resting fMRI and prefrontal cortex for task-based studies. 4 Although the findings from each of these studies are likely to be potentially crucial to the pathophysiology of schizophrenia, our rigorous data-driven discovery approach using a very large sample suggests that the most robust pattern within the 'widespread dysconnectivity' involves the thalamocortical circuitry in schizophrenia.
Several fMRI findings support the idea that the thalamus is a key region in schizophrenia, [29][30][31][32][33][34] but all of these studies focus on thalamus only, thus failing to demonstrate the prominence of thalamic connections among the ubiquitous dysconnectivity seen in patients. Our finding is essentially a strong and data-driven verification of thalamic-cortical dysconnectivity in schizophrenia. Across two different approaches (voxel-wise versus atlas-based) and across five data sets, we not only show that a large effect pertaining to aberrant connectivity of the thalamus, but we also show that the abnormalities involving this structure is more prominent than all other large-scale connectivity aberrations in schizophrenia.
A distinct feature of the identified subcortical (thalamic) connections is that there is a pattern of frontal reduction and sensorimotor enhancement (PreCG, LING, FFG, PCL, and PoCG) in connectivity. This pattern was first reported in a study that included 62 subjects with schizophrenia 29 and has been replicated and shown to be a common feature across the two major psychotic disorders-bipolar disorder and schizophrenia using another seedbased fMRI analyses 32 although these samples were not sufficiently powered to detect correlations with specific clinical symptom scores. Our findings confirm that this apparently hierarchical distribution of thalamocortical abnormalities 32 is related to symptom burden in schizophrenia. This observation raises an interesting question of whether the failure to update/propagate prediction models upon processing sensory information 35,36 could be mapped on to the frontal reduction and sensorimotor enhancement pattern of thalamic dysconnectivity in patients. Voxel-level analysis of functional connectivity changes in schizophrenia W Cheng et al Experimental neuroimaging using magnetic stimulation has provided more direct evidence of thalamic dysconnectivity in schizophrenia. 37 Given the normally high degree of connectivity of thalamus within the brain, we cannot dismiss the possibility that thalamic aberrations are nonspecific in schizophrenia, and only reflect the generalized nature of network-level aberrations that arise in this illness. Nevertheless, the pattern of selective increase in certain thalamic connections while the others show a reduction suggests a putative illness-related imbalance. As highlighted by earlier authors, 29 this pattern is also consistent with the importance of this structure in developmental cortical maturation and the neurodevelopmental theories pertaining to schizophrenia.
Our analysis of correlations between the thalamo-frontal and thalamo-sensorimotor functional connectivities revealed a restricted pattern of associations with specific symptoms. Reduced functional connections between the thalamus and superior frontal gyrus were particularly associated with PANSS delusion scores, in agreement with a previous structural-based study. 38 On the other hand, increased thalamic functional connections with both the pre-and postcentral gyri were particularly associated with difficulty in abstract thinking (see Supplementary Figure S1). This latter finding may reflect a role for these aberrant connections in cenesthopathy and perhaps some disturbed aspects of embodied cognition in schizophrenia, 39 with the sensorimotor system contributing to processing abstract words and conceptual knowledge. 40,41 Several strengths and limitations must be borne in mind while interpreting this study. We examined a large multisite data set, with a meta-analytic approach to address inter-site variations. More commonly used validation approaches with multiple data sets test whether the same observation survives a type 1 error threshold across the data sets, but ignore the distribution of effect-sizes across the samples. Our meta-analytic approach considers the effect-size distribution to be of more importance than whether observations survive arbitrary thresholds set for data sets of varying sample sizes. In addition to voxel-wise search, we also used an anatomical parcellation scheme to obtain functional links and show a convergence. With regard to limitations, in line with a number of other studies, we used a correlation-based approach to infer brain connectivity; no causal influences can be assumed between the linked regions. Although including medication dose has been covaried when computing functional connectivity, we cannot rule out some possible nonlinear effects and effects due to lifetime exposure. Many of our patients were medicated and antipsychotic medications may attenuate functional connectivity patterns in the short term, though to date it is unclear how thalamocortical connectivity is affected by antipsychotics. 42,43 The observed dysconnectivity was consistent across sites, despite some samples being antipsychotic-naive (Huaxi), mixed (Xiangya), and mostly medicated (COBRE, Nottingham), suggesting that the effect of antipsychotics on the 'core' dysconnectivity in schizophrenia, if present, is weak. In addition to antipsychotics, we addressed other potential confounds including movement artefacts carefully (Supplementary Material).
In summary, we identified the probable 'sejunction circuitry' located on thalamic hubs in the largest reported sample of patients with schizophrenia to date. The implications of this finding are three-fold: First, by using a meta-analytic approach on readily available, anatomically parcellated resting-state data, we have provided an initial database upon which further studies could be pooled to continuously shape the search for the consistently dysconnected neural system in schizophrenia. While we do not claim that the frontal reduction and sensorimotor enhancement pattern of thalamocortical dysconnectivity is the core abnormality in schizophrenia, at a more modest level, we suggest that the BWAS approach can facilitate discovery akin to the GWAS, with a greater promise of success due to the proximity of neuronal changes to disease expression (as evidenced by the relationship between dysconnectivity and symptom burden). In addition to the large-scale genomic and enviromic data, such connectomic approaches should finally contribute to our understanding of mechanistic pathways that relate to the expression of the schizophrenia phenotype. Second, our data-driven results challenge some of the existing but inconsistent notions of dysconnectivity, highlighting the importance of an exhaustive bottom-up data mining in the investigation of schizophrenia. Finally, given the potential importance of thalamocortical connectivity to animal models and cognition in schizophrenia, our findings emphasize the role of agnostic data-driven neuroimaging in the pathways of drug discovery. The quest for the core mechanisms of schizophrenia has a long history; our observation clarifies the paths that are likely to be promising in this trail. Figure 3. ROC curve derived from the aberrant connectivity patterns emerging from the brain-wide association study using regional parcellations. Using the connectivity metrics as discriminant markers in a multivariate pattern classification framework, the true-positive rate is plotted against the false-positive rate. The black dots denote the best cut-off values. ROC, receiver operating characteristic.