This study was approved by the Ethics Committee of West China Hospital, Sichuan University and registered on the Chinese Clinical Trial Registration Platform (ChiCTR2000033402). Drug-naïve, first-episode MDD adolescents and healthy controls aged 12–17 years from September 2020 to December 2021. All study subjects were from the inpatient or outpatient departments of West China Hospital of Sichuan University, and healthy controls were recruited through advertising. All participants and their guardians signed written informed consent. All adolescents with MDD were diagnosed according to the diagnostic criteria for MDD in the Diagnostic and Statistical Manual of Mental Disorders-Fourth Edition (DSM-IV). Two senior psychiatrists performed mental health examinations using the Chinese version of Affective Disorders and Schizophrenia-Present and Lifetime Version (KSADS-PL; Kaufman et al., 1997). Thirty-five adolescents with MDD (age 14.40 ± 1.24; 8 males and 27 females) and 24 healthy controls (HCs, age 14.21 ± 1.41; 11 males and 13 females) were included in this study. Inclusion criteria for both MDD patients and HCs were 12–17 years old, right-handed, have at least an elementary school education, have a normal brain structure, be able to understand the contents of the scales, IQ > 85, has not received electroconvulsive therapy, and has not taken any other medications recently. We exclude those subjects with other axis I and axis II mental illness, severe physical illnesses, history of neurologic disease or injury, previous use of illicit substances, pregnancy or breastfeeding. No HCs had psychiatric illness or had attempted suicide.

All adolescents with MDD were treated with antidepressants, either alone or in combination with antipsychotics, sedative-hypnotics or anxiolytics. Antidepressants included agomelatine (43%, 15/35, 50 mg/d), sertraline (37%, 13/35, 50–150 mg/d), escitalopram (14%, 5/35, 10–15 mg/d) and venlafaxine (6%, 2/35, 150–225 mg/d). Antipsychotics included quetiapine (14%, 5/35, 50–100 mg/d) and olanzapine (9%, 3/35, 2.5 mg/d). Sedative-hypnotics included alprazolam (26%, 9/35, 0.4 mg/d), zopiclone (9%, 3/35, 7.5 mg/d) and lorazepam (9%, 3/35, 0.5 mg/d). Anxiolytic was tandospirone (6%, 2/35, 30 mg/d). In subsequent analyses, the doses of antidepressants, antipsychotics and sedative-hypnotics were converted to equivalent doses of fluoxetine (Hayasaka et al., 2015), chlorpromazine (Andreasen et al., 2010) and diazepam (Borrelli et al., 2022), respectively.

Beck Depression Inventory (BDI) was used to evaluate the depressive symptom severity of patients before and after treatment (Jackson-Koku, 2016). Patients were defined as responder group (group A, BDI score reduction ≥50%), non-responder group (group B, 25% ≤ BDI score reduction <50% and group C, BDI score reduction <25%; Reeves et al., 2012). HCs also completed the BDI at baseline.

Resting-state structural and functional MRI data were collected in MDD adolescents before and after treatment and in HCs on the same 3.0 T MRI scanner (uMR790, United-Imaging Healthcare, Shanghai, China). Participants were asked to remain calm, close their eyes, relax, and not think about anything during the scan. Foam padding and earplugs were used to minimize the effects of head movements and scanner noise. After the scan, participants were asked if they fell asleep during the scan. During the scan, images were visually inspected for structural abnormalities, head motion and artifacts. Only images without the above influence were retained for the study.

T1-weighted (T1w) images were acquired using a magnetization-prepared, rapid gradient-echo sequence. The parameters were as follows: repetition time (TR), 8.4 ms; echo time (TE), 3.8 ms; flip angle, 8°; slice thickness, 0.8 mm; total number of sagittal slices, 208; matrix size, 256 × 256; field of view (FOV), 256 × 256 mm²; and voxel size, 0.8 × 0.8 × 0.8 mm³.

Resting-state background echo-planar imaging (EPI) were obtained. The parameters were as follows: TR, 1000 ms; TE, 30 ms; flip angle, 60°; FOV, 210 × 210 mm²; number of slices, 65; voxel size, 2.5 × 2.5 × 2.5 mm³; and layer thickness, 2.5 mm.

The original dicom file was converted to nifty file using MRIcron.¹ Then images of each subject were visually inspected again to ensure the absence of artifacts, deletions and ghosting before preprocessing. The data were preprocessed and analyzed using MATLAB R2022b (MathWorks, Natick, MA, United States), SPM12² and CONN toolbox version 22b.³ The steps for preprocessing the structural MRI data included translation to the (0, 0, 0) coordinates, segmentation and normalization (using gray/white/cerebrospinal fluid segmentation and MNI normalization). T1W image was used for segmentation because it has higher spatial resolution. Since we are ultimately interested in functional image data, we then had to resample these tissue maps onto the functional data. Functional pre-processing involved the following steps: realignment and unwarp (for motion and field map correction), translation of the image center (to the origin 0, 0, 0), slice-timing correction, outlier scan detection and scrubbing (using ART: artifact removal toolbox), spatial normalization to an MNI template (functional target resolution 2 mm) and functional smoothing (full width half maximum of 8 mm). Default preprocessing parameters were used to define and reduce possible confounding due to head movement and blood oxygen level-dependent signals in cerebrospinal fluid and white matter. For noise reduction due to physiological effects (such as respiration and pulsation) and drift caused by scanner noise, the signals for white matter and cerebrospinal fluid were regressed from the functional data, and a bandpass filter of 0.008–0.09 Hz was applied by default (Hallquist et al., 2013). After preprocessing, quality control was performed for each subject.

The ICA can separate signals from noise and increase the sensitivity of detecting inter-individual differences (Monakhova and Rutledge, 2020). Dyn-ICA matrices represent a measure of different modulatory circuits expression and rate of connectivity change between each pair of regions of interest (ROIs), characterized by the strength and sign of connectivity changes covarying with a given component/circuit timeseries. In this research, dynamic connectivity analysis was performed to explore cross-time functional modulation in the matrix of connectivity between ROIs (Calhoun et al., 2014). The dyn-ICA was conducted using the CONN toolbox and performed on the connectivity time-series, calculated the connectivity strength between any pair of ROIs at any given time point. The computation involved ICA decomposition of the ROI-ROI dynamic connectivity timeseries into 20 components with 30 smoothing kernels (selected by default in CONN). CONN’s dyn-ICA implementation follows Calhoun’s group-ICA methodology (Calhoun et al., 2001), with optional subject-level dimensionality reduction, concatenation across subjects, iterative dual regression on group-level data, a fastICA algorithm for group-level independent component definition (with a hyperbolic tangent contrast function), and subject-level back-projection (back-projection of the group-level connectivity matrices into a series of subject-specific connectivity components). Then, T tests (including independent and paired T tests) were used to compare different dFC patterns between MDD and HCs, as well as between before and after treatment. Threshold-free cluster enhancement was used for multiple comparison correction (Smith and Nichols, 2009). These ROIs were defined from the brain atlas file in the CONN toolbox and included 132 brain regions derived from the Harvard-Oxford Cortex Atlas (provided by the Harvard Center of Morphometric Analysis, Cambridge, Massachusetts) and cerebellum parcellations.

We extracted the dFC values between each two brain areas at different time points and calculated the average value of the entire scan time. Pearson correlation was used to investigate the association between mean dFC values (before and after treatment) and BDI score (BDI score reduction, BDI score before and after treatment) in MDD patients. Multivariate linear regression analyses were carried out to determine whether dFC patterns at baseline can predict antidepressant treatment response. Covariates included age, sex, baseline BDI score and medication dosage (including dose of antidepressant, antipsychotics, sedative-hypnotics and anxiolytics) of patients. To validate the prediction analyses, we conducted repeated 3-fold cross-validation (Ash and Hughes-Oliver, 2018) and permutation test. Cross-validation is virtually unbiased, but is known to be variable when applied to small samples (Li et al., 2022). To reduce its variability, the whole process was repeated 10 times. While permutation test was most used statistical analysis for small sample studies (Liu, 2024). We also conducted sensitivity and specificity analyses for each dFC pattern identified as a predictor to identify neuroimaging biomarkers of differential treatment response. An adjusted p value <0.05 [False Discovery Rate (FDR) correction] was used. All statistical analyses were performed using R 4.3.3.⁴

Demographic and clinical characteristics of the participants are listed in Table 1. There were no statistical differences in sex distribution ( x 2 = 2.471, p = 0.116) or age distribution (t = −0.551, p = 0.584) between HC and MDD group. In MDD patients, the BDI scores decreased significantly after treatment (t = 4.226, p < 0.001). 14% of patients achieved treatment response (Group A, 5/35), while 86% of patients failed to achieve response (Group B, 29%, 10/35 and Group C, 57%, 20/35).

In this study, seven brain regions (including right temporal occipital fusiform cortex, bilateral planum polar, bilateral planum temporal and bilateral Heschl’s gyrus) with different dFC patterns between MDD and HC group were found, which were listed in Supplementary Table S1. And several brain regions with different dFC patterns in MDD group between pre-and post-treatment were found (Table 2). These brain regions included right triangular part of inferior frontal gyrus (IFG tri r), bilateral insular cortex (right, IC r; left, IC l), bilateral central opercular cortex (right, CO r; left, CO l), right temporooccipital part of middle temporal gyrus (toMTG r), right supracalcarine cortex (SCC r), left posterior division of supramarginal gyrus (pSMG l), right cuneal cortex (Cuneal r), left frontal pole (FP l), right superior parietal lobule (SPL r; Figure 1).

Before treatment, we found the mean dFC value between IFG tri r and IC l (r = −0.419, p-FDR = 0.022) was negatively correlated with BDI score in MDD group (Figure 2A). There were no correlations between mean dFC values and BDI score after treatment.

The mean dFC value between IFG tri r and IC r (r = −0.461, p-FDR = 0.012) at baseline, as well as which between IFG tri r and IC l (r = −0.518, p-FDR = 0.007), were negatively correlated with BDI score reduction (Figures 2B,C). The mean dFC value between FP l and SPL r after treatment was positively correlated with BDI score reduction (r = 0.442, p-FDR = 0.014; Figure 2D).