First, to evaluate PSD symptom variation, we examined core PSD psychopathology metrics captured by two instruments: the Brief Assessment of Cognition in Schizophrenia (BACS) and the Positive and Negative Syndrome Scale (PANSS). We refer to these 36 items as ‘symptom measures’ throughout the rest of the paper. We observed group mean differences across DSM diagnoses (Figure 1A, p<0.05 Bonferroni corrected); however, symptom measure distributions revealed notable overlap that crossed diagnostic boundaries (Tamminga et al., 2014; Keshavan et al., 2011). Furthermore, we observed marked collinearity between symptom measures across the PSD sample (Figure 1B), indicating that a dimensionality-reduced solution may sufficiently capture meaningful varation in this symptom space. We hypothesized that such a dimensionality-reduced symptom solution may improve PSD symptom-neural mapping as compared to pre-existing symptom scales.

Figure 1

Download asset Open asset

Here, we report results from a principal component analysis (PCA) as it produces a deterministic solution with orthogonal axes (i.e. no a priori number of factors needs to be specified) and explains all variance in symptom measures. Results were highly consistent with prior symptom-reduction studies in PSD: we identified five PCs (Figure 1C), which captured ∼50.93% of all variance (see Materials and methods and Appendix 1—figure 2; Chen et al., 2020).

The key innovation here is the combined analysis across PSD diagnoses of core PSD symptoms and cognitive deficits, a fundamental PSD feature (Barch et al., 2013). The five PCs revealed few distinct boundaries between DSM categories (Figure 1D). Notably, for PCs 2-5 we observed substantial overlap in PC scores between controls and patients for all DSM diagnostic groups. This overlap is not unexpected, given that behaviors measured by the BACS and PANSS (e.g. mood, cognition) occur to some extent in ‘healthy’ individuals not meeting DSM diagnostic criteria (Verdoux and van Os, 2002; Nuevo et al., 2012; Kelleher and Cannon, 2011; Stefanis et al., 2002). In contrast, PC1 showed marked differentiation between PSD and CON, reflecting global functioning which was reduced across DSM diagnostic groups.

Figure 1E highlights loading configurations of symptom measures forming each PC. To aid interpretation, we assigned a name for each PC based on its most strongly weighted symptom measures. This naming is qualitative but informed by the pattern of loadings of the original 36 symptom measures. For example, PC1 was highly consistent with a general impairment dimension (i.e. ‘Global Functioning’); PC2 reflected variation in cognition (i.e. ‘Cognitive Functioning’); PC3 indexed a complex configuration of psychosis-spectrum relevant items (i.e. ‘Psychosis Configuration’); PC4 captured variation in mood and anxiety related items (i.e. ‘Affective Valence’); finally, PC5 reflected variation in arousal and level of excitement (i.e. ‘Agitation/Excitation’). For instance, a generally impaired patient would have a highly negative PC1 score, which would reflect low performance on cognition and elevated scores on most other symptoms. Conversely, an individual with a high positive PC3 score would exhibit delusional, grandiose, and/or hallucinatory behavior, whereas a person with a negative PC3 score would exhibit motor retardation, social avoidance, and possibly a withdrawn emotional state with blunted affect (Gelenberg, 1976). Comprehensive loadings for all five PCs are shown in Appendix 1—figure 1. Figure 1F highlights the mean of each of the three diagnostic groups (colored spheres) and healthy controls (black sphere) projected into a three-dimensional orthogonal coordinate system for PCs 1,2 and 3 (x,y,z axes respectively; alternative views of the three-dimensional coordinate system with all patients projected are shown in Appendix 1—figure 1). Critically, PC axes were not parallel with traditional aggregate symptom scales. For instance, PC3 is angled at ∼45° to the dominant direction of PANSS Positive and Negative symptom variation (purple and blue arrows respectively in Figure 1F).

PC3 loads most strongly onto hallmark symptoms of PSD (including strong positive loadings onto PANSS Positive symptom and strong negative loadings onto most PANSS Negative symptoms). Therefore, we focus on PC3 as an opportunity to quantify a fully data-driven dimension of symptom variation that is highly characteristic of the PSD patient population. Additionally, this bi-directional symptom axis captured shared variance from additional symptoms, such as PANSS General items and cognition. PC3 provides an empirical demonstration that using a data-driven dimensionality-reduced solution (via PCA) can reveal novel symptom patterns underlying PSD psychopathology.

Notably, independent component analysis (ICA), an alternative dimensionality reduction procedure which does not enforce component orthogonality, produced similar effects for this PSD sample (see Appendix 1 - Note 1 and Appendix 1—figure 3A). Certain pairs of components between the PCA and ICA solutions appear to be highly similar and directly comparable (IC5 and PC4; IC4 and PC5) (Appendix 1—figure 3B). On the other hand, PCs 1–3 and ICs 1–3 do not exhibit a one-to-one mapping. For example, PC3 appears to correlate positively with IC2 and equally strongly negatively with IC3, suggesting that these two ICs are oblique to the PC and perhaps reflect symptom variation that is explained by a single PC. The orthogonality of the PCA solution forces the resulting components to capture maximally separated, unique symptom variance, which in turn map robustly on to unique neural circuits. We observed that the data may be distributed in such a way that highly correlated independent components emerge in the ICA, which do not maximally separate the symptom variance associated with neural variance. We demonstrate this by plotting the relationship between parcel beta coefficients for the ${\beta }_{PC3}GBC$ map versus the ${\beta }_{IC2}GBC$ and ${\beta }_{IC3}GBC$ maps (Appendix 1—figure 11G).

Next, we show that the PCA solution was highly stable when tested across sites, and k-fold cross-validations, and was reproducible in independent split-half samples. First, we show that the symptom-derived PCA solution remained highly robust across all sites (Figure 2A) and 5-fold cross-validation iterations (see Materials and methods). The total proportion of variance explained as well as the total number of significant PCs remained stable (Figure 2A). Second, PCA loadings accurately and reliably computed the scores of subjects in a hold-out sample (Figure 2B). Specifically, for each left-out site, we tested if single-subject predicted scores were similar to the originally observed scores of these same subjects in the full PCA (i.e. performed with all 436 subjects). This similarity was high for all sites (Appendix 1—figure 4). Also, we verified that the predicted-versus-observed PCA scores also remained similar via a k-fold cross-validation (for k = 2 to 10, Figure 2C and Appendix 1—figure 3). Finally, across all five PCs, predicted-to-observed similarity of PCA loadings was very high. Moreover, PCA loadings were stable when testing via leave-site-out analysis and 5-fold cross-validation. We furthermore demonstrated the reproducibility of the solution using 1000 independent split-half replications (Figure 2D). For each run of split-half replication, PSD subjects were split into two samples and a PCA was performed independently for each sample. The loadings from both PCAs were then compared using a Pearson’s correlation. Importantly, results were not driven by medication status or dosage (Appendix 1—figure 6). Collectively, these data reduction analyses strongly support a stable and reproducible low-rank PSD symptom geometry.

Figure 2

Download asset Open asset

Next, we tested if the dimensionality-reduced symptom geometry can identify robust and novel patterns of neural variation across PSD patients. We opted to use global brain connectivity (GBC), a summary FC metric, to measure neural covariance because it yields a parsimonious measure reflecting how globally coupled an area is to the rest of the brain (Cole et al., 2010; see Materials and methods). Furthermore, we selected GBC because: (i) the metric is agnostic regarding the location of dysconnectivity as it weights each area equally; (ii) it yields an interpretable dimensionality-reduction of the full FC matrix; (iii) unlike the full FC matrix or other abstracted measures, GBC produces a neural map, which can be related to other independent neural maps (e.g. gene expression or pharmacology maps, discussed below). Furthermore, GBC has been shown to be sensitive to altered patterns of global neural connectivity in PSD cohorts (Anticevic et al., 2013; Fornito et al., 2011; Hahamy et al., 2014; Cole et al., 2011) as well as in models of psychosis (Preller et al., 2018; Driesen et al., 2013).

All five PCs captured unique patterns of GBC variation across the PSD (Appendix 1—figure 9), which were not observed in CON (Appendix 1—figure 10). Again we highlight the hallmark ‘Psychosis Configuration’ dimension (i.e. PC3), here to illustrate the benefit of the low-rank PSD symptom geometry for symptom-neural mapping relative to traditional aggregate PANSS symptom scales. The relationship between total PANSS Positive scores and GBC across N = 436 PSD patients (${\beta }_{Pos}GBC$, Figure 3A) was statistically modest (Figure 3B) and no areas survived whole-brain type-I error protection (Figure 3C, p<0.05). In contrast, regressing PC3 scores onto GBC across N = 436 patients revealed a robust symptom-neural ${\beta }_{PC3}GBC$ map (Figure 3E–F), which survived whole-brain type-I error protection. Of note, the PC3 ‘Psychosis Configuration’ axis is bi-directional whereby individuals who score either highly positively or negatively are symptomatic. Therefore, a high positive PC3 score was associated with both reduced GBC across insular and superior dorsal cingulate cortices, thalamus, and anterior cerebellum and elevated GBC across precuneus, medial prefrontal, inferior parietal, superior temporal cortices, and posterior lateral cerebellum – consistent with the default-mode network (Fox et al., 2005). A high negative PC3 score would exhibit the opposite pattern. Critically, this robust symptom-neural mapping emerged despite no differences between mean diagnostic group PC3 scores (Figure 3D). These two diverging directions may be captured separately in the ICA solution, when orthogonality is not enforced (Appendix 1—figure 12). Moreover, the PC3 symptom-neural map exhibited improved statistical properties relative to other GBC maps computed from traditional aggregate PANSS symptom scales (Figure 3G–J). Notably, the PC2 – Cognitive Functioning dimension, which captured a substantial proportion of cognitive performance-related symptom variance independent of other symptom axes, revealed a circuit that was moderately (anti-)correlated with other PC maps but strongly anti-correlated with the BACS composite cognitive deficit map (r = −0.81, Appendix 1—figure 8O). This implies that the PC2 map reflects unique neural circuit variance that is relevant for cognition, independent of the other PC symptom dimensions.

Figure 3

Download asset Open asset

After observing improved symptom-neural PC3 statistics, we tested if these neural maps replicate. Recent attempts to derive stable symptom-neural mapping using multivariate techniques, while inferentially valid, have not replicated (Dinga et al., 2019). This fundamentally relates to the tradeoff between the sample size needed to resolve multivariate neurobehavioral solutions and the size of the feature space. To mitigate the feature size issue we re-computed the ${\beta }_{PC}GBC$ maps using a functionally-derived whole-brain parcellation via the recently-validated CAB-NP atlas (Glasser et al., 2016; Ji et al., 2019b; Materials and methods). Here, a functional parcellation is a principled way of neural feature reduction (to 718 parcels) that can also appreciably boost signal-to-noise (Glasser et al., 2016; Ji et al., 2019b). Indeed, parcellating the full-resolution ‘dense’ resting-state signal for each subject prior to computing GBC statistically improved the group-level symptom-neural maps compared to parcellating after computing GBC (Figure 4A–C, all maps in Appendix 1—figure 9). Results demonstrate that the univariate symptom-neural mapping was indeed stable across fivefold bootstrapping, leave-site-out, and split-half cross-validations (Figure 4D–N, see Appendix 1 - Note 2), yielding consistent symptom-neural PC3 maps. Importantly, the symptom-neural maps computed using ICA showed comparable results (Appendix 1—figure 12).

Figure 4

Download asset Open asset

Several recent studies have reported ‘latent’ neurobehavioral relationships using multivariate statistics (Xia et al., 2018; Drysdale et al., 2017; Yu et al., 2019), which would be preferable because they simultaneously solve for maximal covariation across neural and behavioral features. Given the possibility of deriving a stable multivariate effect, here we tested if results improve with canonical correlation analysis (CCA) (Hardoon et al., 2004; Figure 5A).

Figure 5

Download asset Open asset

To evaluate if the number of neural features affects the solution, we computed CCA using GBC from: (i) 180 symmetrized cortical parcels; (ii) 359 bilateral subcortex-only parcels; (iii) 192 symmetrized subcortical parcels; (iv) 12 functional networks, all from the brain-wide cortico-subcortical CAB-NP parcellation (Glasser et al., 2016; Ji et al., 2019b). We did not compute a solution using 718 bilateral whole-brain parcels, as this exceeded the number of samples and rendered the CCA insolvable. Notably, the 359 subcortex-only solution did not produce stable results according to any criterion, whereas the 192 symmetrized subcortical features (Appendix 1—figure 13) and 12 network-level features (Appendix 1—figure 14) solutions captured statistically modest effects relative to the 180 symmetrized cortical features (Figure 5B–D). Therefore, we characterized the 180-parcel solution further. We examined two CCA solutions using these 180 neural features in relation to symptom scores: (i) all 36 item-level symptom measure scores from the PANSS/BACS (‘180 vs. 36 CCA’, Figure 5C,F); and (ii) five PC symptom scores (‘180 vs. 5 CCA’, Figure 5D,G, see Materials and methods). Only 2 out of 36 modes for the 180 vs. 36 CCA solution survived permutation testing and false discovery rate (FDR) correction (Figure 5C). In contrast, all 5 modes of the 180 vs. 5 CCA survived (Figure 5D). Critically, we found that no single CCA mode computed on item-level symptom measures captured more variance than the CCA modes derived from PC symptom scores, suggesting that the PCA-derived dimensions capture more neurally relevant variation than any one single clinical item (Figure 5E–G). Additional CCA details are presented in Appendix 1 - Note 3 and Appendix 1—figure 15.

We highlight here an example canonical variate, CV3, across both symptom and neural effects, with data for all CVs shown in Appendix 1—figure 16. CV3 scores across diagnostic groups normalized to controls are shown Figure 5H in addition to how CV3 loads onto each PC (Figure 5I). The negative loadings on to PCs 1, 4, and 5 and the high positive loadings on to PC3 in Figure 5I indicate that CV3 captures some shared variation across symptom PCs. This can also be visualized by computing how CV3 projects onto the original 36 symptom measures (Figure 5J). Finally, the neural loadings for CV3 are shown in Figure 5K.

Lastly, we tested if the 180 vs. 5 CCA solution is stable and reproducible, as done with PC-to-GBC univariate results. The CCA solution was robust when tested with k-fold and leave-site-out cross-validation (Appendix 1—figure 17) likely because these methods use CCA loadings derived from the full sample. However, the CCA loadings did not replicate in non-overlapping split-half samples (Figure 5L, see Appendix 1 - Note 4). Moreover, a leave-one-subject-out cross-validation revealed that removing a single subject from the sample affected the CCA solution such that it did not generalize to the left-out subject (Figure 5M). This is in contrast to the PCA-to-GBC univariate mapping, which was substantially more reproducible for all attempted cross-validations relative to the CCA approach. This is likely because substantially more power is needed to resolve a stable multivariate neurobehavioral effect with this many features (Dinga et al., 2019). Indeed, a multivariate power analysis using 180 neural features and five symptom features and assuming a true canonical correlation of $r=0.3$ suggests that a minimal sample size of $N=8,145$ is needed to sufficiently detect the effect (Helmer et al., 2020), which is an order of magnitude greater than the available sample size (Appendix 1—figure 18). Therefore, we leverage the univariate symptom-neural result for subsequent subject-specific model optimization and comparisons to molecular neuroimaging maps.

Most studies look for differences between clinical and control groups, but to our knowledge no study has tested whether both PSD and healthy controls actually share a major portion of neural variance that may be present across all people. If the bulk of the neural variance is similar across both PSD and CON groups then including this clinically irrelevant neural signal might obscure clinical-relevant neurobehavioral relationships. To test this, we examined the shared variance structure of the neural signal for all PSD patients (N = 436) and all controls (N = 202) independently by conducting a PCA on the GBC maps (see Materials and methods). Patients’ and controls’ neural signals were highly similar for each of the first three neural PCs (>30% of all neural variance in each group)(Appendix 1—figure 19A–J). These neural PCs may reflect a ‘core’ shared symptom-irrelevant neural variance that generalizes across all people. These data suggest that the bulk of neural variance in PSD patients may actually not be symptom-relevant, which highlights the importance of optimizing symptom-neural mapping. Under this assumption, removing the shared neural variance between PSD and CON should not drastically affect the reported symptom-neural univariate mapping solution, because this common variance does not map to clinical features. Thus, we conducted a PCA using the parcellated GBC data from all 436 PSD and 202 CON (see Materials and methods), which we will refer to as ‘GBC-PCA’ to avoid confusion with the symptom/behavioral PCA. This GBC-PCA resulted in 637 independent GBC-PCs. Since PCs are orthogonal to each other, we then partialled out the variance attributable to GBC-PC1 from the PSD data by reconstructing the PSD GBC matrix using only scores and coefficients from the remaining 636 GBC-PCs (${\widehat{GBC}}_{woPC1}$). We then reran the univariate regression as described in Figure 3, using the same five symptom PC scores across 436 PSD (Appendix 1—figure 20). Removing the first PC of shared neural variance (which accounted for about 15.8% of the total GBC variance across CON and PSD) from PSD data attenuated the statistics slightly (not unexpected as the variance was by definition reduced) but otherwise did not strongly affect the univariate mapping solution. We repeated the symptom-neural regression next with the first two GBC-PCs partialled out of the PSD data Appendix 1—figure 21, with the first three PCs parsed out Appendix 1—figure 22, and with the first four neural PCs parsed out Appendix 1—figure 23. The symptom-neural maps remain fairly robust, although the similarity with the original ${\beta }_{PC}GBC$ maps does drop as more common neural variance is parsed out.

Above we demonstrated that symptom PC scores can be reliably computed across sites and cross-validation approaches (Figure 2). Next, we show that leave-one-subject-out cross-validation yields reliable effects for the low-rank symptom PCA solution (Figure 6A). This stable single-subject PC score prediction provides the basis for testing if the derived neural maps can yield an individually reliable subset of features. To this end, we developed a framework for univariate neural feature selection (Appendix 1—figure 23) based on PC scores (e.g. PC3 score). Specifically, we computed a dot product GBC metric ($dpGBC$) that provides an index of similarity between an individual ${\displaystyle \mathrm{\Delta }GBC}$ topography relative to a ‘reference’ group-level ${\beta }_{PC}GBC$ map (see Materials and methods and Appendix 1—figure 23). Using this $dpGBC$ index via a feature selection step-down regression, we found a subset of parcels for which the symptom-neural statistical association was maximal (Figure 6A). For PC3, we found $P=39$ maximally predictive parcels out of the group neural ${\beta }_{PC}GBC$ map. Specifically, the relationship between PC3 symptom scores and $dpGBC$ values across subjects was maximal (Figure 6B, top panel, $r=0.36$) as was the relationship between predicted $dpGB{C}^{pred}$ vs. observed $dpGB{C}^{obs}$ (Figure 6B, bottom panel, $r=0.35$) (see Appendix 1—figure 25 for all PCs). Importantly, the ‘subset’ feature map (i.e. ${[{\beta }_{PC3}GB{C}^{obs}]}^{P=39}$, Figure 6C) exhibited improved statistical properties relative to the full map (i.e. ${[{\beta }_{PC3}GB{C}^{obs}]}^{P=718}$, Figure 6D). Furthermore, the relationship between observed vs. predicted subset feature maps (i.e. $r[dpGB{C}^{obs},dpGB{C}^{pred}]$) was highly consistent across DSM diagnoses (Figure 6E). Finally, a single patient is highlighted for whom the correlation between predicted and observed subset feature maps was high (i.e. ${\displaystyle r[\mathrm{\Delta }GB{C}^{pred},\mathrm{\Delta }GB{C}^{obs}]}$, Figure 6F), demonstrating that the dimensionality-reduced symptom scores can be used to quantitatively optimize symptom-neural map features for individuals.

Figure 6

Download asset Open asset

We showed a quantitative framework for computing a neurobehavioral model at the single-patient level. This brain-behavior space (BBS) model was optimized along a single dimensionality-reduced symptom PC axis. Next, we tested a hybrid patient selection strategy by first imposing a PC-based symptom threshold, followed by a target neural similarity threshold driven by the most highly predictive symptom-neural map features. This is described in Materials and methods, Appendix 1 - Note 5, and Appendix 1—figure 26–27. We found that for patients with a high (either positive or negative) PC symptom score the symptom-neural relationship was robust (Appendix 1—figure 26A-D), and these patients could be reliably selected for that particular dimension (Appendix 1—figure 26E,G). Conversely, patients with a low absolute PC score showed a weak relationship with symptom-relevant neural features. This is intuitive because these individuals do not vary along the selected PC symptom axis. We also tested this patient selection strategy on an independent cross-diagnostic ‘replication’ sample, which yielded consistent results (Appendix 1—figure 26F,H). Collectively, these results show that data-driven symptom scores can pinpoint individual patients for whom neural variation strongly maps onto a target neural reference map. These data also highlight that both symptom and neural information for an independent patient can be quantified in the reference ‘discovery’ BBS using their symptom data alone.

Next, we use the personalized BBS selection in a proof-of-concept framework for informing molecular mechanism of possible treatment response by relating subject-specific BBS to independently acquired pharmacological neuroimaging maps. Here, we examine two mechanisms implicated in PSD neuropathology via ketamine, a N-methyl-D-aspartate (NMDA) receptor antagonist (Krystal et al., 1994), and lysergic acid diethylamide (LSD), primarily a serotonin receptor agonist (Preller et al., 2018; González-Maeso et al., 2007; Egan et al., 1998). We first quantified individual subjects’ BBS ‘locations’ in the established reference neurobehavioral geometry. The radarplot in Figure 7A shows original symptoms whereas Figure 7B shows ${\displaystyle \mathrm{\Delta }GB{C}^{obs}}$ maps for two patients from the replication dataset (denoted here with $X_{PC3}$ and $Y_{PC3}$, see Appendix 1—figure 11 for other example patients). Both of these patients exceeded the neural and behavioral BBS selection indices for PC3 (defined independently in the ‘discovery’ dataset, Appendix 1—figure 26C). Furthermore, both patients exhibited neurobehavioral variation in line with their expected locations in the BBS geometry. Specifically, patient $X_{PC3}$ from the replication dataset scored highly negatively on the PC3 axis defined in the ‘discovery’ PSD sample (Figure 7A). In contrast, patient $Y_{PC3}$ scored positively on the PC3 axis. Importantly, the correlation between the ${\displaystyle \mathrm{\Delta }GB{C}^{obs}}$ map for each patient and the group-reference ${\beta }_{PC3}GBC$ was directionally consistent with their symptom PC score (Figure 7B–C).

Figure 7

Download asset Open asset

We then tested if the single-subject BBS selection could be quantified with respect to a neural map reflecting glutamate receptor manipulation, a hypothesized mechanism underlying PSD symptoms (Moghaddam and Javitt, 2012). Specifically, we used an independently collected ketamine infusion dataset, collected in healthy adult volunteers during resting-state fMRI (Anticevic et al., 2012a). As with the clinical data, here we computed a ${\displaystyle \mathrm{\Delta }GBC}$ map reflecting the effect of ketamine on GBC relative to placebo (Materials and methods). The maximally predictive PC3 parcels exhibited high spatial similarity with the ketamine map (ρ = 0.76, see Materials and methods), indicating that the ${\displaystyle \mathrm{\Delta }GBC}$ pattern induced by ketamine tracks with the GBC pattern reflecting PC3 symptom variation.

Critically, because $X_{PC3}$ is negatively loaded on the PC3 symptom axis, an NMDA receptor antagonist like ketamine may modulate symptom-relevant circuits in a way that reduces similarity with the PC3 map. This may in turn have an impact on the PC3-like symptoms. Consistent with this hypothesis, $X_{PC3}$ expresses predominantly depressive symptoms (Figure 7A), and ketamine has been shown to act as an anti-depressant (Berman et al., 2000). This approach can be applied for patients that load along another axis, such as PC5. Figure 7D–E shows the symptom and neural data for two patients whom met thresholds for PC5 selection (Appendix 1—figure 26C). Notably, the selected PC5 map is anti-correlated with a ${\displaystyle \mathrm{\Delta }GBC}$ map reflecting LSD vs. placebo effects (Preller et al., 2018) (ρ = −0.44, Figure 7F). Hence areas modulated by LSD may map onto behavioral variation along PC5. Consequently, serotonergic modulation may be effective for treating $Q_{PC5}$ and $Z_{PC5}$, via an antagonist or an agonist respectively. These differential similarities between pharmacological response maps and BBS maps (Appendix 1—figure 29) can be refined for quantitative patient segmentation.

To further inform molecular mechanism for the derived BBS results, we compared results with neural gene expression profiles derived from the Allen Human Brain Atlas (AHBA) (Hawrylycz et al., 2012; Burt et al., 2018; Figure 8A, Materials and methods). Specifically, we tested if BBS cortical topographies, which reflect stable symptom-neural mapping along PSD, covary with the expression of genes implicated in PSD neuropathology. We focus here on serotonin receptor subunits (HTR1E, HTR2C, HTR2A), GABA receptor subunits (GABRA1, GABRA5), and the interneuron markers somatostatin (SST) and parvalbumin (PVALB). Serotonin agonists such as LSD have been shown to induce PSD-like symptoms in healthy adults (Preller et al., 2018) and the serotonin antagonism of ‘second-generation’ antipsychotics are thought to contribute to their efficacy in targeting broad PSD symptoms (Geyer and Vollenweider, 2008; Meltzer et al., 2012; Abi-Dargham et al., 1997). Abnormalities in GABAergic interneurons, which provide inhibitory control in neural circuits, may contribute to cognitive deficits in PSD (Benes and Berretta, 2001; Inan et al., 2013; Dienel and Lewis, 2019) and additionally lead to downstream excitatory dysfunction that underlies other PSD symptoms (Lisman et al., 2008; Grace, 2016). In particular, a loss of prefrontal parvalbumin-expression fast-spiking interneurons has been implicated in PSD (Enwright Iii et al., 2018; Lodge et al., 2009; Beasley and Reynolds, 1997; Lewis et al., 2012). Figure 8B shows the distribution of correlations between the PC3 map and the cortical expression patterns of 20,200 available AHBA genes (results for other PCs are shown in Appendix 1—figure 30). Seven genes of interest are highlighted, along with their cortical expression topographies and their similarity with the PC3 BBS map (Figure 8C–E). This BBS-to-gene mapping can potentially reveal novel therapeutic molecular targets for neurobehavioral variation. For example, the HTR1E gene, which encodes the serotonin 5-HT_1E receptor, is highly correlated with the PC3 BBS map. This could drive further development of novel selective ligands for this receptor, which are not currently available (Kitson, 2007).

Figure 8

Download asset Open asset

Psychosis spectrum is associated with notable clinical heterogeneity such deficits in cognition as well as altered beliefs (i.e. delusions), perception (i.e. hallucinations), and affect (i.e. negative symptoms) (Lefort-Besnard et al., 2018). This heterogeneity is captured by clinical instruments that quantify PSD symptoms across dozens of specific questions and ratings. This yields a high-dimensional symptom space that is intractable for reliable mapping of neural circuits (Helmer et al., 2020). Here we show that a low-rank solution captures principal axes of PSD symptom variation, a finding in line with prior work in SZP (van der Gaag et al., 2006a; van der Gaag et al., 2006b; Lindenmayer et al., 1994; Emsley et al., 2003; White et al., 1997; Dollfus et al., 1996; Blanchard and Cohen, 2006; Chen et al., 2020; Lefort-Besnard et al., 2018).

These results highlights two key observations: (i) Existing symptom reduction studies (even those in SZP specifically) have not evaluated solutions that include cognitive impairment – a hallmark deficit across the psychosis spectrum (Barch et al., 2013). Here we show that cognitive performance captures a notable portion of the symptom variance independent of other axes. We observed that cognitive variation captured 10% of PSD sample variance even after accounting for ‘Global Functioning’ psychopathology. (ii) No study has quantified cognitive deficit variation along with core psychosis symptoms via dimensionality reduction across multiple PSD diagnoses. While existing studies have evaluated stability of data-reduced solutions within a single DSM category (van der Gaag et al., 2006a; Lindenmayer et al., 1995; Lefort-Besnard et al., 2018), current results show that dimensionality-reduced PSD symptom solutions can be reproducibly obtained across DSM diagnoses.

Importantly, on each data-reduced symptom axis, some subset of PSD patients received a score near zero. This does not imply that these patients were unimpaired; rather, the symptom configurations for these patients were orthogonal to variation along this specific axis. While in the current paper we demonstrate individual selection based on symptom scores along a single axis, it may also be possible to use a combination of several symptom dimensions to further inform individual treatment selection. For example, for patients with above-threshold symptom scores along several dimensions, a combination of these PCs may be more predictive of their neural features and in turn be a better indication of treatment selection for these individuals. Alternatively, the symptom dimensions can be used as exclusionary criteria for selectively picking out individuals along only one particular axis (e.g. filtering out patients above-threshold on PCs 1, 3, 4, and five to selectively study patients with only severe PC2 cognitive deficits). Additionally, the ICA solution revealed oblique axes which are linear combinations of the PCs, for example IC2 which is oblique to PCs 1, 2, and 3 (Appendix 1—figure 3). However, these ICs do not result in as robust or unique a neural circuit mapping as do the PCs, as shown in Appendix 1—figure 12. The observation that PSD are associated with multiple complex symptom dimensions highlights an important intuition that may extend to other mental health spectra. Additionally, the PSD symptom axes reported here are neither definitive nor exhaustive. In fact, close to 50% of all clinical variance was not captured by the symptom PCA – an observation often overlooked in symptom data-reduction studies, which focus on attaining ‘predictive accuracy’. Such studies rarely consider how much variance remains unexplained in the final data-reduced model and, relatedly, if the proportion of explained variance is reproducible across samples. This is a key property for reliable symptom-to-neural mapping. Thus, we tested if this reproducible low-dimensional PSD symptom space robustly mapped onto neural circuit patterns.

We show that the dimensionality-reduced symptom space improved the mapping onto neural circuit features (i.e. GBC), as compared to a priori item-level clinical symptom measures (Figure 3). This symptom-neural mapping was highly reproducible across various cross-validation procedures and split-half replication (Figure 4). The observed statistical symptom-neural improvement after dimensionality reduction suggests that data-driven clinical variation more robustly covaried with neural features. As noted, the low-rank symptom axes generalized across DSM diagnoses. Consequently, the mapping onto neural features (i.e. GBC) may have been restricted if only a single DSM category or clinical item was used. Importantly, as noted, traditional clinical scales are uni-directional (i.e. zero is asymptomatic, hence there is an explicit floor). Here, we show that data-driven symptom axes (e.g. PC3) were associated with bi-directional variation (i.e. no explicit floor effect). Put differently, patients who score highly on either end of these data-driven axes are severely symptomatic but in very different ways. If these axes reflect clinically meaningful phenomena at both tails then they should more robustly map to neural feature variation, which is in line with reported effects. Therefore, by definition, the derived map for each of the PCs will reflect the neural circuitry that may be modulated by the behaviors that vary along that PC (but not others). For example, we named the PC3 axis ‘Psychosis Configuration’ because of its strong loadings onto conventional ‘positive’ and ‘negative’ PSD symptoms. This PC3 ‘Psychosis Configuration’ showed strong positive variation along neural regions that map onto the canonical default mode network (DMN), which has frequently been implicated in PSD (Fryer et al., 2013; Anticevic et al., 2012b; Ongür et al., 2010; Woodward et al., 2011; Baker et al., 2014; Meda et al., 2014) suggesting that individuals with severe ‘positive’ symptoms exhibit broadly elevated connectivity with regions of the DMN. On the contrary, this bi-directional ‘Psychosis Configuration’ axis also showed strong negative variation along neural regions that map onto the sensory-motor and associative control regions, also strongly implicated in PSD (Ji et al., 2019a; Anticevic et al., 2014). The ‘bi-directionality’ property of the PC symptom-neural maps may thus be desirable for identifying neural features that support individual patient selection. For instance, it may be possible that PC3 reflects residual untreated psychosis symptoms in this chronic PSD sample, which may reveal key treatment neural targets. In support of this circuit being symptom-relevant, it is notable that we observed a mild association between GBC and PC scores in the CON sample (Appendix 1—figure 10). The symptom-neural mapping of dimensionality-reduced symptom scores also produced notable neural maps along other PC axes (Appendix 1—figure 9). Higher PC1 scores, indicating higher general functioning, may be associated with lower global connectivity in sensory/cingulate cortices and thalamus, but higher global connectivity with the temporal lobe. Higher PC2 scores, indicating better cognitive functioning were positively associated with higher medial prefrontal and cerebellar GBC, and negatively with visual cortices and striatum. Importantly, the unique neural circuit variance shown in the ${\beta }_{PC2}GBC$ map suggest that there are PSD patients with more (or less) severe cognitive deficits independent of any other symptom axis, which would be in line with the observation that these symptoms are not treatable with antipsychotic medication (and therefore should not correlate with symptoms that are treatable by such medications; i.e. PC3). Of note, the statistics in the ${\beta }_{PC2}GBC$ map were relatively moderate, opening up the possibility that a more targeted measure of neural connectivity (e.g. seed FC, see below) would reveal a robust circuit map for this particular symptom dimension. Given the key role of cognitive deficits in PSD (Barch et al., 2013), this will be a particularly important direction to pursue, especially in prodromal or early-stage individuals where cognitive dysfunction has been shown to predict illness trajectories (Fusar-Poli et al., 2012; Bora et al., 2014; Bora, 2015; Antshel et al., 2017). The PC4 – Affective Valence axis was positively associated with higher GBC in cingulate and sensori-motor cortices and subcortically with the thalamus, basal ganglia, and anterior cerebellum, which may be implicated in the deficits in social functioning and affective aspects of this symptom dimension. Lastly, positive scores on PC5, which we named ‘Agitation/Excitement’, were associated with broad elevated global connectivity in the frontal associative and sensori-motor cortices, thalamus, and basal ganglia. Interestingly, however, the nucleus accumbens (as well as hippocampus, amygdala, cerebellum, and temporal lobe) appear to be negatively associated with PC5 score (i.e. more severe grandiosity, unusual thought, delusional and disorganized thinking), illustrating another example of a bi-directional axis.

Deriving a neurobehavioral mapping that is resolvable and stable at the individual patient level is a necessary benchmark for deploying symptom-neural ‘biomarkers’ in a clinically useful way. Therefore, there is increasing attention placed on the importance of achieving reproducibility in the psychiatric neuroimaging literature (Noble et al., 2019; Balsters et al., 2016; Cao et al., 2019; Woo et al., 2017), which becomes especially important for individualized symptom-neural mapping. Recently, several efforts have deployed multivariate methods to quantify symptom-neural relationships (Drysdale et al., 2017; Xia et al., 2018; Smith et al., 2015; Moser et al., 2018; Rodrigue et al., 2018; Yu et al., 2019), highlighting how multivariate techniques may perhaps provide clinically innovative insights. However, such methods face the risk of overfitting for high-dimensional but underpowered datasets (Dinga et al., 2019), as recently shown via formal generative modeling (Helmer et al., 2020).

Here, we attempted to use multivariate solutions (i.e. CCA) to quantify symptom and neural feature co-variation. In principle, CCA is well-suited to address the brain-behavioral mapping problem. However, symptom-neural mapping using CCA in our sample was not reproducible even when using a low-dimensional symptom solution and parcellated neural data as a starting point. Therefore, while CCA (and related multivariate methods such as partial least squares) are theoretically appropriate and may be helped by regularization methods such as sparse CCA, in practice many available psychiatric neuroimaging datasets may not provide sufficient power to resolve stable multivariate symptom-neural solutions (Helmer et al., 2020). A key pressing need for forthcoming studies will be to use multivariate power calculators to inform sample sizes needed for resolving stable symptom-neural geometries at the single subject level. Of note, athough we were unable to derive a stable CCA in the present sample, this does not imply that the multivariate neurobehavioral effect may not be reproducible with larger effect sizes and/or sample sizes. Critically, this does highlight the importance of power calculations prior to computing multivariate brain-behavioral solutions (Helmer et al., 2020).

Consequently, we tested if a low-dimensional symptom solution can be used in a univariate symptom-neural model to optimize individually predictive features. Indeed, we found that a univariate brain-behavioral space (BBS) relationship can result in neural features that are stable for individualized prediction. Critically, we found that if a patient exhibited a high PC symptom score, they were more likely to exhibit a topography of neural ΔGBC that was topographically similar to the ${\beta }_{PC}GBC$ map. This suggests that optimizing such symptom-neural mapping solutions (and ultimately extending them to multivariate frameworks) can inform cross-diagnostic patient segmentation with respect to symptom-relevant neural features. Importantly, this could directly inform patient identification based on neural targets that are of direct symptom relevance for clinical trial design.

Selecting single patients via stable symptom-neural mapping of BBS solutions is necessary for individual patient segmentation, which may ultimately inform treatment indication. However, it is also critical to be able to relate the derived symptom-neural maps to a given mechanism. Here we highlight two ways to ‘benchmark’ the derived symptom-neural maps by calculating their similarity against independent pharmacological neuroimaging and gene expression maps. We show a proof-of-principle framework for quantifying derived symptom-neural reference maps with two PSD-relevant neuropharmacological manipulations derived in healthy adults via LSD and ketamine. These analyses revealed that selecting single patients, via the derived symptom-neural mapping solution, can yield above-chance quantitative correspondence to a given molecular target map. These data highlight an important effect: it is possible to construct a ‘strong inference’ (Platt, 1964) evaluation of single patients’ differential similarity to one molecular target map versus another. For instance, this approach could be applied to maps associated with already approved PSD treatments (such as clozapine, olanzapine, or chlorpromazine [Lally and MacCabe, 2015; Miyamoto et al., 2005]) to identify patients with symptom-neural configurations that best capture available treatment-covarying neural targets.

Relatedly, AHBA gene expression maps (Burt et al., 2018) may provide an a priori benchmark for treatment targets that may be associated with a given receptor profile. Here we show that identified BBS maps exhibit spatial correspondence with neural gene expression maps implicated in PSD – namely serotonin, GABA and interneuron gene expression. This gene-to-BBS mapping could be then used to select those patients that exhibit high correspondence to a given gene expression target.

Collectively, this framework could inform empirically testable treatment selection methods (e.g. a patient may benefit from ketamine, but not serotonin agonists such as LSD/psilocybin). In turn, this independent molecular benchmarking framework could be extended to other approaches (e.g. positron emission tomography (PET) maps reflecting specific neural receptor density patterns [Farde et al., 1986; Arakawa et al., 2008]) and iteratively optimized for quantitative patient-specific selection against actionable molecular targets.

There are several constraints of the current result that require future optimization – namely the ability to generalize across time (early course vs. chronic patients), across a range of symptom severity (e.g. severe psychotic episode or persistent low-severity psychosis) and across distinct symptom spectra (e.g. mood). This applies to both the low-rank symptom solution and the resulting symptom-neural mapping. It is possible that the derived lower-dimensional symptom solution, and consequently the symptom-neural mapping solution, exhibits either time-dependent (i.e. state) or severity-dependent (i.e. trait) re-configuration. Relatedly, medication dose, type, and timing may also impact the solution. Another important aspect that will require further characterization is the possibility of oblique axes in the symptom-neural geometry. While orthogonal axes derived via PCA were appropriate here and similar to the ICA-derived axes in this solution, it is possible that oblique dimensions more clearly reflect the geometry of other psychiatric spectra and/or other stages in disease progression. For example, oblique components may better capture dimensions of neurobehavioral variation in a sample of prodromal individuals, as these patients are exhibiting early-stage psychosis-like symptoms and may show signs of diverging along different trajectories.

Critically, these factors should constitute key extensions of an iteratively more robust model for individualized symptom-neural mapping across the PSD and other psychiatric spectra. Relatedly, it will be important to identify the ‘limits’ of a given BBS solution – namely a PSD-derived effect may not generalize into the mood spectrum (i.e. both the symptom space and the resulting symptom-neural mapping is orthogonal). It will be important to evaluate if this framework can be used to initialize symptom-neural mapping across other mental health symptom spectra, such as mood/anxiety disorders.

These types of questions will require longitudinal and clinically diverse studies that start prior to the onset of full-blown symptoms (e.g. the North American Prodrome Longitudinal Study (NAPLS) [Addington et al., 2012; Seidman et al., 2010]). A corollary of this point is that ∼50% of unexplained symptom variance in the current PCA solution necessitates larger samples with adequate power to map this subtle, but perhaps clinically essential, PSD variation.

Notably, the current cohort was adequately powered for symptom data reduction that drove univariate neural mapping. However, this sample was insufficiently powered for resolving stable multivariate symptom-neural relationships even with dimensionality-reduced symptom features. Consequently, the limited sample size necessitated choices for dimensionality reduction of the neural feature space in this study even for univariate analyses. While both parcellation and GBC constitute principled choices, symptom-relevant neural information may have been lost (which may be embedded in a higher dimensional space). One obvious solution is to increase sample sizes (e.g. via datasets such as the UK Biobank [Sudlow et al., 2015]). However, in parallel, it will be critical to develop neurobiologically informed feature space reduction and/or to optimize the stability of multivariate solutions via regularization methods. Another improvement would be to optimize neural data reduction sensitivity for specific symptom variation (Shehzad et al., 2014). We chose to use GBC for our initial geometry characterizations as it is a principled and assumption-free data-reduction metric that captures (dys)connectivity across the whole brain and generates neural maps (where each region can be represented by a value, in contrast to full functional connectivity matrices) that are necessary for benchmarking against molecular imaging maps. However, GBC is a summary measure that by definition does not provide information regarding connectivity between specific pairs of neural regions, which may prove to be highly symptom-relevant and informative. Thus symptom-neural relationships should be further explored with higher-resolution metrics, such as restricted GBC (Anticevic et al., 2013), which can summarize connectivity information for a specific network or subset of neural regions, or seed-based FC using regions implicated in PSD (e.g. thalamus [Anticevic et al., 2014; Woodward et al., 2011]).

Here, we focused on the neural blood oxygen level dependent (BOLD) signal from fMRI. However, other modalities such as diffusion-weighted imaging (DWI), PET imaging, or electroencephalography (EEG) could be leveraged. Additional clinically-relevant information could be derived from genetics (such as polygenic risk scores [Hyman, 2011; O'Donovan et al., 2009; Cao et al., 2020]) or ecological momentary assessment (EMA) (Corcoran et al., 2018; Niendam et al., 2018), especially to parse state vs. trait biomarker variation.

Lastly, building on the proof-of-concept molecular neuroimaging comparisons, it will be imperative to eventually test such predictions in formal clinical trials. An actionable next step would be to optimize patient selection against existing treatments, which could result in higher success rates for drug development trials and potentially have massive impact for developing new interventions. Critically, the opportunity to develop, validate, and refine an individual-level quantitative framework could deliver a more rapid and cost-effective way of pairing patients with the right treatments.

We show that complex and highly heterogeneous PSD symptom variation can be robustly reduced into a low-rank symptom solution that is cross-diagnostic, individually predictive, generalizable and incorporates cognitive deficits. In turn, the derived PSD symptom axes robustly mapped onto distinct yet reproducible neural patterns, which were predictive at the single-patient level. Leveraging these stable results, we show a proof-of-principle framework for relating the derived symptom-relevant neural maps, at the individual patient level, with molecular targets implicated in PSD via LSD and ketamine neuro-pharmacological manipulations. Lastly, we used AHBA gene expression maps to show that identified PSD symptom-relevant neural maps co-vary with serotonin, GABA and interneuron neural gene expression patterns. This set of symptom-neural mapping results can be iteratively and quantitatively optimized for personalized treatment segmentation endpoints.

We used publicly available behavioral and neural data from the Bipolar-Schizophrenia Network on Intermediate Phenotypes (BSNIP) consortium (Tamminga et al., 2014). All data were obtained from the National Data Archive (NDA) repository (https://nda.nih.gov/edit_collection.html?id=2274). Participants were collected at six sites across the United States. Full recruitment details are provided in prior publications (Tamminga et al., 2014; Meda et al., 2015; Sheffield et al., 2017). In brief, participants were recruited through advertising in Baltimore MD, Boston MA, Chicago IL, Dallas TX, and Hartford CT. All assessments were standardized across sites. Participants were excluded if they had (i) a history of seizures or head injury resulting in >10 min loss of consciousness, (ii) positive drug screen on the day of testing, (iii) a diagnosis of substance abuse in the past 30 days or substance dependence in the past 6 months, (iv) history of serious medical or neurological disorder that would likely affect cognitive functioning, (v), history of serious medical or neurological disorder that would likely affect cognitive functioning, (vi) insufficient English proficiency, or (vii) an age-corrected Wide-Range Achievement Test (4th edition) reading test standard score <65. Additionally, participants were required to have had no change in medication and been clinically stable over the past month. Participants completed a SCID interview and were given a diagnosis via consensus from study clinicians; participants with an Axis one clinical psychosis diagnosis were additionally given assessments including the Positive and Negative Symptom Scale (PANSS)(Kay et al., 1987). Note that apart from the measures used in this paper, the full BSNIP dataset includes a rich characterization of measures from multiple modalities, including electrophysiological, eye tracking, structural and diffusion neuroimaging, as well as a number of clinical batteries, which are not quantified in this study. We used data from a total of 638 participants with complete behavioral and neural data after preprocessing and quality control, including 202 healthy controls (CON), 167 patients diagnosed with schizophrenia (SZP), 119 patients with schizoaffective disorder (SADP), and 150 patients with bipolar disorder with psychosis (BPP). For full demographic and clinical characteristics of the sample see Appendix 1—table 2.

Participants completed a neural magnetic resonance imaging (MRI) scan at 3T, including resting-state functional blood-oxygen-level-dependent imaging (BOLD) and a magnetization-prepared rapid gradient-echo (MP-RAGE) sequence for T1 weighted data. Full details on scanners and acquisition parameters used at each of the sites have previously been described (Meda et al., 2012) and are summarized here in Appendix 1—table 3. Neuroimaging data were preprocessed using the Human Connectome Project (HCP) minimal preprocessing pipeline (Glasser et al., 2013), adapted for compatibility with ‘legacy’ data, which are now featured as a standard option in the HCP pipelines provided by our team (https://github.com/Washington-University/HCPpipelines/pull/156, copy archived at swh:1:rev:1334b35ab863540044333bbdec70a68fb19ab611; Repovš et al., 2021). These modifications to the HCP pipelines were necessary as the BSNIP data did not include a standard field map and did not incorporate a T2w high-resolution image. The adaptations for single-band BOLD acquisition have previously been described in detail (Ji et al., 2019a).

A summary of the HCP pipelines is as follows: the T1-weighted structural images were first aligned by warping them to the standard Montreal Neurological Institute-152 (MNI-152) brain template in a single step, through a combination of linear and non-linear transformations via the FMRIB Software Library (FSL) linear image registration tool (FLIRT) and non-linear image registration tool (FNIRT) (Jenkinson et al., 2002). Next, FreeSurfer's recon-all pipeline was used to segment brain-wide gray and white matter to produce individual cortical and subcortical anatomical segmentations (Reuter et al., 2012). Cortical surface models were generated for pial and white matter boundaries as well as segmentation masks for each subcortical gray matter voxel. Using the pial and white matter surface boundaries, a ‘cortical ribbon’ was defined along with corresponding subcortical voxels, which were combined to generate the neural file in the Connectivity Informatics Technology Initiative (CIFTI) volume/surface ‘grayordinate’ space for each individual subject (Glasser et al., 2013). BOLD data were motion-corrected by aligning to the middle frame of every run via FLIRT in the initial NIFTI volume space. In turn, a brain-mask was applied to exclude signal from non-brain tissue. Next, cortical BOLD data were converted to the CIFTI gray matter matrix by sampling from the anatomically defined gray matter cortical ribbon and subsequently aligned to the HCP atlas using surface-based nonlinear deformation (Glasser et al., 2013). Subcortical voxels were aligned to the MNI-152 atlas using whole-brain non-linear registration and then the Freesurfer-defined subcortical segmentation applied to isolate the subcortical grayordinate portion of the CIFTI space.

After the HCP minimal preprocessing pipelines, movement scrubbing was performed (Power et al., 2013). As done in our prior work (Anticevic et al., 2012a), all BOLD image frames with possible movement-induced artifactual fluctuations in intensity were flagged using two criteria: frame displacement (the sum of the displacement across all six rigid body movement correction parameters) exceeding 0.5 mm (assuming 50 mm cortical sphere radius) and/or the normalized root mean square (RMS) (calculated as the RMS of differences in intensity between the current and preceding frame, computed across all voxels and divided by the mean intensity) exceeding 1.6 times the median across scans. Any frame that met one or both of these criteria, as well as the frame immediately preceding and immediately following, were discarded from further preprocessing and analyses. Subjects with more than 50% frames flagged using these criteria were excluded. Next, a high-pass filter (threshold 0.008 Hz) was applied to the BOLD data to remove low-frequency signals due to scanner drift. In-house Matlab code was used to calculate the average variation in BOLD signal in the ventricles, deep white matter, and across the whole gray matter, as well as movement parameters. These signals, as well as their first derivatives to account for delayed effects, were then regressed out of the gray matter BOLD time series as nuisance variables (as any change in the BOLD signal due to these variables are persistent and likely not reflecting neural activity) (Power et al., 2018). It should be noted that using global signal regression to remove spatially persistent artifact is highly controversial in neuroimaging (Power et al., 2017; Yang et al., 2017), but it remains the field-wide gold standard (though see other recent and emerging approaches at [Glasser et al., 2018; Aquino et al., 2019]).

The behavioral measures analyzed included the Positive and Negative Syndrome Scale (PANSS), an assessment of psychosis symptom severity (Kay et al., 1987), and the Brief Assessment of Cognition in Schizophrenia (BACS) battery, which provided an assessment of cognitive functioning (Keefe et al., 2004). Only control subjects with complete BACS measures were used for analyses; PANSS symptom scores were imputed for control subjects for whom the PANSS had not been administered under the assumption that these subjects were asymptomatic. Only patient subjects with complete PANSS and BACS measures were used in analyses. The BACS scores used here are presented as standardized Z-scores normalized to the mean and standard deviation of the control group for each BSNIP site, as done in prior work (Tamminga et al., 2014). The full PANSS battery is conventionally divided into three sub-scales: Positive symptom scale (seven items), Negative symptom scale (seven items), and General Psychopathology symptom scale (16 items). The BACS battery consists of 6 individual sub-scores (Keefe et al., 2004). In total, this yielded 36 symptom variables per participant. Effects of symptom variation across assessment sites have been rigorously characterized in prior publications (Tamminga et al., 2013). Nevertheless, we explicitly tested for site effects in our analyses, described in detail below.

The principal component analysis (PCA) of behavioral data was computed by including all 36 symptom variables across all N = 436 patients. Variables were first scaled to have unit variance across patients before running the PCA. Significance of the derived principal components (PCs) was computed via permutation testing. For each permutation, patient order was randomly shuffled for each symptom variable before re-computing PCA. This permutation was repeated 5000 times to establish the null model. PCs which accounted for a proportion of variance that exceeded chance (p<0.05 across all 5000 permutations) were retained for further analysis.

Stability of the symptom PCA solution was assessed with within-sample cross-validation. To evaluate if there were site differences that uniquely drove the PCA solution, we performed a leave-site-out cross-validation analysis. Specifically, we re-ran the PCA using all patients except those from one site, which was held out. This process was repeated for each of the six sites. To further evaluate the stability of the derived PCA solutions, we performed 1000 runs of k-fold cross-validation for values of k between 2 and 10. For each k-fold run, the full sample of patients was randomly divided into equally-sized k sets and a PCA was re-computed using subjects in k-1 sets (as the left-out set was used as the held-out sample).

For each run of leave-site-out and k-fold cross-validations significance was assessed via permutation testing as described above. The number of significant PCs and the total proportion of variance explained by these significant PCs remained highly stable across all runs (see Figure 2A and Appendix 1—figure 4 and Appendix 1—figure 5). Additionally, we compared observed and predicted PCA scores for the held-out sample of patients. Predicted scores for the held-out sample of patients were computed by weighing their raw symptom scores with the loadings from the PCA computed in all other patients. Observed scores for held-out patients were obtained from the original PCA computed on the full sample presented in the main text. The similarity between predicted and observed scores was high for all five significant PCs across all runs of leave-site-out and k-fold cross-validation, exceeding r = 0.9 in most analyses (see Figure 2B–C and Appendix 1—figure 4 and 5). Notably, the PCA solution did not show medication status or dosage effects (Appendix 1—figure 6). We also examined effects of age, sex, and SES in Appendix 1—figure 7. Briefly, we observed a significant positive relationship between age and PC3 scores, which may be because older patients (whom presumably have been ill for a longer time) exhibit more severe symptoms along the positive PC3 – Psychosis Configuration dimension. We also observed a significant negative relationship between Hollingshead index of SES and PC1 and PC2 scores. Lower PC1 and PC2 scores indicate poorer general functioning and cognitive performance respectively, which is consistent with higher Hollingshead indices (i.e. lower-skilled jobs or unemployment and fewer years of education). We also found significant sex differences in PC2 – Cognitive Functioning, PC4 – Affective Valence, and PC5 – Agitation/Excitement scores.

We also assessed the similarity of the PCA loadings using leave-site-out and 1000 runs of fivefold cross-validation frameworks (Figure 2D). Importantly, this cross-validation was designed to test if the observed loadings remained stable (as opposed to predicted patient-level scores). The loadings for significant PCs from each leave-site-out PCA solution as well as each run of the fivefold cross-validation were highly correlated (Figure 2D).

We also assessed the reproducibility of the PCA solution using independent split-half samples. The leave-site-out and k-fold cross-validation PCA analyses by definition use overlapping samples of patients for each iteration. Therefore, we additionally conducted a full split-half replication using entirely non-overlapping sets of patients in each iteration. For each split-half iteration, the full patient sample was randomly divided into two sets with equal proportions of each of the three diagnostic groups (BPP, SADP, SZP). Then, a PCA was computed using each of the split-half patient samples. The loadings from the two PCA solutions were then evaluated for reproducibility. This process was repeated 1000 times. The loadings for significant PCs were again highly similar even when comparing PCA solutions derived from completely non-overlapping patient samples (Figure 2D).

To predict individual patient PC scores for the leave-one-out analysis (Appendix 1—figure 26A), a PCA was computed using all patients except one held-out patient (N = 435). In turn, the derived loadings were then used to compute the predicted PC scores for the left-out patient. This process was repeated until predicted PC scores were calculated for every patient. Finally, the predicted score for each patient was evaluated for reproducibility relative to the observed score obtained from the PCA solution computed using the full N = 436 sample of patients. In addition, we computed an independent component analysis (ICA) to evaluate the consistency of the behavioral data reduction solution across methods (see Materials and methods).

Following preprocessing, the functional connectivity (FC) matrix was calculated for each participant by computing the Pearson’s correlation between every grayordinate in the brain with all other grayordinates. A Fisher’s r-to-Z transform was then applied. Global brain connectivity (GBC) was calculate by computing every grayordinate’s mean FC strength with all other grayordinates (i.e. the mean, per row, across all columns of the FC matrix). GBC is a data-driven summary measure of connectedness that is unbiased with regards to the location of a possible alteration in connectivity (Cole et al., 2016) and is therefore a principled way for reducing the number of neural features while assessing neural variation across the entire brain.

• where $GBC(x)$ denotes the GBC value at grayordinate $x$;

• where $N$ denotes the total number of grayordinates;

• where ${\displaystyle \sum _{y=1}^N}$ denotes the sum from $y=1$ to $y=N$;

• where $r_{xy}$ denotes the correlation between the time-series of grayordinates $x$ and $y$;

For parcel-wise GBC maps (described below), we first computed the mean BOLD signal within each parcel (see section below for parcellation details) for each participant and then computed the pairwise FC between all parcels. Finally, to obtain the parcellated GBC metric we computed the mean FC for each parcel and all other parcels. This order of operations (first parcellating the dense data series and then computing GBC) was chosen because it resulted in stronger statistical values due to increased within-parcel signal-to-noise of the BOLD data (Figure 4A).

• where $GBC(p)$ denotes the GBC value at parcel $p$;

• where $N$ denotes the total number of parcels;

• where ${\displaystyle \sum _{y=1}^N}$ denotes the sum from $q=1$ to $q=N$;

• where $r_{pq}$ denotes the correlation between the time-series of parcels $p$ and $q$;

Here, we applied a recently developed Cole-Anticevic Brain Network Parcellation (CAB-NP) parcellation (Ji et al., 2019b), which defines functional networks and regions across cortex and subcortex that leveraged the Human Connectome Project's Multi-Modal Parcellation (MMP1.0) (Glasser et al., 2016; Ji et al., 2019b). The final published CAB-NP 1.0 parcellation solution can be visualized via the Brain Analysis Library of Spatial maps and Atlases (BALSA) resource (https://balsa.wustl.edu/rrg5v) and downloaded from the public repository (https://github.com/ColeLab/ColeAnticevicNetPartition, copy archived at swh:1:rev:94772010ac26f487fd6baf1d33d121d57a37e0ed; Cole et al., 2021). The cortex component of the parcellation solution is comprised of 180 bilateral parcels (a total of 360 across both left and right hemispheres), consistent with the Human Connectome Project's Multi-Modal Parcellation (MMP1.0) (Glasser et al., 2016). The subcortex component is comprised of 358 parcels defined using resting-state functional BOLD covariation with the cortical network solution (Ji et al., 2019b).

Behavioral scores (a priori and PCA) were quantified in relation to individual GBC variation (either dense or parcellated) via a mass univariate regression procedure. The resulting map of regression coefficients reflected the strength of the relationship between patients’ behavioral PC score and GBC at every neural location, across all 436 patients. The greater the magnitude of the coefficient for a given location, the stronger the statistical relationship between GBC and the behavioral variation across patients. The coefficients were then Z-scored for each map. Significance of the maps was assessed via nonparametric permutation testing, 2000 random shuffles with TFCE (Smith and Nichols, 2009) type-I error-protection computed via the Permutation Analysis of Linear Models program (Winkler et al., 2014).

Stability of the symptom-neural mapping was assessed following the same cross-validation logic as described for the symptom-driven PCA solutions. Specifically, fivefold cross-validation was performed by first randomly partitioning all patients (N = 436) into five subsets. Regression of the behavioral PC scores onto GBC across patients was performed while holding out 1/5 of the patient sample (N = 349). The correlation between the resulting neural coefficient map was then computed with the neural map obtained from the full sample calculation. For leave-site-out cross-validation, all subjects except those from one site were used when calculating the symptom-neural regression model. The resulting maps were compared to the map from the full sample regression. This process was repeated six times, each time leaving out subjects in one of the six sites.

Reproducibility of the symptom-neural mapping was assessed with 1000 iterations of split-half replication. For each iteration, the full sample of subjects was first randomly split into two halves (referred to as H1 and H2) with the proportion of subjects in each diagnostic group (BPP, SADP, SZP) preserved within each half. For each iteration we used the H1 PCA solution and loadings to compute the predicted PCA scores for H2. In turn, the observed H1 scores were computed from a PCA loadings on the same H1 half-sample of patients. These H1 scores were then regressed against parcellated GBC for patients in H2. This coefficient GBC map reflects the strength of the relationship between the predicted PC score and GBC across H2 patients. Finally, the GBC coefficient maps derived from the H1 observed and H2 predicted PCA scores were then correlated for each PC axis. This process was then repeated 1000 times and evaluated for reproducibility.

To evaluate the consistency of the neural GBC, we computed a PCA solution for the parcellated neural GBC data for all 202 control subjects as well as for all 436 patients (Appendix 1—figure 19). The resulting neural GBC-derived PCs capture the striking consistency of the neural variance components irrespective of clinical status, which highlights that the bulk of the neural variance is not symptom-relevant. As with the behavioral PCA, significance of the neural PCA solution was assessed via permutation testing (1000 random shuffles of parcels within subject). To evaluate the effect of shared neural variance on the symptom-neural solution, we parsed out common neural PCs from the PSD GBC data. We first conducted a PCA using the parcellated GBC data from all 436 PSD and 202 CON (a matrix with dimensions 638 subjects x 718 parcels). This GBC-PCA resulted in 637 independent GBC-PCs. Since PCs are orthogonal to each other, we then partialled out the variance attributable to GBC-PC1 from the PSD data by reconstructing the PSD GBC matrix using only scores and coefficients from the remaining 636 GBC-PCs (${\widehat{GBC}}_{woPC1}$). We then reran the univariate regression as described in Appendix 1—figure 2, using the same five symptom PC scores across 436 PSD. We also repeated the symptom-neural regression with the first two GBC-PCs partialled out of the PSD data (${\widehat{GBC}}_{woPC1-2}$), with the first three PCs parsed out (${\widehat{GBC}}_{woPC1-3}$), and with the first four neural PCs parsed out (${\widehat{GBC}}_{woPC1-4}$).

Canonical correlation analysis (CCA) is a multivariate statistical technique which examines simultaneously the relationships between multiple independent variables and multiple dependent variables by computing linear combinations of each variable set that maximizes the correlations between the two sets (Figure 5A). Here, the two variates used were symptom and neural features across all subjects. Each feature was Z-scored prior to computing the CCA solution. Given the size of the ‘dense’ neural feature space reduced the number of neural features in a principled manner via the described parcellation. This reduced the number of neural features to 180 symmetrized cortical parcels (i.e. GBC was averaged for each pair of analogous parcels in the left and right hemispheres). Critically, corresponding cortical parcels in the left and right hemispheres of this parcellation have been shown to be highly similar (Glasser et al., 2013; Burt et al., 2018). First, we computed the CCA solution using all 36 item-level symptom measures as behavioral features. In turn, we computed an additional CCA solution using scores from the significant five PCs.

Each symptom PC is a weighted linear composite of the original symptom measures and each behavioral CV is a weighted linear composite of the symptom PCs. Therefore, to compute the loadings of the 36 symptom measures on each of the behavioral CVs (computed using the symptom PCs), we multiplied the matrix of loadings from the CCA with the matrix of loadings from the PCA (Figure 5J).

Stability of the CCA solution was assessed via cross-validation. Five-fold cross-validation of the CCA was performed by first randomly partitioning all patients (N = 436) into five subsets. The CCA was then performed between neural and behavioral features using all but one of the subsets. The results of this CCA was then compared to the full sample CCA. Appendix 1—figure 17A-D shows the comparisons across all five five-fold cross-validation runs of the CCA compared to the full model, including neural factor loadings, behavioral factor loadings, and projected symptom measure loadings. For leave-site-out cross-validation, all subjects except those from one site were used in the CCA. The resulting outputs were then compared to those from the full sample CCA. This process was repeated for all six sites. These data are shown in Appendix 1—figure 17E-H.

For leave-one-out cross-validation of the CCA, one subject was held out as a CCA was performed using neural and behavioral features from the other 435 subjects. The loadings matrices $\mathrm{\Psi }$ and $\mathrm{\Theta }$ from the CCA were then used to calculate the ‘predicted’ neural and behavioral latent scores for all five canonical modes for the holdout subject. This was repeated for every subject, such that predicted neural and behavioral latent score matrices (${\displaystyle \hat{N}}$ and ${\displaystyle \hat{B}}$) were computed, of the same dimensions as (${\displaystyle N}$ and ${\displaystyle B}$) respectively. The corresponding CVs in ${\displaystyle \hat{N}}$ and ${\displaystyle \hat{B}}$ were then correlated across subjects, as shown in Figure 5M. If the CCA solution is stable, these correlations should be comparable to the canonical correlations of the full CCA (Figure 5D).

Reproducibility of the CCA solution was assessed via split-half replication. This was performed by first randomly splitting the full sample into two halves (H1 and H2) with the proportion of subjects in each diagnostic group (BPP, SADP, SZP) preserved within each half. A CCA was then performed separately in each half and the resulting outputs were then compared to each other. This process was repeated 1000 times to obtain mean and standard deviation performance metrics. These data are shown in Figure 5L.

Multivariate power analyses to estimate the minimum sample size needed to sufficiently power a CCA were computed using methods described in Helmer et al., 2020, using the Generative Modeling of Multivariate Relationships tool (gemmr, https://github.com/murraylab/gemmr [v0.1.2, copy archived at swh:1:rev:672c2fbb8ff8f5b08afb16fe9b790536c69f2cf3; Helmer, 2021]). Briefly, a model was built by: (1) Generating synthetic datasets for the two input data matrices, by sampling from a multivariate normal distribution with a joint covariance matrix that was structured to encode CCA solutions with specified properties; (2) Performing CCAs on these synthetic datasets. Because the joint covariance matrix is known, the true values of estimated association strength, weights, scores, and loadings of the CCA, as well as the errors for these four metrics, can also be computed. In addition, statistical power that the estimated association strength is different from 0 is determined through permutation testing; (3) Varying parameters of the generative model (number of features, assumed true between-set correlation, within-set variance structure for both datasets) the required sample size n_req is determined in each case such that statistical power reaches 90% and all of the above described error metrics fall to a target level of 10%; and (4) Fitting and validating a linear model to predict the required sample size n_req from parameters of the generative model. This linear model was then used to calculate n_req for CCA in three data scenarios: (i) 718 neural vs. 5 symptom features; (ii) 180 neural vs. 5 symptom features; (iii) 12 neural vs. 5 symptom features.

To illustrate the generalizability of the single-subject prediction results across independent datasets and DSM disorders, we use an independently collected dataset consisting of 30 patients with a formal diagnosis of schizophrenia (SZP) and 39 patients diagnosed with obsessive compulsive disorder (OCD). These patients were recruited via clinician referral and regional advertising and assessed at the Connecticut Mental Health Center in New Haven, CT. Demographic and symptom data are shown in Appendix 1—table 4.

The behavioral assessment for the replication dataset included the PANSS and the Penn Computerized Neurocognitive Battery (PennCNB)(Moore et al., 2015). Items from the PennCNB were matched to items from the BACS as follows: BACS digital sequencing (working memory) to LNB (N-back Total score); BACS symbol coding (attention) -to Penn Continuous Performance Test; BACS token motor task (motor) to CTAP; BACS Tower of London (abstraction and mental flexibility) to PCET or RAVEN (Penn Conditional Exclusion Test); BACS verbal fluency (verbal ability) to SPVRT (Penn Verbal Reasoning Test); BACS verbal memory test to PLLT (Penn List Learning Test) (Gur et al., 2010).

Neural data were collected using a Siemens 3T scanner with a 64-channel head coil at the Yale Center for Biomedical Imaging. Imaging acquisition parameters were aligned with those of the Human Connectome Project (HCP) (Van Essen et al., 2013). High-resolution T1w and T2w structural images were acquired in 224 AC-PC aligned slices, 0.8 mm isotropic voxels. T1w images were collected with a magnetization-prepared rapid gradient-echo (MP-RAGE) pulse sequence (TR = 2400 ms, TE = 2.07 ms, flip angle = 8°, field of view = 256×256 mm). T2w images were collected with a SCP pulse sequence [TR = 3200 ms, TE = 564 ms, flip angle mode = T2 var, field of view = 256×256 mm]. Resting-state BOLD images were collected with a multi-band accelerated fast gradient-echo, echo-planar sequence (acceleration factor = 6, time repetition (TR) = 800 ms, time echo (TE) = 31.0 ms, flip angle = 55°, field of view = 210×210 mm, matrix = 84×84, bandwidth = 2290 Hz); 54 interleaved axial slices aligned to the anterior-posterior commissure (AC-PC) with 2.5 mm isotropic voxels. Additionally, a pair of reverse phase-encoded spin-echo field maps (anterior-to-posterior and posterior-to-anterior) were acquired (voxel size = 2.5 mm isotropic, TR = 7220 ms, TE = 73 ms, flip angle = 90°, field of view = 210×210 mm, bandwidth = 2290 Hz). These neural data were preprocessed using the HCP minimal preprocessing pipeline as described in the section ‘Neural data acquisition and preprocessing’ above, with the same motion scrubbing and nuisance signal regression parameters as were used for the BSNIP dataset.

For the purposes of patient selection, we were focused on individual differences in neural features as they co-vary in relation to individual symptom scores. However, this individual variation in neural features may be small compared to the overall group mean. Hence for each patient, we compute the difference in each brain location (parcel) between the patient’s actual GBC and the group mean GBC for that location. This is denoted by ΔGBC. Importantly, using a de-meaned GBC metric ‘standardizes’ the data and helps to correct for possible differences in scanners/protocols across different datasets. Next, we developed an optimized univariate regression framework leveraging a ‘dot product’ metric to relate a vector of neural features with a single symptom scalar value. This process for neural feature selection (results in Appendix 1—figure 26) is shown as a systems flow diagram in Appendix 1—figure 24. The observed dot product GBC metric ($dpGB{C}^{obs}$) is computed as follows:

• where for a given subject $i$,$dpGB{C}_i^{obs}$ denotes the dot product GBC value of the two vectors ${\displaystyle \mathrm{\Delta }GB{C}_i^{obs}\cdot \beta GB{C}_{N-1}^{obs}}$ across all $P_{select}$ parcels,

• ${\displaystyle \mathrm{\Delta }GB{C}_i^{obs}}$ is a vector of length $P_{select}$ denoting a difference map of subject $i$’s $GB{C}^{obs}$ map relative to the group mean $GB{C}^{obs}$ map, within a given number of parcels $P_{select}$,

• the $\beta GB{C}_{N-i}^{obs}$ vector denotes the PCA-to-GBC statistical group-level β map for a low-dimensional PC symptom score across selected parcels $P_{select}$ for $N$ subjects excluding subject $i$.

This calculation is then repeated for each subject $i$, resulting in a final vector $dpGB{C}^{obs}=[dpGB{C}_1^{obs},\mathrm{\dots },dpGB{C}_N^{obs}]$ for $N$ subjects. There are several key properties of the $dpGB{C}^{obs}$ statistic: (i) it is not inflated by individual GBC map similarity to the group map because each subject’s ${\displaystyle \mathrm{\Delta }GB{C}^{obs}}$ map is demeaned relative to the reference group computed independently of the left-out subject; (ii) this statistic is not biased by the parcel number (which drops with iterative selection) because the resulting $dpGB{C}^{obs}$ value variation is quantified relative to the low-dimensional symptom score across all subjects (see Equation 2). Put differently, the final evaluation considers the relationship between $dpGB{C}^{obs}$ and the low-dimensional PC symptom scores across individuals; (iii) The dot product statistic can yield both positive and negative values – a property which some map similarity measures lack (e.g. ${\eta }^2$); (iv) It is unbounded (unlike a correlation), which is key to maximize co-variation with low-dimensional symptom scores across individuals (see Equation 2); (v) The ${\displaystyle \mathrm{\Delta }GB{C}_i^{obs}}$ map for a given individual is projected onto the basis set of the $\beta GB{C}_{N-i}^{obs}$ map, which is independent of the left-out individual but directly related to the low-dimensional PC symptom score variance, thus maximizing the dot product optimization.

Next, we select $N-i$ individuals and compute a univariate regression where participants’ low dimensional symptom scores are regressed onto the $dpGB{C}_{N-i}^{obs}$ values:

• where for $dpGB{C}^{obs}$ each element denotes the dot product value of the ${\displaystyle \mathrm{\Delta }GB{C}_i^{obs}\cdot \beta GB{C}_{N-1}^{obs}}$ vectors per subject across $P$ parcels,

• α denotes the regression coefficient in the univariate linear model,

• where for $S^{obs}$ each element denotes the observed low-dimensional symptom score (e.g. PC3 score) per subject,

• $ϵ$ denotes the error term in the univariate linear model.

After the $dpGB{C}^{obs}=\alpha S^{obs}+ϵ$ regression is computed on $N-i$ subjects, it is applied to the left-out-subject $i$. This is repeated for all $N$ subjects and the model is evaluated for the number of parcels that maximize two key dot product evaluation metrics (Figure 6B):

Metric A

Metric B

• where $A$ denotes the maximum $r$ correlation value for the two vectors of $dpGB{C}^{pred}$ values and the predicted $S^{pred}$ low-dimensional symptom scores (e.g. PC3 axis) for $N$ subjects from the leave-one-out cross-validation,

• where $B$ denotes the maximum $r$ correlation value for the two vectors of $obs$ and $pred$ $dpGBC$ values for $N$ subjects.

In the initial step in the step-down model, all $P=718$ parcels are retained in the initial dot product calculation. For each iteration of ${\displaystyle P}$ selected parcels, the least predictable parcel ${\displaystyle p}$ (i.e. the parcel with the weakest value in the PC3 map) is eliminated from the map. Then, the step-down regression is repeated until $P=1$.

Methods for the lysergic acid diethylamide (LSD) neuroimaging study are described in detail in prior publications (Preller et al., 2018). The use of LSD in humans was authorized by the Swiss Federal Office of Public Health, Bern, Switzerland. The study protocol was approved by the Cantonal Ethics Committee of Zurich (KEK-ZH_No: 2014_0496). Participants received both written and oral descriptions of the study procedures and the effects and possible risks of LSD. All participants provided written informed consent in accordance with the Declaration of Helsinki. The study employed a fully double-blind, randomized, within-subject cross-over design with three conditions: (1) placebo + placebo (Pla) condition: placebo (179 mg Mannitol and Aerosil 1 mg po) after pretreatment with placebo (179 mg Mannitol and Aerosil 1 mg po); (2) Pla+LSD (LSD) condition: LSD (100 μg po) after pretreatment with placebo (179 mg Mannitol and Aerosil 1 mg po), or (3) Ketanserin+LSD (Ket+LSD) condition: LSD (100 μg po) after pretreatment with the 5-HT2A antagonist Ket (40 mg po). Data were collected for all subjects in a randomized counterbalanced order at three different sessions each two weeks apart. For all conditions, the first substance was administered 60 min before the second substance, and the first neural scan was conducted 75 min after the second administration, with a second scan conducted at 300 min post-administration. In the present study, only data from the two neural scans for the LSD and Pla conditions were evaluated.

Briefly, neuroimaging data acquisition details for the LSD study are as follows. MRI data were acquired on a Philips Achieva 3.0T whole-body scanner (Best, The Netherlands). A 32-channel receiver head coil and MultiTransmit parallel radio frequency transmission was used. Images were acquired using a whole-brain gradient-echo planar imaging (EPI) sequence (repetition time = 2,500 ms; echo time = 27 ms; slice thickness = 3 mm; 45 axial slices; no slice gap; field of view = 240 mm²; in-plane resolution = 3 mm × 3 mm; sensitivity- encoding reduction factor = 2.0). 240 volumes were acquired per resting state scan resulting in a scan duration of 10 mins. Additionally, two high-resolution anatomical images were acquired using T1-weighted and T2-weighted sequences. T1-weighted images were collected via a 3D magnetization-prepared rapid gradient-echo sequence (MP-RAGE) with the following parameters: voxel size = 0.7 mm³, time between two inversion pulses = 3123 ms, inversion time = 1055 ms, inter-echo delay = 12 ms, flip angle = 8°, matrix = 320 × 335, field of view = 224 × 235 mm, 236 sagittal slices. Furthermore T2-weighted images were collected using via a turbo spin-echo sequence with the following parameters: voxel size = 0.7 mm³, repetition time = 2500 ms, echo time = 415 ms, flip angle = 90°, matrix = 320 × 335, field of view = 224 mm × 235 mm, 236 sagittal slices.

Similar to the LSD study, the ketamine pharmacological neuroimaging protocol employed a within-subject design where all healthy volunteer participants underwent a single scanning session consisting of two infusions: (i) placebo (saline solution) followed by (ii) ketamine infusion. The study was approved by the Yale University institutional review board. Participants received both written and oral descriptions of the study procedures and the effects and possible risks of ketamine, and all participants provided written informed consent before beginning the study. Healthy volunteers were informed prior to scanning that they would undergo one placebo run and one ketamine run but were blinded to the order of administration. Because of the sustained effects of ketamine, this infusion was always the second of the two runs, consistent with prior work (Anticevic et al., 2012a). A doctor and two nurses as well as a study coordinator remained present for the duration of the scan. A bolus of ketamine (0.3 mg/kg of bodyweight) or saline were delivered via infusion 5 s after the start of the run and then continuously at a rate of 0.65 mg/kg per hour through the duration of the session. The sequence of scans in each for either the placebo or ketamine infusion was as follows: (i) resting state (4.67 min); (ii) blood draw (sham if saline condition); (iii) a cognitive working memory task (total 14 min); (iv) blood draw (sham if saline condition); (v) a cognitive working memory task (total 14 min); (vii) blood draw (sham if saline condition); (vii) a cognitive working memory task(total 8.63 min). Data from the the cognitive working memory task were not used in the present study and are actively undergoing a distinct analysis. Participants were scanned at Yale University on a Siemens Trio 3T whole-body scanner and 32-channel receiver head coil. High-resolution structural T1-weighted images were acquired using an MP-RAGE sequence and the following parameters: voxel size = 0.8 mm³, time between two inversion pulses = 3123 ms, inversion time = 1055 ms, inter-echo delay = 12 ms, flip angle = 8°, matrix = 320 × 335, field of view = 227 × 272 mm, 227 sagittal slices. T2-weighted images were acquired with the following parameters:voxel size = 0.8 mm³, repetition time = 2500 ms, echo time = 415 ms, flip angle = 90°, matrix = 320 × 335, field of view = 227 × 272 mm, 227 sagittal slices. BOLD images were acquired using a whole-brain gradient-echo planar imaging (EPI) sequence (400 frames at TR = 0.7 ms; TE = 30 ms; slice thickness = 2.5 mm; 54 axial slices; no slice gap; field of view = 250 mm²; in-plane resolution = 2.5 mm × 2.5 mm). In addition, a field map and pair of spin-echo gradients were collected at the end of every scanning session.

All data were preprocessed to be consistent with the BSNIP and replication dataset processing steps. Specifically, we used the HCP minimal preprocessing pipeline and with the same motion scrubbing and nuisance signal regression parameters as were used for the BSNIP dataset as described in the section ‘Neural data acquisition and preprocessing’ above. Parcel-wise GBC maps were computed as described in the section ‘Resting-state functional connectivity and global brain connectivity’ above, by first parcellating the dense time-series for each subject and then computing the parcel-level GBC. As with the BSNIP and replication analyses, parcel GBC was first calculated by parcellating the dense time-series for every subject and then computing GBC across all 718 parcels from the whole-brain functional parcellation (Ji et al., 2019b; Glasser et al., 2016). Group-level GBC maps were computed for every participant for all conditions. We computed the contrast map for ‘LSD-Placebo’ as well as ‘Ketamine-Placebo’ conditions as a t-test between both the pharmacological scans versus placebo scans. The Z-scored t-contrast maps between pharmacological and placebo conditions were used as a pharmacological target map in relation to the PC-derived neural GBC maps (Figure 7 and Appendix 1—figure 29).

The gene mapping analyses in this study utilize the procedure described in Burt et al., 2018. Briefly, we used cortical gene expression data from the publicly available Allen Human Brain Atlas (AHBA, RRID:SCR-007416) , mapped to cortex (Burt et al., 2018). Specifically, the AHBA quantified expression levels across 20,737 genes obtained from six postmortem human brains using DNA microarray probes sampled from hundreds of neuroanatomical loci. Recent studies demonstrated the ability to map expression of each gene onto neuroimaging-compatible templates (Hawrylycz et al., 2012; Burt et al., 2018). Building on these innovations we mapped gene expression on to 180 symmetrized cortical parcels from the HCP atlas (Glasser et al., 2016) in line with recently published methods (Burt et al., 2018). This yielded a group-level map for each gene where the value in each parcel reflected the average expression level of that gene in the AHBA dataset. These group-level maps were in turn as a gene expression target map in relation to the PC-derived neural GBC maps.

Code and core results presented in this paper are available at https://github.com/AnticevicLab/bbs_manuscript/ (copy archived at swh:1:rev:4776b571bbf9a27a2f43146728a6834edbb552ce; Ji, 2021).

We performed an independent component analysis (ICA) on the symptom data as a comparison to our PCA data-reduction solution. Specifically, an ICA was conducted on the behavioral data as an alternative dimensionality reduction method to the PCA. Five components were selected to allow for a direct comparison to be made to the results of the PCA (Appendix 1—figure 3). ICA decomposition was computed on the 436 × 36 (patients x symptom variables) data matrix using Hyvarinen’s (Hyvärinen, 1999) FastICA algorithm, executed via the ‘icafast’ package in R (Helwig and Hong, 2013). Data were mean-centered before decomposition. After decomposition, the source signal estimate (‘ICA score’) for all subjects for each independent component (IC) was regressed on to the neural data across patients, as with a priori and PCA scores, described in the methods section. The proportion of variance explained by each independent component (IC) is shown in Appendix 1—figure 3A. These percentages are comparable to those observed across all the PCs in Figure 1C. Appendix 1—figure 3B shows the correlations between subject scores across all five ICs and five PCs (N = 436). Notably, certain pairs of components between the two solutions appear to be highly similar and exclusively mapped to each other (IC5 and PC4; IC4 and PC5). The correspondence between PCs 1–3 and ICs 1–3 appears to be more distributed as opposed to following a one-to-one mapping, likely because PCs are orthogonal to each other whereas ICs are not necessarily orthogonal. For example, PC3 appears to correspond to IC2, but is negatively related to IC3, suggesting that these two ICs are oblique to the PC and perhaps share behavioral symptom variation in the ICA solution that is explained by a single PC3. Of note, PC3 (the ‘Psychosis Configuration’ axis) is a bi-directional axis whereby individuals who score either highly positively or highly negatively are highly symptomatic in different ways. The diverging nature of this axis may be not be captured well in the ICA solution, when orthogonality is not enforced. Appendix 1—figure 12F shows the correlations across parcels (N = 718 parcels) between the five IC maps (Appendix 1—figure 12A-E) and five PC maps. Again, the neural pattern resembles that in Appendix 1—figure 12, with certain pairs corresponding highly between the the ICA and PCA solutions. The coefficient maps for each of the five ICs are shown in Appendix 1—figure 12. To compare the ICA solution to the PCA solution, we performed a correlation between the individual PCA and ICA scores across subjects, as well as between parcel coefficients of the five PCA maps and five ICA maps across parcels.

Here, we expand on the specific details regarding the cross-validation and reproducibility of the symptom-neural mapping result described in Figure 4. Specifically, Figure 4D shows a summary of a fivefold bootstrapping symptom-neural map validation. Here, the patient sample was randomly split into five folds and each fold was held out once as the regression was performed between symptom scores and parcel-level GBC in the other four folds. The resulting GBC coefficient map was then correlated with the full sample map performed in all 436 patients. Notably, the correlations between maps for each of the k-fold runs and the full model reached effects between r = 0.8–0.9. The coefficient map for one such fold is shown in Figure 4E for and the correlation between the coefficients in this map and the coefficients of the full regression model (performed in all 436 subjects) is shown in Figure 4F. Figure 4G shows a similar summary for leave-one-site-out validation for the regression of all five PC dimensions and traditional symptom scales onto neural GBC. Here, patients from one of the sites were left out while a regression was performed in all remaining patients. In turn, this effect was correlated with the coefficient map from the full patient sample. Again, the maps are highly similar for all sites, showing that the relationships between behavior and neural data are not driven by any one site in particular. Here, the coefficient map for one example is shown in Figure 4H for PC3. Specifically, the results performed without Site three are shown because this is the site with the greatest number of patients and therefore may be the most likely to impact the final effect. The correlation between the coefficients in the resulting map without Site three and the coefficients of the full sample regression model are shown in the scatterplot in Figure 4I. Finally, we performed a stratified split-half replication. Here, the full sample of patients was randomly divided into two halves (H1 and H2) with the proportion of each diagnostic group (BPP, SADP, SZP) preserved within each half. Then, PCA on the symptom data was computed in each half independently. As noted in the main text, this resulted in the ‘observed’ PCA values for H1 and H2. In turn, we used loadings from H1 PCA to compute the ‘predicted’ PC scores for subjects in H2, and vice versa. These ‘predicted’ and ‘observed’ PC scores were then independently regressed against parcellated GBC for subjects a given half. This yielded two maps for each PC: one for the ‘predicted’ PC scores based on H1 data; and one for the ‘observed’ PC scores from H1 data. These two symptom-neural maps were then correlated to evaluate reproducibility of the effects. This process was performed for both H1-to-H2 reproducibility as well as H2-to-H1 reproducibility. We performed 1000 runs of the entire split-half replication analysis for the PC1-5 maps (Figure 4J. Example data for PC3 H1 and H2 are shown in Figure 4K–M). These results highlight that the behavioral PC-to-GBC mapping yields a highly reproducible ‘brain-behavior’ effect.

CCA solves for transformation matrices $\mathrm{\Psi }$ and $\mathrm{\Theta }$ such that the correlations between the linearly transformed ‘latent’' matrices U and V are maximal (see Materials and methods for computational details). Matrix U contains the latent ‘behavioral’ subject scores for all CV and is a linear composite of behavioral feature matrix B scaled by the loadings in $\mathrm{\Psi }$; similarly, V is the latent neural score matrix and is the linear composite of neural features N transformed by $\mathrm{\Theta }$. Each column in U and V is known as a canonical variate (CV); each corresponding pair of canonical variates (e.g. U1 and V1) is referred to as a canonical mode. The correlations between the transformed latent matrices U and V are maximized (Figure 5B), under the constraint that the canonical modes are orthogonal and independent to each other. Hence the correlation between each subsequent pair of CVs is computed from the residuals of the previous pair, until the number of columns in the smaller of the two matrices is reached.

In the main text we show that our CCA solution is not reproducible. We also examined whether the size of the neural feature space affects the CCA solution. To this end, we computed CCA using: (i) 180 symmetrized cortical parcels; (ii) 359 bilateral subcortex-only parcels; (iii) 192 symmetrized subcortical parcels; (iv) 12 functional whole-brain networks, from the recently developed cortico-subcortical CAB-NP parcellation (Glasser et al., 2016; Ji et al., 2019b). We show results for subcortical (192 neural features) in Appendix 1—figure 13, as well as network-level (12 features neural features) in Appendix 1—figure 14. We evaluated CCA significance via permutation testing whereby patient order was randomly shuffled 5000 times (such that the association between behavioral and neural features was violated) and a CCA solution was computed for each shuffle; (see Materials and methods).

For each CCA solution, we examined how much symptom variance can be accounted for by the neural feature variance. Specifically, the correlation between behavioral features B and the latent neural matrix V (i.e. N$\mathrm{\Theta }$) reflects the amount of variance in B that can be explained by the neural CVs in V (Figure 5E). In the main text, Figure 5F shows the proportion of behavioral variance explained by each of the neural CVs in a CCA performed between 180 neural features and all 36 symptom measures. Figure 5G shows the result for the five PC-derived symptom dimensions. Given that the symptom PCs accounted for ∼50.93% of total symptom variance (see Figure 1), the percentages here were scaled to reflect the total proportion of behavioral variance explained by each CV. The amount of neural variance explained by each behavioral CV is shown in Appendix 1—figure 16.

Next, we examined the proportion of symptom variance explained by each neural canonical variate (CV). Because the PCs accounted for ∼50.93% of total symptom variance (Figure 1), the percentages here were scaled accordingly. While the CCA using PC scores has fewer features than the one using item-level symptom measures, each neural CV explained more behavioral variance (Figure 5G versus 5F), suggesting that the PCs may capture neurally-relevant symptom variation with fewer features. Each CV has a pair of associated symptom and neural profiles, which can inform interpretation. This is shown for CV3 in Figure 5H–K. Specifically, the latent scores for an example behavioral canonical variate (CV3) across diagnostic groups normalized to controls are shown in Figure 5H. The profile of PC loadings onto CV3 are shown in Figure 5I. Full results are shown in Appendix 1—figure 16. CCA loadings are computed as the correlation of the latent matrix U (or V) and the input data matrix B (or N). These CCA loadings reflect the amount of variance in the original data that is extracted by the canonical variates. The high negative loadings of CV3 onto PCs 1,4,5 versus the high positive loadings of CV3 onto PC3 in Figure 5I indicate that CV3 captures complex symptom variation. This can also be seen by the projecting 36 symptom measure loadings from a given PC onto CV3 (Figure 5J). Here the most positively loaded items were #7 ‘Delusions’ (PANSS P1), #9 ‘Hallucinations’ (PANSS P3), #21 ‘Somatic Concern’ (PANSS G1) and #29 ‘Unusual Thought Content’ (PANSS G9), while negatively loaded items include #16 ‘Poor Rapport’ (PANSS N3), #19 ‘Lack of Spontaneity’ (PANSS N6), and #28 ‘Uncooperativeness’ (PANSS G8). Conversely, the corresponding neural factor loadings for CV3 are shown in Figure 5K, illustrating GBC patterns associated with the CV3 profile. Data for all five CVs are shown in Appendix 1—figure 15.

Interestingly, the cross-validation of the CCA effect was reproducible when using overlapping cross-validation strategies Appendix 1—figure 17. However, as noted in the main text, the CCA solution was not reproducible with out-sample cross-validation strategies even though the effect passed p<0.05 permutation testing. This highlights the importance of rigorously testing for out-of-sample reproducibility of brain-behavioral effects across independent sub-samples (Helmer et al., 2020).

Here, we outline the details of the split-half replication for the CCA solution. Specifically, the full patient sample was randomly split (referred to as ‘H1’ and ‘H2’ respectively), while preserving the proportion of patients in each diagnostic group. Then, CCA was performed independently for H1 and H2. While the loadings for behavioral PCs and original symptom measures are somewhat similar (mean r$\approx 0.5$) between the two CCAs in each run, the neural loadings were not stable across H1 and H2 CCA solutions. Critically, CCA results did not perform well for leave-one-subject-out cross-validation (Figure 5M). Here, one patient was held out while CCA was performed using all data from the remaining 435 patients. The loadings matrices $\mathrm{\Psi }$ and $\mathrm{\Theta }$ from the CCA were then used to calculate the ‘predicted’ neural and behavioral latent scores for all five CVs for the patient that was held out of the CCA solution. This process was repeated for every patient and the final result was evaluated for reproducibility. As described in the main text, this did not yield reproducible CCA effects (Figure 5M).

Of note, CCA may yield higher reproducibility if the neural feature space were to be further reduced. As noted, our approach was to first parcellate the BOLD signal and then use GBC as a data-driven method to yield a neurobiologically and quantitatively interpretable neural data reduction, and we additionally symmetrized the result across hemispheres. Nevertheless, in sharp contrast to the PCA univariate feature selection approach, the CCA solutions were still not stable in the present sample size of $N=436$. Indeed, a multivariate power analysis (Helmer et al., 2020) estimates that the following sample sizes will be required to sufficiently power a CCA between 180 neural features and five symptom features, at different levels of true canonical correlation ($r_{t}rue$):

• if $r_{true}=0.1$: $n_{required}=103,545$

• if $r_{true}=0.3$: $n_{required}=8,145$

• if $r_{true}=0.5$: $n_{required}=2,497$

To test if further neural feature space reduction may be improve reproducibility, we also evaluated CCA solutions with neural GBC parcellated according to 12 brain-wide functional networks derived from the recent HCP-driven network parcellation (Ji et al., 2019b). Again, we computed the CCA for all 36 symptom measures as well as five PCs (Appendix 1—figure 14). As with the parcel-level effects, the network-level CCA analysis produced significant results (for CV1 when using 36 symptom measure scores and for all 5 CVs when using the five PC-derived scores). Here, the result produced much lower canonical correlations (∼0.3–0.5); however, these effects (for CV1) clearly exceeded the 95% confidence interval generated via random permutations, suggesting that they may reflect the true canonical correlation. We observed a similar result when we evaluated CCAs computed with neural GBC from 192 symmetrized subcortical parcels and 36 symptoms or 5 PCs (Appendix 1—figure 13). In other words, data-reducing the neural signal to 12 functional networks likely averaged out parcel-level information that may carry symptom-relevant variance, but may be closer to capturing the true effect. Indeed, the power analysis suggests that the current sample size is closer to that needed to detect an effect with 12 vs. 5 features:

• if $r_{true}=0.1$: $n_{required}=8,855$

• if $r_{true}=0.3$: $n_{required}=696$

• if $r_{true}=0.5$: $n_{required}=213$

Note that we do not present a CCA conducted with parcels across the whole brain, as the number of variables would exceed the number of observations. However, the multivariate power analysis using 718 neural features and five symptom features estimates that the following sample sizes would be required to detect the following effects:

• if $r_{true}=0.1$: $n_{required}=421,608$

• if $r_{true}=0.3$: $n_{required}=33,166$

• if $r_{true}=0.5$: $n_{required}=10,168$

This analysis suggests that even the lowest bound of ∼10 k samples exceeds the present available sample size by two orders of magnitude.

We tested a hybrid neurobehavioral patient selection strategy for single PC symptom axes by first imposing a PC-based symptom threshold, followed by a target neural similarity threshold driven by the most highly predictive symptom-neural map features. Specifically, the ‘neural similarity prediction index (NSPI)’ computes a patient-specific Spearman’s ρ between that patient’s ${\displaystyle \mathrm{\Delta }GB{C}^{obs}}$ and the group reference ${\beta }_{PC3}GB{C}^{obs}$ map using the maximally predictive $P=39$ parcels (see Appendix 1—figure 27 for whole-brain results and alternative similarity metrics). Appendix 1—figure 26A shows a significant relationship between each patient’s PC3 symptom score (X-axis) and the neural similarity index (Y-axis). In turn, Appendix 1—figure 26B shows binned results, which provides a visual intuition for patient segmentation across both the neural and behavioral indices (Appendix 1—figure 26A, right: binned by $\rho =0.1$ and $PC{3}_{score}=0.5$). For patients at either tail, the neurobehavioral relationship was robust. Conversely, patients with a low absolute PC3 score showed a weak relationship with symptom-relevant neural features. This is intuitive because these individuals do not vary along the selected PC symptom axis.

Appendix 1—figure 26C shows the mean NSPI across subjects within each PC symptom bin along the X-axis. The resulting sigmoid captures that patients exhibit greater neural similarity if their PC symptom scores are more extreme. To evaluate if this relationship can yield a personalized patient selection, we computed the absolute NSPI (Appendix 1—figure 26D). The effect was approximated by a quadratic function, highlighting that patients with extreme PC3 scores (either positive or negative) exhibited a stronger NSPI (i.e. personalized neural effects that strongly resembled the reference ${\displaystyle ({\beta }_{PC3}GB{C}^{obs}{)}^{P=39}}$ map). Appendix 1—figure 26E shows the application of this neurobehavioral selection procedure, demonstrating that PSD patients with extreme PC3 scores (defined at the top/bottom 10th percentile of the ‘discovery’ sample, $+2.17<PC{3}_{score}<-2.41$) exhibit high NSPI values. We observed an inherent trade-off such that if the PC score threshold was raised then neural target similarity confidence goes up, but fewer patients will be selected. In the discovery PSD sample, all patients were accurately selected above the following neurobehavioral thresholds: ${90}^{th}\%tile<PC{3}_{score}<{10}^{th}\%tile$ and $|\rho |>0.4$ (Appendix 1—figure 26E, 34/436 patients selected, green line). We show consistent effects for when the selection was applied to PC5 (Appendix 1—figure 26F; results for all PCs in Appendix 1—figure 28).

To test if the neurobehavioral selection is generalizable, we used an independent cross-diagnostic sample of 30 patients diagnosed with SZP and 39 diagnosed with obsessive-compulsive disorder (OCD) (Appendix 1—figure 26G, see Materials and methods and Appendix 1—table 3 for details). Applying the ‘discovery’ selection thresholds yielded similar results for ∼6% of the cross-diagnostic ‘replication’ sample for PC3 (Appendix 1—figure 26G, full analyses in Appendix 1—figure 27C). Notably, no replication sample patients were selected along the neurobehavioral thresholds for PC5 (Appendix 1—figure 26H). Although there are SZP patients in the replication cross-diagnostic sample, few scored highly on PC5 and none met the neural similarity threshold, emphasizing that not all patients within the same DSM-based diagnosis will exhibit variation along the same neurobehavioral axis. Collectively, these results show that data-driven symptom scores can pinpoint individual patients for whom their neural variation strongly maps onto a target neural reference map. These data also highlight that both symptom and neural information for an independent patient can be quantified in the reference ‘discovery’ BBS using their symptom data alone.

Appendix 1—figure 1

Download asset Open asset

Appendix 1—figure 2

Download asset Open asset

Appendix 1—figure 3

Download asset Open asset

Appendix 1—figure 4

Download asset Open asset

Appendix 1—figure 5

Download asset Open asset

Appendix 1—figure 6

Download asset Open asset

Appendix 1—figure 7

Download asset Open asset

Appendix 1—figure 8

Download asset Open asset

Appendix 1—figure 9

Download asset Open asset

Appendix 1—figure 10

Download asset Open asset

Appendix 1—figure 11

Download asset Open asset

Appendix 1—figure 12

Download asset Open asset

Appendix 1—figure 13

Download asset Open asset

Appendix 1—figure 14

Download asset Open asset

Appendix 1—figure 15

Download asset Open asset

Appendix 1—figure 16

Download asset Open asset

Appendix 1—figure 17

Download asset Open asset

Appendix 1—figure 18

Download asset Open asset

Appendix 1—figure 19

Download asset Open asset

Appendix 1—figure 20

Download asset Open asset

Appendix 1—figure 21

Download asset Open asset

Appendix 1—figure 22

Download asset Open asset

Appendix 1—figure 23

Download asset Open asset

Appendix 1—figure 24

Download asset Open asset

Appendix 1—figure 25

Download asset Open asset

Appendix 1—figure 26

Download asset Open asset

Appendix 1—figure 27

Download asset Open asset

Appendix 1—figure 28

Download asset Open asset

Appendix 1—figure 29

Download asset Open asset

Appendix 1—figure 30

Download asset Open asset