Systematic Analysis of Network-driven Adaptive Resistance to CDK4/6 and Estrogen Receptor Inhibition using Meta-Dynamic Network Modelling

Joint Public Review:

In this manuscript, the authors proposed an approach to systematically characterise how heterogeneity in a protein signalling network affects its emergent dynamics, with particular emphasis on drug-response signalling dynamics in cancer treatments. They named this approach Meta Dynamic Network (MDN) modelling, as it aims to consider the potential dynamic responses globally, varying both initial conditions (i.e., expression levels) and biophysical parameters (i.e., protein interaction parameters). By characterising the "meta" response of the network, the authors propose that the method can provide insights not only into the possible dynamic behaviours of the system of interest but also into the likelihood and frequency of observing these dynamic behaviours in the natural system.

The authors studied the Early Cell Cycle (ECC) network as a proof of concept, specifically focusing on PI3K, EGFR, and CDK4/6, with particular interest in identifying the mechanisms that cancer could potentially exploit to display drug resistance. The biochemical reaction model consists of 50 equations (state variables) with 94 kinetic parameters, described using SBML and computed in Matlab. Based on the simulations, the authors concluded the following main points: a large number of network states can facilitate resistance, the individual biophysical parameters alone are insufficient to predict resistance, and adaptive resistance is an emergent property of the network. Finally, the authors attempt to validate the model's prediction that differential core sub-networks can drive drug resistance by comparing their observations with the knock-out information available in the literature. The authors identified subnetworks potentially responsible for drug resistance through the inhibition of individual pathways. Importantly, some concerns regarding the methodology are discussed below, putting in doubt the validity of the main claims of this work.

While the authors proposed a potentially useful computational approach to better understand the effect of heterogeneity in a system's dynamic response to a drug treatment (i.e., a perturbation), there are important weaknesses in the manuscript in its current form:

(1) It is unclear how the random parameter sets (i.e., model instances) and initial conditions are generated, and how this choice biases or limits the general conclusions for the case studied. Particularly, it is not evident how the kinetic rates are related to any biological data, nor if the parameter distributions used in this study have any biological relevance.
(2) Related to this problem, it is not clear whether the considered 100,000 random parameter samples sufficiently explore parameter space due to the combinatorial explosion that arises from having 94 free parameters, nor 100,000 random initial conditions for a system with 50 species (variables).
(3) Moreover, the authors filter out all the cases with stiff behaviour. This filtering step appears to select model parameters based on computational convenience, rather than biological plausibility.
(4) Also, it is not clear how exactly the drug effect is incorporated into the model (e.g., molecular inhibition?), nor how it is evaluated in the dynamic simulations (e.g., at the beginning of the simulation?). Moreover, in a complex network, the results may differ depending on whether the inhibition is applied from the start or after the network has reached a stable state.
(5) On the same line, the conclusions need to be discussed in the context of stability, particularly when evaluating the role of initial conditions. As stable steady states are determined by the model parameters, once again, the details of how the perturbation effect is evaluated on the simulation dynamics are critical to interpret the results.
(6) The presented validation of the model results (Fig. 7) is only qualitative, and the interpretation is not carefully discussed in the manuscript, particularly considering the comparison between fold-change responses without specifying the baseline states.

Author Response:

We would like to express our gratitude to the reviewers for the time and effort dedicated to evaluating our manuscript. We are committed to addressing each of the comments and recommendations they have presented.

It appears that a majority of the feedback emphasizes the need for clarity and expanded explanations. We acknowledge these points and are confident that offering a clearer exposition and delving into further details will notably enhance the manuscript. In our initial draft, our intention was to ensure accessibility to non-mathematical readers by minimizing technical jargon. However, the feedback underscores the importance of certain details, particularly for those well- versed in ODE modelling, and the need to provide complete information.

While we find the reviewers' feedback invaluable, it is worth noting that none of the critiques suggest a fundamental change in our presented analyses. Below, we offer brief responses to the primary critiques mentioned in the public review:

1. The first notable comment pertains to the selection criteria for parameter and initial condition values. This critique is indeed valid. In brief, parameter values were chosen from a range of 10^- 5 to 10^4, representing rates from 10 femtomolar/minute to 10 micromolar/minute, spanning a biologically plausible spectrum. It is conceivable that values outside this range exist but are exceedingly rare. Similarly, initial conditions were chosen within the range 10^0 to 10^4, typically represented in nM.

2. The second comment highlights the challenges in systematically determining a full spectrum of parameter sets with 94 free parameters. In our observations, as we expanded the number of model instances, the distribution of protein dynamics exhibited minimal variation. A doubling of model instances from 100,000 to 200,000 led to less than a 1% error change. This error was calculated based on the differences across every protein species and dynamic category. These findings suggest that examining more than 100,000 model instances neither shifts the dynamic distributions significantly nor unveils new resistance mechanisms. We are committed to presenting these analyses more comprehensively in the revised manuscript.

3. The query about the appropriateness of filtering our models based on computational feasibility is pertinent. Our contention is that this filter does not exclude a significant number of model instances. Furthermore, stiff ODEs generally arise from extremely high reaction rates, which are exceedingly rare in a biological context. Thus, their exclusion only filters out exceedingly rare biological contexts, and only a small proportion of the time.

4 & 5) Clarifications sought about the simulations will be addressed. Though we feel the details were implicitly incorporated, we will make them explicit in the subsequent version.

1. The final major comment underscores the qualitative nature of our validation, which we agree. Currently, we are exploring experimental techniques or datasets for a more robust validation. In our next revision, we will ensure a more in-depth discussion of the validation in the manuscript's discussion section.

Once again, thank you for your valuable feedback. We look forward to submitting a revised version that addresses all concerns and enhances the manuscript's overall quality.