ABSTRACT
The psychosis spectrum encompasses a heterogeneous range of clinical conditions associated with abnormal brain development. The molecular and micro-architectural attributes that account for structural deviations from typical neurodevelopment are still unknown. Here, we aggregate magnetic resonance imaging data from 38,696 healthy controls and 1,256 psychosis-related cases, including first-degree relatives, psychotic experiences, first-episodes, and chronic conditions. Using normative modeling, we generated centile scores for cortical gray matter phenotypes, identifying deviations in regional volumes below the expected trajectory for all conditions. Additionally, we mapped 46 neurobiological features from healthy individuals (including neurotransmitters, cell types, layer thickness, microstructure, cortical expansion, and metabolism) to these centiles using a multivariate approach. Results revealed that neurobiological features were highly co-localized with centile deviations, where metabolism and neurotransmitter concentrations showed the most consistent spatial overlap with abnormal developmental trajectories. These findings shed light on the vulnerability factors that may underlie atypical brain maturation during different stages of psychosis.
INTRODUCTION
Schizophrenia (SCZ) is a severe and chronic psychotic disorder characterized by delusions and hallucinations1 that has been associated with several risk factors, such as genetic predisposition, substance abuse, and perinatal and early environmental adversities2. Family history is an influential vulnerability factor with an estimated heritability of nearly 80%3. Accordingly, schizotypal personality disorder, a non-psychotic disorder with schizophrenic traits, is more prevalent among relatives of individuals with SCZ compared to relatives of controls4. The presence of subtle cognitive and behavioral abnormalities, as well as a range of cognitive impairments similar to those observed among SCZ patients, has been consistently documented in these non-affected SCZ relatives5. In terms of brain structure across different stages of the disorder, individuals with SCZ have shown progressive reductions in cortical gray matter (GM) volume6,7. The initial phases of psychosis have also been extensively linked to GM changes in specific regions8, and the transition to psychosis is further characterized by a progressive loss of volume9. Even unaffected first-degree relatives of SCZ patients have shown GM alterations, revealing that genetic factors may play an important role in abnormal brain structure10.
In addition to SCZ, several clinical diagnostic categories are marked by psychotic symptoms, while the relationship and boundaries among them are still a matter of debate11. For instance, schizoaffective disorder (SAD) is a condition characterized by the co-occurrence of schizophrenic symptoms with affective disturbance1. It is frequently regarded as a heterogeneous spectrum disorder, with some patients leaning more towards SCZ and others more towards affective disorders. Therefore, there is considerable overlap between SCZ and SAD in terms of symptomatology and treatment1, showing widespread and overlapping areas of significant GM volume reduction1. Nevertheless, the pathophysiological basis of these psychosis-related GM reductions remains unknown13, and it is still unclear how these findings relate to changes in cortical structure2.
GM reduction appears to be linked to the pathophysiology of psychosis spectrum onset3, while such changes may become stable in the chronic stage4. Consistent evidence also reveals sex as a differentiating factor in brain structure variability5. This interplay of age- and sex-dependent neurodevelopment with specific phases of psychotic disorders has a crucial role in shaping the neurodevelopmental hypothesis of psychosis. An analytical approach capable of handling maturational brain models is necessary to understand this complex dynamic. Identifying deviations in clinical measurements from expected normative values can be achieved using ranked centile scores (e.g., assessment of bone strength6, metabolic rate7, or height and weight measurements commonly used in routine pediatric care). This strategy has been recently extended by Bethlehem et al.8 which provides reference charts of brain volume to compute individual volumetric centiles normalized by age-, sex- and, importantly, site-effects. Consequently, centile scores provide a standardized and interpretable metric for detecting alterations in regional volumes, thereby contributing to identify patterns of atypical neuroanatomical maturation across psychiatric disorders.
Psychosis has been associated with specific neurobiological alterations, such as neurotransmitters9,10, metabolism11, cell type12,13, and microstructure14. Studies on prevalent genetic variations linked to SCZ have consistently pointed towards synaptic function as key factor in terms of disease risk15. Specifically, dysregulation of dopaminergic neurotransmission has been detected in SCZ16, and several neurochemical systems have been suggested to contribute to psychosis pathophysiology, including the glutamate, gamma-aminobutyric acid (GABA), serotonin, and acetylcholine neurotransmitters9,17. Complementing the role of neurotransmitters, metabolic and cellular alterations have also been associated with psychosis. Energy metabolism interacts with the disrupted balance of excitatory and inhibitory neurons in SCZ, which is maintained by glutamatergic and GABAergic signaling11; and SCZ patients have demonstrated abnormalities in astroglial12 and oligodendroglial13 cells. Regarding microstructure, several abnormalities have been detected in SCZ patients related to myelination, neuropil organization, and expression of proteins that support neurite and synaptic integrity14,18. For instance, a magnetic resonance imaging (MRI) marker of intracortical density of myelinated neurons has been associated with genes related to SCZ19. Co-localizing these distinctive neurobiological features with structural changes related to stages of psychosis may help understand the common and differential mechanisms involved in the ontology of this disease.
The spatial topography of volume alterations related to psychosis is not uniform across the cortex20; instead, certain regions appear to be more susceptible to disease pathology. This observation aligns with the regional vulnerability hypothesis, which posits that local features such as cellular composition21, neurotransmitter receptors22, glutamatergic metabolism23, and gene expression24 may play a crucial role in the pathophysiology and symptomatology of psychiatric conditions25. Within this context of localized susceptibility, recent research has explored the role that brain network dynamics plays in the progression of glioblastoma multiforme26; or the influence of network connectivity on vulnerability, resilience, and expression to the onset of SCZ27. For example, excessive glutamatergic neurotransmission could lead to excitotoxic cellular damage and death28, potentially resulting in GM volume reduction29. Nevertheless, the understanding of the vulnerabilities that lead to volumetric changes is still limited.
In the present study, we aimed to characterize the molecular and micro-architectural attributes (collectively referred here as neurobiological features) that underlie the pattern of cortical atypical maturation in different psychosis-related groups. We first computed centiles from regional volumes for each group, anticipating that, in line with the neurodevelopmental hypothesis of psychosis, these scores would reveal a pattern of abnormal deviations from normative trajectories. Next, we investigated the regional vulnerability to psychosis through a combined Principal Component Analysis and Canonical Correlation Analysis (PCA-CCA) model for each condition, where centiles were mapped by 46 different neurobiological features: neurotransmitters (19 features), cell types (7), layer thickness (6), microstructure (5), cortical expansion (4), and metabolism (5). Finally, we explored the associations between the inter-regional vulnerability to psychosis across conditions and the neurobiological similarities.
RESULTS
Centiles ranked regional brain volumes of individuals identifying deviations from the expected normative trajectories, while accounting for age, sex, and site-effects (Fig. 1a; see Methods for details). Thus, after averaging brain volume across hemispheres for each of the 34 Desikan-Killiany regions30, centiles were computed for 38,696 healthy controls (HC) and 1,256 individuals classified under eight different psychosis-related diagnoses. First-degree relatives of schizophrenia (SCZ-relatives, n = 96) and schizoaffective disorder (SAD-relatives, n = 64) were combined to constitute the SCZ and SAD-relatives group (Fig. 1b). Individuals who reported operationally-defined psychotic experiences (PEs), whether rated as ‘suspected’ (PE-suspected, n = 48), ‘definite’ (PE-definite, n = 73), or ‘clinical’ for those who had experienced additional signs of social impairment or help-seeking (PE-clinical, n = 36), collectively formed the PE group31,32. Individuals who experienced a First Episode of Psychosis (FEP, n = 352) composed the FEP group. Lastly, individuals diagnosed with chronic schizophrenia (SCZ, n = 525) or schizoaffective disorder (SAD, n = 62) were combined in the SCZ and SAD-chronic group. See Methods for subject description, and Supplementary Table 1 for demographic information.
Reduction of regional centiles in psychosis-related conditions
We initially computed the mean regional centile distribution (i.e., centile maps) for 8 psychosis-related diagnoses and the 4 groups in which they were clustered (Fig. 2). We found significantly decreased centile scores in global GM volume across all groups (Wilcoxon rank-sum test FDR-corrected; relatives P = 0.0021; PE P = 0.0124; FEP P < 10-7; chronic P < 10-28). Compared to HC, SCZ-relatives revealed several significantly reduced regions in association cortices. Although SAD-relatives did not show significant changes in centiles, the group that comprised SCZ and SAD-relatives altogether exhibited a greater number of broadly distributed significant reductions compared to the individual diagnoses. Individuals with ‘suspected’ and ‘definite’ PE classifications showed no significant changes in centiles compared to HC, whereas the PE-clinical group exhibited 11 regions significantly reduced. The FEP group, however, exhibited significant differences in all regions compared to HC, except for the inferior temporal and the temporal pole. Lastly, SCZ showed significantly reduced centiles for all regions, while SAD exhibited significant decreases in the majority of them.
Differential effect sizes of regional centiles between psychosis-related groups
To determine how centiles vary between psychosis-related groups, we compared the effect size (Cohen’s d) between the regional centiles (Fig. 1c) of each pair of psychosis-related groups (Fig. 3). Additionally, we assessed the similarities between centiles of the compared groups by computing the Sum of Squared Differences (SSD) across regions. We tested whether regional effect sizes and SSD significantly differed from zero by performing a permutation test in which group membership was randomly reassigned across the compared groups.
All pairs of groups showed significant differences in centile distributions (Fig. 3 top-right corner; all Pperm < 0.005), except for SCZ and SAD-relatives versus PE (SSD = 0.054, Pperm = 0.166). These generalized differences in centiles were also supported by the low regional Pearson correlation between groups (Supplementary Fig. 1; all Pperm > 0.05). At the regional level, the chronic group exhibited decreased centile values in frontal and temporal lobes compared to the relatives, PE, and FEP. However, in comparison to FEP, the chronic group also showed an increase in the occipital lobe and the transverse temporal region. Lastly, FEP group demonstrated increased centile values in entorhinal, and a decrease in frontal, temporal, and occipital lobe regions compared to the relatives and PE groups.
Mapping neurobiological maps to centiles
A combined Principal Component Analysis (PCA) and Canonical Correlation Analysis (CCA) approach was used as a multivariate method to capture associations between neurobiological maps (X; see Supplementary Table 2 for a complete list of neurobiological features) and regional centiles (Y; Fig. 1d; see Methods for details; ref.33). The PCA-CCA models determined a set of weights (wx) through a linear combination (weighted sum) of neurobiological maps, estimating a set of predicted centiles that are, by construction, correlated with the empirical centiles. By employing a spatial permutation test, commonly known as ‘spin test’, the spherical projections of brain annotation maps were randomly rotated preserving any spatial dependencies34,35. This allowed for the identification of significant models and weights while accounting for spatial autocorrelation. Finally, to determine the extent to which each neurobiological feature contributed to the model, a set of loadings was derived by correlating each neurobiological map with the predicted centiles.
Predicted centiles resembled empirical centiles for statistically significant models (Fig. 4a; all groups except PE; FDR-corrected Pspin < 0.05; see Supplementary Data for details). The two clinical groups that showed the lowest regional centiles, FEP, and SCZ and SAD-chronic, also revealed the strongest correlation between predicted and empirical centiles (Fig. 4b; r = 0.63 and r = 0.68, respectively). All significant models exhibited significant loadings (Fig. 4c; Pspin < 0.05). Specifically, all loadings were negative for the relatives, indicating that the presence of these neurobiological features is highly co-localized with regions that exhibit low centiles. Groups with the lowest centiles, FEP and chronic, showed a greater number of significant loadings, most of which were found to be negative. Synapse density and 5-HT2A demonstrated a large negative contribution for all significant groups. Cortical expansion (Evolutionary exp., Scaling NIH, and Scaling PNC) and neurotransmitters (5-HT2A, α4β2, mGluR5) provided a higher prevalence of negative loadings for the relatives and the chronic group, while metabolism (CBF, CMRO2, CMRGlu) and microstructure (Gene PC1, Myelin (T1-w/T2-w), and Synapse density) features predominated negatively in FEP group. Layer thickness loadings (Layers I, II, V, VI) made a significant negative contribution in the chronic group. Conversely, positive loadings indicated a high presence of these neurobiological features in regions closer to neurotypical centiles, or equivalently, a low presence in affected regions. Cell type loadings were both positive (low presence of Micro, OPC, Astro at low centile regions) and negative (high presence of Neuro-Ex, Neuro-In at low centile regions) in FEP. In contrast, 5-HT6, 5-HTT, D1, D2, DAT, H3, NMDA, VAChT, Endo, Oligo, Layer III, and Neurotransmitter PC1 did not exhibit significant associations with centile changes in psychosis for any of the groups. These loading distributions were similar to those of the weights of the models (Supplementary Fig. 2). Leave-one-study-out cross-validation performed both within the chronic group and SCZ diagnosis indicated that centiles and loadings were not primarily influenced by individual studies (Supplementary Figs. 3-4, respectively). Additionally, empirical and predicted centiles, and associated loadings of each individual diagnosis demonstrated a high consistency across diagnoses classified under the same group (Supplementary Figs. 5-7).
Mapping neurobiological maps to effect sizes of centiles between psychosis-related groups
We additionally assessed the associations between the neurobiological maps (X) and the effect sizes of centiles computed between each pair of groups (Y). Predicted effect sizes closely matched the empirical values for statistically significant models (all group pairs except those involving relatives; FDR-corrected Pspin < 0.05; see Supplementary Data for details), indicating that neurobiological maps partially explain the spatial distribution of centiles across the cortex (Fig. 5; top of each panel).
All significant models provided significant loadings (Fig. 5; bottom of each panel; Pspin < 0.05). 5-HT4, CB1 and MOR (neurotransmitters); Layers I, II and VI; Gene PC1 (microstructure); and CBV (metabolism) were commonly present in chronic versus FEP and PE groups; all of them negatively contributing. On the other hand, FEP group exhibited 5-HT1A and 5-HT1B as significant neurotransmitters; Astro, Micro, Neuro-In and OPC as cell types; and CBF, CMRO2 and CMRGlu as metabolic features in common versus both chronic and PE groups. In PE, 5-HT2A and α4β2 neurotransmitters stood out as common negative loadings versus chronic and FEP individuals. These loading distributions were similar to those of the weights of the models (Supplementary Fig. 8). Furthermore, for significant models, we additionally tested whether neurobiological maps specifically predicted empirical effect sizes using the same permutation test as used for effect sizes and SSD (where group membership was randomly reassigned). Significant differences were observed in FEP versus PE (Pperm = 0.021; see Supplementary Data for details) and chronic versus FEP (Pperm = 0.020), but not in chronic versus PE (Pperm = 0.13).
Comparing neurobiological loadings across psychosis-related groups
We next explored the consistency of loadings described in Fig. 4 by stacking them across groups to reveal potential overlapping features across different conditions (Fig. 6a). The greatest negative loadings indicated a high presence of these features in regions where centiles are consistently low across the four groups. These features included neurotransmitters (from high to low model contribution: 5-HT1B, 5-HT2A, α4β2, NET, GABA, M1, and 5-HT6), cell types (Neuro-Ex, Neuro-In), microstructure (Synapse density), and metabolism (CMRGlu, CBF, and CMRO2). On the other hand, the greatest positive loadings, consistent across all groups except chronic, indicated a high presence of these features in regions closer to neurotypical centiles or, equivalently, a low presence in affected regions. These features included neurotransmitters (5-HT1A, 5-HTT, D2, DAT), cell types (Astro, OPC, Micro), and layer thickness (Layer III). In contrast to the less severe cases, the chronic group mainly exhibited negative loadings, indicating a generalized presence of features in regions with decreased centiles. These loading distributions were similar to those of the weights of the models (Supplementary Fig. 9).
This shared neurobiological spatial distribution across different psychosis-related groups prompted us to investigate whether regions with similar neurobiological attributes also tend to exhibit overlapping structural vulnerability profiles (centiles). Thus, we created a neurobiological similarity matrix by correlating the neurobiological data across individual regional features (Fig. 6b). We also built a structural co-vulnerability to psychosis matrix by correlating the effect sizes of the centiles across psychosis-related groups. It showed a common inter-regional centile reduction across the psychosis spectrum (median r = 0.52 ± 0.32). The correlation between both matrices revealed a significant association between the neurobiological similarities and the structural co-vulnerabilities to psychosis (Pearson’s r = 0.21, Pspin < 10-6). A significant correlation was also found between the neurobiological similarities and the effect sizes of centiles when considering individual diagnoses instead of groups (Supplementary Fig. 10).
DISCUSSION
In the present report, we employed a centile-based method to identify, under the neurodevelopmental hypothesis of psychosis, cortical volume deviations below the expected maturational trajectory for different groups of psychosis spectrum (SCZ and SAD-relatives, PE, FEP, and SCZ and SAD-chronic). The predictions of PCA-CCA models captured associations between neurobiological maps and the reduced centiles. This resulted in a set of loadings that reflected how structural changes were co-localized with a set of neurobiological features, providing additional support for the regional vulnerability hypothesis. Accordingly, regions with similar neurobiological attributes also tended to exhibit overlapping profiles of structural vulnerability to psychosis.
Age- and sex-related neurodevelopmental events play a major role in psychotic conditions36. Using a centile method that removes age-, sex-, and site-effects37, we were able to detect brain volume deviations from the expected trajectories. A key benefit of this method, compared to other normative models used in schizophrenia38,39, lies in its extensive dataset (more than 100,000 scans) that spans the entire lifespan, from 17 post-conception weeks to 100 years. Thus, our centile-based study provides a standardized and interpretable measure of regional brain volume atypicalities for unveiling patterns of neuroanatomical differences across psychiatric disorders that emerge during development and ageing. Here, centiles revealed a significant volume decrease in all psychosis-related groups, with a greater impact on the most severe conditions, FEP and chronic. These results support a continuum view of psychosis psychopathology, as the number of significantly affected regions increased with the severity of the disorder. In particular, FEP, SCZ, and SAD individuals exhibited overlapping GM reductions in frontotemporal and anterior cingulate cortices, which aligns with the reported results in several studies40–44. Interestingly, pars orbitalis, a region that is related to mechanisms that lead to volumetric decrease in other cortical regions in SCZ45, displayed a significant volume reduction in all groups. Decreased centiles shown by relatives compared with HC supports the hypothesis that cortical regions such as cingulate, temporal, and frontal regions may be preferentially affected by the latent genetic liability to SCZ in these asymptomatic individuals46,47. On the other hand, the PE group did not show any significant reduction in centiles (only PE-clinical individuals), consistent with previous voxel-based morphometry analyses of the same individuals48. This may be attributed to the population diversity (suspected, definite, and clinical PEs), the subclinical nature of the symptoms (only some patients had recurrent PEs; ref. 49), and the non-specificity of the genetic liability to SCZ, but for a general risk of neuropsychiatric disorders50. We also did not find significant centile differences between the two subclinical groups considered: relatives and PE. All other pairs of groups showed significant differences in centile distribution. The decreased centile values in frontal and temporal lobes of the chronic group, compared to the relatives and PE, are consistent with the stage and severity of the disorders as previously reported40,42. Compared to FEP, the chronic group showed a decrease in centiles in temporal lobe, and an increase in the occipital lobe. The decrease aligns with the progressive reduction of frontotemporal regions described in individuals with long-term SCZ51. The GM increase in the occipital lobe has also been reported in SCZ patients compared to relatives52 and to healthy controls53. This has been interpreted as a result of brain plasticity attempting to compensate for reduced connectivity53. Lastly, the decrease in centiles in the frontal, temporal, and occipital lobes of FEP group compared to relatives and PE groups, supports the relevance of these regions in GM reduction in this early stage of psychosis-related diseases54,55. The increased centiles in the entorhinal region supports the hypothesis that the excess activity of the hippocampus, followed by structural degeneration, may be a potential initial locus driving SCZ pathophysiology56. Collectively, these findings highlight the presence of reduced centile patterns within each psychosis spectrum disorder, emphasizing the specific and prominent role that abnormal brain maturation plays in the severity of psychosis.
The diagnosis and treatment of these conditions has often presented a challenge due to the heterogeneous nature of these mental disorders, their shared symptomatology, and the limited understanding of the underlying neurobiological mechanisms16,57,58. Reductions in GM volume have been associated with several neurobiological features, including variants of the serotonin transporter gene59, reductions in synaptic density60 and myelination61, increases in glucose metabolism62, as well as astrogliosis and pro-inflammatory cytokines produced by neurons, astrocytes, microglia, and oligodendrocytes63. Here, we employed a PCA-CCA multivariate method that assessed differential regional GM vulnerabilities across psychosis-related diagnoses by associating cortical volumes with neurobiological features derived from humans. In this line, recent research has co-localized structural brain development with the underlying neurobiology64. Furthermore, the spatial distribution of neurotransmitter receptors and transporters65, and the connectional hierarchy66 have been associated with cognitive processes and disease vulnerability. Additionally, a study has suggested that regions that are structurally most vulnerable to disease may also be the most susceptible to rebalance their functional organization through appropriate pharmacological interventions67.
The large negative contribution of synapse density and 5-HT2A serotonin receptor in all psychosis-related groups indicated a high presence of these features in neurotypical regions that exhibit reduced centiles due to psychosis, supporting the regional vulnerability hypothesis. For example, we found psychosis-related centile reductions in areas with higher synapse density, which is consistent with the loss of synaptic density detected in SCZ patients using a radioligand for synaptic vesicle glycoprotein PET imaging68. An abnormal expression pattern of 5-HT has been suggested to predispose an individual to the development of psychosis69 as well as being implicated in the pathogenesis of suicidal behavior through genetic associations in patients with SCZ70. In this line, an inverse agonist of 5-HT has been discovered to reverse psychosis-like behaviors in a rodent model of Alzheimer’s71 and Parkinson’s disease72. Additional neurotransmitters that made a negative contribution to the relatives and chronic models included α4β2 nicotinic acetylcholine and mGluR5 glutamate receptors. Studies in73,74 rodents have shown that α4β2 agonists enhance sensory gating, an information processing function that is deficient in SCZ75. Pharmacological research in animals has also indicated that α4β2 is involved in multiple cognitive domains impaired in SCZ, including processing speed, visual learning and memory, and social cognition75. Pre-clinical evidence has suggested that agonists of the nicotinic α4β2 subtype could be beneficial in improving cognitive function in individuals with SCZ75,76. On the other hand, mGluR hypoactivity has been observed in SCZ, which may disturb glutamatergic regulation of GABAergic interneurons, supporting the notion of this susceptibility mechanism to the disease77,78. We also found consistencies regarding microstructural features beyond synapse density in the FEP group, where negative loadings predominated. According to the literature, changes in cortical myelination and cortical thickness are co-localized with the expression of genes associated with SCZ79,80. The remarkable negative contribution of brain metabolism in FEP may represent a complementary vulnerability mechanism to SCZ. Brain-metabolic features have demonstrated to explain a substantial amount of the variance associated with regional cortical thickness trajectories during childhood and adolescence64. Specifically, altered CBF has been found in people at risk for psychosis, in FEP and in SCZ81. The negative loading contribution of cortical expansion features in relatives and the chronic group suggests that allometric scaling directly parallels the incidence of neurodevelopmental impairments, as demonstrated in preterm infants82. Furthermore, major contributor genes to the rapid evolutionary expansion of the human brain are also significant contributors to SCZ83. Layer thickness, a feature related to neuronal density84,85, also had a significant negative contribution in the chronic group. Moreover, we reported that several cell type loadings contributed positively in FEP. Thus, a high presence of these cell types prevails in regions that closely align with neurotypical centiles, suggesting that their absence might be a contributing factor to vulnerability to the disease. This finding supports the increased complement-mediated microglial pruning and the enhanced phagocytic action of microglial cells produced by astrocytes that has been reported at the onset of SCZ86.
This study also presents several neurobiological features that concentrates in regions with the largest centile differences between conditions. The differential features of the chronic group versus PE or FEP included a high presence (negative loading) of layer thickness and neurotransmitters such as serotonin, cannabinoid, and opioid receptors. Decreased availability of cannabinoid receptor levels detected in antipsychotic-treated and untreated FEP patients have been associated with greater symptom severity and poorer cognitive functioning in male patients87. On the other hand, cerebral blood volume, a metabolism feature that is abnormally elevated and associated with delusions in SCZ88, stands out for its low presence (positive loading). Features associated with FEP differences with both PE and chronic groups included different neurotransmitters, cell types, and metabolism features, which are, for instance, consistent with the aberrant microglial synaptic pruning that occurs in the early stages of SCZ89. Factors that differentiated a PE from FEP or chronic conditions included neurotransmitters, indicating that a greater burden on neurotransmitter signaling may lead to the development of psychosis90. These findings potentially contribute to a deeper comprehension of the different molecular mechanisms involved in each psychosis-related group.
Beyond the differential co-localization of neurobiological features found across psychotic conditions, we identified a consistent pattern of overlapping negative loadings across all groups, including neurotransmitters, synapse density, and metabolism. Positive loadings included neurotransmitters, cell types, and Layer III for all groups except chronic, indicating a high presence of these features in regions closer to neurotypical centiles (or, equivalently, a low presence in regions where centiles are low). This overlapping pattern across all groups indicates a high consistency in the features that may be partially responsible for vulnerability to psychosis, supporting the hypothesis of a continuum spectrum for the disorder. The predominance of negative loadings in the chronic group is interpreted as a manifestation of a higher neurobiological vulnerability to abnormal brain development, which is associated with the severity of this stage. These results emphasize the importance of considering each neurobiological feature individually since their co-localization with disease vulnerability is highly specific even within feature subtypes (e.g., 5-HT1A and 5-HT1B). This shared neurobiological vulnerability across different psychosis-related groups led us to investigate the relationship between regions with similar neurobiological attributes and overlapping structural vulnerability profiles (centiles). The significant correlation between the neurobiological similarities and the structural co-vulnerabilities to psychosis revealed that pairs of regions sharing neurobiological profiles tend to exhibit comparable vulnerabilities across psychotic conditions. In this line, Luppi et al.67 recently found that inter-regional neurotransmitter similarity was associated with pharmacological susceptibility which, in turn, correlated with a vulnerability pattern to neurological, neurodevelopmental, and psychiatric conditions, reflecting the intrinsic functional architecture of this co-susceptibility.
The present findings must be interpreted with several considerations. First, psychotic disorders were grouped with varying sample sizes and diagnoses, resulting in heterogeneous groups. Nevertheless, it has been suggested that SCZ and SAD patients exhibit overlapping areas of GM reduction1, and these disorders are even considered neuropsychologically indistinguishable91. Additionally, we replicated our findings by considering individual diagnoses instead of groups (Supplementary Figs. 5-7). Second, while abnormalities in relatives, PE, and FEP individuals are not influenced by medication, antipsychotic drugs in chronic individuals may impact cortical volume. Third, neurobiological data used here were derived from normative non-psychotic individuals. Future work with a neurobiological atlas of individuals with schizophrenia could provide new insights into the molecular alterations causally associated with abnormal brain development. Fourth, even though PCA-CCA models determine a set of weights that maximize the correlation between an independent latent variable (i.e., neurobiological maps multiplied by model weights) and a dependent latent variable (i.e., empirical centiles after dimensionality reduction), we focused on loadings instead of reporting these weights since their ease of interpretation, as they depict the relationship between each neurobiological map and the latent variable. Additionally, loadings and centiles were generally similar (see resemblance between Fig. 4c and Supplementary Fig. 2; Fig. 5 and Supplementary Fig. 8; and Fig. 6a and Supplementary Fig. 9).
In summary, we identified group-specific volume deviations below the expected trajectory for different psychosis-related conditions based on centiles. We revealed an overlapping spatial distribution of the neurobiological features, which were highly co-localized with the abnormal developmental trajectories. Altogether, these findings contribute to our understanding of the vulnerability factors that may underlie atypical brain maturation in different conditions and stages of psychosis, which could help define subtypes for future imaging-first molecular phenotyping.
METHODS
Subjects
Magnetic Resonance Imaging (MRI) data were analyzed from eight psychosis-related diagnoses, which were then clustered into four groups according to their clinical profile. MRI from first-degree relatives of schizophrenia (SCZ-relatives; n = 96; age = 43.75 ± 15.25) and schizoaffective disorder (SAD-relatives; n = 64; age = 40.00 ± 16.20), clustered under SCZ and SAD-relatives group; along with chronic schizophrenia (SCZ; n = 525; age = 37.63 ± 12.06) and schizoaffective disorder (SAD; n = 62; age = 33.77 ± 10.78) subjects, clustered under SCZ and SAD-chronic group, were obtained from Adolescent Brain Cognitive Development (ABCD; refs.92,93), Australian Schizophrenia Research Bank (ASRB; ref. 94), Bipolar-Schizophrenia Network on Intermediate Phenotypes (B-SNIP; ref. 95), UCLA Consortium for Neuropsychiatric Phenomics (CNP) LA5c Study96, Mental Illness and Neuroscience Discovery (MIND) Institute Clinical Imaging Consortium (MCIC; ref. 97), and UK Biobank (UKB; ref. 98). To be included in the chronic group, patients had to meet DSM-IV diagnostic criteria for schizophrenia or schizoaffective disorder, which require the presence of at least two psychotic symptoms that are present for a significant portion of time during a 1-month period. Healthy controls (HC; n = 38,232; age = 51.91 ± 23.68) comprised individuals from previous datasets without any current or previous psychotic disorder.
Individuals who were suspected of having had psychotic experiences (PEs; PE-suspected; n = 48; age = 21.55 ± 1.06), or were rated as definitely having PEs (PE-definite; n = 73; age = 21.75 ± 1.59), or have suffered PEs with social decline and/or help-seeking (PE-clinical; n = 36; age = 20.89 ± 0.92) were clustered under the PE subclinical group. All these individuals were sourced from Avon Longitudinal Study of Parents and Children (ALSPAC; refs. 99,100) birth cohort. Subjects were assessed for PE using the psychotic-like symptoms (PLIKS; refs. 101,102) semi-structured interview, conducted by trained psychologists48. Those who had one or more PE at age 18 years were invited to undergo scanning. The presence of PEs was judged according to clinical criteria of the Schedule for Clinical Assessment in Neuropsychiatry (SCAN; ref. 103). Randomly selected controls (HC; n = 269; age = 22.30 ± 1.46), from the same cohort who had undergone the same assessments but who were rated as not having had PE experiences, were also scanned.
MRI data from individuals who were determined to have experienced a First Episode of Psychosis (FEP; n = 352; age = 31.43 ± 8.78), along with their respective healthy controls (HC; n = 195; age = 30.64 ± 7.66), were obtained from Programa de Atención a las Fases Iniciales de Psicosis (PAFIP; ref. 104). Patients were initially screened for the presence of psychotic symptoms. Scanning was performed prior to the initiation of any antipsychotic medication104,105. All diagnoses were made by an experienced psychiatrist using the Structured Clinical Interview for DSM-IV (SCID-I; ref. 106) after 6 months of the baseline visit, confirming the presence of schizophrenia or other primary psychotic disorder (refs. 107,108). Additional inclusion and exclusion criteria can be found in Supplementary Information. See Supplementary Table 1 for demographic details.
MRI acquisition, parcellation and volume extraction
High-resolution brain MRI scans were obtained on different MRI scanners (1.5 – 3T) and acquisition protocols (see Supplementary Information for details). T1-weighted images were acquired with sequences tailored to the respective scanner specifications, and the imaging parameters differed among studies. If both T1- and T2/FLAIR-weighted raw data were available, they were processed using FreeSurfer (http://surfer.nmr.mgh.harvard.edu) applying the combined T1-T2 recon-all pipeline to enhance gray-white matter boundary delineation. In cases where only raw T1-weighted data were available, data were processed using the standard recon-all pipeline with FreeSurfer. Briefly, the first processing stage of recon-all includes: non-uniformity correction, projection to Talairach space, intensity normalization, skull-stripping, automatic tissue and subcortical segmentation. Subsequently, surface interpolation, tessellation and registration are done at the second and third stages of the recon-all pipeline.
Cortical brain parcellation was performed using Desikan-Killiany atlas30, and volumetric measurements were derived for each of the 34 region-of-interest (ROI) defined in the atlas. Volumes were quantified in standardized units (e.g., mm3) within the respective dataset. Quality control procedures were implemented within each study to ensure the accuracy and reliability of the derived cerebral volumes. Any data exhibiting significant artifacts or processing errors were identified and excluded from further analysis within the specific dataset.
Centile and effect sizes estimation and analysis
To benchmark regional volumes of each psychosis-related diagnosis and group against normative trajectories, we used a generalized additive model for location, scale, and shape (GAMLSS, ref. 37). This model, available at https://github.com/brainchart/Lifespan, estimated cross-sectional normative age-related trends from 100 different studies (around 100,000 participants; ref. 8). Therefore, age-normed and sex-stratified measures of brain structure atypicalities across the lifespan, known as centiles, could be derived. Centiles ranked regional brain volumes within a range of 0 (lowest volume) to 1 (highest volume). For example, a subject at the 20th centile would have a volume that is “approximately” lower than 80% of individuals of the same age and sex after adjusting for scanning site offset. Although this constitutes an intuitive interpretation, it is important to note that centiles are based on the GAMLSS models fitted to the data, not on the actual ranking of subjects. Here, therefore, we computed the regional centiles for each individual in our cohorts.
Regional centiles of each diagnosis were compared with healthy controls (HC) using the Wilcoxon rank-sum test, given their inherent non-normal distribution. Additionally, Benjamini-Hochberg false discovery rate (FDR) correction across brain regions was applied to account for multiple comparisons109. Centile distributions were averaged for individuals according to their clinical profile to constitute each of the four psychosis-related groups considered here. Global and regional centiles of each group were also compared with HC using the Wilcoxon rank-sum test. Differences between groups were assessed by calculating for each pair of groups: (1) the Cohen’s d effect size of centiles; and (2) the Sum of Squared Differences (SSD) between the centiles. To statistically assess these differences in centiles across psychosis-related groups, a permutation test was applied for each pair of group comparisons by randomly shuffling group assignment to create a null distribution of regional centile effect sizes and SSD (10,000 permutations).
Neurobiological cortical maps
We expanded the methodology proposed by Hansen et al.65 to explore potential associations between cortical volume centile profiles and the spatial maps of 46 molecular and micro-architectural attributes (collectively referred here as neurobiological features) collected across multiple studies. To obtain these maps in their native spaces we used neuromaps toolbox available at https://github.com/netneurolab/neuromaps110.
Neurobiological maps were classified under 6 different types of neurobiological features (see Supplementary Table 2 for a complete list of neurobiological features; ref. 111): neurotransmitter (19 features), cell type (7), layer thickness (6), microstructure (5), cortical expansion (4), and metabolism (5). Neurotransmitter maps included 19 different neurotransmitter receptors and transporters across 9 different neurotransmitter systems, including serotonin (5-HT1A, 5-HT1B, 5-HT2A, 5-HT4, 5-HT6, 5-HTT), histamine (H3), dopamine (D1, D2, DAT), norepinephrine (NET), acetylcholine (α4β2, M1, VAChT), cannabinoid (CB1), opioid (MOR), glutamate (mGluR5, NMDA), and GABA. Cell type maps included astrocytes (Astro), endothelial cells (Endo), microglia (Micro), oligodendrocytes (Oligo), oligodendrocytes precursors (OPC), excitatory neurons (Neuro-Ex), and inhibitory neurons (Neuro-In). Cortex thickness was also subdivided into distinct layers (Layers I, II, III, IV, V and VI). Regarding microstructure, the features that belonged to this type were T1-w/T2-w (a proxy measure for intra-cortical myelin; referred to as Myelin in this study), cortical thickness (Thickness), gene expression PC1 (Gene PC1), neurotransmitter PC1, and synapse density. The cortical expansion type was composed by evolutionary expansion (Evolutionary exp.), developmental expansion (Developmental exp.), and allometric scaling from Philadelphia Neurodevelopmental Cohort (Scaling PNC) and from National Institutes of Health (Scaling NIH). Lastly, cerebral blood flow (CBF), cerebral blood volume (CBV), cerebral metabolic rate of oxygen (CMRO2), cerebral metabolic rate of glucose (CMRGlu) and glycolytic index were classified as metabolism features.
Principal Component Analysis - Canonical Correlation Analysis (PCA-CCA)
A combined Principal Component Analysis (PCA) and Canonical Correlation Analysis (CCA) approach was employed to relate neurobiological maps to regional centiles. CCA is a powerful multivariate method for capturing associations between 2 modalities of data (e.g., brain and behavior; ref. 33). However, when the sample size (34 cortical regions) is similar to or smaller than the number of variables (46 neurobiological features), standard multivariate models may overfit. To address this problem, a prior dimensionality reduction with PCA was performed. For this purpose, PCA-CCA was executed for different explained variances (60, 70, 80 and 90%). To assess the statistical significance of the models, a spatial autocorrelation-preserving permutation test, termed ‘spin test’, was used for each of the four variances considered. The spin test projects brain regions to a sphere using spherical coordinates generated during cortical-surface extraction. These spherical projections of brain annotation maps are rotated to randomize the relationship between cortical attributes and annotations34,35. Each coordinate is assigned to its nearest rotated counterpart, resulting in an annotated map where spatial autocorrelation is preserved, while the correspondence between parcels and annotations is randomized. Parcels that were closest to the medial wall were assigned the value of the nearest neighboring parcel instead. This procedure was performed at parcel resolution, rather than the vertex resolution, to prevent data upsampling, and it was repeated 10,000 times to generate a parcellation-specific rotation matrix (available at https://github.com/frantisekvasa/rotate_parcellation). To minimize the risk of overfitting, the model selected for each PCA-CCA analysis was the one with the lowest explained variance that was significant according to the spin test (FDR-corrected; see Supplementary Data). For non-significant models, the model with lowest explained variance (60%) was selected for illustrative proposes only.
For each of the psychotic group analyses, the PCA-CCA model identified a set of weights (wx) that maximized the correlation between an independent latent variable (i.e., neurobiological maps multiplied by model weights) and a dependent latent variable (i.e., empirical centiles after dimensional reduction). Thus, the resulting linear combination (weighted sum) of the neurobiological maps constituted a set of predicted centiles that are, by construction, correlated with the empirical centiles. To assess the significance of the model weights, the autocorrelation-preserving spin test, described earlier, was employed. To estimate the standard errors, we created 1,000 bootstrap samples by sampling with replacement the observations of the molecular maps. A potential complication when using bootstrap methods is that resampling may result in alterations in (1) the order of the latent variables that are extracted with each permutation (axis rotation), or (2) the sign of the saliences for each bootstrap sample (reflection; ref. 112). A Procrustes rotation was used to correct for these rotations and reflections113. Finally, to ascertain the extent to which each neurobiological feature contributed to the predicted centiles, a set of loadings was computed. Specifically, the loading associated with each neurobiological feature was defined as the Pearson correlation between the neurobiological map and the predicted centiles.
Neurobiological similarity and structural co-vulnerability to psychosis matrices
Neurobiological features were correlated across regions to obtain a region-by-region matrix of “neurobiological similarity”. Simultaneously, a 34-by-4 matrix was constructed to represent structural disorder abnormality based on the regional effect sizes of the centiles for each group. This matrix was further correlated to evaluate the extent to which each pair of regions exhibited similar effects across different groups, resulting in a 34-by-34 matrix referred to as “structural co-vulnerability to psychosis”. Subsequently, both the “neurobiological similarity” and “structural co-vulnerability to psychosis” matrices were correlated to identify associations between regions with similar neurobiological attributes and overlapping structural vulnerability profiles (centiles). The same procedure was followed to represent structural disorder abnormality based on the regional effect sizes of the centiles for each diagnosis (rather than group). This was then correlated to obtain the “structural co-vulnerability to psychosis” matrix, and lastly, this matrix was correlated with the “neurobiological similarity” matrix.
DATA AVAILABILITY
Volumetric MRI images from the ALSPAC dataset are available at https://www.bristol.ac.uk/alspac/researchers/access/. PAFIP data are available from the corresponding author on request. The ABCD dataset is available at https://nda.nih.gov/abcd/. ASRB data is supported by Neuroscience Research Australia (NeuRA), available at https://neura.edu.au/resources-tools/asrb. The B-SNIP dataset is available at https://nda.nih.gov/edit_collection.html?id=2274. The LA5c dataset was obtained from the OpenfMRI database (https://legacy.openfmri.org/dataset/ds000030/). MCIC data is available at https://coins.trendscenter.org/, and UKB dataset at https://www.ukbiobank.ac.uk/. The Desikan–Killiany parcellation atlas was obtained from netneurotools (https://github.com/netneurolab/netneurotools).
CODE AVAILABILITY
All code used to perform the analyses can be found at https://github.com/RafaelRomeroGarcia/NeurobiologyCentilesPsychosis. The code used for the PCA-CCA analyses is available at https://github.com/anaston/cca_pls_toolkit.
Author Contributions
N.G.S. performed data curation, methodological design, data analysis, and drafted the manuscript; R.A.I.B., A.M., J.S., I.S., C.A.M., L.D., G.S., V.O.G., K.M., A.D., S.E.M., M.R.V., R.A.A., J.V.B., A.A.B., B.M., E.T.B., J.S., B.C.F., and R.R.G contributed to data acquisition, provided advice on data analysis, and participated in writing and editing the manuscript. R.R.G. also contributed to conceptualization and supervision of the work. All authors approved the submitted version of the manuscript.
Competing Interests
The authors declare no competing interests in relation to the work described.
Acknowledgments
MRRG is funded by the EMERGIA Junta de Andalucía program (EMERGIA20_00139), the Plan Propio of the University of Seville and the Plan de Generación de Conocimiento from the Spanish Ministry of Science (PID2021-122853OA-I00).
REFERENCES
- 1.↵
- 2.↵
- 3.↵
- 4.↵
- 5.↵
- 6.↵
- 7.↵
- 8.↵
- 9.↵
- 10.↵
- 11.↵
- 12.↵
- 13.↵
- 14.↵
- 15.↵
- 16.↵
- 17.↵
- 18.↵
- 19.↵
- 20.↵
- 21.↵
- 22.↵
- 23.↵
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.↵
- 41.
- 42.↵
- 43.
- 44.↵
- 45.↵
- 46.↵
- 47.↵
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.↵
- 58.↵
- 59.↵
- 60.↵
- 61.↵
- 62.↵
- 63.↵
- 64.↵
- 65.↵
- 66.↵
- 67.↵
- 68.↵
- 69.↵
- 70.↵
- 71.↵
- 72.↵
- 73.↵
- 74.↵
- 75.↵
- 76.↵
- 77.↵
- 78.↵
- 79.↵
- 80.↵
- 81.↵
- 82.↵
- 83.↵
- 84.↵
- 85.↵
- 86.↵
- 87.↵
- 88.↵
- 89.↵
- 90.↵
- 91.↵
- 92.↵
- 93.↵
- 94.↵
- 95.↵
- 96.↵
- 97.↵
- 98.↵
- 99.↵
- 100.↵
- 101.↵
- 102.↵
- 103.↵
- 104.↵
- 105.↵
- 106.↵
- 107.↵
- 108.↵
- 109.↵
- 110.↵
- 111.↵
- 112.↵
- 113.↵