A network modeling approach to elucidate drug resistance mechanisms and predict combinatorial drug treatments in breast cancer

Abstract

Background

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.

Results

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.

Conclusions

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.

Background

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.

Fig. 1

A logical dynamic network model of signal transduction. a Simplified network of PI3K and MAPK signaling initiated by receptor tyrosine kinases (RTKs). b Three possible trajectories (using general asynchronous updating) of the model constructed from the network shown in panel A and the regulatory functions in Eq. (1), with an initial condition in which the only active node is Growth Factor (GF). c Time course of the activity (average node state) of each node using equal update probability for all nodes (left) or using a smaller update probability for the Transcription Factors (TF) node compared to the rest of the nodes (right). Inset shows a zoom in of the time course for the early time steps. Note that the time courses of the activity of several nodes overlap in panels B and C, in particular, RAS and PI3K, RAF and PIP3, and MEK/ERK and AKT

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 f_i 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 f_i, 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

$$ {\displaystyle \begin{array}{l}{\mathrm{f}}_{\mathrm{Growth}\ \mathrm{Factor}}={\upsigma}_{\mathrm{Growth}\kern0.34em \mathrm{Factor}}\\ {}{\mathrm{f}}_{\mathrm{RTK}}={\upsigma}_{\mathrm{Growth}\kern0.34em \mathrm{Factor}}\\ {}{\mathrm{f}}_{\mathrm{RAS}}={\upsigma}_{\mathrm{RTK}}\\ {}{\mathrm{f}}_{\mathrm{PI}3\mathrm{K}}={\upsigma}_{\mathrm{RTK}}\;\mathrm{OR}\;{\upsigma}_{\mathrm{RAS}}\\ {}{\mathrm{f}}_{\mathrm{BRAF}}={\upsigma}_{\mathrm{RAS}}\\ {}{\mathrm{f}}_{\mathrm{PI}\mathrm{P}3}={\upsigma}_{\mathrm{PI}\mathrm{P}3\mathrm{K}}\\ {}{\mathrm{f}}_{\mathrm{MEK}/\mathrm{ERK}}={\upsigma}_{\mathrm{BRAF}}\\ {}{\mathrm{f}}_{\mathrm{AKT}}={\upsigma}_{\mathrm{PI}\mathrm{P}3}\\ {}{\mathrm{f}}_{\mathrm{Transcripton}\kern0.34em \mathrm{Factor}\mathrm{s}}={\upsigma}_{\mathrm{AKT}}\;\mathrm{AND}\;{\upsigma}_{\mathrm{MEK}/\mathrm{ERK},}\end{array}} $$

(1)

which are mathematical statements of the transmission of information between the elements in RTK signaling. For example, f_RTK = σ_{Growth Factor} indicates that the RTK becomes active in the presence of external growth factors, and f_PI3K = σ_RTK OR σ_RAS indicates that PI3K becomes active in the presence of either active RTK or active RAS. Another example is f_{Transcription 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 f_j (σ_j(t) = f_j[Σ(t − 1)], where Σ(t) = (σ₁(t), σ₂(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 f_{Growth 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 (a_i(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⁻³-1 s for signaling events and ~10¹–10² 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 f_RTK = σ_{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 f_RAS = 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.

Fig. 2

Drug inhibition and resistance mechanisms in dynamic network models. a The network model of Fig. 1a with additional nodes denoting a RTK inhibitor (RTKi) and a MEK inhibitor (MEKi). b-d Time courses of node activity (average node state) in response to Growth Factor (GF) in the presence of RTKi, a RAS activating mutation, and MEKi. For panel B, we start with an initial state in which the only active node is Growth Factor and there is no RTK inhibitor (σ_RTKi = OFF); we introduce the RTKi by setting σ_RTKi = ON for time ≥ 20. For panel C, we start with an initial state with σ_GF = σ_RTKi = ON and introduce a RAS activating mutation by setting f_RAS = ON for time ≥ 20. For panel D, we start with an initial state with σ_GF = σ_RTKi = ON and with the modified function f_RAS = ON, and introduce a MEK inhibitor by setting σ_MEKi = ON for time ≥ 20. Note that the time courses of the activity of several nodes overlap in panels c and d, in particular, PI3K and RAF, and PIP3 and MEK/ERK

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.

Results

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.

Fig. 3

Network model of oncogenic signal transduction in ER+ breast cancer. The nodes are colored according to the pathway they are part of: RTK signaling, PI3K signaling, MAPK pathway, AKT pathway, mTORC1 pathway, ER signaling, cell-death signaling (apoptosis) and cell-cycle regulation (proliferation). The network also includes selected drugs of interest in the context of breast cancer: Alpelisib (PI3K inhibitor), Ipatasertib (AKT inhibitor), Fulvestrant (ER inhibitor), Palbociclib (CDK4/6 inhibitor), Everolimus (mTOR inhibitor), Trametinib (MEK inhibitor), and Neratinib (HER1/2 inhibitor). For clarity of the figure, we merge certain nodes into a single node when there is no ambiguity; for example, a node denoting a transcript is merged with the protein it codes, e.g., BCL2 and BCL2_T (for BCL2 transcript) are shown as BCL2. The full list of network nodes is indicated in Additional File 1

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 interest¹: 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.

$$ {\mathrm{f}}_{\mathrm{AKT}}=\left(\mathrm{PIP}3\ \mathrm{or}\ \mathrm{PIP}3\_2\right)\ \mathrm{and}\ \left(\mathrm{PDK}1\_\mathrm{pm}\ \mathrm{or}\ \mathrm{mTORC}2\_\mathrm{pm}\right)\ \mathrm{and}\ \left(\mathrm{not}\ \mathrm{Ipatasertib}\ \mathrm{or}\ \mathrm{PIP}3\_2\right) $$

$$ {\mathrm{f}}_{\mathrm{ER}\_\mathrm{transcription}}=\mathrm{ER}\ \mathrm{and}\ \left(\mathrm{ESR}1\ \mathrm{or}\ \mathrm{ESR}1\_2\right) $$

$$ {\mathrm{f}}_{\mathrm{ER}\_\mathrm{transcription}\_2}=\mathrm{KMT}2\mathrm{D}\ \mathrm{and}\ \mathrm{FOXA}1\ \mathrm{and}\ \mathrm{PBX}1\ \mathrm{and}\ \mathrm{ESR}1\_2\ \mathrm{and}\ \mathrm{ER} $$

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, f_AKT 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 f_{ER_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 f_{ER_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 Apoptosis_norm and Proliferation_norm (Additional File 1), which in the current simulations take the initial values Apoptosis_norm = 0.00 and Proliferation_norm = 0.50, and the final values Apoptosis_norm = 0.70 and Proliferation_norm = 0.25.

Fig. 4

Network model response to PI3K inhibitors. a Time course of node activity (based on 10,000 simulations) in response to PI3K inhibition from time = 2. Only the nodes that change during the time course are shown. For multi-state nodes, we show the node activity for each node state and denote each state of a multi-state node with a “_n”, where n is the state it is referring to (e.g. MAPK_1 refers to state 1 of MAPK and MYC_2 refers to state 2 of MYC). The Apoptosis_norm and Proliferation_norm is a weighted and normalized (between 0 and 1) measure of the state of the node Apoptosis and Proliferation, respectively. b-c Time course of selected nodes, each representative of a different pathway. Panel C zooms in to the early time points of Panel B

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 Apoptosis_norm = 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 Proliferation_norm = 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, Apoptosis_norm and Proliferation_norm values despite PI3K inhibition, in the case of the 90% effective PIM there is a decrease in the average level of mTORC1 and Proliferation_norm and an increase in Apoptosis_norm, closer to the result obtained for SGK1.

Fig. 5

Illustration of PIM and SGK1-mediated resistance to PI3K inhibitors. a The relevant subnetwork of the full network shown in Fig. 3. PIM shares four targets of AKT (compare blue and red edges), while SGK1 shares two (green edges). Their post-translational regulation is different from AKT’s: SGK1 is activated by different pools of PDK1 and mTORC2 than AKT, while PIM is constitutively active. b Time course of node activity in response to PI3K inhibition at time = 2 in cells with full activity of PIM (top) or reduced PIM activity (10% chance of inactive PIM at any time step) (bottom). c Time course of node activity in response to PI3K inhibition at time = 2 in cells with constitutive PDK1 and SGK1 activity. The symbol legend applies to both B and C

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 Apoptosis_norm and Proliferation_norm) 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).

Table 1 Single-node resistance mechanisms to PI3K inhibitors ordered by their effect on Apoptosis

Table 2 Double-node resistance mechanisms to PI3K inhibitors ordered by their effect on Apoptosis

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 (Proliferation_norm = 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 (Proliferation_norm = 0.50) and also lower apoptosis significantly (Apoptosis_norm = 0.13), and FOXO3 = 1 together with MAPK = 2, which blocks apoptosis and restores proliferation (Proliferation_norm = 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 (Proliferation_norm = 0.00 − 0.13 and Apoptosis_norm = 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 (Proliferation_norm = 0.00) and maintain the apoptosis-inducing effect of PI3K inhibition (Apoptosis_norm = 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 (Apoptosis_norm = 0.73) and maintain its anti-proliferative effect (Proliferation_norm = 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

Table 3 Single-node interventions which in combination with PI3K inhibitors yield an increase in Apoptosis_norm from 0.7 and/or a decrease in Proliferation_norm from 0.25

Discussion

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.

Table 4 Illustration of experimental and clinical outcomes in ER+ and HER2+ breast cancer reproduced by the model

Fig. 6

Meta-network illustrating synergistic interventions and resistance mechanisms to PI3K inhibitors. The colored rectangles correspond to the pathways introduced in Fig. 3, and the edges between them represent aggregated regulatory relationships between pathways. In these relationships, and also when referring to a pathway as active or inactive, we focus on the index (named) member of the pathway. Thus, when saying that the mTORC1 pathway is active we mean that mTORC1, EIF4F and S6K are active, TSC and PRAS40 are inactive. Thick continuous lines indicate active pathways/interactions, thick and dashed lines represent partially active pathways/interactions, thin and dashed lines mean inactive pathways/interactions. For the node names indicated inside the colored rectangles, blue indicates inhibition/inactivity and red indicates increased activity. a Signaling pathway activity in response to PI3K inhibition. ER signaling is still active (partly due to the release of its inhibition by AKT), while the apoptosis and proliferation pathways are partially active. Inhibition of the nodes indicated in blue font or constitutive activity of Rb is predicted to have a synergistic effect with PI3K inhibition. b Resistance mechanisms to PI3K inhibitors. Sustained activity of the nodes indicated with red font inside each pathway can (at least partially) restore the pathway’s activation and obstruct the effectiveness of PI3K inhibition. Sustained inactivity of the nodes indicated with blue font can have a similar effect. For simplicity, the HER2/HER3 resistance mechanism is not included in a separate RTK module but as part of the pathways activated by HER2/HER3, namely MAPK and PI3K

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.

Conclusions

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.

Methods

Model simulations

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.