Distinct multivariate brain morphological patterns and their added predictive value with cognitive and polygenic risk scores in mental disorders.

The brain underpinnings of schizophrenia and bipolar disorders are multidimensional, reflecting complex pathological processes and causal pathways, requiring multivariate techniques to disentangle. Furthermore, little is known about the complementary clinical value of brain structural phenotypes when combined with data on cognitive performance and genetic risk. Using data-driven fusion of cortical thickness, surface area, and gray matter density maps (GMD), we found six biologically meaningful patterns showing strong group effects, including four statistically independent multimodal patterns reflecting co-occurring alterations in thickness and GMD in patients, over and above two other independent patterns of widespread thickness and area reduction. Case-control classification using cognitive scores alone revealed high accuracy, and adding imaging features or polygenic risk scores increased performance, suggesting their complementary predictive value with cognitive scores being the most sensitive features. Multivariate pattern analyses reveal distinct patterns of brain morphology in mental disorders, provide insights on the relative importance between brain structure, cognitive and polygenetic risk score in classification of patients, and demonstrate the importance of multivariate approaches in studying the pathophysiological substrate of these complex disorders.

Introduction

Schizophrenia and bipolar disorder are complex disorders with high heritability (van Os and Kapur, 2009), and are along with other mental and behavioral disorders the largest global contributors to years lived with disability (Vos et al., 2013). Whereas recent large-scale collaborative efforts have provided some clues about the genetic underpinnings (Psychiatric Genomics Consortium, 2014), the search for genes or single-nucleotide polymorphisms (SNPs) with strong impact on clinical expression has largely proven unsuccessful. This has reinforced the search for distinct intermediate phenotypes, quantitative biological traits that are reasonably heritable (Preston and Weinberger, 2005), which may be more strongly associated to gene function compared to diagnostic categories (Gottesman and Gould, 2003).

Both clinical expression and empirical evidence suggest molecular, neurochemical, and macrostructural brain abnormalities in severe mental disorders (van Os and Kapur, 2009). In line with the clinical and pathophysiological heterogeneity evident in these disorders, the brain underpinnings are multidimensional, reflecting a myriad of complex pathological processes and environmental effects. Multivariate techniques are well suited to disentangle the different sources of variability in brain imaging data, where each source presumably reflects independent biological pathways with unique genetic and environmental determinants. Using these methods to determine the underlying structure of brain morphology changes in severe mental disorders could provide novel clues about disease mechanisms and candidate avenues for novel and personalized interventions (Insel and Cuthbert, 2015). Partial least squares methods, for instance, which are based on the hypotheses of a linear relationship between the observations (predictor variables) and dependent variables and maximize the covariance between the two set of variables (Sui et al., 2012), have been used for multimodal fusion of function magnetic resonance imaging (MRI) and EEG (Martınez-Montes et al., 2004) as well as PET and structural MRI (Chen et al., 2009). Supervised machine learning regression approach was used to predict brain development index and assess the relevant multiparametric imaging patterns from T1 and diffusion tensor images (Erus et al., 2015). Linked independent component analysis (LICA), a probabilistic technique based on Bayesian framework, is a promising data-driven multimodal fusion technique for modeling co-variances across modalities (Groves et al., 2011). Its potential application has been demonstrated in multimodal studies of healthy subjects (Douaud et al., 2014, Groves et al., 2012), Alzheimer's disease spectrum (Doan et al., 2017) and patients with attention deficit hyperactivity disorder (Francx et al., 2016) and schizophrenia (Brandt et al., 2015).

To identify distinct brain morphology patterns in schizophrenia and bipolar disorder, we performed data-driven fusion of brain imaging phenotypes including cortical thickness, surface area and gray matter density maps using LICA (Groves et al., 2011, Groves et al., 2012). This multivariate technique enables data-driven decomposition of multimodal imaging features into a linear combination of independent components, each of which may represent a biologically interpretable mode of brain variation. Cortical surface area and thickness are important morphological indices of the brain cerebral cortex. Although the exact neurobiological constituents remain unclear, surface area and cortical thickness may be determined by the number and the laminar patterning of cortical columns, respectively (Chen et al., 2013, Rakic, 2009), thus providing a biologically relevant decomposition of cortical volume. Gray matter density map provides a mixed measure of gray matter that may partly reflect interactions between cortical thickness and surface area, but which is also likely reflecting relevant brain morphological variance that is not modeled by the surface-based thickness and area measures (Hutton et al., 2009). Surface-based thickness and area contribute to only a proportion of the gray matter density alterations in SZ and considering both surface-based and voxel-based measures in SZ studies may therefore both increase sensitivity (Palaniyappan and Liddle, 2012), and provide important information about the inter-dependence between the three measures across the brain.

We next applied machine learning, or multivariate pattern analysis, to identify robust combinations of multivariate brain morphology patterns that were most discriminative between groups. Here, our primary aim was not to build a classification model with maximal performance compared to existing models. Instead, by means of multivariate machine learning classification, we aimed to gain further insights regarding the most discriminative combinations of features and the importance of each feature with respect to the others in such combinations. To study their sensitivity to clinical and cognitive traits, we tested for associations with diagnosis, in addition to a range of cognitive domains and clinical variables. Polygenic risk score (PGRS) quantifies the additive effects of a large number of genetic variants with individually weak effects on complex traits. It is computed for each individual by summing up the effect sizes of a set of SNPs, selected based on a p-value threshold and weighted by the number of individual risk alleles (Purcell et al., 2009, Tesli et al., 2014). PGRS has been used as a useful genetic marker for evaluating the effect of cumulative genetic risks on the brain phenotypes (Kanai and Rees, 2011, Kauppi et al., 2015). Therefore, we also assessed to which degree the multivariate brain morphology patterns are modulated by PGRS for schizophrenia and bipolar disorder. Finally, to test the added predictive value provided by multivariate brain imaging features when used in combination with cognitive performance and polygenic risk, we performed and compared group classification using the different combinations of these feature sets. Based on previous reports on brain morphometric alterations in schizophrenia and bipolar disorders (Elvsåshagen et al., 2013, Hibar et al., 2017, Hibar et al., 2016, Kuperberg et al., 2003, Moberget et al., 2017, Rimol et al., 2012, van Erp et al., 2016a, van Erp et al., 2016b), we hypothesized that LICA would reveal morphological patterns that not only converge with univariate findings such as the characteristic pattern of fronto-temporal thinning in the patient groups, but also show novel modes of brain morphometric variability involving multiple morphometric measures sensitive to diagnosis. Based on the previously reported brain structural heterogeneity in patients relative to controls (Kambeitz et al., 2015, Nenadic et al., 2012), we expected to observe moderate discriminative power of the resulting multivariate features at an individual level. Lastly, we hypothesized that combining brain morphology, polygenic and cognitive scores would yield improved predictive value compared to using each of these feature sets alone.

Material and methods

Participants

Participant demographics are summarized in Table 1. Briefly, we included 223 participants with schizophrenia spectrum disorders (age = 32.1 ± 9.3 years, 95 women, 168 schizophrenia, 22 schizophreniform, 33 schizoaffective), 190 with bipolar spectrum disorders (age = 35.0 ± 11.3 years, 111 women, 117 bipolar disorder I, 64 bipolar disorder II, and 9 bipolar disorder not otherwise specified), and 284 healthy controls (age = 35.2 ± 9.6 years, 134 women). All participants were recruited as part of the Thematically Organized Psychosis (TOP) study, which was approved by the Regional Committee for Medical Research Ethics and the Norwegian Data Inspectorate. Written informed consent was obtained from all participants. Patients were recruited from both inpatient and outpatient clinics (but mostly outpatient) at major psychiatric hospitals in the Oslo area with the following criteria: aged between 18 and 65 years, understood and spoke a Scandinavian language, had no history of severe head trauma, and obtained an IQ score of above 70. Healthy controls were excluded if they or any of their first-degree relatives had a lifetime history of a severe psychiatric disorder (schizophrenia, bipolar disorder or major depression). Patients were clinically characterized through a personal interview conducted by trained physicians or clinical psychologists. The assessment covered diagnostics, symptomatology, neurocognition, drug use and medication status (Simonsen et al., 2011). Psychiatric diagnosis was established using the Structured Clinical Interview for DSM-IV Axis I Disorders (SCID) (First et al., 1996). An overview regarding family history of the patients is presented in Supplementary Table 3. The controls were randomly sampled from national registries from the same catchment area and age range as the patients, and had no psychiatric or alcohol/substance use disorder, as well as no cannabis use the last 3 months.

Table 1

Sample characteristics.

Group	SZ	BD	HC	Group comparisons
				SZ vs. BD	SZ vs. HC	BD vs. HC
Demographics, education
Total N	223	190	284
Sex (% female)	95 (42.6)	111 (58.4)	134 (47.2)	χ2 = 9.65, p = 1.89E-03	χ2 = 5.32, p = 2.11E-02	χ2 = 0.88, p = 3.48E-01
Age (years)	32.10 ± 9.31	34.99 ± 11.31	35.22 ± 9.60	t = − 2.73, p = 1.96E-02	t = − 3.42, p = 1.99E-03	t = − 0.37, p = 1
Education (years)	12.75 ± 2.49	13.50 ± 2.27	14.15 ± 2.26	t = − 2.36, p = 5.62E-02	t = − 5.43, p < 1E-04	t = − 2.81, p = 1.54E-02
Ethnicity (% caucasian)	180 (80.7)	171 (90)	280 (98.5)
Cognitive domains mean z ± sd (n)
g	− 1.64 ± 2.45 (181)	− 0.01 ± 2.20 (177)	1.08 ± 1.58 (278)	t = − 7.58, p < 1E-04	t = − 14.59, p < 1E-04	t = − 6.19, p < 1E-04
Processing speed	− 0.63 ± 0.92 (193)	− 0.13 ± 0.92 (177)	0.51 ± 0.82 (283)	t = − 5.76, p < 1E-04	t = − 14.90, p < 1E-04	t = − 8.33, p < 1E-04
Verbal learning/memory	− 0.45 ± 0.83 (194)	0.06 ± 0.90 (178)	0.25 ± 0.72 (284)	t = − 6.07, p < 1E-04	t = − 9.83, p < 1E-04	t = − 2.99, p = 8.56E-03
Executive function	− 0.53 ± 0.87 (193)	− 0.01 ± 0.69 (178)	0.37 ± 0.49 (284)	t = − 7.08, p < 1E-04	t = − 14.41, p < 1E-04	t = − 6.41, p < 1E-04
Working memory/attention	− 0.35 ± 0.89 (194)	− 0.07 ± 1.01 (177)	0.28 ± 0.99 (284)	t = − 3.17, p = 4.79E-03	t = − 7.33, p < 1E-04	t = − 3.72, p = 6.5E-04
Symptomsamean ± sd (n)
Duration of illnessb	10.33 ± 8.15 (220)	13.69 ± 9.91 (188)	na	t = − 3.7, p = 2.4E-04
PANSS negative	15.55 ± 6.59 (220)	10.11 ± 3.60 (186)	na	t = 10.5, p < 1E-04
PANSS positive	14.86 ± 5.38 (219)	10.08 ± 3.64 (187)	na	t = 10.6, p < 1E-04
PANSS total	62.18 ± 17.21 (217)	45.59 ± 10.12 (186)	na	t = 12, p < 1E-04
Drug use n (%)
Alcohol use last month	131 (58.7)	139 (73.2)	213 (75.0)
Amphetamine last month	8 (3.6)	2 (1.1)	1 (0.4)
Cocaine use last month	4 (1.8)	3 (1.6)	2 (0.7)
Cannabis use last month	17 (7.6)	16 (8.4)	3 (1.1)
AUDIT score mean ± sd (n)	6.51 ± 5.87 (198)	8.05 ± 5.83 (178)	6.10 ± 2.93 (247)
Medication n (%)
Antipsychotic	179 (80.3)	74 (38.9)	na
Lithium	2 (0.9)	35 (18.4)	na
Antiepileptic	29 (13.0)	82 (43.2)	na
Antidepressant	57 (25.6)	52 (27.4)	na
Hospitalization (mean ± sd)	2.4 (3.3)	1.6 (2.3)	na
Polygenic risk scores (PGRS) (total N = 505)
N	145	144	216
SZ PGRS (mean ± sd)	0.34 ± 0.94	0.11 ± 1.01	− 0.30 ± 0.94	t = 1.72, p = 2.57E-01	t = 6.05, p < 1E-04	t = 4.13, p = 1.26E-04
BD PGRS (mean ± sd)	0.03 ± 0.99	0.22 ± 0.96	− 0.17 ± 1.01	t = − 1.62, p = 3.16E-01	t = 1.84, p = 2E-01	t = 3.61, p = 1E-03

AUDIT: Alcohol Use Disorders Identification Test. SZ = Schizophrenia, BD = Bipolar Disorder, HC = Controls.

The median time period between the clinical assessment and MRI scan was 155 days.

Defined as years between age at the first occurrence of any psychiatric illness episode (psychotic, affective, or other) and date of MRI scan.

MR acquisition

Magnetic resonance imaging data were obtained on a 1.5 Tesla Siemens MAGNETOM Sonata scanner (Siemens Medical Solutions, Erlangen, Germany) supplied with a standard head coil. For each participant, two T1-weighted images were acquired using a repeated 3D T1-weighted magnetization prepared rapid acquisition gradient echo (MPRAGE) sequence with the following parameters: repetition time (TR) = 2730 ms, echo time (TR) = 3.93 ms, inversion time (TI) = 1000 ms, field of view (FOV) = 240 mm, flip angle (FA) = 7o, matrix = 192 × 256, voxel size = 1.33 × 0.94 × 1 mm, 160 sagittal slices. The two T1-weighted scans obtained for each participant were averaged after rigid registration to improve signal-to-noise ratio (SNR).

Image pre-processing

T1-weighted scans were processed using FreeSurfer (http://surfer.nmr.mgh.harvard.edu) to estimate vertex-wise cortical thickness and surface arealization (Dale et al., 1999). Surface maps were resampled to a common coordinate system (fsaverage5, 10,242 vertices) using a non-rigid high-dimensional spherical averaging method to align cortical folding patterns (Fischl and Dale, 2000). FSL-VBM (Douaud et al., 2007) was used to derive gray matter density (GMD) maps (see Supplementary Information (SI)).

Thickness maps were smoothed using a Gaussian kernel with a commonly used full width of half maximum (FWHM) of 15 mm, resulting in a smoothness of 20.8 mm. To match smoothness across surfaces, the surface area maps were smoothed with different values of FWHM (ranging from 10 mm to 15 mm) and maps having approximately the same smoothness as cortical thickness (corresponding to a FWHM of 11 mm) were used for further processing. In accordance with the preprocessing steps described in (Groves et al., 2012), the GMD maps were resampled to have an isotropic voxel size of 4 × 4 × 4 mm3 using a trilinear interpolator, and smoothed with a FWHM of 10.5 mm (sigma = 4.5 mm).

Linked independent component analysis

We performed a data-driven decomposition of the imaging features obtained from all subjects into independent components using the LICA implementation by FMRIB, University of Oxford (http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FLICA (Groves et al., 2012). LICA was developed based on the conventional ICA technique, which assumes the signal to be a linear mixture of statistically independent spatial patterns that are non-Gaussian. The subject weights, or mixing parameters, are unknown. During the optimization process, ICA searches for maximally non-Gaussian patterns by iteratively updating the subject weights. To model covariance patterns across modalities, LICA allows simultaneous ICA decompositions on different modalities but constrains the subject weights to be the same across modalities. Furthermore, LICA balances the information content from different modalities, allowing a modality to be absent from a component. A LICA component is thus characterized by its spatial maps (one per modality) and the subject weights that are shared across measures, and may involve multiple (hereafter generally referred to as multimodal components) or only one measure (unimodal components). Whereas the group spatial maps indicate the spatial variability at the group level, the subject weights indicate the relative contribution of a subject to the component. For instance, for a given modality in a component, if its spatial map shows positive pseudo z-score, then larger subject weights would reflect relative increases in the corresponding modality; and vice versa, for a spatial map with negative values, larger subject weights would reflect relative decreases. To aid visualization and interpretation of the anatomical distribution of an effect of interest in subsequent analyses on the subject weights, we thresholded the spatial maps (3 < | z | < 10), in line with previous publications (Brandt et al., 2015, Groves et al., 2012), with stronger effects implicated in regions with larger z-scores. Model order selection procedures based on cophenetic coefficient (Ray et al., 2013) suggested that fifty components provided a good fit with the data (see SI), with each component representing a distinct mode of brain variation.

Statistical analyses

Univariate analyses of the LICA subject weights

We applied general linear models (GLM) to test for the main effects of group, age and sex, as well as of the interaction between age and group, on the subject weights of each LICA component (group and sex coded as categorical variables R's lm function). The group pairwise contrasts were computed using the lsmeans R package (Lenth, 2016). We assessed the effect of intracranial volume (ICV) on the main findings by additionally including this measure as a covariate in the group effect analyses, accounting for age and sex. We further checked for the association between the components and total lateral ventricular volume, accounting for age, sex and ICV. We studied the pairwise correlations between all components by means of full and partial correlation.

Further, we assessed the associations between LICA subject weights and cognitive domain scores (g factor, processing speed, verbal learning/memory, executive function, working memory/attention), symptom domain (duration of illness, negative, positive and total PANSS) and polygenic risk scores (Table 1 and SI). We included measures of cognitive domains previously found to be impaired in schizophrenia and bipolar disorder (Bowie and Harvey, 2005, Simonsen et al., 2011). Possible confounding roles of ethnicity and subject motion, quantified using signal-to-noise ratio as a proxy, on the main effect and cognitive analyses were assessed and presented in SI. For the LICA components showing significant group effects, we explored possible effects of medication, substance use (drug, alcohol), and education within either schizophrenia or bipolar disorder group by including each of these variables in the GLM.

Correction for multiple testing across all LICA components was performed using permutation testing. The subject weights were permuted 10,000 times with respect to age, sex, and group, as well as cognitive, polygenic risk, and symptom scores when applicable. In particular, for the main group effect analysis, in each permutation, the maximum value of the f statistics associated to the main group effect across all components was computed and stored. This entire process resulted in a null distribution of the f statistics consisting of 10,000 values and the original f statistic obtained from the non-permuted data. The corrected p value was computed as the ratio between the number of values in the null distribution greater or equal to the original f statistic and the number of values in the null distribution. A similar permutation procedure was applied to obtain the corrected p values for the group pairwise comparisons. Unless specified otherwise, the p values reported throughout the manuscript are corrected for multiple comparisons using permutation testing as this non-parametric procedure only make weak assumptions and can provide the exact control of false positives from the data (Winkler et al., 2014). However, to facilitate comparisons with previous literature in which FDR correction is commonly used, for the main effect analyses, we also report FDR-corrected values (Benjamini and Hochberg, 1995) (Supplementary Table 1).

Multivariate machine learning analyses

Finally, we performed pairwise group classification based on a leave-one-out cross-validation approach using as features the LICA components after regressing out the linear effect of age and sex. We used the random forest classifier as implemented in the randomForest R package (Liaw and Wiener, 2002). For each classification, one subject was left out as a testing set and a classifier was trained on the remaining subjects. During training, the optimal number of features that should be randomly chosen for splitting at each tree was estimated such that it minimized the out-of-bag error (Breiman, 2001) on the training data using the tuneRF function (the testing set remained untouched). The trained classifier was then used to predict the class of the testing set. This process was repeated for each and every subject. Due to the imbalanced sample sizes across groups, to alleviate biased performance towards higher accuracy on the majority class, during training of the random forest classifier, we repeatedly sampled with replacement an equal number of subjects per class, and additionally, we used the area under the receiver-operating-characteristic curve (AUC) as the main performance measure. We repeated this classification using only data from subjects with a diagnosis of schizophrenia (N = 168) in the schizophrenia group, subjects with a diagnosis of bipolar I (N = 117) in the bipolar group, and the controls. On a subset of the sample for which imaging, cognitive and genetic data are available (schizophrenia = 122, bipolar disorder = 132, controls = 212), we performed group classifications using the different combinations of these different feature sets in the same cross-validation framework as described above to assess their complementary predictive value. Significance in terms of classification improvement was evaluated using permutation testing (more details can be found in SI).

LICA decomposition was performed using Matlab (version R2014a). All statistical analyses and classification were performed in R (http://cran.r-project.org, version 3.2.1).

Results

Univariate analyses of the LICA subject weights

Association with diagnosis

Six components (IC1,2,5,8,12,25 explaining 15.5%, 12.3%, 1.7%, 1.6%, 1.3%, 0.95% of the data variance (Supplementary Fig. 1), respectively) showed significant (p < 0.05, corrected) main effects of diagnosis (Fig. 1). Eleven (IC1,2,3,5,8,9,11,12,16,20,22) showed significant main effects of age (see SI, Supplementary Fig. 3), and nine (IC1,5,8,9,11,12,14,25,28) showed significant main effects of sex. Detailed results are presented in Supplementary Table 1.

Fig. 1

(a) Spatial maps of six LICA components showing significant group effects. Volumetric map is presented for GMD, whereas surface maps are presented for cortical thickness and surface area. Only spatial maps of modalities with considerable contribution to the components are shown. All maps were thresholded with | z | > 3. SA = surface area, CT = cortical thickness. Whereas IC1,2,12,25 were unimodal, i.e. dominated by a single modality, IC5,8 were multimodal and driven by variation in both CT and GMD maps. (b) Subject weight distribution of the above six components. The diamond signs represent the group means. The plotted values are the residuals obtained by regressing out age and sex effects. SZ = Schizophrenia, BD = Bipolar Disorder, HC = Controls.

IC1 reflects global cortical surface area with schizophrenia showing smaller subject weights than bipolar disorder (t = − 4.26, p = 3.8E-03) and controls (t = − 4.71, p = 5E-04). This component correlated positively with ICV (Buckner et al., 2004) (t = 31.70, p < 1E-04). Importantly, however, the group effect remained when including ICV in the GLM (F = 10.04, p = 2.6E-03, schizophrenia < bipolar disorder: t = − 4.15, p = 6.2E-03, schizophrenia < controls: t = − 3.61, p = 4.7E-02). IC2 reflects global cortical thickness. Schizophrenia and bipolar disorder patients showed decreased subject weights compared to controls (schizophrenia < controls: t = − 7.31, p = 1E-04, bipolar disorder < controls: t = − 4.88, p = 2E-04), indicating global cortical thinning in both patient groups.

IC5 is a multimodal component mainly driven by thickness (54% weight) and GMD (38% weight), implicating cingulate, insular and medial frontal thickness, as well as global variation with an opposite weighting in GM and WM. Schizophrenia (t = 5.4, p = 1E-04) and bipolar disorder (t = 5.31, p = 1E-04) had larger subject weights than controls, indicating increased thickness and GMD and decreased white matter density (WMD) in the patient groups. This component showed a negative association with total WM volume (t = − 13.88, p < 1E-03).

IC8 reflects multimodal involvement of lateral temporal, parietal and medial frontal thickness (27% weight) as well as a bilateral pattern of temporal occipital fusiform cortex, lingual cortex, precentral, postcentral gyrus and hippocampus GMD (65% weight). Compared to controls, schizophrenia (t = − 4.28, p = 3.7E-03) and bipolar disorder (t = − 3.63, p = 4.35E-02) showed decreased subject weights, indicating reduced thickness and GMD in these regions.

IC12 revealed a distinct gradient in thickness along the anterior-posterior axis with positive and negative weights in the frontal and posterior lobes, respectively. Compared to controls, schizophrenia (t = − 3.84, p = 1.97E-02) and bipolar disorder patients (t = − 3.34, p = 1.2E-01) showed smaller subject weights, indicating a shift in thickness distribution across the brain surface with relatively less anterior compared to posterior thickness in patients relative to controls.

IC25 reveals reduced cerebellar GMD in schizophrenia compared to controls (t = − 3.43, p = 8.75E-02, FDR corrected p = 1.19E-02).

ICV also showed an association on IC5 and IC12 (IC5: t = − 6.8, p < 1E-03, IC12: t = − 3.8, p < 1E-03). However, the results on group differences did not change when including ICV as an additional covariate in the GLM. Ventricular volume showed an association with IC1,2,5,8 (IC1: t = − 6.3, p < 1E-03, IC2: t = − 2.8, p = 5E-03, IC5: t = 2.5, p = 1.3E-02, IC8: t = − 8.6, p < 1E-03, uncorrected), accounting for age, sex and ICV, reflecting larger ventricular GMD regions in the patient groups. Although IC1 and IC2 were mainly dominated by surface area and thickness (Fig. 1), they were also subtly influenced by the GMD measure (Supplemental Fig. 10).

In general, weak correlations were observed between the components (maximum partial correlation = 0.26, full correlation = 0.23), and only a small subset of components was significantly correlated with at least one other component (IC2 with IC5,8,9,20,26; IC4 with IC40; IC0 with IC1,5) (Supplementary Fig. 6).

Associations with cognitive domain scores, PGRS and symptom domain scores

Groups were significantly different in all cognitive domains with lower performance in patients (schizophrenia < bipolar disorder < controls). Bipolar disorder had higher PGRS for bipolar disorder than controls, and both schizophrenia and bipolar disorder had higher PGRS for schizophrenia than controls. PGRS for schizophrenia and bipolar disorder explained ~ 13% and 6% of the variance in diagnostic status (Nagelkerke pseudo R2), respectively (Table 1).

Permutation testing revealed no significant univariate associations between the brain morphology patterns and any of the cognitive, PGRS, or symptom scores. However, our results showed several nominally significant associations with cognitive scores and PGRS, which warrant replication analyses in independent samples. Among others, processing speed showed a positive association with IC8 in schizophrenia (t = 3.31, p = 1.1E-03, uncorrected) and in bipolar disorder (t = 2.1, p = 3.73E-02, uncorrected), indicating that reduced processing speed may be related to reduced lateral temporal and medial frontal cortical thickness and specific local GMD regions in the patient groups (Fig. 4). IC1, reflecting global surface area, and IC5, reflecting cerebellar GDM, also showed a trend association with processing speed and working memory/attention, respectively (see SI for more details and other results).

Fig. 4

Scatter plot of processing speed as a function of IC8's subject weights. Linear effects of age and sex on both IC8 and processing speed were regressed out before plotting. SZ = Schizophrenia, BD = Bipolar Disorder, HC = Controls.

Post-hoc analyses with medication, substance use, education

Post-hoc analyses on the set of six components with significant diagnosis effect showed no significant associations with clinical variables, or with education, within the patient groups.

Multivariate machine learning classification analyses

Classification of schizophrenia and controls using all LICA components yielded an AUC of 0.75 (Fig. 2) (accuracy = 69%, specificity = 74%, sensitivity = 62%). Classification of bipolar disorder and controls, and schizophrenia and bipolar disorder yielded lower performance (bipolar disorder vs. controls: AUC = 0.67, accuracy = 66%, specificity = 72%, sensitivity = 58%; schizophrenia vs. bipolar disorder: AUC = 0.59, accuracy = 58%, specificity = 52% and sensitivity = 62%). Fig. 2B showed that IC1,2,5,8,12,25 are the most important features for both schizophrenia vs. controls and bipolar disorder vs. controls, and IC1 and IC39 (reflecting hippocampal GMD) for schizophrenia vs. bipolar disorder. This performance remained unchanged when classifying only subjects with schizophrenia diagnosis (excluding those with schizoaffective and schizophreniform), bipolar disorder I and controls (AUC = 0.76, 0.67 and 0.53 for schizophrenia vs. controls, bipolar I vs. controls, and schizophrenia vs. bipolar I, respectively).

Fig. 2

Performance of the pairwise group classification on the whole sample: A. ROC plot; B. Spider plot indicating the relative importance of all features included (1 = most important, 0 = least important). SZ = Schizophrenia, BD = Bipolar Disorder, HC = Controls.

Further, we evaluated the complementary predictive value of brain patterns, cognitive scores and PGRS (Fig. 3). In classification of schizophrenia and controls, cognitive scores gave a higher performance than LICA features (AUC = 0.81 vs. 0.67). Using the combined set of cognitive and LICA features yielded significantly higher results (AUCcognitive + LICA = 0.85, p = 7.4E-03, permutation testing keeping the original ordering of cognitive features while permuting randomly LICA features, repeated 5000 times) compared to using cognitive features alone (AUCcognitive only = 0.81) or using LICA alone (AUCLICAonly = 0.67, p = 2E-04, permutation testing keeping the original ordering of LICA features while randomly permuting cognitive features). PGRS when used alone yielded a modest performance (AUC = 0.6), but when used in addition to either LICA features or cognitive scores also led to an improved performance (AUC = 0.73, p = 2E-04; 0.84, p = 1E-02, respectively). However, additionally including PGRS to the combination of LICA and cognitive scores did not improve performance. The relative feature importance plot (Fig. 3D) showed that the cognitive features are more informative than LICA and PGRS features. IC1,2 and 5 are the most informative features among the LICA components. Similar patterns were observed for bipolar disorder vs. controls, although accuracies were generally lower. Classification between schizophrenia and bipolar disorder showed a low performance even when using all feature sets together (AUC < 0.63).

Fig. 3

Performance of the pairwise group classification performed on a subset where imaging, cognitive and polygenic risk score are available: A,B,C. ROC plots for classification schizophrenia (SZ) and controls (HC), bipolar disorder (BD) and HC, SZ and BD, respectively; D. Spider plot indicating the relative importance in case all feature sets were included (1 = most important, 0 = least important). PS = Processing speed, VL/mem = Verbal learning/memory, EF = Executive function, WM/att = Working memory/attention.

Discussion

Using a data-driven multivariate approach for multimodal fusion of brain morphometric properties, we identified six novel distinct multivariate patterns of complex brain variability, showing significant differences in schizophrenia and bipolar disorder compared to controls. Two independent patterns reflected anatomically distributed cortical thickness and surface area, suggesting global cortical thinning in the patient groups, and area reduction in schizophrenia compared to bipolar disorder and controls. Beyond these global effects, four multimodal regional brain patterns also showed strong effects, capturing 1) correlated increases in cingulate thickness and decreases in white matter volume, 2) fronto-parietal and temporal thickness, 3) cerebellar gray matter density, and 4) a distinct anterior-posterior gradient in cortical thickness. Group classification revealed high accuracy for diagnostic prediction using cognitive data alone. Adding polygenic risk and imaging features significantly increased accuracy, suggesting complementary predictive value of brain imaging, cognitive performance and polygenic risk in classifying between patients with severe mental illness and healthy controls.

LICA allowed us to isolate a relatively small number of distinct and biologically interpretable patterns of brain structure variability, each of which may explain just a small portion of the total variation in brain phenotypes underlying severe mental disorders. One important aspect of this approach is that it is fully data-driven and unbiased, i.e. no other information (age, sex or diagnosis) than the imaging data was used to derive the patterns. Furthermore, by means of simultaneous modeling of both global and regional features linking complementary morphological measures (thickness, surface area, GMD maps), we were able to untangle global from regional effects that would have been difficult to disentangle when considering each of the measures separately using univariate approaches. Such a multivariate approach focuses on interrelated patterns across multimodal measures, allowing the identification of complex and weak effects hidden in high-dimensional data not previously shown in studies applying standard univariate methods (Calhoun and Sui, 2016, Douaud et al., 2014).

Case-control difference in LICA subject weights

We identified six distinct multivariate patterns of brain structure variations showing significant group differences. IC1 reflected global variability in surface area, and the strong correlation with ICV supports its global properties. Schizophrenia showed smaller subject weights than bipolar disorder and controls, indicating that this component is specific to schizophrenia. This finding is also consistent with recent meta-analyses confirming reduced ICV in schizophrenia (van Erp et al., 2016a) but not in bipolar disorder (Hibar et al., 2016). Compared to controls, we also observed reduced subject weights in schizophrenia and bipolar disorder in a component reflecting global cortical thickness (IC2), in line with evidence of widespread cortical thinning in schizophrenia (Goldman et al., 2009, Kuperberg et al., 2003, Rimol et al., 2010a), and also thinning in bipolar disorder (Elvsåshagen et al., 2013, Lyoo et al., 2006, Rimol et al., 2010a), although to a lesser extent.

Since IC1 and IC2 modeled the global effects, which accounted for a large amount of data variance, subtler and more anatomically specific variance potentially related to disease processes can be revealed. We identified four other components (each explaining < 2% variance) showing strong main effects of group. Each of these components potentially represents independent pathogenic mechanisms, as supported by the weak pairwise correlation (Supplementary Fig. 6).

IC5 revealed a distinct pattern of distributed GMD and WMD voxels, reflecting decreasing WMD with increasing GMD, in addition to localized and positively weighted lateral fronto-parietal and cingulate thickness variation. IC5 is similar to a component recently suggested to reflect a network of gray matter regions involved in transmodal processing (Douaud et al., 2014), with estimated age-trajectories that mirrored both healthy developmental and aging processes, and with heightened vulnerability to schizophrenia and Alzheimer's disease by spatially recapitulating the pattern of macrostructural abnormalities seen in both disorders. In line with this recent study, we observed increased subject weights in schizophrenia, indicating a correlated pattern of increased cortical GMD, increased cortical thickness along the cingulate gyrus, along with decreased WMD. This specific component captured a source of variation over and above global effects, and the focal pattern of increased GMD and cortical thickness was identifiable since the variability related to global cortical thickness was accounted for in IC2.

The brain pattern observed in IC5 is in line with the balloon model hypothesis (Harasty et al., 2003, Seldon, 2005), which proposes that as the WM grows, the cortical columns stretch and become thinner. This tangential cortical expansion caused by WM growth is hypothesized to increase the capacity of the cortex to differentiate afferent signals (Seldon, 2005). In accordance with this, the increased subject weights in schizophrenia and bipolar disorder compared to controls may reflect decreased WM volume along with decreased cortical column differentiation as captured in increased thickness and GMD.

IC8 implicated the lateral temporal, parietal and medial frontal thickness along with GMD map. The GMD pattern largely comprises temporal and occipito-parietal regions, including the occipital fusiform, lingual cortex, precentral and postcentral gyrus and hippocampus. Schizophrenia and bipolar disorder showed significantly reduced subject weights compared to controls, indicating reduced cortical thickness and GMD in the implicated regions. Our GMD findings in this component are consistent with several previous studies reporting GMD reduction in schizophrenia at the precentral, postcentral, parahippocampus gyrus (Glahn et al., 2008, Gupta et al., 2015), fusiform (Onitsuka et al., 2003, van Erp et al., 2016b) and lingual cortex (Antonova et al., 2005). In our study, all of these structures are implicated in the same independent component, suggesting coordinated volumetric alterations in the patient groups. The thickness pattern reflected in this component is in line with previous reports on cortical thinning in schizophrenia and bipolar disorder (Kuperberg et al., 2003, Lyoo et al., 2006, Rimol et al., 2010a). Beyond findings in univariate studies, IC8 further revealed co-occurring alterations of thickness and GMD at partially non-overlapping brain regions, which may be influenced by the same aspects of the etiology in mental disorder. Furthermore, association results with cognitive scores, although only nominally significant, revealed a positive correlation between IC8 and processing speed (stronger in schizophrenia than in bipolar disorder), suggesting that regional cortical thinning and GMD reduction is associated with reduced processing speed in patients.

IC12 revealed a distinct anterior-posterior gradient in cortical thickness with increasing frontal thickness along with decreasing posterior thickness, thus reflecting the relative thickness increase along the posterior-anterior axis over and above all other components. The observed decreased subject weights in schizophrenia and bipolar disorder compared to controls reflects a shift of the thickness distribution towards a stronger posterior compared to anterior weighting. The decreased frontal distribution corresponds with the characteristic fronto-temporal cortical thinning in schizophrenia and bipolar disorder (Elvsåshagen et al., 2013, Kuperberg et al., 2003, Rimol et al., 2012), but extends previous findings by suggesting that this pattern of effects reflects a shift in the anterior-posterior distribution of cortical thickness, over and above global thinning. The linear combination of IC12 and the global cortical thinning component (IC2), would likely result in a net effect of frontal-temporal thinning, as obtained using mass-univariate analyses on the observed absolute thickness values (SI, Supplementary Fig. 9). An important strength of LICA was the ability to model the global thickness reduction effect into a separate independent component (IC2), and thus more subtle patterns of relative thickness changes as reflected by IC12 can be revealed. The spatial distribution of IC12, which was derived solely from the imaging data in a data-driven manner, is consistent with patterns implicating anterior-posterior gradients of gene expression shown in animal studies (Kudo et al., 2007, O'Leary et al., 2007), and partially overlaps with a gradient pattern of cortical thickness genetic correlation revealed in a seed-based analysis seeding from the primary visual cortex (Rimol et al., 2010b). This pattern is also in line with a genetically informed two-cluster solution of cortical morphology, although more closely resembling an anterior-posterior arealization pattern than a ventral-dorsal cortical thickness gradient (Chen et al., 2013).

IC25 reflected a specific cerebellar pattern. Group analyses revealed reduced subject weights in schizophrenia compared to controls indicating reduced cerebellar GMD in this patient group. While the cerebellum has not been among the brain regions most consistently reported in studies of brain structure in schizophrenia (Honea et al., 2005, Shenton et al., 2001), this subcortical structure is increasingly seen as relevant to several brain disorders with marked cognitive changes and distinct patterns of cerebro-cortical degeneration, such as Alzheimer's disease and frontotemporal dementia (Guo et al., 2016). Moreover, both structural (Bostan et al., 2013, Palesi et al., 2015) and functional (Buckner et al., 2011) data suggest extensive connectivity between the cerebellum and associative cerebral cortex. Importantly, the current well-powered study supports previous reports of altered cerebellar structure in schizophrenia (Okugawa et al., 2007) and bipolar disorder (Johnson et al., 2015, Laidi et al., 2015). Our results for IC25 are also in line with two recent large-scale VBM-based analyses of schizophrenia (Gupta et al., 2015, Moberget et al., 2017), underscoring the relevance of the cerebellum for the pathophysiology of severe mental disorders (Watson et al., 2014). The fact that IC25 is specific to the cerebellum, despite that LICA is able to probe co-occurring variation across multiple measures, and the low correlation between this component and other components, may support the notion of a complex interaction between this structure and other brain regions (Andreasen and Pierson, 2008, Barch, 2014). However, a recent study identified robust associations between regional cerebellar volumes and cortical thickness, in particular in patients with schizophrenia(Moberget et al., 2017). Interestingly, the observed patterns cortical thickness associations with cerebellar volumes mirrored the case-control differences in cortical thickness, suggesting coordinated cortico-cerebellar reductions in schizophrenia (Moberget et al., 2017).

Ventricle volume showed an association with IC1,2,5 and 8, reflecting increased ventricular GMD in the patient groups. Ventricle enlargement is among the most robust brain structural findings in schizophrenia (Moberget et al., 2017, van Erp et al., 2016a). Our results showed that ventricular volumetric alteration is not implicated in one single independent mode of variability, but rather distributed and co-occurring with global alteration pattern of surface area and thickness in schizophrenia.

Associations with cognitive domain and polygenic risk

Whereas the rationale for a close link between the brain structure and inter-individual differences in performance on cognitive tests is sound (Kanai and Rees, 2011), empirical results have been mixed. Previous studies have reported either null association (Killgore et al., 2009) or mixed results (Van Petten, 2004), and have proven difficult to replicate (Boekel et al., 2015). Our results suggest moderate yet potentially interesting associations between the derived multivariate brain morphology patterns and neuropsychological performance. For instance, we found a nominal positive association between IC8 and processing speed within the patient groups. Previous studies have reported an association between processing speed and white matter in schizophrenia involving occipital regions (Antonova et al., 2005, Karbasforoushan et al., 2015), as well as gray matter volume in several frontal, parietal and occipital regions in healthy subjects (Chee et al., 2009). Our results may complement these findings by suggesting that reduced processing speed is associated with a coordinated pattern of reduced medial frontal, lateral temporal, parietal thickness and regional GMD, comprising the precentral gyrus, postcentral gyrus, hippocampus and occipital fusiform, lingual cortices. In accordance with previous findings on a positive correlation between processing speed and total brain volume (Betjemann et al., 2010), the global surface area component (IC1), strongly correlated with ICV, showed a positive association with processing speed. Our results also showed a nominal association between the cerebellar component (IC5) and working memory/attention, in line with studies indicating that the cerebellum underwent rapid phylogenetic increases in size compared to the neocortex (Barton and Venditti, 2013), lesion studies demonstrating that cerebellar alterations produce deficits in verbal working memory (Ravizza et al., 2006) and VBM studies reporting positive correlations between GM volume in cerebellar subregions and working memory (Ding et al., 2012) and IQ (Moberget et al., 2017). It should be noted however that the current associations were non-significant after stringent multiple testing correction and thus replications in future studies are warranted.

Our results revealed no robust associations between brain structure and schizophrenia or bipolar disorder PGRS. This may partly be explained by the fact that PGRS is a cumulative composite score (Psychiatric GWAS Consortium Bipolar Disorder Working Group, 2011, Schizophrenia Psychiatric Genome-Wide Association Study Consortium, 2011) which only explained a moderate amount of variance in diagnostic status in the current sample, and further refinements may increase their sensitivity to brain variability. Although some evidence suggests common genetic variance in brain anatomy and schizophrenia (Lee et al., 2016), our results correspond with a recent large-scale study which reported no evidence for genetic overlaps between schizophrenia risk and volumes of several brain structures (Franke et al., 2016). It is therefore conceivable that the genetic variants influencing individual differences in brain structure is only weakly overlapping with the variants modulating risk for severe mental illness, and thus the hidden heritability of severe mental disorders needs to be revealed elsewhere.

Classification analysis

Pairwise group classifications using LICA features

To gain further insights on the multivariate combinations of the LICA patterns as well as the sensitivity of such combinations to diagnosis at an individual level, we submitted the LICA components' subject weights to multivariate classification using a cross-validated approach. Classification analysis revealed moderate predictive value for the schizophrenia versus controls in a two-class classification. This performance is comparable to previous findings in a sample with equivalent size using VBM and support vector machine (Nieuwenhuis et al., 2012). In this study, the authors applied VBM and feature selection to use the most predictive sets of voxels for optimal group classification. Our study has a slightly different focus. With the primary aim on finding covariance patterns across complementary brain morphological measures, we applied the LICA approach for feature extraction and applied machine learning classification as a means to further evaluate the resulting features in terms of their discriminative power. Classification of schizophrenia from controls showed slightly lower performance compared to other studies (Kambeitz et al., 2015). There are however several factors, including sample characteristics, that may influence the classification performance. In this study, we originally investigated the entire spectrum of schizophrenia and bipolar disorder, which may show increased heterogeneity within the patient groups and phenotypic overlap with the control group than e.g. schizophrenia and bipolar I disorder. The classification performance obtained was very similar when we excluded patients with schizophreniform and schizoaffective disorders. However, the latter two subgroups comprised a very small number of subjects compared to schizophrenia. The heterogeneity within schizophrenia spectrum remains to be explored in future studies with more reasonably sized subgroups.

Compared to the classification between schizophrenia and controls, the performance of the two-class classifications of bipolar disorder was less accurate versus controls and rather low versus schizophrenia. Schnack et al. reported substantially higher accuracy for classification of schizophrenia versus bipolar disorder (88%) and lower for bipolar disorder versus controls (59%) (Schnack et al., 2014). Several factors may explain the differences in performance. In addition to the difference in sample size, where our sample was significantly larger and included a wide range of schizophrenia and bipolar patients, the schizophrenia group included in that study showed longer duration of illness and higher symptom severity, as reflected by the higher PANSS scores. Furthermore, most of the bipolar disorder patients received lithium treatment (68%) compared to 18% in our sample. These differences may in turn imply that the brain phenotypes of our bipolar group showed lower and higher discriminative power compared to schizophrenia and controls, respectively, as compared to those used in (Schnack et al., 2014), leading to lower classification performance of bipolar disorder versus schizophrenia and higher for bipolar disorder versus controls. This is in accordance with a recent systematic assessment on the relation between sample size and classification in psychiatric disorders, where the authors reported that studies with larger sample size are more likely to show lower accuracy than those with smaller sample sizes, possibly attributed to the increased clinical heterogeneity of the sample, but studies with smaller sample size may show lower generalizability to other samples (Schnack and Kahn, 2016). The difference in classification algorithms used, random forest versus support vector machine, may also have an influence, although most likely to a lesser extent than the patient clinical profiles as the utility of both classifiers in clinical studies has been validated (Greenstein et al., 2012, Kambeitz et al., 2015, Schnack et al., 2014). The difference in feature extraction approaches may also partly explain the difference in performance.

Other structural MRI-based studies documented in a recent review (Wolfers et al., 2015) reported comparable performance with our results on classification of bipolar disorder and controls (Serpa et al., 2014). A higher performance was reported by Rocha-Rego et al. (Rocha-Rego et al., 2014) (72% and 73%) and Bansal et al. (Bansal et al., 2012) (98.2%). However, the considerable differences in sample size (much smaller number of patients (N ≤ 26 versus N = 190 in our study)) and sample characteristics (e.g. older patients (Bansal et al., 2012, Rocha-Rego et al., 2014) with our study showing larger sample heterogeneity likely explain the difference in performance. In general, although some previous studies have reported relatively high classification accuracies, our results of a low classification performance between the schizophrenia and bipolar groups are in line with the reported clinical, genetic, and neurobiological overlaps between these disorders (Cross-Disorder Group of the Psychiatric Genomics Consortium et al., 2013, Karege et al., 2012).

Our results documented moderate classification performance, which is accordance with the notion that MRI-based machine learning for prediction of diagnosis is still in its infancy and has not yet reached the level of accuracy required in a clinical setting (Wolfers et al., 2015). Nevertheless, the multivariate structural patterns sensitive to diagnosis as well as the added discriminative value they carried, as reported in our study, motivate further research on the extraction of novel brain phenotypes using multivariate multimodal fusion techniques to improve the prediction performance.

Assessing added discriminative values of LICA, cognitive and polygenic risk features

The cognitive scores showed a higher classification performance than LICA features, suggesting that cognitive performance has higher clinical sensitivity than imaging features, which is in line with the notion of an intimate link between cognitive function and the underlying pathophysiology (Keefe and Harvey, 2012). Adding cognitive to LICA features and vice versa, LICA to cognitive features, yielded a significantly improved performance, indicating that the two feature sets contain complementary predictive information. This increased performance may reflect added specificity provided by the brain features, a hypothesis that should be addressed in future studies including data from patients with a range of different brain disorders. Whereas the neurocognitive test scores may be non-specific to biological pathways underlying psychiatric diseases, neuroimaging may yield novel insights in the neural mechanisms underlying psychiatric symptoms to inform new interventions (Linden, 2012). Our results support that combining complementary feature sets from different domains improves prediction sensitivity. Similarly, PGRS showed complementary predictive value together with LICA or cognitive features. Adding PGRS to the combined sets of LICA and cognitive features did however not improve performance implying that the predictive value of PGRS is captured by a multivariate combination of brain and cognitive variables. For classification between schizophrenia and bipolar disorder, adding the LICA features did not improve the performance, indicating that the multivariate morphological features, as derived using LICA, did not provide added predictive value. More research using different approaches such as partial least squares and on cohorts with different clinical/cognitive profiles is needed to assess the added value of the brain morphological phenotype in distinguishing between these two disorders.

Limitations

LICA performs joint analysis of multiple measures and relies on the assumption that the same subject weights are shared across measures. Alternatively, one could consider running conventional ICA separately for each group and subsequently perform post-hoc analyses. However, this would not allow for a joint modeling of the shared variance across features, and it is not clear how one would couple the components (Groves et al., 2011). LICA has the ability to automatically balance information across measures, allowing a component to be only driven by one measure or jointly driven by multiple measures depending on how strongly related the information of a measure to the underlying pattern of the component. It is possible that LICA would return unimodal components that resemble components obtained using separate ICA results. In this case, LICA is still preferred over separate ICA decompositions since its joint analysis framework allows for better interpretation in terms of specificity. In particular, such unimodal components, when found together with other multimodal components in the LICA context, i.e. those driven by multiple measures, may reflect disease or environmental mechanism that is specific to the relevant measure. This is in principle different from post-hoc correlation analysis of the components resulting from separate ICAs. Nevertheless, like with other multivariate methods, the results obtained with LICA should be interpreted with caution, keeping in mind the inherent limitations of the technique, including the assumption of the similar spatial maps across groups and of shared subject weights across measures (Groves et al., 2011).

We included all three groups in the same decomposition. One could consider performing separate analyses for schizophrenia vs. controls, and bipolar disorder vs. controls, and subsequently compare the two patient groups by means of post-hoc analyses. However, such a comparison requires that the components involved are the same across group pairs, which may not be possible to ensure. To compare between schizophrenia and bipolar disorder in relation to controls, which would allow us to assess the specificity of the patterns, one approach is to include all groups in the same decomposition. Unlike supervised techniques such as partial least squares, which maximize the covariance between the input variables and diagnosis label, LICA is a fully data-driven technique that makes use of no diagnosis information, and decomposes the data into independent modes of brain variability. The resulting pattern may reflect different biological and environmental factors, or a combination of them. Whereas data-driven analysis has the potential to reveal novel patterns best describing the data structure that could be sensitive to diagnosis effect, given the high overlap between schizophrenia and bipolar disorder, it may be more sensitive to treat them as separate groups and search for patterns that best differentiate between the two disorders in a supervised analysis framework.

We consistently used random forest as a tool for validating the discriminative power of the LICA features and their added value when used in combination with cognitive and PGRS scores. Whereas the classification performance was in line with the literature, this classifier may however not represent the optimal choice across all feature sets, particularly the PGRS feature set with only two features included.

Lastly, we included a large sample size covering the entire spectrum of schizophrenia and bipolar disorder, since we aimed to find morphological patterns generalized to the psychosis spectrum instead of specific to a subgroup. Whereas the large sample size is a strength, the clinical heterogeneity of the sample may likely partly explain the relatively low classification performance (Schnack and Kahn, 2016).

Future directions

A possible extension of this work could be to include the genetic data and/or cognitive domain scores as additional “modalities” in LICA decomposition, which may present an effective way to model common variances between genetic, cognitive data and brain morphometry. Furthermore, using multivariate techniques such as LICA to combine structural MRI, functional MRI and diffusion MRI may reveal novel sensitive imaging markers characterizing the potentially complicated relationships among these modalities. Any analysis making use of the diagnosis information may be limited by the potential caveat of the current diagnosis. Another future direction would be to utilize unsupervised clustering or machine learning techniques to identify subgroups that may be more sensitive to clinical or cognitive profiles.

Conclusions

Using a data-driven approach for multimodal fusion of brain morphometric features, we identified six distinct multivariate brain patterns underlying strong effects of schizophrenia and bipolar disorder. We provided evidence of complementary diagnostic predictive value of brain imaging, cognitive and genetic features with cognitive performance measures being most sensitive to classification of patients. Although the classification accuracy does not support a direct clinical utility, the novel multivariate brain morphology patterns are biologically interpretable and may inform models of pathogenic mechanisms in severe mental illness.

Fundings

The study has received funding from the European Commission's Seventh Framework Programme (FP7/2007-2013, #602450, IMAGEMEND), Research Council of Norway (#213837, #223273, #204966/F20), the South-Eastern Norway Regional Health Authority (#2013123, #2014097, #2016083) and Kristian Gerhard Jebsen Foundation.