Treatment Response Prediction in Major Depressive Disorder Using Multimodal MRI and Clinical Data: Secondary Analysis of a Randomized Clinical Trial

The study methods were preregistered prior to analysis at the Open Science Framework register (12). We performed a secondary analysis of data from the Establishing Moderators and Biosignatures of Antidepressant Response for Clinical Care for Depression (EMBARC) study, a multisite double-blind randomized clinical trial (ClinicalTrials.gov identifier, NCT01407094) funded by the National Institute of Mental Health. EMBARC was designed to identify pretreatment and early-treatment (week 1) clinical, neuroimaging, neurophysiological, and behavioral moderators for predicting response (at week 8) to sertraline and placebo in patients with MDD. Patients provided written informed consent after receiving a complete description of the study. Details of the trial design have been reported elsewhere (13).

The primary outcomes of EMBARC were based on symptom severity as measured with the 17-item Hamilton Depression Rating Scale (HAM-D). The trial consisted of two 8-week phases. In phase 1, all patients were randomized to treatment with either sertraline or placebo (see the supplemental methods section of the online supplement for the dosage regimen). In phase 2, nonresponders to placebo were switched to sertraline.

EMBARC included 296 outpatients 18–65 years of age who were diagnosed as having recurrent or chronic MDD and were not taking medication for depression; participants were not required to be medication naive. To test our hypotheses, we defined three population subgroups, for each of which we trained four sets of models. A flow chart of these subgroups is provided in Figure 1. Subgroup A consisted of all 146 patients randomized to the sertraline treatment arm, and subgroup B the 150 patients randomized to the placebo treatment arm in phase 1 of the study. Subgroup C comprised 77 patients from subgroup B who did not respond to placebo treatment in phase 1 of the study, after which they were switched to sertraline in phase 2.

Patients were excluded from the analyses if data were missing for them on the primary outcome for the second to last and last assessments (i.e., weeks 7 and 8 [phase 1] or weeks 12 and 16 [phase 2]), on more than one MRI sequence, or on pharmacological treatment for at least 2 weeks. We tested whether the exclusion of patients significantly affected the populations that were analyzed. We compared population characteristics before and after exclusion using Student’s t test for continuous variables and the chi-square test for categorical variables. The sociodemographic and clinical variables of the analyzed population are summarized in Table 1.

HAM-D symptom severity outcomes were collected at weeks 8 and 16. The first primary clinical outcome is response, defined as a reduction ≥50% in HAM-D score after 8 weeks. The second is remission, defined as a HAM-D score ≤7 after 8 weeks.

Data were acquired before treatment and 1 week after treatment initiation (“early treatment”). The data consisted of clinical assessments and MRI neuroimaging. Clinical assessments consisted of sociodemographic, behavioral, and neuropsychological assessments. See Table S1 in the online supplement for details. Table S2 in the online supplement provides details on the MRI scanner and pulse sequence information by site. In brief, four MRI sequences were acquired: T₁-weighted structural MRI, diffusion-weighted imaging, resting-state functional MRI (rs-fMRI), and perfusion MRI using arterial spin labeling (ASL). Diffusion-weighted imaging was acquired only at baseline. Although task-based fMRI was acquired, it was excluded from the present analysis because of the practical challenges it poses in clinical application.

All software programs used are listed in Table S3 in the online supplement. In short, T₁-weighted scans were standardized using FreeSurfer (14) and segmented based on the Desikan-Killiany atlas (15) with FastSurfer (16). Volumes per anatomical region were computed using PyRadiomics (17). FreeSurfer segmentations underwent systematic quality control (QC) using MRIQC (18). Diffusion-weighted imaging data were preprocessed and analyzed using the JHU ICBM DTI-81 white matter tract atlas (19) and tract-based spatial statistics (20). QC was performed automatically using EDDY QC (21). For rs-fMRI, standardization and QC were performed using the fMRIPrep data preprocessing pipeline (22). Subsequently, dual-regression independent component analysis was performed using FSL (23). The resulting independent components were masked by the corresponding well-matched maps from the BrainMap resting-state networks (24) at Z score >3. ASL coregistration, analysis, and QC were performed using ExploreASL (25). Mean quantitative cerebral blood flow was extracted for regions of interest defined by the Talairach atlas (26). In addition to the various QC tools employed, all segmentations underwent visual inspection.

For pretreatment prediction, pretreatment-acquired predictors were used, and for early-treatment prediction, we calculated relative reductions of predictors. Since MRI data were acquired only during phase 1, we did not perform early-treatment modeling when predicting phase 2 outcomes (i.e., week 16).

The clinical and MRI data contain many predictors. We limited the number of predictors using three predefined subsets (“tiers”) of predictors based on their level of scientific evidence. See Table S4 in the online supplement for a comprehensive list of these predictors.

Tier 1 predictors consisted of six MRI predictors for which level 1 scientific evidence exists, as indicated by meta-analyses and systematic reviews (summarized in reference 8). These predictors are the volumes of the left and right hippocampus, the mean and standard deviation of the frontal-limbic connectivity, and the mean and standard deviation of the quantitative cerebral blood flow in the anterior cingulate cortex (ACC). To these MRI predictors, the following clinical variables were added: age, sex, body mass index, chronicity, employment status, relative reduction in HAM-D score, and scores on the Mood and Anxiety Symptom Questionnaire, the Snaith-Hamilton Pleasure Scale, and the HAM-D. For tier 2 (see Table S4 in the online supplement), 54 predictors were added based on weaker (e.g., single study) evidence. Tier 3 consists of all 240 predictors without selection. For unimodal models, all predictors from other modalities were excluded, and the population was unaffected (Table 1, Figure 1). QC of produced predictors consisted of checking missing and outlier values.

We followed stringent recommendations to reduce the risk of bias and overfitting and implemented a nested cross-validation framework with repeats. Imputation of scan and outcome data, site harmonization, scaling of data values, oversampling of outcomes, and hyperparameter optimization were embedded in the inner nest and tested in the outer cross-validation nest.

Train and test partitions for our primary hypothesis and single-modality models consisted of folds from population subgroup A (sertraline-treated patients). For our second hypothesis on sertraline specificity compared with placebo, we externally validated models that were trained on subgroup A on subgroups B (placebo-treated patients in phase 1) and C (sertraline-treated placebo nonresponders in phase 2). For internal validation, we performed 10-fold, 10-repeat, randomized K-fold cross-validation (RSKF) to improve our performance estimation. For external validation, we also applied 10 repeats since repeats affect the validation subset used during hyperparameter optimization. In addition, we performed a post hoc analysis where we compared RSKF with repeated leave-site-out cross-validation.

Missing sequence-specific predictors were imputed using K-nearest neighbor imputation on data from the same missing sequence to prevent information leakage across MRI sequences. HAM-D outcomes were imputed using iterative multivariate conditional modeling on all available HAM-D scores. To mitigate confounding effects, we used an open-source harmonization method called ComBat (27). Covariates used in harmonization were age, age-squared, sex, brain volume, and site. See the methods section of the online supplement for more details. Response labels were balanced by random oversampling of the minority class without replacement.

We trained XGBoost (28) extreme gradient boosting machine learning models implemented in Python with the area under the receiver operating characteristic curve (AUROC) as the criterion and a logistic objective function. Gradient boosting models have been shown to outperform other standard machine learning methods in classification studies related to MDD (29). In a post hoc analysis, we compared the results of this XGBoost classifier with support vector machine and logistic regression classifiers. See the methods section of the online supplement for details on the implementation of these classifiers. Hyperparameters (see Table S5 in the online supplement) were optimized using Bayesian hyperparameter optimization implemented for 100 epochs (30). To handle the large range of number of predictors provided to the model, feature selection was performed using chi-square-based feature selection (31) with the number of features selected optimized as a hyperparameter. (See Table S3 in the online supplement for details on the software used.)

Additionally, we report the most important predictors using impurity-based feature importance. Importance provides a score indicating how useful or valuable each feature was in constructing the boosted decision trees. For more information on the calculation of predictor importance, see the methods section of the online supplement.

Primary predictive performance metrics are mean balanced accuracy (bAcc) and AUROC across repeats and folds. In addition, we report sensitivity and specificity. Balanced accuracy simplifies model evaluation independently of a priori outcome rates. We also computed the natural logarithm of the diagnostic odds ratio to allow comparison with a previous meta-analysis (10). For our primary hypothesis that our predictive models significantly outperform chance, we calculated the 95% confidence intervals for bAcc by chance, using one-tailed binomial tests using the number of unique test samples (32). Chance is defined as the a priori response rate, the best alternative method. As a supplementary (post hoc) analysis, we predicted continuous clinical outcome measures directly; see the methods section of the online supplement for more details.

To test our second hypothesis on treatment specificity, we tested whether external validation performance in subgroups B (placebo-treated) and C (sertraline-treated placebo nonresponders) was significantly lower than on internal validation in subgroup A (sertraline-treated) by performing a one-sided dependent t test between subgroup B and the independent test set A, and between subgroup C and test set A. The t test samples consisted of the mean performance estimates for all modeling configurations, that is, pretreatment and early-treatment prediction of remission and response. We excluded configurations that did not perform significantly better than chance from this analysis to avoid regression toward the no-information rate.

Finally, we tested whether multimodal models outperformed unimodal models. We tested the null hypothesis that their mean performances were identical. We used one-tailed sign tests between the unimodal and the multimodal AUROC score pairs and a p value of 0.05 as the significance threshold (33).

A total of 229 patients were included in the analysis; their mean age was 38.1 years (SD=13.3), and 65.9% were female (Table 1, Figure 1). Details on causes for missing data are provided in Table S6 in the online supplement. Also, excluding patients with missing primary outcomes, primarily owing to lack of follow-up, did not significantly affect population characteristics. Treatment adherence in the remaining cohort was generally high, since among patients included there were none who stopped medication for more than a week (see Table S6 in the online supplement).

An overview of mean cross-validation results is presented in Table 2. All tier 1 models achieved significantly better performance than chance, in addition to one tier 2–based model and two tier 3–based models. Additional metrics are provided in Table S7 in the online supplement. For the prediction of treatment response with sertraline, best model performance was achieved in the early-treatment prediction of response (bAcc=68% [SD=2.4], AUROC=0.73 [SD=0.03], sensitivity=70% [SD=3.3], specificity=70% [SD=3.3], ln[DOR]=2.46 [SD=1.10], respectively) (Figure 2, panel 1, and Table 2). This performance is significantly better than chance (p<0.001). On average, early-treatment performance (bAcc=66.5% [SD=4.2], AUROC=0.72 [SD=0.11]) was better than pretreatment performance (bAcc=63% [SD=4.2], AUROC=0.72 [SD=4.7]), but the difference was not significant (p=0.06). The performance was higher for models that used predictors with a higher level of scientific evidence; all the best models used tier 1 predictors. Results of the post hoc analyses using leave-site-out cross-validation (see Table S8 in the online supplement), the comparison of machine learning classifiers (see Table S9 in the online supplement), and prediction of continuous clinical outcomes (see Table S10 and Figure S1) are provided in the online supplement. The mean absolute error on tier 1–based regression models achieved a mean absolute error of 5 points on the HAM-D scale. Regression models performed worse on the binary classification task than the corresponding classifiers presented in the main work.

The most important pretreatment predictors were ASL predictors. The median importance of the mean and standard deviation of the quantitative cerebral blood flow in the ACC is 0.13. Note that importance sums to 1. The third most important predictor was the anhedonic depression score on the Mood and Anxiety Symptom Questionnaire (0.11). However, at early treatment, we saw a shift from ASL predictors to clinical predictors. The most important predictors were reduction in HAM-D score (0.17), HAM-D score at week 1, and anhedonic depression score on the Mood and Anxiety Symptom Questionnaire at baseline (0.12). All predictors mentioned were consistently selected through hyperparameter optimization. For more details, see Figure S2 and Table S11 in the online supplement.

For model specificity to treatment, performance was significantly reduced (p=0.01) when tested on placebo-treated patients in subgroup B (Table 2). However, half the models still achieved performance significantly better than chance compared with the primary analysis, being tier 1- and tier 2–based pretreatment and early-treatment prediction of remission and tier 2–based early-treatment prediction of response. An exception to this lower performance was the early-treatment prediction of remission (bAcc=64%, AUROC=0.66) in subgroup B using tier 3 predictors (Figure 2, panel 2), which outperformed the same model in subgroup A (bAcc=59% [SD=2.7], AUROC=0.62 [SD=0.05]).

When tested on the sertraline-treated placebo nonresponders of subgroup C, mean predictive performance was not lower than internal validation (p=0.06). The performance of predicting remission in subgroup C (bAcc=64%, AUROC=0.68 (Figure 2, panel 3) using tier 1 predictors was similar to that of internal validation in subgroup A (bAcc=64%, AUROC=0.68) (Figure 2 panel 1), and both tier 1 and 2 models performed significantly better than chance. However, performance on predicting response in subgroup C using tier 1 predictors (bAcc=56%, AUROC=0.56) was worse compared with subgroup A (bAcc=65% [SD=4.1], AUROC=0.71 [SD=0.06]).

The mean performance of multimodal models significantly outperformed all unimodal models except early-treatment and early-treatment clinical assessment-based prediction of remission (Figure 3). Statistical outcomes are listed in Table S12 in the online supplement. The 11 unimodal models that performed significantly better than chance consisted of six pretreatment models using ASL and five early-treatment models using clinical assessment data (Table 2).

Compared with other multimodal studies, our results confirm the positive conclusions of Patel et al. (34) and Leaver et al. (35). These studies’ small sample sizes of 44 and 19 patients, respectively, raise the possibility that their findings are affected by performance bias (36). More recently, Sajjadian et al. (11) presented an approach using data from the Canadian Biomarker Integration Network in Depression (CAN-BIND-1) study. Compared with that study, our results are an improvement both in terms of bAcc and AUROC on corresponding analyses (pretreatment prediction of response: bAcc, 65% vs. 60%; AUROC, 0.71 vs. 0.60; early-treatment prediction of response: bAcc, 68% vs. 66%; AUROC, 0.73 vs. 0.70). This improvement is larger for pretreatment prediction than for early-treatment prediction. Since CAN-BIND-1 follow-up clinical data were acquired at week 2 instead of week 1 as in EMBARC, we suspect that the CAN-BIND-1 clinical data might present a larger early-treatment effect. Since the two approaches are methodologically similar, another explanation for our improved results is the availability of ASL in EMBARC but not in CAN-BIND-1. Our unimodal results substantiate this hypothesis, indicating significant predictive power for ASL alone.

We performed three post hoc analyses: First, our findings are consistent with leave-site-out cross-validation (see Table S8 in the online supplement), suggesting that site effect is limited. Second, comparing alternative classifiers to the XGBoost classifier shows that these are also capable of predicting treatment response, albeit at lower performance levels (see Table S9 in the online supplement). Lastly, we found that regression models can predict continuous outcome scores (see Table S10 and Figure S1 in the online supplement).

When compared with the predictive performance reported by two recent meta-analyses of treatment response prediction based on MRI for pharmacotherapy, our approach outperformed the weighted estimate of the natural logarithm of the diagnostic odds ratio (2.46 [SD=1.10] and 2.11 [SD=0.56], respectively) (10), as well as the mean bAcc for adequate-quality studies (68% [SD=10] and 63% [SD=7], respectively) (36).

Several considerations should be borne in mind when qualitatively evaluating the utility of the model. First, it should be noted that with the accuracy measure used in this work (balanced accuracy), the actual accuracy will be higher when the response rate is greater than 50%. An AUROC >0.7 can be considered good, depending on the use. In our case, our model outperforms the current trial-and-error standard used in clinical practice. Thus, the benefits of utilizing our treatment planning support tool, which can predict treatment efficacy early on, might quickly outweigh the costs of MRI scanning. Additionally, our work outperforms studies of similar scope. However, the clinical effectiveness of treatment decision support tools should be determined in independent prospective randomized clinical trials.

In clinical practice, it might also be desirable to increase the detection rate of nonresponders at the cost of decreased sensitivity (indicating nonresponse in final responders). This would allow improved identification of true nonresponders to switch to another treatment early at a higher, acceptable rate of false “nonresponder” assignments among patients who might respond to a second selective serotonin reuptake inhibitor (SSRI) or serotonin-norepinephrine reuptake inhibitor (37). Finally, the lower bound at which performance predictive modeling using MRI balances burden and cost should be investigated.

The results for the sertraline-treated placebo nonresponders (subgroup C) were not significantly different from the sertraline-treated training set of subgroup A. In contrast, performance in the placebo-treated subgroup B did decrease significantly. This contrast suggests some specificity for predicting sertraline compared with placebo treatment response. Nevertheless, external validation in subgroup B achieved performance significantly better than chance, indicating some overlap. One explanation suggested by the literature is that, in addition to the pharmacotherapeutic effect, sertraline treatment is partially driven by a placebo response (38). Further studies are required to establish the specificity of our approach for other antidepressants.

Our results indicate consistently significantly higher performance in multimodal compared with unimodal models. This supports the claim that integrating neuroimaging data with clinical data can improve treatment response prediction. However, models of two modalities did perform significantly better than chance: pretreatment models using ASL only and early-treatment models using clinical assessment data only. These two data types are primary candidates for a lean treatment response model. As an exploratory analysis, we report on the performance of these models in Table S13 in the online supplement. Nevertheless, our findings corroborate the call for multimodal data integration to improve predictive performance (9).

Although differences in study design limit comparisons with other studies, the performance of unimodal models was similar to that of previous unimodal studies (this vs. other). This confirms results achieved by previous studies and is a sign that our multimodal results may be similar in other study populations. Bartlett et al. (39) predicted remission based on T₁-weighted MRI (bAcc, 53% vs. 51%, and AUROC, 0.55 vs. 0.59, respectively). Korgaonkar et al. (40) predicted remission using baseline diffusion-weighted imaging (bAcc, 52% vs. 54%). Finally, Chekroud et al. (41) used early-treatment clinical variables for predictions of remission (bAcc, 65% vs. 63%, and AUROC, 0.69 vs. 0.70, respectively).

Balancing prediction performance with interpretability is an ongoing challenge in neuroimaging-based modeling. However, given that even the best-performing neuroimaging-based models show modest prediction performance, more emphasis on interpretability may be warranted (42). Therefore, we performed a post hoc exploration of feature importance, as shown in Figure S2 and Table S11 in the online supplement. Early-treatment response and perfusion in the ACC consistently contributed to our model’s performance. ACC perfusion is a well-replicated biomarker for treatment response (8), and this finding aligns with previous work in MDD. A recent meta-analysis found that the ACC plays a crucial role in treatment outcome because of its role as a hub in the interplay between the ventral emotion-generating and dorsal cognitive-controlling systems. Here, SSRIs are thought to play a role in improving this ventral-dorsal control (43). After 1 week, we observed a high importance of symptom severity reduction, a known predictor (44). Unexpectedly, the hippocampal volume—the most replicated finding concerning treatment response in MDD (8)—did not consistently contribute. A possible explanation is that predictors were selected individually, while a strong effect of left lateralization of hippocampal volume in treatment prediction has been noted (45). Similarly, the connectivity in the default mode network measured with rs-fMRI did not contribute consistently to our models. For more information on predictor importance, see Figure S2 and Table S11 in the online supplement.

The design of the EMBARC trial allowed us to increase the generalizability of our results in three ways. First, EMBARC is the largest multimodal neuroimaging MDD data set available, and the sample size has been shown to reduce bias in performance estimation (9, 36). Second, data were acquired at four sites using 3-T scanners from different manufacturers. We preprocessed and harmonized the data using standard software packages to increase generalizability to new populations (46) and tested our results using leave-site-out cross-validation. Third, we performed validation in a separate study arm, instead of just internal cross-validation as in most studies (10). Performance was not significantly reduced in the separate set of sertraline-treated patients, which is a positive indicator of the generalizability of our method. An additional strength is that the hypotheses and methods in this work were preregistered prior to analysis. Preregistration avoids overestimation of performance due to overfitting, which happens when algorithms are tailored to training data after data analysis. Although preregistration of tier evidence (8) inhibits the exploration of other combinations of predictors, it increases the reproducibility and validity of our results since the prespecified predictors (tier 1) were backed by strong scientific evidence (systematic reviews or meta-analyses).

While the study design provides valuable strengths, several limitations should be acknowledged. First, our results are limited to sertraline treatment from data acquired in the EMBARC study. In the EMBARC study, information about response to previous antidepressant treatment was not available, which might be valuable to add to future predictive tools (47). Thus, further validation of our approach is required on even larger collections of external data, for different antidepressants, and in populations with more clinical heterogeneity than in a randomized clinical trial. Until then, our results should be interpreted with caution. To our knowledge, no data set exists that has collected the same multimodal data at the same time points as the EMBARC trial that would allow such external validation. Second, although excluding patients did not significantly affect population characteristics, excluding patients who did not complete the study’s first phase may have favored prediction results in the direction of adherent participants. Third, the EMBARC study’s design allowed us to validate our models on the same patients treated with placebo (subgroup B) and sertraline (subgroup C). However, subgroup C lacked patients from subgroup B who had responded to placebo treatment in the first phase of the study. This selection can be expected to incur a bias in our results. Since response is known to be partially driven by placebo response (38), lowered treatment response and decreased treatment response prediction performance in subgroup C could have been expected. However, the response rate in subgroup C was not negatively affected by this selection compared with subgroup A (Table 1). Our results also show that performance in subgroup C was not significantly reduced compared with subgroup A. Thus, we can assume that this type of selection bias did not affect our conclusions. Finally, although task-based fMRI might provide promising information to add to multimodal modeling (48–50) and likewise provide important insights into context-dependent neural mechanisms, we excluded this option, given the difficulty of their scalability and replicability in clinical practice (51). Therefore, we could not assess its potential benefits in treatment response prediction (8).

We found no significant difference in performance between pretreatment and early-treatment prediction. If externally validated, early-treatment response prediction would likely not require a second session of MRI scanning, reducing cost and lowering patient burden. Because our results show a performance drop in unimodal ASL models at early treatment but not at pretreatment, and we decided a priori to use relative changes at early treatment, which are sensitive to physiological variability, we suspect that this choice was suboptimal for ASL predictors. To overcome this limitation, we suggest combining absolute ASL predictors with early-treatment clinical predictors to improve performance. Other improvement options include integrating predictors from other sources and utilizing novel analytic methods, such as normative modeling (52).

In summary, our findings on a multimodal machine-learning-based method applied to data from the EMBARC trial show that pretreatment and early-treatment prediction of sertraline treatment response in MDD patients is feasible using brain MRI and clinical data and significantly outperforms chance and most unimodal models. Our results also suggest the specificity of our models for sertraline compared with placebo treatment. We found that ASL was the best unimodal predictor. With additional external validation, these findings will contribute toward the use of predictive modeling in individualizing clinical sertraline treatment of patients with MDD.