Development and Validation of Machine Learning Models in Predicting Prognosis of Breast Cancer Patients with Lymph Nodes Metastasis Following Neoadjuvant Chemotherapy

Introduction

Breast cancer (BC) remains the most prevalent malignancy among women globally, with significant heterogeneity in treatment responses. Neoadjuvant chemotherapy (NAC) is widely used in lymph node (LN)-positive BC patients to downstage tumors and eliminate metastatic LN deposits.^1–3 While 30% of BC patients achieve pathological complete response (pCR) post-NAC, a substantial subset exhibits resistance, underscoring the need for predictive biomarkers. This has spurred increased interest in predicting the effectiveness of neoadjuvant therapy. Previous studies have indicated that approximately 35–68% of patients with positive axillary lymph nodes (ALN) can achieve axillary pCR following NAC,⁴ leading to improved prognosis.^5–7 Efforts have been made to forecast the LN status in BC patients with LN metastasis post-NAC; however, further research is still needed.

Currently, magnetic resonance imaging (MRI) has demonstrated greater accuracy in predicting residual lesions after NAC compared to clinical examination, mammography or ultrasound.⁸ Previous studies have focused on predicting axillary pCR using radiomics or clinical variables alone.^9–14 Although deep learning algorithms have demonstrated efficacy in predicting individual oncological outcomes, few studies have integrated multimodal data to simultaneously predict LN response and achieve prognostic stratification of patient risk.^15,16 Our study addressed this gap through a multimodal framework integrating clinical parameters, radiomic features, and deep learning algorithms to predict post-NAC LN status and enable prognostic risk stratification.

Methods

Population Selection and Study Design

This study was complied with the ethical standards of Helsinki Declaration and was approved by the Ethics Committee of the First Affiliated Hospital of Chongqing Medical University (ID: 2025–286-01). Informed consent was waived by the Ethics Committee due to the retrospective nature of the study and the use of deidentified patient data. This study involved a retrospective collection of clinical and MRI data from breast cancer patients with ALN metastasis who underwent neoadjuvant chemotherapy at the First Hospital of Chongqing Medical University between January 1, 2019, and January 31, 2023. Treatment protocol and pathological evaluation were shown in Supplementary document 1.^17–19 All patients included in the study met the following criteria: 1) Diagnosis of breast cancer confirmed by pathological testing; 2) Completion of MRI prior to any invasive procedures or treatment initiation; 3) Presence of high suspicion of ALN metastasis either through imaging or confirmed by pathological examination; 4) Receipt of NAC; 5) Completion of the treatment procedure with comprehensive follow-up records.

The study retrospectively collected clinicopathological data and MRI images from patients. The MRI image processing steps are detailed in the Supplementary document 2. Patients were divided into two groups based on the response of axillary LN to NAC: those achieving pCR in the axillary LN and those who did not (Non-pCR). First, using three combination of dependent variables—Clinical, Clinical-Radiomics, and Clinical-Deep Learning Radiomics (Clinical-DLR). Machine learning models were trained and validated specifically to predict LN-pCR from the input data. This was its direct predictive endpoint. Then, patients were stratified into risk groups (low-risk vs high-risk) based purely on the model’s predicted LN-pCR probability thresholds. The performance of each model was evaluated and compared. Figure 1 provides a detailed illustration of the study methodology.

Figure 1 Flowchart of the study.

Follow-Up

All patients underwent interviews either through outpatient visits or telephone consultations from discharge until March 31st, 2024. Disease-free survival (DFS) was chosen as the metric for evaluating patient survival duration. DFS was defined as the duration from surgery until the earliest of the following events: 1) the initial recurrence of the tumor locally, regionally, or distantly; 2) the diagnosis of a secondary malignant tumor or contralateral breast cancer; or 3) death from any cause.

Construction and Validation of Machine Learning Models

The integration occurred within the logistic regression or the random forest models during the stepwise enhancement phase: 1) Clinical Model (Step 1): A logistic regression model and a random forest model were built using only selected clinical variables (eg, ER, HER2, T-stage). 2) Clinical-Radiomics Model (Step 2): A new logistic regression model and a random forest model were then built, combining the clinical variables mentioned above with the rad_scores (details were described in the Supplementary document 2). 3) Clinical-DLR Model (Final Step): Similarly, the final logistic regression model and random forest model were built by combining the clinical variables, rad_scores with the ResNet3D50-derived LN-pCR probability scores (treated as a continuous variable). These models were evaluated on stratified 70:30 training-testing splits. Performance metrics (sensitivity, specificity, AUC) were computed. For the Clinical models, clinicopathological variables were evaluated for association with post-NAC LN status by correlation analysis, as well as stepwise logistic regression. For the Clinical-Radiomics model, a radiomics score (rad_score) was calculated based on the selected radiomics features and incorporated into the Clinical models. For the Clinical-DLR model, the probability scores for LN status post-NAC predicted by deep learning model was incorporated in to the Clinical-Radiomics model.

Deep Learning Architecture

Based on the Darwin research platform (http://premium.darwin.yizhun-ai.com/), ResNet-3D-50 was used as the deep learning model. The ResNet-3D-50 model was adapted from the 2D ResNet50 architecture to process volumetric medical imaging data. The architecture comprised 3D convolutional layers, residual blocks with identity skip connections, and global average pooling. Input volumes (64×64×64 voxels, channel-wise normalized) first passed through a 7×7×7 convolutional layer with batch normalization (BN) and ReLU activation, followed by four residual stages containing [3, 4, 6, 3] blocks, respectively. Each block integrated three 3×3×3 convolutional layers with BN/ReLU and skip connections to preserve gradient flow. Spatial features were aggregated via global average pooling before final classification through a softmax-activated fully connected layer. Training used the Adam optimizer (lr=0.001) with categorical cross-entropy loss over 30 epochs (batch size=64, early stopping on validation loss plateau).

Statistical Analysis

Statistical analysis and data processing were conducted by RStudio (version 4.3.2) and Darwin research platform. For the clinical variables, continuous variables were assessed using either the Wilcoxon rank-sum test, Mann–Whitney test, or t-test, contingent upon the distribution’s normality. Meanwhile, categorical variables underwent analysis via the χ2 test or Fisher’s exact test. Two machine learning models—Random Forest and Logistic Regression—were developed for each of the three independent variable combinations: Clinical, Clinical-Radiomics, and Clinical-DLR. To evaluate the predictive performance of these models comprehensively, we employed Area Under the Curve (AUC), Accuracy (ACC), Sensitivity (SEN), Specificity (SPE), and F1-score. Subsequently, Kaplan–Meier analysis was conducted to estimate disease-free survival (DFS) curves, while the Log rank test was used to compare DFS between different groups, providing insights into survival trends associated with HER2 expression. All statistical tests were two-tailed, with a significance threshold of P<0.05 to determine statistical significance.

Results

Baseline of the Clinical Characteristics

In this study, a total of 447 patients were ultimately enrolled. The median age of the participants was 51 years, with an IQR of 45 to 57 years. Notably, more than half of the patients received TEC regimens as part of their NAC. Approximately two-thirds of the patients underwent between 6 to 8 cycles of NAC. Following the completion of NAC, a significant proportion of patients demonstrated favorable outcomes: 171 patients (38.3%), exhibited no residual tumor in the breast tissue. Additionally, 216 patients (48.3%) showed no residual tumor in the ALN. Furthermore, 138 patients (30.9%) achieved pCR. The median DFS duration for the cohort was 23.4 months, with an IQR ranging from 15.9 to 32.5 months. During the follow-up period, 60 patients experienced an event. For the purposes of analysis, the patients were divided into two groups: a training set consisting of 315 cases and a testing set comprising 132 cases, following a 7:3 ratio. The study flowchart was showed in Figure 1 and the relevant information was presented in Table 1.

Table 1 Baseline Characteristics of the Patients

Comparison of Patients with or without Residual Tumor in the Axilla

Compared to the patients with axillary-non-pCR, patients with axillary-pCR were more likely to have undergone more chemotherapy cycles (72.3% vs 60.6, p=0.039); a greater proportion of negative for hormone receptors (62.8% vs 33.8% and 77.5% vs 52.8% for ER and PR, respectively. p<0.001); more likely to be HER2-positive (51.9% vs 29.2%, p<0.001); higher Ki67 index (p=0.041); lower T-stage (in total of T1 and T2, 75.3% vs 68.1%, p=0.014); and more likely to be Breast-pCR post NAC (59.7% vs 15.3%, p<0.001). The relevant information could be found in Table 2. Additionally, survival analysis indicated that patients with axillary-pCR had a significantly better prognosis than those with axillary-non-pCR (p<0.001, Figure 2).

Table 2 Comparison of Characteristics Between Patients with Axillary-pCR and Those Without Axillary-pCR

Figure 2 Comparison of the survival outcome between the patients with axillary-pCR and those without axillary-pCR.

Model Construction and Validation

Clinical Models

Correlation analysis and stepwise logistic regression identified ER status, HER2 status, and T-stage as significantly associated with LN status post-NAC (Table 3). We first developed a logistic regression model using these variables. The model achieved a training set specificity of 0.684, sensitivity of 0.675, and AUC of 0.705, with testing set performance of 0.574, 0.789, and 0.738, respectively. Subsequently, a random forest model was trained using the same variables, with hyperparameters set to mtry = 1 and ntree = 500 (Supplementary Figure 1A and B). This yielded improved training (specificity: 0.665, sensitivity: 0.700, AUC: 0.714) and testing performance (specificity: 0.689, sensitivity: 0.718, AUC: 0.755). ROC curves are illustrated in Figure 3A and B.

Table 3 Stepwise Logistic Regression for Selecting Variables to Construct the Machine Learning Models

Figure 3 ROC curves of the Clinical models (A and B); Clinical-Radiomics models (C and D); Clinical-DLR models (E and F).

Clinical-Radiomics Models

From 1781 radiomic features extracted, six were selected for analysis after rigorous screening (Supplementary Figure 2): exponential_glcm_MCC, square_root_glcm_Imc1, log_sigma_3_0_mm_3D_first_order_TotalEnergy, lbp_3D_m2_glrlm_RunLengthNonUniformity, square_glcm_ClusterProminence, and wavelet_HHL_glcm_JointEnergy. These features were combined with the clinical variables to create a radiomics score (rad_score), which was integrated into the models. For the random forest model (mtry = 1, ntree = 500; Supplementary Figure 1C and D), the training set demonstrated specificity (0.852), sensitivity (0.744), and AUC (0.880), while the testing set showed specificity (0.590), sensitivity (0.887), and AUC (0.801). The logistic regression model achieved specificity (0.761), sensitivity (0.675), and AUC (0.776) in the training set, and specificity (0.852), sensitivity (0.606), and AUC (0.798) in the testing set. ROC curves are shown in Figure 3C and D.

Clinical-DLR Models

Deep learning-derived predictions were incorporated as independent variables to develop Clinical-DLR models. The random forest model (mtry = 1, ntree = 500; Supplementary Figure 1E and F) achieved exceptional performance: training set specificity (0.974), sensitivity (0.925), and AUC (0.988); testing set specificity (0.803), sensitivity (0.930), and AUC (0.923). The logistic regression model also performed well, with training set specificity (0.884), sensitivity (0.869), and AUC (0.928), and testing set specificity (0.820), sensitivity (0.887), and AUC (0.911). ROC curves are illustrated in Figure 3E and F.

Comparison of the Performance Among Each Model and Stratification of Patients

Model performance was evaluated using decision curve analysis (DCA, as illustrated in Figure 4) and summarized metrics (Table 4), with the Clinical-DLR model consistently outperforming alternatives. Subsequent risk stratification into high- and low-risk groups based on model predictions revealed significant prognostic differences via Kaplan–Meier analysis (Figure 5), independent of the modeling approach.

Table 4 Performance Matrix of the Machine Learning Models

Figure 4 The DCA curves of the logistic regression models (A) and the random forest models (B).

Figure 5 Stratification of the patients by the Clinical-DLR logistic regression model (A) and the Clinical-DLR random forest model (B) in the testing set.

Discussion

In this study, we investigated axillary response patterns in breast cancer patients with axillary LN metastasis following NAC. Our findings revealed that achieving pCR in the ALN strongly correlated with breast pCR (p < 0.001). Critically, patients with LN pCR exhibited significantly improved DFS compared to those with non-pCR counterparts, underscoring the prognostic value of LN response assessment. To translate these insights into clinical utility, we hypothesized that predictive modeling of LN response could refine survival risk stratification. Based on clinicopathological variables and pre-treatment MRI scans, three predictive frameworks were developed: Clinical, Clinical-Radiomics and Clinical-DLR. Furthermore, to avoid redundancy and maintain focus, we selected two representative algorithms—the simpler logistic regression and the more complex random forest—for stacked generalization. Our results confirm that integrating radiomics features (especially deep learning-based) significantly improves predictive accuracy. The Clinical-DLR model achieved superior performance, demonstrating that machine learning-driven assessment of LN response enhances prognostic accuracy. These results highlight the potential of multimodal predictive frameworks to guide personalized post-NAC management.

For years, the ability to predict pCR in LN early during treatment has long held significant clinical importance, as it facilitates personalized, evidence-based therapy adjustments through escalation or de-escalation. Li et al have underscored the potential of MRI to provide quantitative biomarkers for effectively predicting pCR in breast cancer patients.²⁰ Zhang et al have demonstrated that a radiomics nomogram, when combined with clinical factors, exhibited promising predictive efficacy for assessing the response of node-positive breast cancer patients undergoing NAC.²¹ In a different approach, Gu et al utilized DLR applied to ultrasonography to forecast tumor response and LN status following NAC.²² Similarly, Zhang et al develop an MRI-based deep learning signature for predicting axillary response after NAC in breast cancer patients.²³ While previous studies mark substantial progress, critical gaps persist: prior investigations have failed to validate prognostic outcomes, and conventional methodologies in this domain often overlook deep learning architectures that could optimize outcome prediction through multimodal data integration.²⁴

Our study highlights several critical considerations. First, we have identified 6 radiomics features, these features might have potential correlations with patients clinical manifestations and survival outcomes. For instance, in several studies, GLCM Cluster Prominence has been identified as a key feature in predictive models for survival outcomes.²⁵ During Clinical-Radiomics model development, direct integration of radiomics features initially provoked severe model overfitting, necessitating adoption of a conventional Rad-score framework. While Rad-scores improve reproducibility, their reliance on dimensionality reduction risks obscuring biologically meaningful tumor heterogeneity encoded in radiomics data.²⁶

Second, due to the unique characteristics of breast tumors, sequences such as T2-weighted imaging frequently exhibit indistinct tumor boundaries that blend with surrounding normal glandular tissue. This ambiguity presents substantial difficulties for target delineation. Our literature review confirmed that no widely accepted clinical guidelines or gold standards currently exist to address this issue. Consequently, we utilized only the T1-DCE sequence. Furthermore, rad_score integration alone failed to confer added predictive benefit, while excessive feature incorporation compromised generalizability by inducing overfitting.^21,27 Third, the random forest algorithm exhibited heightened susceptibility to overfitting when processing high-dimensional radiomics datasets, suggesting tree-based models may be suboptimal for such tasks. Instead, logistic regression or support vector machines (SVM) are recommended for their balance of interpretability and robustness.

To the best of our knowledge, our study confirmed that DCE-MRI images from breast cancer patients with LN metastasis harbor prognostically relevant data capable of effectively predicting clinical outcomes. Deep learning algorithms, particularly when applied to these datasets, demonstrate significant potential to decode complex prognostic biomarkers. Notably, while this work focused exclusively on T1-weighted DCE-MRI, it established a framework for guiding early therapeutic interventions in advanced breast cancer. The model has the potential to aid clinicians in optimizing patient treatment processes, enabling more personalized surgical and radiation pathways. Specifically, it could facilitate axillary surgery de-escalation, radiation therapy optimization, and systemic therapy de-escalation. Ultimately, non-invasive LN-pCR prediction offers a tool for risk-adapted axillary management, transcending the current dependence on invasive nodal staging and surgery-only response assessment. Notably, Grad-CAM visualizations indicated that the peritumoral region warrants focused attention (as revealed in Supplementary Figure 3).

However, there were several limitations. First, this study was a single-center study, which might have potential biases in patient selection and imaging protocols. Second, reliance on T1-DCE sequences alone precludes integration of complementary multi-parametric MRI data, which likely contributes to the suboptimal discriminative performance of the Clinical-Radiomics model.²⁸ Third, The association between these model-derived risk groups and long-term survival outcomes is an indirect assessment performed during the analysis phase. And most critically, the retrospective design necessitates validation in prospective clinical studies.

Although we have performed 5-fold cross-validation during model development to assess internal generalizability. The logistic regression model demonstrated robust discriminative performance, with AUC values of 0.955, 0.913, 0.907, 0.907, and 0.915 across validation cohorts (mean ± SD: 0.919 ± 0.020). And the random forest model was with AUC values of 0.952, 0.932, 0.938, 0.911, and 0.874 across validation cohorts (mean ± SD: 0.921 ± 0.030). Future research should make collaboration with multi-center trials to test our model prospectively. Or leveraging public datasets (eg, TCIA) with compatible clinical/imaging variables for benchmarking. Furthermore, tumor heterogeneity analyses could be integrated with multi-omics frameworks—including pathomics, transcriptomics, and genomics—to elucidate spatially resolved biological mechanisms driving metastasis and treatment resistance. We also advocate for integrating deep learning architectures, prioritizing algorithm-agnostic comparisons to identify optimal feature fusion strategies. Finally, it is essential to rigorously examine the stability of deep learning-derived features across heterogeneous cohorts and conduct reproducibility assessments to bridge the gap between computational outputs and clinically actionable insights.

In conclusion, while traditional radiomics and standalone clinical variables exhibited limited predictive utility, the integration of MRI-based deep learning algorithm, radiomics signatures, and clinicopathological risk factors yielded a robust model for forecasting axillary LN response to NAC. The proposed Clinical-DLR framework achieved superior discriminative accuracy in stratifying LN responders from non-responders. This could not only enhances prognostic precision but also enables dynamic risk stratification, thereby identifying high-risk patients who may benefit from early therapeutic escalation or adaptive treatment regimens.