Predicting individual improvement in schizophrenia symptom severity at 1-year follow-up: Comparison of connectomic, structural, and clinical predictors.

In a machine learning setting, this study aims to compare the prognostic utility of connectomic, brain structural, and clinical/demographic predictors of individual change in symptom severity in individuals with schizophrenia. Symptom severity at baseline and 1-year follow-up was assessed in 30 individuals with a schizophrenia-spectrum disorder using the Brief Psychiatric Rating Scale. Structural and functional neuroimaging was acquired in all individuals at baseline. Machine learning classifiers were trained to predict whether individuals improved or worsened with respect to positive, negative, and overall symptom severity. Classifiers were trained using various combinations of predictors, including regional cortical thickness and gray matter volume, static and dynamic resting-state connectivity, and/or baseline clinical and demographic variables. Relative change in overall symptom severity between baseline and 1-year follow-up varied markedly among individuals (interquartile range: 55%). Dynamic resting-state connectivity measured within the default-mode network was the most accurate single predictor of change in positive (accuracy: 87%), negative (83%), and overall symptom severity (77%) at follow-up. Incorporating predictors based on regional cortical thickness, gray matter volume, and baseline clinical variables did not markedly improve prediction accuracy and the prognostic utility of these predictors in isolation was moderate (<70%). Worsening negative symptoms at 1-year follow-up were predicted by hyper-connectivity and hypo-dynamism within the default-mode network at baseline assessment, while hypo-connectivity and hyper-dynamism predicted worsening positive symptoms. Given the modest sample size investigated, we recommend giving precedence to the relative ranking of the predictors investigated in this study, rather than the prediction accuracy estimates.

INTRODUCTION

While many individuals with schizophrenia can achieve varying degrees of remission, the risk of relapse, recurrence and worsening of symptoms remains significant, particularly after treatment discontinuation (Emsley, Chiliza, Asmal, & Harvey, 2013). Relapse may indicate disease progression (Lieberman et al., 1993), and in many cases, is accompanied by brain changes (Cropley & Pantelis, 2014) or serious consequences, including extreme despair and increased risk of suicide (Pompili et al., 2007). Predicting relapse events and illness course is therefore critical to enable early and targeted initiation of appropriate treatments to individuals with a high risk of relapse. Moreover, individuals predicted to follow a favorable illness course may be spared unnecessary treatment and associated side‐effects (Arana, 2000; Bruijnzeel, Suryadevara, & Tandon, 2014; Cha & McIntyre, 2012; Fusar‐Poli et al., 2013; Vita, De Peri, Deste, & Sacchetti, 2012).

Predicting individual illness course, particularly the likelihood of specific relapse events, is very challenging. Different variables have been identified as predictors of outcome in the past, including: (a) demographics such as age of onset and sex (Eggers & Bunk, 1997; Fleischhaker et al., 2005; Hafner, Maurer, Loffler, & Riecher‐Rossler, 1993; Kydd & Werry, 1982) and family history of psychiatric illness (Fleischhaker et al., 2005); (b) cognitive measures such as IQ (Werry & McClellan, 1992; Werry, Mcclellan, & Chard, 1991) and processing speed (Milev, Ho, Arndt, & Andreasen, 2005; Puig et al., 2014; Remschmidt & Theisen, 2012); (c) baseline clinical symptoms (Maziade et al., 1996; Meng et al., 2006; Milev et al., 2005); (d) biological markers such as baseline antioxidant status (Fraguas et al., 2012; Martínez‐Cengotitabengoa et al., 2014) and neuroimaging markers (Greenstein, Wolfe, Gochman, Rapoport, & Gogtay, 2008; Khodayari‐Rostamabad, Hasey, MacCrimmon, Reilly, & de Bruin, 2010; Koutsouleris et al., 2018). The most accurate predictors of relapse and worsening of symptoms include nonadherence with medication (Robinson et al., 1999), poor premorbid adjustment, and substance abuse (Alvarez‐Jimenez et al., 2012), yet these predictors are unreliable in some cases, potentially difficult to ascertain and at most explain a three to fourfold increase in relapse risk. While biomarker discovery in schizophrenia remains a key research area (Singh & Rose, 2009), which specific types of variables are most useful in predicting outcome remains unclear (Díaz‐Caneja et al., 2015; Menezes & Milovan, 2000). Accurate, objective and reliable predictors of long‐term illness outcome are therefore needed to enable individualized prognostication and treatment.

Functional and structural neuroimaging is a promising candidate for such a predictor that has received significant attention. Several recent studies demonstrate the prognostic utility of training machine learning algorithms to recognize features derived from an individual's neuroimaging data that predict transition to psychosis (Koutsouleris et al., 2009, 2012, 2015; Ramyead et al., 2016), illness course (Kambeitz‐Ilankovic et al., 2016; Koutsouleris et al., 2018; Mourao‐Miranda et al., 2012), treatment response (Khodayari‐Rostamabad et al., 2010; Koutsouleris, Wobrock, et al., 2018), relapse (Nieuwenhuis et al., 2017), and resilience (de Wit et al., 2016). These studies have sought to compute individual prognostications at follow‐ups ranging between 1 and 7 years, primarily using anatomical features derived from structural neuroimaging, such as regional cortical thickness estimates or measures of gray matter volume. Prediction accuracies typically range between 60 and 90% (Janssen, Mourão‐Miranda, & Schnack, 2018).

More recently, the prognostic capability of resting‐state functional connectivity inferred from functional MRI has emerged in the context of predicting antipsychotic treatment response in individuals experiencing a first episode of psychosis (Cao et al., 2018; Sarpal et al., 2016). However, it remains unclear whether resting‐state connectivity enables improved prognostication in a machine learning setting compared to classic neuroimaging measures indexing cortical structure and morphology, such as cortical thickness and gray matter volume, and/or demographic and clinical variables. Furthermore, the dynamic properties of resting‐state connectivity particularly warrant investigation as candidate predictors of illness outcome. When compared to classic time‐averaged resting‐state functional connectivity approaches, dynamic measures of connectivity yield enhanced sensitivity to disease‐related effects (Damaraju et al., 2014) and improved diagnostic accuracy (Kottaram et al., 2018).

Schizophrenia is associated with widespread reductions in resting‐state connectivity strength (Bluhm et al., 2007; Fornito, Zalesky, & Breakspear, 2015; Fox & Greicius, 2010; Ganella et al., 2017; Liu et al., 2008; Salvador et al., 2007; Zalesky, Fornito, & Bullmore, 2010), which are associated with the severity of specific symptom dimensions (Hare et al., 2018; Repovs, Csernansky, & Barch, 2011; Rotarska‐Jagiela et al., 2010; Venkataraman, Whitford, Westin, Golland, & Kubicki, 2012). In particular, the default‐mode network, involved in interoceptive functions (Andrews‐Hanna, Smallwood, & Spreng, 2014; Buckner, Andrews‐Hanna, & Schacter, 2008; Raichle et al., 2001), has been consistently shown to exhibit aberrant connectivity in schizophrenia (Bluhm et al., 2007; Hu et al., 2017; Jafri, Pearlson, Stevens, & Calhoun, 2008; Liemburg et al., 2012; Liu et al., 2012; Shim et al., 2010; Skudlarski et al., 2010; Whitfield‐Gabrieli et al., 2009; Zhou et al., 2007). Assessing the temporal variations in functional connectivity revealed periods of both hypo and hyper‐connectivity with the default‐mode network (Damaraju et al., 2014; Kottaram et al., 2018) and reduced dynamism in patients (Kottaram et al., 2019; Miller et al., 2016).

In the present study, we aim to predict individual illness outcomes (i.e., worsening or improvement in symptoms) at 1‐year follow‐up in a cohort of individuals with schizophrenia. We train a simple (linear) machine learning algorithm to predict individual outcomes based on baseline measures of: (a) static and dynamic resting‐state connectivity (functional neuroimaging), (b) cortical thickness and gray matter volume (structural neuroimaging), and/or (c) clinical and demographic variables. We aim to determine which of these measures, or combinations thereof, provide the most accurate prediction of illness outcome. We hypothesized that dynamic measures of resting‐state connectivity within the default‐mode network would provide the most accurate and objective predictors of outcome. Our hypothesis is motivated by recent studies that suggest aberrant switching and connectivity dynamics within the default‐mode network are core pathophysiological features of schizophrenia (Damaraju et al., 2014; Du et al., 2016; Kottaram et al., 2018, 2019; Miller et al., 2016; Rashid et al., 2016). This emerging literature led us to specifically target the default‐mode network as a predictor, particularly given that: (a) the default‐mode network shows the highest levels of activation during rest (Buckner et al., 2008; Raichle et al., 2001); (b) while wide‐spread connectivity disruptions have been reported in schizophrenia, the most consistent and robust abnormalities have been reported in this network (Bluhm et al., 2007; Hu et al., 2017; Jafri et al., 2008; Liemburg et al., 2012; Liu et al., 2012; Shim et al., 2010; Skudlarski et al., 2010; Whitfield‐Gabrieli et al., 2009; Zhou et al., 2007); (c) from a methodological perspective, the default‐mode network can be more easily and robustly detected compared to other networks (Raichle & Snyder, 2007). This is the first study to evaluate the prognostic capability of dynamic measures in resting‐state connectivity. Given our modest sample size in the context of machine learning (N = 30), we consider this study preliminary and emphasize that our primary goal is to compare and rank the predictive power of various neuroimaging and clinical variables. This study can inform the choice of predictors (i.e., static vis‐a‐vis dynamic, functional vis‐a‐vis structural neuroimaging) investigated in future studies, which should be conducted in larger cohorts that are preferably acquired across multiple sites.

METHODS

Participants

Individuals with schizophrenia (N = 30; age = 27.3 ± 8.6 years; 10 females) were recruited through Orygen Youth Health clinical sites in Victoria, Australia. All participants had a schizophrenia‐spectrum disorder. More specifically, 11 individuals had an established diagnosis of schizophrenia based on DSM‐IV, while the remaining 19 individuals had recently experienced a first psychotic episode but were without an established diagnosis. All individuals were prescribed antipsychotic medications at the time of recruitment, image acquisition and follow‐up assessment. Exclusion criteria included a history of neurological disease or significant brain injury, a documented history of developmental delay or intellectual disability and presence of any contraindications to MRI. All individuals provided written informed consent prior to participation. The study was approved by the Melbourne Health Human Research Ethics Committee (MHREC ID 2012.069).

Baseline and follow‐up assessments of symptom severity, behavior and cognitive ability were completed for all individuals. Baseline assessments were performed within 1 week of image acquisition. Follow‐up assessments were performed after an interval of approximately 1 year (372 ± 12.7 days). Positive and negative symptoms were assessed using the Brief Psychiatric Rating Scale (BPRS; Andersen et al., 1989) and general cognitive ability was assessed using the Wechsler Abbreviated Scale of Intelligence (WASI; Wechsler, 1955). Additional assessments were administered but not analyzed in this study. Table 1 shows summary statistics at baseline and follow‐up. While group‐averaged ratings of symptom severity did not significantly differ between baseline and follow‐up (Table 1), marked heterogeneity was evident across individuals (Figure 1).

TABLE 1

Demographic, behavioral, and clinical characteristics

	Baseline	Follow‐up (1 year)	Baseline vs. follow‐up comparison
Sex (male/female)	20/10
Age (years)	27.3 ± 8.6	28.31 ± 8.6
Illness duration (years)	3 ± 4.6
IQ (WASI)	96 ± 19.2
Education (years)	12.4 ± 3.9
BPRS positive	8.9 ± 5.1	8.2 ± 4.5	p > .1
BPRS negative	5.2 ± 2.3	5.8 ± 2.9	p > .1
BPRS total	42.1 ± 11.6	44.1 ± 12.6	p > .1
Chlorpromazine equivalent dosage (mg/day)	412.9 ± 237.3	411.0 ± 213.5	p > .1

Note: BPRS, Brief Psychiatric Rating Scale; IQ, Intelligence Quotient; WASI, Wechsler Abbreviated Scale of Intelligence.

FIGURE 1

Relative change in symptom severity from baseline to 1‐year follow‐up. (a) Relative change in aggregate scores of positive (red, Δ) and negative (blue, Δ) symptom severity (b) Relative change in overall symptom severity ( Δ). Individuals ordered from most improved (leftmost) to most worsened (rightmost) in overall symptom severity. Symptom severity assessed with the Brief Psychiatry Rating Scale (BPRS). Binary classifiers were trained to predict individuals that worsened ( Δ) or improved ( Δ) with respect to symptom severity , where cut‐off thresholds of (Figure 3) and (Figure S6) were considered. Black dotted lines in Panels a and b correspond to 20% threshold

FIGURE 3

Accuracy of predicting improvement in symptom severity at 1‐year follow‐up using connectivity, structural, and/or clinical features. Individuals were classified as either worsening or improving/stabilizing based on relative change in positive, negative and overall symptom severity. Binary classifiers (linear discriminant analysis) were then trained to predict individual outcome using different combinations of features. (a) Prediction accuracies for classifiers trained using features from only one of the following five feature classes: cortical thickness, gray matter volume, baseline (BL) clinical, and demographic variables, static resting‐state functional connectivity (sFC) and dynamic resting‐state functional connectivity (dFC). A separate classifier was also trained using the top‐8 features selected from each of the above five feature classes as well as the combination of all of them (All—8 features). (b) Prediction accuracies for classifiers trained using 5 (All—40 features), 4 (No Thickness, No Volume, No BL Clin & Dem, No sFC, and No dFC), and 3 (No s and d FC) of the five feature classes. BPRS, Brief Psychiatric Rating Scale; FC, functional connectivity

Image acquisition and preprocessing

Functional and structural magnetic resonance imaging (MRI) was performed in all individuals. Acquisition details are provided in Supporting Information S1. In brief, T1‐weighted structural (TR = 1980 ms, TE = 4.3 ms and voxel resolution of 0.98 mm × 0.98 mm × 1.0 mm) and resting‐state (TR = 2 s, TE = 40 ms and voxel dimensions = 3.3 mm × 3.3 mm × 3 mm, 234 frames over ∼8 min) images were acquired on a 3T Siemens scanner.

Preprocessing was primarily performed using FSL 5.0.9 (FMRIB Software Library, www.fmrib.ox.ac.uk). Each individual's T1‐weighted image was skull‐stripped and spatially normalized via nonlinear registration to the Montreal Neurological Institute (MNI) 152 template. Resting‐state images were slice‐time corrected, realigned to the mean functional image to correct for motion, coregistered to the corresponding T1‐weighted scan via rigid‐body registration and then warped to the normalized structural image using the previously computed warp transformation. Images were sampled at 2 mm isotropic in MNI space. Images were then spatially smoothed using a Gaussian kernel of full width at half maximum of 4 mm. The smoothed images were further processed with ICA‐AROMA (Pruim et al., 2015), an automated tool for motion artifact removal, which has been shown to improve both sensitivity and specificity of connectivity analyses (Parkes, Fulcher, Yücel, & Fornito, 2018; Pruim, Mennes, Buitelaar, & Beckmann, 2015). Signals from white matter and ventricles were then regressed from the voxel time courses and the resulting residuals were used for all subsequent analyses. Finally, any linear trends were removed from each voxel time course and band pass filtering (0.01–0.1 Hz) was performed to suppress the effects of physiological noise (Cordes et al., 2001). Framewise displacement (Power, Barnes, Snyder, Schlaggar, & Petersen, 2012) was computed to quantify the extent of intrascan head movement (0.11 ± 0.04 mm).

Operationalization of clinical outcome

According to a recently proposed scheme (Dazzi, Shafer, & Lauriola, 2016), aggregate scores indexing the severity of positive (+), negative (−), and overall (∑) symptoms were computed by summing across relevant BPRS sub‐items (positive: 7 BPRS sub‐items, negative: 4, total: all sub‐items). Further details about sub‐item allocation are provided in Supplementary Table 1. For each individual, clinical outcome was operationalized as the relative change in symptom severity between baseline and one‐year follow‐up. More specifically, if denotes the value of aggregate score at baseline () and follow‐up (, we defined relative change as follows,

We classified individuals as either improving or worsening with respect to Δ, Δ , and Δ. In particular, worsening of symptoms was operationalized as Δ. Individuals not satisfying this criterion were deemed to improve (or remain stable). This operationalization allowed individuals to improve (worsen) with respect to positive symptoms but worsen (improve) with respect to negative symptoms. While dichotomizing the outcome measure can reduce statistical power, our rationale for doing so is to provide a clear and interpretable measure of outcome that is clinically relevant. Categorizing continuous measure can aid in making sense of unfamiliar measurement scales and treatment effects of uncertain implication (Lewis, 2004).

In this study, the cut‐off thresholds of = 0% (main analysis) and = 20% (supplementary analyses) were investigated. Figure 1 shows the relative change in Δ, Δ , and Δ for each individual.

Resting‐state functional connectivity

Resting‐state functional connectivity was computed between 18 cortical regions comprising the default mode network (DMN; Figure 2a), as defined elsewhere (Andrews‐Hanna, Reidler, Sepulcre, Poulin, & Buckner, 2010; Dodell‐Feder, DeLisi, & Hooker, 2014; Du et al., 2016; Kucyi & Davis, 2014). Refer to Supporting Information (Section S1, Table S2) for MNI coordinates. For each individual, the processed fMRI data were spatially averaged across the voxels comprising each of the 18 regions, yielding a representative time course for each region. A schematic of static and dynamic connectivity estimation is shown in Figure 2b and described below.

FIGURE 2

Schematic of functional connectivity estimation and feature combinations evaluated for individual prediction of change in symptom severity (a) Cortical regions comprising the default‐mode network. Functional connectivity was estimated between all pairs of regions. (b) Schematic of static and dynamic functional connectivity estimation. For each pair of regions, the Pearson correlation coefficient was computed across all time (static) or within overlapping windows (dynamic). (c) List of feature combinations evaluated for individual prediction of change in symptom severity using binary machine learning classifiers. Ticks indicate the presence of a particular class of feature. Feature selection was used to determine the top‐8 features within each feature class. Classifiers were trained on the top‐8 features from individual feature classes (eight features; Rows 1–5), combination of top‐8 features from each of the five feature classes (40 features in total, Row 7) as well as the top‐8 features across all feature classes (8 features in total; Row 6). In addition, combinations of the top‐8 features from different classes were considered, as indicated (Rows 8–13). dFC, dynamic FC; FC, functional connectivity; sFC, static FC

Static connectivity

Static resting‐state functional connectivity was computed for each of the (18 × 17) / 2 = 153 unique pairs of DMN regions. To this end, the Pearson correlation coefficient was computed between the representative time courses of each regional pair (Zalesky, Fornito, & Bullmore, 2012). This was repeated independently for each individual.

Dynamic functional connectivity

Dynamic functional connectivity was computed for each of the 153 unique pairs of DMN regions. Based on our previous study (Kottaram et al., 2018), time was first partitioned into overlapping and contiguous windows of length 20 s, using an overlap of one TR (2 s). The Pearson correlation coefficient was used to compute functional connectivity within each window, yielding a time‐resolved (dynamic) measure of connectivity for each regional pair. To reduce the dimensionality of the resulting data, the mean and standard deviation was computed across the set of windows for each pair of regions (Kottaram et al., 2018). Therefore, each individual's connectivity dynamics were characterized in terms of 153 mean values and 153 standard deviations.

Gray matter volume and cortical thickness estimation

Gray matter volume and cortical thickness were estimated regionally for each individual using T1‐weighted structural images and established software. The Desikan–Killiany atlas (Desikan et al., 2006) was used to define 35 bilateral cortical regions (70 in total). FreeSurfer (version 5.3.0, https://surfer.nmr.mgh.harvard.edu/) was then used to estimate cortical thickness and gray matter volume for each region (Dale, Fischl, & Sereno, 1999; Fischl & Dale, 2000; Fischl et al., 2002; Fischl, Sereno, & Dale, 1999). This involved delineating a three‐dimensional cortical surface model for each individual, which was then inflated and registered to the Desikan‐Killiany brain atlas. All gray/white and gray/pial boundaries were manually examined and corrected by a trained rater to increase the accuracy of the cortical estimates.

Baseline clinical and demographic variables

Clinical and demographic measures assessed at baseline were investigated as predictors of change in symptom severity at follow‐up. A total of eight baseline variables were considered: aggregate scores indexing the severity of positive, negative, and overall symptoms (Section 2.3), illness duration, age of onset of psychosis, sex (binary), intelligence quotient (IQ; WASI), and medication (chlorpromazine equivalent dosage, mg/day).

Classification

Binary machine learning classifiers were trained to predict individual improvement (Δ) or worsening (Δ) with respect to aggregate scores of positive, negative, and overall symptoms, . Under a leave‐one‐out cross‐validation (LOOCV) framework, separate classifiers were trained to predict improvement/worsening in each of the three aggregate scores. Further, classifiers were trained and evaluated using different combinations of features to determine the relative predictive power among different classes of features. We considered five feature classes: (a) cortical thickness (70 regions), (b) gray matter volume (70 regions), (c) baseline clinical and demographic variables (eight variables), (d) static connectivity (153 connections), and (e) dynamic connectivity (153 mean values +153 standard deviations). Feature selection was used to determine the top‐ most predictive features within each of the five feature classes. A two‐sample t‐test was used to perform feature selection within each feature class using only the training data. This involved assessing the between‐class separability of individual features, ranking all features based on the t‐statistic magnitude and then choosing the top‐ features. Among all feature classes, the baseline clinical and demographic variables contained the least number of features (n = 8); hence the number of features chosen was set to = 8 in the main analysis. Since our aim was to compare the predictive power of different feature classes, no further analysis was performed to optimize the value of ; however, all analyses were repeated with = 20 to investigate the impact of variation in the total number of features used. Importantly, in each iteration of the LOOCV, features were selected exclusively based on the training sample (N = 29), leaving the testing sample (N = 1) aside, thereby avoiding cross‐contamination of data. Using a relatively small set of features (i.e., = 8) minimized the risk of over‐fitting, although feature selection per se does not eliminate this risk.

To predict the improvement/worsening of each aggregate measure, 13 different classifiers were trained using different combinations of features, as listed in Figure 2c. A separate classifier was trained and evaluated for each of the five feature classes. Furthermore, classifiers were trained on the top‐ features across all feature classes ( features in total) as well as the top‐ features from each of the five feature classes (5 × features in total). In addition, to quantify the detriment of excluding a particular feature class, classifiers were trained on the top‐ features from 4 of the 5 feature classes (4 × features in total). Finally, a classifier was trained using the top‐ features from all feature classes except the connectivity‐based (static and dynamic) features (3 × features in total).

Given the modest sample size of the cohort studied, linear discriminant analysis (LDA) was employed; LDA models the distribution of features separately for each response class and minimizes over‐fitting, thus being particularly suited to small datasets (Izenman, 2013). Classification results can vary with the choice of machine learning algorithm and for this reason, ancillary analyses were then performed using other linear (support vector machine, SVM) and nonlinear (decision trees) classifiers. For each LOOCV analysis, established measures of classification performance were computed, including prediction accuracy, specificity, sensitivity, precision, area under the curve (AUC), and the F‐measure were calculated.

The feature selection described above is a univariate approach, as individual features are ranked independent of other features. We have also performed the classification experiments using a different feature selection method based on principal component analysis (PCA), which is further described in Supporting Information (Section S3).

RESULTS

Predicting individual change in symptom severity

Relative change in overall symptom severity between baseline and 1‐year follow‐up varied markedly among individuals (interquartile range: 55%; Figure 1), consistent with the heterogeneity in outcomes that is characteristic of schizophrenia. Some individuals worsened with respect to positive symptoms but showed an improvement in negative symptoms, and vice versa. Individuals were classified as either: (a) worsening or (b) improving (or remaining stable) with respect to positive (improving: N = 18), negative (improving: N = 19), and overall symptom severity (improving: N = 12) at 1‐year follow‐up. Classifiers were then trained to predict these individual outcomes using different combinations of features (predictors), with the aim of determining which of functional neuroimaging (static and dynamic resting‐state connectivity), structural neuroimaging (regional gray matter volume and cortical thickness) or baseline demographic and clinical variables, or combinations thereof, provide the most accurate prognostic utility.

Figure 3 shows prediction accuracies for the range of feature combinations indicated in Figure 2c. Combining features across all feature classes (i.e., classifiers labeled “All‐8 features” in Figure 3a) yielded the highest accuracies in predicting change in overall symptoms. Classification accuracies consistently exceeded 80% when features were combined across all feature classes. Remarkably, classifiers trained using only features associated with dynamic resting‐state connectivity yielded accuracies that were comparable to combining features across all feature classes. In contrast, the other single‐class classifiers trained using only one feature class yielded poorer accuracies (<70%, Figure 3a) and often did not yield predictions that exceeded chance level.

The predictive utility of dynamic resting‐state connectivity computed within the default‐mode network is further evidenced by the substantial reduction in prediction accuracies that were apparent when this feature class was excluded from the classifier combining all feature classes (Figure 3b). Excluding dynamic resting‐state connectivity resulted in a 10–20% reduction in accuracy, while excluding both static and dynamic connectivity features resulted in the most substantial accuracy reduction (>30% reduction). In contrast, excluding features related to cortical thickness, gray matter volume or baseline demographic and clinical variables did not affect prediction accuracies by more than 5% (Figure 3b). Table 2 shows additional measures that quantify the accuracy of the classifiers trained using dynamic resting‐state connectivity.

TABLE 2

Performance of predicting improvement in symptom severity at 1‐year follow‐up for classifiers trained using dynamic resting‐state connectivity

Performance measure	Positive symptoms	Negative symptoms	BPRS total
Accuracy (%)	86.67	83.33	76.67
Sensitivity (%)	75	81.82	72.22
Specificity (%)	94.44	84.21	83.33
Precision (%)	90	75	86.67
AUC	0.85	0.83	0.78
F‐measure (%)	81.82	78.26	78.79

Note: AUC—area under the receiver operating characteristic (ROC) curve.

Other performance measures evaluated, such as classification sensitivity (Figure S1), specificity (Figure S2), and AUC (Figure S3) accord with the above findings, suggesting that our classification models were not biased to predicting either improving or worsening outcomes alone. Selecting the top‐20 (instead of top‐8) features within each feature class did not markedly alter prediction accuracies (Figure S4). Compared to LDA (Figure 3), prediction accuracies decreased when using decision trees (Figure S5) or linear support vector machines (Figure S6). Despite this decrease, the relative ranking of features remained consistent, with dynamic resting‐state connectivity consistently providing the most accurate prognostic utility. Further, increasing the cut‐off threshold to = 20% (instead of = 0%) resulted in a marginal reduction (<5%) in the accuracy of predicting change in overall symptom severity, although prediction of change in positive and negative symptoms was unaffected. Also, the analysis with PCA‐based feature selection yielded similar results (Figure S8), consistent with the main analysis reported in Figure 3.

Finally, the variation in prediction accuracy in the case of holding out larger fractions of data for testing was evaluated. While all analyses so far employed a LOOCV approach, the held‐out sample size was systematically increased from 1 to 6 and the accuracies are reported in Figure S8. Only classifiers trained on dFC were considered in this supplementary analysis, given that this was the best performing feature class. While the accuracies reduced in general as the test sample size increased (thus yielding a smaller training sample), the accuracies remained >70%, despite holding out six samples for testing.

Characteristics of connections predicting outcome

Having found that dynamic resting‐state connectivity computed within an individual's default‐mode network yielded the most accurate prognostic utility, we next aimed to determine the characteristics and locations of the specific connections that were most informative with respect to prediction. We assumed a cut‐off threshold of = 0% in this section to operationalize improving and worsening subgroups of individuals with respect to relative change in symptom severity. Given that feature selection varied between iterations of the cross‐validation process, for this analysis, feature selection was performed on the full sample without any individual prediction. Instead, the top‐20% of connections associated with the largest effect sizes (two‐sample t‐test) were identified with respect to the null hypothesis of equality between the improving and worsening subgroups. In particular, a t‐test was computed independently for each pair of regions (i.e., each connection); such that a positive t‐statistic indicated that the worsening group showed increased connectivity at baseline, relative to the improving group.

Figure 4 shows the location of connections with the largest effect sizes (top‐20%), where positive effect sizes (warm colors) indicate increased connectivity strength (left column) and increased connectivity dynamics (standard deviation, right column) in the worsening group at baseline, relative to the improving group. In contrast, negative effect sizes (cool colors) indicate increased connectivity strength and increased connectivity dynamics in the improving group at baseline. It can be seen that worsening positive symptoms are primarily predicted by decreased connectivity strength but increased temporal dynamics among a distributed sub‐network spanning the entire DMN (Figure 4, uppermost row). Notably, connection strengths between ventral medial prefrontal cortex (vmPFC) and bilateral temporal poles, and between left temporal pole and right temporal parietal junction (TPJ) were lower in the worsening group. In addition, connectivity of left parahippocampal cortex (PHC), left temporal pole and left retrosplenial cortex to most other regions of DMN were reduced. In contrast, worsening negative symptoms are predominantly predicted by increased strength but decreased dynamics within a separate DMN sub‐network (Figure 4, center row). Connectivity of anterior medial prefrontal cortex (amPFC) to posterior cingulate cortex (PCC) and vmPFC to PHC was greater in the worsening group. In addition, left lateral temporal cortex (LTC) was strongly connected to other regions in the left hemisphere. The network that predicted worsening of negative symptoms appears to be largely localized to the left hemisphere, while the network predicting worsening of positive symptoms comprised many bilateral regions.

FIGURE 4

Connections most informative in predicting improvement in symptom severity at 1‐year follow‐up. Individuals were classified as either improving or worsening with respect to positive, negative, and overall symptom severity. Two‐sample t‐tests were performed to identify which connections differed most in mean connectivity strength (left column) and connectivity dynamics (standard deviation, right column) between the improving and worsening subgroups. The top‐20% of connections according to t‐statistic magnitude are shown. Note that t‐statistic magnitude is proportional to effect size. Each sphere represents a region within the default mode network (depicted in Figure 2a). The color and thickness of each connection represents the t‐statistic, as per the color bar. Boxplots represent the distribution of mean connection strength (left column) or standard deviation of connection strengths (right column), averaged over all regions within the default‐mode network. “*” indicates significant between‐group differences (p < .05). The top‐20% connections should not be interpreted as significant but are rather visualized here to identify those connection strengths which lead to high classification accuracies

Figure 4 suggests a marked contrast between predictors of change in positive versus negative symptom severity. In particular, worsening positive symptoms (compared to improving positive symptoms) were predicted by hypo‐connectivity and hyper‐dynamism in baseline resting‐state connectivity, while worsening negative symptoms (compared to improving negative symptoms) were predicted by hyper‐connectivity and hypo‐dynamism. We hypothesize that hyper‐connectivity within the DMN predicted worsening negative symptoms because increased engagement of this network has been associated with increased internal, self‐related processing (Molnar‐Szakacs & Uddin, 2013; Qin & Northoff, 2011), which is characteristic of rumination and other internalizing negative symptoms of schizophrenia (see Section 4).

In ancillary analyses, statistical inference was performed using the network‐based statistic (NBS) (Zalesky et al., 2010) to correct for multiple comparisons (i.e., 153 independent t‐tests). This enabled formal identification of the locations of connections that predicted outcome with respect to a family‐wise error corrected p‐value criterion (p < .05), rather than simply investigating the top‐20% as above. The primary t‐statistic threshold for the NBS was chosen as 2 to ensure that any significant effects exceeded an effects size (Cohen's d) of approximately 0.3, which is deemed a medium effect, and a total of 5,000 permutations were generated to estimate p‐values. This was repeated for changes in positive, negative and overall symptom severity. For overall symptom severity, the NBS identified subnetworks of connections that significantly differentiated the worsening and improving groups with respect to connectivity strength (p < .05) and dynamics (p < .05; Figure S7). These subnetworks were highly consistent with the top‐20% of connections (Figure 4; bottommost row). However, no significant effects were found with respect to positive and negative symptoms (p > .05).

DISCUSSION

Mounting evidence suggests that neuroimaging coupled with machine learning can provide utility in guiding long‐term prognostication in schizophrenia, and neuropsychiatry more generally; see Janssen and colleagues (Janssen et al., 2018) for a comprehensive review. However, it remains unclear: (a) whether neuroimaging can provide superior prognostic utility compared to clinical variables such as baseline symptom ratings; and (b) which neuroimaging modalities, or combinations thereof, provide the most accurate prognostic utility. To address these knowledge gaps, the present study aimed to determine which of three broad classes of candidate predictors would yield the most accurate prognostic utility in predicting relative change in schizophrenia symptom severity at 1‐year follow‐up.

We found that the dynamics of resting‐state functional connectivity within the default‐mode network (DMN) provided the most accurate prediction of whether an individual worsened or improved with respect to positive, negative and overall symptom severity at 1‐year follow‐up. Remarkably, while previous studies have reported prognostic utility of gray matter volume and cortical thickness measures (Greenstein et al., 2008; Khodayari‐Rostamabad et al., 2010; Koutsouleris, Wobrock, et al., 2018), prediction accuracies showed minimal improvement when incorporating structural neuroimaging predictors (regional cortical thickness and gray matter volume), static resting‐state connectivity and/or baseline demographic and clinical variables. The top‐8 selected features from each feature group for predicting worsening/improvement of positive, negative and overall symptoms are listed in Tables S3–S5. In particular, for the combination of all feature groups, the top‐8 selected features are mostly those from dFC. This explains why we did not observe any improvement in accuracy for the “All‐8 features.” However, note that even the classifier trained on the combination of top‐8 features from different feature groups (“All‐40 features”) did not provide any significant improvement in classification accuracies (Figure 3b). Similarly removing volume, thickness, and baseline clinical and demographic variables did not affect the classifier performance. In the absence of dynamic resting‐state connectivity information, these predictors alone yielded moderate‐to‐poor prognostic utility (<60% accuracy) that often did not exceed chance‐level predictions. A notable exception was static resting‐state connectivity, which yielded modest predication accuracies (60–70%). The difference in predictive prognostic powers of sFC and dFC measures supports the notion that static versus dynamic connectivity might be capturing different aspects of physiology (Liegeois, Laumann, Snyder, Zhou, & Yeo, 2017). While regional measures of brain structure did not substantially improve prediction performance here, further investigation is needed in larger cohorts to establish whether change in specific symptom dimensions and cognition is best predicted by structure and function.

Our best performing machine classifiers were able to predict whether an individual improved or worsened in terms of symptom severity with accuracies ranging between 75 and90%. However, these accuracies must be interpreted with caution given the modest sample size (N = 30), homogeneity of the sample (i.e., all individuals were recruited from the same clinical service) and general limitations of cross validation, particularly the risk of unreliable estimates with the LOOCV procedure used here (Varoquaux, 2018; Varoquaux et al., 2017). We emphasize that the main goal of the present study was not to provide precise estimates of prediction accuracy, but rather to rank the relative predictive power of various connectomic, structural and clinical predictors of individual outcome. If prediction accuracies were overestimated due to the limited sample size, any potential overestimation was most likely uniform across all feature classes, and thus the relative ranking of features (predictors) would remain unchanged.

While this is the first study to utilize connectivity dynamics to predict outcome, it is important to remark that several recent studies have investigated the prognostic utility of static resting‐state functional connectivity. Most of these studies have considered prediction of response to antipsychotic medication (Cadena et al., 2018; Cadena, White, Kraguljac, Reid, & Lahti, 2018; Cao et al., 2018; Kraguljac et al., 2016; Sarpal et al., 2015). For example, using task‐based fMRI data and a seed‐based analysis, Cadena and colleagues (Cadena, White, Kraguljac, Reid, Jindal, et al., 2018; Cadena, White, Kraguljac, Reid, & Lahti, 2018) found that connectivity between the anterior cingulate cortex and putamen was associated with treatment response over a 6‐week duration. Another study using resting‐state data and seed‐based analysis found that the connectivity of striatal regions to other brain regions was correlated with changes in symptoms over a 12‐week period (Sarpal et al., 2015). In addition, Cao and colleagues (Cao et al., 2018) reported that striatal connectivity was a predictor of response to antipsychotic medication over a period of 10 weeks in a medication‐naïve cohort.

We considered a significantly longer follow‐up period (1 year) than these recent studies; and furthermore, all individuals in our study were prescribed antipsychotic medication at baseline and follow‐up. Hence, our predictions were with respect to relatively long‐term change in symptom severity, as opposed to short‐term response to pharmacological treatment. It may be argued that accurately predicting treatment response is a more beneficial clinical capability than predicting long‐term change in symptom severity. However, it is important to remark that long‐term outcome is an important prognostic consideration given that schizophrenia can often become a chronic and lifelong disorder and such predictive information may orient treatment and follow‐up strategies, limiting relapses.

We found certain differences in dynamic connectivity features that predicted worsening positive versus negative symptoms (Figure 4, Tables S3 and S4). Individuals with worsening positive symptoms were predicted by hypo‐connectivity and hyper‐dynamism within the DMN at baseline, compared to individuals with improving positive symptoms. The converse pattern was evident for negative symptoms (i.e., hyper‐connectivity and hypo‐dynamism at baseline predicted worsening negative symptoms). While the reasons for this dissociation is unclear, it is consistent with the observation that several individuals improved with respect to one symptom dimension but worsened with respect to others (Figure 1). One hypothesis is that positive and negative symptom dimensions tap distinct and possibly mutually exclusive subsystems of the DMN, and thus predictors of change in positive and negative symptom severity circumscribe distinct loci. This notion is in line with the observation that the subnetwork predicting change in positive symptom severity (Figure 4, topmost row) shows little overlap with the connections predicting negative symptoms (center row). Worsening of negative symptoms was predicted by increased connectivity of a left‐lateralized sub‐network with direct involvement of hub regions such as PCC and amPFC, while for positive symptoms, a hypoconnectivity within a bilateral network without the participation of the hub nodes predicted worsening. However, the between‐group comparison of connection strengths using the NBS did not reveal any significant differences between the improving and worsening groups.

We hypothesize that this distinction between the predictors of worsening positive versus negative symptoms may relate to the putative introspective and extrospective modes of the DMN. The DMN toggles between an introspective self‐referential mode and an extrospective mode that remains alert to changes in the external environment (Fransson, 2005, 2006), suggesting a functional segregation within DMN regions (Fransson, 2005; Hearne, Cocchi, Zalesky, & Mattingley, 2015; Uddin, Clare Kelly, Biswal, Xavier Castellanos, & Milham, 2009). Hyperactivity in prefrontal regions characterizes internalization and self‐reflective thinking (Fransson, 2005; Ingvar, 1979), which is frequently interrupted and shifted toward an externalized state vigilant to external stimuli (Fransson, 2005; Laufs et al., 2003). Therefore, frontal DMN hyper‐connectivity may predict worsening of negative symptoms due to increased internalization and fewer shifts to the externalized mode, which are characteristic of the internalizing negative symptoms of schizophrenia. Disengagement among DMN regions (i.e., hypo‐connectivity) on the other hand leads to externalization (Farina et al., 2018), which may be necessary to engender delusions, hallucinations and other positive symptoms.

While we recommend that the dynamics of resting‐state functional connectivity should be given priority when selecting neuroimaging predictors of schizophrenia outcome in future studies, the precise biological and pathological correlates of connectivity dynamics are unknown. In this study, connectivity dynamics were summarized using the standard deviation computed across multiple temporal windows. Therefore, larger temporal fluctuations in functional connectivity strength were indexed by larger standard deviations. Connectivity dynamics could potentially index the rate at which the DMN cycles between putative internalizing and externalizing modes, where an abnormal switching rate may place individuals at greater risk of poor outcome. The hypo‐dynamism associated with worsening of negative symptoms (i.e., low standard deviation) suggests that individuals do not (or are unable to) freely switch between externalizing and internalizing states, consistent with internalization and rumination.

It is worth noting that worsening and improvement were defined as relative deviations from baseline scores in the present study. While this is a commonly followed approach (Cao et al., 2018; Leucht, Davis, Engel, Kane, & Wagenpfeil, 2007; Schennach, Riedel, Musil, & Möller, 2012), this definition does not take into account the actual symptom values. For example, as per this definition, an individual with a low score at baseline and a slightly higher score at follow‐up will be considered as worsening, even though the actual score at follow‐up might still be relatively low. Similarly, a slight improvement in symptoms when the actual scores are high does not imply a good outcome. Note that we performed a supplementary analysis where a threshold of = 20% was considered, such that an individual will be considered as worsening only if the percentage increase in follow‐up symptoms is >20%. This analysis yielded similar results (Figure S7), consistent with the main analysis ( = 0%) reported in Figure 3. However, we considered only percentage changes in symptoms in this study and no absolute changes were considered. Future studies are warranted to assess if the neuroimaging‐based classifiers are predictive of symptom severity (above or below a clinically meaningful cut‐off) or only sensitive to relative change in symptoms.

Several limitations should be considered. First, our sample size is modest for a machine learning study. It is worth mentioning that acquiring longitudinal neuroimaging data with long follow‐up periods is difficult in psychiatric cohorts. While we followed a cross‐validation framework which reduces the risk of over‐fitting, it does not eliminate the risk completely, particularly when all data are acquired from the same site. Intersite cross‐validation in which one acquisition site is used for training and another for testing would have further mitigated this risk. Given this, the reported prediction accuracies must be interpreted cautiously, particularly with respect to the high variability of cross‐validation estimates when applied to typically sized neuroimaging datasets (Schnack & Kahn, 2016; Varoquaux et al., 2017). For this reason, we used a simple linear classifier to minimize the risk of overfitting and principally focused on ranking the feature classes investigated, as opposed to interpreting the prediction accuracies in an absolute sense. We also used a relatively small set of predictors to further reduce the risk of overfitting as well as domain‐specific knowledge to select the regions and networks considered. Second, the feature class that was found to yield the most accurate predictions was also the class with the greatest number of features (153 mean values + 153 standard deviations), and thus the risk of overfitting may have been increased for this class. However, feature selection was performed to choose the same number of features across different classes. Moreover, prediction accuracy did not markedly improve when combining features from the five feature classes investigated, suggesting that the influence of overfitting on accuracy was possibly minimal; otherwise the classifier combining features across all classes would have yielded the highest accuracy. Third, our feature class with clinical and behavioral measures was limited to the information that was collected during data acquisition, which does not include all clinical/behavioral variables. There could be other variables with better predictive power that were not collected and included in this study. For example, incorporating structural connectivity inferred from diffusion‐weighted MRI into the prediction model could potentially enhance accuracy by explaining unique variance that is not captured by connectivity dynamics. Fourth, all individuals were prescribed antipsychotic medication and medication represents a potential confounding effect on structural and functional neuroimaging. The adherence to medication and certain environmental factors such as substance use during the period between baseline and follow‐up examinations affect the outcomes, but these cannot be accounted in a predictive model based on the baseline measures alone. Also, no information regarding adherence to medication was collected during the follow‐up assessment; hence we cannot comment on the extent to which the current findings relate to treatment response. Fifth, all participants were recruited from the same clinical service, resulting in a geographically and ethnically homogeneous cohort. While this homogeneity is beneficial in alleviating confounds inherent to interindividual variation, it is unclear whether the classifiers trained will generalize to other cohorts (Hahn, Ebner‐Priemer, & Meyer‐Lindenberg, 2019). Finally, the use of sliding windows to estimate connectivity dynamics has been criticized (Hindriks et al., 2016; Laumann et al., 2016), although we show that this methodology can yield predictors that outperform cortical thickness and gray matter volume.

CONCLUSIONS

Functional neuroimaging provides improved prognostic utility in predicting long‐term changes in schizophrenia symptom severity, compared to predictors derived from structural neuroimaging and baseline clinical variables. While our results might not necessarily generalize to other cohorts due to the modest sample size, our study can inform and constrain the choice of candidate predictors investigated in future machine learning studies that aim to predict illness outcome using larger cohorts. In particular, we recommend that dynamic fluctuations in resting‐state functional connectivity, particularly within the default‐mode network, are given prime consideration as candidate predictors in future studies.

CONFLICT OF INTEREST

The authors do not have any conflict of interest to declare.

Appendix S1. Supporting Information.

Click here for additional data file.