Data acquisition and preprocessing
The REST-meta-MDD project (http://rfmri.org/REST-meta-MDD) was the primary source of all collected data. Brain imaging data were gathered from 1300 patients with major depressive disorder and 1128 HC subjects, sourced from 25 research groups in 17 hospitals across China17. Subsequently, the data underwent further screening to meet the specific requirements of this study. For detailed methods, please refer to our previous papers18,19. The main processing steps included removing duplicate non-2.0 time point data and data with zero time series for participants. Gender and age matching tests were then conducted for each site's data, and unmatched site data was excluded to ensure that all site data passed the age and gender matching tests. Next, we performed data preprocessing using DPARSF. The primary steps involved global signal regression, extraction of time series for each brain region based on anatomical automatic labeling (AAL) brain atlas, slice timing correction, realignment, spatial normalization, and nuisance regression. To ensure generalizability of results, we chose to define the nodes of the brain network using 116 regions of interest (ROIs) from the widely used AAL atlas20. For detailed information on data preprocessing procedures, please refer to our previous paper17. Finally, functional connectivity networks based on Pearson correlation were generated using BrainNetClass Toolkit21.
Network metrics
After the completion of functional connectivity network construction, significant functional connections were selected and the functional network was reconstructed through two-sample t-tests (p < 0.05 and p < 0.01) to explore whether differential functional connections affect the comparison results of topological brain network indicators between groups. Subsequently, Gretna software (version 2.0, www.nirc.org/projects/Gretna) was used to calculate several common global and nodal topological brain network indicators for two types of networks (original network, reconstructed network)22. Global indicators included global efficiency (Eglob), local efficiency (Eloc), shortest path length (Lp), clustering coefficient (Cp), and small-worldness (σ) based on 100 randomizations. Nodal-level indicators included nodal local efficiency, nodal clustering coefficient, nodal shortest path length, and nodal efficiency. In previous related studies, calculations were typically performed within a sparsity threshold range of 10–34%, with an increment interval of 1%, to ensure estimability of the threshold networks for the above-mentioned metrics while minimizing false edges23–25. In this study, all metrics were calculated within a two-stage sparsity threshold range: first calculating metric values at various thresholds within a wide range of sparsity thresholds (3%-91%, with a step size of 2%), then selecting an effective sparsity threshold range through statistical analysis of metric differences at each threshold(p < 0.05) to ensure valid estimation within the threshold range. Subsequently, within this effective threshold range, the area under the curve (AUC) for each network metric was computed as a summary scalar for global topological features. This approach is independent from single-threshold selection and sensitive to topological changes in brain diseases23–26. Additionally, in order to further delineate differences in the functional brain network topology between the MDD and HC groups, we retained the results of global and nodal structural measurements under single thresholds.
Statistics analysis
In the assessment of demographic and clinical characteristics of patients with severe depression and healthy controls, the normality of data distribution was evaluated using the Shapiro-Wilk test. Differences in age and education level were analyzed using non-parametric Mann-Whitney U test (for non-normally distributed data) or independent samples t-test (for normally distributed data). Gender comparison was conducted using χ2 test. To compare the statistical differences in brain network topological metrics between groups, non-parametric Mann-Whitney U test or independent samples t-test were initially used to examine the differences in global topological properties (including normalized clustering coefficient (γ), normalized shortest path length (λ), and σ values) between MDD and HC groups across a wide range of sparsity threshold levels (3%-91%, with a 2% increment), in order to determine an effective threshold range. Subsequently, within the determined effective threshold range, differences in AUC values for global topological properties (including γ, λ, σ, Eg, Eloc) between the two groups were examined. Age, gender, and education level were included as covariates with Bonferroni correction.
Receiver Operating Characteristic (ROC) analysis
The ROC curve27 is a graphical analysis tool that reflects the comprehensive indicators of sensitivity and specificity. It uses sensitivity or true positive rate (TPR) as the vertical axis and 1-specificity or false positive rate (TNR) as the horizontal axis. In this study, ROC curve area under the curve (AUC) was used to measure the discriminative ability of brain node topological properties between MDD and HC groups. A study emphasized the use of the area under the ROC curve for evaluating the diagnostic capability of a test. Despite its limitation to binary outcomes and disregard for continuity, it can still be utilized in decision theory by converting it into a simplified measure28. Similarly, brain functional networks' topological properties are continuous variables categorized into MDD patients and HC. The AUC calculation of the topological properties of the reconstructed brain FC network post two-sample t-test is depicted under the sparsity threshold in Parts 3 to 5 of Fig. 1. The global attribute part (Part 4 of Fig. 1) divides the AUC value at all thresholds for each sample's global topological property as a sample feature value. In the nodal part (Part 5 of Fig. 1), all average nodal topological property values across sparsity thresholds and functional network nodes serve as sample feature values. These feature values are categorized as 0 for MDD and 1 for HC to distinguish between the two populations. The true positive rate refers to correctly identifying non-diseased subjects among all non-diseased subjects, while the true negative rate refers to correctly identifying diseased subjects among all diseased subjects. Furthermore, there are two methods for selecting significantly different node topological property values and using their mean as sample feature: first, taking averages of node topological features at all sparsity thresholds; secondly, selecting nodes before averaging their values as sample features. Alternatively, conducting inter-group differential analysis on node topological property features at each sparsity threshold level before averaging selected nodes' values; finally taking means across all thresholds to obtain overall sample features.
Functional indicators calculations
ALFF is considered an effective method for detecting the spontaneous oscillation intensity of brain regions and reflecting the spontaneous activity of the brain in rs-fMRI BOLD signals29,30, Meanwhile, fALFF represents the relative contribution of specific low-frequency oscillations to the entire frequency range31, Regional homogeneity (ReHo), as a measure of regional synchronization in functional magnetic resonance imaging (fMRI) time series, has been applied in MDD-related research32, In this study, three functional indicators were calculated: ALFF, fALFF, and Reho. Firstly, each voxel's BOLD time series was transformed into the frequency domain using fast Fourier transform. Then, the square root of power spectrum was computed and averaged within a specified frequency range (0.01–0.08 Hz) for each voxel to obtain ALFF. fALFF represents a normalized measure of ALFF by dividing its value by the total amplitude within the entire frequency range, effectively suppressing physiological artifacts31. Reho measures local synchronous activity by calculating voxel-wise similarity in fluctuation during a given voxel's time course under resting state conditions16, the calculation involves computing Kendall’s coefficient of concordance for brain voxels using a neighborhood size of 26 voxels.
Efficiency testing of topological properties
To further investigate the clinical application efficiency of topological features in identifying MDD, we utilized a support vector regression (SVM) model based on linear kernel33, which is a supervised learning model used for classification and regression analysis. In order to prevent test data leakage during the model training process, the reconstruction of sample functional networks and the selection of metric features of sample functional networks (including local efficiency of brain region nodes, functional indicators such as ALFF, and functional connectivity features) were all conducted within the training dataset. The least absolute shrinkage and selection operator (LASSO)34 was employed as the feature selection method in this study. LASSO is a common linear regression regularization method suitable for feature selection35, After filtering out some brain nodes, the topological features of the training set's brain nodes were input into the SVM model, as shown in part 6 of Fig. 1. To ensure robustness and reduce bias in the experiment, we adopted leave-one-site-out cross validation (LOSCV) strategy to validate the model. The performance measurement index was accuracy (Acc), as shown in Eq. (1):
$$\begin{array}{c}Acc =\frac{TP + TN}{TP + FN + TN + FP},\#\left(\text{1}\right)\end{array}$$
where TP represents true positive samples correctly predicted as MDD patients; TN represents true negative samples correctly predicted as healthy controls; FN represents false negative samples incorrectly predicted as healthy controls; FP represents false positive samples incorrectly predicted as MDD patients.
Analysis of the Importance of Brain Functional Regions
Understanding and determining the importance and impact of each feature value in predictive models is a challenging task36. SHAP analysis can be utilized to assess the significance of features37, Unlike black box models such as deep neural networks or SVM, this model encodes the relationship between input and output as a series of easily understandable if-then rules, providing a globally interpretable method based on optimal Shapley values to measure the contribution of each feature in the model38, and it demonstrates interpretable predictions for machine learning classifiers, addressing deficiencies caused by lack of feature directionality39. The feature selection method based on SHAP analysis determines the importance of brain functional regions carrying topological attribute features in distinguishing MDD from HC groups in a more interpretable manner. In experiments, functional region node features selected by LASSO were used as input features for an XGBoost model to determine the level of importance for distinguishing MDD from HC groups during LOSCV.