This retrospective study was approved by the Ethics Committee of Sun Yat-sen Memorial Hospital, with a waiver for informed consent from all participants. In our institution, a total of 1,757 patients with primary breast invasive cancer were diagnosed between April 16, 2016, and April 30, 2020. The inclusion criteria were as follows: 1) an initial diagnosis of primary invasive breast cancer; 2) DCE-MRI performed before biopsy and within 1 week before NAC; 3) surgical excision of the tumor whether achieving pCR or non-pCR after NAC treatment. The exclusion criteria were distant metastasis (n = 150), another malignant tumor (n = 16), surgery but without NAC (n = 1,187), without any treatment (n = 33), non-standard NAC treatment (n = 8), or tumor progression during NAC (n = 7). The patient enrollment pathway is shown in the consort diagram (Figure 1). Finally, 356 patients were included for analysis. The entire cohorts were split into independent training and validation dataset by 5-fold cross-validation (18). Four-fold data (80% of the tumors) were used as training dataset, and the remaining one-fold data (20% of the tumors) were used as validation dataset. The prediction probabilities of five independent validation sets were collected as a whole set and used to evaluate the model performance. The 5-fold cross-validation procedure is illustrated in Supplementary E1 and Supplementary Figure S1.

Breast MRI was performed on a 1.5T unit (Magnetom Avanto; Siemens Medical Solutions, Erlangen, Germany) with patients in the head-first prone position. The body coil was used as the transmitter, and a dedicated 8-channel phased-array breast coil (Siemens Medical Solutions, Erlangen, Germany) was used as the receiver. MRI sequences consisted of axial T2-weighted turbo spin-echo (TSE) with short tau inversion recovery (STIR) sequence; axial T1-weighted volume interpolated body examination (T1W-VIBE) with Dixon sequence, and axial diffusion-weighted imaging (DWI) with spectral attenuated inversion recovery (SPAIR) fat saturation with 2 b values (b = 0, 800 s/mm²) and axial DCE imaging. DCE images were acquired by using a 3D fat-suppressed T1W-VIBE sequence. The DCE acquisition consisted of 40–70 measurements with a temporal resolution of 8 s and a total of 5–7 min of imaging time. After two consecutive measurements, gadodiamide (Gd-DTPA-BMA) (Omniscan; GE Healthcare, Ireland) was administered via intravenous bolus injection at a dosage of 0.1 mmol/kg and a flow rate of 3.5 ml/s, followed by a 20-ml saline flush. Before DCE acquisition, multiple flip angle images (2°, 4°, 6°, 8°, 10°, and 12°) were obtained for the calculation of T1 maps using the same sequence and parameters except for the flip angle. The details of acquisition parameters of MRI pulse sequence are provided in Supplementary Table S1.

The diagnosis of all patients was established by a core needle biopsy of the primary tumor before NAC. The regimens of NAC, provided in Supplementary E2, were defined according to the National Comprehensive Cancer Network (NCCN) guideline (19). According to the Food and Drug Administration criteria (20), all patients underwent surgical resection of the tumors and sentinel lymph node dissection (SLNB) or axillary lymph node dissection (ALND) after NAC. The resected tumors and lymph nodes were sampled for histologic examination to evaluate the chemotherapeutic response. The pCR (ypT0/Tis-ypN0) was defined as the absence of residual invasive tumor in the breast and axillary lymph nodes on the operative specimen (breast tumor and axillary lymph nodes) following NAC. In contrast, non-pCR was defined as a residual invasive cancer in the breast or axillary nodes.

DCE-MRI data were analyzed independently by two radiologists (ZC and CZ with 10 years and 8 years of experience with breast MRI) using specialized quantitative analysis software (Omni Kinetics, GE Healthcare). The kinetic parameters were calculated using the extended Tofts model. During measurement, the regions of interest (ROIs) were carefully drawn to cover the whole tumor. Necrotic or cystic areas of the lesions, if presented, were excluded from the evaluation. The intraclass correlation coefficient (ICC) of kinetic parameters between the two readers was 0.834–0.977. Data from the two readers were averaged for analysis. The least absolute shrinkage and selection operator (LASSO) regression analysis was applied to select independent predictive kinetic parameters. These selected kinetic parameters were used to construct the kinetic-only RA model using a robust supervised classifier, linear discriminant analysis (LDA) (21), which was employed to classify the NAC treatment efficiency by searching for a linear combination of the independent predictive kinetic parameters. A multilayer perceptron (MLP) neural network (22) was employed to construct the kinetic-only DL model. The structure of the MLP neural network is shown in Supplementary Figure 2A.

Molecular information, including the status of hormone receptor [estrogen receptor (ER), progesterone receptor (PR)], human epidermal growth factor receptor 2 (HER2), and Ki67 expression, was recorded from IHC results. ER/PR negative was defined as <1% of tumor cells with positive nuclear staining and ER/PR positive as ≥1% of tumor cells with positive nuclear staining; the cutoff for Ki67 was 14%; tumors with IHC staining of 0 or 1 were defined as HER2 negative, whereas tumors that either showed 3+ IHC staining or had gene copy number >2.0 were considered HER2 positive (23). The molecular-only LDA and MLP models were constructed by using the molecular information as input. The structure of the MLP neural network is shown in Supplementary Figure 2B.

For RA, the tumors were segmented on DCE-MRI images obtained 88 s after the beginning of the contrast agent injection, as the clinical breast DCE-MRI guideline indicates peak enhancement and obvious conspicuity at this time point in most breast cancers (24). Tumor segmentation was performed using ITK-SNAP software (https://www.itksnap.org) by one radiologist (ZC, with 10 years of experience in breast MRI) who was blinded to the clinical and histopathologic results. Tumors were segmented on a section-by-section basis until the whole tumor volume was captured and a three-dimensional ROI was acquired. A second radiologist (JS, with 21 years of experience in breast MRI) reviewed all the delineations to ensure correct segmentation. The segmented images were processed by using the open-source Python 3.7 (https://www.python.org.) and PyRadiomics toolkit to extract 851 radiomics features, including image intensity statistical, shape, texture, and wavelet features (Supplementary Table S2). A coarse-to-fine feature selection strategy was applied to reduce the dimension and avoid overfitting. Redundant features were removed according to the Spearman correlation coefficient, and then the optimal feature subsets (Supplementary Table S3) were selected using least absolute shrinkage and selection operator (LASSO) regression. The prediction models, based on optimal image features, were built by using the five machine learning classifiers [i.e., LDA, support vector machine (SVM), random forest (RF), AdaBoost, and Naive Bayes] to verify the performance of the classifiers to predict pCR successfully. Then, the optimal classifier was used to build the image-only and image-based RA model.

The integrative image-based RA model was further developed by incorporating kinetic parameters (image-kinetic RA model), molecular information (image-molecular RA model), or both (image-kinetic-molecular RA model) into the image-only model. The optimal feature subsets of integrative image-based RA models are shown in Supplementary Tables S4–S6. The workflow for building RA predictive models is shown in Figure 2. All the RA models were constructed by using Matlab R2018b (MathWorks, Natick, MA, USA).

For DL analysis, a rectangular box of 128 × 128 × 3 pixels in size was used to crop three consecutive slices showing the maximum cross-sectional area of the tumor as input. To ensure comparability of the image signal intensity across patients, image intensity was normalized to a fixed range of 0–1. Random rotation, flip, and translation were used for data augmentation to alleviate the possible overfitting in the training procedure of model development. The image features were extracted by using a deep residual neural network, ResNeXt50 (25), pretrained on a large-scale, well-annotated ImageNet dataset to automatically learn discriminative image features, as illustrated in Supplementary E3 and Supplementary Figure S3. The whole DL structure contained a ResNeXt50 CNN and three fully connected layers, with the probability of pCR as output to build the image-only CNN model. Adam optimizer was used to train all DL models with a learning rate of 0.0001 and a batch size of 32. The triplet loss procedure was introduced to extract more discriminative features using the output of ResNeXt50, and the cross-entropy was introduced as classification loss using the final output of the fully connected layer. Details of the loss function are provided in Supplementary E4.

The integrative image-based DL model was further developed by adding kinetic (image-kinetic DL model), molecular information (image-molecular DL model), or both (image-kinetic-molecular DL model) into the CNN of the image-only model. The kinetic and molecular information was incorporated in the first fully connected layer of DL models. The kinetic and molecular information was incorporated in the first fully connected layer of DL models. The framework for building DL predictive models is shown in Figure 3. All the DL programs were implemented in Pytorch (https://pytorch.org.) on an Intel Core i7-7700 K processor (Intel, Santa Clara, CA, USA) and Nvidia RTX 2080 Ti GPU with 11 GB RAM (Nvidia, Santa Clara, CA, USA).

Data were expressed as mean ± standard deviation for continuous variables, and categorical variables were summarized as frequencies and percentages. The differences in age, molecular information, histopathologic types, tumor number type, clinical T stage, clinical N stage, clinical TNM stage, and treatments between pCR and non-pCR groups were compared by χ² or Wilcoxon rank-sum tests as appropriate. The inter-rater agreement of kinetic parameter evaluation was assessed by using the ICC. An ICC value >0.75 indicates good to excellent agreement. The predictive performance of the models was assessed by the ROC curve analysis. The sensitivity, specificity, positive predictive value (PPV), negative predictive value (NPV), and accuracy of the models were calculated based on a cutoff value determined by the maximum Youden index. And their confidence intervals were calculated by bootstrap analysis with 10,000-fold resampling. DeLong method was used to compare the AUROC between the models. A two-sided P value <0.05 indicated statistical significance. All statistical analyses were performed by using SPSS software (version 21; SPSS, Chicago, IL, USA) and MedCalc software (version 18.9.1; MedCalc, Ostend, Belgium).

A total of 356 female patients (mean age, 46.9 ± 9.4 years) were included in this study. The clinicopathologic characteristics are shown in Table 1. Here, 83 patients (23.3%) achieved pCR (pCR group), while the remaining 273 patients (76.7%) were non-pCR (non-pCR group). pCR group had a higher prevalence of ER-negative, PR-negative, HER2-positive compared with the non-pCR group (all P < 0.001). There was no significant difference in age, Ki67, histological type, tumor number type, clinical T stage, clinical N stage, clinical TNM stage, breast surgery, and axillary surgery between the two groups (all P > 0.05).

The LDA was the most robust classifier across multiple classifiers (Supplementary Table S7). The image-only LDA model had 12 image features selected by LASSO regression (Supplementary Table S3). The image-only CNN models had 1,000 image features extracted by ResNeXt50. The K^trans, K_ep, and MaxSlope were the independent predictors and included in the kinetic-only LDA and MLP models. Their AUROC, sensitivity, specificity, PPV, NPV, accuracy, and corresponding 95% CI are shown in Table 2 and Figures 4A, B. The AUROC of the molecular-only LDA model was 0.744, which was higher than that of the kinetic-only LDA model (0.682, P = 0.012) and image-only LDA model (0.55, P < 0.001). The AUROC of the molecular-only MLP model was 0.752, which was higher than that of the kinetic-only MLP model (0.652, P = 0.007) and image-only CNN model (0.554, P < 0.001). The AUROC of the kinetic-only LDA model was 0.682, which was higher than that of the kinetic-only MLP model (AUROC = 0.652, P = 0.008). There was no significant difference between image-only LDA and image-only CNN models (AUROC = 0.55 and 0.554, P = 0.208), as well as between molecular-only LDA and molecular-only MLP models (AUROC = 0.744 and 0.752, P = 0.33).

The AUROC, sensitivity, specificity, PPV, NPV, accuracy, and corresponding 95% CI of integrative image-based RA and DL models are shown in Table 3 and Figures 5A, B. The AUROC of the image-kinetic-molecular RA model was 0.781, which was higher than that of the image-kinetic RA model (0.629, P < 0.001), while it did not differ from the image-molecular RA model (0.755, P = 0.118). The AUROC of the image-kinetic-molecular DL model was 0.832, which was higher than that of image-kinetic and image-molecular DL models (0.707, 0.79; both P < 0.001). The heatmaps (Figure 6) generated from ResNeXt50 based on the Grad-Cam algorithm (26) indicated that locations were crucial in generating the output.

The AUROC of image-kinetic, image-molecular, and image-kinetic-molecular DL model (0.707, 0.79, and 0.83, respectively) were significantly higher than that of the corresponding image-kinetic, image-molecular, and image-kinetic-molecular RA models (0.629, 0.755, and 0.781, respectively; all P < 0.001). The image-kinetic-molecular DL model had significantly higher AUROC than other integrative models (Table 3).