First-episode psychosis (FEP) is characterised by heterogeneous clinical presentations, ^{Reference Voineskos, Jacobs and Ameis1} temporal progression ^{Reference Prakash, Chatterjee, Srivastava and Chauhan2} and responses to medication. ^{Reference Leucht, Crippa, Siafis, Patel, Orsini and Davis3} Cross-sectional magnetic resonance imaging (MRI) studies consistently report alterations in frontal and temporal cortices, ^{Reference Oertel-Knöchel, Knöchel, Rotarska-Jagiela, Reinke, Prvulovic and Haenschel4} with some identifying changes in insular, ^{Reference Crespo-Facorro5} parietal ^{Reference Yildiz, Borgwardt and Berger6} and occipital regions. ^{Reference Tordesillas-Gutierrez, Ayesa-Arriola, Delgado-Alvarado, Robinson, Lopez-Morinigo and Pujol7} Longitudinal studies have additionally revealed progressive cortical deterioration in FEP, especially in frontal and temporal regions, over 1, ^{Reference Del Re, Stone, Bouix, Seitz, Zeng and Guliano8} 2–3, ^{Reference Cobia, Smith, Wang and Csernansky9} 5 ^{Reference van Haren, Hulshoff Pol, Schnack, Cahn, Mandl and Collins10} and 10 years. ^{Reference Rodriguez-Perez, Ayesa-Arriola, Ortiz-García de la Foz, Setien-Suero, Tordesillas-Gutierrez and Crespo-Facorro11} Despite these findings, the complex interplay between brain development, gender and neurodevelopmental pathophysiology complicates the disentanglement of their individual contributions to brain structure and function. Normative modelling addresses this by creating ‘growth charts of brain development, ^{Reference Bethlehem, Seidlitz, White, Vogel, Anderson and Adamson12} allowing individuals to be ranked and assigned centile scores based on deviations from expected neurotypical trajectories. This approach has been used to identify abnormal structural subtypes in autism ^{Reference Shan, Uddin, Xiao, He, Ling and Li13} and deviations in attention-deficit hyperactivity disorder symptoms. ^{Reference Marquand, Rezek, Buitelaar and Beckmann14} In psychosis, normative modelling reveals deviations in cortical thickness in temporal areas during high-risk states ¹⁵ and across the cortex in FEP, ^{Reference Worker, Berthert, Lawrence, Kia, Arango and Dinga16} with these deviations diminishing during treatment. ^{Reference Berthet, Haatveit, Kjelkenes, Worker, Kia and Wolfers17}

Brain structure and clinical outcomes

MRI studies have also extensively explored the clinical presentation of FEP, linking positive symptoms, such as hallucinations and delusions, to thinning in frontal ^{Reference Padmanabhan, Tandon, Haller, Mathew, Eack and Clementz18} and temporal ^{Reference Walton, Hibar, van Erp, Potkin, Roiz‐Santiañez and Crespo‐Facorro19} cortices. Conversely, negative symptoms, including anhedonia and social withdrawal, are associated with reduced orbitofrontal thickness. ^{Reference Walton, Hibar, van Erp, Potkin, Roiz-Santiañez and Crespo-Facorro20} The effects of antipsychotics on brain structure remain debated, with studies reporting both grey matter reductions ^{Reference Vita, De Peri, Deste, Barlati and Sacchetti21} and increased regional thickness. ^{Reference Goghari, Smith, Honer, Kopala, Thornton and Su22} Cognitive decline in FEP has been observed prior to the onset of psychosis ^{Reference Bora and Murray23} while post-onset trajectories vary, with reports of stabilisation, ^{Reference McCleery and Nuechterlein24} decline ^{Reference Fett, Velthorst, Reichenberg, Ruggero, Callahan and Fochtmann25} and improvement. ^{Reference Burgher, Scott, Cocchi and Breakspear26} Normative modelling has shown neurocognitive delays in youth reporting psychotic symptoms, ^{Reference Gur, Calkins, Satterthwaite, Ruparel, Bilker and Moore27} but the long-term impact of antipsychotics on cognition remains uncertain. ^{Reference Husa, Moilanen, Murray, Marttila, Haapea and Rannikko28}

Regional vulnerability to atypical development

The spatial patterns of cortical deterioration shared across psychiatric and neurological conditions have been interpreted as a reflection of differential regional vulnerability to pathology. ^{Reference Hettwer, Larivière, Park, van den Heuvel, Schmaal and Andreassen29} Under this hypothesis, cortical vulnerability to psychosis-related deterioration may reflect regional differences in cellular composition, neurotransmitter receptors and metabolism, ^{Reference Hettwer, Larivière, Park, van den Heuvel, Schmaal and Andreassen29} with regions of high serotonin and acetylcholine receptor density showing greater atypical deviations. ^{Reference García-San-Martín, Bethlehem, Mihalik, Seidlitz, Sebenius and Alemán-Morrillo30} In this study, we hypothesise that FEP patients exhibit atypical brain development influenced by treatment and medication, which impacts cognitive and symptom progression. We analysed baseline and 10-year follow-up longitudinal data to explore associations among brain deviations, clinical diagnosis, medication, cognition and symptoms. Using cytoarchitectonic and neurobiological atlases, we additionally characterised psychosis-related brain deviations based on neurotransmitters, cell types, microstructure and metabolism.

Subjects

MRI and phenotypic data were collected from healthy controls (n = 195; 120 males, mean age 29.1 ± 7.63 years) and individuals with FEP (n = 357; 213 males, mean age 29.8 ± 8.76 years), who were followed longitudinally. Follow-up assessments were conducted at baseline (193 controls, 333 FEP) and at 1 (62 controls, 96 FEP), 3 (51 controls, 136 FEP), 5 (76 controls, 70 FEP) and 10 (91 controls, 101 FEP) years (see Supplementary Figs 1 and 2 available at https://doi.org/10.1192/bjp.2025.10482 for flow diagram for study participants and alluvial diagram of attrition).

Recruitment was conducted through the Program for Attention to the Initial Phases of Psychoses (PAFIP) at Marqués de Valdecilla University Hospital and Valdecilla Health Research Institute (IDIVAL) in Spain. The authors assert that all procedures contributing to this work comply with the ethical standards of the relevant national and institutional committees on human experimentation, and with the Helsinki Declaration of 1975 as revised in 2013. All procedures involving human subjects were approved by CEIC-Cantabria (clinical trial nos NCT0235832 and NCT02534363) and the Ethics Committee of the University of Seville (no. 2024-2534). Following a comprehensive explanation of the study, each participant provided written informed consent prior to enrolment. Participants in the clinical group met DSM-IV criteria for schizophrenia spectrum disorders, including schizophrenia (n = 156), schizophreniform disorder (n = 99), schizoaffective disorder (n = 6), brief psychotic disorder (n = 48) and psychosis not otherwise specified (n = 24). The full list of inclusion and exclusion criteria is provided in Supplementary Methods. Detailed demographic and clinical information, stratified by gender and diagnosis, are provided in Supplementary Figs 3 and 4 and Supplementary Tables 1–4).

Cognitive and symptom assessment

Cognition was assessed for the following domains: (a) verbal memory: Rey Auditory Verbal Learning Test long-term recall score; (b) visual memory: Rey Complex Figure test long-term recall score; (c) motor dexterity: grooved pegboard, time to complete with dominant hand; (d) executive functions: Trail Making Test part B; (e) working memory: WAIS III-Backward Digits total score; (f) speed of processing: WAIS III-Digit Symbol standard total score; and (g) attention: Continuous Performance Test Degraded-Stimulus, the total number of correct responses. ^{Reference Rodríguez-Sánchez, Setién-Suero, Suárez-Pinilla, Van Son, Vázquez-Bourgon and López31} For each cognitive test, z-scores were computed by subtracting the mean and dividing by the standard deviation of the control group’s scores. Overall cognitive functioning was calculated as the mean z-score across the seven cognitive tests for each individual. FEP participants were also assessed using the Scale for the Assessment of Negative Symptoms (SANS, 21 items), Scale for the Assessment of Positive Symptoms (SAPS, 30 items) and Brief Psychiatric Rating Scale (BPRS, 24 items).

Principal component analysis (PCA) was conducted on cognition and symptom items for comparison with average scores used in the calculations. Sensitivity analyses examined the impact of cannabis use and study dropout on cognitive and symptom measures. Psychological evaluations were performed by experienced psychiatrists using the structured clinical interview for DSM-IV, with diagnoses confirmed at 6 months. All statistical analyses were performed using RStudio version 2023.9.1.494 for Windows (RStudio, Boston, Massachusetts, USA; https://www.rstudio.com/products/rstudio/download/).

MRI acquisition and volume extraction

Structural MRI scans were acquired using a 1.5T General Electric SIGMA System (GE, Milwaukee, USA) and a 3T Philips Medical Systems MRI scanner (Achieva, Best, The Netherlands) equipped with an 8-channel head coil. Baseline and follow-up scans for each individual were consistently conducted using the same machine. The imaging protocol included a T1-weighted image that underwent visual inspection for quality control. The parameters for the 1.5T system were: echo time 5 ms, repetition time 24 ms, NEX 2, rotation angle 45°, FOV 26 × 19.5 cm², slice thickness 1.5 mm, voxel size 1.02 × 1.02 × 1.5 mm³ and matrix 256 × 192; and for 3T: echo time 3.7 ms, repetition time 8.2 ms, flip angle 8°, acquisition matrix 256 × 256, voxel size 0.94 × 0.94 × 1 mm³ and 160 contiguous slices. The longitudinal pipeline of FreeSurfer 7.0.0 (recon-all) was used to estimate the global cortical volume and regional volume of each region defined in the Desikan–Killiany atlas (34 regions per hemisphere). One regional outlier (>5 s.d. from the mean) was excluded. Raw regional volumetric data from both scanners were harmonised using ComBatLS, preserving age, gender and diagnosis variance. ^{Reference Gardner, Shinohara, Bethlehem, Romero-Garcia, Warrier and Dorfschmidt32}

Regional volumetric centile estimation

We leveraged the cross-sectional normative trajectories established by the BrainChart study, which benchmarked cortical volumes from the Desikan–Killiany atlas against data from over 100 000 neurotypical individuals. ^{Reference Bethlehem, Seidlitz, White, Vogel, Anderson and Adamson12} Based on these trajectories, we computed out-of-sample centile estimates for our longitudinal sample, using the control cohort to account for scanner-related offsets (https://github.com/brainchart/Lifespan). Following the pipeline of Bethlehem et al, ^{Reference Bethlehem, Seidlitz, White, Vogel, Anderson and Adamson12} the volumetric values of left and right homologous regions were averaged to obtain a single measure per region, and centiles were subsequently computed from these averaged values. Age-normed, gender-stratified and site-corrected deviations were derived and expressed as centile scores; centile scores represent the relative position of an individual’s grey matter volume with respect to normative lifespan trajectories. Scores below 0.5 indicate volumes lower than the median for age- and gender-matched individuals, whereas scores above 0.5 indicate volumes higher than the median. Residual effects of site or motion (Euler index) on centile scores were assessed using analysis of variance (ANOVA).

Baseline analyses using multiple regression

To evaluate the impact of FEP diagnosis on baseline centiles, a multiple regression was performed for each cortical region using the following model:

(1) $$\small{\rm{regional\ centile\sim 1 + diagnosis + gender+ age\ of\ inclusion + etiv}}$$

where diagnosis is a binary variable (control versus FEP), age of inclusion refers to participants’ age at the first scan and etiv (estimated total intracranial volume) represents cranial volume residuals following age and gender correction (which allowed us to account for anthropometric factors using a single variable). Age of inclusion and etiv were z-scored prior to regression. All P-values were corrected for multiple comparisons across regions using false discovery rate (FDR).

Longitudinal analyses using mixed modelling

The longitudinal effects of disease and medication on centiles and clinical outcomes were analysed using mixed modelling over 1209 assessments (195 controls, 357 FEP). Antipsychotic medication exposure was assessed at each visit (Supplementary Fig. 5) and converted to chlorpromazine equivalents (CPZ-equivalents). ^{Reference Danivas and Venkatasubramanian33}

LMM for prediction of regional centiles as a function of disease progression and treatment

Linear mixed modelling (LMM) was used to assess the impact of time from the first scan (for controls) and time from treatment initiation (for FEP), collectively referred to as time, in regional centiles. The contribution of age was separated between age of inclusion and time to ability, to capture specific interaction effects among treatment time, medication and centiles trajectory in patients:

(2) $$\small{\eqalign{\rm{regional\ centile} & \sim 1+ {\rm diagnosis+ gender + age\ of\ inclusion} \cr& \hskip 11pt + \, {\rm etiv} + {\rm medication+time+ time \times diagnosis} \cr & \,\,\,\,\,\,\,{\rm +\ medication \times time+(1\vert subject)}}}$$

where, in addition to the aforementioned variables, we also considered the following: medication at the time of MRI acquisition (antipsychotic dose in CPZ-equivalents), time × diagnosis (interaction between time and diagnosis) and medication × time (interaction indicating whether medication alters the effect of treatment time on centiles); and ‘1|subject as a random intercept to account for random subject-specific offsets overall. Age of inclusion, etiv, medication and time were z-scored prior to regression.

Sensitivity analyses repeated these models, using raw regional volumes rather than centiles as the dependent variable.

LMM for prediction of cognitive functioning as a function of treatment and regional centiles

After exploring the effects of disease and treatment on centiles, we evaluated how these factors relate to overall cognitive functioning:

(3) $$\small{\eqalign{\rm cognition & \sim 1+ \rm diagnosis+ \rm gender+ \rm age\ \rm of\ \rm inclusion + \rm etiv\cr & \hskip 11pt + \rm medication+\rm time+\rm centile+\rm time \times \rm diagnosis \cr & \hskip 11pt + \rm medication \times \rm time+centile \times diagnosis \cr & \hskip 11pt +\rm centile \times \rm time+(1\vert \rm subject)}}$$

Models were fitted independently for each region, resulting in one parameter estimate per region and covariate. For visualisation, a cortical map was generated for each covariate to summarise the regional patterns of association. All P-values were corrected for multiple comparisons across regions using FDR.

Generalised mixed modelling for prediction of clinical symptomatology as a function of treatment and regional centiles

Symptomatology (BPRS, SAPS and SANS scores) was modelled as:

(4) $$\small{\eqalign{& \rm symptoms \sim 1+ gender+age\ of\ inclusion + etiv+centile\cr& \hskip 52.51pt +\rm time + \rm medication+medication \times time \cr& \hskip 52.5pt +\rm centile\times medication +centile \times time+ (1\vert subject)}} $$

Generalised mixed models with a Poisson distribution were used, due to the count-like nature of symptom scores and their alignment with exponential decay patterns. Separate analyses were conducted for each region and psychometric scale.

Furthermore, for each scale, we conducted a mediation analysis to determine whether the association between medication and symptoms was mediated by the regional centile most strongly associated with symptoms. A bootstrap with 10 000 resamples was used to test for the presence of average direct effect (ADE), average causal mediation effect (ACME) and total effect (ADE + ACME).

Cytoarchitectural characterisation of centile associations

Cortical centiles associated with phenotypic variables (diagnosis, medication, BPRS, SAPS, SANS) were mapped onto the cortical hierarchies defined by the Mesulam atlas (unimodal, paralimbic, idiotypic, heteromodal; https://github.com/ucam-department-of-psychiatry/maps_and_parcs; Supplementary Table 5. ^{Reference Mesulam34} ANOVA tests and Tukey post hoc analyses were used to identify differential distributions of centile associations across Mesulam areas.

Neurobiological mapping of centile associations

Centile–phenotypic associations were also profiled using 46 molecular and microarchitectural maps derived from MRI/positron emission tomography studies, expanding on the methodology presented in ref. ^{Reference Hansen, Shafiei, Markello, Smart, Cox and Nørgaard35} and implemented at https://github.com/netneurolab/neuromaps and https://github.com/RafaelRomeroGarcia/NeurobiologyCentilesPsychosis. These maps included neurotransmitter receptors (5-HT_1A, 5-HT_1B, 5-HT_2A, 5-HT₄, 5-HTT, H₃, D₁, D₂, DAT, NET, α₄β₂, VAChT, CB₁, M₁, MOR, mGluR₅, NMDA, GABA); cell types (astrocytes, endothelial cells, microglia, oligodendrocytes, oligodendrocyte precursors, excitatory neurons, inhibitory neurons); cortical layers (I–VI); microstructural properties (synapse density, neurotransmitter PC1, gene expression PC1 (which is highly driven by cell types ^{Reference Dear, Wagstyl, Seidlitz, Markello, Arnatkevičiūtė and Anderson36} ); cortical thickness, myelin); metabolism (glycolytic index, cerebral blood flow, cerebral blood volume, oxygen metabolism, glucose metabolism); cortical expansion (evolutionary expansion, developmental expansion, allometric scaling from the Philadelphia Neurodevelopmental Cohort; and allometric scaling from the National Institutes of Health). See Supplementary Table 6 for the full list of neurobiological features included in the analysis.

Principal component–canonical correlation analysis (PC-CCA; https://github.com/RafaelRomeroGarcia/cca_pls_toolkit) was used to identify combinations of neurobiological features that best explained centile–phenotypic associations. ^{Reference Mihalik, Chapman, Adams, Winter, Ferreira and Shawe-Taylor37} PC-CCA reduces data dimensionality using PCA to minimise overfitting and then identifies the linear combination (weighted sum of loadings) of the neurobiological maps that best predicts the inter-regional variance of each centile–phenotypic association map. Statistical significance was evaluated using spatial autocorrelation-preserving permutation tests (spin tests; https://github.com/frantisekvasa/rotate_parcellation), with 10 000 parcellation-specific rotations. To validate the multiple regression patterns identified with PC-CCA, we conducted a sensitivity analysis with partial least squares (PLS). PLS maximises the covariance between the two variable sets, providing a robust alternative when collinearity might affect the canonical solution. ^{Reference Mihalik, Chapman, Adams, Winter, Ferreira and Shawe-Taylor37}

Limitations

Several limitations must be acknowledged. First, the primary inclusion criterion was the diagnosis of FEP at baseline, regardless of subsequent follow-up diagnoses, which predominantly included schizophrenia but also other psychotic disorders. Second, the use of CPZ-equivalents provides a simplified representation of the heterogeneous antipsychotic treatments across the sample, allowing us to capture the overall medication effect with a single parameter. However, different classes of antipsychotics may have differential effects on centiles and clinical outcomes. Inclusion of all drug-specific doses was not feasible, because this would have substantially increased the number of variables and the risk of overfitting. Finally, higher medication doses were generally prescribed to individuals with more resistant positive symptoms. However, these prescriptions were influenced by factors such as age of onset and treatment duration, complicating the disentanglement of individual contributions in the models.

This Brain-Wide Association Study (BWAS) demonstrates that centiles are associated with clinical outcomes in FEP. Nevertheless, growing evidence highlights important concerns: BWAS findings are susceptible to inflated effect sizes and low replication rates unless samples include thousands of individuals. ^{Reference Marek, Tervo-Clemmens, Calabro, Montez, Kay and Hatoum52} Recent findings by Kang et al have shown, however, that longitudinal designs can substantially improve replicability. ^{Reference Kang, Seidlitz, Bethlehem, Xiong, Jones and Mehta53} In that study, participants were comprehensively assessed over a 10-year period. While replication gains plateau after the first follow-up, particularly for structural markers, which exhibit limited within-subject variability, longitudinal assessments remain valuable for capturing the pronounced phenotypic decline in cognition and symptoms over time in FEP. Importantly, this high within-subject phenotypic variability is critical for enhancement of replicability, especially when modelled separately from between-subject variability, ^{Reference Kang, Seidlitz, Bethlehem, Xiong, Jones and Mehta53} as implemented here.

We aimed to maintain the models as parsimonious as possible. With 552 participants, up to 11 fixed variables and at least 449 degrees of freedom, the resulting balance was appropriate for robust and replicable estimation of model parameters. ^{Reference Maas and Hox54} Nevertheless, aspects such as the inclusion of moderate covariates, potential underlying collinearity or missing data further contribute to the complexity of assessing the actual statistical power and the potential risk of overfitting of the models.

Centiles derived from cross-sectional brain charts have demonstrated underestimation of longitudinal brain changes, resulting in inaccurate predictors of longitudinal individual changes. ^{Reference Di Biase, Tian, Bethlehem, Seidlitz, Alexander-Bloch and Yeo55} In the present study, which includes extensive follow-up data, our objective was not to perform individualised age-related prediction but to identify volumetric deviations associated with long-term cognitive and symptomatic outcomes in FEP. To account for residual age effects retained by the centiles, both age of inclusion and time were included as main effects in the models. However, longitudinal normative modelling was not performed due to the limited sample size. Future work will benefit from calibrating brain reference charts with longitudinal data-sets. ^{Reference Di Biase, Tian, Bethlehem, Seidlitz, Alexander-Bloch and Yeo55} This will be crucial for characterisation of individual-level trajectories, and for identification of subgroups that may account for clinical heterogeneity and differences in treatment response.

This study applied a normative modelling approach to demonstrate that FEP is characterised by regional decreases in cortical centiles. Centiles recovered to normative values during treatment, but were negatively affected by high doses of long-term medication. Similarly, improvements in cognitive abilities and symptomatology were positively associated with time but negatively influenced by prolonged medication use. Regions sensitive to symptomatology co-located with higher serotonin receptor densities and other neurobiological features. Collectively, these findings provide critical insights into the long-term interactions among cortical structure, treatment, medication, cognitive function and symptom profiles.