# A comparative machine learning study of schizophrenia biomarkers derived from functional connectivity

Victoria Shevchenko, R Austin Benn, Robert Scholz, Wei Wei, Carla Pallavicini, Ulysse Klatzmann, Francesco Alberti, Theodore D Satterthwaite, Demian Wassermann, Pierre-Louis Bazin, et al.

► To cite this version:

Victoria Shevchenko, R Austin Benn, Robert Scholz, Wei Wei, Carla Pallavicini, et al. A comparative machine learning study of schizophrenia biomarkers derived from functional connectivity. Scientific Reports, 2025, 15 (2849), 10.1038/s41598-024-84152-2. hal-04911960

HAL Id: hal-04911960

https://hal.science/hal-04911960v1

Submitted on 25 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

OPEN

# A comparative machine learning study of schizophrenia biomarkers derived from functional connectivity

Victoria Shevchenko<fnref n="1" />,²,³,⁴✉, R. Austin Benn<fnref n="1" />,², Robert Scholz<fnref n="1" />,²,⁵,⁶ Wei Wei<fnref n="1" />,², Carla Pallavicini<fnref n="1" />,⁷,⁸ Ulysse Klatzmann<fnref n="1" />,², Francesco Alberti<fnref n="1" />,², Theodore D. Satterthwaite⁹, Demian Wassermann³,⁴, Pierre-Louis Bazin¹⁰ & Daniel S. Margulies<fnref n="1" />,²✉

Functional connectivity holds promise as a biomarker of schizophrenia. Yet, the high dimensionality of predictive models trained on functional connectomes, combined with small sample sizes in clinical research, increases the risk of overfitting. Recently, low-dimensional representations of the connectome such as macroscale cortical gradients and gradient dispersion have been proposed, with studies noting consistent gradient and dispersion differences in psychiatric conditions. However, it is unknown which of these derived measures has the highest predictive capacity and how they compare to raw functional connectivity specifically in the case of schizophrenia. Our study evaluates which connectome features derived from resting state functional MRI — functional connectivity, gradients, or gradient dispersion — best identify schizophrenia. To this end, we leveraged data of 936 individuals from three large open-access datasets: COBRE, LA5c, and SRPBS-1600. We developed a pipeline which allows us to aggregate over a million different features and assess their predictive potential in a single, computationally efficient experiment. We selected top 1% of features with the largest permutation feature importance and trained 13 classifiers on them using 10-fold cross-validation. Our findings indicate that functional connectivity outperforms its low-dimensional derivatives such as cortical gradients and gradient dispersion in identifying schizophrenia (Mann-Whitney test conducted on test accuracy: connectivity vs. 1st gradient: $U = 142, p < 0.003$; connectivity vs. neighborhood dispersion: $U = 141, p = 0.004$). Additionally, we demonstrated that the edges which contribute the most to classification performance are the ones connecting primary sensory regions. Functional connectivity within the primary sensory regions showed the highest discrimination capabilities between subjects with schizophrenia and neurotypical controls. These findings along with the feature selection pipeline proposed here will facilitate future inquiries into the prediction of schizophrenia subtypes and transdiagnostic phenomena.

Functional connectivity holds promise as a potential biomarker for schizophrenia<fnref n="1" />⁻⁵, as evidenced by a robust body of fMRI literature that highlights distinct functional profiles between people with schizophrenia and neurotypical individuals. Prior studies have reported lower connectivity across regions, reduced small-worldness of the resting state networks, and lower functional network segregation⁵⁻⁷. However, when the goal is to predict clinical status using the entire functional connectome as features for model training, the resulting model becomes high-dimensional and overly complex. This high dimensionality, coupled with small sample sizes in clinical research, increases the risk of overfitting⁸⁻¹⁰.

Recently, low-dimensional representations of the connectome such as macroscale cortical gradients¹¹,¹² and gradient dispersion¹³,¹⁴ have been proposed. The gradients are derived from functional connectivity matrices

<footnote n="1">¹Cognitive Neuroanatomy Lab, INCC UMR 8002, CNRS, Université Paris Cité, Paris, France. ²Wellcome Centre for Integrative Neuroimaging, Nuffield Department of Clinical Neurosciences, FMRIB Centre, University of Oxford, Oxford, UK. ³MIND Team, Inria Saclay, Université Paris-Saclay, Palaiseau, France. ⁴Neurospin, CEA, Gif-Sur-Yvette, France. ⁵Max Planck School of Cognition, Leipzig, Germany. ⁶Wilhelm Wundt Institute for Psychology, Leipzig University, Leipzig, Germany. ⁷National Scientific and Technical Research Council (CONICET), Buenos Aires, Argentina. ⁸Department of Physics, Institute of Applied and Interdisciplinary Physics, University of Buenos Aires, Buenos Aires, Argentina. ⁹University of Pennsylvania, Perelman School of Medicine, Philadelphia, PA, USA. ¹⁰Full Brain Picture Analytics, Leiden, The Netherlands. ✉email: victoria.shevchenko@inria.fr; daniel.margulies@cnrs.fr</footnote>
through dimensionality reduction algorithms such as principal component analysis (PCA) or diffusion map embedding. The aim of this computation is to maximize the cumulative amount of variance explained by the resulting components. The first component (also known as the *principal gradient*) explains the largest fraction of connectivity variance. It reflects the functional hierarchy of the cortex<sup>11,12,15</sup>, spanning from the primary sensory (*unimodal*) regions to higher-order (*transmodal*) regions. The principal gradient has been demonstrated to be consistent across individuals<sup>16,17</sup>.

Dong et al.<sup>18</sup> revealed that the principal gradient is contracted in schizophrenia. That is, the primary sensory regions were reported to be closer to the higher-order regions in terms of their functional connectivity profile, as indicated by their principal gradient values. This finding indicates lower functional differentiation between uni- and transmodal regions. In other words, the primary sensory and higher-order processing regions are more functionally similar in subjects with schizophrenia compared to neurotypical individuals. In addition, Holmes et al.<sup>19</sup> reported significant differences in the second gradient which spans the hierarchy of primary sensory areas. Based on the gradient framework, other neurodevelopmental, psychiatric and neurodegenerative disorders have been reported to manifest changes in the cortical functional hierarchy<sup>14,16,20-28</sup>.

As an extension of the gradient framework, gradient dispersion further quantifies the density of local connectivity between or within clusters of areas. This density characterizes the extent to which the areas are functionally segregated or integrated. Specifically, higher dispersion would indicate higher functional segregation within a cluster of regions, whereas lower dispersion would mean higher functional integration<sup>13,27</sup>. Thus, gradient dispersion characterizes functional differentiation across the cortical hierarchy: high dispersion of a cluster of regions points to their functional dissimilarity, placing them farther apart along the functional hierarchy of the cortex. Gradient dispersion is of interest to us since changes in the cortical hierarchy have been reported to be idiosyncratic to schizophrenia<sup>1,18,29</sup>. However, gradient dispersion has not yet been tested as a potential predictor of schizophrenia.

Overall, there is a growing interest in the macroscale cortical hierarchy as a predictor of various phenotypes, be it age<sup>13,30</sup>, psychosis<sup>18,19</sup>, or task performance<sup>31,32</sup>. However, it is unknown if cortical gradients or gradient dispersion can discriminate between people with schizophrenia and neurotypical individuals, or if they outperform raw connectivity. Given this gap and the recent prominence of the gradient framework, we selected gradients and gradient dispersion as biomarker candidates for this study, in addition to raw functional connectivity. Here, we attempt to identify the features with the largest biomarker potential from a large set of features including the aforementioned measures. In addition, we explore the impact of the number of features on the choice of classifier. We also seek to address the question germane to neuroscience and computational psychiatry: when one has a limited number of subjects and a disproportionately rich set of independent variables, how does one justify the choice of features? To this end, we develop a pipeline based on permutation feature importance which allows us to assess the predictive power of all features at once in a computationally efficient manner (Fig. 1). Additionally, we elaborate on putative functional underpinnings of schizophrenia based on the features with the highest predictive potential.

# Methods

## Data and preprocessing

The present study's sample was derived from three publicly available datasets: COBRE<sup>33</sup>, LA5c study from UCLA Consortium for Neuropsychiatric Phenomics<sup>34</sup>, and the SRPBS-1600 multidisorder MRI dataset<sup>35</sup>. All data were acquired in accordance with the Declaration of Helsinki. For each dataset, the acquisition was approved by the local institutional review board (IRB). Informed consent was obtained from all subjects and/or their legal guardian(s) to the collection of imaging and behavioral data. The links to the datasets are provided in the data availability statement below.

The initial cohort for our investigation consisted of 996 individuals; subsequently, 60 subjects were excluded from analysis due to substantial motion artifacts (mean framewise displacement (FD) > 0.5 mm). The analyzed sample comprised 248 subjects with schizophrenia (SCZ) and 688 neurotypical controls (NC). We conducted a Chi-squared test for sex and Mann-Whitney U tests for age and mean FD. Some datasets were significantly different in terms of sex, age and mean FD distributions (see Table 1). We address this limitation by including these variables as covariates in all analyses. The scanning parameters of each dataset are reported in Supplementary Table 1.

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Overview of the methods and main outcome of the paper. Schematic images: Flaticon.com. NC: neurotypical controls, SCZ: individuals with schizophrenia.</figcaption></figure>
<table><caption>**Table 1.** Demographic statistics of the sample. Significant comparisons are in bold and underlined. *These datasets are part of a larger dataset, SPRBS-1600. FD: framewise displacement.</caption><thead><tr><th>Site</th><th>Group</th><th>N</th><th>Mean age (SD)</th><th>U<sub>AGE</sub></th><th>N female (N male)</th><th>χ<sup>2</sup><sub>SEX</sub></th><th>Mean FD (SD) (mm)</th><th>U<sub>FD</sub></th></tr></thead><tbody><tr><td rowspan="2">COBRE</td><td>NC</td><td>81</td><td>37.9 (11.9)</td><td rowspan="2">2598</td><td>23 (58)</td><td rowspan="2">1.87</td><td>0.25 (0.1)</td><td rowspan="2"><strong>1860</strong></td></tr><tr><td>SCZ</td><td>59</td><td>36.4 (13.2)</td><td>10 (49)</td><td>0.29 (0.1)</td></tr><tr><td rowspan="2">LA5c</td><td>NC</td><td>102</td><td>30.6 (8.3)</td><td rowspan="2"><strong>1456</strong></td><td>49 (53)</td><td rowspan="2"><strong>7.62</strong></td><td>0.21 (0.07)</td><td rowspan="2"><strong>1232</strong></td></tr><tr><td>SCZ</td><td>45</td><td>36.3 (8.9)</td><td>10 (35)</td><td>0.21 (0.09)</td></tr><tr><td rowspan="2">KTT*</td><td>NC</td><td>75</td><td>28.9 (9)</td><td rowspan="2"><strong>879</strong></td><td>27 (48)</td><td rowspan="2">0.39</td><td>0.1 (0.04)</td><td rowspan="2"><strong>1260</strong></td></tr><tr><td>SCZ</td><td>46</td><td>37.6 (9.7)</td><td>20 (26)</td><td>0.13 (0.05)</td></tr><tr><td rowspan="2">KUT*</td><td>NC</td><td>159</td><td>36.5 (13.6)</td><td rowspan="2"><strong>2497</strong></td><td>66 (93)</td><td rowspan="2">1.51</td><td>0.15 (0.07)</td><td rowspan="2">3112</td></tr><tr><td>SCZ</td><td>43</td><td>42.1 (10.4)</td><td>23 (20)</td><td>0.15 (0.07)</td></tr><tr><td rowspan="2">SWA*</td><td>NC</td><td>101</td><td>28.4 (7.9)</td><td rowspan="2"><strong>205</strong></td><td>15 (86)</td><td rowspan="2">0.11</td><td>0.15 (0.07)</td><td rowspan="2"><strong>472</strong></td></tr><tr><td>SCZ</td><td>19</td><td>42.9 (8.4)</td><td>4 (15)</td><td>0.22 (0.09)</td></tr><tr><td rowspan="2">UTO*</td><td>NC</td><td>170</td><td>35.6 (17.5)</td><td rowspan="2">3194</td><td>92 (78)</td><td rowspan="2"><strong>4.34</strong></td><td>0.12 (0.07)</td><td rowspan="2"><strong>2810</strong></td></tr><tr><td>SCZ</td><td>36</td><td>31.4 (10.3)</td><td>12 (24)</td><td>0.13 (0.06)</td></tr><tr><td rowspan="2">All data</td><td>NC</td><td>688</td><td>33.5 (13.2)</td><td rowspan="2"><strong>65103</strong></td><td>272 (416)</td><td rowspan="2"><strong>4.34</strong></td><td>0.15 (0.08)</td><td rowspan="2"><strong>62110</strong></td></tr><tr><td>SCZ</td><td>248</td><td>37.4 (11.1)</td><td>79 (169)</td><td>0.19 (0.1)</td></tr></tbody></table>

Preprocessing of MRI data was done using fMRIPrep 20.2.1³⁶ which is based on Nipype 1.5.1³⁷ (Supplementary Methods 1). The preprocessed BOLD time series were parcellated with the Schaefer parcellation (1000 parcels, 7 Yeo networks)³⁸. Then, we computed a connectivity matrix (Pearson correlation) for each subject.

## Macroscale cortical gradients

For each subject, we computed cortical gradients by applying PCA to the Fisher z-transformed and thresholded connectivity matrix (Fig. 2A). Hong et al.¹⁶ showed that PCA, when applied to thresholded connectivity matrices, yields more reliable gradients compared to the other dimensionality reduction techniques frequently featured in the gradient literature¹²,³⁹,⁴⁰. We thresholded the connectivity matrices by discarding 90% of the lowest correlation values including negative values. We used Procrustes alignment to align the gradients of all subjects⁴¹. To avoid introducing dataset-specific bias to the alignment, we used the gradients computed from the group connectivity matrix of the Human Connectome Project (HCP)⁴² as reference gradients. Figure 2B displays mean variance explained across all subjects for 200 gradients. On average, the principal gradient accounted for ~6% of variance of thresholded connectivity matrices. Collectively, 200 gradients accounted for ~80% of variance. Given that 1000 gradients cumulatively explained all variance (Supplementary Fig. 1), we deemed 200 gradients computationally optimal as 20% of gradient values in this case account for 80% of variance.

## Centroid (gradient) dispersion

Gradient dispersion was first introduced by Bethlehem et al.¹³ and has since been employed in other studies¹⁴,¹⁹,⁴³. Prior investigations, akin to the present study, used gradient dispersion to operationalize functional modularity— i.e., functional (dis)similarity—across the cortex. However, these methods necessitated the identification of centroids of networks for which dispersion was computed relative to the regions encompassed within the network¹³. Thus, dispersion quantifies functional modularity within a network or between networks.

Following this procedure, we computed within- and between-network dispersion for seven Yeo networks⁶⁵ in the 3-dimensional gradient space constituted by the first three principal components of functional connectivity. Within-network dispersion was quantified as the sum of squares of Euclidean distances from the centroid of the network to the rest of its regions. For the gradient values belonging to a given network, the centroid was defined as three median values of the first three cortical gradients, as in Bethlehem et al.¹³; the position of the centroid in the 3D latent gradient space is defined by these three values. Thus, given a network *k* consisting of *P* parcels with each parcel *p* having gradient values *g*, within-network centroid dispersion is computed as follows:

$$ \text{Within - network CD}_k = \sum_{p=1}^{P} \left\| g_{p,k} - \text{median}(G_k) \right\|^2 \quad (1) $$

where $G_k$ is a matrix of all gradient values of all parcels in network *k*, $g_{p,k} \in \mathbb{R}^3$.

Between-network dispersion was defined as the Euclidean distance between the network centroids as follows:

$$ \text{Between - network CD}_{kl} = \|\text{median}(G_k) - \text{median}(G_l)\| \quad (2) $$

where $G_k$ and $G_l$ are matrices containing all gradient values of all parcels in networks *k* and *l*.
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. (A) Parcel-wise time series (Schaefer atlas, 1000 parcels, 7 Yeo networks³⁸) of each subject were correlated to produce a 1000 × 1000 connectivity matrix. Principal component analysis (PCA) was applied to the thresholded matrix to extract 200 gradients. (B) Variance explained by 200 gradients, mean across subjects ± 1 s.d.</figcaption></figure>

As a result of these computations, centroid dispersion amounted to 28 values per subject: 7 values for within-network and 21 values for between-network dispersion. Centroid dispersion is schematically illustrated in Fig. 3A.

## Neighborhood (gradient) dispersion

Centroid dispersion can vary depending on how the networks are delineated (e.g., if a different network parcellation is used). Hence, we computed neighborhood dispersion with the aim to circumvent this potential confounding factor. Specifically, we calculated dispersion for individual regions which enabled us to maintain the spatial resolution congruent with the gradients (1000 regions per measure).

To compute neighborhood dispersion for every region, we identified K closest neighboring regions (Fig. 3B,C) via the K-Nearest Neighbors (KNN) algorithm. Then, we computed the mean Euclidean distance between the focal region and its designated neighboring regions in the gradient space. The resulting value quantified dispersion for the focal region. Neighborhood dispersion was computed for combinations of gradients spanning from 1 to 200.

Unlike previous investigations where the primary source of variability in dispersion originated from network delineation, our study has its own unique challenge in determining the number of the nearest neighbors (K). To address this issue, we included all dispersion values computed based on a range of nearest neighbors from 10 to 170 with a step size of 40. We chose this step size to reduce the computational load of this operation, i.e., we did not perform any optimization for K. Mathematically, neighborhood dispersion can be formulated as follows:

$$ND_f(K) = \frac{1}{K} \sum_{j=1}^{K} \|x_f - x_j\| \quad (3)$$

where $K \in [10, 50, 110, 140, 170]$ (number of neighbors), $x_f, x_j \in \mathbb{R}^d$ with $x_f$ being the gradient values of the focal region and $x_j$ being the gradient values of its $j^{th}$ neighbor in a $d$-dimensional gradient space ($d \in [1: 200]$).
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Illustration of the methods used to compute gradient dispersion. (A) Centroid dispersion. The sum of squares of distances between the centroid of a network and its regions (black dashed lines) quantifies within-network dispersion. The distance between the centroids of networks quantifies between-network dispersion. (B) Neighborhood dispersion. In a multidimensional gradient embedding, for a given region (red) K nearest neighbors are identified (blue). These regions are shown within the black circle. (C) Neighborhood dispersion of a given region i is the mean distance between said region and its K closest neighbors. The same operation is done for every region (N regions = 1000).</figcaption></figure>

Basing the calculation of dispersion on differing sets of gradients (each including up to 200 gradients), we derived 1000 × 200 × 5 = 1,000,000 neighborhood dispersion values, which we then used as input to our analytic workflow (along with flattened connectivity matrices and gradients).

In summary, centroid dispersion is computed based on the Yeo networks and the first three gradients as in Bethlehem et al.¹³ (Eqs. 1, 2). Conversely, neighborhood dispersion is computed using a range of the nearest neighbors and combinations of gradients from 1 to 200 (Eq. 3).

## Analytic workflow

The objective of our workflow was to identify the features with the largest predictive capacity from an extensive array of connectivity-based features. Our dataset included vectorized connectivity matrices (N=499,500), 200 gradients (N=200,000; 1000 values per gradient; 200 values per region), neighborhood dispersion (N=1,000,000; 1000 values per region) and centroid dispersion (N=28) (Fig. 4A). Since the number of candidate features is too large for every feature combination to be tested separately (for a fixed $N_{features}$ = 1000: $N_{combs}$ = $\frac{1699528!}{1000!(1699528-1000)!} \approx 3.95 * 10^{3010}$), we designed a custom feature selection pipeline (Fig. 4). For each participant, we combined all features into one flat vector of size [1,699,528]. Then, we concatenated the vectors for all 936 participants, resulting in a matrix of size [936 × 1,699,528] (Fig. 4B, 1). Next, we applied PCA to each type of features separately and retained 20% of variance for each type (effectively compressing the feature dimension), except for the centroid dispersion for which all variance was included (Fig. 4B, 2). The aim of the decomposition was i. to retain the same amount of variance for all feature types, ii. to ensure that for each feature type more than one component is extracted when applying PCA, and iii. to alleviate the imbalance of the number of features available for each feature type (feature type ratio in the full feature matrix: 50:20:100:0.0028 and in the PCA feature matrix: 7:72:42:28). Thus, we retained 149 components in total ($N_{conn}$=7, $N_{grad}$=72, $N_{centroid Disp}$=28, $N_{cortex Disp}$=42; Fig. 4B, 2). The decomposed dataset was divided into the train (75% of participants) and holdout (25%) sets, with the ratio of $N_{SCZ}/N_{NC}$ = 0.35 in both.

The purpose of the following steps was to quantify the importance of each component for classification performance. We fitted an L2-regularized logistic regression on the train set (Fig. 4B, 3). Logistic regression was determined as the best model to compute component importance since the resulting coefficients are interpretable and the logic behind their computation is well understood. We used permutation feature importance (Fig. 4C, 1) as the measure of the contribution of each component and feature to classification performance. Initially conceived for random forests⁴⁴,⁴⁵, it allows us to estimate the importance of the features in a classifier-agnostic way using the holdout set. First, the classifier is trained on the train set and the baseline accuracy on the holdout set is obtained. Second, each feature of the holdout set is randomly shuffled $N_{perm}$ times and at each shuffle the permutation accuracy is computed. Permutation feature importance is the difference between the baseline accuracy and the permutation accuracy at each permutation. This feature importance metric can be estimated for any classifier; it quantifies the extent to which classification performance deteriorates or improves as every feature is shuffled in the holdout set. The feature with the largest permutation importance contributes the most to classification performance (Fig. 4C, 1). We computed mean permutation importance for each component ($N_{perm}$ = 10,000) and inverse transformed it based on the components' projection matrices to obtain feature importance in the feature space (Fig. 4C, 2). Permutation feature importance enabled the selection of features for further assessment of their utility for classification.
A All features

<figure><img src="image_7.png" /></figure>

B Modeling Pipeline

<figure><img src="image_8.png" /></figure>

C Permutation Feature Importance

<figure><img src="image_9.png" /></figure>

$$
\begin{align*}
I_{\text{PCA}} \leftarrow & I_{\text{conn}} = P_{\text{conn}}^{-1} I_{\text{PCA}_{\text{conn}}} \\
& \cdot \\
& \cdot \\
& I_{\text{centroid\_disp}} = P_{\text{centroid\_disp}}^{-1} I_{\text{PCA}_{\text{centroid\_disp}}}
\end{align*}
\quad \rightarrow \quad
\begin{tabular}{c}
\text{Feature Importance} \\
$\mathbf{[0.004 \quad 0.0001 \quad \cdot \cdot \cdot \cdot]}$ \\
$N = 1 \ 699 \ 528$
\end{tabular}
$$

**Fig. 4.** (A) The types of predictors tested in this work (left to right): connectivity matrices (vectorized), macroscale cortical gradients, neighborhood, and centroid dispersion. (B) All four types of features were concatenated together (2) and decomposed using group PCA (2) (each feature group is decomposed separately). The resulting dataset, along with covariates, was divided into the train and holdout set; 10-fold cross-validation (CV) was used to assess the performance of L2-regularized logistic regression on the PCA dataset (3). (C) Permutation component importance was computed for each component using the holdout set (1). For each feature type, component importance was inverse transformed to obtain feature importance (2). COVARS: covariates (age, sex, framewise displacement, site), CONN: connectivity, CV: cross-validation, DISPcnr: centroid dispersion, DISPnbr: neighborhood dispersion, GRAD: cortical gradients, I<sub>PCA</sub>: component permutation importance, L2-Log Reg: L2-regularized logistic regression, PCA: principal component analysis.

Classifier analysis

Finally, we assessed the predictive capacity of each feature type in a classifier-agnostic manner to prevent the
results from being driven by the choice of a specific classifier. To this end, we selected the top 1% of features with
the largest permutation feature importance and we trained and tested 13 distinct classifiers on them:
* Logistic regression (L2-regularized, LR).

* K-Neighbors Classifier (KN).

* Naïve Bayes (NB).

* Decision Tree Classifier (DT).

* Support Vector Machine (SVM).

* Ridge Classifier (Ridge).

* Random Forest Classifier (RF).

* Ada Boost Classifier (AB).

* Gradient Boosting Classifier (GB).

* Light Gradient Boosting Machine (LGB).

* Linear Discriminant Analysis (LDA).

* Extra Trees Classifier (ET).

* Quadratic Discriminant Analysis (QDA).

To verify that the difference in classification performance between feature types persists regardless of the number of features, we repeated this analysis for a range of 100 to 10,000 features with the largest permutation feature importance from each type.

For each feature count, we identified the best classifier based on its mean cross-validation (CV) accuracy across 10 folds. Next, we tested the best classifier on the holdout set. All data transformations were done within Scikit-Learn pipelines, i.e., they were performed separately on the test and holdout sets. The multi-classifier analysis was done using Pycaret (https://github.com/pycaret/pycaret). Age, sex, site and framewise displacement (FD) were always included as covariates.

For classifier performance, we report both accuracy and the F1-score. The F1-score, calculated using the Scikit-Learn package in Python, represents the harmonic mean of precision and recall. This metric is particularly useful in cases where class distribution is imbalanced, as it accounts for both false positives and false negatives. We report accuracy and F1-score for cross-validation (CV) folds and the held out sets. In CV, each fold consists of a training and validation set. Thus, CV metrics were computed on the validation set of each fold. For instance, mean CV accuracy is the mean across accuracies computed on validation sets of 10 folds. Regarding test performance, whenever we refer to test F1-score or accuracy, we refer to the performance on the subjects that the model has not seen during training, i.e., these subjects were not present in any of the 10 folds the model was trained and validated on.

The CV and test performance across all classifiers and feature subsets was compared to two baselines, namely:

1) Test performance of the logistic regression trained on the PCA dataset (all 149 components) on the hold-out set.

2) The dummy classifier which randomly picks the class for each sample. It is frequently used as a baseline in machine learning research$^{46,47}$. Note that in this study the accuracy of the dummy classifier consistently remained at 73.7%. In contrast, the F1 score for this classifier is always 0, meaning that it does not differentiate between the two classes. These values constitute our chance reference. We also conducted the multi-classifier analysis for the principal components of the feature types (for the full ranking of the models trained on the components of features see Supplementary Table 2). Logistic regression was ranked 2nd for accuracy and 1st for F1-score.

## Localization of the best features in the brain: weighted degree centrality

The analyses described above allow us to identify the features exerting the largest influence on model predictions. However, it is also necessary to establish their links to specific brain areas. To this end, we conducted an additional exploratory analysis. Since raw connectivity was identified as the best feature type, this analysis was tailored for connectivity edges.

For three feature subsets including 500, 1000 and 5000 connectivity edges with the largest permutation feature importance, we computed weighted degree centrality (WDC). For each subset, the selected edges were transformed back to a 1000 × 1000 connectivity matrix and for each region the sum of edges was computed, each edge weighted by its correlation value. The computation of WDC can be mathematically formulated in what follows. Given a connectivity matrix $C_{subj}$ of a specific subject, let $E_{selected}$ be a set of indices of the edges with the highest permutation feature importance. In $C_{subj}$ we mask the edges whose indices are not in $E_{selected}$:

$$C_{subj,selected} = M \odot C_{subj} \text{ where } M = \begin{cases} 1, & \text{if } (i, j) \in E_{selected} \text{ or } (j, i) \in E_{selected} \\ 0, & \text{otherwise} \end{cases} \quad (4)$$

Next, for neurotypical and individuals with schizophrenia, we compute mean WDC for every region using the unmasked correlation values:

$$WDC_{group}(C_{subj,selected}) = \frac{1}{S} \sum_{subj=1}^{S} \sum_{j=1}^{1000} c_{i,j}; \quad c_{i,j} \in C_{subj,selected} \quad (5)$$

where S is the number of subjects in the group. Of note, those edges which were not selected based on their feature importance ($(h, m) \notin E_{selected}$) were set to 0 in $C_{subj,selected}$ as per Eq. 4. In summary, for each region, we sought to quantify its connectivity strength to the other regions given the selected edges.
# Results

## Permutation feature importance and classification performance

Upon visual inspection, we observed that the principal components of functional connectivity had the largest permutation feature importance, followed by gradients, centroid dispersion, and neighborhood dispersion (Fig. 5A). We sought to verify that permutation feature importance indeed reflects an advantage in classification performance regardless of classifier. To this end, we selected the top 1% of features with the largest permutation feature importance for each type and trained and tested 13 classifiers on them (Fig. 5). Connectivity outperformed the other feature types. This conclusion was supported by the Mann-Whitney U test (Fig. 5B, Table 2) for both accuracy and F1-score. To limit the number of comparisons, we only performed tests for connectivity vs. all other types which amounted to 5 tests for accuracy and F1-score. The results of all comparisons are displayed in Table 2.

Multi-classifier analyses further emphasized the superior predictive capacity of functional connectivity compared to the other feature types. Here we report CV and test accuracy and F1-score across all classifiers and for the best classifier (Fig. 5C). Connectivity edges with the largest permutation feature importance consistently outperformed logistic regression fitted on all feature components (Supplementary Table 2). The other feature types performed substantially worse.

In addition, we examined the impact of the number of best features on the choice of the best classifier (fits on all feature types are considered). Figure 5D illustrates the evolution of the best classifier with the increasing number of features as the change in relative density of instances where the classifiers were identified as best. Three main patterns can be noted. Firstly, several classifiers clearly performed better with N_features < 3000: GB, SVM, RF, KN and Ridge. Secondly, ET and LGB had a relatively stable winning rate across all feature subsets. Thirdly, LR and LDA had an increased performance when N_features > 3000. However, LDA rarely outperformed the other classifiers, whereas LR and LGB often emerged as optimal.

## The features with the largest permutation feature importance

For this analysis, we selected 500, 1000 and 5000 connectivity edges with the largest permutation feature importance and computed weighted degree centrality (WDC, Eq. 5. For each edge subset, we plotted the difference in group WDC between SCZ and NC (Fig. 6):

$$WDC_{diff} = WDC_{SCZ} - WDC_{NC} \qquad (6)$$

WDC was overall lower in subjects with schizophrenia which is in line with previous accounts of lower overall functional connectivity characteristic of this disease<sup>48–53</sup>. This finding lends support to the dysconnectivity hypothesis of schizophrenia<sup>54,55</sup>. Furthermore, the spatial pattern illustrated in Fig. 6 indicated that the significant differences in WDC were concentrated in the primary regions for all edge subsets displayed here. Put differently, the edges with the largest permutation feature importance reflect connectivity strength in the primary regions, indicating that their connectivity profiles contribute the most to classification performance. The edges selected for this analysis can be viewed in Supplementary Material 1, 2 and 3 for 500, 1000 and 5000 most important edges respectively (the edges that were not selected are set to 0).

Additionally, we conducted a Mann-Whitney U test on WDC. The results align with the observations outlined above and can be viewed in Supplementary Fig. 2. However, we acknowledge that in the context of our feature selection procedure this test may be viewed as circular. Therefore, its results are to be interpreted with caution.

# Discussion

In this paper, we extend the effort to determine the optimal connectivity-based predictors of schizophrenia by attempting to benchmark connectivity-based features against each other for diagnosis prediction. To this end, we applied a feature selection workflow based on permutation feature importance to a large dataset comprising connectivity, macroscale cortical gradients and gradient dispersion. Our analysis revealed that, despite growing interest in cortical gradients and gradient dispersion as potential biomarkers of schizophrenia and psychosis<sup>18,19,56</sup>, functional connectivity holds superior predictive potential over its low-dimensional derivatives (Fig. 5B,C). Furthermore, we showed that this result remains consistent even when compared to the performance of a larger number of features of a different type. Specifically, the top 1% of connectivity features with the highest permutation feature importance performed better than the top 1% of all features and neighborhood dispersion features. Recently, an investigation similar to ours has been conducted for subject task performance<sup>31</sup>. Therein, the authors reported that individual parcellations achieved better predictions compared to cortical gradients. Although we employed a group parcellation, our findings align with this result.

Additionally, we demonstrated that the connectivity edges connecting the primary sensory regions have the largest permutation feature importance. This result indicates that variations in connectivity strength of these regions encapsulate critical information for distinguishing people with schizophrenia from neurotypical controls. This finding, however, does not discard the relevance of the gradients for the studies looking specifically into functional hierarchical variations. Regarding our use case, the diminished predictive capacity of the gradients could stem from the conservative matrix thresholding applied prior to dimensionality reduction. Future studies are needed to test this hypothesis.

To gain a deeper understanding of the connectivity patterns associated with schizophrenia, we conducted an exploratory analysis whereby we compared weighted degree centrality for the regions linked by the 500, 1000 and 5000 most important edges. Firstly, degree centrality appeared to be lower in individuals with schizophrenia,
<figure><img src="image_10.png" /><figcaption>**B** Feature Type Comparison</figcaption></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**C** Model Performance</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**D** The Best Model vs. Number of Features</figcaption></figure>

<figure><img src="image_14.png" /></figure>

**Fig. 5.** (A) Permutation importance across feature types. (B) Accuracy and F1 score across 13 classifiers (mean cross CV folds; dummy classifier was excluded) fit on top 1% best features from each feature type, the principal gradient, the 28 values of centroid dispersion and the top 1% best features from the whole feature set (mixed: the inset shows the number of features from each feature type that were included in the top 1%). P-values indicate significant difference as per Mann-Whitney U test, $α ≤ 0.01$ (connectivity vs. all: Bonferroni-corrected). The stars denote the performance of the best classifier as identified based on the mean accuracy across 10 CV folds. (C) Mean ± s.e.m. CV and test performance across classifiers for N features 100–10,000 for connectivity (left), gradients (middle), and neighborhood dispersion (right). Horizontal lines represent test performance of the logistic regression on all principal components (blue), and the performance of the dummy classifier (brown). The shading indicates s.e.m. (D) Relative density of fits where the corresponding classifier was identified as best. Larger area indicates that the corresponding classifier had the highest CV accuracy more often. The legend features all classifiers that were tested in this study; the classifiers in black were never identified as the best. CV: cross-validation; PCA: principal component analysis; SVM: support vector machine.
<table><caption>**Table 2.** Results of Mann-Whitney U test. We compared the performance of the top 1% of the most important features from connectivity to all the other types. Significant results are marked with an asterisk. The results were corrected for multiple comparisons (Bonferroni: α = 0.01). These results are also represented graphically in a boxplot in Fig. 5B.</caption><thead><tr><th>Top 1%</th><th>Gradient 1</th><th>Neighborhood dispersion</th><th>Centroid dispersion</th><th>All gradients</th><th>All features</th><th>Metric</th></tr></thead><tbody><tr><td rowspan="2">Connectivity</td><td>U = 142,<br/>p = 0.003*</td><td>U = 141,<br/>p = 0.004*</td><td>U = 128,<br/>p = 0.027</td><td>U = 127,<br/>p = 0.03</td><td>U = 83,<br/>p = 0.95</td><td>Accuracy</td></tr><tr><td>U = 154,<br/>p = 0.0004*</td><td>U = 153,<br/>p = 0.0004*</td><td>U = 137,<br/>p = 0.008*</td><td>U = 126,<br/>p = 0.04</td><td>U = 85,<br/>p = 1.0</td><td>F1-score</td></tr></tbody></table>

Weighted Degree Centrality: SCZ vs. NC

<figure><img src="image_15.png" /><figcaption>**Fig. 6.** The difference in mean weighted degree centrality (WDC, averaged across subjects, Eq. 5) between the two groups for 500, 1000 and 5000 connectivity edges with the largest permutation feature importance. Inset violinplots display WDC averaged across regions for the two groups. Color bars denote the difference in WDC between SCZ and NC. SCZ: subjects diagnosed with schizophrenia, NC: neurotypical controls, WDC: weighted degree centrality.</figcaption></figure>

corroborating previous evidence of overall hypoconnectivity typical for the disorder$^{48–53}$. Secondly, for the edges with the largest importance, these differences were predominantly concentrated in the sensorimotor, auditory, and visual cortex. These findings appear to contrast with the studies highlighting discrepancies predominantly in higher-order areas, such as the default mode network (DMN)$^{57–61}$. However, the differences we observed need to be contextualized as the ones relevant for classification performance: while they may or may not constitute the neurological basis of schizophrenia, they are most informative for the differentiation of the two groups. Viewed from this perspective, our results resonate with the research demonstrating that across individuals the measurements are most reliable in the primary, unimodal areas$^{62,63}$. In addition, one study reported more accurate surface registration for the primary sensory areas$^{64}$. It has been postulated that this stability can be attributed to the fact that the primary areas are phylogenetically the most ancient cortical areas and are therefore considered as evolutionary anchors around which most of the cortical expansion in humans unfolded$^{63,65–68}$. The largest feature importance in the primary and most stable regions may indicate that the classifiers rely on those individual differences which are replicable across individuals of the same group. In such a setting, intra-individual variability cannot be accounted for. This shortcoming can be potentially addressed by studies focusing on datasets featuring hours of scanning time per individual$^{69–71}$.

We were also interested in how the ranking of the 13 classifiers tested in this work evolved with the number of best features selected based on permutation feature importance. As depicted in Fig. 5D, some classifiers' performance relative to each other fluctuated considerably depending on the size of the feature subset. Overall, no classifier outperformed others at all times for any number of features. It is a clear indication of the necessity of empirical testing of the candidate classifiers. However, it appears that as the number of features increases, the performance of linear classifiers — such as logistic regression and LDA—improves. Prior work offers insight into why linear classifiers would outperform non-linear approaches when the number of features disproportionately exceeds the number of training samples. Ardeshir et al.$^{72}$ benchmarked the performance of support vector machine (SVM) and logistic regression in a high-dimensional (overparametrized) setting. The authors demonstrate that the solutions provided by SVM (which involves non-linear operations) and a linear
classifier such as logistic regression coincide when the number of model parameters is large. This trend is the result of *support vector proliferation*<sup>73</sup>. In this scenario, SVM does not provide a performance boost compared to logistic regression. This phenomenon is in line with the evolution of our model ranking as a function of the number of predictors included in the model. Specifically, it is possible that adding more parameters to the model offsets the benefit of non-linear classifiers compared to their linear analogues. This account could explain the increase in the number of instances for higher numbers of features when linear models reached the top of the performance leaderboard.

Clinical neuroscientists frequently face the imbalance between sample size and the number of features. The choice of features with the highest predictive potential is critical in this line of work. With this problem in mind, we developed a pipeline which allowed us to select the most promising features from over one million candidates for diagnosis prediction. Based on PCA applied separately to each feature group and permutation feature importance, our approach enabled a systematic evaluation of the predictive capacity of each group, identifying connectivity as the winner. Notably, this result was achieved without testing all possible combinations of feature types which would be computationally challenging.

A number of limitations need to be noted. First, most subjects in our sample underwent scans lasting less than 6 minutes. Short scan time such as this prevented us from deriving individual parcellations<sup>74</sup> which, as it has been demonstrated<sup>74,75</sup>, account of individual neuroanatomy and enrich functional connectivity with intra-individual variance. This variance was unavailable in our case since we relied on a group parcellation. In addition, prior investigations have articulated that cortical gradients are most discriminable across subjects when scan duration exceeds 30 minutes<sup>16</sup>. Given that our circumstance did not align with this recommendation, a certain degree of caution is warranted when interpreting our findings.

Second, we did not account for the effect of medication in our study since the medication data were missing for a large number of participants diagnosed with schizophrenia. Excluding these subjects would have resulted in a drastically decreased sample size and, more importantly, in an increased imbalance of classes. We note however that even without these covariates, classification accuracy reached 80% and beyond for the features with the highest permutation feature importance.

Third, there were significant differences between both diagnosis groups and datasets in terms of age, sex distribution and framewise displacement. To address this limitation, we included all these variables as covariates in all models. Alternatively, cross-site harmonization could be applied to remove site-specific variance from the data. Various harmonization methods are available, ranging from linear regression to techniques more specific to MRI data such as ComBat<sup>76,77</sup>. However, all harmonization approaches known to date present a significant risk of reduction in sensitivity. That is, biologically relevant information may be regressed out along with site-related signals. In this work, we prioritized avoiding this risk.

Fourth, most classifiers in this study were fitted on datasets with more features than observations which increases the risk of overfitting. Nonetheless, we believe to have addressed this issue by employing a 10-fold CV for each classifier. In addition, we have also witnessed that for some classifiers test performance exceeded CV performance (Fig. 5C) which is not characteristic of overfitting.

Fifth, while permutation feature importance allowed us to select the most important features from multiple types, it was computed solely for logistic regression. Although feature rankings might vary depending on the model used to compute feature importances, we deem the results of this study robust for two reasons. First, based on cross-validation, logistic regression was rarely identified as the best model for feature sets containing fewer than 3000 features. That is, even though the features were ranked using logistic regression importances, this model did not outperform others for smaller feature subsets. Second, the usage of logistic regression for feature selection is a well-known strategy and proven effective in several domains<sup>78,79</sup>. Nonetheless, a comprehensive comparison of feature rankings generated by different models would be valuable and should be explored in future work.

In addition, we did not consider the issue of comorbidity, transdiagnostic phenomena across psychiatric disorders or subtypes of schizophrenia. While these aspects lie beyond the purview of the present study, it is imperative to acknowledge their critical significance. Inquiries into transdiagnostic classification, symptom prediction and subtyping hold potential to profoundly transform the paradigms governing the diagnosis and treatment of psychiatric and neurodevelopmental disorders. Indeed, as of late, the research looking into transdiagnostic effects<sup>7,56,80-84</sup> and schizophrenia subtypes<sup>85-87</sup> has garnered substantial momentum. We expect our work to assist future endeavors investigating these outstanding questions in benchmarking biomarker candidates.

The emergence of novel connectivity-based methods broadens our toolkit for predicting psychiatric disorders, introducing a necessity for empirical validation. Our findings indicate that functional connectivity outperforms its more recent, low-dimensional derivatives such as cortical gradients and gradient dispersion in predicting schizophrenia. Additionally, in this study, the connectivity within the primary sensory regions showed the highest discrimination capabilities, possibly due to the reduced anatomical and functional variability of those regions. We anticipate, however, that it is also informative for a broad spectrum of major psychiatric disorders. The exploration of this latter possibility warrants thorough examination in future work.

## Data availability

The data used in this work are in public access; for COBRE: http://schizconnect.org/, SPBRS-1600: https://bicr-resource.atr.jp/srpbopen/, and LA5c: https://openfmri.org/dataset/ds000030/. A user account on http://schizconnect.org/ is required to download COBRE. Access to SPBRS-1600 is conditioned upon completing the access application (https://bicr-resource.atr.jp/bicr_add/jsp/appForm.jsp) and signing a data sharing agreement. For questions about access to the SPBRS-1600 dataset, please contact deconf-db-admin@atr.jp or the corresponding author of this article: https://doi.org/https://doi.org/10.1038/s41597-021-01004-8.
## Code availability

The code used to perform the analyses and to produce the figures featured in this study is available Github: https://github.com/victoris93/feature-selection-scz. All figures are generated using Seaborn (v. 0.13.2: https://seaborn.pydata.org/), Surfplot (v. 0.2.0: https://surfplot.readthedocs.io/) and Matplotlib (v. 3.9.3: https://matplotlib.org/) packages in Python. When necessary, the figures were compiled in panels using Affinity Designer (https://affinity.serif.com/, not open source).

Received: 16 July 2024; Accepted: 20 December 2024

Published online: 22 January 2025

## References

1. Yang, G. J. et al. Functional hierarchy underlies preferential connectivity disturbances in schizophrenia. Proc. Natl. Acad. Sci. U. S. A. 113, E219-E228 (2016).

2. Chan, Y. H., Yew, W. C., Chew, Q. H., Sim, K. & Rajapakse, J. C. Elucidating salient site-specific functional connectivity features and site-invariant biomarkers in schizophrenia via deep neural networks. Sci. Rep. 13, 21047 (2023).

3. Gifford, G. et al. Resting state fMRI based multilayer network configuration in patients with schizophrenia. Neuroimage Clin. 25, 102169 (2020).

4. Sendi, M. S. E. et al. Aberrant dynamic functional connectivity of default mode network in schizophrenia and links to symptom severity. Front. Neural Circuits 15, 649417 (2021).

5. Li, S. et al. Dysconnectivity of multiple brain networks in schizophrenia: A meta-analysis of resting-state functional connectivity. Front. Psychiatry 10, 482 (2019).

6. Lynall, M.-E. et al. Functional connectivity and brain networks in schizophrenia. J. Neurosci. 30, 9477-9487 (2010).

7. Ma, Q. et al. Transdiagnostic dysfunctions in brain modules across patients with schizophrenia, bipolar disorder, and major depressive disorder: A connectome-based study. Schizophr. Bull. 46, 699-712 (2020).

8. Cearns, M., Hahn, T. & Baune, B. T. Recommendations and future directions for supervised machine learning in psychiatry. Transl. Psychiatry 9, 271 (2019).

9. Cui, Z. & Gong, G. The effect of machine learning regression algorithms and sample size on individualized behavioral prediction with functional connectivity features. Neuroimage 178, 622-637 (2018).

10. Serin, E., Zalesky, A., Matory, A., Walter, H. & Kruschwitz, J. D. NBS-Predict: A prediction-based extension of the network-based statistic. Neuroimage 244, 118625 (2021).

11. Jung, H., Wager, T. D. & Carter, R. M. Novel cognitive functions arise at the convergence of macroscale gradients. J. Cogn. Neurosci. 34, 381-396 (2022).

12. Margulies, D. S. et al. Situating the default-mode network along a principal gradient of macroscale cortical organization. Proc. Natl. Acad. Sci. U. S. A. 113, 12574-12579 (2016).

13. Bethlehem, R. A. I. et al. Dispersion of functional gradients across the adult lifespan. Neuroimage 222, 117299 (2020).

14. Pasquini, L. et al. Dysfunctional cortical gradient topography in treatment-resistant major depressive disorder. Biol. Psychiatry Cogn. Neurosci. Neuroimaging (2022).

15. Paquola, C. et al. Shifts in myeloarchitecture characterise adolescent development of cortical gradients. Elife (2019).

16. Hong, S.-J. et al. Toward a connectivity gradient-based framework for reproducible biomarker discovery. Neuroimage 223, 117322 (2020).

17. Knodt, A. R. et al. Test-retest reliability and predictive utility of a macroscale principal functional connectivity gradient. Hum. Brain Mapp. 44, 6399-6417 (2023).

18. Dong, D. et al. Compressed sensorimotor-to-transmodal hierarchical organization in schizophrenia. Psychol. Med. (2021).

19. Holmes, A. et al. Disruptions of hierarchical cortical organization in early psychosis and schizophrenia. Biol. Psychiatry Cogn. Neurosci. Neuroimaging 8, 1240-1250 (2023).

20. Wang, J., Zhou, Y., Ding, J. & Xiao, J. Functional gradient alteration in individuals with cognitive vulnerability to depression. J. Psychiatr. Res. 144, 338-344 (2021).

21. Xia, M. et al. Connectome gradient dysfunction in major depression and its association with gene expression profiles and treatment outcomes. Mol. Psychiatry 27, 1384-1393 (2022).

22. Xiao, Y. et al. Charting the dorsal-medial functional gradient of the default mode network in major depressive disorder. J. Psychiatr. Res. 153, 1-10 (2022).

23. Wan, B. et al. Diverging asymmetry of intrinsic functional organization in autism. Mol. Psychiatry 28, 4331-4341 (2023).

24. Yang, Y. et al. Altered asymmetry of functional connectome gradients in major depressive disorder. Front. Neurosci. 18, (2024).

25. Park, B.-Y. et al. Connectome-wide structure-function coupling models implicate polysynaptic alterations in autism. Neuroimage 285, 120481 (2024).

26. Nguyen, T. T. et al. Variations in cortical functional gradients relate to dimensions of psychopathology in preschool children. J. Am. Acad. Child Adolesc. Psychiatry 63, 80-89 (2024).

27. He, Y. et al. Functional gradients reveal altered functional segregation in patients with amnestic mild cognitive impairment and Alzheimer's disease. Cereb. Cortex 33, 10836-10847 (2023).

28. Hu, Q., Li, Y., Wu, Y., Lin, X. & Zhao, X. Brain network hierarchy reorganization in Alzheimer's disease: A resting-state functional magnetic resonance imaging study. Hum. Brain Mapp. 43, 3498-3507 (2022).

29. Bassett, D. S. et al. Hierarchical organization of human cortical networks in health and schizophrenia. J. Neurosci. 28, 9239-9248 (2008).

30. Setton, R. et al. Age differences in the functional architecture of the human brain. Cereb. Cortex 33, 114-134 (2022).

31. Kong, R. et al. Comparison between gradients and parcellations for functional connectivity prediction of behavior. Neuroimage 273, 120044 (2023).

32. Gao, Z. et al. Flexing the principal gradient of the cerebral cortex to suit changing semantic task demands. Elife 11, (2022).

33. Aine, C. J. et al. Multimodal neuroimaging in schizophrenia: Description and dissemination. Neuroinformatics 15, 343-364 (2017).

34. Poldrack, R. A. et al. A phenome-wide examination of neural and cognitive function. Sci. Data 3, 160110 (2016).

35. Tanaka, S. C. et al. A multi-site, multi-disorder resting-state magnetic resonance image database. Sci Data 8, 227 (2021).

36. Esteban, O. et al. fMRIPrep: A robust preprocessing pipeline for functional MRI. Nat. Methods 16, 111-116 (2019).

37. Gorgolewski, K. et al. Nipype: A flexible, lightweight and extensible neuroimaging data processing framework in python. Front. Neuroinform. 5, 13 (2011).

38. Schaefer, A. et al. Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity MRI. Cereb. Cortex 28, 3095-3114 (2018).

39. Girn, M. et al. Serotonergic psychedelic drugs LSD and psilocybin reduce the hierarchical differentiation of unimodal and transmodal cortex. Neuroimage 256, 119220 (2022).

40. Huntenburg, J. M., Bazin, P.-L. & Margulies, D. S. Large-scale gradients in human cortical organization. Trends Cogn. Sci. 22, 21-31 (2018).
41. Wang, C. & Mahadevan, S. Manifold alignment using Procrustes analysis. In *Proceedings of the 25th international conference on Machine learning - ICML '08* (ACM Press, New York, New York, USA, 2008).

42. Van Essen, D. C. et al. The human connectome project: A data acquisition perspective. *Neuroimage* **62**, 2222–2231 (2012).

43. De Rosa, A. P. et al. Functional gradients reveal cortical hierarchy changes in multiple sclerosis. *Hum. Brain Mapp.* **45**, e26678 (2024).

44. Breiman, L. Random forests. *Mach. Learn.* **45**, 5–32 (2001).

45. Gordon, A. D., Breiman, L., Friedman, J. H., Olshen, R. A. & Stone, C. J. Classification and regression trees. *Biometrics* **40**, 874 (1984).

46. Alzaid, M. A comparative study between classical feature engineering and RNNs for seizure detection in imbalanced data. In *2020 IEEE Signal Processing in Medicine and Biology Symposium (SPMB)* (IEEE, 2020).

47. Viglione, A. et al. Behavioral impulsivity is associated with pupillary alterations and hyperactivity in CDKL5 mutant mice. *Hum. Mol. Genet.* **31**, 4107–4120 (2022).

48. Cai, M. et al. Disrupted local functional connectivity in schizophrenia: An updated and extended meta-analysis. *Schizophrenia (Heidelb)* **8**, 93 (2022).

49. Du, X. et al. Unraveling schizophrenia replicable functional connectivity disruption patterns across sites. *Hum. Brain Mapp.* **44**, 156–169 (2023).

50. Liu, Y. et al. Disrupted small-world networks in schizophrenia. *Brain* **131**, 945–961 (2008).

51. van Dellen, E. et al. Functional brain networks in the schizophrenia spectrum and bipolar disorder with psychosis. *NPJ Schizophr.* **6**, 22 (2020).

52. Zalesky, A., Fornito, A. & Bullmore, E. T. Network-based statistic: identifying differences in brain networks. *Neuroimage* **53**, 1197–1207 (2010).

53. Zhuo, C. et al. Brain structural and functional dissociated patterns in schizophrenia. *BMC Psychiatry* **17**, 45 (2017).

54. Friston, K., Brown, H. R., Siemerkus, J. & Stephan, K. E. The dysconnection hypothesis (2016). *Schizophr. Res.* **176**, 83–94 (2016).

55. Zhou, Y., Fan, L., Qiu, C. & Jiang, T. Prefrontal cortex and the dysconnectivity hypothesis of schizophrenia. *Neurosci. Bull.* **31**, 207–219 (2015).

56. Park, B.-Y. et al. Multiscale neural gradients reflect transdiagnostic effects of major psychiatric conditions on cortical morphology. *Commun Biol* **5**, 1024 (2022).

57. Krishnan, R. R., Keefe, R. & Kraus, M. Schizophrenia is a disorder of higher order hierarchical processing. *Med. Hypotheses* **72**, 740–744 (2009).

58. Kong, M. et al. Abnormal network homogeneity of default-mode network and its relationships with clinical symptoms in antipsychotic-naïve first-diagnosis schizophrenia. *Front. Neurosci.* **16**, 921547 (2022).

59. Roiser, J. P. et al. Dysconnectivity in the frontoparietal attention network in schizophrenia. *Front. Psychiatry* **4**, (2013).

60. Tu, P.-C., Lee, Y.-C., Chen, Y.-S., Li, C.-T. & Su, T.-P. Schizophrenia and the brain's control network: aberrant within- and between-network connectivity of the frontoparietal network in schizophrenia. *Schizophr. Res.* **147**, 339–347 (2013).

61. Wang, H. et al. Evidence of a dissociation pattern in default mode subnetwork functional connectivity in schizophrenia. *Sci. Rep.* **5**, 14655 (2015).

62. Gratton, C. et al. Functional brain networks are dominated by stable group and individual factors, not cognitive or daily variation. *Neuron* **98**, 439–452.e5 (2018).

63. Mueller, S. et al. Individual variability in functional connectivity architecture of the human brain. *Neuron* **77**, 586–595 (2013).

64. Tardif, C. L. et al. Multi-contrast multi-scale surface registration for improved alignment of cortical areas. *Neuroimage* **111**, 107–122 (2015).

65. Yeo, B. T. T. et al. The organization of the human cerebral cortex estimated by intrinsic functional connectivity. *J. Neurophysiol.* **106**, 1125–1165 (2011).

66. Buckner, R. L. & Krienen, F. M. The evolution of distributed association networks in the human brain. *Trends Cogn. Sci.* **17**, 648–665 (2013).

67. Preuss, T. M. & Wise, S. P. Evolution of prefrontal cortex. *Neuropsychopharmacology* **47**, 3–19 (2022).

68. Heuer, K. et al. Evolution of neocortical folding: A phylogenetic comparative analysis of MRI from 34 primate species. *Cortex* **118**, 275–291 (2019).

69. Paugam, F., Pinsard, B., Lajoie, G. & Bellec, P. A benchmark of individual auto-regressive models in a massive fMRI dataset. *Imaging Neuroscience* **2**, 1–23 (2024).

70. Rastegarnia, S., St-Laurent, M., DuPre, E., Pinsard, B. & Bellec, P. Brain decoding of the human connectome project tasks in a dense individual fMRI dataset. *Neuroimage* **283**, 120395 (2023).

71. Pinho, A. L. et al. Individual brain charting dataset extension, second release of high-resolution fMRI data for cognitive mapping. *Sci. Data* **7**, 353 (2020).

72. Ardeshir, N., Sanford, C. & Hsu, D. Support vector machines and linear regression coincide with very high-dimensional features. *Adv. Neural Inf. Process. Syst.* 4907–4918 (2021).

73. Hsu, D., Muthukumar, V. & Xu, J. On the proliferation of support vectors in high dimensions. *J. Stat. Mech.* **2022**, 114011 (2022).

74. Gordon, E. M. et al. Precision functional mapping of individual human brains. *Neuron* **95**, 791–807.e7 (2017).

75. Braga, R. M. & Buckner, R. L. Parallel interdigitated distributed networks within the individual estimated by intrinsic functional connectivity. *Neuron* **95**, 457–471.e5 (2017).

76. Carré, A. et al. AutoComBat: A generic method for harmonizing MRI-based radiomic features. *Sci. Rep.* **12**, 12762 (2022).

77. Wang, Y.-W., Chen, X. & Yan, C.-G. Comprehensive evaluation of harmonization on functional brain imaging for multisite data-fusion. *Neuroimage* **274**, 120089 (2023).

78. Wang, J., Wang, H., Nie, F. & Li, X. Feature selection with multi-class logistic regression. *Neurocomputing* **543**, 126268 (2023).

79. Khandezamin, Z., Naderan, M. & Rashti, M. J. Detection and classification of breast cancer using logistic regression feature selection and GMDH classifier. *J. Biomed. Inform.* **111**, 103591 (2020).

80. Buckholtz, J. W. & Meyer-Lindenberg, A. Psychopathology and the human connectome: toward a transdiagnostic model of risk for mental illness. *Neuron* **74**, 990–1004 (2012).

81. Hettwer, M. D. et al. Coordinated cortical thickness alterations across six neurodevelopmental and psychiatric disorders. *Nat. Commun.* **13**, 6851 (2022).

82. Elliott, M. L., Romer, A., Knodt, A. R. & Hariri, A. R. A connectome-wide functional signature of transdiagnostic risk for mental illness. *Biol. Psychiatry* **84**, 452–459 (2018).

83. Segal, A. et al. Regional, circuit and network heterogeneity of brain abnormalities in psychiatric disorders. *Nat. Neurosci.* **26**, 1613–1629 (2023).

84. Xia, C. H. et al. Linked dimensions of psychopathology and connectivity in functional brain networks. *Nat. Commun.* **9**, 3003 (2018).

85. Ivleva, E. I., Turkozer, H. B. & Sweeney, J. A. Imaging-based subtyping for psychiatric syndromes. *Neuroimaging Clin. N. Am.* **30**, 35–44 (2020).

86. Jiang, Y. et al. Neuroimaging biomarkers define neurophysiological subtypes with distinct trajectories in schizophrenia. *Nat. Mental Health* **1**, 186–199 (2023).
87. Chen, J. et al. Neurobiological divergence of the positive and negative schizophrenia subtypes identified on a new factor structure of psychopathology using non-negative factorization: An international machine learning study. *Biol. Psychiatry* **87**, 282–293 (2020).

**Acknowledgements**

We thank Joshua Vogelstein for methodological guidance and Vincent Kreft for contributing to the conceptual development of the study. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 866533) awarded to D.S.M. This project was also supported by the contrat doctoral spécifique normalien (CDSN) awarded to V.S. by the Ministry of Higher Education and Research of France. The computational aspects of this research were supported by the Wellcome Trust Core Award Grant Number 203141/Z/16/Z and the NIHR Oxford BRC. The views expressed are those of the author(s) and not necessarily those of the NHS, the NIHR or the Department of Health.

**Author contributions**

V.S. researched the data, conducted the analyses, generated the figures and wrote the manuscript. C.P., D.W., P.L.B., and D.S.M provided methodological expertise. C.P., R.A.B, R.S., W.W., U.K., F.A., D.W., T.D.S., P.L.B., and D.S.M. assisted in interpretation of findings and manuscript revision.

**Declarations**

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-024-84152-2.

**Correspondence** and requests for materials should be addressed to V.S. or D.S.M.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2025