Biologically-informed integration of drug representations for breast cancer treatment using deep learning

Domain: neuroscience
Published: 2025-12-11
Editorial outcome: major_revision

Abstract
Approximately 50% of breast cancer patients receiving neoadjuvant therapy do not achieve pathological complete response (pCR). Identifying personalized optimal therapeutic strategy is an unmet major challenge. Here, based on 4371 eligible patients from 31 datasets, we develop GDnet, an interpretable deep learning model integrating drug representations and tumor transcriptome profiles to predict the response to neoadjuvant therapies and aid in selecting optimal treatment strategy. We demonstrate that GDnet significantly outperforms transcriptome-only model in predicting treatment response. Then we conduct two series of in-silico simulated clinical trials based on I-SPY2 trial and test datasets respectively and show that GDnet can significantly increase the pCR rate in all trials. Moreover, the odds ratios for pCR increase from 1.6 to 2.5 linearly as optimization intensifies. Overall, GDnet can function as a digital drug-testing surrogate to optimize treatment decision-making. This approach may have broader applications across various treatment settings and cancer types.

Introduction
Breast cancer treatment currently faces many challenges. Approximately 50% of patients receiving neoadjuvant therapy do not achieve pathological complete response (pCR)1,2, and approximately 21.4% patients will experience local recurrence and 38.2% will experience distant recurrence within 15 years3. However, current methods for predicting the efficacy of neoadjuvant therapy in patients with breast cancer are inadequate. Previous models based on genomics4, transcriptomics5, radiomics6, pathology7, or even multiomics8 have been developed based on small sample sizes and cohorts with simplex treatment strategy. However, diverse treatment options, which are critical for cancer response, were not integrated into these models. Although these models can achieve high predictive performance in the treatment scenarios similar to the ones with which the models were built, they cannot be generalized to more complex real-world scenarios well and directly guide treatment regimen selection for patients. Unsurprisingly, physicians continue to select neoadjuvant treatment regimens only using empirical clinical factors such as molecular subtype and tumor stage.

Theoretically, accurately predicting the efficacy of cancer treatment requires addressing two aspects. First, a more comprehensive characterization of the tumor hallmarks is essential to better assess its functional status9. This task requires incorporating as many omics data types as possible, including transcriptomic, genomic, proteomic, and pathological information, to achieve a more comprehensive characterization of tumors and improve the predictive performance of models8,10. However, as the number of features increases, the sample size required for training also grows significantly11. Currently, aside from transcriptomic data, there is limited availability of patients with biologically multi-omics information, and existing large sample-based models are primarily based on transcriptomic data12. Therefore, a comprehensive characterization of transcriptomic status is crucial. Deep learning methods allow for the use of as many genes as possible to optimally characterize the tumor’s state due to the flexibility of neural networks, capturing not only tumor heterogeneity13 but also the tumor microenvironment14,15 information. However, deep learning methods require a large sample size, and their results are more difficult to interpret. Integrating prior biological knowledge can facilitate the construction of sparse, high-performance, and interpretable deep learning models16,17.

The second aspect involves appropriately characterizing drugs to accurately simulate the perturbations in tumors induced by treatments, which is essential for predicting the tumor’s response to treatment. To develop a model that can better guide clinical drug decision-making in precision medicine, incorporating as many diverse treatment types as possible is necessary to learn complex drug‒tumor interaction patterns. To incorporate therapeutic drugs into the model, three important factors need to be considered: (1) The types of anticancer drugs used in clinical practice are relatively limited (fewer than one hundred) and structurally diverse (e.g., small molecules, antibodies). Given the insufficient number of breast cancer therapeutic drug types to learn compound structural characterizations and the difficulty of integrating compound structural data with transcriptomic states, previous strategies that directly characterize compounds on the basis of their structure18 are not suitable for our task. (2) Many chemotherapeutic drugs lack specific targets, and the mechanisms of action for targeted therapies often extend beyond the primary target (e.g., antibody-dependent cell-mediated cytotoxicity related genes in anti-HER2 therapies). Thus, previous characterization of drugs based on specific gene targets cannot fully represent all types of drugs or comprehensively characterize their mechanisms19–21. (3) Clinically, drugs are often used in combination, leading to complex synergistic or antagonistic interactions that require flexible model structures (such as deep learning models) to simulate these interactions. Therefore, existing frameworks based solely on late fusion ensembles cannot achieve optimal efficacy predictions22,23. In general, A unified, genome-wide approach that can be applied to characterize different types of drugs and account for the complex interactions in combination therapies is needed.

In this work, we assemble a multi-dataset cohort of breast cancer transcriptomic data paired with diverse treatment records. Through systematic exploration of drug and transcriptome integration strategies to simulate drug-tumor interactions, we develop GDnet, a deep learning model for predicting neoadjuvant treatment response. GDnet demonstrates robust performance across various treatment regimens and breast cancer subtypes, serving as a digital drug-testing surrogate to optimize clinical treatment decision-making. This study will significantly advance precision neoadjuvant treatment selection and improve treatment efficacy in breast cancer.

Results
Schematic of GD network flows
To sufficiently fuse drugs with transcriptomic data, mapping drugs as a series of gene–value pairs is one of the most suitable approaches. The selection of gene–value pairs is the key consideration. The GDSC and CTRP datasets, which contain the transcriptome status of more than 800 cell lines of various cancer types (pan-cancer) and their drug sensitivities to more than 500 molecular drugs, are great reference datasets. The heterogeneous responses of different pan-cancer cell lines in GDSC or CTRP to drugs can be attributed to the heterogeneity of tumor transcriptome states; thus, genes that are significantly correlated with drug sensitivity are more likely to be response- or resistance-related genes and therefore more likely to be genes that are perturbed when tumors are exposed to a specific drug. Therefore, we used the correlation of genome-wide expression with each drug’s half maximal inhibitory concentration (IC50) value (where a lower value indicates greater sensitivity) across different cell lines as the gene–value pairs to represent each drug (Fig. 1a). Notably, antibody or peptide drugs such as anti-HER2 drugs are not included in these drug sensitivity datasets, and IC50 values are not traditionally applied to measure their sensitivity due to complex mechanisms of action24–26. Given that the sensitivity of these targeted drugs is strongly positively correlated with targeted gene expression, we used the negative value of the genome-wide expression correlation with drug target expression as the gene–value pairs to represent antibody drug (Fig. 1a). For regimens with multiple drugs, a gene-wise summation of each drug representation in the regimen was performed (Fig. 1b).

To better simulate the perturbations in tumors induced by treatments, we input rank-normalized transcriptome and regimen representations in parallel into the gene level layer of the subsequent neural network, which we refer to as GDnet (Gene-Drug net) (Fig. 1b). We further explore three additional model frameworks: (1) GAnet (Gene-Attention net), which integrate drugs and transcriptome utilizing a method imitating key-value attention mechanism (Supplementary Fig. 1a); (2) Gnetens (Gene net ensemble), which employs an ensemble learning method to aggregate predictions across all drugs administered to a patient (Supplementary Fig. 1b); and (3) Gnet (Gene net), which relies solely on transcriptome information (Supplementary Fig. 1c). Next, the gene layer is fed into the biologically informed neural network proposed by Elmarakeby et al.16, which sets nodes, layers, and connections by incorporating the affiliation (child‒parent) relationships between genes, fine pathways, and more complex pathways on the basis of the Reactome biological knowledge dataset (version 2022) (Fig. 1b). In other words, the flow of information from the input to hidden layers are guided by prior biological knowledge. The finally built prediction model is designed to precisely predict the response of breast cancer patients to neoadjuvant therapy, aiding in the selection of regimens among multiple candidates.

Drug representations based on the GDSC and CTRP datasets
Applying the drug representation strategies described in Fig. 1a, we successfully generated the gene-correlation pairs for all drugs used in the included breast cancer neoadjuvant therapy datasets. The representations of all drugs were aggregated (Fig. 2 and Supplementary Data 1). Methotrexate has the largest number of represented gene-value pairs with 5378 pairs in total. Among 22 drug representations, 17 have larger than 1000 gene-value pairs, 5 (ANG, lapatinib, carboplatin, epothilone B, and doxorubicin) have less than 1000 gene-value pairs. Consensus clustering of drugs clustered these 5 drugs with fewer number of pairs together (Fig. 2). For drugs with larger number of pairs, three main clusters were generated, the first cluster contains anti angiogenic therapies, anti HER2 therapies, and anti IGF1R therapies, the second cluster contains chemotherapies, AKT inhibitors, and PARP inhibitors, and the last cluster contains anti PD-1 therapies and HSP90 inhibitors (Fig. 2). Then the gene-correlation pairs of drug representations and gene-expression pairs of the tumor transcriptome were fed into the gene layer of GDnet network framework described in Fig. 1b.

GDnet showed superiority in predicting the outcome of neoadjuvant therapy
We collected 9550 breast cancer samples from 68 independent public datasets. After removing duplicates and filtering out samples with non-pretreatment status or missing key information (outcome and treatment), 31 datasets with 4371 samples and 16,133 genes were filtered out (Supplementary Fig. 2a). All breast cancer subtypes in populations from Northern America, Europe, and East Aisa were included (Supplementary Fig. 3a, b). Seventeen datasets containing the largest number of cases, totaling 3459 samples, were used to train the model, the remaining 14 datasets with the least number of cases, totaling 912 samples, were employed as external test sets, and 2 datasets with 134 patients published after February 15, 2024 were employed as prospective test sets (Supplementary Data 2).

For benchmark comparison, we trained the biologically informed models GDnet, GAnet, Gnetens, and Gnet (informed group), where Gnet was the baseline model with only transcriptome data; the dense neural network (dnet) variants (without biologically informed knowledge designed connections) of GDnet, GAnet, Gnetens, and Gnet, named GDdnet, GAdnet, Gdnetens, and Gdnet (dense group); dense deep neural network variants with hyperparameter optimized nodes and layers (without designed nodes, layers, and connections), named GDdnet-o, GAdnet-o, Gdnetens-o, and Gdnet-o (optimized group); deep learning models using previously applied drug representation strategies with optimized hyperparameters, named GSnet (Gene-SMILES net) and GTnet (Gene Target net) (previous group); and machine learning variants of GDnet with optimized hyperparameters (machine learning group). GDdnet, Gnet, GTnet, GSnet, and the drug-only model Dnet (Drugnet) were used for ablation study for GDnet. Tenfold cross-validation was applied for hyperparameter tuning and internal validation using the training set. The results revealed that models with multiple optimized hyperparameters had higher AUROC values in tenfold cross validation than informed and dense group models (Supplementary Fig. 4).

The deep learning models were then retrained with the whole training set 10 times and tested with the external test set, facilitating a robust comparison across different models. In the external test, GDnet significantly performed better (median AUROC=0.712, median AUPRC=0.691) than all the other models did (e.g., for the transcriptome-only model Gnet, the median AUROC=0.676; median AUPRC=0.632) (Fig. 3a; Supplementary Fig. 5). Notably, GDdnet and GDdnet-o, which have the largest number of parameters (168,977,476 and 200,048,640), exhibited worse performance (AUROC=0.678, AUPRC=0.661) than GDnet, which has only 37,127 parameters, demonstrating the importance of prior knowledge in guiding training and increasing performance. The ablation study further showed that deleting either component of GDnet decreased performance (Supplementary Fig. 6).

To generate the final robust deep learning model for future application, we performed ensemble learning by averaging the predictions from the 10-times retrained models for the final prediction. Similarly, GDnet (AUROC=0.725) achieved the highest AUROC (Fig. 3b), and its ROC curve significantly differed from that of transcriptome model Gnet (AUROC=0.683) (Delong’s Z=2.734, P=0.006, ΔAUROC=0.042, 95% CI [0.012, 0.073]) (Fig. 3c). The AUPRC results were similar to the AUROC results (Supplementary Fig. 7a). GDnet also showed generally high performance in terms of accuracy, F1-score, precision, recall, and specificity when cutoff value was set as 0.5 (Supplementary Fig. 7d–f). Moreover, GDnet demonstrated robust performance across all breast cancer subtypes (Fig. 3d; Supplementary Fig. 8) and outperformed Gnet in most cases. The only exception was TNBC, where GDnet showed slightly worse performance (Fig. 3d; Supplementary Fig. 8). We then tested the performance of GDnet in 13 individual datasets in the test datasets via the 10-times repetitions ensemble models to determine whether GDnet can perform well in limited regimen settings compared with Gnet. The GSE191127 dataset was not included in individual testing because it contains only residual disease (RD) patients. The results showed that GDnet outperformed Gnet, achieving higher AUROCs in 10 out of 13 datasets, while exhibiting only marginally lower AUROCs (difference<0.02) in the remaining three datasets (Supplementary Fig. 9a). Similar results were observed for the AUPRC values (Supplementary Fig. 9b). The comparison among all models is displayed in Supplementary Fig. 10. These results suggested that GDnet successfully captures the patterns reflecting how treatments perturbate tumor transcriptomes rather than simply overfitting the drug distribution bias among datasets. To further support the reliability of the model, we tested the model performance in the prospective test sets. The comparison results also showed that GDnet had higher performance than most other models (Fig. 3e; Supplementary Fig. 11), with a significantly better ROC curve than Gnet (Delong’s Z=2.76, P=0.006, ΔAUROC=0.059, 95% CI [0.017, 0.1]) (Fig. 3f).

GDnet can optimize the regimen selection for better outcomes
To assess the ability of GDnet to guide regimen selection for clinical use, we then conducted in-silico clinical trials to inspect if GDnet can improve breast cancer neoadjuvant treatment outcome (Fig. 4a). Given the consideration of multiple potential regimens for a patient, we input tumor transcriptome data and drug representation of each regimen into GDnet to computationally predict the pathological complete response (pCR) probability for each regimen. We selected the best regimens with the top N% highest predicted pCR probabilities, which are likely to benefit the patient, and the worst regimens with the top N% lowest predicted pCR probabilities, which are not suitable for patients. If the regimen actually used by the patient was included among the best regimens, we assigned the patient to the optimized group, and if it was included among the worst regimens, we assigned the patient to the control group. To inspect the baseline characteristics of the two groups, we also calculated the predictive score via the transcriptome-only model to evaluate the degree of tumor malignancy between them, and Gnet is priorly considered because its similar structure with GDnet and its more robust performance than GNNet regardless of AUROC (Levene’s W(1, 18)=5.9, P=0.026) or AUPRC (Levene’s W(1, 18)=6.9, P=0.017) (Fig. 3a, b).

We first conducted a series of in-silico trials based on the I-SPY2 trial, one of the datasets used to train our model. A total of 653 patients who received one of 11 regimens in I-SPY2 were enrolled in our simulation study, where these regimens were set as optional choices. We set varying optimization thresholds, starting at 50% and decreasing in 5% increments (e.g., 50%, 45%, 40%, 35%) until either group size dropped below 30 (Supplementary Fig. 12a), to categorize the best and worst regimens, thereby ensuring the robustness of the results. Lower thresholds indicate increased optimization stringency. The results revealed that the pCR rates (Supplementary Fig. 12c) were significantly greater in the optimized group than in the control group in all optimize thresholds. Moreover, the odds ratio (OR) for pCR increased linearly (r(6)=−0.849, P=0.008, 95% CI [−0.972, −0.358]) as the optimization intensifies (optimize threshold decrease) (Supplementary Fig. 12d). However, a higher Gnet prediction score was observed in the optimized group (Supplementary Fig. 12b), suggesting that tumors in this group were less malignant. This difference may bias the effect of our model. Thus, we performed Propensity Score Matching (PSM) to balance transcriptome profile based on Gnet prediction score (Supplementary Fig. 12e, Supplementary Fig. 12f). Similarly, pCR rates were significantly higher in the optimized group than in the control group (Fig. 4b), with OR increased linearly (r(6)=−0.932, P<0.001, 95% CI [−0.988, −0.661]) as the optimization intensifies (Fig. 4c), suggesting that GDnet has the ability to optimize treatment regimen choices.

We further performed a series of simulation trials with the external test datasets. Given the limited variety of regimens in each individual test dataset, we combined all available test datasets together, yielding 912 patients and 19 regimens. Similar to what was observed in the above I-SPY2 simulation, the results for the external test set showed that the pCR rates (Supplementary Fig. 13a, Supplementary Fig. 13c) were significantly greater in the optimized group than in the control group in all optimize thresholds, with OR increased linearly (r(6)=−0.824, P=0.003, 95% CI [−0.957, −0.404]) as the optimization intensifies (Supplementary Fig. 13d), even a lower Gnet prediction score was observed in the optimized group (Supplementary Fig. 13b). PSM analysis based on Gnet prediction score (Supplementary Fig. 13e and Supplementary Fig. 13f) also revealed that pCR rates were significantly higher in the optimized group than in the control group (Fig. 4d), with OR increased linearly (by PSM, r(6)=−0.971, P<0.001, 95% CI [−0.993, −0.879]) as the optimization intensifies (Fig. 4e), further indicating the excellent ability of GDnet to guide neoadjuvant treatment selection for breast cancer.

Interpretation of GDnet reveals important regimens and genes involved in the treatment response
We then sought to investigate which treatment types contribute to response importantly. With respect to the treatment optimization strategy shown in Fig. 4a, we further inspected the distribution of regimens and drugs used in each group in PSM analysis at different optimization thresholds; the differences in the proportions of regimens and drugs in different groups can reveal those with high contributions to the model. The results revealed that the paclitaxel|doxorubicin|cyclophosphamide|olaparib|PDCD1 and paclitaxel|doxorubicin|cyclophosphamide|ERBB2|ERBB2 regimens were used much more frequently in the optimized group than in the control group, with an increasing tendency as the optimization intensifies (Fig. 5a). Docetaxel|cyclophosphamide, doxorubicin, and 5-fluorouracil|epirubicin|cyclophosphamide|lapatinib were regimens that were mostly or even exclusively used in the control group, with an increasing tendency as the optimization intensifies (Fig. 5a). Specifically, drugs such as PDCD1 (PD-1 inhibitors), olaparib, and ERBB2 (anti-HER2 antibodies) may play important roles (Supplementary Fig. 14) in improving the pCR rate.

To understand the important genes and pathways that influence the response to the promising regimens, paclitaxel|doxorubicin|cyclophosphamide|olaparib|PDCD1 (indicated as the immunotherapy regimen below) and paclitaxel|doxorubicin|cyclophosphamide|ERBB2|ERBB2 (indicated as the dual anti-HER2 regimen below), in the optimized group, we visualized the layers of GDnet based on the patients who applying the corresponding regimen in the external test dataset (Fig. 5b and Supplementary Fig. 15). An integrated gradient attribution method27 was applied to obtain the importance score of each node in each layer for each patient in the patients set. The mean importance score across all patients was used as the final score for each node, which was displayed as the node length in the Sankey diagram. For the dual anti-HER2 regimen (Fig. 5b), high-scoring genes included PSMB3 (reported to be co-expressed with ERBB2 in breast cancer28), PSMB7 (reported associated with anthracycline resistance in breast cancer29), PIK3CA (reported associated with anti-HER2 treatment resistance in breast cancer30,31), YWHAH (reported associated with cell migration in breast cancer32), and SRC (reported associated with adhesion, and cell migration in breast cancer33). Pathways at various levels, such as ER-Phagosome pathway in hidden layer 2, which is reportedly linked with innate immune system and phagocytosis34, RHO GTPase cycle in hidden layer 3, which is reportedly linked with tumorigenesis, invasion, and metastasis35, and immune system in hidden layer 6, were identified, highlighting their critical role in anti-HER2 treatment. For the immunotherapy regimen (Supplementary Fig. 15), high-scoring genes included CCND1 (reported to be linked to enhanced invasiveness, a poorer prognosis, and resistance to drug therapies in ER-positive breast cancer36), SRC, JUN (reported to promote breast cancer progression via the AKT signalling pathway37), TRIO (reported to facilitate breast cancer cell migration38), and DUSP6 (reported to be associated with poor survival in HER2 positive breast cancer patients39). Collectively, we identified multiple key genes implicated in the response or resistance mechanisms to specific treatment regimens.

Methods
Publicly available data collection and processing
Transcriptome data related to breast cancer neoadjuvant therapy were systematically retrieved from PubMed and the Gene Expression Omnibus (GEO) database on 2024.2.15. The search terms used were (“Neoadjuvant Therapy”[MeSH Terms] OR neoadjuvant) AND (“Breast Neoplasms”[MeSH Terms] OR ((breast OR mammary) AND (cancer* OR tumor* OR tumour* OR carcinoma* OR neoplasm* OR oncology OR malignanc*))). “*” is a wildcard, which is used to match any sequence of zero or more characters. For example, tumor* can match tumor and tumors. High-throughput sequencing such as next-generation sequencing (NGS) and microarrays were included, and only datasets containing no less than 20 samples were selected. In total, 9550 breast cancer samples from 68 independent public datasets were accessed.

We obtained transcriptomic expression data in both raw and processed forms, encompassing NGS-derived count matrices or normalized (such as FPKM, TPM, or RPKM) value matrices and microarray raw intensities or normalized (such as RMA) values. Probesets of microarrays are annotated for genes referring to the reference table provided in GEO. For quality control of genes and samples, 16,133 genes occurred in larger than 50% samples and 4580 samples with less than 50% missing genes were retained. Among these datasets, 4371 patients from 31 datasets (GSE194040, GSE25066, GSE16716, GSE41998, GSE180962, GSE20271, GSE34138, GSE50948, GSE149322, GSE22226, GSE32603, GSE22358, GSE32646, GSE16446, GSE130788, GSE231629, GSE123845, GSE4779, GSE173839, GSE22093, GSE21997, GSE42822, GSE66399, GSE181574, GSE23988, GSE41656, GSE8465, GSE122630, GSE21974, GSE207248, GSE191127), covering data on complete pCR, treatment regimens, and sufficient transcriptome information were used for model construction, validation, and test.

Moreover, we collected the breast cancer neoadjuvant therapy transcriptome datasets obtained after February 15, 2024, as prospective test. A total of 134 patients out of 195 patients from two independent datasets (GSE240671, GSE260693) were included.

Rank value encoding of transcriptomes
RNA-seq count or patient-level normalization data (e.g., transcripts per kilobase million (TPM) or fragments per kilobase million (FPKM)), microarray patient-level normalized data (e.g., robust multichip average (RMA) or other method) from GEO were collected. A greater upstream level is maintained if multilevel data are accessed. To facilitate translation among gene expression testing technical platforms, we applied a rank transformation (rank value encoding) to normalize the data. For this purpose, gene expression values were transformed from microarray intensities or RNA-seq counts or their normalization format to their respective ranks. This transformation was conducted in a gene wise manner. Therefore, all gene expression values per gene were assigned a rank on the basis of their order from lowest to highest. The rankings were then scaled into quantiles, ranging from 0 (lowest gene expression) to 1 (highest gene expression). Missing values are filled with 0. This rank-based approach may be more robust against technical artifacts that may systematically bias the absolute transcript count value, whereas the overall relative ranking of genes within each cell remains more stable.

Treatment drug and regimen representation
To facilitate the fusion of drug treatments and transcriptomes, we aimed to represent drugs by combining related genes and correlation values. To build a sufficiently large pharmacogenomic dataset for drug representation, raw drug sensitivity data of cell lines to molecular drugs and RMA-normalized or TPM-normalized followed by log-transformed transcriptomes of the cell lines were retrieved from the Genomics of Drug Sensitivity in Cancer database v2 (GDSC) and the Cancer Therapeutics Response Portal v2 (CTRP)51–56. We used the IC50 metric, where a high IC50 value corresponds to low cell-killing ability and a low IC50 value corresponds to high cell-killing ability, to represent drug sensitivity. For each drug, we calculated the correlation between the gene expression levels of all the cell lines and the log2(IC50+1) of the response of the cell lines to each drug. Due to the fact that the normalized gene expression levels and log2(IC50+1) values showed a generally normal distribution, we applied the Pearson correlation method. The significance cutoff of correlation was set to 1*10−6 according to the Bonferroni correction to reduce the false discovery rate, considering nearly 5*104 genes to be tested. Genes identified as significant are more likely to play an important role in predicting the drug response of cancer treatment, and the correlation values are prior knowledge that can inform the importance degree of genes to drug response and thus can reflect the perturbation of drugs to genes. The correlated genes were removed if they were not included in the retained set of 16,133 genes in the patient transcriptome. Therefore, each drug was ultimately represented by its significant genes and correlation values (ranging from −1 to 1). If a drug appeared in both the GDSC and CTRP datasets, the GDSC dataset was considered first.

There are several special cases in drug representation. Some types of drugs are not found in the GSDC or CTRP datasets. Therefore, we chose the most relevant drug with similar mechanism and within the same anti-tumor drug type to represent it. Examples include capecitabine, which is represented by 5-fluorouracil; ganetespib, which is represented by luminespib; and ixabepilone, which is represented by epothilone B. Some drug descriptions in GEO datasets, such as taxane or anthracycline, lacks sufficient drug details; therefore, we chose an averaged representation of the drugs within the taxane or anthracycline categories. Some types of drugs are antibody- or peptide-based drugs, which are not included in either dataset and are not commonly represented by the IC50 to measure drug sensitivity due to complex mechanisms of action24–26. For these drugs, we used the negative correlation of their target genes with other genes in cell lines because the response of these drugs is strongly associated with target expression, including PDCD1 (pembrolizumab and durvalumab), VEGFA (bevacizumab), and IGF1R (ganitumab). If multiple targets exist for a drug, we performed a gene wise summation of the correlations among all the targets as the final representation, including ANG, which contains the ANGPT1 and ANGPT2 genes and represents an ANG peptide inhibitor (AMG-386). In addition, drugs such as T-DM1 (which is an antibody-drug conjugate) and endocrine therapy have no corresponding drug sensitivity data in the GDSC or CTRP. Considering the small number of patients treated with these drugs in our datasets, patients treated with regimens including T-DM1 or endocrine therapy were excluded. The representations of all drugs are listed in Supplementary Data 1. For patients receiving multiple drugs in a regimen, we performed a gene-wise summation of the gene correlation representations of each drug to represent the regimen.

Fusion of drug representations and transcriptome
Three approaches to integrate regimens and transcriptome information are tested. The first approach involves the use of both gene expression (represented as Ge) and gene-sensitivity correlation (represented as Gcd) as inputs, which we refer to as GDnet (Gene Drug net) (Fig. 1b). The second approach uses Gatt = Ge × exp(Gcd) as the final input, utilizing a key-value attention mechanism like approach (Supplementary Fig. 1a), and named this model GAnet (Gene-Attention net). In this mechanism, for each query (Ge), we find the corresponding key (gene) in the drug representation, and the attention (Gatt) is calculated as attention = query × exp(value), where the value represents the Gcd of the drug representation. The attention calculation is set to exp(value) to ensure that Gatt can simulate the tumor’s response under drug exposure conditions; in other words, Ge is expected to decrease or increase relative to the response (Gcd < 0) or resistance (Gcd > 0) gene of the drug when no drug is present. For example, ERBB2 is a poor-prognosis gene, but when exposed to an anti-HER2 drug with ERBB2 as the response gene (Gcd < 0), the expression of ERBB2 is expected to decrease. In addition, we explored a third approach, which involves the use of an ensemble learning approach. This method focuses on building different transcriptome-only models for each drug and combining the predictions of all transcriptome-only models corresponding to patient-applied drugs by averaging them for the final prediction. This model is referred to as Gnetens (Gene net ensemble) (Supplementary Fig. 1b). In Gnetens, only drugs with available data for more than 10% of the patients were used to construct individual models. Five drug types (anthracycline, microtubule inhibitors, cyclophosphamide, pyrimidine analogs, and anti-HER2 antibodies) were chosen (Supplementary Fig. 3a). Transcriptome-only model named as Gnet (Gene net) is also explored (Supplementary Fig. 1c).

Design of the biologically informed model
We introduced prior knowledge into the training process, drawing inspiration from the P-NET framework built by Elmarakeby HA et al.16 to accelerate model training and enhance both performance and explainability. The complete Reactome dataset57 was downloaded and processed into a layered network. This network includes one layer for features, one layer for genes, and five layers for pathways. The affiliation relationships between nodes in adjacent layers (representing genes or pathways) are represented by six 2D binary matrixes. The constructed layered network originated from 8,377 genes (the first hidden layer) that are available both in our datasets and the Reactome dataset. Pathways that are directly related to these genes form the second layer (846 initial pathways). Subsequently, pathways directly related to those in the second layer constitute the third layer (218 superior pathways), and this layering process continues analogously to the fourth layer (108 pathways), the fifth layer (51 pathways), until the sixth layer with the top 18 pathways representing biological processes is built.

The deep learning models were implemented in PyTorch (Python version 3.12.4; torch version 2.4.0; cuda version 12.2), with the layers and nodes numbers from the prior knowledge network defining the structure of a base feedforward dense neural network model. The six 2D binary matrixes containing affiliation relationships between genes or pathways in adjacent layers were fed into the neural network model as the masks between adjacent layers. This mask is multiplied by weights after each epoch during weight updates, imposing constraints on nodes and edges. This configuration yields a sparse neural network with six layers, comprising 8377-gene nodes (hidden layer 1), 846 (hidden layer 2), 218 (hidden layer 3), 108 (hidden layer 4), 51 (hidden layer 5), and 18 pathway nodes (hidden layer 6) in each hidden layer. Specifically, for GDnet and GDdnet, the input layer before the 8377-gene layer comprised gene expression and drug representation data encompassing 19,300 nodes (half gene expression nodes and half drug representation nodes). L2 regulation was set to 0 since the model is sufficiently sparse. Each node encodes a biological entity (for example, a gene or a pathway), while each edge represents a known relationship between the corresponding entities. The node constraints allow for a better understanding of the state of different biological components. The edge constraints lead to a smaller number of parameters than those of fully connected networks with the same number of nodes and thus potentially fewer computations. Since the dataset is unbalanced, we weighted the classes differently to reduce the network bias toward one class on the basis of the bias in the training set. The model was trained 1000 epochs via the Adam optimizer with default parameters to reduce the binary cross-entropy loss functions for imbalance. To allow each layer to be useful by itself, we added a predictive layer with sigmoid activation after each hidden layer. Since it is more challenging to fit the data using a smaller number of weights in the later layers, we used higher loss weights for later layer outcomes (from the gene & drug layer to the later layers, the weights are 1, 1.2, 1.4, 1.6, 1.8, 2, and 2.5) during the optimization process. The final prediction of the network was calculated by averaging all layer outcomes.

Models applying drug representation strategies in previous studies
Using the drug representation strategy based on SMILES (Simplified Molecular-Input Line-Entry System) molecule structure expressions proposed by Kuenzi et al.18, we introduced GSnet (Gene-SMILES net). In this approach, drugs are hashed into 2048-dimensional vectors to represent molecular drugs. Antibody-based drug representations were not generated due to the lack of SMILES entries. Regimens were represented as the vector summation of involved drugs. The SMILES-based 2048-dimensional vectors and gene expression vectors were subsequently used as input for GSnet.

Using the drug representation strategy based on the targets of drugs proposed by Zhang et al.19, we introduced GTnet (Gene-Target net). Drug target information was obtained from the DrugBank database58. Chemotherapeutic agents with nonspecific targets (e.g., DNA, RNA) that cannot be mapped to specific genes were excluded from drug representation. This approach incorporated the following 10 drugs with clearly mappable gene targets: trastuzumab (ERBB2), pertuzumab (ERBB2), bevacizumab (VEGFA), AMG-386 (ANGPT1, ANGPT2), MK-2206 (AKT1), olaparib (PARP1, PARP2), durvalumab (PDCD1), pembrolizumab (PDCD1), ganitumab (IGF1R), luminespib (HSP90AA1). Regimens were represented as the vector summation of involved drugs. The target-based 10-dimensional vectors and gene expression vectors were subsequently used as input for GTnet.

Hyperparameters optimization
Hyperparameter optimization was performed using 10-fold cross-validation on the training set to maximize the AUROC. The selected hyperparameters and their tuning ranges are described below:
(1) Hyperparameter settings for the informed group models: The GDnet, GAnet, Gnet, and Gnetens models incorporate biological prior knowledge into their architecture. Since their neural network structures (including nodes, layers, and connections) were predefined by this prior knowledge ([8388, 846, 218, 108, 51, 18] for GDnet, and [846, 218, 108, 51, 18] for GAnet and Gnet), no further optimization was required for these parameters. Thus, only the learning rate was tuned with candidates [0.1, 0.01, 0.001, 0.0005, 0.0001, 0.00005, 0.00001]. A learning rate of 0.0005 had the highest mean AUROC across all informed models and thus was set as the optimized learning rate.
(2) Hyperparameter settings for the dense group models: For GDdnet, GAdnet, Gdnet, and Gdnetens, to ensure comparability with prior knowledge models, the same node and layer structures were retained, but connections were converted to fully connected. Thus, only the learning rate was tuned, and the candidate values were the same as described above. A learning rate of 0.0005 was selected as the optimized learning rate.
(3) Hyperparameter settings for the optimized group models and previous group models: The GDdnet-o, GAdnet-o, Gdnet-o, Gdnetens-o, GSnet, and GTnet models employ fully connected architectures, requiring the tuning of layers, nodes, and learning rate. Learning rate used the same candidate values as above. For hidden layers, we used a range of 3 to 6 layers. For nodes, we used GDdnet-o as an example, the candidate architecture combinations include: [[8192, 4096, 2048], [4096, 1024, 256], [8192, 4096, 2048, 1024], [4096, 1024, 256, 64], [8192, 4096, 2048, 1024, 512], [4096, 1024, 256, 64, 16], [8192, 2048, 1024, 256, 128], [4096, 2048, 512, 256, 64], [8192, 4096, 2048, 1024, 512, 256], [4096, 1024, 256, 64, 16, 4], [8192, 2048, 1024, 256, 128, 32], [4096, 2048, 512, 256, 64, 32]]. The optimized learning rate was 0.00001 for GDdnet-o, and 0.0001 for others. The optimized hidden layers and nodes were [8192, 4096, 2048] for GDdnet-o, and [2048, 512, 128] for others.
(4) Hyperparameter settings for machine learning group models: We applied 4 machine learning models including L2-regularized logistic regression, random forest, support vector machine (SVM), and XGBoost in this study. Hyperparameters such as regularization strength (C) in L2 logistic regression, n_estimators, max_depth, and min_samples_split in random forest, regularization strength (C) and kernel in SVM, and n_estimators, max_depth, and learning_rate in XGBoost were tuned with empirical candidate values.

Ablation study for GDnet
The proposed GDnet model comprises three critical modules: the transcriptomic module, drug representation module, and prior knowledge module. To demonstrate the significance of these modules, we designed corresponding ablation experiments: (1) For the transcriptomic module, we developed a drug-representation-only model (Dnet) to validate the importance of transcriptomic data. (2) For the drug representation module, we used the transcriptome-only model (Gnet) and alternative drug representation-based models (GSnet and GTnet) to verify the superiority of our proposed drug characterization approach. (3) For the prior knowledge module, we used GDdnet, a fully connected feedforward neural network model without prior knowledge integration.

Model training and testing
The analysis included datasets that met specific criteria, providing a total of 4371 available patients. The 17 largest datasets, comprising approximately 3459 patients (79% of the total), were designated as the training set. The 14 remaining datasets, comprising approximately 912 (21% of the total), were used as the external test set, assigned without subjective assignment. The model was trained with the training set, and hyperparameters such as the learning rate were optimized via tenfold validation. The trained models based on the whole training set were finally validated with the external test set. The prediction performance of different deep learning models was compared using the test set by calculating the median AUROC and AUPRC across 10 rounds of replicated training using roc_auc_score and average_precision_score implemented in Python sklearn.metrics library separately. Then, we ensembled the 10 repetitions of the deep learning models to derive the final prediction for stability and robustness. Machine learning models were trained only once due to the same outcome across repetitions. The final models were evaluated according to the AUROC and AUPRC, and further evaluated using metrics such as accuracy, F1-score, precision, recall, and specificity with the cutoff value set as 0.5.

Model interpretation
We applied Integrated Gradients27, a gradient-based attribution method implemented in captum.attr library, to calculate the sample-level importance of all nodes in all layers. From 10-repetition GDnet models, we chose a GDnet model with a similar AUROC value (0.727) to that (0.725) of the 10-repetition ensemble GDnet for interpretation. The absolute importance score (always positive) of each layer’s nodes was normalized to facilitate comparability between layers via the following formula:
I j = A j / P n i=1 A i
where I j is the normalized absolute importance score for each node in a layer, A j is the absolute importance score, and n is the total number of nodes in a layer.

To calculate the total node-level importance, we aggregated the sample-level importance scores across all samples in the test set with a specific regimen (the dataset-level interpretation is represented as the mean value of all patient-level weights in the test set with a specific regimen). The importance score of each node is represented as the node length in the Sankey diagram (Fig. 5a). The importance score of links between the source node (in the former layer) and target node (in the latter layer) was represented by multiplying the normalized (within the absolute weight of all source nodes that point to the same target node) absolute weight of the source node that points to the target node in the model by the importance score of the target node, as shown in the following formula:
L (s→t) j = I t × P n s W j / P n s W j
L (s→t) j is the importance score of the link from the j source node (refer to as s) to the target node (refer to as t), which is the normalized absolute importance score of the target node, W s is the weight of the j source node pointing to the target node in the deep learning model, and n is the number of source nodes that point to a target node.

Statistical analyses
The change in the area under the ROC curve between GDnet and other models was tested via the DeLong test59. The AUROCs derived from tenfold cross-validation and tenfold repeated training and testing for GDnet and other models were compared via the Wilcoxon rank-sum test. The variance of the AUROCs and AUPRCs between Gnet and GNNet was compared via the Levene’s test to assess the robust of performance. The 95% confidence intervals for AUROC, AUPRC, accuracy, F1-score, precision, recall, and specificity were calculated using 1,000 times of bootstrapping. The odds ratio (OR) was calculated to determine preferences for pathological complete response (pCR) in the simulated trials, and differences were compared via Fisher’s exact test. A linear regression model was fitted to examine the relationships between the degree of optimization and OR. The significance of the regression coefficients was evaluated via t tests, and the R-squared value was used to quantify the proportion of variance explained by the model. To address potential confounding caused by the imbalance of transcriptome profile (inferred by the Gnet prediction score) between optimized and control groups, we employed Propensity Score Matching (PSM). PSM was conducted by estimating propensity scores using the transcriptome-only model Gnet prediction score. Matching was performed using nearest neighbor matching (1:1 ratio) with a caliper of 0.1, without replacement, to minimize imbalance. All tests were two-sided. P values < 0.05 were considered statistically significant unless otherwise specified. All the statistical analyses were performed via R statistical software version 4.3.

Statistics and reproducibility
No statistical method was used to predetermine sample size. We systematically retrieved patient data from the GEO database to maximize cohort size and statistical power, initially accessing 9745 patients. Patients were included regardless of sex or gender. Based on pre-established exclusion criteria, we excluded: (1) 628 duplicate samples to prevent data leakage and over-representation of individual patients; (2) 3805 samples collected after treatment initiation or lacking outcomes, as the model predicts treatment response from pre-treatment data; (3) 209 samples with >50% missing genes to maintain transcriptomic data integrity; and (4) 598 samples without treatment details or treatment representations, which are essential model inputs. The final training and validation cohort included 3459 patients from 17 datasets, all with complete pCR status, treatment regimens, and sufficient transcriptomic data. The testing cohort comprised 912 patients from 14 datasets, with an additional 134 patients from two subsequent datasets for prospective testing. Patient samples were allocated into experimental cohorts based on their source dataset, not random allocation. This strategy was chosen to evaluate the model’s generalizability across independent studies and different patient populations, which is a key strength of our study design. Blinding of investigators to cohort allocation during data collection and outcome assessment was not applicable to this computational study. The analysis was performed using pre-existing, de-identified data from the GEO database, and all model predictions were generated algorithmically, eliminating the potential for subjective bias.

To ensure computational reproducibility, we implemented: (1) tenfold cross-validation during training for hyperparameter optimization and internal validation; (2) 10 rounds of replicated training across 17 datasets with external testing in 14 datasets for robust comparison; and (3) prospective testing on an independent cohort of 134 patients from two subsequent datasets. All replication attempts successfully confirmed our initial findings, demonstrating model robustness.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_2": "Thank you for the opportunity to review the manuscript titled \"Biologically informed deep learning integrating drug representations for breast cancer treatment selection\". I have carefully evaluated the content of this paper. Below are my detailed comments:\nMajor Comments:\nStudy Design and Methodology:\nThe manuscript presents an interpretable deep learning model to predict treatment responses in breast cancer patients using drug representations and tumor transcriptome profiles. While the approach is technically sound, there is a critical limitation in the validation strategy. The authors rely entirely on publicly available datasets for model validation. I strongly recommend that the authors use their own institutional or cohort data for validation rather than relying solely on public datasets. Internal validation with an external dataset from their own hospital or research group would greatly enhance the credibility and clinical applicability of the findings. Public datasets, while valuable, do not always reflect the heterogeneous nature of patient populations across different institutions.\nClinical Relevance and Impact:\nWhile the model shows promise, the clinical relevance is underemphasized. The manuscript would benefit from a more thorough discussion of how this predictive model could be integrated into daily clinical practice. Specifically, the authors should explain how this model might influence clinical decision-making in real-world settings, especially in diverse patient populations. For instance, the model should be validated on their own patient data to ensure that the findings can be applied to their specific cohort and institution. The broader implications for patient management should be discussed in more depth.\nStatistical Analysis:\nThere are concerns regarding the statistical analysis, particularly in terms of model validation. The authors do not provide sufficient explanation for the choice of validation metrics, nor do they include an external validation cohort to confirm the model's generalizability. Without validating this model using data from the authors' own cohort, the model's true clinical utility cannot be accurately assessed. A more comprehensive validation strategy using both internal and external datasets would significantly strengthen the findings.\nInterpretation of Results:\nWhile the results are presented clearly, the interpretation lacks depth in terms of clinical impact. The relationship between the model outputs and actual clinical outcomes such as treatment responses, survival, and recurrence rates is not sufficiently explored. Furthermore, external validation using data from the authors' own institution would provide stronger evidence that the model can be effectively used in practice.\nModel Generalizability and Transparency:\nThe manuscript lacks transparency in terms of the model training process, especially regarding data preprocessing, feature selection, and model evaluation. While the authors use publicly available data, they should also describe how they would apply this model to their own institution's data. This would not only improve transparency but also help assess the model's generalizability across different clinical environments.\nEthical Considerations:\nGiven the reliance on public datasets and the potential for deploying this model in clinical settings, the manuscript should include a discussion on ethical issues such as data privacy, biases in the model, and the potential for inequitable outcomes in diverse patient populations. This is particularly important when using predictive models in clinical decision-making.\nMinor Comments:\nThe introduction would benefit from a more detailed review of the current literature, specifically in relation to predictive models used in breast cancer. A stronger connection to the existing body of research would provide better context for the model presented.\nSome sections of the manuscript are repetitive, especially the explanations of the methodology. A more concise and focused presentation would improve readability.\nA number of grammatical and typographical errors need to be addressed for clarity and flow.\nConclusion:\nWhile the authors have made a commendable attempt to develop a predictive model for treatment responses in breast cancer, the manuscript has significant limitations. The most critical issue is the lack of validation using the authors' own cohort data. Public datasets alone do not provide sufficient evidence for the model's robustness or clinical applicability. I strongly recommend that the authors validate their model using data from their own institution to enhance the model’s clinical relevance and generalizability. Without such validation, I do not believe this manuscript is ready for publication.", "reviewer_3": "The authors employed an attention mechanism-based deep learning approach, utilizing transcriptomic data from breast cancer patients receiving neoadjuvant therapy to predict therapeutic efficacy. Additionally, they explored the potential of AI-driven methodologies for optimizing treatment strategies. Despite their findings, questions persist concerning the model’s robustness and reliability when translated into real-world clinical applications. Below are a collection of comments and suggestions for your review:\n(1) Line 103: The sensitivity of targeted therapy is not necessarily correlated with the ex-pression levels of targeted genes. For instance, ERBB2 mutations can occur in HER2-low breast cancer and lead to the activation of downstream pathways, which can be targeted by anti-HER2 TKIs such as pyrotinib (mentioned in Yi Z. et al., Signal Transduction and Targeted Therapy, 2024). Therefore, would it be possible to integrate somatic mutations and tran-scriptomic profiles to predict therapeutic responses ?\n(2) Line 135: The sequencing data in GDSC and CTRP datasets derived solely from tumor cell lines. The therapeutic responses to anti-angiogenic therapies and anti-HER2 therapies depend on microenvironmental factors, more specifically, the abundance of tumor vascula-ture and the antibody-dependent cell-mediated cytotoxicity (ADCC) effect of T lymphocytes.\n(3) Line 139: The authors should clarify how the paired data, comprising drug sensitivity and expression profiles, are leveraged as prior knowledge to inform and enhance subsequent modeling.\n(4) Line 160: It is suggested that the authors include additional metrics, such as precision and F1-score, to assess the accuracy of model predictions and highlight the superiority of GDnet.\n(5) Line 190: The authors should verify whether the definition of pathological complete re-sponse (pCR) is consistent across different public datasets. The definition of pCR is partic-ularly stringent, requiring the absence of invasive carcinoma in both the primary lesion and lymph nodes, as well as the absence of vascular tumor emboli. However, some public da-tabases may interpret pCR solely as the complete response in primary lesion. Furthermore, in HR+HER2- breast cancer, pCR is not a strong prognostic indicator, which may affect the practical utility of this model.\n(6) Line 215: The translational research team of the I-SPY2 clinical trial utilized immune cell-related biomarkers to distinguish between durvalumab-sensitive and durval-umab-resistant populations (mentioned in Pusztai, L. et al. Cancer Cell, 2021). It is recom-mended that authors compare the predictive power of these biomarkers and their deep learning model in assessing durvalumab efficacy.\n(7) Line 244: The authors listed several pathways potentially suggesting therapeutic re-sponse to neoadjuvant therapy. However, as with the public datasets involved in this study, real-world neoadjuvant therapy regimens vary significantly and exhibit unique sensitivity and resistance pathways. Therefore, it would be wiser to select a specific regimen and discuss the regimen-specific sensitivity or resistance mechanisms. These insights could ultimately enhance the effectiveness of the model's interpretability analysis.", "reviewer_4": "This article developed a deep learning model called GDnet, which integrates drug representations and tumor transcriptomic data to predict breast cancer patients' responses to neoadjuvant therapy. By simulating drug-tumor interactions, GDnet outperforms transcriptome-only models in selecting optimal treatment strategies. Ultimately, the study demonstrates how biologically informed deep learning can optimize neoadjuvant treatment selection for breast cancer and potentially have broader applications across other treatment settings and cancer types. Detailed comments are listed here:\n1. Some descriptions are unclear: “To build a model with such functions, learning the interactions between tumor and treatments is necessary.” What do the “interactions between tumor and treatment” refer to? Please provide a more specific explanation.\n2. “We demonstrate that GDnet (area under the curve (AUC) = 0.725) significantly outperforms the transcriptome-only model (AUC = 0.683) in predicting treatment response.” Is the improvement in AUC statistically significant? A statistical test should be performed to justify this claim.\n3. What does the term “pan-cancer cell line” in Figure 1 mean? There is no explanation of this term in the text. Please clarify.\n4. An ablation study of the network structure should be conducted to demonstrate the importance of each component of the model.\n5. Typo corrections:\n• In Figure 1: “Gnome-wide” should be corrected to “Genome-wide.”\n• On page 13: “malignanc” should be corrected to “malignancy.”\n6. The second paragraph in the section “GDnet showed superiority in predicting the outcome of neoadjuvant therapy” should be rewritten. The illustrations of the results should be reorganized in a more coherent and logical way. Currently, too many methods and results are listed without clear organization.\n7. How were the numbers of nodes in each layer, such as the gene layer and pathway layers, determined? Please provide justification or additional details.\n8. How were the connections between nodes in adjacent layers determined? Clarify the rationale or methodology behind these connections.\n9. The biological meanings of “layer1,” “layer2,” and “layer3” are not clear. Additional explanations should be provided.\n10. The potential applications of the proposed method beyond breast cancer should be elaborated upon. More discussion is needed on how this method can support more informed decision-making in selecting neoadjuvant treatments.", "reviewer_5": "Summary of the work\nIn this work, Ge et al. presented GDnet, a deep learning model integrating drug representations with bulk transcriptomics to predict neoadjuvant therapy response in breast cancer. The approach of representing each drug as a set of genes with significant correlation to the corresponding cell line response is a resourceful and innovative approach that facilitates incorporating a large number of drugs into a single model effectively when using biologically informed DL. The validation of GDnet on a large number of patients highlights the potential of such approaches, with promising results on optimizing drug regimen selection. However, I have a couple of concerns that the authors should address to further improve the manuscript.\nComments to the authors\nOverall, the manuscript is concise and easy-to-follow with sufficient details for the most part, albeit with a number of typos and image quality issues. Here are the specific comments for further improvement:\n1. The authors have compared their GDnet model against a number of different DL models incorporating different architectures or data modalities. However, as a baseline, they should also compare their performance against a model based on known clinical factors such as tumor grade, ER status, HER2 status etc.\n2. Breast cancer has multiple established subtypes – clinical (based on ER/HER2 status) and molecular (based on PAM50 assay) with implications to patient response and survival. The authors should perform a subtype-specific to evaluate whether GDnet performs similar across different subtypes.\n3. There are a couple of missing technical info from the main text:\na. The authors mentioned using the genome-wide gene expression correlation with IC50 to represent different drugs. Is it Pearson or Spearman correlation?\nb. They also mentioned applying rank normalization on the transcriptomics data. is this normalization in a descending or ascending order?\n4. In the performance evaluation section (line 179-180), the authors mentioned “Though GDnet showed lower AUCs in the left 5 datasets, all the differences were less than 0.03”. Looking at Figure 3- it is not readily clear what they are referring to, especially with the difference being <0.03. The authors should revise it to clearly state what the intended message was.\n5. The interpretability analysis identified a few genes with high importance in response prediction. I was wondering if the authors could add the roles/biological significance of these genes in breast cancer (or cancer in general) from relevant literature.\n6. Figure 1, Extended Data Figure 1-2, 8-10 are too blurry and low-quality to properly examine. The authors should replace these figures with high quality versions.\n7. There are multiple typos across the manuscript that the authors should rectify, e.g.,\na. On page 5, line 116, they mentioned HA et al. which should be Elmarakeby et al.\nb. On page 9, line 244-248, the authors mentioned different pathways in different layers. However, the numbering of the layers is not consistent e.g., “ER-Phagosome pathway in layer 1” while it is actually on layer 2 of the network. They should resolve these discrepancies.\nc. On page 10, line 280, they mentioned RNA-sequencing methods which should be transcriptomics technologies (to include microarray techniques)\nd. On page 11, line 298, they mentioned “the 21-gene recurrence score”, which should be OncotypeDX score.", "reviewer_6": "This manuscript presents GDnet, a biologically-informed deep learning approach that integrates drug representations with transcriptomic data to predict and optimise breast cancer treatment response. The work makes several notable contributions: assembly of a large dataset (4,371 patients across 31 datasets), an innovative approach to drug representation using genome-wide expression correlations, and demonstration of potential clinical utility through in-silico trials. The study is methodologically rigorous with detailed model specifications, thorough ablation studies, and proper validation procedures. However, several aspects require attention. The benchmark comparison needs substantial expansion to include existing approaches from the literature. Some analysis results show counterintuitive patterns that require explanation, such as the unexpected relationship between optimisation thresholds and pCR rates, and the surprisingly better performance on test data compared to training data. The manuscript would also benefit from improved clarity in methodology description and figure quality. With these improvements, particularly in comparative evaluation and results analysis, this work could make a significant contribution to precision oncology.\nMajor:\nThe paper's core contribution needs better clarification and stronger empirical validation. While the authors have assembled an impressively large dataset (4,371 patients across 31 datasets), this alone isn't sufficient as a key contribution. The proposed GDnet model, which aims to integrate drug representations with transcriptomic data for treatment response prediction, requires more rigorous validation against existing approaches.\nThe authors describe three main approaches to drug characterisation in previous work: 1.Structure-based drug representation, 2.Target-based drug representation, 3.Genome-wide correlation-based representation (their approach).However, simply stating that previous approaches are unsuitable for this specific task isn't convincing without direct empirical comparison. The benchmark section should include implementations of at least 2-3 representative methods from these different approaches. While the extensive comparison of GDnet variants is valuable for ablation studies, it doesn't demonstrate advancement over the state-of-the-art.\nThe paper lacks crucial details and justification regarding the architecture of the dense neural network baselines. While the layer sizes for the biologically-informed models are well documented, the corresponding information for the dense models has to be inferred from the code. It appears the authors simply matched the layer sizes of their biologically-informed models, but this may not be optimal given the sample size and task complexity. For a fair comparison, the dense neural networks should be evaluated with different architectural configurations - varying both the number of layers and nodes per layer - along with proper hyperparameter optimisation.\nFurthermore, the comparison would be more complete by including traditional machine learning approaches that often perform well on tasks with this scale of data. Models like logistic regression, random forest, and XGBoost could serve as strong baselines, especially since they tend to be more sample-efficient than deep learning models. These simpler models might actually perform competitively given the relatively modest dataset size of 4,371 samples, and their inclusion would provide better context for evaluating the necessity of a deep learning approach. This additional benchmarking would strengthen the authors' claims about the advantages of their biologically-informed architecture.\nWhile the authors emphasise how biological knowledge improved performance in their GDnet architecture, they don't address why incorporating the same biological knowledge actually decreased performance in all other model variants (GAnet, Gnetens, and Gnet compared to their respective dense versions). This inconsistency raises questions about the general effectiveness of their biological knowledge integration approach and needs to be thoroughly discussed. Are there specific aspects of the GDnet architecture that make it uniquely suitable for biological knowledge integration? Or could this pattern suggest potential issues with the current implementation of biological knowledge in the other architectures? Without addressing these questions, the claimed benefits of biological knowledge integration remain uncertain.\nIn the PSM analysis with Fig 4, the odds ratios for optimised vs control groups show a stronger linear trend (R² = 0.94, p < 0.001) in the external validation compared to the training data (R² = 0.87, p < 0.001). The pCR rates also appear to show more significant differences in the validation set results. This pattern contradicts standard machine learning principles - a model should typically perform better on data it was trained on compared to unseen test data, unless there's underfitting. However, the overall performance metrics shown in Figure 3 don't suggest underfitting, as the model achieves reasonable AUC values. The authors need to investigate and explain this unusual pattern, as it raises questions about the robustness of their validation approach and the model's generalisation properties.\nThe relationships between optimisation threshold and pCR rates shown in Figure 4 display some counterintuitive patterns that require explanation. Theoretically, we would expect two clear trends: the optimised group should show decreasing pCR rates as the threshold becomes less stringent (includes more cases), while the control group should show increasing pCR rates. However, this expected pattern isn't consistently observed, particularly in the I-SPY2 training data results (4b-c). Most concerningly, the control group shows an illogical pattern where the worst 20% of cases actually have a higher pCR rate than the worst 50%. This makes little sense, especially in training data where the model should have learned to rank cases more reliably. The authors should explain this unexpected behaviour, as it raises questions about the robustness of their treatment optimisation approach.\nThe visualisation could also be improved by including more threshold steps between the current intervals. This would help readers better understand if these unexpected patterns are genuine trends or possibly artifacts of the chosen threshold points. A smoother curve with more data points would provide a clearer picture of how the optimisation truly performs across different stringency levels.\nMinor:\nThe description of \"gene and value pairs\" on lines 91-92 needs clarification. The concept is central to the method but isn't clearly explained.\nFigure 2 would be more informative if annotated with the number of gene-correlation pairs for each drug.\n\"AUC\" should be labelled as \"AUROC\" to explicitly distinguish it from AUPRC.\nSince Figure 3c and 3d present results averaged from ten runs, confidence intervals should be calculated and displayed to show the stability of these results.\nAll bar plots, including Figure 3d, should be annotated with sample sizes to help readers assess the reliability of results across different datasets.\nThe current Extended Data Figure 5a, which shows a more complete view than Figure 3d, would be more appropriate as a main figure panel. The current arrangement buries important information in supplementary materials.\nThere's an interesting pattern in Extended Data Figure 5 where GANNet outperforms GDnet on several datasets. What characteristics of these datasets might explain this performance difference?\nThe Figure 4 legend needs reorganisation and better descriptions to match the panel order and clearly explain the experimental setup.\nThere's a typo on line 221 where \"detests\" should presumably be \"datasets\".\nMost figures look blurred when zoomed in."}, "author_response": "Reviewer #2\nReviewer #2 (Remarks to the Author): Expert in breast cancer neoadjuvant\nchemotherapy, clinical trials, and statistics\nThank you for the opportunity to review the manuscript titled \"Biologically informed\ndeep learning integrating drug representations for breast cancer treatment selection\". I\nhave carefully evaluated the content of this paper. Below are my detailed comments:\nMajor Comments:\nStudy Design and Methodology:\nThe manuscript presents an interpretable deep learning model to predict treatment\nresponses in breast cancer patients using drug representations and tumor transcriptome\nprofiles. While the approach is technically sound, there is a critical limitation in the\nvalidation strategy. The authors rely entirely on publicly available datasets for model\nvalidation. I strongly recommend that the authors use their own institutional or cohort\ndata for validation rather than relying solely on public datasets. Internal validation with\nan external dataset from their own hospital or research group would greatly enhance the\ncredibility and clinical applicability of the findings. Public datasets, while valuable, do\nnot always reflect the heterogeneous nature of patient populations across different\ninstitutions.\nResponse: Thank you for pointing out validating the model using patient populations\nfrom diverse institutions with high heterogeneity to enhance its credibility and clinical\napplicability. We apologize for the lack of clarity in our validation strategy, which may\nhave given the impression that only internal validation was performed. We have now\nrevised the flowchart illustrating the model training and validation strategy\n(Supplementary Fig. 2). As shown, we utilized 17 independent datasets from different\ninstitutions, clinical trial cohorts, and transcriptomic profiling methods (Supplementary\nTable 2) for hyperparameter tuning and 10-fold interval cross-validation, followed by\nretraining the model with the optimal hyperparameters. Subsequently, we\nconducted external multicenter test based on 14 additional independent datasets from\ndifferent institutions, encompassing varied clinical trial cohorts and transcriptomic\nplatforms (Supplementary Table 2). Thus, our validation set meets the criteria of\nbeing external multicenter and highly heterogeneous.\nAdditionally, we noted that after our initial submission, two new breast cancer\nneoadjuvant therapy transcriptomic datasets from different centers were publicly\nreleased. To further strengthen our model’s validation, we performed a prospective\nexternal test using these two datasets (Supplementary Fig. 2b, Supplementary Table 2),\nfurther confirming the model’s credibility and clinical applicability.\nRegarding validation within our own institution, our group has not previously conducted\nneoadjuvant therapy studies, so no corresponding archived specimens were available. We\nattempted to retrospectively collect Formalin Fixed Paraffin Embedded (FFPE) tissues\nfrom neoadjuvant-treated cases for transcriptomic sequencing, but this was declined by\nthe pathology department, as it would have nearly exhausted the patients' specimens. On\nthe other hand, we are currently conducting a breast cancer neoadjuvant therapy clinical\ntrial (NCT05983094), in which we have planned prospective collection of pretreatment\nbiopsy specimens from enrolled patients. However, as the trial is still in the active\nrecruitment phase, model validation using this cohort will not be feasible in the short\nterm.\nGiven that our model has already undergone external validation based on 16 independent\nmulticenter patient cohorts, we believe its credibility and clinical applicability have been\nrobustly demonstrated. The absence of validation at our institute does not diminish the\nevidence level supporting the model’s performance.\nThank you again for your valuable feedback, which has helped us improve the clarity and\nrigor of our validation approach.\nClinical Relevance and Impact:\nWhile the model shows promise, the clinical relevance is underemphasized. The\nmanuscript would benefit from a more thorough discussion of how this predictive model\ncould be integrated into daily clinical practice. Specifically, the authors should explain\nhow this model might influence clinical decision-making in real-world settings,\nespecially in diverse patient populations. For instance, the model should be validated on\ntheir own patient data to ensure that the findings can be applied to their specific cohort\nand institution. The broader implications for patient management should be discussed in\nmore depth.\nResponse: Thank you for raising the important question regarding the practical\napplication value of our model. In response, we have expanded the discussion in our\nmanuscript (Lines 366-385) to include a detailed analysis of the model’s clinical utility,\nfocusing on the following three key applications:\n(1) In real-world clinical practice, for a breast cancer patient scheduled for neoadjuvant\ntherapy, we can first list all eligible treatment options based on guidelines and clinical\ntrial evidence. The patient’s tumor transcriptome data and treatment options can then be\ninput into GDnet to predict the probability of achieving a pCR. The regimen with the\nhighest pCR probability can then be recommended as the optimal treatment choice for the\npatient. The efficacy of this treatment selection strategy is demonstrated by in-silico\nsimulation trials based on ISPY-2 and external test dataset, showing significantly\nimproved pCR rate compared to the control group in our study.\nFurthermore, to further validate this treatment selection strategy with higher-level\nevidence, we are planning a randomized clinical trial based on GDnet. The trial will\ncompare GDnet-guided therapy selection (intervention arm) against standard-of-care\ntreatment (control arm), with the primary endpoint being the difference in pCR rates. The\ntrial protocol and ethical approvals are currently underway.\n(2) In clinical trial design, GDnet can assist in optimizing treatment strategies.\nSpecifically, given a set of candidate monotherapies or combination regimens, GDnet can\nsimulate and predict the pCR probabilities for each regimen option in patient\nsubpopulations that meet the inclusion criteria. The regimen demonstrating the highest\noverall pCR probability across the patient subpopulations would be the most promising\ncandidate for clinical trial. Even for novel drug combinations not included in GDnet’s\ntraining set, the model can still guide trial strategy selection. This is evidenced by\nGDnet’s strong performance in the external test set, where it successfully predicted the\nefficacy of drugs that were not included in model training (e.g., durvalumab,\nbevacizumab, cisplatin, olaparib, and gemcitabine).\n(3) Beyond trial design, GDnet can further optimize patient selection for clinical trials.\nBy estimating individual patients’ likelihood of responding to a trial regimen, GDnet can\nidentify the subgroup most likely to benefit, enabling selective enrollment and thereby\nimproving trial success rates.\nRegarding validation on our institute patient data, as mentioned earlier, we have\nconducted external validation using a total of 14 external datasets from different centers\nplus 2 prospective datasets, which fully demonstrates our model's credibility and clinical\napplicability.\nWe sincerely appreciate your insightful feedback, which has helped us better articulate\nthe clinical translation potential of GDnet.\nStatistical Analysis:\nThere are concerns regarding the statistical analysis, particularly in terms of model\nvalidation. The authors do not provide sufficient explanation for the choice of validation\nmetrics, nor do they include an external validation cohort to confirm the model's\ngeneralizability. Without validating this model using data from the authors' own cohort,\nthe model's true clinical utility cannot be accurately assessed. A more comprehensive\nvalidation strategy using both internal and external datasets would significantly\nstrengthen the findings.\nResponse: We appreciate your suggestion regarding the use of more objective and\ncomprehensive evaluation metrics. For model assessment, we selected AUROC and\nAUPRC as our primary metrics since they evaluate the model’s performance across all\npossible threshold settings, providing a more comprehensive evaluation than other\nmetrics like F1-score or accuracy. These metrics are well-established in neoadjuvant\ntherapy modeling studies1-5.\nTo complement these primary metrics, we systematically incorporated five additional\nevaluation measures (F1-score, accuracy, recall, precision, and specificity) using a cutoff\nthreshold of 0.5. We recognize that these single-threshold metrics are less comprehensive\nthan AUROC/AUPRC. Especially, in clinical practice, the optimal cutoff threshold is\ninherently flexible and must be carefully tailored to specific clinical scenarios. Among\nthese metrics, we prioritized F1-score (harmonic mean of precision and recall) and\naccuracy (proportion of correct predictions) as they provide balanced assessments of both\npositive and negative predictions. The remaining metrics (recall, precision, specificity),\nwhile valuable, offer more limited perspectives as they focus on either positive or\nnegative predictions alone. As demonstrated in Supplementary Fig. 7, GDnet achieved\nbalanced performance across all these metrics in external test datasets, outperforming\nmost comparative models. These results further validate the reliability and clinical\napplicability of GDnet.\nRegarding validation with our institutional data, as previously mentioned, the model has\nbeen externally validated using datasets from 14 independent centers along with 2\nprospective datasets. This extensive multi-center validation provides robust evidence\nsupporting the model's reliability and clinical utility.\nInterpretation of Results:\nWhile the results are presented clearly, the interpretation lacks depth in terms of clinical\nimpact. The relationship between the model outputs and actual clinical outcomes such as\ntreatment responses, survival, and recurrence rates is not sufficiently explored.\nFurthermore, external validation using data from the authors' own institution would\nprovide stronger evidence that the model can be effectively used in practice.\nResponse: We sincerely appreciate your valuable suggestions regarding additional\nclinically important endpoints worthy of prediction in breast cancer beyond pCR,\nincluding ORR, EFS, and OS. In our current study, we selected pCR as the primary\nprediction target and trained our model using pCR as the label. This decision was based\non several considerations: pCR is a well-established and widely used surrogate endpoint\nfor neoadjuvant therapy in breast cancer1, 3-8, with numerous studies demonstrating its\nstrong correlation with various long-term outcomes9. Furthermore, pCR represents the\nmost generally reported outcome measure across most neoadjuvant therapy studies and\npublicly available datasets (Supplementary Table 2).\nWe fully acknowledge the significant clinical relevance of other endpoints such as ORR\nand survival outcomes (EFS, OS), which may not be completely captured by pCR alone.\nTo address this, we have conducted a comprehensive review of available outcome\ninformation (ORR, EFS, DRFS, OS) across all included datasets, as presented in\nSupplementary Table 2. However, our analysis revealed that most datasets lacked these\nadditional outcome measures - only 2 out of 34 datasets reported ORR, 4 reported EFS or\nDRFS, and 2 reported OS. This limited availability currently precludes the development\nof robust predictive models for these endpoints. To advance predictive modeling in this\nfield, future efforts should prioritize the building of datasets encompassing complete\noutcome measures - including pCR, ORR, and comprehensive survival metrics - to better\ncapture the full spectrum of patient prognosis.\nRegarding validation using our own patient data, as previously detailed, our model has\nundergone extensive external validation utilizing datasets from 14 independent centers\nalong with 2 prospective datasets. We believe this external multi-center validation\nstrategy effectively addresses concerns regarding model generalizability in practice.\nModel Generalizability and Transparency:\nThe manuscript lacks transparency in terms of the model training process, especially\nregarding data preprocessing, feature selection, and model evaluation. While the authors\nuse publicly available data, they should also describe how they would apply this model to\ntheir own institution's data. This would not only improve transparency but also help\nassess the model's generalizability across different clinical environments.\nResponse: Thank you for pointing out the need for a more detailed description of our\ndata preprocessing methods. We have now expanded the technical details regarding data\npreprocessing, feature selection, and model evaluation in the Methods section (Lines 422-\n460).\nRegarding data preprocessing, we have provided a more comprehensive description of\nour procedures. We obtained transcriptomic expression data in both raw and processed\nforms, encompassing NGS-derived count matrices or normalized (such as FPKM, TPM,\nor RPKM) value matrices and microarray raw intensities or normalized (such as RMA)\nvalues (Lines 431-435). For NGS data specifically, since GEO datasets already contained\ngene expression level data, we utilized these pre-processed data directly rather than\nprocessing from raw fastq files.\nFor feature engineering, we retained 16,133 genes that appeared in more than 50% of\nsamples while removing genes with higher missing rates to ensure proper model training\n(Lines 435-437). It's important to note that this study employs a deep learning-based\nmodeling strategy. Unlike simple machine learning models, deep learning models are\nbuilt by incorporating all available features to the greatest extent2, 10-12. Therefore, we did\nnot perform feature selection steps such as collinearity removal or k-best feature selection\nbut rather included all features as much as possible to train the deep learning model.\nConcerning validation with our own institute’s data, as previously stated, we have\nconducted extensive external validation using datasets from 14 independent centers along\nwith 2 prospective datasets. This rigorous validation process robustly demonstrates our\nmodel's generalizability across diverse clinical environments. Moreover, we validated our\nmodel's robust performance across diverse breast cancer subtypes in external test datasets\n(Figure 3d; Supplementary Fig. 8), highlighting its complementary role alongside\nestablished clinical factors in practical settings.\nEthical Considerations:\nGiven the reliance on public datasets and the potential for deploying this model in clinical\nsettings, the manuscript should include a discussion on ethical issues such as data\nprivacy, biases in the model, and the potential for inequitable outcomes in diverse patient\npopulations. This is particularly important when using predictive models in clinical\ndecision-making.\nResponse: Thank you for raising the important ethical issues regarding the application of\nour model to clinical practice.\nRegarding data privacy, the dataset our model relies on consists of de-identified, secure\npublic data, so there are no data security concerns. We have also added a discussion in the\npaper about data privacy for future clinical applications of the model (Lines 385-388). We\nrecommend storing all data and deploying the model exclusively within hospital firewalls\nto mitigate breach risks.\nAs for biases in the model, we conducted subgroup analyses on patients with different\nmolecular subtypes of breast cancer and found that the model performs well across all\nsubgroups, indicating that our model covers all molecular subtypes with low bias (Figure\n3d). For other potential biases, such as age, menopausal status, tumor stage, and grade,\nfurther validation is needed in the future due to the lack of such information in existing\npublic databases.\nRegarding the potential for inequitable outcomes in diverse patient populations, our\nmodeling dataset includes individuals from multiple regions, such as the Americas,\nEurope, and Asia (Supplementary Table 2), so our model is applicable to diverse patient\npopulations. Additionally, GDnet is based on transcriptomic data, which differs from\ngenomic data that may involve issues like SNP polymorphisms across ethnicities.\nTranscriptomics reflects biological functional states, transcending genotypic differences\namong ethnicities. Therefore, we believe our model is minimally affected by ethnic\nvariations, and this discussion has also been added to the paper (Lines 388-392).\nMinor Comments:\nThe introduction would benefit from a more detailed review of the current literature,\nspecifically in relation to predictive models used in breast cancer. A stronger connection\nto the existing body of research would provide better context for the model presented.\nResponse: Thank you for your suggestion. We have supplemented and refined this\nsection (Lines 37-46), making the connection between the background and our study\nclearer.\nSome sections of the manuscript are repetitive, especially the explanations of the\nmethodology. A more concise and focused presentation would improve readability.\nResponse: Thank you for pointing out the issue of redundancy in our article. We have\ncarefully polished the entire manuscript and removed the repetitive sections, particularly\nby condensing redundant methodological details in the Results, thereby improving the\nreadability of the paper.\nA number of grammatical and typographical errors need to be addressed for clarity and\nflow.\nResponse: Thank you for your careful review. We have thoroughly examined and\npolished the entire manuscript, correcting grammatical and typographical errors.\nConclusion:\nWhile the authors have made a commendable attempt to develop a predictive model for\ntreatment responses in breast cancer, the manuscript has significant limitations. The most\ncritical issue is the lack of validation using the authors' own cohort data. Public datasets\nalone do not provide sufficient evidence for the model's robustness or clinical\napplicability. I strongly recommend that the authors validate their model using data from\ntheir own institution to enhance the model’s clinical relevance and generalizability.\nWithout such validation, I do not believe this manuscript is ready for publication.\nResponse: Thank you for raising a series of critical questions. We have addressed each of\nyour suggestions and revised the manuscript accordingly. In particular, we introduced two\nadditional test datasets from different centers to validate the model. These suggestions\nhave significantly enhanced the quality of our manuscript.\nReviewer #3\nReviewer #3 (Remarks to the Author): Clinical expert in breast cancer neoadjuvant\nchemotherapy and clinical trials\nThe authors employed an attention mechanism-based deep learning approach, utilizing\ntranscriptomic data from breast cancer patients receiving neoadjuvant therapy to predict\ntherapeutic efficacy. Additionally, they explored the potential of AI-driven methodologies\nfor optimizing treatment strategies. Despite their findings, questions persist concerning\nthe model’s robustness and reliability when translated into real-world clinical\napplications. Below are a collection of comments and suggestions for your review:\nResponse: Thank you for raising a series of important questions. We have responded to\neach of your questions and revised the manuscript accordingly. Regarding the model’s\nrobustness and reliability when translated into real-world clinical applications, we further\nvalidated the model by introducing additional two independent test datasets,\ndemonstrating its strong generalization capability. We also validated GDnet's robust\nperformance across diverse breast cancer subtypes in external test datasets (Figure 3d;\nSupplementary Fig. 8), highlighting its complementary role alongside established clinical\nfactors in practical settings. Moreover, in durvalumab-treated patients, GDnet achieved\npredictive performance comparable to or better than the top three immune-related\nsignatures mentioned in Pusztai, L. et al. (Cancer Cell, 2021) for distinguishing between\ndurvalumab-sensitive and durvalumab-resistant populations. These results further\nvalidate GDnet's robustness and reliability.\n(1) Line 103: The sensitivity of targeted therapy is not necessarily correlated with the ex-\npression levels of targeted genes. For instance, ERBB2 mutations can occur in HER2-low\nbreast cancer and lead to the activation of downstream pathways, which can be targeted\nby anti-HER2 TKIs such as pyrotinib (mentioned in Yi Z. et al., Signal Transduction and\nTargeted Therapy, 2024). Therefore, would it be possible to integrate somatic mutations\nand tran-scriptomic profiles to predict therapeutic responses ?\nResponse: Thank you for your suggestion. You raised an important point: the\nmechanisms underlying targeted therapy sensitivity depend not only on gene expression\nlevels but also on factors such as genetic mutations. Therefore, integrating transcriptomic\nand genomic data would provide a more comprehensive characterization of tumor states,\nenabling more accurate prediction of targeted therapy drug sensitivity1.\nDuring our study design, we also considered multi-omics fusion models. While\ncompiling the dataset, we concurrently collected available genomic data from the\nincluded datasets (Supplementary Table 2). However, only 2 out of 34 datasets\n(GSE123845 with 105 samples and GSE191127 with 20 samples) contained paired\ntranscriptomic and genomic data, totaling just 125 patients, and the genomic data was not\npublicly accessible. With such a limited sample size, constructing a robust, integrated\ntranscriptomic-genomic predictive model for treatment efficacy was unfeasible, which is\nwhy genomic data were not included in this study. Future efforts should prioritize the\nbuilding of multi-omics datasets to facilitate the development of more precise multi-\nomics predictive models.\nHowever, as you pointed out, while genetic mutations may not be directly reflected in the\nexpression levels of the mutated genes themselves, they can influence the expression\npatterns of upstream and downstream genes. Since our predictive model is built using\ndeep learning and incorporates expression levels of all available genes as input features, it\nis theoretically capable of learning the altered expression patterns of related genes caused\nby HER2 mutations based on transcriptomic data. This allows our model to partially\ncompensate for the missing genomic information, which was one of the key\nconsiderations in our choice of deep learning for this study.\n(2) Line 135: The sequencing data in GDSC and CTRP datasets derived solely from\ntumor cell lines. The therapeutic responses to anti-angiogenic therapies and anti-HER2\ntherapies depend on microenvironmental factors, more specifically, the abundance of\ntumor vascula-ture and the antibody-dependent cell-mediated cytotoxicity (ADCC) effect\nof T lymphocytes.\nResponse: We sincerely appreciate your valuable feedback regarding the drug\nrepresentation approach. You have pointed out that the therapeutic mechanisms of anti-\nangiogenic agents, anti-HER2 targeted therapies, and immunotherapies extend beyond\nthe expression levels of their direct targets to encompass the tumor immune\nmicroenvironment state. Therefore, the representation of these targeted drugs should\nincorporate feature dimensions that simulate drug-induced perturbations of tumor\nvascular status, ADCC effects, and immune microenvironment modulation.\nIn this study, in addition to target genes, we employed over 1,000 genes whose expression\ncorrelates with therapeutic targets in tumor cell lines to represent targeted therapies.\nAlthough these correlations were derived from tumor cell line data, the actual tumor and\nmicroenvironment interact and influence each other, resulting in strong correlations\nbetween tumor state and microenvironment status. The expanded representation\ndimensions provide the training foundation (and parameter space) for modeling drug\nperturbation mechanisms beyond direct target effects, including microenvironmental\ninfluences. Specifically, for anti-angiogenic drugs, VEGFA and its correlated 2,028 genes\ncan influence angiogenesis in the microenvironment, enabling partial simulation of\nvascular perturbations under such high-dimensional drug representation. For\nimmunotherapies, PDCD1 and its correlated 3,697 genes also impact the immune\nmicroenvironment, allowing partial simulation of immunotherapy-induced immune\nmicroenvironment perturbations (as shown in Supplementary Fig. 15 where drug\nresponse-related pathways included immune-associated pathways such as Signaling by\nInterleukins and Immune System). For anti-HER2 therapies, HER2 and its correlated\n2,382 genes affect ADCC-related components in the microenvironment like NK cells,\nmacrophages, and neutrophils, enabling partial simulation of ADCC effect perturbations\n(as demonstrated in Figure 5b where anti-HER2 regimen response-related pathways\ninclude innate immune system-related pathways such as ER-Phagosome pathway and\nSARS-CoV Infection pathway).\nWhile our current drug representation approach may not be perfect, it represents the\noptimal available method we have identified that simultaneously satisfies three critical\ncriteria: unified representation across all types of targeted therapies, maximal simulation\nof drug efficacy and resistance mechanisms, and seamless integration of drug\ncharacterization with tumor transcriptomes. Undoubtedly, future approaches will enable\nmore comprehensive characterization of targeted therapies across all mechanistic levels\nand more complete representation of patient tumor and microenvironment states,\ntheoretically achieving better simulation of drug-induced tumor perturbations and\ntreatment response prediction. However, this will require larger-scale single-cell multi-\nomics datasets that are currently unavailable. With the advancement of single-cell omics\ntechnologies and the accumulation of relevant datasets, we will pursue further\nimprovements in this area in the future.\n(3) Line 139: The authors should clarify how the paired data, comprising drug sensitivity\nand expression profiles, are leveraged as prior knowledge to inform and enhance\nsubsequent modeling.\nResponse: Thank you for your suggestion. We apologize for the unclear wording at the\nend of this paragraph, which caused confusion. The gene-correlation pairs of drug\nrepresentations and gene-expression pairs of tumor transcriptome were fed into the gene\nlayer of GDnet described in Figure 1b. Since the nodes and connections in GDnet's deep\nlearning network are designed based on biological prior knowledge of gene-pathway\nrelationships from databases like Reactome, the information flows from inputs (drug\nrepresentations and transcriptomes) to deep learning hidden layers (gene layer and\nmultiple pathway layers) is guided by these biological knowledges. We have made\ncorresponding revisions in the manuscript (Lines 121-126; Lines 140-142).\n(4) Line 160: It is suggested that the authors include additional metrics, such as precision\nand F1-score, to assess the accuracy of model predictions and highlight the superiority of\nGDnet.\nResponse: Thank you for your valuable suggestion regarding the use of additional\nevaluation metrics. In our study, we selected AUROC and AUPRC as the primary\nevaluation metrics because they assess the model’s performance across all possible\nthreshold settings, providing a more comprehensive evaluation than other metrics like F1-\nscore or accuracy. These metrics are well-established standards in neoadjuvant therapy\nmodeling studies1-5.\nTo complement these primary metrics, we systematically incorporated five additional\nevaluation measures (F1-score, accuracy, recall, precision, and specificity) using a cutoff\nthreshold of 0.5. We recognize that these single-threshold metrics are less comprehensive\nthan AUROC/AUPRC. Especially, in clinical practice, the optimal cutoff threshold is\ninherently flexible and must be carefully tailored to specific clinical scenarios. Among\nthese metrics, we prioritized F1-score (harmonic mean of precision and recall) and\naccuracy (proportion of correct predictions) as they provide balanced assessments of both\npositive and negative predictions. The remaining metrics (recall, precision, specificity),\nwhile valuable, offer more limited perspectives as they focus on either positive or\nnegative predictions alone. As demonstrated in Supplementary Fig. 7, GDnet achieved\nbalanced performance across all these metrics in external test datasets with a cutoff\nthreshold of 0.5, outperforming most comparative models. These results further validate\nthe reliability and clinical applicability of our approach.\n(5) Line 190: The authors should verify whether the definition of pathological complete\nre-sponse (pCR) is consistent across different public datasets. The definition of pCR is\npartic-ularly stringent, requiring the absence of invasive carcinoma in both the primary\nlesion and lymph nodes, as well as the absence of vascular tumor emboli. However, some\npublic da-tabases may interpret pCR solely as the complete response in primary lesion.\nFurthermore, in HR+HER2- breast cancer, pCR is not a strong prognostic indicator,\nwhich may affect the practical utility of this model.\nResponse: We sincerely appreciate your reminders. We have carefully reviewed the GEO\ndataset descriptions and corresponding publications to verify the definition of pCR across\nall datasets. We confirmed that pCR consistently refers to complete response in both\nprimary tumor and regional lymph nodes without residual lesions (Supplementary Table\n2).\nRegarding HR+HER2- breast cancer, your question is exceptionally insightful. Compared\nto other subtypes, HR+HER2- breast cancer exhibits a higher incidence of late recurrence\n(beyond 5 years), making pCR - as a short-term efficacy indicator - potentially less\npredictive of long-term outcomes. Additionally, these patients typically require\npostoperative endocrine therapy, which further influences long-term prognosis and\ncomplicates the correlation between pCR and survival outcomes. However, considering\nthat neoadjuvant therapy primarily aims to achieve tumor volume reduction in the short\nterm and guide postoperative adjuvant or intensified treatment strategies, using pCR as\nthe primary endpoint aligns with current clinical objectives of neoadjuvant treatment.\nWhile other endpoints like EFS and OS indeed provide more prognostic value,\ndeveloping predictive models incorporating neoadjuvant therapy, surgery, and adjuvant\ntreatment strategies to forecast EFS/OS represents an important future direction. As\nshown in Supplementary Table 2, most datasets in this study lack survival information -\nonly 4/34 datasets include EFS or DRFS data, and 2 datasets contain OS information -\nwhich is currently insufficient for building robust predictive models. We will monitor\nupdates to these datasets and pursue development of more prognostically relevant models\nwhen adequate data becomes available.\n(6) Line 215: The translational research team of the I-SPY2 clinical trial utilized immune\ncell-related biomarkers to distinguish between durvalumab-sensitive and durvalumab-\nresistant populations (mentioned in Pusztai, L. et al. Cancer Cell, 2021). It is\nrecommended that authors compare the predictive power of these biomarkers and their\ndeep learning model in assessing durvalumab efficacy.\nResponse: We sincerely appreciate your suggestion. For comparison with GDnet model,\nwe selected three signatures from Pusztai, L. et al. (Cancer Cell, 2021) that showed\nsignificant differences between responders and non-responders across all molecular\nsubtypes in durvalumab-treated patients: the STAT1 signature, TAMsurr_TcClassII ratio\nsignature, and Mitotic signature.\nThe results demonstrated that in durvalumab-treated patients, GDnet achieved superior\nperformance compared to both the TAMsurr_TcClassII ratio signature and Mitotic\nsignature, while performing comparably to the STAT1 signature. This validates the high\npredictive capability of GDnet for evaluating durvalumab treatment response.\nNotably, in patients who did not receive durvalumab treatment, GDnet showed\nsubstantially better performance than all three signatures.\n(7) Line 244: The authors listed several pathways potentially suggesting therapeutic re-\nsponse to neoadjuvant therapy. However, as with the public datasets involved in this\nstudy, real-world neoadjuvant therapy regimens vary significantly and exhibit unique\nsensitivity and resistance pathways. Therefore, it would be wiser to select a specific\nregimen and discuss the regimen-specific sensitivity or resistance mechanisms. These\ninsights could ultimately enhance the effectiveness of the model's interpretability\nanalysis.\nResponse: Thank you for your suggestion, this is a very good point. We fully agree with\nyour perspective that molecular mechanisms of drug response or resistance are treatment\nspecific. Therefore, it makes more sense to interpret models for specific treatment\nregimens.\nFor our analysis, we selected treatment regimens that appeared significantly more\nfrequently in the optimized group compared to the control group in our simulated trials:\nthe dual anti-HER2 regimen (paclitaxel|doxorubicin|cyclophosphamide|ERBB2|ERBB2)\nand the immunotherapy regimen\n(paclitaxel|doxorubicin|cyclophosphamide|olaparib|PDCD1) (Figure 5a). We then\nperformed model interpretation on the biological entities (genes and pathways) for each\nregimen separately.\nFor the dual anti-HER2 regimen, gene-level analysis identified PIK3CA, PSMB3,\nYWHAH, PSMB7, and SRC as key genes associated with treatment response (Figure 5b).\nFor the immunotherapy regimen, TRIO, JUN, DUSP6, CCND1, and SRC were found to\nbe important for response (Supplementary Fig. 15). Brief functional descriptions of these\ngenes are also provided in the manuscript (Lines 272–286).\nReviewer #4\nReviewer #4 (Remarks to the Author): Expert in machine learning, deep learning, cancer\nomics, and drug response\nThis article developed a deep learning model called GDnet, which integrates drug\nrepresentations and tumor transcriptomic data to predict breast cancer patients' responses\nto neoadjuvant therapy. By simulating drug-tumor interactions, GDnet outperforms\ntranscriptome-only models in selecting optimal treatment strategies. Ultimately, the study\ndemonstrates how biologically informed deep learning can optimize neoadjuvant\ntreatment selection for breast cancer and potentially have broader applications across\nother treatment settings and cancer types. Detailed comments are listed here:\n1. Some descriptions are unclear: “To build a model with such functions, learning the\ninteractions between tumor and treatments is necessary.” What do the “interactions\nbetween tumor and treatment” refer to? Please provide a more specific explanation.\nResponse: We sincerely appreciate your careful review. Indeed, the original phrasing of\nthis statement was not clear. It is well established that drugs can induce molecular\nperturbations in tumors that disrupt critical biological functions (e.g., microtubule\ndynamics, DNA replication), ultimately leading to tumor growth inhibition or cell death.\nIn our study, the term \"interaction\" specifically refers to these drug-induced perturbation\nprocesses that mediate therapeutic efficacy. Our neural network architecture models these\nperturbations through its hidden layers, thereby enabling treatment response prediction.\nWe have carefully revised the relevant sections to clarify this concept (Line 17; Lines 62-\n63; Line 114).\n2. “We demonstrate that GDnet (area under the curve (AUROC) = 0.725) significantly\noutperforms the transcriptome-only model (AUROC = 0.683) in predicting treatment\nresponse.” Is the improvement in AUROC statistically significant? A statistical test\nshould be performed to justify this claim.\nResponse: Thank you for your valuable suggestion. We performed DeLong's test, which\nyielded a significant result (p = 0.006) as shown in Figure 3c and Line 25.\n3. What does the term “pan-cancer cell line” in Figure 1 mean? There is no explanation of\nthis term in the text. Please clarify.\nResponse: Thank you for your careful review. We apologize for not providing this\nexplanation in the original manuscript. The term \"pan-cancer cell lines\" refers to the 800\ncell lines from the GDSC and CTRP datasets, which include multiple cancer types such\nas breast, lung, colorectal, and pancreatic cancers - hence their designation as pan-cancer.\nWe have now added this clarification in Figure 1a’s legend and the corresponding text\nsection (Line 99).\n4. An ablation study of the network structure should be conducted to demonstrate the\nimportance of each component of the model.\nResponse: We sincerely appreciate your suggestion. The proposed GDnet model\ncomprises three critical modules: transcriptomic module, drug representation module, and\nprior knowledge module. To demonstrate the significance of these modules, we designed\ncorresponding ablation experiments: (1) For the transcriptomic module: We developed a\ndrug-representation-only model (Dnet) to validate the importance of transcriptomic data.\n(2) For the drug representation module: We used the transcriptome-only model (Gnet)\nand alternative drug representation-based models (GSnet and GTnet) to verify the\nsuperiority of our proposed drug representation approach. (3) For the prior knowledge\nmodule: We used GDdnet, a fully connected feedforward neural network model without\nprior knowledge integration. The results showed that removing any component impaired\nmodel performance (Supplementary Fig. 6), conclusively demonstrating that tumor\ntranscriptome profiles, optimized drug representations, and biological prior knowledge\nintegration are all indispensable for GDnet's high performance. We have incorporated\nthese analyses in Supplementary Fig. 6 and detailed them in both the Results (Lines 177-\n178) and Methods (Lines 616-624) sections.\n5. Typo corrections:\n• In Figure 1: “Gnome-wide” should be corrected to “Genome-wide.”\nResponse: Thank you for your careful review. We have made the corresponding\ncorrection (Figure 1a).\n• On page 13: “malignanc” should be corrected to “malignancy.”\nResponse: Thank you for your careful review. We apologize for the lack of detail in the\narticle regarding this part. Here, we used \"*\" as a wildcard, meaning \"malignanc*\"\nincludes all expressions with the prefix \"malignanc,\" such as \"malignancy,\"\n\"malignancies,\" etc. The scope is thus broader than the single term \"malignancy.\"\nWe have added a detailed explanation in the Methods section of the article (Lines 427-429).\n6. The second paragraph in the section “GDnet showed superiority in predicting the\noutcome of neoadjuvant therapy” should be rewritten. The illustrations of the results\nshould be reorganized in a more coherent and logical way. Currently, too many methods\nand results are listed without clear organization.\nResponse: Thank you for your suggestion. We have reorganized this section, presenting\nthe results sequentially as follows: the training and testing dataset configuration (Lines\n145-153, Supplementary Fig. 2), the benchmarking comparison design of GDnet (Lines\n154-164), internal 10-fold cross-validation results (Lines 165-168, Figure 3a,\nSupplementary Fig. 5), external test results from 10 independent training repetitions\n(Lines 169-178, Figure 3b, Supplementary Fig. 6), external test results after ensemble of\nthe 10 training repetitions (Lines 179-186, Figure 3c, Supplementary Fig. 7), subgroup\nanalysis of external test including molecular subtypes of breast cancer (Lines 186-188,\nFigure 3d, Supplementary Fig. 8) and 13 individual external test datasets (Lines 188-198,\nSupplementary Figs. 9-10), and prospective test results (Lines 198-202, Figure 3e-f,\nSupplementary Fig. 11).\n7. How were the numbers of nodes in each layer, such as the gene layer and pathway\nlayers, determined? Please provide justification or additional details.\n8. How were the connections between nodes in adjacent layers determined? Clarify the\nrationale or methodology behind these connections.\n9. The biological meanings of “layer1,” “layer2,” and “layer3” are not clear. Additional\nexplanations should be provided.\nResponse: We address these three model structure-related questions collectively. First,\nwe appreciate your scrutiny of GDNet's architecture and apologize for the insufficient\nclarity in our original description.\nThe constructed layered network originated from 8,377 genes that are available both in\nour patient transcriptomic datasets and the Reactome dataset. These 8,377 genes\nconstructed the first hidden layer (gene layer).\nThen we selected all the pathways that contain at least one of the 8,377 genes in the first\nlayer as the second hidden layer. As a result, 846 pathways were selected, thus the\nnumber of nodes in the second layer was determined as 846. Specifically, the affiliation\nrelationships between genes and pathways (the strings starting with ‘R-HSA’) or between\npathways and pathways were obtained from the Reactome dataset and listed as below:\nNext, the affiliation relationships between 846 pathways and 8,377 genes are reorganized\nas a 2D binary sparse matrix with genes as rows and pathways as columns as displayed\nbelow. 0 means no affiliation relationship between a gene and a pathway, and 1 means the\ngene is a member of the pathway. The 2D matrix served as a mask between the first\nhidden layer and the second hidden layer, multiplied by weights after each epoch during\nweight updates, imposing constraints on connections between the two layers. This is the\napproach we determined the connections between nodes in adjacent layers.\nSubsequently, pathways directly related to those in the second layer constitute the third\nlayer (218 superior pathways), and this layering process continues analogously to the\nfourth layer (108 pathways), the fifth layer (51 pathways), until the sixth layers with 18\ntop pathways representing broad biological processes are built. The 2D matrixes defining\nthe connections between nodes in adjacent layers are also established. We have added a\ndetailed explanation in the Methods section of the article (Lines 527-562). The detailed\ncodes are uploaded to CodeOcean (the\n‘2_Masks_construction_from_biologically_knowledge.ipynb’ file).\nRegarding the biological interpretation of each hidden layer: Hidden layer 1 represents\ngenes, layer 2 corresponds to primary signaling pathways with singular functions, while\nlayers 3 through 6 progressively represent higher-order pathways with increasingly\ncomplex functionalities. Specifically, layer 6 encompasses the most broadly functional\npathways, such as \"immune system\" and \"cell cycle\" (as illustrated in Figure 5b). Much\nlike the hierarchical taxonomy in the Linnaean classification system (Kingdom →\nPhylum → Class → Order → Family → Genus → Species), GDNet's hidden layers 1\nthrough 6 exhibit a reverse progressive relationship - evolving from specific biological\nelements in the initial layers to more comprehensive biological functions in deeper layers\n(with node numbers gradually decreasing). This architecture ultimately predicts patient\ntreatment response, which itself represents the most complex biological event in this task.\nThis design philosophy aligns with fundamental deep learning principles where early\nlayers capture low-level features while subsequent layers progressively integrate them\ninto higher-level abstractions.\n10. The potential applications of the proposed method beyond breast cancer should be\nelaborated upon. More discussion is needed on how this method can support more\ninformed decision-making in selecting neoadjuvant treatments.\nResponse: We appreciate your suggestions and have expanded the discussion in our\nmanuscript to address both the potential applications beyond breast cancer neoadjuvant\ntherapy (Lines 406-414) and how this method can support more informed decision-\nmaking in selecting neoadjuvant treatments (Lines 365-385).\nRegarding the potential application beyond breast cancer neoadjuvant therapy (Lines\n406-414), we believe this framework can be extended to a broader range of diseases and\ntherapeutic scenarios. While this study primarily focuses on the neoadjuvant treatment\nsetting for breast cancer, the same framework can be adapted to other treatment stages\n(e.g., adjuvant therapy, advanced-stage therapy) and other cancer types (e.g., lung cancer,\npancreatic cancer, prostate cancer) by constructing corresponding digital drug-testing\nsurrogate models based on their respective efficacy evaluation systems. Moreover, the\nframework can also be tailored to non-oncological diseases, such as diabetes and\npsychiatric disorders, by developing disease-specific drug representation methods that\nalign with their unique omics profiles and building corresponding models for treatment\nselection.\nConcerning how this method can support more informed decision-making in selecting\nneoadjuvant treatments, we listed three potential scenarios in which GDnet can play an\nimportant role (Lines 365-385): (1) In real-world clinical practice, for a breast cancer\npatient scheduled for neoadjuvant therapy, we can first list all eligible treatment options\nbased on guidelines and clinical trial evidence. The patient’s tumor transcriptome data\nand treatment options can then be input into GDnet to predict the probability of achieving\na pCR. The regimen with the highest pCR probability can then be recommended as the\noptimal treatment choice for the patient. The efficacy of this treatment selection strategy\nis demonstrated by in-silico simulation trials based on ISPY-2 and external test dataset,\nshowing significantly improved pCR rate compared to the control group in our study. (2)\nIn clinical trial design, GDnet can assist in optimizing treatment strategies. Specifically,\ngiven a set of candidate monotherapies or combination regimens, GDnet can simulate and\npredict the pCR probabilities for each regimen option in patient subpopulations that meet\nthe inclusion criteria. The regimen demonstrating the highest overall pCR probability\nacross the patient subpopulations would be the most promising candidate for clinical trial.\nEven for novel drug combinations not included in GDnet’s training set, the model can\nstill guide trial strategy selection. This is evidenced by GDnet’s strong performance in the\nexternal test set, where it successfully predicted the efficacy of drugs that were not\nincluded in model training (e.g., durvalumab, bevacizumab, cisplatin, olaparib, and\ngemcitabine). (3) Beyond trial design, GDnet can further optimize patient selection for clinical trials.\nBy estimating individual patients’ likelihood of responding to a trial regimen, GDnet can\nidentify the subgroup most likely to benefit, enabling selective enrollment and thereby\nimproving trial success rates.\nReviewer #5\nReviewer #5 (Remarks to the Author): Expert in machine learning, deep learning, cancer\nomics, and drug response\nSummary of the work\nIn this work, Ge et al. presented GDnet, a deep learning model integrating drug\nrepresentations with bulk transcriptomics to predict neoadjuvant therapy response in\nbreast cancer. The approach of representing each drug as a set of genes with\nsignificant correlation to the corresponding cell line response is a resourceful and\ninnovative approach that facilitates incorporating a large number of drugs into a single\nmodel effectively when using biologically informed DL. The validation of GDnet on a\nlarge number of patients highlights the potential of such approaches, with promising results on\noptimizing drug regimen selection. However, I have a couple of concerns that the authors\nshould address to further improve the manuscript.\nComments to the authors\nOverall, the manuscript is concise and easy-to-follow with sufficient details for the most\npart, albeit with a number of typos and image quality issues. Here are the specific\ncomments for further improvement:\nResponse: Thank you for your positive feedback on our work! We have carefully\naddressed each of your comments and made corresponding revisions to the manuscript,\nwhich have further improved the quality of the paper.\n1. The authors have compared their GDnet model against a number of different DL\nmodels incorporating different architectures or data modalities. However, as a baseline,\nthey should also compare their performance against a model based on known clinical\nfactors such as tumor grade, ER status, HER2 status etc.\nResponse: We appreciate your suggestion regarding the importance of comparison with\nbaseline models. After careful consideration, we identified three key aspects requiring\ndirect baseline comparisons. For clinical factors—which differ in modality from\ntranscriptomic data—we treated them as complementary (synergistic) factors rather than\ncompetitive benchmarks. Below, we outline the rationale for these four aspects:\n(1) Transcriptome-only baseline (Gnet): Since one of our model's key innovations\ninvolves incorporating treatment representation into transcriptomic data, the\ntranscriptome-only model (Gnet in our paper) serves as a primary baseline.\nDemonstrating GDnet's superiority over Gnet validates the importance of integrating\ndrug representation.\n(2) Alternative drug representation methods: As our drug representation strategy itself\nrepresents an innovation, we included additional control models (GSNet and GTNet)\nemploying previously published drug representation methods to demonstrate the\nadvantages of our novel approach.\n(3) Traditional machine learning models: Given our use of deep learning, we included\ncomparisons with conventional machine learning approaches to highlight the benefits of\nour architecture.\n(4) Clinical factors (synergy factors): In current clinical practice where treatment\ndecisions are primarily guided by clinical factors, we need to clarify the relationship\nbetween GDnet and these well-established clinical parameters. Fundamentally, we view\nthis relationship as complementary (demonstrating that GDnet's functionality does not\nentirely overlap with clinical factors and serves to supplement them) rather than\ncompetitive (where GDnet's performance surpassing clinical factors leads to their\nsubstitution). Specifically, existing clinical factors (demonstrating robust predictive\nperformance with AUROC=0.7 for neoadjuvant therapy response in Sammut et al.'s\nstudy1) are evidence-based and clinically validated, however, integrating other modality\nfeatures like transcriptomic, genomic, and pathologic data can further enhance predictive\naccuracy1. This establishes GDnet and clinical factors as complementary dimensions of\npatient characterization, rather than competitive and replaceable information.\nTo empirically demonstrate this complementarity, we can adopt the following two\napproaches: (1) developing a combined clinical-GDnet multimodal model to show\nsuperior performance over clinical-factor-only model, and (2) conducting clinical-factor-\nstratified analyses to prove GDnet's consistent performance across subgroups. Through\nsystematic clinical data collection of included datasets (Supplementary Table 2), we\nfound only ER/HER2 status was relatively complete across datasets, while other factors\nlike tumor stage/grade were largely unavailable, precluding comprehensive multimodal\nmodeling. Focusing on the available ER/HER2 status, our stratified analyses confirmed\nGDnet maintains strong predictive performance across all subtypes (Figure 3d;\nSupplementary Fig. 8), demonstrating its complementary value to these clinical factors.\nIn summary, we conducted a comprehensive evaluation of baseline models, including a\ntranscriptome-only model, two models based on previous drug representation strategies,\nand simple machine learning models, to demonstrate the superior performance of GDnet.\nRegarding clinical factors, since the datasets we included only had relatively complete\nmolecular subtyping data, we performed a subgroup analysis based on molecular\nsubtypes, confirming GDnet’s complementary value to these clinical factors.\n2. Breast cancer has multiple established subtypes – clinical (based on ER/HER2 status)\nand molecular (based on PAM50 assay) with implications to patient response and\nsurvival. The authors should perform a subtype-specific to evaluate whether GDnet\nperforms similar across different subtypes.\nResponse: Thank you for your suggestion, this point is indeed crucial for the clinical\napplication of GDnet in breast cancer. As mentioned above, we collected ER and HER2\nstatus data from the included datasets (Supplementary Table 2) and stratified patients into\nfive groups: ER+HER2+, ER-HER2+, ER+HER2-, TNBC, and unknown (incomplete ER\nor HER2 data). The results demonstrate that GDnet maintains strong predictive\nperformance across all molecular subtypes (Figure 3d; Supplementary Fig. 8; all\nAUROC > 0.65). Furthermore, GDnet outperforms the transcriptome-only model (Gnet)\nin all subtypes except TNBC, where its performance is only slightly lower than Gnet\n(Figure 3d; Supplementary Fig. 8). This highlights the superiority of GDnet's modeling\nstrategy over transcriptome-only approaches.\n3. There are a couple of missing technical info from the main text:\na. The authors mentioned using the genome-wide gene expression correlation with IC50\nto represent different drugs. Is it Pearson or Spearman correlation?\nb. They also mentioned applying rank normalization on the transcriptomics data. is this\nnormalization in a descending or ascending order?\nResponse: Thank you for your careful review, and we apologize for omitting these\ndetails in the manuscript.\nFor the correlation analysis, since the normalized gene expression data and log-\ntransformed IC50 values generally followed a normal distribution, we used Pearson\ncorrelation for the calculations. These details have now been added to the manuscript\n(Lines 471–472).\nRegarding rank normalization, we employed ascending order, ensuring that the ranked\norder preserved the original expression-level hierarchy of the genes. Specifically, we\nsorted each gene’s original expression levels, assigning a value of 0 to the lowest\nexpression level and 1 to the highest, with intermediate values scaled proportionally\nbetween 0 and 1. This means the rank-encoded order exactly matches the pre-ranking\nexpression-level order. These methodological details have also been incorporated into the\nmanuscript (Lines 455–458).\n4. In the performance evaluation section (line 179-180), the authors mentioned “Though\nGDnet showed lower AUROCs in the left 5 datasets, all the differences were less than\n0.03”. Looking at Figure 3- it is not readily clear what they are referring to, especially\nwith the difference being <0.03. The authors should revise it to clearly state what the\nintended message was.\nResponse: Thank you for your suggestion. This statement described the content from\nFigure 3d in our previous manuscript, and we sincerely apologize for not having directly\nvisualized these results in the original figure. In our revised manuscript, we have\nreorganized Figure 3's content, with the original Figure 3d now relocated to\nSupplementary Fig. 9a. For improved clarity, we've annotated all bars in the barplots\n(including Supplementary Fig. 9a) with their respective values to facilitate observation of\nAUROC value differences. Furthermore, due to corrections made to previous bugs in\nAUROC calculations, these results have been updated. The manuscript now states (Lines\n192-195): \"The results showed that GDnet outperformed Gnet, achieving higher\nAUROCs in 10 out of 13 datasets, while exhibiting only marginally lower AUROCs\n(difference < 0.02) in the remaining three datasets (Supplementary Fig. 9a).\" Specifically,\nGDnet showed slightly lower AUROC values than Gnet in GSE21997 (0.796 vs 0.811,\ndifference = 0.015), GSE42822 (0.78 vs 0.784, difference = 0.004), and GSE23988 (0.94\nvs 0.945, difference = 0.005), though all differences remained below 0.02.\nRegarding the aforementioned bugs in previous AUROC calculations, the details are as\nfollows:\nPreviously, for ease of plotting in R, we exported all model prediction values and\ncalculated AUROC/AUPRC values using pROC package in R. In this revision, we\nstandardized our analysis by exclusively using Python's sklearn.metrics functions. This\nrevealed discrepancies between the two calculation methods. Upon investigation, we\nconfirmed this was caused by a bug in the pROC package.\nFor example, in GSE8465:\npROC calculation: AUROC=0.55 (with aberrant ROC curve shape shown below)\nsklearn.metrics calculation: AUROC=0.45 (with proper ROC curve shown below)\nFurther exploration of GDnet prediction values showed (below Figure) pCR group had a\nlower median prediction value than the RD group (0.038 vs 0.05); this typically\ncorresponding to AUROC<0.5.\nThe aberrant ROC curve from pROC and inverted AUROC values suggest the package\nmay flip results (calculating 1-correct AUROC) under specific data distributions. We\nhave therefore updated all AUROC/AUPRC values using Python.\nRegarding the observed AUROC values below 0.5 in GSE8465 and GSE122630, our\ninvestigation of these datasets revealed no special characteristics. Considering that GDnet\nstill outperformed Gnet in the two datasets, and only these two small sample datasets\nshowed AUROC < 0.5, while the model performed well in other external test datasets and\nprospective test datasets, we believe it may be an incidental result due to the small sample\nsize.\n5. The interpretability analysis identified a few genes with high importance in response\nprediction. I was wondering if the authors could add the roles/biological significance of\nthese genes in breast cancer (or cancer in general) from relevant literature.\nResponse: Thank you for your suggestion, this can significantly improve the clarity of\nour model interpretation. After reviewing relevant literature, we have annotated each\ngene's association with breast cancer (Lines 272-286).\nNotably, while our previous model interpretation examined all treatment strategies in the\ntest set, we recognized that molecular mechanisms of response or resistance are\ninherently treatment specific. We therefore refined our interpretation approach by\nfocusing on the most promising treatment regimens from the simulated trial optimization\ngroup: the dual anti-HER2 regimen\n(paclitaxel|doxorubicin|cyclophosphamide|ERBB2|ERBB2) and the immunotherapy\nregimen (paclitaxel|doxorubicin|cyclophosphamide|olaparib|PDCD1), conducting\nseparate model interpretations for the biological entities (genes and pathways) in each\ncase.\nFor the dual anti-HER2 regimen, gene-level analysis identified PIK3CA, PSMB3,\nYWHAH, PSMB7, and SRC as important for treatment response (Figure 5b). For the\nimmunotherapy regimen, key response-associated genes included TRIO, JUN, DUSP6,\nCCND1, and SRC (Supplementary Fig. 15). Brief breast cancer-related functional\ndescriptions of these genes have been incorporated into the manuscript (Lines 272-286).\n6. Figure 1, Supplementary Figure 1-2, 8-10 are too blurry and low-quality to properly\nexamine. The authors should replace these figures with high quality versions.\nResponse: Thank you for your careful review. We have improved the resolution of the\nfigures accordingly.\n7. There are multiple typos across the manuscript that the authors should rectify, e.g.,\na. On page 5, line 116, they mentioned HA et al. which should be Elmarakeby et al.\nResponse: Thank you for your careful review. We have incorporated this change into the\nrevised manuscript (Line 122).\nb. On page 9, line 244-248, the authors mentioned different pathways in different layers.\nHowever, the numbering of the layers is not consistent e.g., “ER-Phagosome pathway in\nlayer 1” while it is actually on layer 2 of the network. They should resolve these\ndiscrepancies.\nResponse: We sincerely appreciate you pointing out this issue, and we apologize for the\nambiguity in our previous description. Originally, we referred to the first hidden layer as\nthe \"gene layer\" and the second hidden layer as the \"first pathway layer,\" which was then\nabbreviated as \"layer 1.\" We acknowledge this naming approach was indeed confusing.\nTo address this, we have standardized our terminology to consistently use \"hidden layer\"\ndesignations throughout the manuscript, changing previous \"layer 1\" to current \"hidden\nlayer 2,\" (Lines 277-281) and have implemented these revisions across the entire text.\nc. On page 10, line 280, they mentioned RNA-sequencing methods which should be\ntranscriptomics technologies (to include microarray techniques)\nResponse: Thank you for your careful review. We have made the corresponding revision\nin the revised manuscript (Line 300).\nd. On page 11, line 298, they mentioned “the 21-gene recurrence score”, which should be\nOncotypeDX score.\nResponse: We sincerely appreciate your expert feedback and apologize for our incorrect\nphrasing. We have made the corresponding corrections in the revised manuscript (Line\n318)\nReviewer #5 (Remarks on code availability):\nDid a quick glance on the code. It seems to be in order and well organized.\nReviewer #6\nReviewer #6 (Remarks to the Author): Expert in machine learning, deep learning, cancer\nomics, and drug response\nThis manuscript presents GDnet, a biologically-informed deep learning approach that\nintegrates drug representations with transcriptomic data to predict and optimise breast\ncancer treatment response. The work makes several notable contributions: assembly of a\nlarge dataset (4,371 patients across 31 datasets), an innovative approach to drug\nrepresentation using genome-wide expression correlations, and demonstration of\npotential clinical utility through in-silico trials. The study is methodologically rigorous\nwith detailed model specifications, thorough ablation studies, and proper validation\nprocedures.\nHowever, several aspects require attention. The benchmark comparison needs substantial\nexpansion to include existing approaches from the literature. Some analysis results show\ncounterintuitive patterns that require explanation, such as the unexpected relationship\nbetween optimisation thresholds and pCR rates, and the surprisingly better performance\non test data compared to training data. The manuscript would also benefit from improved\nclarity in methodology description and figure quality. With these improvements,\nparticularly in comparative evaluation and results analysis, this work could make a\nsignificant contribution to precision oncology.\nResponse: We sincerely appreciate your valuable suggestions! In response, we have (1)\nadditionally included benchmark comparisons for GDnet against two models based on\nprior drug representation strategies, dense models with optimized nodes and layers\nhyperparameters, and simple machine learning models; (2) provided explanations and\ndiscussions regarding the counterintuitive patterns in results; (3) improved the\nmethodological descriptions; and (4) enhanced the quality of figures. Your insightful\ncomments have not only significantly improved the quality of our manuscript but have\nalso deepened our understanding of the study results.\nMajor:\nThe paper's core contribution needs better clarification and stronger empirical validation.\nWhile the authors have assembled an impressively large dataset (4,371 patients across 31\ndatasets), this alone isn't sufficient as a key contribution. The proposed GDnet model,\nwhich aims to integrate drug representations with transcriptomic data for treatment\nresponse prediction, requires more rigorous validation against existing approaches.\nThe authors describe three main approaches to drug characterisation in previous work:\n1.Structure-based drug representation, 2.Target-based drug representation, 3.Genome-\nwide correlation-based representation (their approach).However, simply stating that\nprevious approaches are unsuitable for this specific task isn't convincing without direct\nempirical comparison. The benchmark section should include implementations of at least\n2-3 representative methods from these different approaches. While the extensive\ncomparison of GDnet variants is valuable for ablation studies, it doesn't demonstrate\nadvancement over the state-of-the-art.\nThe paper lacks crucial details and justification regarding the architecture of the dense\nneural network baselines. While the layer sizes for the biologically-informed models are\nwell documented, the corresponding information for the dense models has to be inferred\nfrom the code. It appears the authors simply matched the layer sizes of their biologically-\ninformed models, but this may not be optimal given the sample size and task complexity.\nFor a fair comparison, the dense neural networks should be evaluated with different\narchitectural configurations - varying both the number of layers and nodes per layer -\nalong with proper hyperparameter optimisation.\nFurthermore, the comparison would be more complete by including traditional machine\nlearning approaches that often perform well on tasks with this scale of data. Models like\nlogistic regression, random forest, and XGBoost could serve as strong baselines,\nespecially since they tend to be more sample-efficient than deep learning models. These\nsimpler models might actually perform competitively given the relatively modest dataset\nsize of 4,371 samples, and their inclusion would provide better context for evaluating the\nnecessity of a deep learning approach. This additional benchmarking would strengthen\nthe authors' claims about the advantages of their biologically-informed architecture.\nResponse: We appreciate your three suggestions regarding the benchmark comparison\nsetup for GDnet, and we address them collectively in this response. We fully agree with\nyour proposed comparative models and have implemented the following additions to our\nanalysis:\n(1) Models using drug representation strategies in previous studies: GSnet and GTnet\nPrevious drug representation strategies were primarily developed for cancer cell\nmodeling rather than patient-level applications. To our knowledge, no existing models for\nbreast cancer neoadjuvant therapy (except ours) have incorporated drug representations.\nTherefore, we established two additional models based on previous drug representation\nstrategies for comparison with GDnet (Lines 563-580), including hyperparameter tuning\n(Lines 598-609):\n GSnet (Gene-SMILES net): we use the drug representation strategy based on\nSMILES (Simplified Molecular-Input Line-Entry System) molecule structure\nexpressions proposed by Kuenzi et al.10. In this approach, drugs are hashed into\n2048-dimensional vectors to represent molecular drugs. Antibody-based drug\nrepresentations were not generated due to the lack of SMILES entries. Regimens\nwere represented as the vector summation of involved drugs. The SMILES-based\n2048-dimensional vectors and gene expression vectors were subsequently used as\ninput for GSnet.\n GTnet (Gene-Target net): we use the drug representation strategy based on the\ntarget of drugs proposed by Zhang et al.13. Drug target information was obtained\nfrom the DrugBank database14. Chemotherapeutic agents with nonspecific targets\n(e.g., DNA, RNA) that cannot be mapped to specific genes were excluded from\ndrug representation. This approach incorporated the following 10 drugs with\nclearly mappable gene targets: trastuzumab (ERBB2), pertuzumab (ERBB2),\nbevacizumab (VEGFA), AMG-386 (ANGPT1, ANGPT2), MK-2206 (AKT1),\nolaparib (PARP1, PARP2), durvalumab (PDCD1), pembrolizumab (PDCD1),\nganitumab (IGF1R), luminespib (HSP90AA1). Regimens were represented as the\nvector summation of involved drugs. The target-based 10-dimensional vectors and\ngene expression vectors were subsequently used as input for GTnet.\n(2) Dense models with optimized nodes/layers: GDdnet-o, GAdnet-o, Gdnetens-o, and\nGdnet-o\nWe revised fixed-architecture dense model nomenclature (GDNNet→GDdnet,\nGANNet→GAdnet, GNNetens→Gdnetens, GNNet→Gdnet). While maintaining the\ndense models with identical node and layer configurations as our informed models\n(shown in the \"dense\" group of Figure 3a), we have further introduced optimized versions\nwith tuned hyperparameters for nodes and layers. These optimized models are denoted by\nthe \"-o\" suffix (presented in the \"optimized\" group of Figure 3a). For hyperparameter\noptimization, we explored hidden layers ranging from 3 to 6, with node configurations\nfollowing geometric decay of around 2 or 4 folds (Lines 598-609).\n(3) Machine learning models\nWe trained and hyperparameter-tuned four models: L2-regularized logistic regression,\nrandom forest, SVM, and XGBoost, using optimal parameters for final training (Lines\n610-615).\nThe results demonstrate that GDnet consistently outperformed all comparator models in\nthe external test set (Figure 3b). Notably, in our newly introduced prospective test set\ncomprising 134 cases, GDnet maintained superior performance over the majority of\ncomparison models (Figure 3e). These findings collectively validate the advantages of\nGDnet's modeling strategy.\nWhile the authors emphasise how biological knowledge improved performance in their\nGDnet architecture, they don't address why incorporating the same biological knowledge\nactually decreased performance in all other model variants (GAnet, Gnetens, and Gnet\ncompared to their respective dense versions). This inconsistency raises questions about\nthe general effectiveness of their biological knowledge integration approach and needs to\nbe thoroughly discussed. Are there specific aspects of the GDnet architecture that make it\nuniquely suitable for biological knowledge integration? Or could this pattern suggest\npotential issues with the current implementation of biological knowledge in the other\narchitectures? Without addressing these questions, the claimed benefits of biological\nknowledge integration remain uncertain.\nResponse: We appreciate your critical question; this is an excellent point for discussion.\nWe will analyze this issue by collectively discussing the informed, dense, and optimized\ngroup models.\nFor GDnet, the incorporation of prior knowledge resulted in a total of 37,127 trainable\nparameters, whereas GDdnet (with node architecture [19300, 8377, 846, 218, 108, 51, 18,\n2]) and GDdnet-o (nodes: [19300, 8192, 4096, 2048, 2]) contained 168,977,476 and\n200,052,736 parameters, respectively. The results showed that GDnet demonstrated\nsignificantly superior performance over both dense and optimized variants (Figure 3a).\nFor GAnet, Gnet, and Gnetens, the informed models contained 20,373, 20,373, and\n101,865 parameters, respectively. Their dense variants—GAdnet (with node structure\n[8377, 846, 218, 108, 51, 18, 2]), Gdnet (with node structure [8377, 846, 218, 108, 51,\n18, 2]), and Gdnetens (5-model ensemble, each with node structure [8377, 846, 218, 108,\n51, 18, 2])—utilized 7,301,376, 7,301,376, and 36,506,880 parameters, while the\noptimized versions—GAnet-o (with node structure [8377, 2048, 512, 128, 2]), Gnet-o\n(with node structure [8377, 2048, 512, 128, 2]), Gnetens-o (5-model ensemble, each with\nnode structure [8377, 2048, 512, 128, 2])—contained 18,270,464, 18,270,464, and\n91,352,320 parameters, respectively. Here, the informed models exhibited slightly lower\nperformance than dense models but marginally outperformed optimized models (Figure\n3a).\nSome possible explanations are as follows: The prior knowledge in this study was\nderived from biological knowledge database, yet current research on genes, pathways,\nfunctions, etc., remains largely incomplete. Thus, existing biological knowledge\ndatabases do not fully reflect all actual biological behaviors. Assuming a fixed sample\nsize, in cases with an extremely large number of parameters, prior knowledge represents a\nguide for the identification of the correct and generalizable learning direction for deep\nlearning models under the limited sample size, thereby achieving better external test\nperformance compared to models not incorporating prior knowledge. However, when the\nnumber of parameters becomes smaller (while still sufficient for learning), deep learning\ncan accurately reveal more potential actual biological functions beyond the existing prior\nknowledge under the same sample size, leading to high external test performance. Thus,\nthere is a trade-off to incorporating prior knowledge. In the scenario of GAnet, Gnet,\nand Gnetens, the imperfect prior knowledge instead becomes a constraint in the training of\ndeep learning models. Specifically, prior knowledge integration improved GAnet, Gnet,\nand Gnetens performance over parameter-heavy optimized models (GAnet-o, Gnet-o,\nGnetens-o), yet these still underperformed relative to the relatively parameter-efficient\ndense models (GAnet, Gnet, Gnetens). In other words, dense models and optimized\nmodels represent the opposite sides of the balanced benefit-drawback point of\nincorporating prior knowledge (informed models) for GAnet, Gnet, and Gnetens. This\nphenomenon highlights the importance of selecting appropriate scenarios for applying\nprior knowledge. When the parameter count is large and the sample size is small, prior\nknowledge can significantly aid deep learning model training. However, this advantage\ngradually diminishes as the number of parameters decreases or the sample size increases.\nWe have incorporated this important discussion into our manuscript (Lines 329-357), as it\nprovides valuable insights for incorporating prior knowledge into deep learning models.\nIn the PSM analysis with Fig 4, the odds ratios for optimised vs control groups show a\nstronger linear trend (R² = 0.94, p < 0.001) in the external validation compared to the\ntraining data (R² = 0.87, p < 0.001). The pCR rates also appear to show more significant\ndifferences in the validation set results. This pattern contradicts standard machine\nlearning principles - a model should typically perform better on data it was trained on\ncompared to unseen test data, unless there's underfitting. However, the overall\nperformance metrics shown in Figure 3 don't suggest underfitting, as the model achieves\nreasonable AUROC values. The authors need to investigate and explain this unusual\npattern, as it raises questions about the robustness of their validation approach and the\nmodel's generalisation properties.\nResponse: Thank you for your question; it’s an excellent point. This part of the analytical\nprocess is original to our study, and your question has greatly enhanced our\nunderstanding of it. The issue is relatively complex, so please allow us to address it step\nby step.\nFirst, let’s examine the simulation trials based on ISPY-2 and external validation before\nPSM. For GDnet, the AUROC was 0.901 and AUPRC was 0.88 in ISPY-2, while in\nexternal test, the AUROC was 0.725 and AUPRC was 0.707. From the trial results, we\ncan see that in the ISPY-2 based trials, the difference in pCR between the optimized\n(opt) and control (ctl) groups was substantial, with more significant p-values and a higher\noverall OR, increasing from 3 to >7 (below Figure; Supplementary Fig. 12c-d).\nIn contrast, in the external test-based trials, the difference in pCR between the optimized\nand control groups was smaller, with less significant p-values and a lower overall OR,\nincreasing from 1.3 to 3 (below Figure; Supplementary Fig. 13c-d).\nWe believe this aligns with the results that GDnet performs better in the ISPY-2\ndataset (one of the internal validation sets) but shows lower performance in the external\ntest set. \nSecond, let’s look at the simulation trials after Gnet-based PSM. For Gnet in ISPY-2,\nthe AUROC was 0.949 and AUPRC was 0.934, while in external test, the AUROC was\n0.683 and AUPRC was 0.639. The trial results show that, compared to pre-PSM, the\ndifference in pCR between the optimized and control groups in ISPY-2 based trials\ndecreased a lot, and the overall OR also decreased significantly to less than 2.5 (below\nFigure; Figure 4b-c).\nIn the external test-based trials, however, the difference in pCR between the optimized\nand control groups changed little compared to pre-PSM, and the overall OR did not\ndecrease noticeably (below Figure; Figure 4d-e).\nWe observed that after PSM, the OR in the ISPY-2 based trials decreased substantially,\nwhile the OR in the external test set based group showed no significant reduction. We\nattribute this phenomenon to two main reasons:\n The Generally Better Treatment Regimens in ISPY-2: All ISPY-2 regimens\nconsisted of at least three-drug combinations with better efficacy (Supplementary\nTable 2, GSE194040). When tumor malignancy was balanced between the\noptimized and control groups, the theoretical difference in treatment efficacy was\ninherently smaller, resulting in a lower OR. Thus, after PSM using Gnet to\nbalance tumor malignancy, the reduction in OR between the opt and ctl groups\ncan be explained. In contrast, the external test set included a wide range of\ntreatments, from single-drug (less effective) to dual-drug, triple-drug, and more\ndrug combinations (Supplementary Table 2). After PSM balanced tumor\nmalignancy, the higher OR in the external test set can be reasonable.\n Overfitting of Gnet in ISPY-2: Gnet achieved an AUROC of 0.949 and AUPRC of\n0.934 in the ISPY-2 dataset, indicating higher overfitting compared to GDnet.\nConsequently, after PSM using Gnet, the predicted pCR rates in the opt and ctl\ngroups became very similar in ISPY-2 based trials, leaving little space for GDnet\nto demonstrate its effect, which explains the reduced OR.\nAs for why the p-value significances for pCR differences decreased overall in the ISPY-2\nbased trials after PSM, aside from the two factors mentioned above, we believe sample\nsize also plays an important role. As shown in the below figure, after PSM, the sample\nsize in the ISPY-2 based trials (below lest; Supplementary Fig. 12e) was significantly\nsmaller than in the external test-based trials (below right; Supplementary Fig. 13e). This\nsmaller sample size would lead to less significant p-values for pCR differences in the\nISPY-2 based trials.\nFinally, regarding why the R² of the OR fitting curve in the ISPY-2 based trials after PSM\nwas lower than in the external test-based trials, we believe this can also be explained by\nthe three factors mentioned above: (1) the generally better treatment regimens in ISPY-2,\n(2) the higher overfitting of Gnet compared to GDnet in ISPY-2, and (3) the smaller\nsample size in the ISPY-2 group after PSM.\n Treatment Regimen Distribution: Assuming tumor malignancy is balanced across\noptimization and control groups, R² is primarily influenced by the treatment\nregimen distribution in the dataset. The external test set includes a balanced\ndistribution of treatments ranging from single drug (less effective) to multi-drug\ncombinations (more effective) (Supplementary Table 2), creating a linear gradient\nfrom poor to good regimens. Under such a distribution, the expected OR values\nunder different optimization intensities would be more likely to follow a linear\ngradient, resulting in an R² closer to 1. In contrast, ISPY-2 trials exclusively used\nat least triple-drug combinations—all relatively effective—meaning there was no\nclear linear gradient from poor to good regimens. Thus, the expected OR values\nunder different optimization intensities would not like to show a strong linear\ntrend, leading to an R² farther from 1.\n Overfitting of Gnet: As mentioned earlier, Gnet exhibited higher overfitting in the\nISPY-2 dataset, leaving little space for GDnet to demonstrate its effect. This not\nonly reduced the OR but also made the assessed OR values less stable, pushing R²\nfarther from 1.\n Sample Size Reduction After PSM: The smaller sample size in the ISPY-2 based\ntrials after PSM made OR estimation less robust, further reducing R².\nIn summary, the lower R² of the OR fitting curve in the ISPY-2 based trials after PSM\ncompared to the external test-based trials is theoretically consistent.\nThe relationships between optimisation threshold and pCR rates shown in Figure 4\ndisplay some counterintuitive patterns that require explanation. Theoretically, we would\nexpect two clear trends: the optimised group should show decreasing pCR rates as the\nthreshold becomes less stringent (includes more cases), while the control group should\nshow increasing pCR rates. However, this expected pattern isn't consistently observed,\nparticularly in the I-SPY2 training data results (4b-c). Most concerningly, the control\ngroup shows an illogical pattern where the worst 20% of cases actually have a higher\npCR rate than the worst 50%. This makes little sense, especially in training data where\nthe model should have learned to rank cases more reliably. The authors should explain\nthis unexpected behaviour, as it raises questions about the robustness of their treatment\noptimisation approach.\nResponse: Thank you for your good question. We agree with this conclusion: the\noptimized group should show decreasing pCR rates as the threshold becomes less\nstringent (including more cases), while the control group should show increasing pCR\nrates. However, this conclusion should have a prerequisite: the malignancy level of\ntumors remains constant as the threshold changes. As shown in the figure below, in the\nISPY-2 based trials, we ensured balanced tumor malignancy between the optimized (opt)\nand control (ctl) groups through Gnet-based PSM, but the malignancy level (reflected by\nGnet score) actually varies across different optimization thresholds (Supplementary Fig.\n12f).\nThis is precisely why using the odds ratio (OR), which calculates the ratio between opt\nand ctl groups, provides a more reliable comparison of GDnet's optimization capability\nacross different thresholds. The rationale is as follows: as tumor malignancy\nincreases/decreases, the pCR rates in both opt and ctl groups will correspondingly\ndecrease/increase, while the ratio between their pCR rates remains relatively stable.\nRegarding why the pCR rate in the ctl group in Figure 4b-c is higher at 20% threshold\nthan at 50% threshold, our explanation is based on the theory that the final pCR rate is\njointly determined by tumor malignancy level and treatment optimization effectiveness -\nin other words, the final pCR rate represents the combined effect of these two factors. As\ndemonstrated in Supplementary Fig. 12f, where Gnet scores represent tumor malignancy\n(with lower Gnet scores indicating lower pCR rates and higher malignant tumors), we can\nobserve that tumors at the 20% threshold are relatively more benign, while those at the\n50% threshold are much more malignant. This explains why the pCR rate is higher at the\n20% threshold compared to the 50% threshold. To better illustrate this phenomenon, we\ncreated a conceptual diagram for the opt group based on the \"tumor Gnet + treatment\neffect = final pCR prediction\" theory. It shows that despite better treatment regimens at\n50% threshold, the higher tumor malignancy results in lower pCR rates compared to the\n20% threshold group where lower malignancy compensates for less optimized treatment.\nThis same theory can also explain the pCR rate changes in the opt group in Figure 4b,\nas well as the pCR changes in both opt and ctl groups in the external test-based trials shown\nin Figure 4d.\nThe visualisation could also be improved by including more threshold steps between the\ncurrent intervals. This would help readers better understand if these unexpected patterns\nare genuine trends or possibly artifacts of the chosen threshold points. A smoother curve\nwith more data points would provide a clearer picture of how the optimisation truly\nperforms across different stringency levels.\nResponse: We sincerely appreciate your thoughtful suggestion. As shown in\nSupplementary Table 2 (GSE_ID: GSE194040), the ISPY-2 cohort contains 11 distinct\ntreatment regimens. This means any optimization threshold interval smaller than 1/11\n(9.1%) would capture all potential treatment decision change points. Similarly, the\nexternal test set's 19 regimens (Supplementary Table 2) require threshold intervals below\n1/19 (5.3%) to encompass all decision-altering points.\nOur study employed a threshold interval of 1/20 (5%), which not only covers all decision\nchange points but provides additional resolution (that’s why Figure 4b shows identical\npCR rates and OR values between 50% and 45% threshold, and between 40% and 35%\nthreshold).\nIn summary, our selected 5% threshold interval can capture all potential treatment\ndecision change points for both simulation trials, yielding optimally smooth curves.\nMinor:\nThe description of \"gene and value pairs\" on lines 91-92 needs clarification. The concept\nis central to the method but isn't clearly explained.\nResponse: Thank you for pointing out the incomplete description of \"gene and value\npairs\" in our original manuscript. We have revised the relevant section in the manuscript\n(Lines 96-113) to provide a clearer explanation of the gene and value pairs.\nFigure 2 would be more informative if annotated with the number of gene-correlation\npairs for each drug.\nResponse: Thank you for your suggestion. We have color-coded the number of pairs and\nincluded corresponding annotations in parentheses after drug names (Figure 2).\n\"AUC\" should be labelled as \"AUROC\" to explicitly distinguish it from AUPRC.\nResponse: Thank you for your suggestion. We have revised the text accordingly, as this\nformulation is indeed more consistent and standardized.\nSince Figure 3c and 3d present results averaged from ten runs, confidence intervals\nshould be calculated and displayed to show the stability of these results.\nAll bar plots, including Figure 3d, should be annotated with sample sizes to help readers\nassess the reliability of results across different datasets.\nResponse: Thank you for your suggestion. We have calculated confidence intervals for\nall bars based on 1000-times bootstrapping (Lines 670-672) and have annotated them in\nthe bar plots (e.g., Figure 4b, Figure 4d, Figure 4e). Additionally, we have included the\nsample sizes in parentheses following each dataset for all performance evaluations\nconducted on external test subsets (e.g., Figure 4d, Supplementary Fig. 8-9).\nThe current Supplementary Figure 5a, which shows a more complete view than Figure\n3d, would be more appropriate as a main figure panel. The current arrangement buries\nimportant information in supplementary materials.\nResponse: Thank you for your suggestion. Indeed, the previous Supplementary Figure 5a\ncontained more substantial information and carried greater significance than Figure 3d.\nHowever, following our manuscript revisions - which incorporated both a new\nprospective test set (now presented in the new Figure 3e-f) and molecular subtype\nanalyses of breast cancer (the new Figure 3d) - we recognized several considerations\nnecessitating reorganization. First, the previous Figure 3d primarily represents subgroup\nanalyses of previous Figure 3c, resulting in functional overlap. Second, the inclusion of\nmultiple new comparator models had substantially expanded this figure's complexity.\nConsequently, we have relocated previous Figure 3d to the supplementary materials\n(Supplementary Figure 10), while specifically highlighting comparisons between GDnet\nand transcriptome-only Gnet models in another supplementary figure (Supplementary\nFigure 9).\nThere's an interesting pattern in Supplementary Figure 5 where GANNet outperforms\nGDnet on several datasets. What characteristics of these datasets might explain this\nperformance difference?\nResponse: Thank you for your question. While preparing our response, we identified\nseveral package-related bugs that have affected these results. We apologize for any\nconfusion caused and have thoroughly confirmed all recalculated metrics. The specific\nissues are as follows:\nPreviously, for ease of plotting in R, we exported all model prediction values and\ncalculated AUROC/AUPRC values using pROC package in R. In this revision, we\nstandardized our analysis by exclusively using Python's sklearn.metrics functions. This\nrevealed discrepancies between the two calculation methods. Upon investigation, we\nconfirmed this was caused by a bug in the pROC package.\nFor example, in GSE8465:\npROC calculation: AUROC=0.55 (with aberrant ROC curve shape shown below)\nsklearn.metrics calculation: AUROC=0.45 (with proper ROC curve shown below)\nFurther exploration of GDnet prediction values showed (below Figure) pCR group had a\nlower median prediction value than the RD group (0.038 vs 0.05); this typically\ncorresponding to AUROC<0.5.\nThe aberrant ROC curve from pROC and inverted AUROC values suggest the package\nmay flip results (calculating 1-correct AUROC) under specific data distributions. We\nhave therefore updated all AUROC/AUPRC values using Python.\nThe updated results (as shown in the figures below) indicate that GANNet (renamed as\nGAdnet) demonstrated the significant performance advantage over GDnet in the\nGSE207248 dataset. Upon investigation, we found that GSE207248 involved dual anti-\nHER2 therapy in HER2-positive inflammatory breast cancer cases. The superior\nperformance of GAdnet in this dataset may be attributed to either the unique\ncharacteristics of inflammatory breast cancer or potential random errors due to the\nrelatively small sample size. In datasets GSE21974 and GSE23988, GAdnet showed\nmarginally better performance than GDnet. GSE21974 represented chemotherapy across\nall breast cancer subtypes, while GSE23988 involved chemotherapy specifically in\nHER2-negative breast cancer patients. Both datasets employed conventional\nchemotherapy regimens with no distinctive patient population characteristics, suggesting\nthat the observed performance differences were likely due to random variations. In all\nother datasets examined, the performance between the two models was either comparable\nor favored GDnet, and thus no further investigation was conducted.\nRegarding the observed AUROC values below 0.5 in GSE8465 and GSE122630, our\ninvestigation of these datasets revealed no special characteristics. Considering that GDnet\nstill outperformed Gnet in the two datasets, and only these two small sample datasets\nshowed AUROC < 0.5, while the model performed well in other external test sets and\nprospective test sets, we believe it may be an incidental result due to the small sample\nsize.\nThe Figure 4 legend needs reorganisation and better descriptions to match the panel order\nand clearly explain the experimental setup.\nResponse: Thank you for your suggestion. We sincerely apologize for the lack of clarity\nin our previous description. We have now rephrased the text to present the figures in a\nmore sequential manner and have enhanced the experimental setup details (Lines 761-\n781, Figure 4).\nThere's a typo on line 221 where \"detests\" should presumably be \"datasets\".\nMost figures look blurred when zoomed in.\nResponse: Thank you for your careful review. We have corrected the typographical errors\n(Line 238) and regenerated higher-resolution images.\nReviewer #6 (Remarks on code availability):\nthe code repository is easy to follow.\nReferences\n1. Sammut SJ, Crispin-Ortuzar M, Chin SF, et al. Multi-omic machine learning\npredictor of breast cancer therapy response. Nature. Jan 2022;601(7894):623-629.\ndoi:10.1038/s41586-021-04278-5\n2. Elmarakeby HA, Hwang J, Arafeh R, et al. Biologically informed deep neural\nnetwork for prostate cancer discovery. Nature. Oct 2021;598(7880):348-352.\ndoi:10.1038/s41586-021-03922-4\n3. Ogier du Terrail J, Leopold A, Joly C, et al. Federated learning for predicting\nhistological response to neoadjuvant chemotherapy in triple-negative breast cancer.\nNat Med. Jan 2023;29(1):135-146. doi:10.1038/s41591-022-02155-w\n4. Shi Z, Huang X, Cheng Z, et al. MRI-based Quantification of Intratumoral\nHeterogeneity for Predicting Treatment Response to Neoadjuvant Chemotherapy in\nBreast Cancer. Radiology. Jul 2023;308(1):e222830. doi:10.1148/radiol.222830\n5. Tanioka M, Fan C, Parker JS, et al. Integrated Analysis of RNA and DNA from\nthe Phase III Trial CALGB 40601 Identifies Predictors of Response to Trastuzumab-\nBased Neoadjuvant Chemotherapy in HER2-Positive Breast Cancer. Clin Cancer Res.\nNov 1 2018;24(21):5292-5304. doi:10.1158/1078-0432.Ccr-17-3431\n6. Wang Y, Liu X, Yu K, et al. A generalized non-linear model predicting efficacy of\nneoadjuvant therapy in HER2+ breast cancer. iScience. Apr 21 2023;26(4):106330.\ndoi:10.1016/j.isci.2023.106330\n7. Mittendorf EA, Zhang H, Barrios CH, et al. Neoadjuvant atezolizumab in\ncombination with sequential nab-paclitaxel and anthracycline-based chemotherapy versus\nplacebo and chemotherapy in patients with early-stage triple-negative breast cancer\n(IMpassion031): a randomised, double-blind, phase 3 trial. Lancet. Oct 10\n2020;396(10257):1090-1100. doi:10.1016/s0140-6736(20)31953-x\n8. Chen L, Li H, Zhang H, et al. Camrelizumab vs Placebo in Combination With\nChemotherapy as Neoadjuvant Treatment in Patients With Early or Locally Advanced\nTriple-Negative Breast Cancer: The CamRelief Randomized Clinical Trial. Jama. Feb 25\n2025;333(8):673-681. doi:10.1001/jama.2024.23560\n9. Antonini M, Mattar A, Bauk Richter FG, et al. Real-world evidence of\nneoadjuvant chemotherapy for breast cancer treatment in a Brazilian multicenter cohort:\nCorrelation of pathological complete response with overall survival. Breast. Dec\n2023;72:103577. doi:10.1016/j.breast.2023.103577\n10. Kuenzi BM, Park J, Fong SH, et al. Predicting Drug Response and Synergy Using\na Deep Learning Model of Human Cancer Cells. Cancer Cell. Nov 9 2020;38(5):672-\n684.e6. doi:10.1016/j.ccell.2020.09.014\n11. Lipkova J, Chen RJ, Chen B, et al. Artificial intelligence for multimodal data\nintegration in oncology. Cancer Cell. Oct 10 2022;40(10):1095-1110.\ndoi:10.1016/j.ccell.2022.09.012\n12. Park S, Silva E, Singhal A, et al. A deep learning model of tumor cell architecture\nelucidates response and resistance to CDK4/6 inhibitors. Nat Cancer. Mar 5\n2024;doi:10.1038/s43018-024-00740-1\n13. Zhang H, Chen Y, Li F. Predicting Anticancer Drug Response With Deep\nLearning Constrained by Signaling Pathways. Front Bioinform. 2021;1:639349.\ndoi:10.3389/fbinf.2021.639349\n14. Wishart DS, Feunang YD, Guo AC, et al. DrugBank 5.0: a major update to the\nDrugBank database for 2018. Nucleic Acids Res. Jan 4 2018;46(D1):D1074-d1082.\ndoi:10.1093/nar/gkx1037\nReviewer #2\nReviewer #2 (Remarks to the Author):\nAlthough the researchers were still unable to provide an external validation cohort for\nindependent testing, their explanation has addressed the vast majority of my concerns.\nResponse: Thank you for your acknowledgment and supportive feedback.\nReviewer #4\nReviewer #4 (Remarks to the Author):\nThe authors have addressed all the comments.\nResponse: Thank you for your positive feedback.\nReviewer #5\nReviewer #5 (Remarks to the Author):\nThe authors have addressed the concerns/comments raised and made relevant changes\nthat have significantly strengthened the manuscript.\nResponse: Thank you for your positive assessment.\nReviewer #6\nReviewer #6 (Remarks to the Author):\nWe are pleased to see that the authors have completed a comprehensive revision that has\nsignificantly enhanced the quality of the manuscript. It is particularly important to see\nhow thoroughly the authors have investigated the questions and performed the analyses\nwe suggested in our initial review. Based on these substantial improvements, we have no\nfurther comments and believe the manuscript now meets the criteria for publication in\nNature Communications.\nResponse: Thank you for your positive and encouraging feedback.\nReviewer #7\nReviewer #7 (Remarks to the Author): Expert in breast cancer neoadjuvant chemotherapy\nand genomics; replaces Reviewer #3\nthis is an interesting manuscript which develops a model to predict the best treatment\nregimens to use for neoadjuvant treatment for breast cancer.\nThe authors have addressed many of the concerns however, there a 2 that remain\n1. pCR is most predictive of DFS for triple negative breast cancer followed by Her 2. The\nmodel appears to have the lowest AUROC for these 2 cancer subtypes. Is it possible that\nthis is due to the bias of the datasets which usually contain mostly ER+ cancers since\nthey are most common? Should this be corrected?\nResponse: Thank you for pointing out the phenomenon that the predictive performance\nof pCR for DFS varies across different subtypes, as well as the issue of the relatively\nlower performance of our model (GDnet) in the TNBC and ER-HER2+ subtypes in the\ntest set. This has deepened our understanding of the differences in the model's predictive\nperformance across various subtypes.\nWe would first like to clarify that the ability of pCR to predict DFS in different subtypes\nis not directly related to GDnet's performance in predicting pCR across these subtypes.\nSpecifically: (1) The predictive ability of GDnet for pCR is primarily associated with\ntranscriptomic indicators of tumor response and factors such as the neoadjuvant treatment\nregimen before surgery. (2) In contrast, the predictive accuracy of pCR for DFS is\nstrongly influenced by other factors, including the type of surgery performed, the\npresence of minimal residual disease (MRD) after surgery, and the efficacy of adjuvant\ntherapies—such as systemic treatment and/or radiotherapy—in eradicating any remaining\nMRD. Therefore, the prediction of pCR by GDnet and the prediction of DFS by pCR\nrepresent two independent steps. Higher performance (such as TNBC) in step 2 (pCR\npredicting DFS) does not necessarily imply high performance in step 1 (GDnet predicting\npCR).\nThen we tested patients with reported DFS or DRFS information in the included datasets,\nand as shown in the figure below, pCR indeed demonstrated strong predictive value for\nDFS/DRFS in TNBC patients (GSE32063 and GSE22226 reported DFS; GSE16446 and\nGSE25066 reported DRFS). In contrast, its predictive value was less significant in\nER+HER2- patients and appeared lowest in HER2+ patients (though the sample size\navailable for HER2+ analysis was limited, making this conclusion less robust).\nHowever, as we initially explained, the prediction of DFS/DRFS by pCR and the\nprediction of pCR by GDnet are two relatively independent steps, with no direct\nrelationship between them. Therefore, although pCR has the highest predictive value for\nDFS in TNBC, this does not necessarily imply that GDnet will also exhibit high\nperformance in predicting pCR for TNBC.\nAs for why the performance of GDnet is relatively low in the TNBC and ER-HER2+\nsubtypes in test set, we first carefully considered the distribution of molecular subtypes in\nthe training set of this study and have included this information in Supplementary Fig. 3b.\nAmong the 3,459 cases in our training set, 3,249 had complete ER and HER2\ninformation. This included 1,169 TNBC cases, 1,439 ER+HER2- cases, 274 ER+HER2+\ncases, and 367 ER-HER2+ cases (Supplementary Fig. 3b). It can be observed that\nER+HER2+ patients are the least numerous, yet the predictive performance for this\nsubtype is high. The sample size for ER-HER2+ is slightly larger than that of\nER+HER2+, but its performance is relatively low. Meanwhile, TNBC has a relatively\nlarge number of cases. Therefore, the lower pCR predictive performance observed in the\nTNBC and ER-HER2+ subtypes within the test set is unlikely to be attributable to the\ndistribution of sample sizes across subtypes.\nWe believe the reasons for the lower predictive performance in TNBC and ER-HER2+\nmay include:\n(1) Biological factors: Subtypes like TNBC are generally characterized by poorer\ntreatment response, greater challenges in predicting therapeutic efficacy, and a lack of\nreliable predictive biomarkers. Although TNBC accounts for a large proportion of the\ntraining set, its biological nature still makes treatment response challenging to predict.\n(2) Random factors in the dataset: The sample sizes for each molecular subtype in the test\nset are not very large, so random variations could contribute to the lower performance.\n(3) It is important to emphasize that we believe GDnet's performance in the ER-HER2+\nsubtype in the test set of this study is actually not poor. This is because GDnet shows\nsignificant improvement compared to the baseline transcriptome-only model (Gnet) for\nER-HER2+ subtype. Specifically, the low predictive performance of the transcriptome-\nonly Gnet in ER-HER2+ suggests that this subtype in our test set is inherently difficult to\npredict. The improved performance of GDnet demonstrates its relatively high\nperformance and its value in enhancing precision treatment.\nRegarding the question of whether sample size correction across subtypes should be\nperformed, machine learning techniques such as upsampling or downsampling (i.e.,\nresampling with replacement to increase samples in minority subtypes or reduce those in\nmajority subtypes) do exist to balance the number of patients across subtypes to a specific\nratio. However, as mentioned above, it is unlikely that the distribution of subtype sample\nsizes is the reason for GDnet’s lower performance in ER-HER2+ or TNBC. Therefore,\nthere is insufficient justification for correcting the sample size proportions across\nsubtypes. Moreover, if we were to adjust the sample size distribution, we believe it could\nactually harm GDnet’s performance. This is because molecular subtypes are not the\ninputs to GDnet—instead, transcriptome data and treatments serve as the inputs. Our\npriority is to ensure the diversity of transcriptome profiles and treatments, which the\ncurrent model already achieves. Applying additional upsampling or downsampling based\non molecular subtypes, which are not model inputs, could introduce bias. Specifically: (1)\nDownsampling would reduce the overall sample size, which is detrimental to model\ntraining. (2) Upsampling would cause certain patients’ transcriptomic profiles to appear\nmultiple times, potentially leading the model to learn individual-specific expression\npatterns and introducing additional bias.\nIn summary, while the predictive power of pCR for DFS indeed varies across subtypes,\nthis does not necessarily imply that GDnet should exhibit a high AUROC for pCR\nprediction in subtypes like TNBC. Furthermore, it is unlikely that the relatively lower\nproportion of certain subtypes in the train set is the primary cause of GDnet’s lower\nperformance in ER-HER2+ or TNBC. We attribute the observed differences to the\ndistinct biological characteristics of each subtype and potential random variations in the\ntest set. Therefore, there is insufficient justification for correcting the sample size\ndistribution across subtypes. Artificially adjusting the subtype proportions could\nintroduce additional biases. Consequently, we maintain the conventional strategy of not\napplying sample size correction.\nAlso, as mentioned previously, the manuscript would be strengthened by the inclusion of\ninformation from the institutions own patients. While it may not be feasible at this time to\nobtain the data from their own intuition, it may be possible to obtain this information\nfrom another individual institution, rather than a publicly available dataset.\nResponse: Thank you for your suggestion regarding alternative test sets. In fact, we had\nalready attempted to obtain an additional dataset from a collaborating institution during\nthe previous revision. However, due to data security concerns and potential conflicts of\ninterest related to unpublished data, we were ultimately unable to acquire a suitable\ndataset for an additional test. Furthermore, the dataset from the collaborating institution\nwas limited in sample size, subtype diversity, and treatment regimens, which would have\nbeen insufficient for comprehensively validating our model. As a next step, we are\ncurrently designing and seeking funding for a prospective interventional clinical trial.\nThis trial will enroll breast cancer patients across multiple subtypes and treatment\nregimens, utilizing GDnet to optimize clinical decision-making in neoadjuvant therapy.\nWe hope to provide higher-level evidence to validate our model's ability to improve\nclinical decisions and benefit patient outcomes in the future.\nI believe that all of the concerns of reviewer 3 were addressed.\nResponse: Thank you for your assessment and for taking the time to evaluate our revised\nmanuscript!\nReviewer #7\nReviewer #7 (Remarks to the Author):\nno new recommendations\nResponse: Thank you for your assessment and supportive feedback."}, {"round": 2, "reviewer_comments": {"reviewer_2": "Although the researchers were still unable to provide an external validation cohort for independent testing, their explanation has addressed the vast majority of my concerns.", "reviewer_4": "The authors have addressed all the comments.", "reviewer_5": "The authors have addressed the concerns/comments raised and made relevant changes that have significantly strengthened the manuscript.", "reviewer_6": "We are pleased to see that the authors have completed a comprehensive revision that has significantly enhanced the quality of the manuscript. It is particularly important to see how thoroughly the authors have investigated the questions and performed the analyses we suggested in our initial review. Based on these substantial improvements, we have no further comments and believe the manuscript now meets the criteria for publication in Nature Communications.", "reviewer_7": "this is an interesting manuscript which develops a model to predict the best treatment regimens to use for neoadjuvant treatment for breast cancer.\nThe authors have addressed many of the concerns however, there a 2 that remain\n1. pCR is most predictive of DFS for triple negative breast cancer followed by Her 2. The model appears to have the lowest AUROC for these 2 cancer subtypes. Is it possible that this is due to the bias of the datasets which usually contain mostly ER+ cancers since they are most common? Should this be corrected?\nAlso, as mentioned previously, the manuscript would be strengthened by the inclusion of information from the institutions own patients. While it may not be feasible at this time to obtain the data from their own intuition, it may be possible to obtain this information from another individual institution, rather than a publicly available dataset.\nI believe that all of the concerns of reviewer 3 were addressed."}, "author_response": "Response: Thank you for your acknowledgment and supportive feedback.\nResponse: Thank you for your positive feedback.\nResponse: Thank you for your positive assessment.\nResponse: Thank you for your positive and encouraging feedback.\nResponse: Thank you for pointing out the phenomenon that the predictive performance of pCR for DFS varies across different subtypes, as well as the issue of the relatively lower performance of our model (GDnet) in the TNBC and ER-HER2+ subtypes in the test set. This has deepened our understanding of the differences in the model's predictive performance across various subtypes.\nWe would first like to clarify that the ability of pCR to predict DFS in different subtypes is not directly related to GDnet's performance in predicting pCR across these subtypes. Specifically: (1) The predictive ability of GDnet for pCR is primarily associated with transcriptomic indicators of tumor response and factors such as the neoadjuvant treatment regimen before surgery. (2) In contrast, the predictive accuracy of pCR for DFS is strongly influenced by other factors, including the type of surgery performed, the presence of minimal residual disease (MRD) after surgery, and the efficacy of adjuvant therapies—such as systemic treatment and/or radiotherapy—in eradicating any remaining MRD. Therefore, the prediction of pCR by GDnet and the prediction of DFS by pCR represent two independent steps. Higher performance (such as TNBC) in step 2 (pCR predicting DFS) does not necessarily imply high performance in step 1 (GDnet predicting pCR).\nThen we tested patients with reported DFS or DRFS information in the included datasets, and as shown in the figure below, pCR indeed demonstrated strong predictive value for DFS/DRFS in TNBC patients (GSE32063 and GSE22226 reported DFS; GSE16446 and GSE25066 reported DRFS). In contrast, its predictive value was less significant in ER+HER2- patients and appeared lowest in HER2+ patients (though the sample size available for HER2+ analysis was limited, making this conclusion less robust).\nHowever, as we initially explained, the prediction of DFS/DRFS by pCR and the prediction of pCR by GDnet are two relatively independent steps, with no direct relationship between them. Therefore, although pCR has the highest predictive value for DFS in TNBC, this does not necessarily imply that GDnet will also exhibit high performance in predicting pCR for TNBC.\nAs for why the performance of GDnet is relatively low in the TNBC and ER-HER2+ subtypes in test set, we first carefully considered the distribution of molecular subtypes in the training set of this study and have included this information in Supplementary Fig. 3b. Among the 3,459 cases in our training set, 3,249 had complete ER and HER2 information. This included 1,169 TNBC cases, 1,439 ER+HER2- cases, 274 ER+HER2+ cases, and 367 ER-HER2+ cases (Supplementary Fig. 3b). It can be observed that ER+HER2+ patients are the least numerous, yet the predictive performance for this subtype is high. The sample size for ER-HER2+ is slightly larger than that of ER+HER2+, but its performance is relatively low. Meanwhile, TNBC has a relatively large number of cases. Therefore, the lower pCR predictive performance observed in the TNBC and ER-HER2+ subtypes within the test set is unlikely to be attributable to the distribution of sample sizes across subtypes.\nWe believe the reasons for the lower predictive performance in TNBC and ER-HER2+ may include:\n(1) Biological factors: Subtypes like TNBC are generally characterized by poorer treatment response, greater challenges in predicting therapeutic efficacy, and a lack of reliable predictive biomarkers. Although TNBC accounts for a large proportion of the training set, its biological nature still makes treatment response challenging to predict.\n(2) Random factors in the dataset: The sample sizes for each molecular subtype in the test set are not very large, so random variations could contribute to the lower performance.\n(3) It is important to emphasize that we believe GDnet's performance in the ER-HER2+ subtype in the test set of this study is actually not poor. This is because GDnet shows significant improvement compared to the baseline transcriptome-only model (Gnet) for ER-HER2+ subtype. Specifically, the low predictive performance of the transcriptome-only Gnet in ER-HER2+ suggests that this subtype in our test set is inherently difficult to predict. The improved performance of GDnet demonstrates its relatively high performance and its value in enhancing precision treatment.\nRegarding the question of whether sample size correction across subtypes should be performed, machine learning techniques such as upsampling or downsampling (i.e., resampling with replacement to increase samples in minority subtypes or reduce those in majority subtypes) do exist to balance the number of patients across subtypes to a specific ratio. However, as mentioned above, it is unlikely that the distribution of subtype sample sizes is the reason for GDnet’s lower performance in ER-HER2+ or TNBC. Therefore, there is insufficient justification for correcting the sample size proportions across subtypes. Moreover, if we were to adjust the sample size distribution, we believe it could actually harm GDnet’s performance. This is because molecular subtypes are not the inputs to GDnet—instead, transcriptome data and treatments serve as the inputs. Our priority is to ensure the diversity of transcriptome profiles and treatments, which the current model already achieves. Applying additional upsampling or downsampling based on molecular subtypes, which are not model inputs, could introduce bias. Specifically: (1) Downsampling would reduce the overall sample size, which is detrimental to model training. (2) Upsampling would cause certain patients’ transcriptomic profiles to appear multiple times, potentially leading the model to learn individual-specific expression patterns and introducing additional bias.\nIn summary, while the predictive power of pCR for DFS indeed varies across subtypes, this does not necessarily imply that GDnet should exhibit a high AUROC for pCR prediction in subtypes like TNBC. Furthermore, it is unlikely that the relatively lower proportion of certain subtypes in the train set is the primary cause of GDnet’s lower performance in ER-HER2+ or TNBC. We attribute the observed differences to the distinct biological characteristics of each subtype and potential random variations in the test set. Therefore, there is insufficient justification for correcting the sample size distribution across subtypes. Artificially adjusting the subtype proportions could introduce additional biases. Consequently, we maintain the conventional strategy of not applying sample size correction.\nAlso, as mentioned previously, the manuscript would be strengthened by the inclusion of information from the institutions own patients. While it may not be feasible at this time to obtain the data from their own intuition, it may be possible to obtain this information from another individual institution, rather than a publicly available dataset.\nResponse: Thank you for your suggestion regarding alternative test sets. In fact, we had already attempted to obtain an additional dataset from a collaborating institution during the previous revision. However, due to data security concerns and potential conflicts of interest related to unpublished data, we were ultimately unable to acquire a suitable dataset for an additional test. Furthermore, the dataset from the collaborating institution was limited in sample size, subtype diversity, and treatment regimens, which would have been insufficient for comprehensively validating our model. As a next step, we are currently designing and seeking funding for a prospective interventional clinical trial. This trial will enroll breast cancer patients across multiple subtypes and treatment regimens, utilizing GDnet to optimize clinical decision-making in neoadjuvant therapy. We hope to provide higher-level evidence to validate our model's ability to improve clinical decisions and benefit patient outcomes in the future.\nI believe that all of the concerns of reviewer 3 were addressed.\nResponse: Thank you for your assessment and for taking the time to evaluate our revised manuscript!"}, {"round": 3, "reviewer_comments": {"reviewer_7": "no new recommendations"}, "author_response": "Response: Thank you for your assessment and supportive feedback."}]

Decision letter