Introduction

Drug resistance is a widespread phenomenon across all cancer types and is a major obstacle to the development of curative therapeutic strategies [1–3]. Cellular heterogeneity is a known driver of drug resistance, wherein the treatment of a heterogenous population of cancer cells with cytostatic or cytotoxic drugs creates a selective pressure that results in the survival and expansion of any cells that are capable of overcoming the effects of said drugs [4–6]. While this general phenomenon is well established [7–9], it is usually less clear exactly how and why some cells are able to overcome a drug treatment and others are not. Heterogeneity is a catch-all term used to describe any and all differences between cells, however the majority of the differences between cells ultimately converge on differences in how their constituent proteins behave, i.e. their protein dynamics. To deepen our understanding of how tumour heterogeneity drives drug resistance, we must explore the relationship between cellular heterogeneity, protein dynamics and drug resistance.

The efficacy of a targeted therapy largely comes down to how well and how long the targeted protein is suppressed, and how frequently and reliably this suppression results in either cytostasis or apoptosis. A cell can therefore be considered resistant if the target protein is insufficiently suppressed or is initially suppressed but later recovers, or if the suppression of the protein is insufficient to stimulate cytostasis or apoptosis [10, 11]. Proteins are embedded in complex networks and protein dynamics are dictated by the properties of the networks in which they reside. While resistance is often due to direct effects, such as reduced drug-target binding affinity [12, 13] or excessive drug efflux from tumour cells [14, 15], it can also occur due to changes in the state of a cells biochemical networks that counteract the effects of the drug; a phenomenon referred to as adaptive resistance [16–18]. PI3K, EGFR and CDK4/6 are just few prime examples of proteins that have been targeted in the treatment of cancer that display acute and robust adaptive resistance [19–21]. These studies also support the notion that it is rarely the behaviour of a single ‘gatekeeper’ protein that drives adaptive resistance, but an entire protein network that acts and is acted upon by the target protein [22, 23]. Frequently however, we possess a limited knowledge and appreciation of the network-level mechanistic relationships that underpin adaptative resistance.

Understanding the relationship between cellular heterogeneity and adaptive drug resistance requires an ability to explore and characterise the heterogeneity that exists both within a tumour and between patients. This can be achieved in silico by investigating the behaviour of a network, i.e. the dynamics of its constituent proteins, over a broad range of network conditions, allowing us to observe the full spectrum of possible dynamics that can be displayed by a given network topology. A particularly useful technique for investigating protein signalling networks is Ordinary Differential Equation (ODE) modelling [24]. An ODE model is a mathematical representation of how the protein species within a network interact and evolve over time. ODE models have been extensively used to simulate and predict network-level responses to perturbations, such as growth-factor stimulation or drug treatment [25–28]. In an ODE model, the interactions between the network’s constituent proteins are converted into mathematical formulations using well-established biochemical rate laws [29–33]. The end result is a set of ODEs – the model – that can be numerically solved using specialised ODE solvers, allowing one to simulate, and thus predict, how the concentrations of the protein species will change over time or in response to perturbations [29–33].

An ODE model is comprised of state variables, initial conditions, model parameters, and defined inputs and outputs [29–31]. Typically, the parameters represent the strength of the protein interactions and state variables correspond to the concentrations of the protein species. Initial conditions (ICs) are the values of the state variables at the starting point of the solution and represent the basal expression levels of the protein species within the model. Model input(s) are usually a perturbation to the modelled network, such as a growth factor stimulation or a drug inhibition. We can model network heterogeneity by producing a suite of ODE models, each with their own set of parameter values and/or initial conditions. In this manner, each model ‘instance’ represents a unique cellular context and produces a unique set of dynamic behaviours for its constituent proteins. By varying the parameters and initial conditions over wide ranges, we can delineate an extreme upper limit to the heterogeneity in dynamics facilitated by a network topology.

The process of creating models almost always involves a degree of abstraction, either explicitly, due to conscious decisions by the creator, or implicitly, due to imperfect knowledge of the system being modelled. This abstraction forces a model’s parameters to capture information beyond that which is being explicitly modelled and can potentially render the concept of a ‘true’ or ‘biologically accurate’ parameter value less meaningful. There is an array of factors that can influence protein interactions: post-translational modifications, catalysts/scaffolds, pH, etc., and they can all individually influence protein interactions over several orders of magnitude. In the context of cancer, this variability in protein interactions can become extreme, with mutations altering the properties of proteins and dramatically changing how strongly they interact with other proteins. By varying the parameters and initial conditions over wide ranges, we also potentially capture protein dynamics driven by factors we haven’t explicitly modelled, or we aren’t aware exist.

Drug resistance is a problem almost as old as drug treatments; and many computational modelling frameworks have been developed to investigate the emergence of drug resistance during cancer treatment [34–36]. A comprehensive review of the modelling frameworks developed to explore the relationship between cellular heterogeneity and drug resistance, in particular, was written by Chisholm et al. in 2016 [35]. Most of the models covered within focused on population dynamics and were used to identify key processes that govern how resistance emerges within a heterogenous cell population. Notably, none of the models in this review attempted to link heterogenous intracellular processes with adaptive resistance mechanisms, or resistance in general. Perhaps most similar to the ensemble approach developed within this project, He et al. recently developed an ODE-based modelling technique termed RACIPE (random circuit perturbation) [37]. In contrast to our work, RACIPE was a quasi-Monte Carlo method developed to reduce the necessity of identifying precise kinetic parameter values on which model predictions are supposedly dependent. While their technique does involve the exploration of a range of parameter values, it does so to make predictions about broadly robust outcomes, not to link parameter variation to heterogenous network dynamics. Thus, there is a need for a novel computational framework capable of systematically characterising network-dynamic heterogeneity and its relationship with adaptive resistance.

To address this need, we have developed a new modelling framework, coined Meta-Dynamic Network (MDN) modelling, which enables us to study the effects of network heterogeneity on protein signalling dynamics. We then applied MDN modelling to an Early Cell Cycle (ECC) network model to systematically explore and illuminate adaptive resistance mechanisms that arise in response to targeted CDK4/6 and Estrogen Receptor (ER) inhibition. The ECC network model developed herein incorporates the two major mitogenic signalling pathways that are largely responsible for cell cycle initiation, and the G1-S phase cell cycle regulatory network.

Using MDN and the ECC network model, we have demonstrated that variation in both the strength of protein interactions (parameters) and the basal expression of proteins (state variable values) can affect the qualitative shape and features of a network’s protein dynamics. However, we then provide evidence that parameter variation induces a much greater variety in the qualitative shape and features of the network’s protein dynamics and is therefore much more likely to give rise to adaptive resistance dynamics. We further show that there is a surprising large number of contexts in which the network topology of the ECC network can facilitate resistance dynamics, and that we can dissect these contexts to systematically characterise their underpinning kinetic drivers. Finally, we analysed existing publicly available biological data to validate the existence of a spectrum of signalling dynamics in response to CDK4/6 inhibitors and propose a novel relationship between population signalling dynamics and population level drug sensitivity.

Materials and methods

ODE model construction, modelling, and simulations

To simulate the change in protein concentrations over time in response to drug perturbations the ECC network (Figure 2) was converted into a series of ODEs. In this process, protein relationships are converted into mathematical equations that describe the rate of change of a protein species using well-established biochemical rate laws. In this model, catalysed chemical reactions, such as kinase-mediated phosphorylations, were modelled using simple second order kinetics, complex formations were modelled using the law of Mass Action, and protein transcription, translation and degradation were modelled as single processes, using first-order kinetics. The final model consisted of 50 state variables, 94 parameters, 2 mitogenic inputs and 2 drug inputs. A more detailed description of model scope is given in the Supplementary Information; and the Biochemical Reaction, ODE and SBML format model files are given in supplementary Model Files S1 - S3.

All model simulations were implemented in MATLAB (The MathWorks. Inc. 2019a). The IQM toolbox (http://www.intiquan.com/intiquan-tools/) was used for two primary purposes: to generate and convert the reaction schematic to a series of ODEs and to compile the IQM file into a MEX file to accelerate simulation speed. IQM enables the generation of ODEs via the ability to describe the reaction schematic more intuitively by allowing the entry of biochemical reactions that describe how each species evolves over time, which are generally the rate laws mentioned above, and automatically generating the ODEs from the resulting file. This file is an IQM model structure that can be used for simulation and analysis and can also be converted into a MEX file so that the ODEs can be solved using the SUNDIALS ODE solver suite, dramatically speeding up solving time compared to the in-house MATLAB solvers. Each simulation of our model consisted of three phases. The first was a starvation phase to ensure a steady state equilibrium was reached before mitogenic stimulation. In the second phase the model was stimulated with an in silico dose of 100 nM of each mitogenic input (insulin and FGF). Again, this phase was allowed to proceed until all protein dynamics reached a steady state equilibrium. In the third and final phase the model was stimulated with 10 nM of each drug input (CDK4/6i and Ei) and the response of each network protein species was stored.

Filtering of randomly sampled parameters

Part of the MDN modelling process involves the random generation of parameter sets that are used to create unique model instances. To improve the likelihood that these parameter sets were representative of possible biological contexts, a number of filters were applied to exclude inappropriate parameter sets. However, we tried to keep the filters as minimal as possible to avoid excluding rare but possible network states. The primary filter was the ability of our ODE solver to actually solve model simulations, as the vast majority of randomly generated parameter sets produced systems that were too stiff, i.e. models that were too difficult for the software to solve numerically. Modelling multi-time-scale biological processes often produces ODE models that are relatively stiff, however, this filter acts to exclude biologically implausible parameter sets that result in overly stiff ODE systems [38, 39]. The additional secondary filters were included to ensure that the solutions to ODEs reached an equilibrium state in simulation of phases 1 and 2, proteins responded to mitogenic stimulations and drug perturbations, and that protein species never possessed negative concentration values. Critically, the filtering process required a complex of cyclin D and CDK4 (CDK46cycD), and active ER (aERa) - the direct targets of CDK4/6 and ER inhibitors, respectively, to decrease by at least 5% of their respective initial values. This filter ensured that the direct drug targets were at least minimally responsive to drug perturbation, without imposing too strict a requirement of strong drug-mediated suppression.

Quantitative definitions of dynamic patterns

A major component of MDN modelling is the analysis of the qualitative shape of individual protein dynamics. This enables us to group together model instances that produce broadly similar behaviours for any proteins being investigated. We can then analyse these groups to identify and elucidate shared network structures that drive the behaviour. To this end, we defined six qualitative protein dynamic categories: increasing, decreasing, biphasic, rebound, no-response and other.

To assess the qualitative shape of each protein dynamic, we developed a series of mathematical restraints that were capable of sorting protein simulations into their appropriate categories. First, each protein dynamic was normalised to its initial value. A protein’s dynamic was categorised as increasing (INC) if (i) the final concentration exceeded 20% of its initial concentration, (ii) over the entire simulated time period it never dropped below 10% of its initial concentration, and (iii) the final concentration was within 10% of the maximum concentration. It was categorised as decreasing (DEC) if (i) the final concentration was at least 20% lower than the initial concentration, (ii) never went above 10% of its initial concentration, and (iii) the final concentration was within 10% of the minimum concentration. It was categorised as biphasic (BIP) if (i) the maximum concentration exceeded the initial concentration by 20% and (ii) the final concentration was at least 10% lower than the maximum concentration. It was categorised as rebound (REB) if (i) the minimum concentration was at least 20% lower than the initial concentration and (ii) the final concentration was at least 10% bigger than the minimum concentration. It was categorised as no-response (NRP) if the concentration never exceeded more or less than 10% of the initial concentration. Finally, it was categorised as other (ETC) if the protein’s dynamic did not fit one of the previously described categories.

Key Model Output Proteins

There is strong evidence to suggest that the linear cascade of cell cycle regulator progression, wherein the activation of CDK4/6 leads to the activation of CDK2-cyclin E complex, which in turn leads to the activation of CDK2-cyclin A, is overly simplistic [40, 41]. If this traditional understanding was correct, the inhibition of CDK4/6 activity would robustly inhibit cell cycle progression, yet this is not the case [42]. A number of studies have suggested that the activation or reactivation of proteins downstream of CDK4/6 can enable a cell to overcome the cell cycle inhibitory effects of CDK4/6 targeting drugs and enable cell cycle progression. Proteins that have been shown to demonstrate this phenomenon include phosphorylated Rb and the active CDK2-cyclin E complex [21, 43, 44]. As such, we have focused on the dynamics of monophosphorylated Rb (pRb), hyperphosphorylated Rb (pppRb) and the CDK2-cyclin E complex (CDK2cycE) throughout this project, as key mediators of adaptive resistance to CDK4/6 and ER inhibitors.

Results

An Overview of Meta-Dynamic Network (MDN) Modelling

To explore the link between cellular heterogeneity, protein dynamics and drug resistance, we first developed a novel modelling framework coined Meta-Dynamic Network (MDN) modelling. The principal purpose of the MDN modelling framework is to explore how parametric and state variable variation influence the qualitative shape of protein dynamics within a given network. By varying model parameters and initial state variable values (i.e. initial conditions, or ICs) we are able capture many sources of cellular heterogeneity and characterise their influence on protein signalling dynamics. Varying parameter values captures heterogeneity in the interactions between proteins that arises due to phenomena such as cell-specific mutations, post-translational modifications, and epigenetic modifications. Varying initial state variable values captures heterogeneity in protein concentration that arises due to phenomena such as stochastic protein expression, copy number variation and tissue-specific expression programmes.

Once a large number of unique model instances have been generated, each with their own unique parameter and/or state variable values, we can simulate each model’s response to drug treatment and analyse the resulting shape of each protein’s dynamic. From this analysis we can determine if there are contexts where key output proteins, such as downstream targets of the drug perturbation, demonstrate protein dynamics associated with adaptive resistance. Finally, we can group model instances that display similar dynamics together and perform further analyses to systematically identify the mechanistic causes capable of driving adaptive resistance dynamics. An overview of the MDN modelling process can be seen in Figure 1.

Figure 1.

Protein dynamics that are associated with adaptive resistance are those where the protein is insufficiently suppressed over a sustained time period, i.e. increasing, biphasic and rebound [18–20, 22, 45]. Increasing is the most dramatic adaptive-resistance response, where the drug perturbation not only fails to suppress key output proteins, but actually increases their activity/expression. Biphasic is similar to increasing in its initial phase, where the outputs actually increase, however the degree of the second decreasing phase likely influences the degree of resistance conferred by this particular behaviour. Rebound is similar to biphasic wherein the degree of resistance conferred by this dynamic is likely influenced by how strongly the protein recovers. While rebound is probably the most commonly observed adaptive-resistance driving dynamic [18–23], for the purpose of this project we consider all three of these behaviours as adaptive-resistance associated behaviours and consider decreasing as the only sensitivity associated behaviour.

Construction of a mechanistic model of the Early Cell Cycle (ECC) network

The ECC network model generated in this project was built to capture the signalling events driving entry into G1 and progression through the early cell cycle events, up to the G1-S phase transition. A detailed network schematic of the ECC’s biochemical reactions can be seen in Figure 2. We nominally include two Tyrosine Kinase Receptors (TKRs), the Insulin Receptor (INSR) and the Fibroblast Growth Factor Receptor (FGFR), as they strongly stimulate the PI3K and MAPK pathways respectively but are both capable of stimulating each pathway.

Figure 2:

The two major mitogenic signalling pathways, PI3K and MAPK, were included for three reasons. The first is that they are known to potently stimulate the production of cyclin D, the primary activator of CDK4/6 [46, 47]. The second is that activating mutations in these pathways are some of the most common mutations across all cancer types and the final reason is that these same mutations have been shown to facilitate resistance to CDK4/6 inhibitors [48, 49]. The inclusion of proteins within this network is not exhaustive, due to computational limitations, but our model captures the core signalling relationships and network structure. The final component included in this model was the network regulating G1 entry and progression into S phase. This includes the relevant CDKs, cyclins and their inhibitors, as well as the transcription machinery regulating cyclin expression. Clinical approval for the use of CDK4/6 inhibitors has traditionally been limited to patients with ER+ (estrogen receptor-positive) breast cancer who were concurrently receiving ER (estrogen receptor) antagonists but were still experiencing disease progression [50–52]. Due to the ER’s ability to promote CDK4/6 inhibition efficacy and the clinical practice of simultaneously prescribing ER antagonists with CDK4/6 inhibitors [50–52], the ER and its downstream signalling events were also included in the model.

The ECC network model was constructed using ODEs that represent biochemical interactions as a series of ordinary differential equations based on established kinetic laws. More detailed description of the model is given in the Methods and Materials. The full ODE equations and reaction rates are given in supplementary Model Files S1 - S3.

The ECC network robustly facilitates resistance-associated protein dynamics

To investigate the meta-dynamics of the ECC network, i.e. the range of dynamics that can be facilitated by the topology of the ECC network, we first generated model instances with varying parameter values. Specifically, we generated 100,000 model instances, each with the same set of ICs but a unique set of randomly generated parameter values. Each model instance was simulated as laid out in the ODE model construction, modelling, and simulations (see Methods and Materials), and the resulting time-course drug response of each of the network’s proteins was stored. This was repeated for all model instances and the entire process was repeated in triplicate.

We found that parameter variation resulted in a wide range of drug perturbation responses for many of the ECC network proteins. For example, Figure 3A displays a clustered subset of monophosphorylated Rb (pRb) time-course drug responses, showing a fairly extreme variety of drug-response dynamics that are possible in response to parameter variation. pRb not only possesses a quantitative spectrum of intra-category responses, which is probably expected, but it can also display a range of qualitative dynamics, and in fact, can demonstrate all six categories of qualitative drug response dynamics (see Methods and Materials).

Figure 3:

Categorising the dynamics of each protein across all of the model instances allowed us to calculate the distribution/frequency of each dynamic for each protein. The distribution of behaviours for the ‘active’ forms of the model’s protein species can be seen in the heatmap in Figure 3B. This overall distribution of protein behaviours is akin to a map that shows how the network state can affect the response of individual proteins to a drug treatment, in this case, the response to simultaneous CDK4/6 and ER inhibition. For many proteins, particularly proteins far upstream of the drug targets, the network state has only a little influence on drug response. This is particularly well demonstrated by the proteins in cluster 1, which shows close to 100% no-response (NRP) dynamics.

Other proteins are much more strongly influenced by network state, where, under the right conditions, a protein’s response can flip from sensitivity to resistance. The proteins in clusters 4 and 5 show this phenomenon to a moderate degree but also show a predilection towards decreasing dynamics. The predilection to a decreasing dynamic may be due to their proximity to the drug targets CDK4/6 and ER and/or the linearity of the interactions between them and the direct drug targets. It is particularly interesting to observe that most of the major promoters of cell cycle progression downstream of CDK4/6 cluster together (cluster 3, Figure 3B) and are all strongly influenced by the network state. This suggests that heterogeneity can strongly influence whether CDK4/6 inhibition will translate to an effective network-level response and prevent cell cycle progression.

To further assess the capability of the ECC network to facilitate adaptive resistance, we dug into the dynamic distributions of the key output proteins: pRb, pppRb (hyper-phosphorylated Rb) and CDK2cycE. Recall that we have defined resistance dynamics as increasing, biphasic and rebound, and sensitive behaviours as decreasing. While the frequency of sensitive dynamics was much higher than even the sum of resistance dynamics for all three proteins, the ratio of model instances displaying sensitivity to those displaying resistance was much lower than we expected. The ratios of sensitive to resistant behaviours, as percentages of 100,000 model instances, were 30: 14, 22: 8, 19: 11 (rounded) for pRb, pppRb and CDK2cycE, respectively (see supplementary Figure S1A). These results suggest that for every 2-3 model instances/network states that facilitate sensitivity, there is 1 that facilitates resistance.

We then analysed the model instances further to identify how frequently sensitive and resistant behaviours occurred simultaneously across the three key output proteins. Our results suggest that broad sensitivity across all three output proteins is quite rare, and that resistance is remarkably common (Figure 3C). We observed that only 2.79% of model instances show simultaneous sensitivity across all three output proteins (Figure 3C). In contrast, when evaluating model instances where at least 1 of the output proteins displays a resistance associated dynamic, we observe that this occurs in just over a quarter, 26.04%, of all model instances (Figure 3C).

Given that CDK46cycD is only strongly suppressed in just under 60% of the model instances (Figure 3C), we hypothesised that the resistance displayed by the key output proteins may be largely due to insufficient suppression of CDK46cycD. To explore this, we filtered out the model instances that showed strong suppression of CDK46cycD and re-generated the distribution of protein dynamics. The resulting heatmap can be seen in supplementary Figure S1B. Surprisingly, the distribution of behaviours did not appear to change significantly, with the key output proteins showing very similar levels of adaptive resistance dynamics. This result suggests that the resistance dynamics of the key output proteins are bona fide adaptive resistance mechanisms and not just a result of insufficient suppression of the drug target.

Taken together, these results show that in the face of heterogeneity, the ECC network is robustly capable of facilitating adaptive resistance and that key output proteins display adaptive resistance even when the drug targets are being robustly suppressed.

Protein interactions are a stronger driver of adaptive-resistance dynamics than protein expression

Having investigated the influence of parametric variation on drug response, we then wanted to explore the effect of IC variation (basal expression/activity) in a similar fashion. To this end, we generated 100,000 model instances (in triplicate), each with the same set of ‘nominal’ parameter values (see supplementary Table S1), and a unique set of ICs, and repeated the analysis as performed previously. The resulting distribution of dynamics can be seen in supplementary Figure S2A, and the clustering of this distribution can be seen in Figure S2B.

Interestingly, we found that compared to parametric variation, IC variation induced much less variety in protein response dynamics. Most proteins demonstrated a singular dominant behaviour and very few proteins exhibited both sensitive and resistance dynamics. This seems to suggest that IC variation is less capable of influencing a proteins drug response. However, we thought that this may have been due to the selected nominal parameter set and may not be a universal phenomenon. To explore this hypothesis, we compared the ability of parametric variation and IC variation to induce adaptive resistance. We did this by simulating a large number of model instances while hold the parameter values constant and measuring how frequently randomly generated ICs induced a switch from a decreasing dynamic to a resistant dynamic, focusing on our key output proteins. We then repeated this process but held the ICs constant and randomly generated parameter values.

Specifically, we first selected a parameter set that produced a decreasing dynamic of pppRb and then generated 1000 random IC sets. We then combined the 1000 IC sets with each parameter set to create 1000 model instances, each with a constant parameter set and unique set of ICs. We then simulated each model instance to measure how frequently the behaviour of the key output protein switched from a sensitive dynamic to a resistant dynamic. This entire process was then repeated for 1000 parameter sets, allowing us to observe how frequently IC variation induced adaptive resistance independent of the chosen parameter set. We then inverted this process, generating random parameter sets for sensitivity-displaying IC sets, to observe how frequently parameter variation induced adaptive resistance independent of the chosen IC set. An overview of this process can be seen in Figure 4A.

Figure 4:

We found that for the vast majority of parameter sets, varying the initial conditions never induced resistance in any of the three output proteins (Figure 4B, left). There was however, a very small number of parameter sets where altering the initial conditions frequently induced resistance. In contrast, altering the parameter values always induced resistance in a small number of the IC sets; but there were no IC sets wherein parameter variation strongly induced resistance, independent of the protein chosen (Figure 4B, right).

To further characterise the difference between parameter and IC variation, we then investigated the variety of dynamics induced by either parameter or IC variation. For each master parameter or IC set, we measured how many unique dynamic categories were produced across their respective 1000 model instances. We found that varying the ICs still induced a moderate amount of variety in protein dynamics (Figure 4C, left). However, parameter variation almost always produced all six dynamic categories (Figure 4C, right).

These results provide clear evidence that changes in both parameter values and protein expression can facilitate the emergence of adaptive resistance. However, changes in parameter values appear to be much more capable of inducing large, qualitative shifts in protein dynamics compared to changes in protein expression. The ability of expression variation to induce resistance also seems to be dependent on the master parameter set and there also appears to be no IC set that can induce resistance independent of the parameter values. These results imply that heterogeneity in protein interactions is much more likely, and independently able to facilitate the emergence of adaptive resistance compared to heterogeneity in protein expression. Due to the independent ability of parameter variation to strongly induce adaptive resistance, we chose to focus on this particular form of heterogeneity in the following studies.

Protein interactions co-ordinate to drive adaptive resistance, independent of their individual strengths

The data produced by our initial MDN pipeline enabled us to identify thousands of model instances where the key output proteins displayed adaptive resistance dynamics. Using these model instances, we next wanted to explore whether there were shared network features driving resistance, or if each model instance was unique in its ability to generate resistance. To explore this in a systematic manner, we extracted subsets of model instances, one for each key output protein and resistance behaviour combination (i.e. protein-dynamic combination) and performed a series of analyses on their parameters.

First, we investigated the mean and standard deviation of each parameter in each protein-dynamic combination to see if there were any individual parameter values underpinning resistance. We found that the mean of most parameters was between 0.1 and 1, but the standard deviation usually spanned 6 orders of magnitude (supplementary Figure S3). We did observe that all three receptor-activation associated parameters (kc1f1-INSR, kc10f1-FGFR, kc15f1-ER) were generally much lower in value than the rest, but further investigation ruled that this was generally due to higher values causing the model’s system of ODEs to become too stiff, and therefore unsolvable, and so were filtered out. We then undertook hierarchical clustering of each protein-dynamic combination’s parameter values to try and identify if there were groups of parameter values driving resistance. However, we were unable to find any meaningful clusters (supplementary Figure S4). Together, these results led us to conclude that there were no particular individual parameters or groups of parameter values that were responsible for driving resistance.

Having found no patterns in the raw, or absolute, parameter values, we decided to investigate how frequently each parameter contributed to the increase or recovery of protein activity following drug perturbation, for each output protein. To this end, we investigated the effects of parameter knockdowns on the following dynamic features: maximum concentration, final concentration and/or rebound, hereafter collectively referred to as ‘resistance features’. These features were chosen as they represent either an increase or recovery of protein activity, post drug perturbation, that could potentially enable a cell to overcome the drug insult. The relevant resistance features for biphasic and increasing dynamics were maximum concentration and final concentration, and the relevant features for rebound dynamics were rebound and final concentration. Rebound is not applicable to biphasic and increasing dynamics as it simply hasn’t occurred, and the maximum concentration is only applicable to rebounding dynamics if the rebound overshoots the initial concentration, which is then covered by final concentration.

To measure the influence of parameters on the resistance features of the key output proteins we performed a high-throughput perturbation analysis wherein we reduced the value of each parameter by 20%, at the time of the drug treatment, and measured the effect on the output protein’s resistance features. A breakdown of the measured dynamic features, including the resistance features just described, is given in supplementary Figure S5.

If the knockdown of a parameter resulted in a reduction in the resistance features of the output proteins’ dynamic, it scored a 1; otherwise it scored a 0. Repeating this process for each parameter and model instance resulted in an M x N matrix, where M is the number of model instances (protein and dynamic dependent), and N is the number of perturbed parameters (94), for each protein-dynamic combination. See Figure 5A for a visual representation of this, and the succeeding overall process.

Figure 5:

Initially, we wanted to determine how frequently each parameter contributed to the resistance features of each output protein. To achieve this, we summed up the number of times a parameter scored a 1 in the preceding analysis across all model instances, for each protein-dynamic combination. The results were plotted as bar graphs (top and bottom 10 parameters), which can be seen in supplementary Figure S6. The graphs showed clear rankings in how frequently each parameter contributed to the adaptive resistance features of each protein dynamic. Each protein-dynamic combination displayed a unique ranking but there also appeared to be a significant overlap in the top-scoring parameters. There also appeared to be more intra-protein overlap than inter-protein overlap, and pppRb and CDK2cycE shared much more similarity with each other than either with pRb. It was also interesting to observe that for biphasic pppRb and CDK2cycE, their top-scoring parameters accounted for approximately 90% of all model instances, indicating that these parameters must be critical in driving this particular dynamic for these proteins. The remaining top scoring parameters for each protein-dynamic subset only accounted for between 50-80% of model instances, suggesting each protein-dynamic combination possesses a degree of heterogeneity in the parameters that drive their adaptive resistance dynamics.

Together, the above results show that it is not the value of individual parameters that drive adaptive resistance, but rather it is the manner in which parameters co-ordinate to shape the quantitative and qualitative features of a protein’s dynamic that drives adaptive resistance.

MDN can identify a full spectrum of core sub-networks that facilitate resistance

The observation that the top-scoring parameters generally do not broadly drive adaptive resistance dynamics suggests that there may be a number of different network states that are capable of driving a given resistance dynamic. To investigate this possibility, we subjected the matrices produced in the previous analysis to hierarchical clustering to try and find groups of model instances with shared resistance-driving parameter signatures. As a prime example, we first focused on rebounding hyperphosphorylation of Rb; both as the hyperphosphorylation of Rb strongly promotes cell cycle progression, and because previous studies have demonstrated the role of this protein-dynamic combination in CDK4/6 inhibitor resistance [21, 53, 54]. Figure 5B displays a heatmap showing how rebounding pppRb model instances cluster, with respect to the contribution of their parameters to the rebounding dynamic.

The clustering of the pppRb-rebound model instances produced 9 robust clusters, where robust clusters are defined as those that contain at least 1% of the total protein-dynamic combination model instances, and the top-scoring parameter in the cluster accounts for at least 80% of its constituent model instances. Simply put, the clusters accounted for a high enough proportion of the model instances, and the top scoring parameter was robustly representative of the model instances within the cluster. The left-most column of the heatmap in Figure 5C shows the overall ranking for how frequently each parameter contributes to the rebound dynamic of pppRb, and the remaining 9 columns represent the 9 clusters, aligned with the overall ranking. We observed a fairly strong consensus amongst the top-scoring parameters, with the majority of the cluster differences coming from the middle and lower parameters (Figure 5C). Of note, this pattern appeared to be broadly true when we analysed all the 9 protein-dynamic groups (supplementary Figure S7).

Despite broad similarities between the clusters, many of the individual clusters did not correlate particularly well with the overall ranking and possessed clear differences in their resistance-driving parameter signatures (Figure 5C). These results suggest that the same drug response dynamic for the same protein in two different cells can be driven by different network states. To highlight the similarities and differences between the subnetworks underpinning resistance driven by pppRb rebound, we overlayed the parameter signatures of the five largest clusters on the ECC network schematic (Figure 6A). This revealed significant overlap between the clusters around the formation of the CDK2-cyclin E complex and the hyperphosphorylation of Rb itself. These results were to be expected as the CDK2-cyclin E complex is directly responsible for the hyperphosphorylation of Rb and confirms the logical hypothesis that CDK2 activation can drive the reactivation of Rb phosphorylation following CDK4/6 and ER inhibition. However, that almost all of the clusters shared this core network suggests that this subnetwork alone is likely sufficient to enable the reactivation of Rb phosphorylation on its own and would be strongly selected for by cancer during treatment. Similarly, but less expectedly, we also found strong overlap at the level of Myc (Figure 6A), suggesting a strong likelihood for its involvement in the development of resistance.

Figure 6:

On the other hand, clusters that have less overlap highlight potential differences in underlying resistance mechanisms between patients experiencing seemingly similar drug responses. For the pppRb-rebound group, these clusters included those driven by the PI3K and MAPK pathways, protein degradation regulation (GSK3B), and p21/27 inhibition of the CDK2-cyclin E complex (Figure 6A). Consistently, when we targeted (inhibited) parameters related to these signalling nodes across different clusters we did not observe universal suppression of pppRb; instead, we saw cluster-specific suppression of pppRb (Figure 6B).

Expanding this analysis to include all 9 protein-dynamic combinations, we observed that resistance-promoting subnetworks were widespread throughout the ECC network (supplementary Figure S8 and Figures S9A-I). Further, we identified just under 100 robust clusters/subnetworks capable of facilitating resistance dynamics (supplementary Table S2A – J). Finally, to highlight the overlap between protein-dynamic clusters, we combined all of the individual clusters of all of the protein-dynamic combinations (supplementary Figure S10).

This revealed resistance most frequently converges on the activation of CDK2. Mechanisms driving this activation include: the promotion of cyclin E synthesis, reduction in cyclin E degradation, promotion of the formation of the CDK2cyclin E complex and suppression of p21 and p27 mediated CDK2 inhibition. There were also moderate contributions to resistance spread throughout both upstream mitogenic pathways. The contribution of the ER to resistance seemed to be context dependent. A strong deactivation rate of the ER contributed to resistance quite frequently, yet paradoxically strong activation of the ER also somewhat frequently contributed.

The above results suggest three things. First, while there is a moderately large number of network states that can drive resistance to dual CDK4/6-ER inhibition, there is very likely only a finite number of network states that cancer can exploit to overcome drug perturbation. Second, the moderately large number of network states that facilitate resistance and the observation that any given resistance dynamic can be driven by many different mechanisms highlight the low likelihood of finding broadly-applicable resistance mechanisms. Third, some protein nodes and interactions contribute far more frequently to resistance than others. It is quite possible that these frequencies are related to the likelihood and tumorigenicity of mutations affecting these nodes/interactions.

The literature provides strong support for the predictions made by our MDN analysis

Monoculture cell populations demonstrate a significant degree of heterogeneous signalling dynamics in response to CDK4/6 inhibitors

This study provides a theoretical framework that connects heterogeneity, protein signalling dynamics and adaptive resistance mechanisms. To validate this framework, we analysed the literature to pinpoint biological evidence for a number of our key predictions. First, we wanted to observe if and how much drug-response signalling heterogeneity exists in cellular monocultures. Second, we wanted to find evidence that cells possess differential resistance-mediating subnetworks. And third, we wanted to determine if the resistance-mediating subnetworks we identified agree with known resistance mechanisms.

Recently, Xe et al. constructed a novel reporter system that enabled them to investigate the activity of CDK4/6 and CDK2 in individual MCF10A cells in response to three clinically approved CDK4/6 inhibitors: abemaciclib, palbociclib and ribociclib [55]. Plotting the single-cell time-course data for each drug and protein revealed that at the single-cell level there is a wide variety of protein signalling dynamics (Figure 7A). We then subjected this time-course data to our category analysis to calculate the distribution of each protein’s drug response dynamics (Figure 7B). We found that at the level of CDK4/6, the direct target of the three drugs, there was a moderate degree of signalling heterogeneity. Additionally, even the most potent drug (abemaciclib) failed to induce sustained suppression in approximately 10% of the cell population. The activity levels of CDK2, however, showed a much greater degree of signalling heterogeneity; and all three of the drugs only managed to strongly suppress CDK2 activity in a sustained manner in 40% of the population (Figure 7B). We further observed that there were many individual cells that actually demonstrated increasing CDK2 activity in response to CDK4/6 inhibition, a particularly unintuitive result, predicted by our MDN analysis. We then extracted a random sample of model instances from our MDN analysis (parametric variation) and plotted the drug response dynamics of CDK46cycD and CDK2cycE, our model’s equivalent active forms of CDK4/6 and CDK2 (Figure 7C). The resemblance between the distribution of dynamics produced by our model simulations and the single-cell data is striking, both showing very similar proportions of resistant and sensitive dynamics.

Figure 7.

Finally, we examined how the distribution of specific dynamics changed when comparing CDK4/6 activity with CDK2 activity. Figure 7D shows that the single-cell data and our MDN analysis exhibit the same general trends: there are more resistant dynamics at the level of CDK2 than at the level of CDK4/6; there are fewer sensitive dynamics at the level of CDK2 than CDK4/6; and there are more increasing dynamics at the level of CDK2. Perhaps most unintuitively, there are fewer rebound dynamics at the level of CDK2 in the single cell data, which was also predicted by our MDN analysis. Overall, these single-cell analyses support the model predictions of drug perturbation dynamics.

The identified resistance-mediating core subnetworks align with known resistance mechanisms

Previous studies have identified a wealth of mechanisms responsible for driving resistance to CDK4/6 and ER inhibitors. Focusing on mechanisms that have overlap with our model, the vast majority of identified mechanisms align quite well with the predictions made by our MDN analysis, validating the ability of this technique to identify a full spectrum of resistance-driving mechanisms. Moreover, our analysis also makes predictions about additional protein nodes/interactions that have not yet been identified by the literature.

Known mechanisms of resistance to CDK4/6 inhibitors include: Rb loss/mutations [56, 57], CDK4/6 overexpression [58–60], loss of ER expression [58, 61, 62], loss of PTEN [63], reduced expression/activity of p21 and p27 [64–66], activation of mTOR complexes [67, 68], activation of PI3K and PI3K signalling pathways [61, 69], activation of PDK1 [70], upregulation of FGFR [71, 72], AKT amplification and/or over-activation [73], E2F amplification and/or over-activation [74], MAPK pathway activation [75–77], c-Myc activity [78, 79], cyclin E overexpression and increased CDK2-cyclin E complex formation [43, 44, 80, 81]. Known mechanisms of resistance to ER inhibitors include mutations that induce the constitutive activation of ER [82, 83], loss of ER activity [84, 85], activation of PI3K and MAPK signalling [86, 87], overexpression of c-Myc [88, 89] and loss of p21 and p27 [90, 91].

The preceding list represents a wealth of known resistance mechanisms that are widely distributed throughout the ECC network. Every single one of the above mechanisms is represented in at least one of our identified resistance-driving subnetworks. To clearly demonstrate how well our MDN analysis captured the known resistance mechanisms we created a table containing the known resistance mechanisms and the respective resistance-driving parameters identified by our MDN analysis (supplementary Table S3). The ability of MDN to largely capture the known resistance mechanisms highlights its potential usefulness in predicting resistance mechanisms to novel drugs and drug targets. Apart from capturing the known resistance mechanisms, our analysis really extends the observation that resistance can emerge from many different nodes within the ECC network and further adds support to the idea that it is critical to treat each patient individually.

Discussion

Despite incredible advances in our understanding of cancer development and the paradigm shift towards targeted therapies, cancer remains a major cause of death worldwide. Although it is clear that resistance to current cancer therapies is a major barrier to achieving successful outcomes in cancer patients, we are still unable to reliably predict why and when resistance will occur during cancer treatment. In this study, we have focused on a systems-level characterisation of the relationship between heterogeneity, protein signalling dynamics and adaptive resistance.

The inability to reliably predict drug response can be partially explained by the complexity introduced by inter-patient and intra-tumoral heterogeneities; but exactly how heterogeneity translates to drug resistance is still poorly understood. Examples of biological phenomena that can alter the strength of protein-protein interactions include base-pair mutations and tissue/patient-specific post-translational modifications [92, 93]. Examples of cellular phenomena that can cause heterogeneity in the basal expression of proteins include copy number variation and the inherent stochasticity of gene expression [94, 95]. Throughout this investigation, we have argued that a major underlying link between heterogeneity and drug resistance is the ability of heterogeneity to influence protein signalling dynamics. Specifically, we have demonstrated that differences in the network state of a protein signalling network can drive shifts in the quantitative and qualitative response of a protein to a drug and can often facilitate the reversal of a sensitivity associated response to that of resistance.

Adaptive resistance is a drug-induced compensatory shift in signal transduction that enables cellular networks to partially, or even fully, restore a functional output, e.g. proliferation. These dynamic shifts are primarily driven by signalling feedback and crosstalk mechanisms that introduce non-linearity to the system and enable complex protein dynamics such as rebounding and biphasic responses.[18]. The intricate dynamics of these feedback mechanisms are regulated by the intricate interplay of protein interaction strengths (the parameters of the network) and basal protein expression levels (the initial conditions of the network). Such feedback mechanisms provide the system with critical features such as rapid response rates and protection from environmental noise, but can also be exploited by cancer cells to overcome drug perturbations, as observed in the phenomenon of adaptive resistance [18, 23].

In this study we developed a novel modelling technique, Meta-Dynamic Network modelling, to demonstrate how heterogeneity in protein interactions and protein expression translates to changes in protein signalling dynamics. By simulating many thousands of variations of a base network model we can gain insight into not only the types of protein dynamics that are possible for each protein species, but also gain an understanding of how frequently, or how likely these dynamics are to emerge. The generation of meta dynamics within a network can be conceptualized as the production of a map, offering insight into the specific drug responses of proteins, as well as the overall behaviour of the network. Furthermore, this approach allows for the identification of the multitude of mechanisms by which cancer can exploit network state to evade drug perturbation.

Throughout this investigation we modelled heterogeneity in the network properties of the ECC network to determine if there were any network states capable of facilitating adaptive resistance. We found that the ECC network was robustly capable of facilitating adaptive resistance dynamics in the face of broad heterogeneity. Even though all of the key output proteins (i.e. pRB, pppRB and CDK2cycE) frequently displayed drug-sensitive dynamics, the proportion of resistant dynamics was far higher than we expected (Figure 3B and Figure S1A). We also observed that simultaneous sensitivity across key downstream cell cycle regulators was also exceedingly rare. Both of these results were supported by single cell data, which showed that even within a supposed genetically identical monoculture there was a large proportion of cells demonstrating adaptive-resistance associated signalling dynamics for both CDKs, and that robust suppression at the level of CDK4/6 did not necessarily translate to robust suppression of CDK2 activity.

The surprisingly large number of network states that can facilitate resistance may explain the wide range of responses and broad insensitivity to CDK4/6 inhibition seen across hundreds of different cell lines as seen in the Genomics of Drug Screening in Cancer (GDSC) database [42]. The large number of possible resistant network states means that any changes to a cell’s network state have a statistically high likelihood of causing it to become resistant. Thus, across a large enough population of cells, possessing a broad enough range of intra-cellular heterogeneity, we would expect CDK4/6 inhibitor resistance to be relatively common.

We also observed that increasing the number of model instances beyond a certain point did not result in any appreciable changes in the distribution of protein dynamics. Doubling the number of model instances from 100,000 to 200,000 resulted in less than a 1% shift across all 50 protein species (0.02% per species, supplementary Figure S12). This indicates that there is an upper limit to the degree to which parameter sets can influence the qualitative shape of a protein’s dynamic within a given network topology. Combining this with the observation that most proteins displayed a dominant dynamic, seen in Figure 3B and Figure S1A, it appears that protein dynamics are primarily regulated by a network’s topology, i.e. the edges between proteins and the classes of those edges i.e. the rate law used to model the interaction. Similarly, when we compared the effects of parameter variation with variation in the initial conditions, we found that parameter variation was significantly more capable of influencing protein dynamics. Overall, these results suggest there exists three tiers of protein-dynamic regulation: topology, parameters, and initial conditions, or alternatively, protein interaction classes, protein interaction strengths and protein expression.

In the healthy context, this arrangement of tiered regulation likely acts to ensure that cells can continue to function appropriately even after numerous small changes that cause shifts in protein interaction strengths and protein expression. The average cell experiences somewhere in the order of 10⁴ – 10⁵ DNA damage events per day [96] and even with an extremely successful repair rate the number of mutations inevitably accumulates and causes cells to diverge from their original state [97]. Our results show that even in the face of pretty extreme parametric heterogeneity, as long as the overall topology of a cell’s protein network remains intact the majority of cells within a population are still likely to continue to respond to network perturbations appropriately. However, the inherent plasticity of cellular systems, which enables them to adapt to a wide range of physiological conditions, also confers the potential for the emergence of drug-resistant subpopulations during cancer therapy.

A hierarchical clustering analysis of the absolute parameter values was performed to investigate the underlying patterns that give rise to drug resistance. The results of this analysis did not reveal any significant patterns (Figure S4). This suggests that knowledge of the strengths of interaction between individual protein pairs alone is not sufficient to predict resistance, and that adaptive resistance is an emergent property of the network. Supporting this, when we investigated the relative contributions of parameters to resistance, we were able to identify clusters of model instances that possessed common sub-networks driving the resistance-associated dynamics (Figure 6, Figure S8 and Figure S9A-I). This suggests that there is an upper limit to the adaptive mechanisms that can be exploited by cancer cells to circumvent a specific drug therapy. In the future it may be possible to design therapeutic strategies that target both the primary cancer-driving protein activity as well as consequential secondary, resistance-conferring activity, to robustly increase a treatments long-term efficacy.

We also attempted to validate the model prediction that differential core sub-networks can drive drug resistance (Figure 6) by utilizing the Cancer Dependency Map (DepMap) [98]. The DepMap database contains information about the degree to which gene knockouts affect proliferation in over 1000 cell lines. Note that cell lines and gene knockouts phenomenologically correspond to a series of biological model instances and targeted drug treatments, respectively. This allowed us to compare the effects of parameter knockdown on the dynamics of our key output proteins, across various model instances, to the dependency of cell lines on the genes that encode the same proteins. We extracted the dependency data for the genes of the proteins within the CDK4/6-centred network, selected only the genes that have a moderate to strong negative effect on proliferation, removed any genes upon which cells were universally dependent, and performed hierarchical clustering on the resulting gene-dependency matrix. Similar to the results shown in Figure 5, this analysis produced a heatmap that compares the overall gene-dependency ranking with cluster-specific gene dependencies (Figure S11). The heatmap produced by this analysis is not specific to CDK4/6 inhibitors and so it is not surprising that the subnetworks produced do not match those found in our more focused analysis. However, the fact that clusters in the DEPMAP data exist strongly suggests that different cell lines do indeed possess differential vulnerabilities to drug treatments and that cells do indeed possess differential resistance-mediating protein signalling subnetworks.

The identification of parameter patterns that can drive resistance dynamics also provides a link between adaptive resistance and acquired resistance. If a parameter influences an adaptive resistance mechanic, it seems very likely that mutations affecting that parameter can also confer acquired genetic resistance. In our analyses the phosphorylation rate of Rb by the CDK2-cyclin E complex is consistently associated with a rebounding dynamic. In the literature, the overexpression of cyclin E is frequently reported as facilitating resistance to CDK4/6 inhibitors. This demonstrates that in the future meta-dynamic analysis may be able to predict potential resistance mechanisms before they emerge during drug treatment.

Our work demonstrates drug resistance can be driven by individual mutations that strongly alter protein dynamics; however, it also demonstrates that resistance can be driven by the accumulation of the right collection of less significant mutations. Individually, a single mutation may not facilitate resistance, but together they can significantly alter protein dynamics and confer drug resistance that can enable it to survive until it acquires more potent resistance and cause disease relapse. It seems plausible that ‘lesser’, passenger mutations are far more likely to occur than more dominant driver mutations, but when the dominant mutations arise, the clones containing them expand. In this way the lesser mutations enable cancer to initially survive a drug treatment and provide a baseline from which cancer can evolve much more robust resistance mechanisms.

Cells are dynamic entities that are in constant flux, and this is exemplified by the protein signalling networks that drive cellular behaviours. Using a novel technique coined meta-dynamic network modelling, we have demonstrated how heterogeneity influences protein signalling dynamics, and thus how cancer can exploit heterogeneity to overcome a drug treatment. Although we have focused on the ECC network and CDK4/6 therapy as prime examples due to their clinical importance, our work outlines a broadly-applicable theoretical framework that links heterogeneity, protein signalling network dynamics and drug resistance. In the future, it may be possible to use MDN modelling to predict how cancer will escape drug treatments, allowing oncologists to respond accordingly, and robustly maintain cancer remission.

Acknowledgements

L.K.N was funded by a Victorian Cancer Agency Mid-Career Research Fellowship (MCRF18026); a Venture Grant from Cancer Council Victoria, Australia; and an Investigator Initiated Research grant from the National Breast Cancer Foundation and Love Your Sister, Australia (IIRS-20–094). A.H was This research is/was supported by an Australian Government Research Training Program (RTP) Scholarship.

Author contributions

Conceptualization: Anthony Hart, Lan K. Nguyen.

Formal analysis: Anthony Hart, Sung-Young Shin, Lan K. Nguyen.

Funding acquisition: Lan K. Nguyen.

Project administration: Lan K. Nguyen.

Supervision: Lan K. Nguyen.

Writing – original draft: Anthony Hart, Lan K. Nguyen, Sung-Young Shin.