Systematic analysis of network-driven adaptive resistance to CDK4/6 and oestrogen receptor inhibition using meta-dynamic network modelling
- Department of Biochemistry and Molecular Biology, Faculty of Medicine, Nursing and Health Sciences, Monash University, Australia
- Biomedicine Discovery Institute, Monash University, Australia
- Computational Systems Oncology Program, South Australian immunoGENomics Cancer Institute (SAiGENCI), The University of Adelaide, Australia
- Australian Research Council Centre of Excellence for the Mathematical Analysis of Cellular Systems (MACSYS), Australia
Figures
Figure 1
A workflow of the Meta-Dynamic Network (MDN) modelling pipeline.
As with all modelling, the first step is asking a useful question. In our case, this is usually related to understanding complex protein signalling dynamics. The second step is defining the network; this involves selecting nodes and identifying relationships between them. The third step is translating the network schematic to a mathematical representation, which is typically a combination of ordinary differential equations and biochemical rate laws. The fourth step is the generation of thousands of model instances, each with a unique set of properties, for example parameter and initial condition values. The fifth step is the simulation of each model instance and the bulk analysis of the simulations. In this work, we analyse the distribution of the qualitative dynamics of the network’s protein species.
Figure 2
Network schematic of the early cell cycle (ECC) network.
Network schematic showing the detailed biochemical reactions that regulate the initiation of the G1 phase and transition through to the S phase of the cell cycle. The network includes three mitogenic inputs, two mitogenic signalling pathways, PI3K and MAPK, key transcription factors that promote G1 phase transition, key CDK-cyclin complexes and their regulators, a CDK4/6 inhibitor, and an estrogen receptor (ER) inhibitor. Canonically, mitogenic stimulation leads to the cascading activation of the mitogenic signalling pathways, which converge on the synthesis of cyclin D. The synthesis of cyclin D then promotes the activation of CDK4/6, which consequently drives the synthesis of cyclin E. Cyclin E activates CDK2, leading to the synthesis of cyclin A and the formation of CDK2-cyclin A complexes signals the initiation of S phase. ER and CDK4/6 inhibitors are believed to disrupt this linear chain of events and thus prevent proliferation.
Figure 3 with 3 supplements
The ECC network robustly facilitates resistance-associated protein dynamics.
(A) Clustered heatmap of the time-course profiles of a representative subgroup model instances, normalised to maximum value. Time-course profiles are of mono-phosphorylated Rb (pRb), the downstream target of active CDK4/6, and in response to CDK4/6 and ER inhibition. Blue represents low concentration and red represents high concentration. (B) The frequency of the six dynamic categories for the model’s active protein forms, across 100,000 model instances. Blue represents low frequency and red represents high frequency. Species have been clustered to highlight species that have highly correlated distributions of dynamics. (C) The frequency of model instances displaying simultaneous sensitivity (blue), or resistance in at least one of the key output proteins (red). Only 2.79% of model instances show simultaneous sensitivity across all three key output proteins. The percentage of resistant model instances was calculated by measuring the number of unique model instances that show resistance-associated dynamics for at least one of the key output proteins.
Figure 3—figure supplement 1
Methods for generating groups of 100,000 unique model instances.
(A) Model instances are produced by combining a base set of initial conditions (ICs) with a unique set of parameters. Parameter values in each unique parameter set are generated by randomly selecting values for parameters from between the range 10–5 and 104. (B) Model instances are produced by combining a base parameter set with a unique set of ICs. IC values in each unique IC set are produced by randomly selecting values for ICs between the range 100 and 104. (C) Measuring the ability of parameter variation versus IC variation to induce resistance. First, we selected 1000 model instances from the group of model instances produced in (A) that demonstrated decreasing dynamics for one of the key output proteins, and 1000 model instances produced in (B) that show the same. We then generated 1000 new, unique model instances per model instance in group A by taking the parameter set of each group A model instance as the base parameter set and combining it with all of the IC sets from group B, see left side of panel. This produced a total of 1,000,000 unique model instances. We then simulated each of these new model instances and determined the protein dynamic of each of the key output proteins. We then calculated the distribution of observed resistance-associated dynamics for the rows and columns. The distribution of the rows is representative of IC variation and the distribution of the columns is representative of parameter variation.
Figure 3—figure supplement 2
Full repertoire of signalling dynamics identified by MDN, including filtering for strong CDK4/6-Cyclin D suppression.
(A) Heatmap showing the frequency of each dynamic for a selection of active protein species, across 100,000 model instances (parametric variation), that is the network’s meta-dynamic map. Note that CDK4/6cycD shows a sustained decreasing behaviour in only 59.93% of model instances. (B) Heatmap of dynamic frequencies, filtered for model instances that display strong suppression of CDK46cycD. Filtering in this way has very little effect on the distribution of protein dynamics.
Figure 3—figure supplement 3
Bar graph showing the total change in the frequency of dynamics when the number of model instances is doubled.
Specifically, we calculated the Mean Squared Error between the distributions of protein dynamics for model instance populations of size N and 2 N. Y-axis is the calculated MSE and is in log10 scale, and the X axis shows the size of population N.
Figure 4 with 2 supplements
Reaction coefficients are a stronger driver of adaptive-resistance dynamics than protein abundance.
(A) Overview of the computational pipeline undertaken to compare the effects of parametric variation with initial condition (IC) variation. 1000 unique IC sets are applied to a base parameter set that demonstrates a decreasing dynamic for the key output protein being analysed. Each of the new model instances are simulated and the category of the output protein’s dynamic is determined. Then the number of times the model instances produce a resistance-associated dynamic is measured. This entire process is repeated for 1000 base parameter sets. The final step is measuring the distribution of resistance behaviours when the parameter values are changed and comparing the distribution when the IC values are changed. (B) The frequency of which IC variation (left) induces resistance, in each key output protein. Varying IC values usually does not induce resistance, however, for select parameter sets, most ICs produce resistance associated dynamics. This analysis is repeated for parameter variation (right). Contrary to varying IC values, varying parameters usually induces resistance to a small degree, but there are no sets of ICs that are universally resistant. (C) The number of unique dynamic categories produced by IC variation (left), versus parametric variation (right). Varying ICs produces a moderate variety of qualitative dynamics, and can produce all 6 dynamic categories, but varying parameters almost always induces every single possible dynamic category for the key output proteins.
Figure 4—figure supplement 1
MDN modelling exploring initial condition variation.
(A) Heatmap showing the frequency of each dynamic for a selection of active protein species, across 100,000 model instances when initial condition values are varied. The base parameter set used was initially randomly selected. Blue represents low frequency and red represents high frequency. (B) Clustering of the heatmap in A to highlight proteins with similar distributions of dynamic categories.
Figure 4—figure supplement 2
Varying conserved totals shifts the dynamic equilibria of the relevant protein species.
50 different model instances (colours); each with an identical base parameter set but different, randomly selected conserved totals for their constituent protein species. Each model instance was simulated as laid out in the methods section of the main text. The first column demonstrates the shifted equilibria with no stimulation, the second under (constant) mitogenic stimulation and the third after drug treatment. Each row is derived from the labelled protein species. Variety can be seen in both steady state equilibria and in the qualitative dynamics before steady state is reached.
Figure 5 with 5 supplements
MDN modelling identifies core sub-networks that facilitate resistance.
(A) Overview of the process undertaken to investigate the relationship between parameter knockdown and resistance-dynamic features to identify resistance driving parameter signatures. Each parameter is knocked down by 20%, one at a time, and the effect on the protein dynamics is measured. If the knockdown decreases resistance it scores a 1, otherwise it scores a 0. This process is repeated for each parameter, one at a time, for each model instance. The end result is a matrix where the columns represent parameters, and the rows represent the model instances. This matrix can then be analysed using clustering methods to identify parametric resistance signatures. (B) Heatmap produced by hierarchical clustering of the parameters that contribute to rebounding pppRb, for model instances that display rebounding pppRb. (C) Comparing parametric resistance signatures with the parameters that most frequently contribute to resistance overall. Parameters are ranked by how frequently they contribute to rebounding pppRb. Clusters are aligned with overall ranking. The lower bar graph represents the correlation between the overall ranking and each cluster specific ranking. This heatmap shows that some parametric resistance signatures align closely with the overall ranking, but some clusters are quite different, emphasising the context specificity of resistance.
Figure 5—figure supplement 1
Mean and SD of the absolute log10 transformed parameter values for each protein-dynamic combination.
Y-axis represents the range of possible log10 parameter values, the X-axis represents each individual parameter (in the order they appear in the IQM reaction file).
Figure 5—figure supplement 2
Patterns in the absolute values of parameters driving resistance associated dynamics for key output proteins.
Heatmaps showing the results of the hierarchical clustering of the raw parameter values for each protein-dynamic combination. Very little clustering can be observed, suggesting that there are no patterns in the absolute values of parameters with respect to resistance associated dynamics.
Figure 5—figure supplement 3
Breakdown of the various measurable features of a protein’s dynamic response following drug perturbation.
Resistance features (red) are a subset of these dynamic features and are defined by the recovery or increase in the activity of a protein following drug perturbation, as represented by the increase in concentration of the active form of the protein. The transition from black to red lines represent an increase in the resistance being displayed.
Figure 5—figure supplement 4
Top and bottom 10 ranked parameters for each protein-dynamic combination.
Ranking is based on the overall contribution of each parameter to the resistance features for dynamic of each protein-dynamic combination. For biphasic and increasing dynamics, the ability of each parameter to contribute to the maximum concentration and final concentration was measured, and for the rebounding dynamics, the contribution to rebound and final concentration was measured. The y-axis refers to the percentage of the model instances in which the parameters along the x-axis contribute to the relevant resistance features, for that particular protein dynamic-output protein category.
Figure 5—figure supplement 5
Heatmaps comparing the overall ranking of each protein-dynamic combination to the cluster-specific rankings.
Rank is determined by how frequently each parameter contributes to the resistance features of the dynamic of each protein-dynamic combination.
Figure 6 with 12 supplements
Identification of subnetworks driving pppRb rebound-mediated resistance within the ECC network.
(A) Cluster-specific parameter signatures overlaid with the ECC network to highlight subnetworks that drive resistance through rebound of pppRb. Only the five largest clusters (out of nine) are displayed. (B) Parameter knockdowns for the top-scoring parameters in three select clusters that display divergent resistance driving parameter signatures. Black represents treatment with CDK4/6 and ER inhibitors alone, decreasingly-red lines represent the addition of an increasingly potent parameter knockdown. These results highlight the context specificity of parameter knockdown when attempting to prevent adaptive resistance.
Figure 6—figure supplement 1
Parameter-signature ECC network overlays for each of the nine protein-dynamic combinations, presented side-by-side for easier comparison of protein-dynamic combinations.
Note that Figure 5—figure supplement 5A-I displays the results for the individual protein-dynamic combination in more detail.
Figure 6—figure supplement 2
Cluster-specific parameter signatures with the ECC network, shown here for the protein-dynamic pRb-INC combination.
Each figure (Figure 5—figure supplement 5A-I) is the overlay of the respective cluster-specific parameter signatures with the ECC network, for each protein-dynamic combination. Each figure displays the five largest clusters for each protein-dynamic combination.
Figure 6—figure supplement 3
Cluster-specific parameter signatures with the ECC network, shown here for the protein-dynamic pRb-BIP combination.
Each figure (Figure 5—figure supplement 5A-I) is the overlay of the respective cluster-specific parameter signatures with the ECC network, for each protein-dynamic combination. Each figure displays the five largest clusters for each protein-dynamic combination.
Figure 6—figure supplement 4
Cluster-specific parameter signatures of the ECC network, shown here for the protein-dynamic pRb-REB combination.
Each figure (Figure 5—figure supplement 5A-I) is the overlay of the respective cluster-specific parameter signatures with the ECC network, for each protein-dynamic combination. Each figure displays the five largest clusters for each protein-dynamic combination.
Figure 6—figure supplement 5
Cluster-specific parameter signatures with the ECC network, shown here for the protein-dynamic CDK2cycE-INC combination.
Each figure (Figure 5—figure supplement 5A-I) is the overlay of the respective cluster-specific parameter signatures with the ECC network, for each protein-dynamic combination. Each figure displays the five largest clusters for each protein-dynamic combination.
Figure 6—figure supplement 6
Cluster-specific parameter signatures with the ECC network, shown here for the protein-dynamic CDK2cycE-BIP combination.
Each figure (Figure 5—figure supplement 5A-I) is the overlay of the respective cluster-specific parameter signatures with the ECC network, for each protein-dynamic combination. Each figure displays the five largest clusters for each protein-dynamic combination.
Figure 6—figure supplement 7
Cluster-specific parameter signatures with the ECC network, shown here for the protein-dynamic CDK2cycE-REB combination.
Each figure (Figure 5—figure supplement 5A-I) is the overlay of the respective cluster-specific parameter signatures with the ECC network, for each protein-dynamic combination. Each figure displays the five largest clusters for each protein-dynamic combination.
Figure 6—figure supplement 8
Cluster-specific parameter signatures with the ECC network, shown here for the protein-dynamic pppRb-INC combination.
Each figure (Figure 5—figure supplement 5A-I) is the overlay of the respective cluster-specific parameter signatures with the ECC network, for each protein-dynamic combination. Each figure displays the five largest clusters for each protein-dynamic combination.
Figure 6—figure supplement 9
Cluster-specific parameter signatures with the ECC network, shown here for the protein-dynamic pppRb-BIP combination.
Each figure (Figure 5—figure supplement 5A-I) is the overlay of the respective cluster-specific parameter signatures with the ECC network, for each protein-dynamic combination. Each figure displays the five largest clusters for each protein-dynamic combination.
Figure 6—figure supplement 10
Cluster-specific parameter signatures with the ECC network, shown here for the protein-dynamic pppRb-REB combination.
Each figure (Figure 5—figure supplement 5A-I) is the overlay of the respective cluster-specific parameter signatures with the ECC network, for each protein-dynamic combination. Each figure displays the five largest clusters for each protein-dynamic combination.
Figure 6—figure supplement 11
Summary of all the parameter-signature clusters (subnetworks) that drive resistance identified for nine protein-dynamic combinations.
This was derived from the individual analysis results in Figure 5—figure supplements 4 and 5A-I.
Figure 6—figure supplement 12
Heatmap showing the clustered differential gene-dependencies of hundreds of cancer cell lines.
Genes were filtered for moderate dependency and then subject to hierarchical clustering to identify common signatures of gene dependency across approximately 1000 cancer cell lines. The gene signatures identified in this analysis are quasi-similar to the idea of drug resistance signatures identified using our MDN analysis.
Figure 7
Validation of MDN-based predictions.
Qualitative comparisons between single-cell data and our MDN modelling analyses. (A) Single-cell time-course responses to CDK4/6 inhibition from Yang et al., 2020. The reporter readout reflects phosphorylation-dependent nucleo-cytoplasmic shuttling; the y-axis is the nuclear/cytoplasmic fluorescence ratio (not fold-change to t=0). Baseline values at t=0 therefore vary across cells (typically ~1–1.8) due to inherent differences in reporter localisation prior to treatment. Red traces: resistance-associated dynamics; blue: sensitivity-associated; grey: no/limited response. (B) Frequency of each dynamic category in panel A. (C) Representative MDN modelling simulations for CDK4/6·CycD and CDK2·CycE drawn from the model ensemble (parametric variation). Colours denote the same dynamic classes as in A. The distribution of resistant and sensitive dynamics from our simulations are highly correlated with the single-cell data. (D) Agreement between single-cell data and MDN modelling simulations quantified by the distribution of dynamic classes across both datasets.
Tables
Key resources table
| Reagent type (species) or resource | Designation | Source or reference | Identifiers | Additional information |
|---|---|---|---|---|
| Software, algorithm | MATLAB | https://au.mathworks.com/ | N/A | |
| Software, algorithm | IQM | https://iqmtools.intiquan.com/ | N/A |
Additional files
-
Supplementary file 1
Early cell cycle (ECC) model in biochemical reaction format.
Model file describing the ECC network in biochemical reaction format, including model species, initial conditions, parameters, variables, and reaction-rate expressions.
- https://cdn.elifesciences.org/articles/87710/elife-87710-supp1-v1.txt
-
Supplementary file 2
Early cell cycle (ECC) model ordinary differential equations (ODEs).
Plain-text representation of the ECC model containing the system of ODEs, initial conditions, parameter values, model variables, and reaction-rate equations.
- https://cdn.elifesciences.org/articles/87710/elife-87710-supp2-v1.txt
-
Supplementary file 3
Early cell cycle (ECC) model in Systems Biology Markup Language (SBML) format.
SBML representation of the ECC model, including model species, parameters, assignment rules, and biochemical reactions.
- https://cdn.elifesciences.org/articles/87710/elife-87710-supp3-v1.pdf
-
Supplementary file 4
Typeset ordinary differential equations (ODEs) for the early cell cycle (ECC) model.
Typeset representation of the complete system of ODEs and reaction-rate equations defining the ECC model, provided to facilitate interpretation and reproducibility.
- https://cdn.elifesciences.org/articles/87710/elife-87710-supp4-v1.pdf
-
Supplementary file 5
Nominal parameter units and values.
Nominal parameter values and units used for the ECC model. These values provide the base parameter set used in analyses examining variation in protein abundance/initial conditions.
- https://cdn.elifesciences.org/articles/87710/elife-87710-supp5-v1.xlsx
-
Supplementary file 6
Robust parameter signatures capable of facilitating resistance-associated dynamics.
Tables summarising the robust parameter-signature clusters identified by MDN modelling analysis across the resistance-associated protein–dynamic combinations. For each cluster, the tables report its frequency and the parameters contributing to the corresponding resistance signature.
- https://cdn.elifesciences.org/articles/87710/elife-87710-supp6-v1.xlsx
-
Supplementary file 7
Known resistance mechanisms and corresponding resistance-driving model parameters.
Comparison of experimentally reported mechanisms of resistance to CDK4/6 and oestrogen receptor inhibition with the corresponding resistance-driving parameters identified by MDN modelling analysis.
- https://cdn.elifesciences.org/articles/87710/elife-87710-supp7-v1.xlsx
-
MDAR checklist
- https://cdn.elifesciences.org/articles/87710/elife-87710-mdarchecklist1-v1.docx
Download links
A two-part list of links to download the article, or parts of the article, in various formats.
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Systematic analysis of network-driven adaptive resistance to CDK4/6 and oestrogen receptor inhibition using meta-dynamic network modelling
eLife 12:RP87710.
https://doi.org/10.7554/eLife.87710.3
