Author Correction: Functional connectome differences in individuals with hallucinations across the psychosis continuum

Hallucinations may arise from an imbalance between sensory and higher cognitive brain regions, reflected by alterations in functional connectivity. It is unknown whether hallucinations across the psychosis continuum exhibit similar alterations in functional connectivity, suggesting a common neural mechanism, or whether different mechanisms link to hallucinations across phenotypes. We acquired resting-state functional MRI scans of 483 participants, including 40 non-clinical individuals with hallucinations, 99 schizophrenia patients with hallucinations, 74 bipolar-I disorder patients with hallucinations, 42 bipolar-I disorder patients without hallucinations, and 228 healthy controls. The weighted connectivity matrices were compared using network-based statistics. Non-clinical individuals with hallucinations and schizophrenia patients with hallucinations exhibited increased connectivity, mainly among fronto-temporal and fronto-insula/cingulate areas compared to controls (P < 0.001 adjusted). Differential effects were observed for bipolar-I disorder patients with hallucinations versus controls, mainly characterized by decreased connectivity between fronto-temporal and fronto-striatal areas (P = 0.012 adjusted). No connectivity alterations were found between bipolar-I disorder patients without hallucinations and controls. Our results support the notion that hallucinations in non-clinical individuals and schizophrenia patients are related to altered interactions between sensory and higher-order cognitive brain regions. However, a different dysconnectivity pattern was observed for bipolar-I disorder patients with hallucinations, which implies a different neural mechanism across the psychosis continuum.


Supplementary Content
Functional connectome differences in individuals with hallucinations across the psychosis continuum. Sommer.

Supplementary Discussion
Supplementary    Table 6 Results of the network based statistics for non-clinical individuals with hallucinations versus healthy controls using the AAL atlas Supplementary Table 7 Results of the network based statistics for schizophrenia patients with hallucinations versus healthy controls using the AAL atlas Supplementary Table 8 Results of the network based statistics for bipolar-I disorder patients with hallucinations versus healthy controls using the AAL atlas Supplementary Table 9 Results of the network based statistics across all three hallucination groups using the AAL atlas Supplementary Table 10 Results of the network based statistics for non-clinical individuals with hallucinations versus healthy controls using the Power atlas.

Supplementary Table 11
Results of the network based statistics for schizophrenia patients with hallucinations versus healthy controls using the Power atlas.

Supplementary Table 12
Results of the network based statistics for bipolar-I disorder patients with hallucinations versus healthy controls using the Power atlas.

Supplementary Table 13
List of features for random forest classifier in non-clinical individuals with hallucinations.

Supplementary Table 14
List of features for random forest classifier in schizophrenia patients with hallucinations.

Supplementary Table 15
List of features for random forest classifier in bipolar-I disorder patients with hallucinations.

Supplementary Table 16
List of abbreviations for brain areas AAL-atlas

Data acquisition
All participants were acquired on the same 3 tesla Achieva Philips clinical MRI scanner equipped with an 8-channel SENSE head coil at the University Medical Center Utrecht (Philips, Best, The Netherlands). Anatomical 3D highresolution T1-weighted images were obtained for anatomical reference with the following parameter settings: echo time [TE] = 4.6 ms, repetition time [TR] = 10 ms, flip angle = 90º, Field of View (FOV) = 240 mm/100%, voxel size 0.75 × 0.75 × 0.80 mm, reconstruction matrix = 200 x 320 x 320. For 49 participants included in the Spectrum study [4], T1-weighted images with a lower resolution were available: 160 contiguous sagittal slices (TE = 4.6 ms, TR = 10 ms, flip angle = 8°, FOV = 224 mm, 1 × 1 × 1 mm voxels). Task-free resting state functional data were acquired using a 3D PRESTO sequence that allowed for fast functional brain coverage every 609 ms [11]. All participants were instructed to keep their eyes closed. Parameters settings were: (40 (coronal) slices, TE = 32.4 ms, TR = 21.75 ms, flip angle = 10°, FOV = 224 x 256 x 160, 4 x 4 x 4 mm voxels). Using this sequence, between 600-1000 images were obtained depending on the study for which the data was acquired, after which all resting state scans were resized to the first 600 images (~6 minutes). All scans were checked for radiological abnormalities.

Data preprocessing
Preprocessing was done using the FMRIB Software Library (FSL) version 5.0.4 [12] with FEAT's default settings, including skull stripping (BET), motion correction with MCFLIRT, spatial smoothing (5mm kernel at full width at half maximum) and high-pass filtering (100-second cut-off). The mean global signal was not included, as this could potentially bias group differences or increase the risk of removing a valuable signal [13][14][15]. In-scanner motion can influence functional connectivity measures and is known to disproportionally affect the functional connectome of patients versus controls [16][17]. Several steps were taken to prevent a systematic motion-related bias in this study following recent developments in the literature. If the relative root mean square displacement over all frames exceeded 0.2mm; or if 20 individual frames exceeded the threshold of 0.25mm, these subjects were excluded from further analysis [18]. ICA-AROMA was subsequently applied for denoising [19][20]. More information on exclusion of scans due to missing clinical data, poor imaging, data quality, processing errors, motion artefacts, or brain abnormalities reported by the radiologist can be found in the Supplementary material.

Functional network reconstruction
The average time series for each participant were extracted from N=90 functional regions (nodes) of the automated anatomical labeling (AAL) atlas [21] and N=260 nodes of the Power atlas [22] covering the whole brain, including cortical and subcortical regions. We applied wavelet decomposition on the raw time series to extract a signal in the frequency domain of 0.05 -0.1 Hz (Scale 4), as resting state studies commonly apply bandpass filtering over the frequency domain between 0.01 -0.1 Hz [23-25]. Wavelet-based methods have several advantages in terms of denoising [26], robustness to outliers [27], autocorrelations and has the advantage of bypassing the role of positive and negative edge weights compared the traditional band-pass filtering and Pearson correlation coefficients [28]. Because a coherence value is always between 0 and 1, this facilitates possible problems of how to deal with negative edge weights. The coherence was calculated as the mean squared coherence between time series over the chosen wavelet scale in Matlab, using the mscohere function.

Global connectome disturbances
We analyzed disturbances in global network organization related to hallucinations. Two key properties were investigated to characterize network topology. The first property is the weighted global efficiency, which provides information on overall communication efficiency throughout the network [29] and is hence a measure of network integration. It was calculated as the average inverse shortest path length between all node pairs [29]. The second property is the weighted clustering coefficient, which reflects the degree to which nodes in the network tend to cluster together, and is a metric of network segregation. It was calculated as the weighted level of functional connectivity that exists between each nodes' neighbours, proportional to the maximum degree of connectivity [29]. The graph metrics were computed using the Brain Connectivity Toolbox [29]. Thresholding connectivity matrices can bias graph metrics due to differences in edge density and functional connectivity strength between patient and control groups [30-31]. We therefore used unthresholded connectivity matrices and tested if the connectivity strength differed across groups.
Minimum Spanning Tree (MST) analysis provides information on a network's topology and is therefore insensitive to connection strength [32][33]. The Minimum Spanning Tree (MST) is a subnetwork of the original weighted connectivity matrix, connecting all nodes in the network without forming loops, and therefore represents the strongest connections in the network (e.g. a maximum spanning tree) [32]. For each subject, the MST network was constructed by consecutively adding connections with the highest-ranking weight. When adding a new connection resulted in the formation of a loop, this connection was omitted and the next connection in order of weight was added to the network. This procedure was repeated until all nodes were connected, resulting in a unique MST for each subject with a fixed number of (N) 90 nodes and (M = N-1) 89 links. Hence, no bias was induced due to differences in connection density or strength between subjects. The MST therefore represents the strongest connections in the network (e.g. a maximum spanning tree), and network topology was compared using the MST leaf fraction and MST diameter [32].
A combination of two MST measures, namely the diameter and the leaf fraction can be used to determine if the network's topology is more line-like (i.e., less integrated) or more star-like (i.e., more integrated) network [32]. The MST diameter is calculated as the longest distance (in number of edges) between any two nodes of the MST. The MST leaf fraction, reflecting the number of nodes that only link to one other node (i.e., degree 1). A longer diameter and lower leaf number will reflect a more line-like network, whereas a more integrated star-like network will have a shorter diameter and larger leaf number [32].

K-means clustering
We computed an averaged connectivity matrix per participant group and subtracted the average matrix of healthy controls to obtain three difference matrices. These differences matrices were converted into three vectors, thereby only extracting those connections that were part of the NBS component. The three vectors were entered into a k-means clustering algorithm. The k-means clustering method is an unsupervised machine learning algorithm that groups data into a pre-specified number of k clusters. We determined the optimal number of according to Madhulatha and colleagues [34], by running the k-means clustering algorithm for k range of values (e.g. 1-10 clusters) and subsequently calculating the sum of squared errors (SSE) at each k. The most favorable balance between a low SSE and low k determines the optimum number of clusters [34]. In this study, the elbow method validated an optimum number of six clusters (Supplementary Figure 4).

Effects of age, sex, motion and modality
Due to the relatively large differences in age and sex between the groups, we tested for possible collinearity. To identify collinearity among variables the variance inflation factor (VIF) was calculated in SPSS 22.0 using linear regression. As a rule of thumb, a VIF > 3 warrants extra checking, and a VIF > 5 means that collinearity is highly likely. To further explore possible effects of age and sex as covariates, we replicated the Network Based Statistics analyses in smaller matched subgroups (without including age and sex as covariates). These smaller subgroups were matched using exact case-control matching in SPSS 20.0.
The residual relationship between subject movement and connectivity was evaluated according to two benchmarks; 1) assess residual relationship between motion and connectivity by estimating the residual QC-FC (quality control / functional connectivity) correlation; 2) estimating distance-dependent effects of motion on connectivity. The QC-FC relationship was computed by correlating the network edge weight of the 90-nodes with the mean relative displacement motion as calculated by MCFLIRT function in FSL [35]. To account for the participant's age and sex, a partial correlation was calculated. This was subsequently used to calculate the number of edges that significantly correlated with subject motion after using the false discovery rate (FDR) [36]. Secondly, distance dependent effects were estimated by correlating the Euclidean distance between the center of each node with the QC-FC correlation. Additional data on the residual relation between motion and functional connectivity is provided in Supplementary Figure 1.

Random forest classification algorithm
The random forest classifier has the advantage that cross-validation is done internally, meaning that no separate test set of participants is needed to estimate the generalized error of the training set [39]. The cross-validation is done internally by out-of-bag (OOB) error estimate. In each bootstrap training sample, one-third of the features are left out and not used in assembling the decision tree. As the random forest is built, each tree is tested on features not used in building that tree, termed OOB examples. Each feature that is left out is put down the tree to obtain a classification. At the end of the run, the model predicts the class that received most of the votes every time feature x was OOB, as well as the proportion of time the class is not equal to the original class of feature x. This is averaged over all features resulting in the OOB error estimate [39].
Each feature per classification receives a variable importance (VIMP) score between 0 and 1. This score is computed by the mean decrease in the Gini impurity criterion. This Gini impurity specifies how often a specific feature was selected for a split, and how large its overall discriminative value was for the classification. By adding up the Gini impurity decreases for each feature across all trees in the random forest, the feature importance (i.e. variable) is obtained.
The classification metrics sensitivity, specificity, and accuracy were determined as follows:

Subjects
Of all scans collected in the five studies, a subset of 749 participants was selected based on presence of a resting state and anatomical scan. Of these, n = 196 participants were excluded due to missing clinical data, poor imaging, data quality, processing errors, effects of motion, or brain abnormalities reported by the radiologist. See Supplementary  Table 1 for a detailed overview of excluded scans per study, and Supplementary Table 2 for an overview of excluded scans per clinical status (e.g. NC-H, SCZ, BD, CTRL). Next, we selected only those participants who experienced lifetime hallucinations, and bipolar patients without any lifetime history of psychotic experiences. We therefore excluded another 69 patients with bipolar-I disorder, and 6 patients with schizophrenia. Exclusion of these scans resulted in a final eligible pool of n = 483 subjects. On average, sex and age differed across the participant groups. Tukey post-hoc tests were conducted to provide more information on the exact differences (see Supplementary Table  3 for post-hoc tests on age, and Supplementary Table 4 for post-hoc tests on sex). Because age and sex differed across the groups, these variables were entered as covariates in the subsequent analyses.

Effects of age and sex
No significant relationship was present between age and sex was t = -0.055, P = 0.145. The collinearity diagnostics between age and sex was VIF = 1.000; tolerance = 1.000, indicating a low possibility of collinearity (VIF < 3). Age and sex were thus entered as covariates in the analysis. To further explore the effects of age and sex, we also replicated the NBS analyses in smaller matched subgroups, see Supplementary Fig 2. A similar pattern of connectivity disturbances was demonstrated for the non-clinical individuals and patients with schizophrenia. The analysis for patients with bipolar disorder-I was not significant, although a trend was observed (P = 0.097). This is probably due to a power problem as this analysis was done in a smaller matched subgroup. We therefore repeated this analysis at a more lenient significance level (P < 0.10) which yielded a similar pattern of connectivity disturbances as in the original analyses. The bipolar patients without hallucinations did not show any significant differences with healthy controls, neither at a more lenient significance threshold (P > 0.10).

Effects of motion
The residual relationship between motion and connectivity strength was estimated, revealing that none of the edges was significantly related to motion after FDR correction (P < 0.001, for all tests). The median QC-FC for healthy controls was 0.07; for patients with schizophrenia 0.06; for patients with bipolar disorder 0.05; and for non-clinical individuals 0.11. These numbers are in correspondence with previous findings on ICA-AROMA by Ciric and colleagues [35]. Distance dependent effects of motion are represented by the correlation displayed in each of the graphs in Supplementary Figure 1. The distance dependent effects. The results demonstrate that the distance dependent effect on motion was minimal for all groups.

Global connectome disturbances
The groups did not differ in functional connectivity strength (F (4) = 1.0; P = 0.434), meaning its effects on graph theoretical measures can be considered minimal [30]. None of the hallucination groups differed in terms of global network organization relative to healthy controls, neither in graph theoretical or MST measures (P > 0.05, FDR corrected).

Network Based Statistics
The NBS contrasts schizophrenia patients vs healthy controls, and non-clinical individuals vs healthy controls were both tested using a F-test at T = 13.0; P < 0.05 at 10,000 permutations. The contrast bipolar-I disorder patients with hallucinations vs healthy controls was tested at T = 8.5; p < 0.05 at 10,000 permutations. Differences between bipolar-I disorder patients vs healthy controls was tested at a range of T-thresholds, but did not result in any significant differences at either of the thresholds. According to Zalesky and colleagues [40], the choice of the T-statistic is somewhat arbitrary, but does not affect the specificity of the results. We therefore choose T-statistics that provided a roughly equal number of connections across the contrasts to facilitate visual comparison of results. All significant connections in the NBS network are displayed in Supplementary Table 6 for patients with bipolar disorder, Supplementary Table 7 for patients with schizophrenia, and Supplementary Table 8 for non-clinical individuals with hallucinations. The results of the overall group comparison of all hallucinating individuals compared to healthy controls are displayed in Supplementary  Table 9.

Conjunction analysis
To test for precise overlap of connections in each of the NBS networks, a conjunction analysis was applied. The connections in the NBS networks of the non-clinical individuals and schizophrenia patients were found to overlap with eight connections, while the non-clinical individuals and schizophrenia patients showed overlap with bipolar-I patients with hallucinations in only one connection (see Supplemental Figure 4 and Supplemental Table 5).

Supplementary Discussion
Implication for future research Connectivity disturbances in areas such as the cingulate gyrus, insula, hippocampal complex, precuneus and language areas suggests involvement of several large-scale networks, including the auditory, attention and task-control networks in the experience of hallucinations. This typical pattern could suggest that sensory areas may not be adequately mediated by higher-order cognitive regions [41]. These findings can be interpreted in light of the bottom-up/top-down explanatory model which proposes that hallucinations may arise from an imbalance between sensory input (bottomup) and top-down inhibitory control on sensory areas [42][43][44][45]. Disruption of bottom-up processing networks could lead to disturbed interpretation of incoming sensory input, causing patients to hallucinate.
For example, the circular inference theory explains the experience of hallucinations in terms of an imbalance between bottom-up and top-down information [46]. Reduced inhibitory control is hypothesized to lead to top-down and bottom-up information being reverberated (i.e. prior beliefs are misinterpreted as sensory observations) [46]. This may fit with the observed increased connectivity between the frontal and temporal area in both schizophrenia patients and non-clinical individuals in this study, as multiple overcounting of the same redundant sensory and prior information means that both areas will be more active at the same time (e.g. resulting in a higher correlation of BOLD signals) than to be expected in controls.
In support of the reduced inhibitory control theory, previous studies demonstrated that patients with schizophrenia have difficulties to suppress irrelevant information as established by various cognitive tasks measuring inhibitory control [47][48][49][50]. Similar deficits in inhibitory control were reported in non-clinical individuals with hallucinations [51-52]. Reduced inhibitory control was previously reported for bipolar patients in a manic episode [53], but to a lesser extent in euthymic states [54]. The bipolar-I disorder patients in the current study were all in euthymic phase at time of scanning. Euthymic bipolar patients have been found to perform intermediate to schizophrenia patients and controls on executive functioning tasks [55]. In all, these findings might point to mild inhibitory problems related to hallucinations, albeit to a lesser extent in bipolar-I disorder patients than schizophrenia patients and non-clinical individuals with hallucinations.
The present findings are correlational, and can therefore not confirm our hypothesis that connectivity between frontal, parietal and cingulate areas relates to inhibitory dysfunction. Future studies could investigate whether these disturbed connections point to inhibitory and excitatory dysfunction by studying global glutamate (Glu) receptor hypofunction (or the N-methyl-d-aspartate receptor (NMDA-R)) and gamma-aminobutyric acid (     Six sets of clusters were validated by the elbow method and entered into the k-means clustering algorithm. After a NBS analysis across the three hallucinating groups (NC-H, SCZ-H, BD-H), those edges that were found to significantly differ in functional connectivity across the three groups, were entered into the k-means clustering algorithm. A list of these edges can be found in Supplementary Table 9

Supplementary Table 6. Results of the network based statistics for non-clinical individuals with hallucinations versus healthy controls using the AAL atlas a .
a NBS F-test at T = 13.0; P < 0.05 with 10,000 permutations, adjusted for multiple comparisons. b Stats indicates the T-statistic of significant edges.