Regional, circuit, and network heterogeneity of brain abnormalities in psychiatric disorders

Anatomical Data

We estimated regional GMV using voxel-based morphometry⁸⁷ (VBM) of T1-weighted anatomical MRI scans^52,87. We focused on GMV because it is one of the most frequently studied neural phenotypes in psychiatry. VBM is also arguably the most widely used tool for measuring GMV in this context. Image acquisition parameters for each dataset are provided in Table S1. The following outlines our quality assurance procedures and data processing pipeline, which was applied to all raw T1-weighted images obtained for each dataset.

Quality Control

All T1-weighted images were visually inspected and evaluated for the presence of artifacts^88,89, resulting in the exclusion of 53 images with gross artifacts or abnormalities. Next, we used the Computational Anatomy Toolbox⁵² (CAT12 r113, http://dbm.neuro.uni-jena.de/cat/) to generate a weighted overall image quality rating (IQR) for every scan. This metric combines ratings of basic image properties, including the level of noise and geometric distortions, into a single score that quantifies the overall image quality of a participant’s T1-weighted scan (see here for more information www.neuro.uni-jena.de/cat/index.html#QA). On this metric, lower scores denote higher image quality. As per previous work²¹, we excluded 153 images with an IQR > 2.8. An additional 63 images were excluded due to a failure of our processing pipeline. Collectively, this quality assurance process resulted in the exclusion of a total of 269 scans.

MRI preprocessing

Regional GMV was estimated using the CAT12 VBM pipeline, which is included as an extension of Statistical Parametric Mapping software (SPM12, http://www.fil.ion.ucl.ac.uk/spm/software/spm12) in MATLAB v9.8⁹⁰. Briefly, the T1-weighted images were first corrected for intensity nonuniformities, and segmented into GM, white matter (WM), and cerebrospinal fluid (CSF) tissue probability maps. Then, using the high- dimensional Diffeomorphic Anatomical Registration Exponentiated Lie Algebra⁵³ (DARTEL), the segmented scans were normalized into standard Montreal Neurological Institute (MNI) space. Lastly, the images were bias-field corrected and modulated by the linear and nonlinear components of the Jacobian determinant obtained from the DARTEL deformation fields to obtain voxel-wise estimates of GMV. To constrain our analyses to GM voxels, we generated a mean image from all the normalized GM maps and retained voxels with a tissue probability ≥0.2.

Brain parcellation

We employed hierarchical Bayesian regression (HBR) to estimate normative models of regional GMV (see Normative Model section, below). To limit computational burden, we parcellated the brain into 1032 cortical and subcortical regions by combining well-validated parcellations of the cortex⁵⁴ comprising, 1000 regions, and of the subcortex⁵⁵, comprising 32 regions. The cortical parcels have been mapped to well-described, canonical functional networks of the brain⁵⁸, facilitating our analysis of network-level overlap in deviations. Voxels that overlapped between atlases were assigned to the corresponding subcortical region. Regional GMV estimates were obtained using freely available code (http://www0.cs.ucl.ac.uk/staff/gridgway/vbm/get_totals.m).

Note that the choice of parcellation will necessarily affect the magnitude of overlap that can be observed in any given brain area, such that coarser parcellations will be associated with higher levels of overlap. Comparison to the control data thus provides a critical normative benchmark against which to evaluate the levels of overlap observed in each clinical group. We chose our 1032-region parcellation to offer as much spatial resolution as possible while still ensuring that our analyses were computationally feasible. The same parcellation was used in regional, circuit, and network-level analyses, facilitating direct comparison across scales.

Normative Modelling

We used normative modeling to obtain person-specific GMV deviation maps in relation to an underlying model of normative expectations for regional GMV variations (see ²⁷ for details; code available at https://github.com/amarquand/PCNtoolkit, version=0.16). Normative models estimate the mean and variance (referred to as a normative range) of a response variable (e.g., GMV) from a set of clinically relevant covariates (e.g., age and sex) across a large healthy sample, referred to as a training set. These estimates are then used to quantify the deviations of samples in the test subset, which typically consists of cases, from the normative range.

When using multi-site data, scanner- and site-related variability introduce artefactual variance that confounds the result of any subsequent analyses⁹¹. These confounding effects present as site-correlated biases that cannot be explained by biological heterogeneity between samples. To account for these site-related effects, we built our normative model using hierarchical Bayesian regression (HBR), which successfully accommodates signal and noise variance in multisite data by estimating different but connected mean and variance components through shared prior distributions across sites⁵⁶. Prior to modelling, regional GMV estimates were first subjected to a Box-Cox transformation to ensure normality⁹². The optimal lambda parameter for minimizing skewness was estimated for each brain region independently using maximum likelihood.

For the model to accurately parse variance attributable to age, sex, and site, a sufficient number of observations is required for any given combination of these variables. We therefore excluded data cells containing less than 10 HCs for a particular sex at a given site. If this entire exclusion procedure resulted in less than 10 cases in total from any given site, the entire site was excluded, resulting in the exclusion of 217 scans.

The training data for the normative model (HC_train) was created by randomly selecting 90% of HC individuals from each site, provided that the site included data for 230 HCs. For sites with smaller HC samples, all HC data were included in the training set. The test data (HC_test), which was completely independent of the normative model, comprised the 10% of HCs from sites with 230 HCs and all case scans. The HC_test data, comprising 269 individuals (140 male), offered a normative benchmark for assessing case-specific model deviations, as outlined below.

Following stratification of our sample into training and test subsets, for each brain region, we fitted separate HBRs⁵⁶, modeling site and sex effects with random slopes, intercepts, and noise, to each of our 1032 parcellated brain regions to model GM volume as a function of age, sex, and site in the training data, yielding estimates of normative regional GMV variance and predictive uncertainty. Importantly, based on a partial-pooling approach, shared prior distributions were imposed over site- specific and sex-specific model parameters. These shared priors assume that while the model parameters for each site and sex are different, they are drawn from a common distribution. These shared priors thus regularize the model parameters and prevent the model from overfitting small batches for a given sex or site.

To quantify deviations from the normative model predictions for each subject in the test data, we generated deviation z-maps for each individual. Specifically, for each subject, i, at each brain region, j, we combined the predicted GM volume, , true GM volume, y_ij the predictive uncertainty, σ_ij and the normative variance, σ_nj, to calculate a z-score, z_ij as which quantifies the extent to which an individual’s regional GMV estimate deviates from the model prediction, given the uncertainty of the model.

For any given brain region, we are interested in participant values that deviate to a large extent from the normative expectations set by the model, under the assumption that these deviations represent pathophysiologically relevant features of the disease²⁷. We used two complementary approaches to characterize positive and negative deviations from the normative model at each region and for each individual. First, we thresholded the deviation maps at |z| > 2.6 (i.e., p < 0.005) to identify regions showing extreme deviations from model predictions. We used this approach to be consistent with prior work^21,93 and because alternative methods, such as the False Discovery Rate⁹⁴ (FDR), rely on adaptive thresholds that can vary between individuals. Applying a fixed z-threshold allows for an absolute definition of extreme deviations. Second, to ensure that our findings were not driven by this specific choice of threshold, we adapted a threshold-weighted approach⁵⁷, as detailed below.

Evaluation of model performance

To assess model generalizability, we used 5-fold cross-validation (CV) applied to the HC_train cohort (n = 1196 HC individuals, 55% male). Specifically, we split this group into 5 folds. Within each fold, we estimated an HBR model on 80% of participants, using age, sex, and site as covariates, to predict regional GMV in the remaining 20% of people. This procedure was repeated 5 times so that regional GMV values for all participants in the HC_train group were predicted once. As above, we estimated deviation z-maps for every individual and identified extreme deviations.

For each fold in the CV, we assessed model fit for each brain region by evaluating three performance metrics: (i) explained variance (EXPV); (ii) the mean standardized log-loss (MSLL); and (iii) the standardized mean-squared error (SMSE)⁵⁶, as shown in Figure S3. We also evaluated the model’s efficacy in partitioning site-related variance in the data using linear support vector machines (SVM). Specifically, we used a series of one-versus-all linear support vector machines (LSVM with default slack parameter = 1) trained separately on the z-maps from the HC_train and the HC_test subsets to classify scan sites. For each site, we ran a 2-fold SVM classifier to obtain the mean balanced accuracy score for the given site. Here, a balanced accuracy at chance-level (50%) indicated that the resulting deviations were not contaminated by residual site effects, as confirmed in Table S2. To assess whether individual variations in scan quality affected the deviation z-maps, we calculated the Pearson’s correlation between the total number of extrema and CAT12’s IQR rating in each group.

Characterizing the regional heterogeneity of extreme deviations

We characterized the heterogeneity of case-specific thresholded deviation maps for each disorder using a non-parametric approach. Specifically, for each clinical group and the HC_test cohort, we computed the proportion of individuals in each group showing an extreme deviation within each region, resulting in group-specific overlap maps estimated separately for positive and negative extrema (Figure 2f). We used the proportion of individuals rather than raw counts to account for sample size differences between groups. Next, we subtracted the HC_test overlap map from each disorder’s overlap map, resulting in an overlap difference map for each disorder (Figure 2g). We then permuted group labels (i.e.,HC, case) and repeated the procedure 10,000 times to derive an empirical distribution of overlap difference maps under the null hypothesis of random group assignment (Figure S6a-c). For each brain region, we obtained p-values as the proportion of null values that exceeded the observed difference (Figure S6d). The tails of the null distribution (i.e., values associated with p< 0.10) were approximated using a generalized Pareto distribution⁹⁵, as implemented in the Permutation Analysis of Linear Models software package⁹⁶ (PALM alpha116), to allow inference at arbitrarily high levels of precision. Statistically significant effects were identified using an FDR-corrected⁹⁴ threshold of p_FDR<.05, two-tailed.

Threshold-weighted deviation mapping

The procedure to characterize the regional heterogeneity of extreme deviations above used an arbitrary threshold (Z< |2.6|) on each individual z-map. The practical benefit of this approach is that it allows calculation of intuitive metrics, such as the proportion of cases within a group showing a supra-threshold deviation within a given brain region. We complemented this analysis with an alternative approach, which does not yield similarly intuitive metrics but which integrates results across a range of thresholds⁵⁷, thus allowing us to determine the degree to which our findings depend on a specific threshold choice. First, for each diagnostic group, we thresholded z-maps across the threshold range 1.64 < Z < 3.10, in 100 equal log-space increments, separately for positive and negative extrema. Second, we obtained a region- and threshold-specific percentage overlap map quantifying the proportion of individuals within each group showing an extreme deviation for any given region and threshold. Third, we applied a weighted function to penalize less conservative thresholds. Fourth, we obtained a final threshold-weighted overlap map by taking the area under the curve of the cumulative histogram of the proportion of participants at each threshold for each region across participants. We used the following linear weighted function proposed by Seghier and Price⁵⁷, which ensures that overlap map values range between 0 and 1, with 1 indicating an effect in a region is present in each subject at each threshold: where Z_max and Z_min respectively correspond to the maximum and minimum z-value thresholds, and th is the given statistical threshold. We then calculated the difference between the threshold- weighted percentage overlaps between cases and the held-out HC cohort, resulting in an overlap difference map for each disorder. We performed inference on these differences using the same group- based permutation procedure described above with statistically significance defined using a threshold of p_FDR<.05, two tailed.

Functional magnetic resonance imaging (fMRI) data and processing

The functional circuitry of deviant regions was mapped using an independent cohort of 150 resting-state functional MRI (rs-fMRI) scans (71 males, age 21- 35 years) from the S900 release of the Human Connectome Project⁹⁷ (HCP). These individuals corresponded to the 150 people with the lowest total head motion, as estimated using framewise displacement⁹⁸, from the broader set of 282 unrelated participants with the same fMRI reconstruction in the S900 release. We focused on this sub- sample to minimize the effects of head motion in our FC maps and to minimize computational burden due to the large number of analyses required. Mapping deviation-related functional circuity in an independent cohort is essential because we seek to understand the circuit context of deviant regions under normative conditions. Our approach aligns precisely with the logic used in lesion network mapping of neurological patients^31,32.

The rs-fMRI data for each subject consists of four runs (rfMRI_REST1_LR, rfMRI_REST1_RL, rfMRI_REST2_LR, rfMRI_REST2_RL) acquired at two different sessions (REST1 and REST2) using two different directions of phase coding (LR: left to right, and RL: right to left). Here, we mapped deviation-related FC using the rfMRI_REST1_LR data. All rs-fMRI data were obtained using a 32-channel Siemens 3T connectome-Skyra scanner. The imaging parameters for rs- fMRI were as follows: repetition time (TR) = 720 ms, echo time (TE) = 33.1 ms, flip angle (FA) = 52°, field of view (FOV) = 208 × 180 mm², matrix = 104 × 90, slice number = 72, slice thickness = 2 mm, voxel size = 2 × 2 × 2 mm³, multiband factor = 8, and 1200 volumes. Participants were required to stay awake, relaxed, and to keep their eyes open and fixed on a bright cross-hair projected on a dark background. A detailed description of the HCP data is available in previous work⁹⁷.

The rs-fMRI data used underwent the HCP’s minimal pre-processing pipeline⁹⁹, which includes gradient-nonlinearity-induced distortion, motion correction to the single-band reference image using FLIRT, EPI image distortion correction using TOPUP¹⁰⁰, registration into standard space using a customized boundary-based-registration (BBR) algorithm, and single step spline interpolation using all transforms, intensity normalization and bias field removal to resample the original EPI into MNI space. Minimal high pass filtering was applied with a cutoff of 2000ms. Artifacts were then removed using ICA-FIX¹⁰¹. This involves employing an automatic classifier, specifically trained for HCP data, to identify ICA components due to measurement noise, additional motion or physiological artifacts like cardiac pulsation and respiration. Next, the volume timeseries were mapped into the standard CIFTI grayordinate space⁹⁷ and smoothed to 2mm FWHM (where the smoothing was on the surface for the cortex and in volume space for subcortex). This results in a standard set of grayordinates in every subject, with surface vertex data and subcortical volume voxel data. The mean GM signal was then removed each grayordinate’s time series to remove residual widespread signal deflections that are not removed by ICA-FIX^102–104

Mapping the functional circuitry of extreme deviations

We used seed-based FC analysis to map the functional circuitry of each region showing an extreme deviation in any individual within the clinical and HC_test groups. Specifically, each deviant region was used as a seed, from which the average time course was extracted for each individual in the HCP₁₅₀ sample (Figure 4b). These seed time courses were then correlated with all other brain grayordinates and the resulting correlation maps were subjected to Fisher’s r-to-z transformation¹⁰⁵. The transformed correlation maps for each HCP₁₅₀ subject were then aggregated using a one-way t- test at each grayordinate, as implemented in PALM. Thresholded maps representing the functional circuitry of the seed were obtained using permutation testing (500 permutations) with a generalized Pareto approximation of the tail of the null distribution⁹⁵ and threshold-free cluster enhancement¹⁰⁶ (TFCE), run separately for the cortical and sub-cortical areas, with a threshold of p_FWE<.025, family- wise error corrected (i.e., .05/2 to account for cortical and subcortical analyses). This procedure was repeated for each region in which at least one extreme deviation across individuals was observed (Figure 3b).

To quantify circuit-level overlap, we developed an approach that parallels traditional lesion network mapping methodology but which allows statistical inference on observed case-control differences in overlap scores at each region. Specifically, we thresholded each deviant-related FC map by applying TFCE to the grayordinate maps, as described above. We then classified an area in our 1032-region parcellation as showing significant FC with the seed if more than 50% of its grayordinates survived the TFCE threshold. We repeated our analysis using a mapping threshold of 75% to ensure that our findings were not driven by this specific choice (see Figure S13-15). We mapped FC at the grayordinate level and thresholded the data in this way to leverage the superior statistical sensitivity of TFCE and to map the spatial architecture of the seed-related FC patterns more accurately. Our procedure resulted in a binary map representing the specific brain regions that comprise the functional circuitry of each seed (Figure 3c). Then, for each individual, we took the union of the binary FC maps across that person’s set of extrema, resulting in a single map identifying areas showing significant FC to any extreme deviation expressed by the person (Figure 3d). Separate union maps were obtained for positive and negative deviations. The overlap of these union maps was then taken for each clinical group and the HC_test sample, representing, for each region, the proportion of individuals in that group for which significant deviant-related FC was identified. We refer to this image as a circuit-level overlap map (Figure 3e).

Note that our approach focuses on FC between deviant and other areas as derived from a normative sample. It is possible that patterns of FC are altered in diagnosed individuals, and that the underlying network architecture differs from normative expectations. The net effect of such alterations would be to affect coupling between deviant and other areas, and thus our normative benchmark still offers an important reference point for understanding the distributed effects of regional deviations. Further investigation of how altered patterns of FC may moderate these relationships remains an important topic for future investigation.

Statistical inference on group differences in circuit-level overlap

We evaluated group differences in the circuit-level overlap maps of each clinical group relative to the HC_test group in two ways (see Figure S6 for a schematic overview). First, following the analysis of regional overlaps, we computed the difference in the case and control overlap maps (Δ overlap map, Figure 3e) and evaluated the statistical significance of these differences with respect to an empirical null distribution generated by permuting the group labels of the individual-specific FC union maps 10,000 times, with a generalized Pareto Tail approximation⁹⁵. Statistically significant differences in deviation FC overlap were identifying using a threshold of p_FDR < 0.05, two-tailed (Figure S6a-d).

This group-based permutation test identifies regions where the overlap in deviant-related functional circuitry differs between cases and controls, but it does not preserve group differences in the total number of deviations. It therefore cannot distinguish whether any observed differences arise from a preferential aggregation of deviations within specific circuits in one group relative to the other, or whether the group differences are driven by variations in total deviation burden. Note that differences in deviation burden should not be thought of as a confound in this context, since they are an intrinsic feature of the disease that will have real phenotypic consequences. Nonetheless, it is important to determine which findings might be driven by differences in deviation burden compared to a preferential accumulation of pathology within a specific circuit for a given disorder.

To disentangle these possibilities, we considered a second, spatial null model that preserves group differences in deviation burden. Specifically, we generated an ensemble of null cortical deviation maps for each individual in the test data by spatially rotating their empirical, unthresholded deviation maps (Figure S6e) using Hungarian spherical spin tests^107,108 (Figure S6f) and then thresholding these rotated maps in the same way as the observed data (i.e., Z> |2.6|; Figure S6g). The Hungarian method was used as it preserves the original values and spatial autocorrelation of the original map, which ensures that the surrogate data yields the same number of extrema as originally observed. This procedure could only be applied to the cortex, which is topologically equivalent to a sphere. The choice of an appropriate null model for subcortical areas is more complicated. While model-based procedures (e.g.,¹⁰⁹) can generate surrogates with comparable spatial autocorrelation, they require parameter tuning, can show variable fits across individuals, and do not preserve the exact values of the original data, making thresholding challenging. We therefore randomly shuffled deviation values across subcortical areas before thresholding them. Note that this null model is more lenient than a spatially constrained test since the model preserves fewer features of the data. As we found no significant differences using this test, our results will not change under a more stringent, spatially constrained null model.

We then obtained individual-specific surrogate FC union maps using the same procedure described above (Figure S6h), and calculated surrogate within-group overlap maps (Figure S6i) and between-group overlap difference maps (Figure S6j). This procedure was repeated 10,000 times to generate a null distribution of overlap difference maps for each disorder, again using the generalized Pareto tail approximation (Figure S6k). Statistically significant differences were identified using a threshold of p_FDR < 0.05, two-tailed.

To summarize, group differences in circuit-level overlap were assessed with respect to two null models, one based on permutation of group labels and one based on spatial permutation of individual deviation maps. The group-based permutation test identifies regions showing differences in circuit-level overlap regardless of deviation burden. The spatial permutation test can be used to assess the degree to which such differences are driven by variations in deviation burden. More specifically, if a given brain region shows a significant difference for both tests, we have evidence to indicate that extreme deviations within the disorder preferentially aggregate within the functional circuitry of that region. If, on the other hand, a region shows a difference with the group-based permutation test but not the spatial permutation test, then the observed group differences are consistent with expectations from randomly selecting of the same number of seeds in each group, indicating that the effects observed under group-based permutation can be attributed to the fact that one group has a higher deviation burden than the other. Note that this scenario still yields important information about the extent of circuit-level heterogeneity in the disorder, since we generally expect cases to have a higher deviation burden and any resulting circuit-level overlap will still have phenotypic consequences regardless of whether such overlaps are explained by deviation burden or not. The spatial permutation test simply offers additional insights into the potential mechanisms that may drive group differences in circuit-level overlap.

Assigning regions to networks

We assigned each cortical region to one of 7 canonical functional cortical networks using a well-validated network parcellation^54,58 and assigned each subcortical region to either the medial temporal lobe (amygdala and hippocampus), thalamus, or basal ganglia (nucleus accumbens, globus pallidus, putamen, caudate nucleus), as done previously⁵⁵, resulting in a total of 10 distinct functional networks (Figure 4a-c). The subcortical regions were assigned to coarse anatomical areas that were not aligned with the finely mapped cortical function networks because the correspondence between these subcortical nuclei, as parcellated here, and canonical cortical functional networks has not been extensively investigated. An alternative network assignment may influence network-level subcortical overlap values.

Characterizing network-level heterogeneity

To characterize network-level heterogeneity, we examined the degree to which GMV deviations aggregated within each network. To this end, for each network, we estimated the proportion of individuals within each diagnostic group that showed at least one extreme deviation in a region assigned to that network, separately for positive and negative deviations (Figure 4d). We then computed clinical-control differences in overlap proportions for each network (Figure 4e).

Statistical inference on group differences in network-level overlap

As in the circuit-level analysis, these observed group differences in network-level overlap were evaluated with respect to two complementary null models. The first involved a group-basedrandomization test in which a null distribution was generated by permuting group labels 10,000 times (Figure S6a-d). The second null model used the 10,000 spatially rotated maps generated in the circuit- level analysis to obtain surrogate estimates of group differences in network-level overlap while preserving variation in deviation burden (Figure S6e-g,l-o). In both cases, the tail of the null distribution was approximated using a generalized Pareto distribution and we used a significance threshold of p_FDR < 0.05, two-tailed. As in the circuit-level analysis, these two null models allowed us to distinguish preferential network involvement from the effects of total deviation burden. Note that comparisons in overlap values between networks are complicated by differences in network size. Our analysis therefore focuses primarily on comparisons between cases and controls within networks since the control data provide a critical normative benchmark against which cases can be compared.