- Open Access
A network modeling approach to elucidate drug resistance mechanisms and predict combinatorial drug treatments in breast cancer
Cancer Convergence volume 1, Article number: 5 (2017)
Mechanistic models of within-cell signal transduction networks can explain how these networks integrate internal and external inputs to give rise to the appropriate cellular response. These models can be fruitfully used in cancer cells, whose aberrant decision-making regarding their survival or death, proliferation or quiescence can be connected to errors in the state of nodes or edges of the signal transduction network.
Here we present a comprehensive network, and discrete dynamic model, of signal transduction in ER+ breast cancer based on the literature of ER+, HER2+, and PIK3CA-mutant breast cancers. The network model recapitulates known resistance mechanisms to PI3K inhibitors and suggests other possibilities for resistance. The model also reveals known and novel combinatorial interventions that are more effective than PI3K inhibition alone.
The use of a logic-based, discrete dynamic model enables the identification of results that are mainly due to the organization of the signaling network, and those that also depend on the kinetics of individual events. Network-based models such as this will play an increasing role in the rational design of high-order therapeutic combinations.
Decades of cancer research and clinical practice have showed that durable treatment of metastatic solid tumors is limited by the acquisition of resistance to the treatment (Holohan et al., 2013; Garraway & Jänne, 2012; Cree & Charlton, 2017). Attaining durable control of these tumors will likely require therapeutic combinations; i.e. combinations of drugs that target different key pathways within cancer cells. Our current knowledge of drug resistance mechanisms is based on resistance to single-agent treatments in cancer models and patients. The effective drug combinations employed in the clinic today, such as the ones used in chemotherapies and other notable success stories, have been mainly derived through empirical testing and following many failures. Yet the prediction of drug resistance mechanisms and design of therapeutic combinations based on scientific rationales is still an unmet need. The methods to reach these goals will have to take into account the genomic and phenotypic diversity of tumors, the variety of resistance mechanisms, and the intrinsically combinatorial nature of the problem (Higgins & Baselga, 2011; Friedman et al., 2015; Johannessen & Boehm, 2017; Meric-Bernstam & Mills, 2012). This makes the currently used strategies ineffective and calls for new approaches that fall under the broad umbrella of the systems biology paradigm (Werner et al., 2014; Archer et al., 2016).
Mechanistic network models of the signal transduction pathways underlying cancer cells are one of the pillars of systems biology research because of their ability to explain how these signaling cascades integrate internal and external inputs to give rise to a cellular response (Kumar Jolly & Levine, 2017; Tyson et al., 2011; Aldridge et al., 2006; Wang et al., 2012; Alon, 2006; Zhang et al., 2014; Tian et al., 2017). These properties make mechanistic network models ideally suited to approach the problems of identifying drug resistance mechanism and designing effective hypothesis-based drug combinations. In particular, we propose using a subtype of network models, known as discrete dynamic models, which have been shown to reproduce the qualitative behavior of cancer signaling networks and are constructed solely from the regulatory interactions among the signaling proteins and the combinatorial effect of these regulatory interactions (e.g. positive or negative, additive or multiplicative) (Wang et al., 2012; Morris et al., 2010; Steinway et al., 2014; Udyavar et al., 2017; Méndez-López et al., 2017; Collombet et al., 2017).
Network models of signal transduction pathways and discrete dynamics
A signal transduction pathway consists of enzymes (e.g. kinases and phosphatases), adaptors, and signaling molecules that integrate extracellular and intracellular information and relay it to the transcription factors responsible for the required cellular response. Signal transduction pathways can be represented as a network, where each network node denotes an element of the signaling cascade (e.g., a signaling protein) and a directed edge between two nodes means that the first node regulates the activity of the second (target) node.
As an example, consider the simplified version of signaling through receptor tyrosine kinases (RTKs) shown in Fig. 1. In signaling through RTKs, binding of growth factors to the extracellular domain of RTKs induces a conformational change in the RTK, which promotes the recruitment and binding of several signaling proteins and kinases to its intracellular domain. Among the recruited signaling proteins are RAS and PI3K, which are activated by the RTK, and recruit other signaling molecules. RAS activates the kinase BRAF, which phosphorylates and activates the MAPK cascade (MEK/ERK), which then elicits a transcriptional response. Similarly, PI3K phosphorylates the phospholipid PIP2 into PIP3, which in turn leads to the phosphorylation of the kinase AKT, which then activates several transcription factors. Fig. 1 shows the directed interactions involved in the described sequence of events in RTK signaling, and additionally, a directed interaction between RAS and PI3K that represents the RTK-independent activation of PI3K by RAS.
A network representation of a signaling pathway, like the one in Fig. 1, can be converted into a discrete dynamic network model by assigning a state variable σi and a regulatory function fi to each node i. Each state variable σi can take a discrete number of states which denote the level of activity of the signaling element represented by node i, and where each level of activity is defined by its regulatory effect on the state variables of its target nodes. The existing experimental evidence on the number of protein conformations or post-translational modifications that yield different outcomes points to the sufficiency of assuming a small number of states, e.g. two or three (Kapuy et al., 2009; Burra et al., 2009). Each regulatory function fi, which can be represented using the logical operators OR, AND, and NOT, encodes the combinatorial effect on σi of the directed interactions acting on node i and thus depends on the state of the regulators of i.
As an example, we convert the network in Fig. 1a into the simplest type of discrete dynamic model, a logical (or Boolean) dynamic model, in which each node state variable σi can have two states: ON (active) or OFF (inactive). We note that the OFF state does not mean the complete absence of activity but a level of activity that is not sufficient to regulate target nodes. For the regulatory functions of this model we use the following logical rules
which are mathematical statements of the transmission of information between the elements in RTK signaling. For example, fRTK = σGrowth Factor indicates that the RTK becomes active in the presence of external growth factors, and fPI3K = σRTK OR σRAS indicates that PI3K becomes active in the presence of either active RTK or active RAS. Another example is fTranscription Factors = σAKT AND σMEK/ERK, which indicates that the activation of the transcription factors we are considering (activation that may include transcriptional as well as post-translational regulation) requires both active AKT and active MEK/ERK. These latter rules are consistent with signaling in certain types of lung cancer (Castellano & Downward, 2011; Lim & Counter, 2005).
In addition to regulatory functions like those in Eq. (1), a logical model must also specify how the node state variables change with time based on these functions, that is, we need to specify an updating scheme. Here we use the general asynchronous (GA) updating scheme (Steinway et al., 2014; Garg et al., 2008; Saadatpour et al., 2010), which updates the node state variables in discrete time units by the following two steps: (i) choosing one randomly selected node j at each time step t and updating its node state σj(t) by plugging the node states of the previous time step in the regulatory function fj (σj(t) = fj[Σ(t − 1)], where Σ(t) = (σ1(t), σ2(t), …, σn(t)) is the network state and encodes the state of all nodes at time t), and (ii) transferring the node state from the previous time step of the nodes not selected in step (i) (σi(t) = σi(t − 1), i ≠ j).
To illustrate the GA updating scheme, consider the network model in Fig. 1a, the regulatory function Eq. (1), and an initial node state Σ(t = 0) in which σGrowth Factor = ON and the rest of node states are OFF. Note that the rule fGrowth Factor = σGrowth Factor indicates that the state of the growth factor is sustained, which means that σGrowth Factor will stay in its initial state (σGrowth Factor = ON in this case). Three sequences of network states Σ(t) using GA updating, which we refer to as trajectories, are shown in Fig. 1b; note that, because we randomly choose one node at each time step, there are many possible trajectories. In the left-most trajectory in Fig. 1b, MAPK signaling activates before PI3K signaling, while in the right-most trajectory the activation order is reversed. The middle trajectory shows both PI3K and MAPK signaling activating concurrently. In all three cases the long-term behavior is the same: PI3K and MAPK signaling and the target transcription factors are activated, that is, all the node states in Σ are ON. Discrete dynamic models always display patterns of long-term behavior, known as dynamical attractors (e.g. steady states, such as the state Σ with all nodes active in this example), which have been found to be identifiable with stable cell fates, cell states, or stable patterns of intracellular activity. The long-term behaviors (e.g. steady states) of discrete dynamic models can be identified not only by simulations, but also by alternative methods, including stable motif analysis (Zañudo & Albert, 2013; Zañudo & Albert, 2015), network reduction (Klamt et al., 2006; Saadatpour et al., 2013), and algebra-based methods (Veliz-Cuba et al., 2014).
The behavior of a population of cells governed by the same underlying intracellular network can be captured by the model by performing multiple simulations and interpreting each trajectory as the dynamics of a single cell. The simulated population can represent multiple types of heterogeneity by having a constitutive (in)activity of certain nodes in certain simulations, different starting states, or different kinetic parameters (this latter is implicitly captured by using stochastic update). To illustrate this latter type of heterogeneity, we use the network model of Fig. 1a and Eq. (1) with the initial state Σ(t = 0) of Fig. 1b and perform 10,000 simulations wherein we randomly select a node with equal probability and update that node only at each time step. To capture the population-level behavior, we use the average state of node i at time t in a set of trajectories to define a quantity called the activity of the node (ai(t)). The activity of each node is shown in Fig. 1c left. In addition, Fig. 1c right shows how the node activity changes when incorporating the biological constraint that signaling events are faster than transcriptional events, which we do by making the probability of choosing the Transcription Factors node be smaller than that of the rest of nodes. We choose the probability to be 5 times smaller for illustration purposes, even though the difference in time scales is significantly larger; ~10−3-1 s for signaling events and ~101–102 min for transcriptional and translational events (Milo et al., 2010; Milo & Phillips, 2015). Both timecourses in Fig. 1c show how the network elements downstream of GF are activated sequentially in the cell population (RTK first, followed by PI3K and RAS, followed by the elements downstream of them), and how PI3K signaling and MAPK signaling are activated at the same time, on average, in the cell population. The activity of the outcome node of this simple network, the node Transcription Factors, lags behind the activity of its two regulators, as it can only activate when both AKT and MEK/ERK are active. The assumed slower timescale (lower update probability) assumed in Fig. 1c further adds to the delay of the activation of Transcription Factors (TF).
A network model can also be used to simulate the effect of drug inhibition and to identify potential resistance mechanisms. For example, the addition of an RTK inhibitor in the model of Fig. 1a can be simulated by adding a node to the network denoting the RTK inhibitor (RTKi), setting the logical rule of RTK to fRTK = σGrowth Factor and not RTKi, and setting the state of the inhibitor to σRTKi = ON (either initially or at a certain time). Adding RTKi at time = 20 causes the reversal of the increase in the activity of both branches of signaling cascades (Fig. 2b). Ultimately, all the nodes downstream of the RTK become inactive in all the simulated cells, yielding a steady state identical to the initial state (t = 0 in Fig. 2b). A putative resistance mechanism can be evaluated in the model by changing the logical rule of a node suspected to be responsible for the observed resistance (e.g., an activating RAS mutation can be introduced by setting fRAS = ON), and testing its effect on the rest of the network. As shown on Fig. 2c, the activating RAS mutation leads to the reactivation of both signaling pathways despite the continued presence of the RTKi, and yields a steady state that differs from the steady state of Fig. 1b in the state of RTK only. In other words, activating RAS mutation causes resistance to RTK inhibitors in the model. Note that in this toy model no other activating mutation (except RAS activation) would result in resistance to RTK inhibitors because the outcome node TF requires both AKT and MAPK activity for its activation. The model can be used to identify the inhibitor combinations that are able to overcome resistance. For example, the introduction of a MEK inhibitor (MEKi) at time 20 stops the continued activation of the outcome node TF (after a time delay) and leads to its inactivation in all the simulations (Fig. 2d). Thus, although the PI3K pathway is still active under this condition, from the point of view of the outcome node Transcription Factors, the combined application of RTKi and MEKi has overcome the resistance.
Resistance mechanisms to PI3K inhibitors in breast cancer
The PI3K/AKT/mTOR signaling pathway is one of the most important regulatory pathways of cell growth and survival in healthy and cancerous cells, as evidenced by the finding that alterations in this pathway are one of the most common in human cancers (Mayer & Arteaga, 2016; Zhang et al., 2017). In particular, PI3KCA (the gene coding for the isoform α of the catalytic subunit of PI3K) is the most common altered gene in this pathway (mutated in ~15% of human cancers and having copy number amplifications in ~5% (Zhang et al., 2017; Zehir et al., 2017)) and is particularly important in the context of breast cancer (mutated in ~35% and copy number amplifications in ~5% (Ciriello et al., 2015; Koboldt et al., 2012; Stephens et al., 2012; Pereira et al., 2016)). The importance of PI3K in cancer has led to the development of drugs that target it. A variety of targeted drugs against PI3K are currently in clinical trials for breast cancer (Mayer & Arteaga, 2016); they range in specificity from dual PI3K/mTOR inhibitors, pan-PI3K inhibitors, to isoform specific inhibitors of PI3K (e.g. Alpelisib or BYL719, a p110-alpha/PIK3CA specific inhibitor).
As a result of the development of PI3K inhibitors, there has been an increased interest in investigating the resistance mechanisms to PI3K inhibitors in the context of breast cancer, and several studies have been done in this direction (Costa et al., 2015; Castel et al., 2016; Toska et al., 2017; Bosch et al., 2015; Elkabets et al., 2013; Le et al., 2016; Vora et al., 2014; Kodack et al., 2017; Zwang et al., 2017). These studies have elucidated several resistance mechanisms to PI3K inhibitors such as PIK3β signaling (an alternative PI3K isoform), HER3 (ERBB3) receptor activity (which is upstream of PI3K, and strongly activates the MAPK and PI3K pathway), mTORC1 signaling (which would otherwise be activated by the PI3K pathway), estrogen receptor (ER) transcriptional regulatory activity (which provides PI3K-independent means of promoting proliferation), and signaling through the PIM (PIM1, PIM2, and PIM3), SGK (SGK1, SGK2, and SGK3), and PDK1 protein kinases (which act independently of PI3K, and have functions similar to AKT). Importantly, these resistance mechanisms have been found to be dependent on each other in some cases. For example, evidence suggests that mTORC1 signaling in BYL719 resistant breast cancer cell lines HCC1954 and JIMT1 is a consequence of the higher activity of SGK and PDK1 in these cell lines, which is sufficient to activate mTORC1 through the phosphorylation of TSC2 by SGK. This dependence between resistance mechanisms suggests that an integrative approach that fully elucidates their joint and separate mechanism of action on cell signaling and cell survival is needed for a complete understanding and to make predictions of drug interventions that overcome the observed resistance mechanisms.
A network model of oncogenic signal transduction in ER+ breast cancer
We constructed a comprehensive discrete dynamic network model of signal transduction in ER+ breast cancer based on the literature of ER+, HER2+, and PIK3CA-mutant breast cancers (Fig. 3). The construction of the model follows a methodology that has been repeatedly used to model several other oncogenic and biological processes (Wang et al., 2012; Morris et al., 2010). In brief, we perform a comprehensive review of this literature and identify the pathways, molecular components, and interactions that have been mechanistically linked to the response or resistance to several targeted drugs. In particular, the model incorporates the findings of resistance studies in the context of PI3K inhibitors, mTORC inhibitors, AKT inhibitors, MAPK inhibitors, RTK inhibitors, CDK4/6 inhibitors, and ER inhibitors/degraders, the feedback mechanisms and adaptive cellular responses identified during these studies, and also includes the recent results of unbiased genome-wide screens for resistance mechanisms to PI3K inhibitors (Costa et al., 2015; Toska et al., 2017; Elkabets et al., 2013; Le et al., 2016; Vora et al., 2014; Kodack et al., 2017; Zwang et al., 2017; O’Reilly et al., 2006; Chandarlapaty et al., 2011; Zhang et al., 2011; Anderson et al., 2016; Miller et al., 2010; Nahta et al., 2005; Muellner et al., 2011; Turke et al., 2012; Will et al., 2014; Serra et al., 2011; Rodrik-Outmezguine et al., 2011; Vasudevan et al., 2009; Chakrabarty et al., 2012; Carracedo et al., 2008; Massarweh et al., 2008; Miller et al., 2011; Finn et al., 2009). Most of the signaling proteins and interactions in the network have been consistently identified as essential markers of the response or resistance to the selected targeted drugs in multiple cell lines, in vivo mouse models, and patient tumors. For signaling proteins and interactions that are less studied or that were more recently identified as key players in these pathways, we also use interactions and information from other cancers and biological processes, and when available, focus on the consensus between experiments done on canonical cell lines, in particular, MCF7, T47D, and MDA-MB-415 for ER+ breast cancer, and BT474, JIMT1, HCC1954, SKBR3, and HER2-overexpressing MCF7 for HER2+ breast cancer.
The model consists of 51 nodes (34 Boolean and 16 multi-state nodes), of which 13 are nodes with no regulators (source nodes) that encode the initial transcriptional state of the cell (e.g. ER, HER2) or the state of nodes which are not regulated by other elements in the model (e.g. PIM and mTORC2). The nodes correspond to proteins, transcripts (8 nodes) as well as the biological outcomes proliferation and apoptosis (see Additional File 1). The edges correspond to transcriptional regulation, epigenetic mechanisms, post-translational and signaling processes. The model incorporates elements of the main signaling pathways involved in breast cancer: RTK signaling (e.g. IGF1R and HER2/HER3), PI3K signaling (e.g. PI3K and PTEN), MAPK signaling (e.g. RAS and MAPK), AKT signaling (e.g. AKT, PDK1, and FOXO3), mTORC1 signaling (e.g. mTORC1, TSC, and S6K), and ER signaling (e.g. ESR1 and MYC). The model incorporates multiple negative feedback loops through which the PI3K/AKT pathway leads to the negative regulation of RTK signaling. These six pathways converge in the survival signaling proteins that control apoptosis (e.g. BIM, BAD, and MCL1) and proliferation (e.g. cyclins, RB, and E2F, which form a positive feedback loop). The model describes two biological outcomes with multi-state nodes: Proliferation (a 4-state node) and Apoptosis (a 3-state node). In addition to the 51 nodes, we also include 7 nodes that denote inhibitors of specific targets of interestFootnote 1: Alpelisib or BYL719 (PI3K inhibitor – p110-alpha isoform specific), Ipatasertib (AKT inhibitor), Fulvestrant (ER inhibitor – selective estrogen receptor degrader (SERD)), Palbociclib (CDK4 and CDK6 inhibitor), Everolimus or Sirolimus (mTOR inhibitor), Trametinib (MEK1 and MEK2 inhibitor), and Neratinib (HER2 and EGFR inhibitor). To our knowledge, this the first comprehensive network model of its kind in breast cancer.
In Additional File 1 we indicate the full name of each network node and support each relationship and regulatory function with references. To construct the regulatory functions, we start with the assumption that regulators are independent from each other and that negative regulators are dominant. Then incorporate any available conditional knowledge (e.g. that two regulators need to work together in order to be effective); this information is usually related to the biology of the interaction and is distilled from the same literature source as the interaction. Each node variable is initially assumed to have two states (OFF/0 and ON/1), and additional states (e.g., 2) are added if justified by the current knowledge. As an illustrative example, here we explain the regulatory functions of AKT and ER_transcription. For simplicity, the state of each node is described using the node name, thus AKT stands for AKT = 1, PIP3 stands for PIP3 = 1, PIP3_2 stands for PIP3 = 2.
The regulatory function of AKT encodes the facts that PIP3 recruits AKT to the membrane (Castel et al., 2016; CURRIE et al., 1999) and that membrane-bound PDK1 and mTORC2 phosphorylate AKT (Castel et al., 2016; Alessi et al., 1997; Sarbassov et al., 2005). We assume that PIP3-mediated recruitment together with phosphorylation by either PDK1 or mTORC2 is sufficient for AKT activation. Additionally, fAKT encodes that the drug Ipatasertib inhibits AKT activity and that this inactivation can be overcome by a high level of PIP3 (denoted PIP3_2), an assumption that is consistent with the AKT-inhibitor literature (Chandarlapaty et al., 2011; Will et al., 2014). The model describes the gene encoding the estrogen receptor (ER) with two levels of activity (ESR1 and ESR1_2). Similarly, the transcriptional regulatory activity of ER has two levels as well (ER_transcription and ER_transcription_2). The ER_transcription rule indicates that baseline transcriptional regulatory activity of ER (ER_transcription = 1, given by fER_transcription) requires an ER+ cell and a baseline (ESR1 = 1) or upregulated (ESR1 = 2) expression of the ER transcription factor (the ESR1 gene codes for ER). Thus, the downregulation of ESR1 (ESR1 = 0, e.g., by the effect of the drug Fulvestrant) would result in below-threshold transcriptional regulatory activity of ER (ER_transcription = 0). Enhanced ER transcriptional regulatory activity (ER_transcription = 2, given by fER_transcription_2) requires high expression of ESR1, a KMT2D-mediated open chromatin state, and the participation of the co-activators FOXA1 and PBX1 (Toska et al., 2017; Bosch et al., 2015).
We focus on the context of ER+/HER2- breast cancer, which we encode in the model by setting the node ER to ON and HER2 and HER3_T to OFF. In addition, we start by considering a cell state in which the source nodes IGF1R_T and PBX1 are ON (IGF1R is a common RTK in breast cancer signaling, and the subscript T denotes the intrinsic transcript level of IGF1R; PBX1 is a co-factor required for ER-dependent transcription). The source nodes (i.e. nodes with no regulators) PTEN, SGK1_T, PIM1, PDK1_T, and mTORC2, which act as resistance mechanisms to PI3K inhibitors, are OFF, and the source nodes BIM_T and BCL2_T can be ON or OFF (BIM and BCL2 are pro- and anti-apoptotic proteins, respectively).
In the absence of drugs, the model recapitulates a cancerous state (see Additional File 1) in which RTK, PI3K, MAPK, AKT, mTORC1, and ER signaling are active, which results in high survivability: Proliferation is high (Proliferation = 3 or 4) and Apoptosis is low (Apoptosis = 0). In this cancerous state, high survivability is possible even if the apoptotic proteins are active: anti-apoptotic protein BCL2 can be either ON or OFF, and pro-apoptotic protein BIM can be ON as long as BCL2 is also ON to counteract its effect. Thus, this state corresponds to a set of six steady states: Proliferation = 3 or 4 (caused by E2F = 2 or 3), with BCL2 = BIM = OFF, BCL2 = BIM = ON, or BCL2 = ON and BIM = OFF. This indicates high survivability states that can be either primed (BIM = ON) or unprimed for cell death (BIM = OFF), a commonly observed feature of cancer cells (Sarosiek et al., 2017; Lee et al., 2012).
The network model recapitulates the response to PI3K inhibitors and predicts the degree of survivability of different resistance mechanisms
We simulate the effect of a PI3K inhibitor on a population of cancer cells by starting from a combination of steady states corresponding to the cancerous state and setting Alpelisib = ON at time = 2 and maintaining it until the end of the simulation. In order to simulate the dynamics of the network model, we use general asynchronous updating, categorize nodes into fast or slow depending on whether the node is activated by a (fast) signaling event or a (slow) transcriptional/translational event, and set the update probability of fast nodes to be 5 times higher than that of slow nodes. The resulting time course of node activities is shown in Fig. 4, where time is scaled so that the time unit is equal to the average time needed to update a slow node. The time course recapitulates the experimentally observed response to PI3K inhibitors, in which PI3K inhibition has a quick and direct attenuating effect on MAPK, AKT, and mTORC1 signaling, followed by the nuclear localization of FOXO3, which transcriptionally upregulates the transcription factor ER (coded by the gene ESR1) and the pioneer factor FOXA1, which increases ER transcriptional regulatory activity. The PI3K inhibition-induced fast signal transduction events converge with the slow transcriptional events triggered by cell signaling and regulate both apoptosis and proliferation. For example, AKT and mTORC1 are quickly inhibited following PI3K inhibition, which results in the dephosphorylation and activation the pro-apoptotic protein BAD, and in the attenuation of the translational machinery. Meanwhile, pro-apoptotic protein BIM is transcriptionally up-regulated by FOXO3, and cell cycle protein cyclin D is transcriptionally upregulated due to the increased ER transcriptional regulatory activity, both of which occur later in the response to PI3K inhibition. The end result is a marked decrease in survivability: an increase in apoptosis (from Apoptosis = 0 to Apoptosis = 2 or 3, depending on whether BCL2 was initially active) and a decrease in proliferation (from Proliferation = 3 to Proliferation = 1, due to the early downregulation of MAPK, AKT, and mTORC1 activity) followed by an increase (caused by the late upregulation of ER transcriptional activity) and then stabilization at Proliferation = 2. We summarize the apoptosis and proliferation propensity with the normalized and averaged values Apoptosisnorm and Proliferationnorm (Additional File 1), which in the current simulations take the initial values Apoptosisnorm = 0.00 and Proliferationnorm = 0.50, and the final values Apoptosisnorm = 0.70 and Proliferationnorm = 0.25.
We next test whether two recently discovered resistance mechanisms to PI3K inhibitors, PIM1/2/3 and SGK1/PDK1, increase survivability in response to PI3K inhibition in our network model. We start with an initial population of cells in the cancerous state and set either PIM = ON (which stands for any of the PIM family members) or PDK1 = SGK1_T = SGK1 = ON, and simulate the system as in the previous case (Alpelisib = ON at time = 2). The resulting time course of node activities is shown in Fig. 5b-c. Both PIM and SGK1 act as resistance mechanisms to PI3K inhibitors in the model, as evidenced by a decrease in Apoptosis (from Apoptosisnorm = 0.70 in case of PI3K inhibitors only to 0.00/0.25 in the PIM/SGK1 cases) and an increase/lack of change in Proliferation (from Proliferationnorm = 0.25 to 0.50/0.25 in the PIM/SGK1 cases). A closer look at the network and the interactions of PIM and SGK1 (Fig. 5a) shows that they share most of the downstream targets of AKT, and thus, can compensate for the loss of AKT activity due to PI3K inhibition. In particular, PIM shares four out of the six AKT targets in the model (PIM does not phosphorylate TSC nor KMT2D) while SGK1 shares two AKT targets. The fact that SGK1 does not regulate the activity of pro-apoptotic protein BAD and the cyclin dependent kinase inhibitors p21/p27 is the reason why the model predicts that the PIM proteins are a stronger resistance mechanism to PI3K inhibitors compared to PDK1/SGK1. We note that this prediction depends on the relative ability of PIM and SGK1 to phosphorylate their downstream targets. To illustrate this point, Fig. 5b bottom shows the resulting time course for PIM if it is 10% less efficient than AKT on its downstream targets, (which we implement by setting PIM = OFF with a probability of 10% at every time step). While fully effective PIM maintained the cancerous mTORC1, Apoptosisnorm and Proliferationnorm values despite PI3K inhibition, in the case of the 90% effective PIM there is a decrease in the average level of mTORC1 and Proliferationnorm and an increase in Apoptosisnorm, closer to the result obtained for SGK1.
The network model predicts MAPK, FOXO3, AKT, MYC, and cell cycle proteins as resistance mechanisms to PI3K inhibitors
In order to identify new resistance mechanism to PI3K inhibitors, we test every possible single and double node constitutive activation or inactivation (used in conjunction with PI3K inhibition), using an analogous procedure as in the case of PIM and SGK1/PDK1. Table 1 and Table 2 show the top node interventions that increase survivability (as measured by Apoptosisnorm and Proliferationnorm) compared to the control case of PI3K inhibition with no node interventions. For the case of single node interventions, the model recapitulates the known resistance mechanisms to PI3K inhibitors: PIM, SGK1, mTORC1 (mTORC1 = ON, TSC = OFF, PRAS40 = OFF, or translation = ON), and HER2/HER3 (HER2/HER3 = 2), which lead to a decrease in the apoptosis propensity and increase in the proliferation propensity. We identify several additional resistance mechanisms: MAPK (MAPK = 1 or 2, where level 2 is the state associated with HER2/HER3 activity), which partially or fully block apoptosis, AKT (AKT = ON), which fully blocks apoptosis and restores the PI3K-inhibitor-free proliferation levels, and FOXO3 (FOXO3 = OFF or FOXO3_Ub = ON), which leads to a decrease in both the apoptosis and proliferation propensity. We also identify several resistance mechanisms that involve cell cycle proteins, namely, cyclin E and CDK 2 (cycE/CDK2 = ON), p21/p27 (p21/p27 = OFF), E2F (E2F = 3), Rb (pRb = 3), or proteins of the mitochondrial apoptosis pathway, namely BIM (BIM = OFF), BAD (BAD = OFF), MCL1 (MCL1 = ON), and BCL2 (BCL2 = ON). Several of these resistance mechanisms are supported by experimental evidence (Rb, MCL1, and BAD) or are consistent with clinical observations (AKT) (Le et al., 2016; Vora et al., 2014; Zwang et al., 2017; Anderson et al., 2016).
For the case of double-node resistance mechanisms, and excluding those that target the apoptosis and proliferation pathways, we identify several new resistance mechanisms involving the AKT, MAPK, mTORC1, or ER pathways (Table 2). For example, MAPK = 1 combined with SGK1 = 1 fully blocks apoptosis in the model, and MAPK = 2 together with high ER activity (ER_transcription = 2 or MYC = 2) restores proliferation to its PI3K-inhibitor-free level (Proliferationnorm = 0.50) and fully blocks apoptosis. Other examples are MAPK = 1 combined with mTORC1-activating elements (mTORC1 = 1, TSC = 0, PRAS40 = 0), which restore proliferation to its original level (Proliferationnorm = 0.50) and also lower apoptosis significantly (Apoptosisnorm = 0.13), and FOXO3 = 1 together with MAPK = 2, which blocks apoptosis and restores proliferation (Proliferationnorm = 0.50).
The network model predicts that the inhibition of the MYC-CDK4/6 axis of cell-cycle regulation and of mTORC1 synergizes with PI3K inhibitors
We next asked what single or combinatorial interventions would further sensitize cells to PI3K inhibition, i.e. yield an increased apoptosis propensity or decreased proliferation propensity compared to PI3K inhibition alone. Table 3 shows the top interventions that synergize with PI3K inhibition, which in this case are all single-node interventions. Interventions that involve the inhibition of ER activity (e.g. Fulvestrant = 1, ER_transcription = 0, FOXA1 = 0, PBX1 = 0, KMT2D = 0) have a high anti-proliferative and apoptotic effect (Proliferationnorm = 0.00 − 0.13 and Apoptosisnorm = 0.83). Indeed, ER activity is up-regulated in response to PI3K inhibition and attenuates drug response. The synergistic effect of PI3K and ER inhibition has been previously reported (Bosch et al., 2015) and is currently being explored in multiple clinical trials (Mayer & Arteaga, 2016). The model predicts a set of combinatorial interventions that involve inhibition of PI3K and the MYC-CDK4/6 axis of cell-cycle regulation (e.g. Palbociclib = 1, MYC = 0, cyclinD = 0, CDK4/6 = 0, pRb = 0), which completely block proliferation (Proliferationnorm = 0.00) and maintain the apoptosis-inducing effect of PI3K inhibition (Apoptosisnorm = 0.70). The synergistic effect of PI3K and CDK4/6 inhibition was previously reported (Vora et al., 2014); the rest of the predictions are novel. Mechanistically, these interventions act by blocking the proliferative effect of ER, which makes their anti-proliferative effect as potent as the combination of PI3K and ER inhibitors. A second novel set of combinatorial interventions involve inhibition of PI3K and mTORC1 (e.g. Everolimus = 1, mTORC1 = 0, S6K = 0, and EIF4F = 0), which is predicted to modestly increase the pro-apoptotic effect of PI3K inhibition (Apoptosisnorm = 0.73) and maintain its anti-proliferative effect (Proliferationnorm = 0.25). The synergistic effect of PI3K and mTORC1 inhibition has been documented in MCF7-derived xenografts (Elkabets et al., 2013), but the mechanism has not been identified. Indeed, PI3K inhibition results in mTORC1 downregulation. We find that this combinatorial effect on apoptosis depends on the relative timing of the start of the mTORC1 and PI3K inhibition. Early addition of mTORC1 leads to the inhibition of MCL1, which primes the cells for PI3K-inhibitor-induced apoptosis (see Additional file 2: Table S1), and the maximum apoptosis is the same as when MCL1 is initially set to OFF. Thus, the model predicts that MCL1 inhibition is the mechanism through which PI3K and mTORC1 inhibition are synergistic, and combined PI3K and MCL1 inhibition can mimic the effect of combined PI3K and mTORC1 inhibitors
Network models excel in both aspects of model utility: the integration and interpretation of existing knowledge, and the generation of novel predictions. Our network model (Fig. 3) unites information from numerous studies, reproduces several key experimental and clinical outcomes (Table 4), and visualizes the inter-relationships among various pathways and processes. The overlay of the usually-defined pathways (marked by separate colors) on the signal transduction network that starts with receptor tyrosine kinases and ends with two phenotypic outcomes reveals that these pathways cover a variety of subgraphs of the full network, from linear cascades (such as RAS-MAPK) to bow-tie structured neighborhoods of a node (such as ER signaling) and to subgraphs wherein negative regulation or feedback plays an important role. The inter-regulation among subgraphs is also substantial, and biologically significant. To better illustrate this point, in Fig. 6a we represent each the seven pathways relevant to ER+, PI3K mutant breast cancer as single nodes and indicate the aggregated relationships between them. These relationships summarize one or multiple logically consistent paths between the pathways. For example, mTORC1-induced protein translation, which leads to the increased activity of the anti-apoptotic protein MCL1, yields an overall negative regulation between the mTORC1 pathway (orange rectangle) and the apoptosis pathway (green rectangle). The positive effect of the AKT pathway on the mTORC1 pathway summarizes AKT and PIM’s inhibition of PRAS40, as well as AKT and SGK1’s inhibition of TSC; both PRAS40 and TSC are inhibitors of mTORC1.
AKT inhibits ER signaling through downregulating the histone methyltransferase KMTD and by its inhibition of FOXO3, which would otherwise activate ER. This negative edge stands out from and opposes an otherwise sign-consistent meta-network, wherein the five upstream pathways have positive inter-regulation and all favor proliferation and/or disfavor apoptosis. In ER+, PI3K mutant breast cancer cells, this negative edge dampens (but does not block) ER signaling, and all four other pathways are active; yielding the collective effect of a significant proliferation propensity and lack of apoptosis propensity. In case of drug inhibition of PI3K, four pathways (PI3K, MAPK, AKT, mTORC1) are inactivated, and consequently the break on ER signaling is released (Fig. 6a). The overall effect is a significant apoptosis propensity and a low (but non-zero) proliferation propensity. While the quantification of the two biological outcomes depends on specific model and implementation details, the main message is clearly encapsulated in the network: PI3K inhibition does not eliminate all the proliferation-inducing, apoptosis-resisting activity in the network. Our model provides specific predictions on what additional interventions would yield a significant improvement over PI3K inhibition alone. The targets of these predicted interventions lie in the ER signaling, mTORC1, cell cycle and apoptosis pathways; their names and the nature of their control (inhibition or activation) is also indicated in Fig. 6a. Our finding that multiple combinatorial interventions are effective enables the selection of those that are most effective drug targets and minimize toxicity and side effects.
A novel prediction of the model is that PI3K + CDK4/6 inhibition is a very effective combination treatment because of its ability to both induce cancer cell death and cell cycle arrest by suppressing two parallel proliferation regulator complexes (cyclin E and CDK2, and cyclin D and CDK4/6). So far there are few studies of the combined effect of PI3K inhibitors and CDK4/6 inhibitors (Vora et al., 2014; O’Leary et al., 2016), and the ongoing clinical trials all include ER inhibitors simultaneously with both PI3K and CDK4/6 inhibitors (Mayer & Arteaga, 2016; O’Leary et al., 2016). The model predicts that the combination of PI3K and CDK4/6 inhibitors can be as effective as the combination of PI3K and ER inhibitors, and that the addition of CDK4/6 inhibitors to the latter combination does not further increase its effectiveness. Given that resistance to ER inhibitors can be overcome by CDK4/6 inhibitors (Finn et al., 2009) and that targets of CDK4/6 inhibitors are known resistance mechanisms to ER inhibitors (Hui et al., 2002; Musgrove & Sutherland, 2009), the model predictions suggests that PI3K inhibitors + ER inhibitor followed by PI3K inhibitors + CDK4/6 inhibitors after acquisition of resistance is a better strategy than combined PI3K + ER + CDK4/6 inhibitors.
The network of inter-relationships among pathways can also be used to interpret the existing information and new predictions on potential resistance mechanisms to PI3K inhibition in PI3KCA mutant, ER+ breast cancer (Fig. 6b). Broadly speaking, any mechanism that yields the restoration of activity in the PI3K, MAPK, AKT or mTORC1 pathways, or increased activity of the ER pathway, will restore the proliferation-inducing and/or apoptosis-opposing effects of these pathways, and will yield a decrease in the effectiveness of PI3K inhibitors. For example, a mechanism that would partially restore PI3K activity or PIP3 levels (for example by a loss of function alterations in PTEN (Juric et al., 2015), could lead to the restoration of the PI3K → AKT and PI3K → MAPK edges in Fig. 6b and thus reverse the effects of PI3K inhibition. Constitutive activity of AKT, PIM or SGK1, or inhibition of FOXO3, would at least partially restore the four outgoing edges of the AKT pathway. While one of these edges is to dampen ER signaling, the other three will lead to a decreased apoptosis propensity and increased proliferation propensity. Inspecting the multitude of potential resistance mechanisms (indicated by color-coded node names inside each pathway symbol in Fig. 6b), those in the PI3K, AKT and mTORC1 pathways may be categorized as pathway reactivation, if we consider the union of these three, i.e. PI3K/AKT/mTORC1, as the index pathway. Constitutive ER transcriptional regulatory activity is an example of pathway bypass: it leads to cell cycle progression and activates the anti-apoptotic protein BCL2. Constitutive activity of the MAPK pathway, a model-predicted resistance mechanism, resembles pathway bypass in that it inhibits FOXO3, which would otherwise be accomplished by AKT, but it also overlaps the index pathway through its activation of mTORC1. The model predicts that two different combinations of MAPK and FOXO3 activity (FOXO3 = ON and MAPK = 2 or FOXO3 = OFF and MAPK = 1) can both act as resistance mechanisms to PI3K inhibitors. An analysis of the elements regulated by MAPK and FOXO3 reveals that this happens because two normally opposing effects are allowed to co-occur. In the FOXO3 = ON and MAPK = 2 case, MAPK’s regular inhibition of FOXO3 is blocked, thus this combination yields the proliferative effect of FOXO3 = 1 but a lesser pro-apoptotic effect (due to MAPK = 2). In the FOXO3 = OFF, MAPK = 1 case the apoptosis propensity is decreased because MAPK = 1 inhibits BAD.
Here we focused on PI3KCA mutant breast cancer and targeted PI3K inhibition, which is showing promising results in clinical trials. Our network modeling framework can be used or adapted to answer a broader set of questions. For example, we can consider mutants that have one of the model-identified resistance mechanisms, determine the drugs that overcome the resistance, and identify the most effective combinatorial therapies. An example of such a prediction is that patients with activating genetic alterations in PIM would greatly benefit from the combination of PI3K inhibitors with PIM inhibitors (as one would expect) or the combination of ER and mTOR inhibitors. Although we focused on ER+/HER2- breast cancer, the pathways and mechanisms of the HER2+ subtype are included in the model. Indeed, HER2/HER3 appears as a resistance mechanism to PI3K inhibition and the model recapitulates several resistance mechanisms observed in HER2+ breast cancer (Table 4). The model can be expanded to incorporate additional resistance mechanisms relevant to HER2+ breast cancer (e.g. the FGFR signaling pathway in the context of estrogen receptor degraders (Turner et al., 2010; André & Cortés, 2015; Mao et al., 2017)).
Certain predictions of the model rely on considerations that go beyond the network structure, for example timing. The model predicts a non-monotonic decrease in proliferation in response to PI3K inhibition (Fig. 3). This is due to the convergence of fast signal transduction events that decrease proliferation with the slow ER-driven transcriptional events that increase it. Timing also plays a key role in the predicted synergistic effect on apoptosis induction of mTORC1 inhibition followed by PI3K inhibition. This is because mTORC1 inhibition leads to the inhibition of MCL1, which primes the cells for PI3K-inhibitor-induced apoptosis (see Additional file 2: Table S1). The observation that the timing of drugs can prime cells for apoptosis and yield drug synergy is consistent with previous work showing a similar effect in triple-negative breast cancer (Lee et al., 2012).
Even though the model includes several of the pathways and signaling proteins important in ER+, HER2+, and PIK3CA-mutant breast cancer, it is not complete. The model, for example, does not include the DNA damage pathway, signaling through other RTKs such as FGFR and EGFR, the Wnt pathway, and the TGFβ pathway. Despite the role of these pathways in apoptosis and proliferation of breast cancer cells, we did not include them in the model either because the literature we explored pointed to them behaving very similarly to a pathway included in the model (e.g., in the case of EGFR and FGFR), or because we were not able to find strong evidence linking them to the response/resistance to the targeted drugs studied (e.g., in the case of the DNA damage pathway). We expect that extending the model to account for the effect of other inhibitors (e.g. PARP inhibitors) or other oncogenic processes (e.g. the epithelial-to-mesenchymal transition) would necessitate the inclusion of additional pathways. Note that not including these pathways in the current model is not equivalent to the assumption that they do not play a role in the resistance to the studied targeted therapies; rather, it is a reflection of the expectation that their potential role is mediated through one of the included signaling proteins (or an element of their pathways).
The network model we present in this work, just like any mathematical model, is not final and definitive (Box, 1976; Box, 1979). For several regulatory functions there was insufficient evidence regarding the aggregated effect of multiple regulators; in these cases, we tested several alternatives before settling on the function that most faithfully recapitulated biological results. The model can be improved by experimental elucidation of these regulatory functions and by experimental testing of the model’s predictions. Any discrepancies between the model and experiments would lead us to test changes to the model’s assumptions that resolve the discrepancies while keeping the cases of agreement intact. The improved model would give further predictions that could be tested experimentally, and so on, thus completing the model/experiment cycle inherent to any modeling approach.
The breast cancer network model we present in this work integrates the current knowledge of PIK3CA-mutant, ER+ breast cancers, and uses it to identify a set of elements that may eventually be exploited in high-order therapeutic combinations to achieve a more durable control of breast cancer. The model’s predictions will serve as a basis for guiding and interpreting drug resistance and drug combination studies in ER+ breast cancer. The model can be straightforwardly adapted to HER2+ breast cancer; it already recapitulates multiple outcomes in this setting. The model can be expanded to incorporate multiple additional genetic alterations observed in breast cancer patient cohorts (Koboldt et al., 2012; Pereira et al., 2016; Wagle et al., 2017; Cohen et al., 2017) by appropriately introducing these alterations into the model (e.g., as a node activation or inactivation). The inclusion of the most probable intrinsic or acquired resistance mechanisms to a treatment, informed by pre-, on- and post-treatment genetic characterization of tumors (Cohen et al., 2017), will allow the identification and ranking of the combinatorial interventions that are effective even in the presence of tumor drug resistance. We expect that experimentally and clinically validated network models similar to the one presented here will become an integral part of precision medicine, and will be able to identify successful combinatorial therapies in tumor types and subtypes of interest.
The simulations of the discrete network models were done using the BooleanDynamicModeling Java library, which is freely available on GitHub (https://github.com/jgtz/BooleanDynamicModeling). To simulate multi-level nodes, we use a Boolean variable to denote each level greater than 1. For example, for a 3-level node with states 0, 1, and 2, we have 2 variables (Node and Node_2), and for a 4-level node we have 3 variables (Node, Node_2, and Node_3). The regulatory functions of all the nodes are indicated and explained in Additional File 1. We perform 10,000 simulations in each modeled scenario. The number of time steps are 75 (for the simulations in Figs. Fig. 4 and Fig. 5), 100 (for the simulations in Tables 1 and 3), and up to 120 (for the simulations in Additional file 2: Table S1) depending on the scenario. The code used to simulate the model is available on GitHub (https://github.com/jgtz/BreastCancerModel).
Attractor-finding in discrete network models
To attractors of the model were identified using the StableMotifs Java library, which is freely available on Github (https://github.com/jgtz/StableMotifs) and implements the attractor-finding method based on stable motif analysis (Zañudo & Albert, 2013; Zañudo & Albert, 2015), as has been previously described (Steinway et al., 2014). Stable motif analysis can find the attractors of a logical model by identifying the model’s stable motifs, a set of nodes and their node states with certain identifiable topological (intersecting directed cycles) and dynamical properties (partial steady states), which uniquely determine the attractors of the model (Zañudo & Albert, 2013; Zañudo & Albert, 2015).
The full names of the abbreviated node names in Fig. 3 and thereafter are indicated in Additional File 1.
Here we use the term “inhibitor” to refer to drugs that target and inhibit an element regardless of the specific mechanism of action.
Protein kinase B
Human epidermal growth factor receptor 2
Mitogen-activated protein kinase kinase
Phosphatidylinositol – 4,5- biphosphate 3-kinase
Receptor tyrosine kinase
Aldridge BB, Burke JM, Lauffenburger D a, Sorger PK. Physicochemical modelling of cell signalling pathways. Nat Cell Biol 2006;8(11):1195–1203.
Alessi DR, James SR, Downes CP, Holmes AB, Gaffney PR, Reese CB, et al. Characterization of a 3-phosphoinositide-dependent protein kinase which phosphorylates and activates protein kinase Balpha. Curr Biol. 1997;7(4):261–9.
Alon U. An introduction to systems biology: design principles of biological circuits: CRC press; 2006. https://www.crcpress.com/An-Introduction-to-Systems-Biology-Design-Principles-of-Biological-Circuits/Alon/p/book/9781584886426
Anderson GR, Wardell SE, Cakir M, Crawford L, Leeds JC, Nussbaum DP, et al. PIK3CA mutations enable targeting of a breast tumor dependency through mTOR-mediated MCL-1 translation 2016;175(December):1–15.
André F, Cortés J. Rationale for targeting fibroblast growth factor receptor signaling in breast cancer. Breast Cancer Res Treat. 2015;150:1–8.
Archer TC, Fertig EJ, Gosline SJC, Hafner M, Hughes SK, Joughin BA, et al. Systems Approaches to Cancer Biology. Cancer Res. 2016;76(23):6774 LP–6777.
Bosch A, Li Z, Bergamaschi A, Ellis H, Toska E, Prat A, et al. PI3K inhibition results in enhanced estrogen receptor function and dependence in hormone receptor – positive breast cancer. Sci Transl. 2015;7(283):283ra51.
Box GEP. Science and statistics. J Am Stat Assoc. 1976;71(356):791–9.
Box GEP. Robustness in the strategy of scientific model building. Robustness. Stat. 1979;1:201–36.
Burra PV, Zhang Y, Godzik A, Stec B. Global distribution of conformational states derived from redundant models in the PDB points to non-uniqueness of the protein structure. Proc Natl Acad Sci U S A. 2009;106(26):10505–10.
Carracedo A, Ma L, Teruya-feldstein J, Rojo F, Salmena L. Inhibition of mTORC1 leads to MAPK pathway activation through a PI3K- dependent feedback loop in human cancer. J Clin Invest. 2008;118(9):3065–74.
Castel P, Ellis H, Bago R, Toska E, Razavi P, Carmona FJ, et al. PDK1-SGK1 signaling sustains AKT-independent mTORC1 activation and confers resistance to PI3Kα inhibition. Cancer Cell. 2016;30(2):229–42.
Castellano E, Downward JRAS. Interaction with PI3K: more than just another effector pathway. Genes Cancer. 2011;2(3):261–74.
Chakrabarty A, Sanchez V, Kuba MG, Rinehart C, Arteaga CL. Feedback upregulation of HER3 (ErbB3) expression and activity attenuates antitumor effect of PI3K inhibitors. Proc Natl Acad Sci. 2012;109(8):2718–23.
Chandarlapaty S, Sawai A, Scaltriti M, Rodrik-Outmezguine V, Grbovic-Huezo O, Serra V, et al. AKT inhibition relieves feedback suppression of receptor tyrosine kinase expression and activity. Cancer Cell. 2011;19(1):58–71.
Ciriello G, Gatza ML, Beck AH, Wilkerson MD, Rhie SK, Pastore A, et al. Comprehensive molecular portraits of invasive lobular breast cancer. Cell. 2015;163(2):506–19.
Cohen O, Kim D, Oh C, Waks A, Oliver N, Helvie K, et al. Abstract S1-01: whole exome and transcriptome sequencing of resistant ER+ metastatic breast cancer. Cancer Res 2017 ;77(4 Supplement):S1-1-S1-1.
Collombet S, van Oevelen C, Sardina Ortega JL, Abou-Jaoudé W, Di Stefano B, Thomas-Chollier M, et al. Logical modeling of lymphoid and myeloid cell specification and transdifferentiation. Proc Natl Acad Sci. 2017;114(23):201610622.
Costa C, Ebi H, Martini M, Beausoleil SA, Faber AC, Jakubik CT, et al. Measurement of PIP3 levels reveals an unexpected role for p110β in early adaptive responses to p110α-specific inhibitors in luminal breast cancer. Cancer Cell. 2015;27(1):97–108.
Cree IA, Charlton P. Molecular chess? Hallmarks of anti-cancer drug resistance. BMC Cancer. 2017;17(1):10.
CURRIE RA, WALKER KS, GRAY A, DEAK M, CASAMAYOR A, DOWNES CP, et al. Role of phosphatidylinositol 3,4,5-trisphosphate in regulating the activity and localization of 3-phosphoinositide-dependent protein kinase-1. Biochem J. 1999;337(3):575 LP–583.
Ebi H, Costa C, Faber AC, Nishtala M, Kotani H, Juric D, et al. PI3K regulates MEK/ERK signaling in breast cancer via the Rac-GEF, P-Rex1. Proc Natl Acad Sci U S A. 2013;110(52):21124–9.
Elkabets M, Vora S, Juric D, Morse N, Mino-Kenudson M, Muranen T, et al. mTORC1 inhibition is required for sensitivity to PI3K p110alpha inhibitors in PIK3CA-mutant breast cancer. Sci Transl Med. 2013;5(196):196ra99.
Finn RS, Dering J, Conklin D, Kalous O, Cohen DJ, Desai AJ, et al. PD 0332991, a selective cyclin D kinase 4/6 inhibitor, preferentially inhibits proliferation of luminal estrogen receptor-positive human breast cancer cell lines in vitro. Breast Cancer Res. 2009;11(5):R77.
Friedman AA, Letai A, Fisher DE, Flaherty KT. Precision medicine for cancer with next-generation functional diagnostics. Nat Rev Cancer. 2015;15(12):747–56.
Garg A, Di Cara A, Xenarios I, Mendoza L, De Micheli G. Synchronous versus asynchronous modeling of gene regulatory networks. Bioinformatics. 2008;24(17):1917–25.
Garraway LA, Jänne PA. Circumventing cancer drug resistance in the era of personalized medicine. Cancer Discovery. 2012;2:214–26.
Higgins MJ, Baselga J. Targeted therapies for breast cancer. J Clin Invest. 2011;121(10):3797–803.
Holohan C, Van Schaeybroeck S, Longley DB, Johnston PG. Cancer drug resistance: an evolving paradigm. Nat Rev Cancer. 2013;13(10):714–26.
Hui R, Finney GL, Carroll JS, Lee CSL, Musgrove EA, Sutherland RL. Constitutive overexpression of cyclin D1 but not cyclin E confers acute resistance to antiestrogens in T-47D breast cancer cells. Cancer Res. 2002;62(23):6916–23.
Johannessen CM, Boehm JS. Progress towards precision functional genomics. Curr Opin Syst Biol. 2017;2:74–83.
Juric D, Castel P, Griffith M, Griffith OL, Won HH, Ellis H, et al. Convergent loss of PTEN leads to clinical resistance to a PI(3)Kalpha inhibitor. Nature. 2015;518(7538):240–4.
Kapuy O, Barik D, Domingo Sananes MR, Tyson JJ, Novák B. Bistability by multiple phosphorylation of regulatory proteins. Prog Biophys Mol Biol. 2009;100:47–56.
Klamt S, Saez-Rodriguez J, Lindquist JA, Simeoni L, Gilles EDA. Methodology for the structural and functional analysis of signaling and regulatory networks. BMC Bioinformatics. 2006;7(1):56.
Koboldt DC, Fulton RS, McLellan MD, Schmidt H, Kalicki-Veizer J, McMichael JF, et al. Comprehensive molecular portraits of human breast tumours. Nature. 2012;490(7418):61–70.
Kodack DP, Askoxylakis V, Ferraro GB, Sheng Q, Badeaux M, Goel S, et al. The brain microenvironment mediates resistance in luminal breast cancer to PI3K inhibition through HER3 activation. Sci Transl Med. 2017;9(391):eaal4682.
Kumar Jolly M, Levine H. Computational systems biology of epithelial-hybrid- mesenchymal transitions. Curr Opin. Syst Biol. 2017;3:1–6.
Le X, Antony R, Razavi P, Treacy DJ, Luo F, Ghandi M, et al. Systematic functional characterization of resistance to PI3K inhibition in breast cancer. Cancer Discov. 2016;6(10):1134–47.
Lee MJ, Ye AS, Gardino AK, Heijink AM, Sorger PK, MacBeath G, et al. Sequential application of anticancer drugs enhances cell death by rewiring apoptotic signaling networks. Cell. 2012;149(4):780–94.
Lim KH, Counter CM. Reduction in the requirement of oncogenic Ras signaling to activation of PI3K/AKT pathway during tumor maintenance. Cancer Cell. 2005;8(5):381–92.
Mao P, Quartey Q, Cohen O, Piccioni F, Wagle N. Abstract P3-03-08: a large-scale functional screen to identify resistance mechanisms to selective estrogen receptor degraders fulvestrant and GDC-810 in ER+ breast cancer. Cancer Res 2017 ;77(4 Supplement):P3-3-8-P3-03–8.
Massarweh S, Osborne CK, Creighton CJ, Qin L, Tsimelzon A, Huang S, et al. Tamoxifen resistance in breast tumors is driven by growth factor receptor signaling with repression of classic estrogen receptor genomic function. Cancer Res. 2008;68(3):826–33.
Mayer IA, Arteaga CL. The PI3K/AKT pathway as a target for cancer treatment. Annu Rev Med. 2016;67(1):11–28.
Méndez-López LF, Davila-Velderrain J, Domínguez-Hüttinger E, Enríquez-Olguín C, Martínez-García JC, Alvarez-Buylla ER. Gene regulatory network underlying the immortalization of epithelial cells. BMC Syst Biol. 2017;11(1):24.
Meric-Bernstam F, Mills GB. Overcoming implementation challenges of personalized cancer therapy. Nat Rev Clin Oncol. 2012;9(9):542–8.
Miller TW, Balko JM, Fox EM, Ghazoui Z, Dunbier A, Anderson H, et al. ERα-dependent E2F transcription can mediate resistance to estrogen deprivation in human breast cancer. Cancer Discov. 2011;1(4):338–51.
Miller TW, Hennessy BT, González-Angulo AM, Fox EM, Mills GB, Chen H, et al. Hyperactivation of phosphatidylinositol-3 kinase promotes escape from hormone dependence in estrogen receptor-positive human breast cancer. J Clin Invest. 2010;120(7):2406–13.
Milo R, Jorgensen P, Moran U, Weber G, Springer M. BioNumbers—the database of key numbers in molecular and cell biology. Nucleic Acids Res. 2010;38(suppl_1):D750–3.
Milo R, Phillips R. Cell biology by the numbers: Garland Science; 2015. http://garlandscience.com/product/isbn/9780815345374
Morris MK, Saez-Rodriguez J, Sorger PK, Lauffenburger DA. Logic-based models for the analysis of cell signaling networks. Biochemistry. 2010;49:3216–24.
Muellner MK, Uras IZ, Gapp BV, Kerzendorfer C, Smida M, Lechtermann H, et al. A chemical-genetic screen reveals a mechanism of resistance to PI3K inhibitors in cancer. Nat Chem Biol. 2011;7(11):787–93.
Musgrove EA, Sutherland RL. Biological determinants of endocrine resistance in breast cancer. Nat Rev Cancer. 2009;9(9):631–43.
Nahta R, Yuan LXH, Zhang B, Kobayashi R, Esteva FJ. Insulin-like growth factor-I receptor/human epidermal growth factor receptor 2 heterodimerization contributes to trastuzumab resistance of breast cancer cells. Cancer Res. 2005;65(23):11118–28.
O’Leary B, Finn RS, Turner NC. Treating cancer with selective CDK4/6 inhibitors. Nat Rev Clin Oncol. 2016;13(7):417–30.
O’Reilly KE, Rojo F, She QB, Solit D, Mills GB, Smith D, et al. mTOR inhibition induces upstream receptor tyrosine kinase signaling and activates Akt. Cancer Res. 2006;66(3):1500–8.
Pereira B, Chin S-F, Rueda OM, H-KM V, Provenzano E, Bardwell HA, et al. The somatic mutation profiles of 2,433 breast cancers refines their genomic and transcriptomic landscapes. Nat Commun. 2016;7:11479.
Rodrik-Outmezguine VS, Chandarlapaty S, Pagano NC, Poulikakos PI, Scaltriti M, Moskatel E, et al. mTOR kinase inhibition causes feedback-dependent biphasic regulation of AKT signaling. Cancer Discov. 2011;1(3):248–59.
Saadatpour A, Albert I, Albert R. Attractor analysis of asynchronous Boolean models of signal transduction networks. J Theor Biol. 2010;266(4):641–56.
Saadatpour A, Albert R, Reluga TC. A reduction method for Boolean network models proven to conserve attractors. SIAM J Appl Dyn Syst. 2013;12(4):1997–2011.
Sarbassov DD, Guertin DA, Ali SM, Sabatini DM. Phosphorylation and Regulation of Akt/PKB by the Rictor-mTOR Complex. Science. 2005;307(5712):1098 LP–1101.
Sarosiek KA, Fraser C, Muthalagu N, Bhola PD, Chang W, McBrayer SK, et al. Developmental regulation of mitochondrial apoptosis by c-Myc governs age- and tissue-specific sensitivity to cancer therapeutics. Cancer Cell. 2017;31(1):142–56.
Serra V, Scaltriti M, Prudkin L, Eichhorn PJ, Ibrahim YH, Chandarlapaty S, et al. PI3K inhibition results in enhanced HER signaling and acquired ERK dependency in HER2-overexpressing breast cancer. Oncogene. 2011;30(22):2547–57.
Steinway SN, Zanudo JGT, Ding W, Rountree CB, Feith DJ, Loughran TP, et al. Network modeling of TGFβ signaling in hepatocellular carcinoma epithelial-to-mesenchymal transition reveals joint sonic hedgehog and Wnt pathway activation. Cancer Res. 2014;74(21):5963–77.
Stephens PJ, Tarpey PS, Davies H, Van Loo P, Greenman C, Wedge DC, et al. The landscape of cancer genes and mutational processes in breast cancer. Nature. 2012;486(7403):400–4.
Tian X, Huang B, Zhang X-P, Lu M, Liu F, Onuchic JN, et al. Modeling the response of a tumor-suppressive network to mitogenic and oncogenic signals. Proc Natl Acad Sci. 2017;114(21):5337–42.
Toska E, Osmanbeyoglu HU, Castel P, Chan C, Hendrickson RC, Elkabets M, et al. PI3K pathway regulates ER-dependent transcription in breast cancer through the epigenetic regulator KMT2D. 2017;355(6331):1324–30.
Turke AB, Song Y, Costa C, Cook R, Arteaga CL, Asara JM, et al. MEK inhibition leads to PI3K/AKT activation by relieving a negative feedback on ERBB receptors. Cancer Res. 2012;72(13):3228–37.
Turner N, Pearson A, Sharpe R, Lambros M, Geyer F, Lopez-Garcia MA, et al. FGFR1 amplification drives endocrine therapy resistance and is a therapeutic target in breast cancer. Cancer Res. 2010;70(5):2085–94.
Tyson JJ, Baumann WT, Chen C, Verdugo A, Tavassoly I, Wang Y, et al. Dynamic modelling of oestrogen signalling and cell fate in breast cancer cells. Nat Rev Cancer. 2011;11(7):523–32.
Udyavar AR, Wooten DJ, Hoeksema M, Bansal M, Califano A, Estrada L, et al. Novel hybrid phenotype revealed in small cell lung cancer by a transcription factor network model that can explain tumor heterogeneity. Cancer Res. 2017;77(5):1063–74.
Vasudevan KM, Barbie DA, Davies MA, Rabinovsky R, McNear CJ, Kim JJ, et al. AKT-independent signaling downstream of oncogenic PIK3CA mutations in human cancer. Cancer Cell. 2009;16(1):21–32.
Veliz-Cuba A, Aguilar B, Hinkelmann F, Laubenbacher R. Steady state analysis of Boolean molecular network models via model reduction and computational algebra. BMC Bioinformatics. 2014;15(1):221.
Vora S, Juric D, Kim N, Mino-Kenudson M, Huynh T, Costa C, et al. CDK 4/6 inhibitors sensitize PIK3CA mutant breast cancer to PI3K inhibitors. Cancer Cell. 2014;26(1):136–49.
Wagle N, Painter C, Anastasio E, Dunphy M, McGillicuddy M, Kim D, et al. The Metastatic Breast Cancer (MBC) project: Accelerating translational research through direct patient engagement. J Clin Oncol. 2017;35(15_suppl):1076.
Wang R-S, Saadatpour A, Albert R. Boolean modeling in systems biology: an overview of methodology and applications. Phys Biol. 2012;9(5):55001.
Werner HMJ, Mills GB, Ram PT. Cancer systems biology: a peek into the future of patient care? Nat Rev Clin Oncol. 2014;11(3):167–76.
Will M, Qin ACR, Toy W, Yao Z, Rodrik-Outmezguine V, Schneider C, et al. Rapid induction of apoptosis by PI3K inhibitors is dependent upon their transient inhibition of RAS-ERK signaling. Cancer Discov. 2014;4(3):334–48.
Zañudo JGT, Albert R. An effective network reduction approach to find the dynamical repertoire of discrete dynamic networks. Chaos. 2013;23(2).
Zañudo JGT, Albert R. Cell fate reprogramming by control of intracellular network dynamics. PLoS Comput Biol. 2015;11(4):e1004193.
Zehir A, Benayed R, Shah RH, Syed A, Middha S, Kim HR, et al. Mutational landscape of metastatic cancer revealed from prospective clinical sequencing of 10,000 patients. Nat Med. 2017;23(6):703–13.
Zhang J, Tian X-J, Zhang H, Teng Y, Li R, Bai F, et al. TGF-β-induced epithelial-to-mesenchymal transition proceeds through stepwise activation of multiple feedback loops. Sci Signal. 2014;7(345):ra91.
Zhang S, Huang W-C, Li P, Guo H, Poh S-B, Brady SW, et al. Combating trastuzumab resistance by targeting SRC, a common node downstream of multiple resistance pathways. Nat Med. 2011;17(4):461–9.
Zhang Y, Kwok-Shing Ng P, Kucherlapati M, Chen F, Liu Y, Tsang YH, et al. A Pan-Cancer Proteogenomic Atlas of PI3K/AKT/mTOR Pathway Alterations. Cancer Cell. 2017;31(6):820–32. e3
Zwang Y, Jonas O, Chen C, Rinne ML, Doench JG, Piccioni F, et al. Synergistic interactions with PI3K inhibition that induce apoptosis. elife. 2017;6:e24523.
J.G.T.Z. would like to thank Levi Garraway and Nikhil Wagle for hosting him as a visiting Research Fellow at the Dana-Farber Cancer Institute and Broad Institute of Harvard and MIT. We would also like to thank Pingping Mao and Anthony Letai for their helpful suggestions on the manuscript. This project is part of a group effort by the Stand Up to Cancer Drug Combination Convergence Team; we thank the members of the team for constructive discussions and feedback.
Availability of data and materials
All the data are contained in the text and additional files.
This work was funded by the National Science foundation grant NSF PHY 1545832, the Stand Up to Cancer Foundation (to R.A. and M.S.), a Stand Up to Cancer Foundation/The V Foundation Convergence Scholar Award (D2015–039) to J.G.T.Z, and National Institutes of Health grant P30CA008748 and the Breast Cancer Research Foundation (to M.S.).
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file 1:
Details of the model: full node names, regulatory functions and update information. (DOCX 38 kb)
Additional file 2: Table S1.
Timing-dependent combinatorial effect of PI3K inhibitors and mTOR inhibitors. (DOCX 11 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Gómez Tejeda Zañudo, J., Scaltriti, M. & Albert, R. A network modeling approach to elucidate drug resistance mechanisms and predict combinatorial drug treatments in breast cancer. Cancer Converg 1, 5 (2017). https://doi.org/10.1186/s41236-017-0007-6
- Breast cancer
- Signal transduction networks
- Network model
- Dynamic model
- Combination therapy