Identifying functional network changing patterns in individuals at clinical high-risk for psychosis and patients with early illness schizophrenia: A group ICA study.

Although individuals at clinical high risk (CHR) for psychosis exhibit a psychosis-risk syndrome involving attenuated forms of the positive symptoms typical of schizophrenia (SZ), it remains unclear whether their resting-state brain intrinsic functional networks (INs) show attenuated or qualitatively distinct patterns of functional dysconnectivity relative to SZ patients. Based on resting-state functional magnetic imaging data from 70 healthy controls (HCs), 53 CHR individuals (among which 41 subjects were antipsychotic medication-naive), and 58 early illness SZ (ESZ) patients (among which 53 patients took antipsychotic medication) within five years of illness onset, we estimated subject-specific INs using a novel group information guided independent component analysis (GIG-ICA) and investigated group differences in INs. We found that when compared to HCs, both CHR and ESZ groups showed significant differences, primarily in default mode, salience, auditory-related, visuospatial, sensory-motor, and parietal INs. Our findings suggest that widespread INs were diversely impacted. More than 25% of voxels in the identified significant discriminative regions (obtained using all 19 possible changing patterns excepting the no-difference pattern) from six of the 15 interrogated INs exhibited monotonically decreasing Z-scores (in INs) from the HC to CHR to ESZ, and the related regions included the left lingual gyrus of two vision-related networks, the right postcentral cortex of the visuospatial network, the left thalamus region of the salience network, the left calcarine region of the fronto-occipital network and fronto-parieto-occipital network. Compared to HCs and CHR individuals, ESZ patients showed both increasing and decreasing connectivity, mainly hypo-connectivity involving 15% of the altered voxels from four INs. The left supplementary motor area from the sensory-motor network and the right inferior occipital gyrus in the vision-related network showed a common abnormality in CHR and ESZ groups. Some brain regions also showed a CHR-unique alteration (primarily the CHR-increasing connectivity). In summary, CHR individuals generally showed intermediate connectivity between HCs and ESZ patients across multiple INs, suggesting that some dysconnectivity patterns evident in ESZ predate psychosis in attenuated form during the psychosis risk stage. Hence, these connectivity measures may serve as possible biomarkers to predict schizophrenia progression.

Introduction

Schizophrenia (SZ) is a severe and disabling mental disorder, characterized by positive symptoms (including hallucinations, delusions, and thought disorders), negative symptoms (including poor motivation, anhedonia, and social withdrawal), and cognitive impairments. Prior to the onset of psychosis, a prodromal period lasting from a few weeks to several years typically occurs for most SZ patients (Cannon, 2015). Individuals at clinical high-risk (CHR) for psychosis (Klosterkotter et al., 2001, Miller et al., 2003, Yung et al., 2005) exhibit a psychosis-risk syndrome principally defined by the presence of attenuated forms of the positive symptoms characteristic of SZ (McGlashan et al., 2010). Approximately 35% of individuals meeting CHR criteria transitions to psychosis within 2.5 years of initial ascertainment in longitudinal studies (Cannon et al., 2008, Fusar-Poli et al., 2012). The majority of CHR individuals who convert to psychosis develop a disorder in the SZ spectrum (Fusar-Poli et al., 2013, Woods et al., 2009). Yet, relatively few brain imaging studies directly compare CHR individuals with early illness SZ (ESZ) patients to determine whether functional brain abnormalities evident in ESZ are present in attenuated form, or in qualitatively distinct patterns, in the psychosis-risk syndrome. This knowledge gap is evident in the realm of brain intrinsic functional networks (INs) that reflect correlated activity among distributed brain regions. INs have shown promise as neuropsychiatric biomarkers for mental disorder diagnoses in general (Lee et al., 2013, Sporns, 2014), and specific patterns of IN dysconnectivity have been reported in SZ (Calhoun et al., 2009, Du et al., 2015b, Du et al., 2016b, Fitzsimmons et al., 2013, Sheffield and Barch, 2016). Although the heterogeneity of clinical outcomes for CHR individuals might lead to patterns of IN dysconnectivity that are qualitatively distinct from those observed in ESZ, the fact that the psychosis-risk syndrome is largely defined by attenuated forms of SZ positive symptoms leads us to hypothesize that CHR individuals may primarily show similar but attenuated IN dysconnectivity relative to ESZ patients. If true, this would suggest that the IN dysconnectivity patterns seen in SZ predate the onset of full-blown psychosis, albeit in a less severe manner. Moreover, given that most CHR individuals have had minimal or no prior exposure to antipsychotic medication, the presence of attenuated IN dysconnectivity in CHR individuals would provide strong evidence that the dysconnectivity findings in SZ patients are not secondary to the confounding effects of antipsychotic drugs. Accordingly, the current study sought to directly compare IN dysconnectivity patterns in CHR individuals and ESZ patients.

In order to extract INs from resting-state functional magnetic resonance imaging (fMRI) data, researchers have used one of two general analytic approaches: model-based and data-driven analyses (Li et al., 2009, van den Heuvel and Hulshoff Pol, 2010). Among model-based approaches, the most widely used is the region of interest (ROI)-based method (Biswal et al., 1995), which measures functional connectivity (correlation or coherence) among different ROIs or between a specified “seed” ROI and all other voxels in brain. Most previous IN studies of CHR individuals employed an ROI-based method seeded by thalamus (Anticevic et al., 2015), superior temporal gyrus (Yoon et al., 2015), posterior cingulate cortex (Shim et al., 2010), medial prefrontal cortex and anterior insula (Wotruba et al., 2014), Broca's area (Jung et al., 2012), or cerebellum (Wang et al., 2016). These studies computed correlations between the mean time series in a specific ROI and the time series of other voxels in brain. The identified IN abnormalities of CHR individuals were primarily located in the frontal lobes (Anticevic et al., 2015, Jung et al., 2012, Wang et al., 2016, Wotruba et al., 2014, Yoon et al., 2015), temporal lobes (Yoon et al., 2015), thalamus (Anticevic et al., 2015), sensory motor cortex (Anticevic et al., 2015) and Heschl's gyrus (Anticevic et al., 2015), some of which also showed an intermediate degree of abnormality in the CHR syndrome compared to first episode psychosis (Jung et al., 2012, Yoon et al., 2015). One limitation of such studies is that results depend on a somewhat arbitrary definition of a specific ROI, including delineation of its shape, extent, and precise location (Du et al., 2012, Liu, 2011). Furthermore, such studies typically only focus on several predefined ROIs and thus do not evaluate IN patterns across the whole brain.

In contrast to ROI-based method, data-driven approaches estimating INs do not require specification of predefined ROIs. These increasingly popular approaches include spatial independent component analysis (ICA) (Calhoun and Adali, 2012, Calhoun et al., 2001, Du et al., 2016a, Du and Fan, 2013), principle component analysis (PCA), and clustering methods (Du et al., 2014, van den Heuvel et al., 2008). In particular, ICA is a widely used approach that has shown great promise in identifying network-based biomarkers of psychiatric disorders such as SZ (Calhoun et al., 2011, Du et al., 2015b, Garrity et al., 2007, Khadka et al., 2013, Meda et al., 2014, Ongur et al., 2010). Spatial ICA on individual-subject's fMRI data decomposes an fMRI data-converted matrix (time points × voxels) as a linear combination of multiple independent components (ICs), of which meaningful ICs can be regarded as INs. Advantages of ICA, relative to ROI-based methods, include that it does not require selecting prior ROIs and that it can simultaneously estimate multiple INs from whole-brain data. Different from ROI-based method which represents each IN as functional connectivity between different ROIs, ICA regards each spatial IC as an IN where voxels with higher Z-scores tend to have higher intra-connectivity (or co-activation). For multi-subject studies, group ICA approaches (Calhoun and Adali, 2012) enable estimation of individual-subject components, some of which are identified as meaningful INs, based on group-level ICs. Traditional group ICA typically used either PCA-based or regression-based (e.g., dual regression) back-reconstruction (Calhoun et al., 2001, Erhardt et al., 2011). However, these methods cannot guarantee independence of individual-subject components. More recently, we have proposed a new back-reconstruction method called group information guided ICA (GIG-ICA) (Du et al., 2016a, Du and Fan, 2013, Du et al., 2015b). GIG-ICA utilizes a multi-objective function optimization algorithm to simultaneously optimize the independence among each subject's components as well as the correspondence between each group-level component and its associated subject-specific component, consequently resulting in more accurate networks.

In the present study, using resting-state fMRI data from CHR individuals, ESZ patients, and healthy controls (HCs), we estimate whole-brain multiple INs for each subject using the GIG-ICA method and then investigate all possible network changing patterns in each IN. We hypothesize that ESZ patients would show dysconnectivity (either hypo- or hyperconnectivity) across multiple INs relative to HC individuals. We further hypothesize that CHR individuals would generally show similar, but attenuated, dysconnectivity across many INs in which ESZ patients show dysconnectivity, consistent with their attenuated positive symptoms relative to ESZ patients. However, given that dysconnectivity in some INs may emerge after the onset of schizophrenia, we expect to find some INs that showed dysconnectivity in ESZ patients but normal connectivity in CHR individuals. Our overarching goal is to perform a comprehensive whole-brain investigation of resting-state connectivity across multiple INs in order to identify and directly compare the resting-state dysconnectivity patterns in ESZ and in the psychosis-risk syndrome.

Materials and methods

Subjects

We analyzed resting-state fMRI data from 70 HCs, 53 CHR individuals, and 58 ESZ patients. Subject demographic and clinical information are presented in Table 1. There were no significant differences among the three groups in age or gender (see Table 1). CHR individuals were recruited from a psychosis-risk research clinic at the University of California, San Francisco (UCSF). Individuals recruited to the CHR group met the Criteria of Prodromal Syndromes (COPS) based on a Structured Interview for Prodromal Syndromes (SIPS) (Miller et al., 2003) administered by trained clinicians. The COPS comprises three non-mutually exclusive syndromes (McGlashan et al., 2010): Attenuated Positive Symptom Syndrome (APSS), Brief Intermittent Psychotic Syndrome (BIPS), and/or Genetic Risk and Deterioration Syndrome (GRD). In our study, the majority (49/53, 92.5%) of CHR subjects met COPS criteria for APSS. Supplementary Table S1 shows the number and percentage of CHR subjects with each syndrome. Ratings of symptom severity in CHR individuals were obtained using the Scale of Prodromal Symptoms (SOPS) (McGlashan et al., 2010), an embedded scale within the SIPS. Most of the CHR subjects (41/53; 77%) were antipsychotic medication-naive at the time of scanning. ESZ patients within five (mean ± standard deviation = 2.08 ± 1.37) years of illness onset were recruited from an early psychosis clinic at UCSF and from community clinics. Diagnosis of schizophrenia or schizoaffective disorder in ESZ subjects was verified using the Structured Clinical Interview for DSM-IV (SCID) (First et al., 1995) and symptom severity was assessed using the Positive and Negative Syndrome Scale (PANSS) (Kay et al., 1987). Most of the ESZ patients (53/58; 91%) were taking antipsychotic medication at the time of testing. HC subjects were recruited from the community and did not meet current or lifetime criteria for major Axis I psychiatric disorders based on the SCID. The same data were analyzed using dynamic connectivity analysis in our previous study (Du et al., 2017a).

Table 1

Subject demographic and clinical characteristics.

	HCs(n = 70)		CHR individuals(n = 53)		ESZ patients(n = 58)
	Mean	SD	Mean	SD	Mean	SD
Age (years)	21.9	5.6	20.4	4.5	21.8	3.8
PANSS positive symptoms	–	–	–	–	13.7	4.8
PANSS negative symptoms	–	–	–	–	17.4	6.6
PANSS general symptoms	–	–	–	–	32.9	8.9
SOPS positive symptoms	–	–	9.4	4.5	–	–
SOPS negative symptoms	–	–	12.2	5.8	–	–
SOPS general symptoms	–	–	8.1	4.5	–	–
SOPS disorganization symptoms	–	–	5.4	3.4	–	–
Maximum translation motion displacement (mm)	0.8	0.6	1.1	1.4	1.2	1.3
Maximum rotation motion displacement (degree)	0.8	1.0	1.0	0.9	0.9	0.8
	n	%	n	%	n	%
Male	41	59%	32	62%	38	65%
Subjects taking antipsychotic medication	–	–	12	23%	53	91%
Names of primary antipsychotic drugs	–		Abilify, Seroquel, Risperdal, Zyprexa, Risperdal		Abilify, Seroquel, Clozapine, Risperdal, Zyprexa, Risperdal
Name of primary antidepressant drugs	–		Prozac		Prozac

SD, standard deviation; HCs, healthy controls; CHR, clinical high-risk; ESZ, early illness schizophrenia; PANSS, Positive and Negative Syndrome Scale; SOPS, Scale of Prodromal Symptoms. p-value = 0.7 for gender examined by Chi Square test; p-value = 0.2 for age examined by analysis of variance. The maximum translation motion displacement was computed as the maximum translation across all axes (x-axis, y-axis and z-axis) and the whole scanning. The maximum rotation motion displacement was computed as the maximum rotation across the pitch, roll and yaw and the whole scanning.

MRI scan acquisition

All brain images were acquired on a 3 T Siemens TIM Trio scanner at the UCSF Neuroimaging Center. Resting-state scans were acquired using whole-brain echo-planar imaging (EPI) sequences. Rest scans lasted 6 min, during which 180 functional images were obtained (32 axial slices, 3.5 mm slice thickness, 1.05 mm inter-slice gap, TR = 2 s, TE = 29 ms, flip angle = 75°, FOV = 24 cm, 64 × 64 matrix). Subjects were instructed to rest with their eyes closed and to stay awake.

MRI data preprocessing

FMRI data from each subject were preprocessed using Statistical Parametric Mapping (SPM8) (http://www.fil.ion.ucl.ac.uk/spm). The first ten volumes were discarded to allow for T1 equilibration effects, and the remaining volumes in the run were then slice-time corrected and realigned to the first volume to correct for head motion. In terms of head motion, translations were < 4 mm and rotations did not exceed 4° in all axes over the entire fMRI run for all included subjects. There were no significant group differences (p > 0.05, tested by analysis of variance) in the motion parameters displayed in Table 1. Subsequently, the images were spatially normalized to the Montreal Neurological Institute (MNI) EPI template, resliced to 3 mm × 3 mm × 3 mm voxels, and smoothed with a Gaussian kernel with a full-width at half-maximum (FWHM) of 6 mm.

Estimating intrinsic functional networks using GIG-ICA

We applied GIG-ICA (Du et al., 2016a, Du and Fan, 2013, Du et al., 2017b, Du et al., 2015b) to the preprocessed fMRI data of all subjects to estimate each subject's INs, resulting in INs with correspondence across subjects while also retaining subject-specific network features. Our method primarily involved the following three steps. First, based on the temporally-concatenated fMRI data of 181 subjects, ICA with the Infomax algorithm (Bell and Sejnowski, 1995) was used to estimate the group-level ICs. In order to decrease the influence of random initializations, we performed ICASSO (Himberg et al., 2004, Ma et al., 2011) with 20 ICA runs followed by selection of the most reliable ICA run to obtain reliable group-level ICs. Before the group-level ICA, subject-level and group-level PCAs (Calhoun et al., 2001, Erhardt et al., 2011) were applied to reduce the dimensionality of the fMRI data. In this study, the number of ICs was set to 30 based on previous work (Abou-Elseoud et al., 2010, Du et al., 2015b), and the number of principle components used in the subject-level PCAs was specified as a greater number (i.e., 60) as recommended (Erhardt et al., 2011). Thus, 30 group -level ICs were obtained and Z-scored. It is worth noting that the sign of each resulting group-level IC is arbitrary. In our study, we computed the skewness of each IC and flipped the IC if its skewness was negative. Consequently, the positive Z-scores in each IC represented the voxels contributing the most to the associated network. Second, we manually identified the meaningful group-level networks after carefully inspecting all group-level ICs' spatial maps. If the primary activation of one IC's spatial map is located in gray matter, smooth and continuous (Allen et al., 2014, Du et al., 2016a), we preferred to regard it as an interesting functional network. Also, we referred to networks reported in previous studies (Du and Fan, 2013, Du et al., 2015b, Smith et al., 2009, Zuo et al., 2010) to further confirm. The removed noise-related group-level ICs reflected head motion, physiological noise, and influence of scanner. Third, the remaining meaningful group-level ICs were characterized as the group-level INs and then used to calculate the subject-specific INs based on a multi-objective function optimization algorithm (Du and Fan, 2013). GIG-ICA automatically yields Z-scored the subject-specific INs. Finally, a regression model (Du et al., 2016a) was utilized to estimate the corresponding time courses (reflecting fluctuations) of the individual INs for each subject's fMRI data.

It is known that for each voxel in a given individual IN, a positive Z-score indicates that the time series (after preprocessing) in the voxel has a positive correlation with the time course of the IN; in contrast a negative Z-score represents that the time series in the voxel has a negative correlation with the time course of the IN. Therefore, for one specific IN, a region containing voxels with high positive Z-scores (taking a region only including two voxels as an example, Z-scores of two voxels = 10 and 9) tends to have stronger co-activation (i.e., positive correlation) within the IN compared to a region containing voxels with low positive Z-scores (e.g., Z-scores of two voxels = 2 and 1). In such case, the mean of the Z-scores in one region can reflect the co-activation extent of the region. In fact, if a region contains voxels with greater negative Z-scores (e.g., Z-scores of two voxels = − 10 and − 9) relative to the other region (e.g., Z-scores of two voxels = − 2 and − 1), the former also tends to have stronger co-activation than the later. In this study, considering that the positive Z-scores contributing the most include important spatial maps in each IN, we focused on investigating group differences in voxels primarily associated with positive Z-scores to simplify the subsequent analyses. It is also worth pointing out that if one region includes both higher and lower positive Z-scores (e.g., Z-scores of two voxels = 10 and 1) compared to the other region (e.g., Z-scores of two voxels = 2 and 8), we cannot determine which region has greater co-activation using Z-scores. Therefore, in this work, we applied two-sample t-tests (rather than analysis of variance) to identify brain regions including voxels with consistent changing patterns between groups.

Identifying significant discriminative regions in intrinsic functional networks

Based on the estimated subject-specific INs, we investigated differences in each corresponding IN across the three groups in order to examine the network abnormalities of the ESZ and CHR group relative to each other and to the HC group. For each IN, we performed three analysis steps that include (1) extracting important voxels positively contributing to the network (i.e., voxels primarily involving positive Z-scores) and generating a network mask, (2) identifying regions based on all possible 19 inter-group changing patterns in the network mask, and (3) investigating group differences in the co-activation of each significant discriminative region. We describe the detailed processing on each IC as follows.

In the first step, for each voxel, a right-tailed one-sample t-test (p < 0.01 with Bonferroni correction, i.e., p < 0.01/the number of voxels in a brain mask = 0.01/69078 = 1.45e − 7) was applied to the voxel's Z-scores from the corresponding subject-specific INs of all 181 subjects to identify voxels whose Z-scores were positive in most of the subjects, resulting a network mask. The subsequent second and third steps were performed on the voxels in the network mask. Additionally, we tested if our strategy may fail to find group differences in some important voxels (or regions) showing positive Z-scores in one group (e.g., HC) but not showing positive Z-scores in other groups (e.g., CHR and ESZ). So, we examined whether performing one-sample t-tests on each group separately may generate significantly different network masks compared to performing one-sample t-tests on all the three groups. For each group (e.g., HC), voxel-wise right-tailed one-sample t-tests were applied on the corresponding individual-subject's networks. The parameter (p < 0.01 with Bonferroni correction, i.e., p < 0.01/the number of voxels in brain mask = 0.01/69078 = 1.45e − 7) used in one-sample t-tests on each group's subjects was as same as the parameter of one-sample t-tests on all subjects of three groups. We then combined the three masks obtained separately from the three groups in order to compare them with the mask obtained using the whole dataset.

In the second step, we defined 19 non-overlapping patterns (see Fig. 1) that reflect all possible significant differences among the three groups based on statistical analysis results from two-tailed two-sample t-tests, and then identified regions including voxels coinciding with each pattern. Specifically, for each significant voxel in the network mask identified, we performed two-tailed two-sample t-tests (p < 0.05) on its Z-score for three pairs of groups (HC vs. CHR, CHR vs. ESZ, and HC vs. ESZ) and then examined whether its Z-score exhibited one of 19 characteristic patterns. Among the 19 patterns, two patterns (Fig. 1A–B) corresponded to strictly decreasing and strictly increasing Z-scores from HC to CHR to ESZ. For example, for the strictly decreasing pattern, t-value > 0 was required for both the HC vs. CHR test and the CHR vs. ESZ test, reflecting significant decreases in Z-scores from the healthy condition to the psychosis-risk state and from the risk-state to full-blown schizophrenia. Four patterns (Fig. 1C–F) involved situations where the ESZ or HC group was significantly lower or significantly higher than the other two groups. For example, in the ESZ-decreasing pattern, t-value > 0 was required for both the HC vs. ESZ and CHR vs. ESZ comparisons, with no significant difference between the HC and CHR groups. Six additional patterns (Fig. 2G–L) were CHR-unique alteration related patterns in which the CHR group had significantly lower or higher Z-scores than the other two groups. For example, in terms of the CHR-decreasing pattern 1 (Fig. 2G), we required t-value > 0 for the HC vs. CHR comparison as well as t-value < 0 for the CHR vs. ESZ comparison, with no significant difference between the ESZ and HC groups. There were also six patterns (Fig. 1M–R) in which only one pair of groups showed a group difference. Taking the HC > ESZ pattern for an example, t-value > 0 for the HC vs. ESZ comparison and no group differences for other comparisons were required. Finally, a no-difference pattern (Fig. 2S) was possible, where all three groups showed no significant differences (i.e., p > 0.05 for all comparisons). Regarding each pattern, the associated voxels were identified and then corrected (p < 0.05) for multiple comparisons using Monte Carlo simulation (N = 1000) (Ledberg et al., 1998) to obtain regions each of which were spatially continuous. Thus, each remaining region was required to have a size greater than a given voxel number that was determined by the threshold p < 0.05 and the network mask from the voxel-wise one-sample t-tests using 1000 Monte Carlo simulations. Except for the last pattern representing no significant difference across groups, we named the identified regions as significant discriminative regions (SDRs) within the associated INs. In our work, we applied two-tailed two-sample t-tests to define the changing patterns because we cannot make a hypothesis that one group (e.g., HC) shows higher (or lower) Z-scores in network than the other group (e.g., CHR). We are testing for the possibility of the relationship in both directions (i.e., the possibility of one group shows higher and lower Z-scores than the other group).

Fig. 1

The 19 possible non-overlapping changing patterns defined for each voxel within each intrinsic functional network (IN) based on the results of two-tailed two-sample t-tests on Z-scores (p < 0.05) for three pairs of groups (i.e., HC vs. CHR, HC vs. ESZ, and CHR vs. ESZ comparisons). The requirements used to identify each pattern are displayed. Group difference in Z-scores identified by a two-tailed two-sample t-test between a pair of groups (p < 0.05) is denoted by a red line.

Fig. 2

The whole analysis framework. The steps include A: estimating subject-specific networks using GIG-ICA, B: performing voxel-wise right-tailed one-sample t-tests to identify a mask for each network, C: performing two-tailed two-sample t-tests for three pairs (HC vs. CHR, HC vs. ESZ, and CHR vs. ESZ) on each voxel within the network mask to identify significant discriminative region (SDR) with the same changing pattern, and D: investigating group difference in the co-activation of each SDR.

In the third step, for each SDR including voxels with the same pattern in a specific IN, we tested whether its co-activation differed significantly among the HC, CHR and ESZ groups. Based on the results from the above mentioned voxel-wise two-sample t-tests of different voxels in one SDR, we applied a Fisher's combined probability test (Fisher, 1925) to measure the group difference of each SDR's co-activation, resulting in a combined p-value for any comparison between two groups (e.g., HC vs. CHR). Since Fisher's combined probability test was proposed originally under an assumption of independence among separate tests, we also performed a non-parametric combination method including 10,000 permutations as a supplement. The non-parametric combination procedure works even when independence is untenable (Winkler et al., 2016). The detailed steps for the non-parametric combination method are described in the Supplementary materials. If the combined p-value of one comparison between two groups was smaller than 0.05/3, we regarded the co-activation in SDR significantly different between the two groups. Additionally, since the mean Z-score within voxels of each SDR can be used to reflect co-activation strength in the SDR, we calculated the Pearson correlation between the co-activation strengths of each SDR within one specific IN and the symptom severity ratings (in Table 1) for the CHR group and the ESZ group to explore the association between network measures and symptoms. The significance level was set to p < 0.05 for the correlation analyses.

The whole analysis framework is shown in Fig. 2.

We also tested if medication has significant effects on the identified measures showing group differences. First, we converted all anti-psychotic data to their respective chlorpromazine (CPZ) dosage equivalents (Danivas and Venkatasubramanian, 2013, Woods, 2003) for ESZ patients with available dose-level medication data. Only one ESZ patient had no CPZ data. We then used a non-parametric Spearman's rank correlation method to evaluate the association between the CPZ equivalent and the co-activation (reflected by the mean Z-score) of each SDR (p < 0.05 with Bonferroni correction, i.e., p-value threshold = 0.05/the number of all SDRs) for ESZ patients. Spearman's rank correlation analysis was applied due to the non-normality of the CPZ variable. For CHR subjects, we did not perform correlation analysis on the CPZ equivalents due to the small number of the medicated CHR subjects and the relatively noisy property of CPZ equivalents. Instead, we examined difference in the co-activation of each SDR between the CHR individuals without medication treatment and the CHR individuals taking medication using a two-tailed two-sample t-test (p < 0.05 with Bonferroni correction, i.e., p-value threshold = 0.05/the number of all SDRs).

Results

Brain intrinsic functional networks

Thirty reliable group-level ICs were estimated from the fMRI data of 181 subjects. Fifteen group-level ICs were then identified as noise-related components and removed (see Supplementary Fig. S1), leaving 15 group-level INs that were employed to guide the computation of subject-specific INs. The remaining 15 group-level INs are shown in the Supplementary Fig. S2. It is seen that for those networks, the spatial maps with positive Z-scores provided more important information than the spatial maps with negative Z-scores due to that the skewness of each IN was positive. As stated in the method section, for each group-level IN, the corresponding subject-specific IN was estimated for each subject using GIG-ICA. Regarding each of the 15 INs, we show the t-value map obtained by performing voxel-wise right-tailed one-sample t-tests on all subjects' individual networks in Fig. 3. These 15 INs primarily consisted of two default mode networks (DMNs) (IN 5 and IN 10), two fronto-parieto-cerebellar networks (IN 2 and IN 9), three vision-related networks (IN 1, IN 3 and IN 4), a sensory-motor network (IN 12), a salience network (IN 11), an auditory-related network (IN 6), a visuospatial network (IN 7), an anterior insula network (IN 8), a fronto-occipital network (IN 13), a fronto-parietal network (IN 14), and a fronto-parieto-occipital network (IN 15).

Fig. 3

A: t-value map of each of the 15 intrinsic functional networks (INs). Each t-value map was obtained by performing voxel-wise right-tailed one-sample t-tests (p < 0.01 with Bonferroni correction) on the corresponding subject-specific INs from the 181 subjects. B: t-value maps of all 15 INs, shown together. Different INs are shown using different colors.

As stated in the method section, for each IN, the voxels passing multiple comparison correction of the one-sample t-tests on all subjects' individual networks were included in its network mask. It is worth noting that the brain regions (shown in the Section 3.2) identified using the specific changing patterns in each IN located within the relevant network mask. In the Supplementary Fig. S3, we show both the mask obtained from the one-sample t-tests on all groups' individual networks and the combined mask obtained from the one-sample t-tests on each group's individual networks. For all INs, the two kinds of mask were very similar while the slight difference was only on the edge of maps. That indicates if voxels showed significant positive Z-scores in one group, they also tended to show positive Z-scores in other groups. The mask we used had slightly more voxels than the union mask of three separate groups. Therefore, our analyses did not miss any important voxels with positive Z-scores.

Significant discriminative regions in intrinsic functional networks

Given the specified 19 changing patterns shown in Fig. 1, the associated voxels in each IN were identified (see Supplementary Table S2). We found that 12 voxels in the salience network (IN 11) exhibited strictly decreasing trends, and only two voxels in the visuospatial network (IN 7) showed strictly increasing trends, in Z-scores from the HC group to the CHR group to the ESZ group. Few or no voxels belonged to the CHR-decreasing pattern 2 (and 3) and the CHR-increasing pattern 2 (and 3). Excepting for the no-difference pattern that most of the voxels were assigned to, the mostly occupied pattern was the HC > ESZ pattern. After multiple comparison correction, 24 SDRs were identified to be relevant to the 18 discriminative changing patterns (excluding the no-difference pattern). Detailed information of the 24 SDRs can be found in Table 2 and Supplementary Table S3.

Table 2

The 24 significant discriminative regions (SDRs) extracted from the 15 intrinsic networks (INs). For each SDR, the associated network, its index (ID), the voxels' Z-score changing pattern across groups, the number of included voxels, the volume in cubic millimeters, and the related brain regions are listed.

IN ID and IN name	ID of SDR identified in INs	Voxels' Z-score changing pattern	Number of voxels	Region volume (mm3)	Brain region name
IN 1: Vision-related network	SDR-1-1	HC > ESZ	104	2808	Lingual gyrus (L)
IN 3: Vision-related network	SDR-3-1	HC < ESZ	80	2160	Inferior occipital gyrus (R)
IN 3: Vision-related network	SDR-3-2	HC > CHR	116	3132	Lingual gyrus (R)
IN 4: Vision-related network	SDR-4-1	HC > ESZ	125	3375	Lingual gyrus (L)
IN 4: Vision-related network	SDR-4-2	HC < CHR	120	3240	Fusiform gyrus (R)
IN 5: Default mode network	SDR-5-1	CHR > ESZ	140	3780	Precuneus and posterior cingulate cortex (R)
IN 6: Auditory-related network	SDR-6-1	ESZ-decreasing	112	3024	Heschl's gyrus (L)
IN 6: Auditory-related network	SDR-6-2	CHR > ESZ	185	4995	Insula and rolandic operculum (L)
IN 7: Visuospatial network	SDR-7-1	ESZ-decreasing	74	1998	Precentral gyrus (R)
IN 7: Visuospatial network	SDR-7-2	HC > ESZ	106	2862	Postcentral gyrus (R)
IN 8: Anterior insula network	SDR-8-1	HC < CHR	118	3186	Superior temporal gyrus (L)
IN 9: Fronto-parieto-cerebellar network	SDR-9-1	HC < ESZ	98	2646	Superior frontal gyrus (R)
IN 10: Default mode network	SDR-10-1	CHR > ESZ	154	4158	Anterior cingulate gyrus (R)
IN 10: Default mode network	SDR-10-2	CHR < ESZ	127	3429	Superior frontal gyrus, medial (L)
IN 11: Salience network	SDR-11-1	ESZ-decreasing	171	4617	Middle cingulum (L)
IN 11: Salience network	SDR-11-2	ESZ-increasing	77	2079	Middle frontal gyrus (R)
IN 11: Salience network	SDR-11-3	HC > ESZ	161	4347	Thalamus (L)
IN 11: Salience network	SDR-11-4	CHR > ESZ	236	6372	Middle cingulum (L)
IN 12: Sensory-motor network	SDR-12-1	ESZ-decreasing	81	2187	Paracentral lobule (L)
IN 12: Sensory-motor network	SDR-12-2	HC-decreasing	135	3645	Supplementary motor area (L)
IN 12: Sensory-motor network	SDR-12-3	HC < CHR	82	2214	Precentral and supplementary motor area (L)
IN 13: Fronto-occipital network	SDR-13-1	HC > ESZ	125	3375	Calcarine (L)
IN 14: Fronto-parietal network	SDR-14-1	HC < CHR	142	3834	Precuneus (R)
IN 15: Fronto-parieto-occipital network	SDR-15-1	HC > ESZ	142	3834	Calcarine (L)

Note: L, left; R, right. Regarding the SDR's ID, SDR-M-N means the Nth SDR in the Mth IN.

Fig. 4A shows the identified 24 SDRs and the results obtained from comparing the co-activation of each SDR among groups. We also report the relevant co-activation comparison results using the non-parameter combination method in the Supplementary Fig. S4. It is shown that the results in Figs. 4 and S4 were similar in terms of the combined p-values and consistent after the multiple comparison correction. Regarding the ESZ-decreasing pattern, there were four SDRs identified, which included left Heschl's gyrus from the auditory-related network, right precentral gyrus from the visuospatial network, left median cingulate of the salience network, and left paracentral lobule from the sensory-motor network. Evaluated by Fisher's combined probability test, we found that the co-activation in each of the four SDRs was significantly reduced in the ESZ group compared to both the HC and CHR groups, suggesting that ESZ patients had greater hypo-connectivity than CHR individuals in these SDRs. Regarding the ESZ-increasing pattern, we found one SDR, i.e., right middle frontal gyrus from the salience network, that showed relatively higher co-activation in the ESZ population than the other groups. In addition, the HC-decreasing pattern was evident in an SDR comprising left supplementary motor area from the sensory-motor network that persisted after correction for multiple comparisons. Its co-activation had enhanced strength for both CHR and ESZ groups, relative to the HC group.

Fig. 4

A: The significant discriminative regions (SDRs) identified from the 15 intrinsic functional networks (INs) and the co-activation (measured by the mean Z-score) of each SDR in the associated IN for healthy control (HC), clinical high-risk (CHR), and early illness schizophrenia (ESZ) groups, separately. SDR-M-N means the Nth SDR in the Mth IN. Table 2 includes more detailed information about these SDRs. SDRs identified using different changing patterns are shown using different colors, and the associated INs are shown in yellow color. The Fisher's combined p-values corresponding to all comparisons between any pair of groups (HC vs. CHR, HC vs. ESZ, and CHR vs. ESZ) are shown in each subfigure's title. Each SDR's mean Z-scores of all subjects in one group are shown using a boxplot. In each boxplot, the central line is the median; the square is the mean; and the edges of the box are the 25th and 75th percentiles. The whiskers extend to 1 inter-quartile range, and each outlier is displayed with a “+” sign. B: Significant associations between the clinical symptom scores (in Table 1) and the co-activation of SDRs.

Numerous voxels in the INs indicated progressive Z-scores from the healthy condition to psychosis risk to schizophrenia. Regarding the HC > ESZ pattern, six SDRs from different networks were found. These SDRs involved the left lingual gyrus of two vision-related networks, the right postcentral cortex of the visuospatial network, the left thalamus region of the salience network, left calcarine region of the fronto-occipital network and fronto-parieto-occipital network. Interestingly, each of the six SDRs tended to display a strictly decreasing change in its co-activation from HCs to CHR individuals to ESZ patients, assessed using the Fisher's combined probability test. In terms of the HC < ESZ pattern, two SDRs were extracted. Evaluated by combined p-values, the right inferior occipital gyrus in the vision-related network showed a common increased co-activation in ESZ and CHR, while the right superior frontal gyrus in the fronto-parieto-cerebellar network had a progressive increasing alteration from HC to CHR and from CHR to ESZ.

In our study, we also found that using the HC > CHR pattern, the HC < CHR pattern, the CHR > ESZ pattern and the CHR < ESZ pattern, 10 SDRs were identified, reflecting that the co-activation in these regions of the CHR group was not intermediate between the HC group and the ESZ group. Among the 10 SDRs, 8 SDRs in the CHR individuals showed enhanced co-activation compared to the other two groups. Precuneus, posterior cingulate cortex and anterior cingulate gyrus in the DMNs showed increased co-activation in the CHR group compared to the other two groups, indicating the complexity of CHR abnormality.

In summary, widespread INs were diversely impacted. Among the voxels in the identified significant discriminative regions (obtained using all 19 possible changing patterns excepting the no-difference pattern), more than 25% of the voxels from six INs were identified with the HC > ESZ pattern. In addition, the CHR > ESZ pattern was found in 24% of the significant voxels from four INs, and the ESZ-decreasing pattern was found in 15% of the significant voxels from four INs. Our findings thus suggest that many brain regions of CHR individuals showed an intermediate connectivity trend between the HC and ESZ groups. Compared to the HC individuals, the reduced connectivity was the primary type of dysconnectivity evident in ESZ individuals. Some brain regions also showed a CHR-unique connectivity (primarily CHR-increasing alterations).

Furthermore, significant associations were identified between SDRs' co-activation strengths and the symptom scores. As shown in Fig. 4B, the co-activation of superior temporal gyrus region in the anterior insula network was positively correlated with SOPS Disorganization scores in the CHR group (r = 0.42, p-value = 0.01). The co-activation of left supplementary motor cortex in the sensory-motor network was negatively correlated with both PANSS Positive Symptom scores (r = − 0.32; p-value = 0.02) and Negative Symptom scores (r = − 0.27; p-value = 0.04) in the ESZ group; and the co-activation of left calcarine region of the fronto-parieto-occipital network was negatively correlated with SOPS Disorganization scores in the CHR group (r = − 0.37; p-value = 0.02). Strikingly, each of the three SDRs with negative correlations with symptoms demonstrated a decreasing co-activation trend from the HC to the CHR to the ESZ groups (Fig. 4A). Regarding the ESZ group, the CPZ dosage equivalent did not have significant association with the co-activation (measured by the mean Z-score) in any of the identified SDRs showing a main effect of group differences. The details can be found in Supplementary Table S4. For the CHR group, our results (shown in Supplementary Table S5) suggest that most (all but one) of the SDRs showed no significant differences in the co-activation between the medicated and unmedicated CHR subgroups.

Discussion and conclusions

Individuals who meet clinical criteria for the psychosis-risk syndrome are neurobiologically vulnerable to developing a full-blown psychotic disorder, particularly schizophrenia (Cannon et al., 2008, Fusar-Poli et al., 2013, Fusar-Poli et al., 2012, Woods et al., 2009). Abnormalities in several ICA-derived INs have been reported in schizophrenia patients (Calhoun et al., 2009, Du et al., 2015b, Mattiaccio et al., 2016, Ongur et al., 2010). In this study, our purpose was to explore which brain regions in which networks were influenced by psychosis risk; further, we sought to determine if CHR individuals have similar network impairments as seen in SZ patients early in their illness course.

Using the GIG-ICA method, our results revealed multiple altered INs in the ESZ and CHR groups, which primarily included the DMN, salience, auditory-related, visuospatial, sensory-motor and parietal networks. Most of the identified INs have been previously associated with cognitive, emotional, executive, sensory and motor functions (Bassett and Bullmore, 2009, Bressler and Menon, 2010). Regarding different INs, various brain regions exhibited abnormally reduced or enhanced co-activation in both the ESZ and CHR groups. In total, 24 brain regions showing group differences were identified from all INs.

Compared to HCs, ESZ patients exhibited diminished co-activation strength in the majority (17) of these regions, although some regions with increased co-activation strength in the ESZ patients were also found. Our findings generally support prior studies that have predominantly shown reduced connectivity in SZ (Pettersson-Yeo et al., 2011, Zhou et al., 2015), but are also in agreement with some previous work reporting increased connectivity in SZ (Whitfield-Gabrieli et al., 2009, Zhou et al., 2007).

There were some brain regions that showed a significant step-wise worsening of dysconnectivity from HC to CHR to ESZ groups. In particular, the left lingual gyrus of two vision-related networks, the right postcentral cortex of the visuospatial network, the left thalamus region of the salience network, left calcarine region of the fronto-occipital network and fronto-parieto-occipital network displayed a strictly decreasing change in its co-activation from HCs to CHR individuals to ESZ patients, while the right superior frontal gyrus in the fronto-parieto-cerebellar network had a progressive increasing alteration from HC to CHR and from CHR to ESZ. Previous work has also found group differences in thalamic connectivity (Anticevic et al., 2015, Seiferth et al., 2008). While this pattern of results is consistent with the possibility of progressive dysconnectivity as individuals develop the psychosis-risk syndrome and then transition to full-blown psychosis, this pattern could also arise if only the small subgroup of CHR individuals who eventually transition to psychosis exhibit abnormal dysconnectivity on par with the ESZ patients. Unfortunately, sufficient clinical follow-up data were not available to help distinguish between these possibilities.

Compared to both HC and CHR individuals, ESZ patients demonstrated hypo-connectivity mainly in left Heschl's gyrus from the auditory-related network, right precentral gyrus from the visuospatial network, left median cingulate of the salience network, and left paracentral lobule from the sensory-motor network, but showed stronger co-activation in the right middle frontal gyrus from the salience network. Prior schizophrenia studies have similarly shown abnormalities in left Heschl's gyrus (Shinn et al., 2013) associated with auditory hallucinations (Javitt and Sweet, 2015). Our current results indicate that these changes are not prominent in the psychosis-risk syndrome, although it remains a possibility that they would be evident in the subset of CHR individuals who eventually progress to full psychosis.

There were additional brain regions that exhibited similarly abnormal co-activation in CHR and ESZ groups, consistent with the emergence of dysconnectivity in these regions prior to the onset of psychosis and the possibility that they reflect the underlying vulnerability for schizophrenia rather than a mechanism involved in the pathogenesis of the full-blown clinical syndrome. Both CHR and ESZ subjects showed increased co-activation in the left supplementary motor area of the sensory-motor network and the right inferior occipital gyrus of the vision-related network, compared to HCs. Our findings are consistent with prior work showing abnormal connectivity in CHR and SZ groups in the sensory motor cortex (Anticevic et al., 2015) and occipital gyrus (Fryer et al., 2016).

It is worth noting that CHR-unique alterations (primarily the CHR-increasing connectivity) were also evident in some brain regions, reflecting that the co-activation in these regions of the CHR group was not intermediate between the HC group and the ESZ group. Our previous study (Du et al., 2017a) using dynamic connectivity analysis on the same data also found CHR-specific alterations relative to HCs and ESZ patients. A Meta-analysis (Fusar-Poli et al., 2013) reported that although the majority (73%) of CHR individuals who convert to psychosis develops schizophrenic psychoses (including schizophrenia, schizophreniform disorder and schizoaffective disorder), 11% of the converted CHR individuals develops affective psychoses (including depression with psychotic features and bipolar disorder with psychotic features). In our study, ESZ patients included patients with schizophrenia or schizoaffective disorder. So, these CHR-unique alterations that were not present in ESZ patients may reflect brain function alterations that are closer to abnormalities in affective psychoses (e.g., bipolar disorder with psychotic symptoms). Our previous studies (Du et al., 2017c, Du et al., 2015a, Du et al., 2015b) have found that there are significant differences in functional connectivity among schizophrenia, schizoaffective disorder and psychotic bipolar disorder. Unfortunately, in the study, we cannot verify this point due to lacking of data from patients with affective psychoses.

In addition to the results of the group comparisons, clinical symptom correlations were also found for some of the regions. In ESZ patients, the co-activation of left supplementary motor cortex in the sensory-motor network was negatively correlated with both PANSS Positive and Negative Symptom scores. In CHR individuals, the co-activation of the left calcarine sulcus in the fronto-parieto-occipital network was negatively correlated with SOPS Disorganization Symptom scores. Indeed, the co-activation of these regions also showed a decreasing trend from HCs to CHR individuals to ESZ patients, supporting the possible role of these regions as biomarkers of clinical severity. We also found that there was no significant association between the CPZ dose equivalents of the ESZ patients and the co-activation of any SDR. For the CHR subjects, most of the SDRs showed no significant differences in the co-activation between the medicated and unmedicated CHR subgroups.

In summary, multiple INs showed significant changes in ESZ and CHR individuals. Among all regions identified, the most prominent group difference pattern involved reduced co-activation in the ESZ patients. Furthermore, similar network abnormalities were evident in individuals exhibiting the psychosis risk syndrome, although ESZ patients generally showed more severe changes. Many brain regions showed a progressive change from HC to CHR individuals to ESZ patients. Some regions also presented CHR-unique alterations. Additionally, co-activation strengths in some regions were associated with symptom severity ratings. In general, CHR individuals had INs with intermediate connectivity, falling between the normal connectivity present in HC individuals and the dysconnectivity evident in ESZ patients. Taken together, the IN abnormalities found in our study warrant further study in longitudinal designs in order to more precisely determine their roles in the pathogenesis of schizophrenia and other psychotic disorders.

This study had several limitations, including some that could limit the generalizability of our findings. First, a few parameters (e.g., the number of ICs used in GIG-ICA and the thresholds used for multiple comparisons correction) are adjustable and the specific settings we chose may have influenced the pattern of results we observed. Second, insufficient clinical follow-up data from our CHR participants prevented us from examining whether any of the IN abnormalities identified are predictive of conversion to psychosis. We will investigate this issue in the future when more data are available. Third, due to the relatively small sample sizes available for this study, we included a small number of CHR individuals who were treated with antipsychotic medication. Most of the ESZ patients analyzed took antipsychotic medication. Although we did not find significant associations between the medication and the identified measures showing group differences, influence of medication could be complex and deserves further investigation in future. In addition, we did not assess and cannot address whether treatment response or treatment resistance moderated functional connectivity in the ESZ patients, since we did not assess the treatment response information and also did not require stable doses of medication for entry into the study. However, to the extent that many of the IN abnormalities observed were evident in both ESZ and CHR individuals, there is a diminished likelihood that the network abnormalities observed in ESZ patients are the result of antipsychotic medication. Finally, since we flipped the estimated group-level ICs to ensure that their skewness was positive, in order to simplify the analyses for each IN, we only focused on the voxels with positive Z-scores in most subjects, i.e., those voxels that contributed the most to the networks (see Fig. S2). We ignored the voxels with negative Z-scores, thus we did not compute the co-activation within regions with negative Z-scores. We also did not investigate the anti-correlation between regions with positive Z-scores and regions with negative Z-scores. As shown in Fig. S2, the anti-correlation (e.g., the anti-correlation between the DMN and the salience network in IN 11) may deserve further study.