Network dynamics-based cancer panel stratification for systemic prediction of anticancer drug response

Cancer is a complex disease involving multiple genomic alterations that disrupt the dynamic response of signaling networks. The heterogeneous nature of cancer, which results in highly variable drug response, is a major obstacle to developing effective cancer therapy. Previous studies of cancer therapeutic response mostly focus on static analysis of genome-wide alterations, thus they are unable to unravel the dynamic, network-specific origin of variation. Here we present a network dynamics-based approach to integrate cancer genomics with dynamics of biological network for drug response prediction and design of drug combination. We select the p53 network as an example and analyze its cancer-specific state transition dynamics under distinct anticancer drug treatments by attractor landscape analysis. Our results not only enable stratification of cancer into distinct drug response groups, but also reveal network-specific drug targets that maximize p53 network-mediated cell death, providing a basis to design combinatorial therapeutic strategies for distinct cancer genomic subtypes.

In the present manuscript, Choi et al. use the p53 network as a paradigm to study how genomic alterations in cancer cells rewire the topology of a signaling network and thereby change its dynamics upon stimulation as well as the response to pharmacological perturbations. It represents an extension of their previous work investigating the p53 network using attractor landscape analysis (Choi et al. 2012 Science Signaling). In this previous study, the authors established a boolean model that allows them to predict different states of the p53 network and associated phenotypes (proliferation, cell cycle arrest or cell death) in normal cells in presence or absence of DNA damage. This model is then adapted to represent known alterations in the cancer cell line MCF7 and to predict the effect of combination treatments such as DNA damage (etoposide) and feedback inhibition (Mdm2 inhibitor Nutlin and Wip1 siRNA).
The authors now use publicly available genome-wide datasets of somatic mutations and copy number variations to determine genomic alterations within the p53 network for 83 cancer cell lines (including MCF7) that have wild-type p53 and functional caspases. The identified alterations were discretised to represent either constantly activated or inactivated nodes and summarized in 45 differentially wired network. Next, the authors systematically investigated how the rewired networks respond to perturbation of five drugable network nodes / links in the presence or absence of DNA damage, focusing on network states that represent the phenotypes mentioned above. Using a drug synergy score, they discriminate between additive, synergistic and antagonistic combination treatments. Finally, they validate model predictions in three cell lines (U-2 OS, MCF7, A549) by determining cell death upon combinations of BCL2, Mdm2-p53 and Wip1 inhibition with DNA damage induction.
The authors present an interesting approach to systematically address how genomic alterations in cancer cells affect the function of signaling networks and contribute to the heterogeneity of the disease. The mapping and classification of different mutation types is a valuable resource for the field. However, the paper is not easily accessible. The authors introduce many new terms such as "non-intuitive target" and "critical determinant" that are not commonly used, not intuitive and not well explained. Moreover, the authors seem to rely on the reader being familiar with their previous publication. I would have appreciated a more detailed description of modelling and analysis methods to be able to understand how data is generated and conclusions are reached.
At the same time, the study seems like a direct extension from previous work, were the authors already generated a cancer cell specific network (MCF7) and predicted and validated drug combinations (Etoposide, Nutlin, Wip1 -kd). Experimental validation of new predictions is rather limited. They test two additional cell lines compared to the previous publication and add a Bcl2 inhibitor. The validated cell lines share three alterations, which means that a large space of alterations, e.g. cell lines with wild type p14/Arf, remains unexplored. A more systematic validation would be needed to complement the systematic model analysis. Moreover, the phenotypic read out relies on counting dead cells at 48h by microscopy. No errors bars are provided and no indication of cell numbers or how often experiments were repeated can be found. Here the authors should use additional quantitative methods with high statistical power for exampled based on flow cytometry or time-lapse imaging. It would also be helpful if the deviation of model predictions and experimental data would be quantified and put into context.
There are also some restrictions inherent to the approach of the authors. As they focus on an abstracted p53 network in p53 wt cells, alterations that are not part of the initial network and could represent truly "non-intuitive" targets are neglected. Also, p53 mutations that often acquire new functions, are also not included in the analysis. The use of a boolean network together with discretised genomic alterations limit the predictive power. For example, the authors can't distinguish between weakly and strongly activating mutations or consider varying efficiencies of inhibitors. While these limitations are inherent to the approach, the authors should at least discuss them openly.
Minor points: * Figure legends are extremely brief and lack detail. There are many small issues with the figures: * Figure 1: no explanation of e.g. node color after mapping of alterations is given, panel b is not described at all. * Figure 2: abbreviations such as HOMDEL, AMP etc should be spelled out at least once; heat maps in panel a should contain labels and there are typos in y-axis label for the bar graph next to pie charts ("Altered Smaples"). * Figure 3: there are no numbers on scale for heat map * The naming switches between genes and nodes, e.g. CDKN2A -p14ARF, PPM1D -Wip1 in Figure  2. Here a consistent labelling would be easier for the reader.
They experimentally validate some key findings in cell lines as well. This part could be expanded.
The manuscript is poorly written -it needs a lot more work to make it easy to read and assess; In particular the authors tend to repeat themselves quite a bit.
Thus the paper should be completely rewritten.
But globally I don't have any major concern at this point in terms of the scientific impact."

[RESPONSE]
We thank the reviewer for the very positive comments regarding the novelty and significance of our work. In light of the reviewer's comment, we have restructured and completely rewritten large sections of the manuscript, in particular the Results and Discussion section, to clarify the methods of our network dynamics-based analysis and points/terms that were difficult to understand in the original manuscript. For instance, we added an overview of the modeling framework and computational procedures in the beginning of the Results section to illustrate how data were generated and results were analyzed. We also removed repeated sentences and ambiguous terms throughout the manuscript. Moreover, we extended the experimental validation from 3 to 8 cancer cell lines that now encompass all three network subgroups in terms of network topology. In general, the simulation and experimental results showed similar drug response characteristics, indicating the effectiveness of our modeling approach to predict drug response profiles from genomic data. We hope the very significant revision of the manuscript and addition of new experimental data now makes the manuscript up-to-the-standard for the reviewer to assess.

Response to the specific comments of Reviewer 2 [COMMENT #1].
"The authors present an interesting approach to systematically address how genomic alterations in cancer cells affect the function of signaling networks and contribute to the heterogeneity of the disease. The mapping and classification of different mutation types is a valuable resource for the field. However, the paper is not easily accessible. The authors introduce many new terms such as "non-intuitive target" and "critical determinant" that are not commonly used, not intuitive and not well explained. Moreover, the authors seem to rely on the reader being familiar with their previous publication. I would have appreciated a more detailed description of modelling and analysis methods to be able to understand how data is generated and conclusions are reached."

[RESPONSE]
We again appreciate the reviewer's very positive comment on the significance and general interest of our work. Following the reviewer's helpful suggestions, we have rewritten the manuscript comprehensively, for instance, either removing the unclear terms, such as "non-intuitive target", or adding more detailed description to clarify the quantitative and biological meaning of the terms, such as "critical determinant" (refer to page 12 in the revised manuscript). Briefly, if an inhibitory perturbation results in different change of major cellular response phenotype in cancer cell vs. normal cell (i.e., p53 network with no genomic alteration), we considered the particular genomic alterations present in the cancer cell as the determinant of drug response for the network subtype. We then were able to identify a minimal subset of such genomic alterations for a given inhibitory treatment, which we termed "critical determinant". Intuitively, if a drug-induced perturbation results in the same major cellular response phenotype in cancer and normal cell, there is no critical determinant in the cancer cell for this drug treatment.
"At the same time, the study seems like a direct extension from previous work, were the authors already generated a cancer cell specific network (MCF7) and predicted and validated "There are also some restrictions inherent to the approach of the authors. As they focus on an abstracted p53 network in p53 wt cells, alterations that are not part of the initial network and could represent truly "non-intuitive" targets are neglected. Also, p53 mutations that often acquire new functions, are also not included in the analysis. The use of a boolean network together with discretised genomic alterations limit the predictive power. For example, the authors can't distinguish between weakly and strongly activating mutations or consider varying efficiencies of inhibitors. While these limitations are inherent to the approach, the authors should at least discuss them openly." [RESPONSE] We agree the simplified p53 network model that we used as an example of the network dynamics-based approach is limited in the aspects that the reviewers mentioned.
However, our computational approach is still widely applicable to model complex regulatory networks, where continuous modeling is not appropriate or possible, due to insufficient quantitative information on kinetic parameters 1, 2 . And we do think our approach provides a foundational framework that can be further developed to address the limitations pointed out by the reviewer. For instance, we can easily expand the simplified p53 network to a larger, more comprehensive functional network that incorporates additional components crucial to oncogenesis, metastasis and/or tumor response. We could also expand the Boolean network model to multi-valued logical model or use fuzzy logic, instead of Boolean logic, to describe varying degrees of activation/inhibition and drug efficacy 2, 3 . As for cancer cells with mutant p53, the drug effect is likely mediated by alternative signaling pathway so the network dynamics-based analysis should probably be developed and performed other signaling pathway, beyond the p53 network. We added the above discussion regarding the limitation of our approach in page 20 of the revised manuscript. In Figure 1a, the node color represents the status of the node activity. A black (white) node means that the node is constantly activated (inactivated) and a gray node means that the status is dependent on the activity of a given input. Figure 1b illustrates the general framework of our method to predict drug response using attractor landscape analysis of network dynamics.
In the revised figure legend, we describe the p53 network model, selection of drug target and definition of response phenotype based on attractor states. "* Figure 2: abbreviations such as HOMDEL, AMP etc should be spelled out at least once; heat maps in panel a should contain labels and there are typos in y-axis label for the bar graph next to pie charts ("Altered Smaples")."

[RESPONSE]
We corrected the gene names in Figure 2 and also added a full description of "HOMDEL" and "AMP" in the revised manuscript and "* The naming switches between genes and nodes, e.g. CDKN2A -p14ARF, PPM1D -Wip1 in Figure 2. Here a consistent labelling would be easier for the reader."

[RESPONSE]
In the revised manuscript, we consistently used protein names to describe the nodes. For example, CDKN2A and PPM1D have been changed to p14ARF and Wip1 throughout the manuscript. As mentioned in my initial comment, conclusions from the experimental validation would be stronger if the authors find a way to quantify the deviation of model predictions and experimental data and, if possible, provide some insight into how much better a cell line specific prediction is compared to "random" predictions, e.g. generated by shuffling predictions for the corresponding perturbation for each cell line. This would also emphasize "surprising" cell line specific prediction. I also need to mention that the authors still limit experimental validation to only one of the predicted outcomes (cell death), and ignore the other two (cell cycle arrest and proliferation)." [RESPONSE]: We again thank the reviewer for the very helpful comments. Regarding the reviewer's major concern of a lack of quantitative analysis of the deviation of modeling Regarding the point raised by the reviewer that our simulation results showed very similar perturbation responses for A2780 and CAL51, we think this is probably due to the fact that the only difference between A2780 and CAL51 under the simplified p53 network that we employed is MDMX, which is constantly activated in A2780, but not in CAL51. And both networks share two key alterations, including constantly inactivated PTEN and constantly activated AKT (i.e., PTEN(I) and AKT(A)) (refer to Fig. R1 (Fig. R1).
In addition, MDMX is always "ON" in all attractor states of CAL51 under the distinct perturbations, even though MDMX is not constantly activated in the control CAL51 network without inhibitory perturbation. This indicates that AKT(A) is a dominant input that keeps MDMX "ON" in the CAL51-specific network model that we curated, despite that MDMX has multiple upstream regulators (Fig. R2).  predictions, we observed relatively weak correlations between the experimentally observed and randomly predicted responses, compared to that between the experimental data and predicted responses by our approach (p < 0.001, Wilcoxon rank sum test). In addition, each RMSE of our cell line specific predictions is significantly smaller than that of the random prediction (p < 0.001, Wilcoxon rank sum test). Overall, these results demonstrate that our network dynamics analysis performs substantially better than random prediction. We added this discussion in page 19 of the revised manuscript.
In this study, we chose to focus on performing detailed analysis of the cell death response because clinical response of tumor regression or delay in tumor growth often correlates with the cell death response and variation in cell death response to cancer chemotherapy is likely the most crucial factor that distinguishes sensitive and resistant tumors.