Predicting nodal response to neoadjuvant treatment in breast cancer with core biopsy biomarkers of tumor microenvironment using data mining
Neoadjuvant systemic treatment (NAST) has become a standard treatment strategy for patients with stage II-III breast cancer . It allows downstaging of the disease and assessment of response to therapy, which facilitates post-neoadjuvant treatment planning. Response to NAST is not only a favorable prognostic factor, but also allows for less aggressive surgical management of the axilla in this setting . For initially node-positive patients (cN+), axillary lymph node dissection (ALND) used to be the standard surgical procedure. However, in the absence of nodal metastases after NAST, sentinel lymph node biopsy (SLNB) with 3 or more nodes retrieved or other limited surgical options have been proven as feasible alternatives with acceptably low false-negative rates . Performing SLNB alone did not translate into worse outcome even after 10 years of follow-up . Pathologic complete response (pCR) in the axilla is expected in up to 70% of cases, but not all patients are successfully spared ALND . Although imaging has improved the preoperative nodal assessment compared to clinical assessment alone, its conventional use is still inferior to surgical staging . With an accuracy of up to 70%, axillary ultrasound may lead to an incorrect surgical strategy in about one third of cases . Nodal pCR is associated with the biological subtype of the tumor as well as other clinical and pathologic characteristics . In addition to tumor cell characteristics, the tumor microenvironment (TME) has gained attention in recent years to investigate both prognostic and predictive value in several cancers . Tumor-infiltrating lymphocytes (TIL) and the tumor stroma are major constituents of TME and allow for an overall stromal landscape assessment . While the prognostic and predictive value of TILs in breast cancer has been thoroughly studied , the content of tumor stroma has been shown to be not only a prognostic but also a predictive biomarker . The aim of the present study was to use data mining to create a model for predicting nodal response to NAST in cN+ patients that incorporates TME characteristics and could be used to optimize surgical planning. 

Data collection The study was approved by the national ethics committee (approval number 0120–178/2022/3). Due to the retrospective nature of the study, informed consent was waived. We retrospectively reviewed the records of patients diagnosed with breast cancer at our institution between 2008 and 2021. Female patients who were treated with NAST followed by surgery were eligible for the study. Patients with inflammatory breast cancer, primary metastatic breast cancer, bilateral breast cancer, synchronous malignancy at other site, history of invasive or in situ carcinoma of the breast, or breast cancer diagnosed during pregnancy were excluded. Only node-positive (cN+) patients were included in the study, which was defined as biopsy-proven nodal involvement at presentation. The following clinical and pathologic characteristics were recorded retrospectively: Age at diagnosis, menopausal status, tumor and nodal stage at diagnosis (according to TNM), pre-treatment core-needle biopsy (CB) tumor characteristics: ER %, PR %, Her-2 status, Ki-67%, grade, histological subtype, radiological response on breast imaging after NAST and axillary ultrasound findings after NAST. Subtype was determined based on CB characteristics. Axillary staging procedure We employ dual-tracer technique with a recommended sampling of at least 3 SLNs for patients that had biopsy-proven node-positive disease before neoadjuvant treatment, as previously described . Overall, 342/437 (78.3%) patients ultimately underwent ALND and the remaining 95/437 patients (21.7%) had sentinel lymph node biopsy (SNLB) only. Of these, 71 had no evidence of disease in the SLNs (neither ITC, micro- nor macrometastases). If the patient was rendered clinically and radiologically (using axillary US) node-negative after treatment, they were planned for SLNB. If there was a high suspicion of residual disease, they were proceeded directly to ALND. However, when less than 3 SLNs were retrieved, the patients would often undergo ALND directly. Core biopsy revision Because TILs and tumor stroma are not routinely reported on CB reports, all pre-treatment hematoxylin–eosin (H&E) slides were reviewed by an experienced breast pathologist. We had access to a single core biopsy in 30% of cases. For the remaining, multiple cores were available, and we followed the approach outlined by Le et al. This involved screening all available cores to identify the most stroma-rich area of the specimen. The selected cores needed to have invasive tumor cells present at both ends of the specimen, as well as along the side margins, or throughout the entire specimen. The tumor stroma ratio (TSR) was evaluated by tenfold percentage of stroma content. While 50% was established as the cut-off value for stroma-rich vs. stroma-poor tumors in resection samples, a cut-off value for CB specimen that classifies stroma-rich vs. stroma-poor tumors with optimal specificity and sensitivity (maximal Youden index) was set at 40% . TILs were evaluated on CB and reported as percentage values in accordance with previously published reports . For analysis, TILs were considered as three predefined groups of low (0–10%), intermediate (11–59%) and high TILs (≥ 60%) . The primary end-point was axillary pCR, defined as ypN0 in the final histopathology report. Statistical analysis and data mining Clinical and pathologic features were reported as absolute and relative frequencies. Features were compared between patients who achieved nodal pCR and those who did not using the Pearson chi-square test in IBM SPSS Statistics (IBM corp., Armonk, NY, USA). The Orange Datamining Toolbox (Laboratory for bioinformatics, University of Ljubljana, Slovenia) was used for model generation . The five most influential features were selected by backward elimination. Naïve Bayes, logistic regression, neural network, k-nearest neighbors and random forest algorithms were used to train the prediction model: Naive Bayes: a probabilistic model that is simple yet robust and offers valuable insights into effects of attributes on class probabilities. Its outputs can be visualized using nomograms for clear, intuitive support tools, for example in clinical decision-making process. Logistic Regression: we included this for its simplicity and interpretability. It is particularly useful for understanding the contribution of individual predictors in binary classification tasks like ours. k-Nearest neighbors (kNN): a non-parametric model that captures non-linear relationships and serves as a contrast to parametric methods. It predicts outcomes based on similarity by measuring distance between data points. Random forest: an ensemble method that builds multiple decision trees using random subsets of data and features, with predictions based on majority voting. While highly accurate, Random Forests can be less interpretable due to the complexity of multiple trees. They are useful in handling high-dimensional data and are robust against overfitting. Neural Networks: modeled after biological neural networks, they are highly effective for complex supervised learning tasks such as disease prediction. They consist of layers of interconnected nodes (neurons) that assign weights to inputs, adjusting through backpropagation. When multiple layers are used, this becomes "deep learning," which is particularly useful for complex modeling problems. Stratified tenfold cross-validation was used for model validation. The dataset was divided into ten random subsets (with preserving class distribution). The model was trained and evaluated ten times, each time using a different fold as the testing set and the remaining nine-folds as training set. Cross-validation helps in mitigating bias by using different subsets for training and testing in each iteration, ensuring that the model is exposed to various parts of the dataset. It also helps control variance by averaging the performance over multiple runs, making the evaluation more stable and less sensitive to the particularities of a single random split (e.g., manually splitting the data in a training and testing subset in advance). This is crucial for obtaining a robust estimate of the model's generalization performance across different data instances. Since the dataset is not extremely large, it also maximizes the use of all available data for both training and testing purposes, contributing to a more robust assessment of the model's capabilities. The area under the ROC curve (AUC), accuracy, sensitivity, specificity and F1 score were used to evaluate the performance of the models. The F1 score is a measure of accuracy that combines precision and recall. It is calculated as follows: F1 = ((precision x recall)/(precision + recall)) × 2. A nomogram was used for visualization . The study is reported according to REMARK guidelines . 

The study was approved by the national ethics committee (approval number 0120–178/2022/3). Due to the retrospective nature of the study, informed consent was waived. We retrospectively reviewed the records of patients diagnosed with breast cancer at our institution between 2008 and 2021. Female patients who were treated with NAST followed by surgery were eligible for the study. Patients with inflammatory breast cancer, primary metastatic breast cancer, bilateral breast cancer, synchronous malignancy at other site, history of invasive or in situ carcinoma of the breast, or breast cancer diagnosed during pregnancy were excluded. Only node-positive (cN+) patients were included in the study, which was defined as biopsy-proven nodal involvement at presentation. The following clinical and pathologic characteristics were recorded retrospectively: Age at diagnosis, menopausal status, tumor and nodal stage at diagnosis (according to TNM), pre-treatment core-needle biopsy (CB) tumor characteristics: ER %, PR %, Her-2 status, Ki-67%, grade, histological subtype, radiological response on breast imaging after NAST and axillary ultrasound findings after NAST. Subtype was determined based on CB characteristics. 

We employ dual-tracer technique with a recommended sampling of at least 3 SLNs for patients that had biopsy-proven node-positive disease before neoadjuvant treatment, as previously described . Overall, 342/437 (78.3%) patients ultimately underwent ALND and the remaining 95/437 patients (21.7%) had sentinel lymph node biopsy (SNLB) only. Of these, 71 had no evidence of disease in the SLNs (neither ITC, micro- nor macrometastases). If the patient was rendered clinically and radiologically (using axillary US) node-negative after treatment, they were planned for SLNB. If there was a high suspicion of residual disease, they were proceeded directly to ALND. However, when less than 3 SLNs were retrieved, the patients would often undergo ALND directly. 

Because TILs and tumor stroma are not routinely reported on CB reports, all pre-treatment hematoxylin–eosin (H&E) slides were reviewed by an experienced breast pathologist. We had access to a single core biopsy in 30% of cases. For the remaining, multiple cores were available, and we followed the approach outlined by Le et al. This involved screening all available cores to identify the most stroma-rich area of the specimen. The selected cores needed to have invasive tumor cells present at both ends of the specimen, as well as along the side margins, or throughout the entire specimen. The tumor stroma ratio (TSR) was evaluated by tenfold percentage of stroma content. While 50% was established as the cut-off value for stroma-rich vs. stroma-poor tumors in resection samples, a cut-off value for CB specimen that classifies stroma-rich vs. stroma-poor tumors with optimal specificity and sensitivity (maximal Youden index) was set at 40% . TILs were evaluated on CB and reported as percentage values in accordance with previously published reports . For analysis, TILs were considered as three predefined groups of low (0–10%), intermediate (11–59%) and high TILs (≥ 60%) . The primary end-point was axillary pCR, defined as ypN0 in the final histopathology report. 

Clinical and pathologic features were reported as absolute and relative frequencies. Features were compared between patients who achieved nodal pCR and those who did not using the Pearson chi-square test in IBM SPSS Statistics (IBM corp., Armonk, NY, USA). The Orange Datamining Toolbox (Laboratory for bioinformatics, University of Ljubljana, Slovenia) was used for model generation . The five most influential features were selected by backward elimination. Naïve Bayes, logistic regression, neural network, k-nearest neighbors and random forest algorithms were used to train the prediction model: Naive Bayes: a probabilistic model that is simple yet robust and offers valuable insights into effects of attributes on class probabilities. Its outputs can be visualized using nomograms for clear, intuitive support tools, for example in clinical decision-making process. Logistic Regression: we included this for its simplicity and interpretability. It is particularly useful for understanding the contribution of individual predictors in binary classification tasks like ours. k-Nearest neighbors (kNN): a non-parametric model that captures non-linear relationships and serves as a contrast to parametric methods. It predicts outcomes based on similarity by measuring distance between data points. Random forest: an ensemble method that builds multiple decision trees using random subsets of data and features, with predictions based on majority voting. While highly accurate, Random Forests can be less interpretable due to the complexity of multiple trees. They are useful in handling high-dimensional data and are robust against overfitting. Neural Networks: modeled after biological neural networks, they are highly effective for complex supervised learning tasks such as disease prediction. They consist of layers of interconnected nodes (neurons) that assign weights to inputs, adjusting through backpropagation. When multiple layers are used, this becomes "deep learning," which is particularly useful for complex modeling problems. Stratified tenfold cross-validation was used for model validation. The dataset was divided into ten random subsets (with preserving class distribution). The model was trained and evaluated ten times, each time using a different fold as the testing set and the remaining nine-folds as training set. Cross-validation helps in mitigating bias by using different subsets for training and testing in each iteration, ensuring that the model is exposed to various parts of the dataset. It also helps control variance by averaging the performance over multiple runs, making the evaluation more stable and less sensitive to the particularities of a single random split (e.g., manually splitting the data in a training and testing subset in advance). This is crucial for obtaining a robust estimate of the model's generalization performance across different data instances. Since the dataset is not extremely large, it also maximizes the use of all available data for both training and testing purposes, contributing to a more robust assessment of the model's capabilities. The area under the ROC curve (AUC), accuracy, sensitivity, specificity and F1 score were used to evaluate the performance of the models. The F1 score is a measure of accuracy that combines precision and recall. It is calculated as follows: F1 = ((precision x recall)/(precision + recall)) × 2. A nomogram was used for visualization . The study is reported according to REMARK guidelines . 

Overall, 437 cN+ eligible patients were identified. 151/437 (34.6%) achieved nodal pCR (ypN0). Patients’ clinical and pathologic characteristics are summarized in Table . Features with more than 25% missing values were excluded. After feature reduction, the following 5 variables were included in the model: ER, Her2, grade, stroma content and TILs. Stratified tenfold cross-validation was used for internal validation. The performance of the models for predicting pCR is evaluated in Table . The logistic regression model (with ridge regularization) showed the best performance: AUC 0.86, F1 0.72. The ROC curve and calibration curve are shown in Figs. and , respectively. A nomogram was created for visualization (Fig. ). 

Using machine learning methods, we developed a clinical tool to predict nodal pCR in biopsy-proven cN+ patients who underwent NAST. Our model achieved an AUC of 0.86. TILs, Her-2 status, stromal content, ER status and tumor grade were included in the model. It is worth noting that several prediction models have been previously proposed to estimate nodal pCR rates, taking into account different clinical and pathologic characteristics. Corsi et al. constructed a nomogram on a cohort of nearly 2000 patients; although their model yielded an AUC of 0.77, which is considerably lower compared to our tool, the AUC remained at 0.77 after validation on external cohort . Independent predictors included breast cancer subtype, clinical response in the breast, Ki-67 status and tumor grade. Similarly, in the prediction model developed by Guo et al. Her-2 status, ER expression level, response in the breast, and clinical stage at diagnosis were recognized as independent predictors of nodal pCR in hormone receptor positive breast cancer patients . Kim et al. also developed a nomogram on a cohort of 415 patients with a comparable pCR rate to our cohort . After validation on an independent testing cohort, the model yielded an AUC of 0.8. Of note, all patients underwent ALND to determine the final nodal status, which would no longer be acceptable in clinical practice. We developed the model using several commonly used algorithms and selected logistic regression as the model with the best performance characteristics. The high and comparable performance across different algorithms further strengthens the predictive value of the input parameters. In addition, Orange facilitates visualization of the logistic regression model through a nomogram, which is a valuable asset for the clinical setting . To reiterate, high TILs, Her-2-positive status, low stromal content, ER below 85% and high grade were predictors of nodal pCR in our model. It has been previously reported that stromal content of CB can predict response in the breast in Her-2-negative tumors by assessing the reduction of tumor cellularity . To our knowledge, this is the first study to report stroma content from pre-treatment CB of primary breast cancer as a predictor of nodal pCR in cytologically proven cN+ patients. In addition, our model includes TILs, which along with TSR provide a comprehensive assessment of TME . TILs serve as a surrogate biomarker of anti-tumor T-cell-mediated immunity. The extent of TIL infiltration within the TME has been shown to correlate with the tumor’s mutational burden . Notably, high TIL levels are associated with more aggressive clinicopathologic characteristics, such as higher cT stage and Ki-67 levels . While high TIL levels are more commonly found in TNBC and Her-2-positive subtypes, they have been associated with pCR across all breast cancer subtypes. This is consistent with our findings, suggesting a uniform response of immunological infiltrates to chemotherapeutic agents, regardless of breast cancer subtype. Furthermore, high TIL levels (60% or more, consistent with our reporting) have been associated with improved prognosis in TNBC and Her-2-positive subtypes. However, in luminal Her-2-negative cancers, high TILs have been identified as a negative prognostic marker for overall survival . This may be partially attributed to differences in the composition of immune cells within the infiltrate and it may indicate that, despite a favorable response to chemotherapy, other–possibly ER-related pathways–are contributing to the poor prognosis in this subgroup of luminal Her-2-negative cancers. Further research is needed to better understand both the role of peritumoral immune response and ER-related pathways in this subtype . Beyond inflammatory cell infiltrates, TME consists of extracellular matrix (ECM), cancer-associated fibroblasts (CAFs), and vascular endothelial cells. Since fibroblasts in the stroma have both tumor-promoting and tumor-suppressing action, the role of CAFs has not yet been fully comprehended . They are thought to promote angiogenesis, modulate the peritumoral immune response and restruct ECM components. They regulate cancer and other stromal cells, and potentially influence not only tumor growth and progression, but also modulate drug resistance mechanisms and immune evasion . High stromal content, characterized by an abundance of CAFs, has been previously associated with poor response to both neoadjuvant and adjuvant therapy in breast and colorectal cancers . Our findings are in line with the study of Li et al., who showed that a combination of TSR and TIL improves the prediction of ypT0/is ypN0 outcome (pCR in both breast and axilla), particularly in Her-2 negative tumors . While responses to NAST in the breast and axilla often overlap , our study focused specifically on predicting nodal response, which is key in planning the axillary surgical procedure. While it is not meant to replace the decision to opt for SLNB as initial surgical procedure in clinically negative axilla, an accurate nodal response prediction tool is beneficial in different borderline clinical scenarios. For example, when residual nodal disease is highly likely, SLNB procedure may be technically challenging with lower identification and higher false-negative rates . In such cases, proceeding to ALND when in doubt, reduces the risk of potential undertreatment. Conversely, if there is a high probability of nodal pCR, and we attempt a limited surgical staging procedure (SLNB or targeted axillary dissection (TAD)), but fail to identify a minimum of three lymph nodes or targeted node, we can opt for ALND omission. This is because we are already aware of the high probability of nodal pCR. The EUBREAST-01 trial aims to clarify the question of omitting axillary surgery in Her-2-positive and triple-negative patients who achieve a radiologic and pathologic complete response in the breast, but this is currently limited to cN0 patients. In cN+ patients, axillary surgical staging remains essential, as residual nodal disease reflects the chemoresistance of the disease and provides insights for further systemic treatment planning. The current prediction model may have to be adjusted in the future to better recognize subgroups of patients that despite not achieving nodal pCR will be candidates for ALND omission (e.g., OPBC-07 Study), but according to current guidelines, omission of ALND is of now only advised within clinical studies (Alliance A01120202, ADARNAT, OPBC-03 Trials). The potential advantage of the model is also to help guide adjuvant RT/systemic treatment decisions in these cases. Including radiologic assessment of response to NAST in the breast and axilla in the model would likely further improve the prediction tool. This has been previously noted by Kantor et al. who suggested that response in the breast is associated with response in the nodes, whether clinical, radiological or pathological . However, in our cohort of 437 patients, only 45% underwent MRI limiting the inclusion of MRI assessment in the model. Nonetheless, MRI and other imaging modalities open the door to radiomics, which is becoming increasingly popular in cancer imaging analysis and could also be used in the post-NAST setting to assess nodal response . The study has several limitations. The retrospective nature of the analysis harbors potential pitfalls. The stromal content from CBs was assessed by a single pathologist, reducing inter-observer variability but potentially introducing bias. The study spans over a relatively long period of time during which NAST regimens, as well as indications for NAST have changed. This most likely resulted in improved pCR rates over time. The trend may further continue (e.g., immunotherapy in the neoadjuvant setting for triple-negative cancers improves pCR rates). One fifth of the patients in our cohort had SLNB only as staging procedure, of which 75% were ypN0. While we aim to minimize the false-negative rates by dual-tracer technique and sampling of at least 3 SLNs, SLNB is not as accurate method of axillary staging as ALND. As a consequence, the proportion of nodal pCR in our study cohort may have been slightly overestimated. Nevertheless, it would not be acceptable to perform ALND for the sole purpose of axillary staging. Our data analysis was performed using machine learning techniques rather than traditional statistical approaches and lacks conventional null hypothesis testing. Although cross-validation ensures reproducibility of the model and improves generalizability, the most significant limitation of our study remains the lack of external validation .While our findings align with previous literature on predicting pCR, which enhances current study’s relevance and significance, we acknowledge that without validation of the model on an independent external cohort it is impossible to generalize our findings. Specifically, we aim to collaborate with other institutions within the EUBREAST Network to obtain external data sets, which will allow us to assess the robustness and applicability of our findings across different populations. We also encourage other centers to use similar tools on their own data to provide accurate risk estimation tools for individual patients and ultimately improve the outcome and cost-effectiveness of care . 

In summary, we have developed a clinical tool to predict nodal pCR for cN+ patients after NAST that includes biomarkers of TME and achieves an AUC of 0.86 after tenfold cross-validation. Similar tools are needed and we encourage the use of data mining techniques in clinical practice to improve outcomes. 

Below is the link to the electronic supplementary material. Supplementary file1 (DOCX 23 KB) Supplementary file2 (JPG 53 KB) 

