Most cancer cells harbor multiple drivers whose epistasis and interactions with expression context clouds drug and drug combination sensitivity prediction. We constructed a mechanistic computational model that is context-tailored by omics data to capture regulation of stochastic proliferation and death by pan-cancer driver pathways. Simulations and experiments explore how the coordinated dynamics of RAF/MEK/ERK and PI-3K/AKT kinase activities in response to synergistic mitogen or drug combinations control cell fate in a specific cellular context. In this MCF10A cell context, simulations suggest that synergistic ERK and AKT inhibitor-induced death is likely mediated by BIM rather than BAD, which is supported by prior experimental studies. AKT dynamics explain S-phase entry synergy between EGF and insulin, but simulations suggest that stochastic ERK, and not AKT, dynamics seem to drive cell-to-cell proliferation variability, which in simulations is predictable from pre-stimulus fluctuations in C-Raf/B-Raf levels. Simulations suggest MEK alteration negligibly influences transformation, consistent with clinical data. Tailoring the model to an alternate cell expression and mutation context, a glioma cell line, allows prediction of increased sensitivity of cell death to AKT inhibition. Our model mechanistically interprets context-specific landscapes between driver pathways and cell fates, providing a framework for designing more rational cancer combination therapy.
Oncogene-targeted small molecule kinase inhibitors, like imatinib for BCR-ABL [1], have transformed chemotherapy. However, such precision medicine approaches are not always efficacious. In some cases, mutation-matched patients do not respond to the drug [2], or alternatively, resistance stochastically develops [3]. Monotherapy can even activate the target pathway, depending on cellular context [4]. Combination therapy is a logical and clinically-proven path forward [5], but rationalizing even just the choice of combinations from among the at least 28 FDA-approved [6] targeted small molecule kinase inhibitors, notwithstanding important questions related to the dozens of traditional chemotherapeutics, monoclonal antibodies, dose, timing and sequence, remains challenging for basic and clinical research. Approaches are needed that consider quantitative, dynamic, and stochastic properties of cancer cells given a context. Computational modeling can help fulfill this need, since simulations often are much quicker and less expensive than the explosion of experimental conditions one would need to assay the drug combination space. Some of the first approaches define transcriptomic signatures that specify tumor subtypes and suggest drug vulnerabilities [7]. Big data and bioinformatic statistical approaches are at the forefront of predicting drug sensitivity, with penalized regression and other machine learning methods linking mutation or expression biomarkers to drug sensitivity [8,9]. Such statistical modeling approaches largely cannot extrapolate predictions into untrained regimes with confidence. However, this is often a major task-of-interest; for example, drug combination predictions usually only have data for drugs used alone. They also usually cannot account for dose and dynamics, central to combinations. Alternatively, mechanistic computational models based on physicochemical representations of cell signaling have an inherent ability to account for dose and dynamics. They also have greater potential for extrapolation, since they represent physical processes. Multiple mechanistic models exist for almost every major pan-cancer driver pathway identified by The Cancer Genome Atlas (TCGA): receptor tyrosine kinases (RTKs), Ras/ERK, PI3K/AKT, Rb/CDK, and p53/MDM2 [10]. An integrative model that accounts a multi-driver context has not yet been built, but will likely be needed to address drug combination prediction problems with such modeling approaches. Here, we constructed the first mechanistic mathematical model integrating these commonly mutated signaling pathways. We use the MCF10A cell line–a non-transformed mammalian cell line with predictable phenotypic behaviors. We tailor the model to genomic, transcriptomic, and proteomic data from MCF10A cells and train the model using a wealth of literature resources as well as our own microwestern blot and flow cytometry experiments to refine biochemical parameters and phenotypic predictions. We use the model to explore stochastic proliferation and death responses to pro- and anti-cancer perturbations. This mechanistic, biological context-tailored mathematical model depicting several major cancer signaling pathways allows us to probe the mechanisms that underlie how noisy signaling dynamics drive stochastic proliferation and death fates in response to various perturbations, and gain insight into their biochemical mechanisms.
An overarching goal of this work is a foundation that integrates canonical biochemical knowledge of cancer driver pathways with cell-specific context to make quantitative predictions about how drugs, microenvironment signals, or their varied combinations influence cell proliferation and death (Figs 1A and S1; S1 and S2 Tables). Three main considerations guided us. First is scope. TCGA identified pan-cancer driver pathways [10]: multiple RTKs; proliferation and growth (Ras/Raf/MEK/ERK; PI-3K/PTEN/AKT/mTOR); cell cycle; DNA damage. These interface with expression and apoptosis. See Fig 1A for an overview of the model (see S1 Fig for a detailed kinetic scheme) and lists of proteins included in each submodel. Second is inter- and intra-tumoral heterogeneity. Experimental accessibility defines each, which predominantly consists of copy number variation, mRNA expression and mutations, which the foundation integrates. Intra-tumoral heterogeneity is genetic [11] and non-genetic [12], so we require accounting for such noise. Third is formalism. Signaling, cell fate, and drug action are all inherently quantitative, stochastic and dynamic phenomena. Predictions for a multitude of contexts and drugs, many of which would not have yet been calibrated for, is facilitated by mechanistic as opposed to empirical models [13], so we take a mechanistic chemical kinetics approach. We present the model as follows (Fig 1B). In “unit testing” we adapt and develop submodels for individual pathways (see Supplement) in isolation, requiring that each satisfies a set of experimentally observed biological behavior. We focus on a widely-studied non-transformed cell line, MCF10A [14] to leverage literature data and establish a baseline for “normal” cell behavior before describing transformed contexts with genetic instability and multiple poorly understood mutations. We then connect these submodels together via established signaling mechanisms to perform “integrated unit testing”, which is iterative. Finally, we analyze the model to reason about emergent signaling mechanisms underlying biological observations. Each submodel unit is now constrained by a variety of biological observations and data, so we move to integrated unit testing. We do not imply that submodels are “correct”. Rather, the presented characterization simply increases confidence that submodels reproduce expected cellular behavior, which of course is predicated upon the considered data and assumptions. We first evaluate integrated behavior in the reference state—serum starvation. The chemotherapeutic etoposide induces DNA damage. We treated serum-starved MCF10A cells with etoposide for 24 or 48 hours, but found negligible induction of cell death over control (Fig 4A-top). This is consistent with etoposide-induced cell death being amplified by DNA replication [42], which we also observed (EGF + insulin, Fig 4A-bottom). We postulated a pharmacodynamic model where etoposide-induced DNA damage is strongly dependent on S-phase cyclins (see S1 Supplemental Methods) that reproduced the observed cell death results (Fig 4A–4C). Very few simulated cells without growth factors show a p53 response or active cycling (Cyclin A), and therefore very little cell death (cleaved PARP). With EGF and insulin most simulated cells show a robust and sustained p53 pulsing response, which leads to cell cycle arrest (through p21 and Chk1) and cell death. We noticed that etoposide-induced cell death in the presence of EGF and insulin plateaus over time in experiments and simulations at around 60%, as opposed to increased cell killing over time. We reasoned that this was due to p21/Chk1-induced cycle arrest. In simulations where p53 has limited ability to upregulate p21, etoposide is predicted to cause greater cell death (Fig 4D). This result is supported by studies showing that cells with mutated p21 are more sensitive to DNA damaging agents [43–45]. Furthermore, several cancer cell lines from the Genomics of Drug Sensitivity in Cancer Project with p21 mutations exhibit increased sensitivity to etoposide [9]. Thus, our model recapitulates the complex dynamic response to a common chemotherapy and predicts features of its context-dependent efficacy. The mechanisms by which oncogene-targeted drugs cause cell death through loss of survival signaling are not as well understood as their cytostatic effects. We measured cell death dynamics (assayed via Annexin-V/PI flow cytometry) in response to small molecule inhibitors of ERK and AKT pathways (downstream mediators of many oncogene-targeted drugs) alone or in combination, with serum-starvation or EGF + insulin (Fig 5A). As before, serum-starvation induces cell death, which is attenuated by EGF + insulin. Inhibition of either ERK or the AKT pathway alone increases cell death 72 hours post treatment, but negligibly at 48 hours. Combining the two inhibitors revealed synergy in cell death induction at 48 hours, and much more cell death at 72 hours. The integrated model captures these cell death dynamics reasonably well (Fig 5A—compare far right to far left). We used the integrated model to reason about what mechanisms may underlie the synergy between and overall impacts of ERK and AKT inhibition for cell death. Both ERK and AKT can phosphorylate and inactivate the pro-apoptotic function of BAD, ERK can phosphorylate BIM causing similar inhibition of pro-apoptotic function, and AKT-mediated deactivation of FOXO downregulates BIM levels (Fig 5B). Simulations where either BAD- (blue arrows) or BIM-dependent (red arrows) mechanisms were computationally knocked-out predicted BIM to account for most ERK and AKT inhibition effects (Fig 5C). Independent experimental data corroborate this prediction; MCF10A cell detachment-induced death, which depends on such survival signaling mechanisms, has been observed to be predominantly controlled by BIM rather than BAD [46,47]. The model suggests several contributing mechanisms. One reason is expression levels; BAD is lowly expressed compared to anti-apoptotic BCL2 family proteins (~14,000 vs. ~61,000 molecules per cell), and thus is unlikely to exert significant control over cell death. Another is the cascade network structure for BIM regulation. If only AKT is inhibited, phosphorylation by ERK blocks the accumulation of active BIM. If only ERK is inhibited, the amount of BIM in the cell, even if completely unphosphorylated, is largely insufficient to drive cell death in the absence of further strong death signals. However, when both pathways are inhibited, the production and accumulation of active BIM proceeds uninhibited, thus creating a potent apoptotic signal. Furthermore, the inhibition of these activities must be coordinated over relatively long periods of time, with duration of AKT inhibition having enough overlap with that of ERK to cause both the accumulation and activation of BIM. Why do some cells die early and others die late following ERK and AKT inhibition? We used simulation data to gain insight into this question. First, we asked whether initial total protein levels (sum across all forms of a protein) were correlated with time-to-death after ERK and AKT inhibition. In line with conclusions from prior studies focused on extrinsic death pathways [12], time-to-death was not significantly correlated with the initial levels of any single protein (Fig 5D). This suggests that stochastic processes post-inhibitor treatment largely dictate variable cell death fate. We performed lasso regression using time-averaged protein levels as predictors of early or late death (pre/post 40 hours) over 400 individual cell simulations. This analysis identified time-averaged BIM and BCL2 of the top explanatory variables. We note that the temporal dependence of explanatory variables reveals a dynamic signature of the pathway. Next, we trained a support vector machine (200 cells for training, 200 for validation) and found that average BIM and BCL2 levels over a 40-hour time course (or as long as a cell lived) are highly predictive for death timing (Fig 5E). Initial levels, or even those averaged over the first 8 hours post-stimulus, were moderately predictive at best. As one might expect based on the respective pro- and anti-apoptotic functions of BIM and BCL2, respectively, stochastic expression time courses that tend towards high BCL2 and low BIM are the most protective (Fig 5F). Interestingly, these predictors of stochastic death phenomena were not equally applicable to different death stimuli such as the extrinsic ligand TRAIL (S6A Fig), thought to depend more on stochastic activation dynamics of Caspase 8 [12], which highlights the treatment-specificity of stochastic cell fate determinants. Overall, these results suggest that intrinsic death of individual cells induced by inhibition of survival signals may be inherently unpredictable prior to treatment. Such a claim is inherently difficult to prove conclusively, since there always might exist alternative measurable predictors and/or methods of analysis that render such prediction possible. Moreover, in our specific case, we have not yet incorporated detailed stochastic mechanisms of protein degradation, which when included may provide additional informative co-variates. EGF and insulin regulate MCF10A proliferation. We measured (BrdU incorporation/flow cytometry) how these two ligands influence cell cycle progression of serum-starved MCF10A cells and how the ERK and AKT pathways were involved (Fig 6A). Insulin induces negligible cell cycle entry on its own but synergizes with EGF (Fig 6A), which is also seen on the level of cyclin D expression (Fig 6C). The ERK and AKT pathways are both essential to drive EGF + insulin-induced cell cycle progression (Fig 6A). After tuning the dependencies of transcriptional processes downstream of the ERK and AKT pathways, the integrated model reproduces these data (Fig 6B and 6C). What underlies the synergy between EGF and insulin with regards to cell cycle progression? Prior work has suggested that the dynamics of the ERK pathway can determine cell proliferation fate [48]. However, data and simulations suggest that ERK pathway activation dynamics are essentially identical for EGF or EGF + insulin (Fig 6D). Insulin is a very strong activator of AKT but not ERK signaling (Fig 6D). When cells are treated with EGF + insulin, there are negligible differences in acute (<30 min) AKT activation dynamics, and significant differences only develop over long (hours) time scales. Thus, the mitogenic synergy between EGF and insulin seems associated with prolonged AKT pathway activity over long time scales. Simulations suggest that increasing mitogen-induced cyclin D expression cannot be accomplished by only increasing the activity of one pathway, but rather roughly equal amounts of both time-integrated ERK and AKT activity are needed to drive more robust mitogenic responses (Fig 6E). This interpretation cannot rule out coincident factors such as increased insulin-induced glucose uptake, but does suggest that the coordinated activation of both ERK and AKT for several hours post-mitogen treatment seem important to drive cell cycle entry, rather than acute activation dynamics following mitogenic stimulus. These observations suggested to us that stochastic differences from cell-to-cell in AKT activity dynamics might be predictive for EGF + insulin-induced cell cycle entry. Surprisingly, in simulation results, neither the initial AKT activity nor dynamics post-stimulus contained significant discriminatory information related to subsequent cell cycle entry decisions of individual simulated cells (Fig 6F, bottom). Rather, the ERK activity dynamics were much more predictive of subsequent stochastic cell cycle entry (Fig 6F, top). This is consistent with prior live-cell imaging work that showed a strong correlation between time-integrated ERK activity and resultant S-phase entry decisions [34], albeit in a setting of much higher cell confluence than our experiments here entail and with additional mitogenic factors present (hydrocortisone and cholera toxin). Why would AKT activity dynamics control synergy between EGF and insulin but not be predictive of stochastic cell fate? AP1 and cMyc control cyclin D expression in this model. AP1 (cFos-cJun) activity has a sharp response to time-integrated ERK activity, relative to the less sharp dependence of cMyc expression on AKT activity. This sharpness difference is due to positive autoregulation of cJun expression by AP1 [49], and is likely also influenced by the already high expression of cMyc in serum-starved MCF10A (~48,000 molecules/cell, compared to 0 and ~3,000 for cFos and cJun), which is one of the few MCF10A alterations (they also have CDKN2A/ARF loss, which is captured in our gene copy number data) [50]. We also wondered whether higher AKT activity lowers the threshold of ERK activity needed to cause cell cycling; however, we found this not to be the case (S7D Fig). Could stochastic cell cycle entry be predicted by protein abundance fluctuations prior to EGF + insulin treatment? We again applied lasso regression using initial conditions for total protein levels in individual simulated cells to identify explanatory variables for cell cycle entry, and then trained a support vector machine classifier. In contrast to predicting death in individual cells, here we found that initial CRaf and BRaf levels were fairly predictive of whether a cell would enter the cell cycle by 24 hours post EGF + insulin stimulus (Fig 6G). Initial total ERK and AKT (active + inactive) levels were not predictive at all (Fig 6G). We conclude that so long as AKT activity is prolonged, ERK activity fluctuations likely control stochastic cell proliferation fate in response to acute growth factor stimulus. These ERK activity fluctuations across a cell population seem heavily influenced by initial Raf protein levels. The model contains multiple oncogenes and tumor suppressors, and is calibrated to a non-transformed epithelial cell state. We wondered what simulations would predict the transformation potential to be for each gene. To test this, we systematically altered expression of each RTK, proliferation and growth submodel gene (10-fold up for oncogenes, 10-fold down for tumor suppressors), simulated proliferation response to EGF + insulin (emulate a growth factor-containing microenvironment), and ranked nodes by this metric (Fig 7A). Several expected high-ranking results are observed, such as Ras, Raf, ErbB2/HER2 and EGFR, those with well-established transforming viral analogs such as cFos, cJun and cCbl, and more recently recognized translation control with EIF4E and EIF4E-BP1 [51]. Gene products along the PI3K axis (PIK3Cα, PTEN, β-Catenin, cMyc) are lower ranked, which seems counter-intuitive since these are also frequent cancer drivers. However, MCF10A cells, while non-transformed, are immortal and fast growing in culture indefinitely, and one of the reasons (as stated above) is cMyc overexpression. The cMyc protein is downstream biochemically of PI3K/AKT, its overexpression occurs as a resistance mechanism in response to PI3K inhibitor therapy, and PI3K and cMyc alterations tend to be mutually exclusive [52–54]. Thus, we interpret these results as being consistent with the MCF10A context being naturally less reliant on these PI3K axis gene products to drive increased proliferation. One paradoxical result is the low ranking of MEK (MAP2K1/2), despite being functionally surrounded by highly ranked genes. This result is actually consistent with widespread clinical data that seldom report MEK alterations (~<5% TCGA pan-cancer in cBioPortal) as compared to Ras/Raf alterations, the explanation for which remains unclear. What mechanisms might explain MEK’s poor simulated transforming potential? In simulations, free C/B-Raf—those that are available to transmit signals from active Ras—are reduced, and most of this difference ends up sequestered by MEK itself, which seems to reduce active ERK (Fig 7B–7D). Simulations where MEK has artificially high affinity for MEK phenocopy this result. We interpret this as a consequence of stoichiometry and expression levels, since both C-Raf and B-Raf are relatively lowly expressed (~29,000 mpc and ~2,800 mpc, respectively) as compared to MEK (MAP2K1~147,000 mpc; MAP2K2~378,000 mpc), so although this sequestration effect might be relatively small, it can account for a non-negligible proportion of the total Raf pool. Thus, the prediction from these simulations is that MEK typically has poor transforming potential because increasing its levels it tends to sequester Raf from increased signaling as opposed to allowing for increased signal transmission. This of course does not fully explain observations related to the very low frequency of MEK activating mutations, which convolve protein structural robustness, and also potentially genome organization factors. Nevertheless, such a mechanism would not be possible to uncover from genomics or network model-based approaches because it is inherently based in quantitative biochemistry that our model captures from its formulation. Future work will have to perform a thorough genetic screen to experimentally validate these simulation-derived predictions regarding the transformation potential of single gene modifications. Thus far we have focused on tailoring the model to MCF10A cells, and have parameterized the model to reflect expression levels, protein dynamics, cell death and cell cycling percentages following various stimuli. Does our model possess the ability to be predictive of cell fate outcomes when tailored to an alternative cellular context? To test this, we used data from U87 cells, a glioma cell line. Flow cytometry experiments of U87 cells treated with a MEK inhibitor, an AKT inhibitor, or dual kinase inhibition indicated pronounced sensitivity of U87 cells to AKT, but not MEK, inhibition (Fig 8A); this result was not found for MCF10A cells, which were found to be minimally sensitive to either inhibitor used alone (Fig 5A, grey bars). We thus sought to test whether the model tailored to U87 expression data could predict this increased sensitivity to AKT inhibition. To test this, we quantified mRNA levels in serum-starved U87 cells and tailored the model with the same setup pipeline procedure as for the MCF10A-tailored model (Fig 2D). In addition, because U87 cells are PTEN-deficient [55–57], we set the PTEN translation rate constant to zero, essentially removing PTEN from the system. U87 cells are also thought to possess a mutation in CDKN2A; however, since the incorporation of this mutation would negligibly affect cell death responses to kinase inhibition, we excluded it from this analysis. However, future iterations of this model application to U87 should certainly incorporate CDKN2A mutation, in addition to more generally capturing the spectrum of mutations in this and other cell lines. We then simulated our U87 cell model under the above experimental conditions, recording percentage of cell death at 48 hours post stimulus (Fig 8B, red bars) and compared this to flow cytometry data (Fig 8A). Remarkably, the U87-tailored model predicted with reasonable accuracy cell death responses to the different stimuli, capturing the increased sensitivity to AKT inhibition compared to our MCF10A-tailored model simulations. We noted that the U87-tailored model is unable to capture basal cell death under the “No Stim” condition (Fig 8B). We suspect this highlights a limitation in the current initialization procedure, specifically in the tuning of basal caspase 8 activity to reflect balanced levels of intrinsic death signals in the absence of any stimuli. Future work must more carefully consider the phenotypic requirements of initialization across cell types. Simulated time courses suggest that the increased sensitivity to AKT inhibition in U87 cells may be driven by increased accumulation of FOXO (Fig 8C, left), and subsequently the pro-apoptotic BIM (Fig 8C, right), when compared to MCF10A cells. A proposed model-derived mechanism behind this observation lies in the faster degradation rate of phosphorylated FOXO [58–60]. In U87 cells, the presence of the PTEN-null mutation, in addition to higher basal AKT levels, results in higher basal AKT activity. This results in increased FOXO phosphorylation, which causes it to be more rapidly degraded. In turn, when AKT is inhibited and FOXO becomes unphosphorylated, its degradation is slowed and it accumulates, driving BIM upregulation and cell death. Future experimental work will have to test these mechanistic interpretations.
We developed the first pan-cancer driver network model of cell proliferation and death regulation that has a biological context mechanistically parameterized by multi-omics data. The model describes fractional cell kill and partial cytostasis in response to treatment with a variety of anti-cancer drugs, conditioned upon an expression and microenvironmental context. We first developed sets of unit test operations that build confidence in individual submodels, and then evaluated the ability of the model to reproduce cell biological behavior as an integrated whole, which we subsequently analyzed to gain insight into biochemical mechanisms of stochastic cell proliferation, death and transformation. We then applied the model to a glioma cell context, which was able to predict differential sensitivity to AKT inhibition. While this model is complex and larger than any previous such signaling model, some may argue it still does not account for enough biology. We have predominantly focused on a single non-transformed cell line with no mutations. Impeding future work here is that most mutations are not functionally well-understood. We do not incorporate detailed mechanisms of post-translational Myc regulation, pathways outside of pan-cancer scope, tumor metabolism, hypoxia signaling or immune response. We do not include VEGFR, which although an important receptor tyrosine kinase in cancer, predominantly plays a tissue level role in angiogenesis as opposed to a signaling role inside tumor cells themselves, the focus of this study [61,62]. We also have not considered cell-to-cell communication from physical contact, autocrine, or paracrine signaling—so called population effects—which will be important to consider moving forward. In our experiments, we minimized cell confluence to reduce such effects. Nevertheless, MCF10A are strongly adherent to one another, and the presented data, results and phenotypes, and therefore our model implicitly account for these population effects. Future experiments may focus on better teasing these components apart. Yet others may argue that our first model draft is too large to be useful. Symmetric critiques often reflect striking an appropriate balance, but both views have validity. Rule-based modeling [63] and algorithmic model building [64] could help modelers capture more biology. Open and FAIR [65] data that is well-annotated and available by computational query (e. g. LINCS consortium) will be needed to drive such modeling forward. We certainly do not claim that our present model is unique, neither in terms of structure nor parameterization. Several variants of submodels exist [22,29,30,66], as well as efforts to integrate them with expression data [67–70], most notable the whole cell model of mycoplasma [71]. These classes of models are well-known to exhibit parametric non-identifiability but yet are still able to make valid and precise biological predictions [72]. The fact that these models stand on the shoulders of decades of cell and molecular biology research and are formulated to the extent possible based on physicochemical principles facilitates predictive potential despite the inherent uncertainty. Increased ability to algorithmically compare model structure and predictions to a wide net of publicly available experimental data will help reduce model uncertainty. Regardless, model complexity is driven by the questions and data at hand, and here the ultimate goal of cancer precision medicine coupled with omics data for tailoring gives a rational basis for the model size and complexity. The fact that when we tuned the model to the U87 glioma cell context, it could predict increased sensitivity of cell death to AKT inhibition, suggests that the level of detail may be an appropriate foundation. To study stochastic cell fate with our model we combined penalized regression (identify predictive features) with machine learning (relate predictive features to cell fates). Like prior work [12], we found that ERK+AKT inhibitor-induced apoptosis was not predictable from initial protein levels. Rather, the couplet of time-averaged total BIM and BCL2 levels following death stimulus were highly predictive. Thus, apoptosis seems inherently unpredictable at the time of the death stimulus, and rather evolves stochastically over the treatment time course. It is difficult to prove such a claim conclusively, as a combination of different initial measurements with non-linear models could provide such predictive power. However, the results do suggest a strong stochastic element of cell death that is hard to predict from initial protein levels alone. To test the BIM and BCL2 dynamics predictions experimentally, the ability to track cells over time from the onset of perturbation to the final cell fate is needed. If one could engineer cell lines to express fluorescent protein tagged (yet still functional) copies of all BIM and BCL2 isoforms from each allele, and monitor their responses by live-cell imaging, then this prediction could be tested. However, construction of such lines remains a challenging proposition. In contrast to the inherent unpredictability of apoptosis based on initial protein levels, S-phase entry in simulations was predictable from initial C-Raf and B-Raf levels. Why is this the case for S-phase entry and not apoptosis? There are several potential factors in the model. First is the relatively low expression levels of both Raf isoforms, giving larger expression noise. Second is the bottleneck role of Raf in propagating signals from RTKs to S-phase entry, coupled with enzymatic signal amplification as opposed to stoichiometric titration as in apoptosis signaling. Third is the predominant stochastic control of S-phase entry by ERK signaling dynamics, as opposed to AKT signaling dynamics. This may be due in part to the downstream positive feedback of cJun expression, which amplifies noise. Small fluctuations in Raf levels therefore might be amplified to all or nothing cJun and consequently cell cycle entry responses. Our model can naturally account for drug dose, promiscuity, scheduling, kinetics and pharmacodynamics. These hallmarks of pharmacology are as yet difficult to rationalize with genomic approaches. One potential application of ours and similar models is improving upon genomic-based approaches to cancer precision medicine. For example, patient transcriptomic data could be used to tailor the model as described here, giving an in silico test bed to explore a plethora of single or combination therapy options, including dosing and timing. Given advances in patient-derived experimental systems from biopsies [73], perturbation data can be used to further refine patient-specific mechanistic models and test drug predictions prior to clinical administration. Such predictions would have engrained within important properties of clinical response such as the promiscuity of many targeted kinase inhibitors [74] and dynamic feedback mechanisms in cancer cells important for resistance [75]. Moreover, given inference of subclonal evolution patterns of a patient’s tumor [11], multiple instantiations of the model could be generated to find a cross-section of therapies, or perhaps even order of therapies, that could most efficiently control an (epi) genomically heterogeneous tumor. Alternatively, virtual clinical trials could be conducted for a particular drug or combination, given a cohort of patients defined by multi-omics data (e. g. transcriptomics). Simulation results could stratify patients for inclusion or exclusion in trials—an exceedingly important task given the niche efficacy of targeted drugs from single genetic biomarkers. Further simulations could optimize dosing and timing strategies. With now hundreds of FDA-approved anti-cancer therapeutics (including dozens of targeted drugs), and the realization that combination therapy is almost always required for clinical response, computational approaches such as this that help prioritize drugs, their doses, and their scheduling for particular patients will become increasingly important in drug development and personalized oncology.
Our model is the amalgamation of several already-published models from the literature into an integrated system of ordinary differential equations (ODEs). The model can be subdivided into six principal submodels: (1) receptor tyrosine kinase (RTK) [22,23,76–78], (2) proliferation and growth (which includes Ras-Raf-MAPK plus PI3K-AKT-mTOR pathways) [23,29,30], (3) cell cycle [40], (4) apoptosis [37], (5) DNA damage response [35], and (6) gene expression (the genes/proteins involved in each submodel are noted in main text Fig 1). The model is composed of ~1200 total species, which includes all genes, mRNAs, lipids, proteins, and post-translationally modified protein/protein complexes. There are a total of 141 genes across all submodels that can be reduced into 102 functionally unique proteins, which we term “protein conglomerates” (see S1 Table for mapping). These 102 immediate gene products can become post-translationally modified in various ways to make up the remaining ~1100 species in the model. We opted to represent the majority of reactions as elementary steps to the extent possible, using mass action kinetics. However, sometimes we employ Michaelis-Menten or Hill representations to capture events where the mechanisms are not sufficiently elucidated. Kinetic parameters are taken from original source models, previous models, experimental studies, or estimated to fit empirical observations (see S2 Table for rate laws and rate constant sources). See Supplementary Methods for a more detailed description of the model structure and formulation.
Cancer is a complex and diverse disease. Two people with the same cancer type often respond differently to the same treatment. These differences are primarily driven by the fact that two type-matched tumors can possess distinct sets of mutations and gene expression profiles, provoking differential sensitivity to drugs. Over the past few decades, we have seen a shift away from more broadly cytotoxic drugs to more targeted molecules therapies; but how to match a patient with a specific drug or drug cocktail remains a difficult problem. Here, we build a mechanistic ordinary differential equation model describing the interactions between commonly mutated pan-cancer signaling pathways—receptor tyrosine kinases, Ras/RAF/ERK, PI3K/AKT, mTOR, cell cycle, DNA damage, and apoptosis. We develop methods for how to tailor the model to multi-omics data from a specific biological context, devise a novel stochastic algorithm to induce non-genetic cell-to-cell fluctuations in mRNA and protein quantities over time, and train the model against a wealth of biochemical and cell fate data to gain insight into the systems-level, context-specific control of proliferation and death. One day, we hope models of this kind could be tailored to patient-derived tumor mRNA sequencing data and used to prioritize patient-specific drug regimens.
cell death
medicine and health sciences
diabetic endocrinology
cell cycle and cell division
cell processes
messenger rna
hormones
erk signaling cascade
insulin
endocrinology
biochemistry
rna
signal transduction
biochemical simulations
cell biology
nucleic acids
apoptosis
biology and life sciences
computational biology
cell signaling
signaling cascades
cyclins
A mechanistic pan-cancer pathway model informed by multi-omics data interprets stochastic cell fate responses to drugs and mitogens