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 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. Cylin 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.

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 e.g. 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.

The ECC network robustly facilitates resistance-associated protein dynamics.

A) Clustered heatmap of a 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 6 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.

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 variation. 1000 unique initial condition 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 initial condition values are changed. B) The frequency of which initial condition variation (left) induces resistance, in each key output protein. Varying initial condition values usually does not induce resistance, however, for select parameter sets, most initial conditions produce resistance associated dynamics. This analysis is repeated for parameter variation (right). Contrary to varying initial condition values, varying parameters usually induces resistance to a small degree, but there are no sets of initial conditions that are universally resistant. C) The number of unique dynamic categories produced by initial condition variation (left), versus parametric variation (right). Varying initial conditions 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.

MDN analysis 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.

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 9) 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.

Validation of MDN-based predictions.

Qualitative comparisons between single-cell data and our MDN analyses. A) Single-cell time-course responses to CDK4/6 inhibition from Yang et al. (55)). 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 localization 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 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 simulations quantified by the distribution of dynamic classes across both datasets.