Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Boolean modeling of breast cancer signaling pathways uncovers mechanisms of drug synergy

  • Kittisak Taoma,

    Roles Formal analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing

    Affiliations Bioinformatics and Systems Biology Program, School of Bioresources and Technology, King Mongkut’s University of Technology Thonburi, Bangkok, Thailand, School of Information Technology, King Mongkut’s University of Technology Thonburi, Bangkok, Thailand

  • Marasri Ruengjitchatchawalya,

    Roles Conceptualization, Formal analysis, Methodology, Writing – review & editing

    Affiliations Bioinformatics and Systems Biology Program, School of Bioresources and Technology, King Mongkut’s University of Technology Thonburi, Bangkok, Thailand, Biotechnology Program, School of Bioresources and Technology, King Mongkut’s University of Technology Thonburi, Bangkok, Thailand

  • Monrudee Liangruksa ,

    Roles Conceptualization, Formal analysis, Funding acquisition, Methodology, Writing – review & editing

    monrudee@nanotec.or.th (ML); teeraphan.lao@kmutt.ac.th (TL)

    Affiliation National Nanotechnology Center, National Science and Technology Development Agency (NSTDA), Pathum Thani, Thailand

  • Teeraphan Laomettachit

    Roles Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Supervision, Writing – original draft, Writing – review & editing

    monrudee@nanotec.or.th (ML); teeraphan.lao@kmutt.ac.th (TL)

    Affiliations Bioinformatics and Systems Biology Program, School of Bioresources and Technology, King Mongkut’s University of Technology Thonburi, Bangkok, Thailand, Theoretical and Computational Physics Group, Center of Excellence in Theoretical and Computational Science, King Mongkut’s University of Technology Thonburi, Bangkok, Thailand

Abstract

Breast cancer is one of the most common types of cancer in females. While drug combinations have shown potential in breast cancer treatments, identifying new effective drug pairs is challenging due to the vast number of possible combinations among available compounds. Efforts have been made to accelerate the process with in silico predictions. Here, we developed a Boolean model of signaling pathways in breast cancer. The model was tailored to represent five breast cancer cell lines by integrating information about cell-line specific mutations, gene expression, and drug treatments. The models reproduced cell-line specific protein activities and drug-response behaviors in agreement with experimental data. Next, we proposed a calculation of protein synergy scores (PSSs), determining the effect of drug combinations on individual proteins’ activities. The PSSs of selected proteins were used to investigate the synergistic effects of 150 drug combinations across five cancer cell lines. The comparison of the highest single agent (HSA) synergy scores between experiments and model predictions from the MDA-MB-231 cell line achieved the highest Pearson’s correlation coefficient of 0.58 with a great balance among the classification metrics (AUC = 0.74, sensitivity = 0.63, and specificity = 0.64). Finally, we clustered drug pairs into groups based on the selected PSSs to gain further insights into the mechanisms underlying the observed synergistic effects of drug pairs. Clustering analysis allowed us to identify distinct patterns in the protein activities that correspond to five different modes of synergy: 1) synergistic activation of FADD and BID (extrinsic apoptosis pathway), 2) synergistic inhibition of BCL2 (intrinsic apoptosis pathway), 3) synergistic inhibition of MTORC1, 4) synergistic inhibition of ESR1, and 5) synergistic inhibition of CYCLIN D. Our approach offers a mechanistic understanding of the efficacy of drug combinations and provides direction for selecting potential drug pairs worthy of further laboratory investigation.

Introduction

Breast cancer ranks among the most prevalent forms of cancer, with the number of new cases estimated to become approximately 3.19 million worldwide in 2040 [1]. Various breast cancer treatments are currently available, such as targeted-drug therapy, hormone therapy, and chemotherapy [2]. However, the development of drug resistance hinders the effectiveness of these treatments. For example, targeted drug treatments can trigger a rapidly adaptive resistance mechanism of cancer cells by rewiring feedback loops in the protein signaling pathways [3]. In addition, chemotherapy can introduce genomic instability, leading to new resistance clones of cancer cells [4, 5].

Promisingly, drug combination therapy has been shown to mitigate drug resistance by simultaneously targeting more than one signaling pathway, effectively diminishing proliferation or inducing apoptosis [610]. Additionally, combined drugs with diverse mechanisms of action allow for lower individual drug doses than when used as single agents, resulting in fewer side effects while improving patient survival [11, 12].

When two drugs are combined, their effect may interact to enhance or reduce their effectiveness compared to when used separately, referred to as synergistic and antagonistic effects, respectively [13]. Several models have been proposed to quantify the degree of synergy, such as Highest Single Agent (HSA), Loewe, ZIP, and Bliss [14]. These models share a common approach to calculating the synergy score S by determining the deviation of drug combination effect yc from the expected effect ye, S = ycye. Each model uses a different hypothesis for the expected effect calculation. For example, the HSA model assumes that the expected effect (ye) of drug combination is equal to the effect of either drug A (yA) or drug B (yB), whichever is the greater. Then, the HSA synergy score (SHSA) is calculated as SHSA = yc − max(yA, yB), where yc is the effect of the drug combination. Positive and negative values of SHSA indicate synergistic and antagonistic effects, respectively.

Despite the availability of high-throughput assays for drug combination screening [1518], discovering effective drug combinations remains challenging. As the number of potential targeted drugs grows, ample combinatorial space of available drugs is largely unexplored. Therefore, there is a need for computational tools to facilitate the rational selection of potential drug pairs for further investigation in the laboratory.

Many machine-learning algorithms have been proposed for drug synergy prediction in pan-cancer cell lines [1923]. While the prediction performances were satisfactory, much was still to be learned about the underlying mechanism of the synergy between the predicted drug pairs. The ’black-box’ nature of the machine learning models made it difficult to fully understand how the predictions were made. Complementary to machine learning, mechanistic modeling is an alternative approach to understanding the mechanism of cancer’s drug resistance and identifying potential drug combination intervention. Boolean modeling is one of the simplest mechanistic approaches, where kinetic parameters are not required. Thus, creating a larger protein signaling network with Boolean modeling is more feasible compared to kinetic-based models. As a result, the Boolean model can include more relevant proteins, such as drug target proteins [24]. Multiple Boolean modeling frameworks have been proposed to investigate drug resistance mechanisms [2527] and drug combination effects [2831] in various types of cancer. However, these models did not quantify the degree of synergy, which can be quantified experimentally as a synergy score on a continuous scale (e.g., higher scores mean higher synergy). Providing such a metric is essential to drug combination therapy research [32].

In this study, we modeled signaling pathways in breast cancer with a Boolean model, which was tailored to represent five specific breast cancer cell lines by incorporating genetic mutation and gene expression profiles. The cell-line specific models were calibrated to reproduce gene expression profiles and drug-treatment experiments. Then, the models were used to simulate a dose-response profile of 150 drug pairs retrieved from the DrugComb database [33], and the activities of representative proteins selected by a genetic algorithm were used to calculate the synergy score. Finally, the activities of the selected proteins were used to gain insight into the mechanisms of drug combinations with both synergistic and antagonistic interactions.

Methods

Our approach shares similarities with the methodology employed by Tsirvouli et al. (2020) [28]. The study predicted drug synergy in various cancer cell lines, including one breast cancer cell line (MDA-MB-468). The present work focused on drug synergistic effects on breast cancer, extending the analysis to five breast cancer cell lines (MCF-7, T-47D, BT-549, MDA-MB-231, and MDA-MB-468). These five cell lines are commonly used in breast cancer research to represent molecular subtypes of the ER-positive (MCF-7 and T-47D) and triple-negative (BT-549, MDA-MB-231, and MDA-MB-468) cancer (for example, see [3437]). The five cell lines are also experimentally reported in the DrugComb database with the highest numbers of tested drug pairs among other breast cancer cell lines [33].

Fig 1 summarizes the workflow in the present study with five steps. Briefly, the first step involved reconstructing a protein signaling network that comprehensively represents breast cancer proliferation pathways by integrating information from various databases, including KEGG, SIGNOR, and SignaLink (Step 1). Next, the network was translated to a Boolean model (Step 2). The model was further tailored into five cell-line specific Boolean models by integrating gene expression and genetic mutation profiles of each cell line. Then, each cell-line specific model was simulated to obtain steady-state protein activities under unperturbed and drug-perturbed conditions. The protein activities of each cell line under unperturbed conditions were compared to their respective gene expression data. The results of the drug perturbation simulations were compared to the observed responses of the cell lines to drug treatments, which were collected from the literature. The Boolean rules of each cell-line specific model were manually modified during this step to match model simulations with gene expression data and observed behaviors from drug-perturbed experiments (Step 3). Then, the resulting Boolean models were used to simulate the effects of drug combinations on protein activities. We proposed a calculation of protein synergy scores (PSSs), determining the effect of drug combinations on individual proteins’ activities. The PSSs of selected proteins were summed together and compared to the synergistic scores from the experiments (Step 4). The PSS profiles were also used to investigate the mechanisms by which drug pairs exert their synergistic or antagonistic effects on the cell lines (Step 5). Below, we describe each step in detail.

1. Protein signaling network reconstruction

Breast cancer-related protein signaling pathways were retrieved from the KEGG database (https://www.genome.jp/pathway/hsa05224), which contains 76 proteins and 60 interactions (accessed on 1 March 2020). Another set of 41 proteins and an additional 176 interactions from the SIGNOR2.0 (https://signor.uniroma2.it/) and 12 interactions from SignaLink (http://signalink.org/) databases (accessed on 1 January 2022) were combined. A network was used to represent the signaling pathways, where nodes represent the proteins and edges represent the interactions. Proteins forming a complex and proteins with multiple isoforms were merged into a single entity. For example, the CYCLIN_D_c node represents the CDK4-CCND1 complex, and the AKT_i node accounts for AKT1, AKT2, and AKT3. In addition, self-interactions (auto-activation) of nine ligand nodes (e.g., estrogen; ES and progesterone; PG) were included by assuming that these ligands are produced by cancer cells in an autocrine manner [38]. The resulting protein signaling network comprises 117 nodes and 257 edges (Fig 2). The complete node and interaction list from KEGG, SIGNOR2.0, and SignaLink can be found in S1 File.

thumbnail
Fig 2. A generic protein signaling network of breast cancer pathways.

The network was reconstructed using information from KEGG, SIGNOR, and SignaLink databases. Nodes represent protein entities and edges represent the interactions between them. Proteins forming a complex were merged into a single node (names ending with ’_c’, e.g., CYCLIN_D_c represents the CDK4-CCND1 complex). Proteins with multiple isoforms were also merged into a single node (names ending with ’_i’, e.g., AKT_i represents AKT1, AKT2, and AKT3). The black and red arrows represent activating and inhibitory regulation, respectively. The complete node list can be found in S1 File.

https://doi.org/10.1371/journal.pone.0298788.g002

2. Generic Boolean model

We implemented the Boolean modeling approach to represent the network, where each node was associated with a value of either 0 or 1 to represent the protein activity/level (0 = off/low and 1 = on/high). The Boolean values were updated over discrete time steps according to logical rules derived from protein interactions as defined in the databases. The logical rules were initially assigned as follows.

  1. An or function was used to describe the updating rules of nodes with multiple activating regulators, e.g., C(t+1) = A(t) or B(t), where C(t+1) is the Boolean value of node C in the next time step, and A(t) and B(t) are the Boolean values of the activator nodes A and B in the current time step.
  2. A not or function was used to describe the updating rules of nodes with multiple inhibitory regulators, e.g., C(t+1) = not (A(t) or B(t)). Here, node C turns on in the next time step only if nodes A and B are off in the current time step.
  3. An inhibitory regulator wins over all activating regulators, e.g., D(t+1) = (A(t) or B(t)) and not(C(t)), where nodes A and B are activating regulators while node C is an inhibitory regulator of node D.

These assumptions were initially made for the generic model in a manner similar to the previous studies [2830]. However, more specific rules were derived when biological support was present. For instance, AKT is recruited by PIP3 at the plasma membrane. Subsequently, AKT is phosphorylated twice by MTORC2 and PDPK1 at the S473 and T308 sites, respectively, to be fully activated [39]. These interactions were described as

The complete Boolean rules of the generic model (Fig 2) can be found in S1 File.

3. Cell-line specific models

The generic Boolean model describes common signaling pathways related to cell growth and proliferation in breast cancer. However, the same combination of drugs can produce different synergistic effects on different cell lines. To account for these heterogenic responses, the generic Boolean model was tailored into cell-line specific models to capture each cancer cell line’s unique characteristics and responses to drug perturbation. We created five cell-line specific models for MCF-7, T-47D, BT-549, MDA-MB-231, and MDA-MB-468 by

  1. integrating cell line mutation profiles into the Boolean functions,
  2. setting initial states of the cell lines based on gene expression profiles and
  3. modifying the Boolean rules to match model simulations with gene expression data and observed behaviors from drug-perturbed experiments. These modifications are explained in detail below.

3.1 Integrating cell line mutation profiles into the Boolean functions.

The missense mutations for the five cell lines were identified from the Cell Model Passports (https://cellmodelpassports.sanger.ac.uk/) and DepMap (https://depmap.org/portal/) databases. Subsequently, the mutations were annotated for the biological function using the OncoKB database (https://www.oncokb.org/). We assigned a fixed Boolean value of 0 or 1 to nodes whose mutations were annotated as loss-of-function or gain-of-function, respectively. S1 Table displays the mutation profile of seven genes across the five breast cancer cell lines utilized in the present study.

3.2 Setting initial states of the cell lines based on gene expression profiles.

The gene expression data normalized with the transcripts per million (TPM) method was retrieved from Cell Model Passports. The normalized gene expression values were further scaled across five cancer cell lines by the Min-Max normalization method (Eq 1) to obtain values between 0 and 1: (1) where and xjk are the scaled and original expression values of gene j from cell line k, respectively. max(xj) and min(xj) are the maximum and minimum, respectively, of the expression values of gene j among five breast cancer cell lines. The scaled values (S1 Fig) were then used as a probability that the corresponding nodes were on at the beginning of each simulation. For example, if the scaled value of gene J of cell line K is 0.7, the initial state of protein J will be set to 1 (on) with a 70% chance at the beginning of the simulation of cell line K.

3.3 Modifying the Boolean rules to match model simulations with gene expression data and observed behaviors from drug-perturbed experiments.

Next, each cell-line specific model (i.e., the generic model that incorporated genetic mutation and gene expression profiles) was simulated to generate protein activities at the steady state. The simulation results were compared to gene expression data and observed behaviors from drug-perturbed experiments of each respective cell line. During this step, we manually modified the Boolean rules of each cell-line specific model to capture observed behaviors derived from the experimental data.

To compare protein activities simulated from the Boolean model (values of 0 or 1) to gene expression data (values obtained from the Cell Model Passports database and scaled to range between 0 and 1), we performed model simulation for 5000 time steps using a uniform asynchronous update scheme. In this scheme, only one node from all nodes that could change their state was selected with equal probability among others to change per time step. One hundred repeats were simulated for each cell-line specific model. The resulting Boolean values (0 or 1) of each protein from all repeats were averaged over the same time steps to obtain the protein activities (values ranging between 0 and 1), and the protein activities at the steady state were compared to the retrieved gene expression data from the Cell Model Passports database. The time step at which the steady state of the protein activities was reached was determined by entropy H(t) [40]: (2) where Pt(S) is the probability the state S is observed over time step 0 to t. The entropy was calculated across the simulation time step (t = 0 to 5000) to obtain the entropy profile. A higher entropy value indicates a higher impurity level in the system (i.e., various states of protein activities are observed over time t). When a particular state of protein activities becomes dominant over time, indicating the steady state, the entropy level approaches 0 as the probability of observing the steady state is near 1, i.e., . However, it may take a long computational time to simulate the model until the system is completely homogenous with a steady state. Therefore, we assume that the protein activities have reached a steady state when the slope of the entropy profile is at its minimum. The protein activities from the time step the steady state was reached to the final simulation time step (t = 5000) were averaged to represent the steady-state protein activities.

Next, we collected 42 drug-treatment experiments on the five cell lines from published literature. We conducted a corresponding simulation for each experiment by fixing the Boolean values of the drug target nodes to 0 or 1, depending on whether the drug was an antagonist or agonist. We then compared the simulated response of the model to the observed behavior of the cell lines in the experiments.

During this step, inconsistency between simulations and experimental observations was addressed by manually modifying logical rules (e.g., changing the rule from and to or, and vice versa) and by adding/removing edges until the Boolean model of each cell line reproduced the gene expression and drug-treatment experiment observations. The adjustments were informed based on the literature as listed in S1 File.

For example, the generic Boolean function of MTORC1 (3) was modified into the ER-positive specific Boolean function (4) because, in the generic Boolean model, MTORC1 was activated only by RHEB, which is a downstream target of AKT (Eq 3). However, the regulatory function of MTORC1 of T-47D and MCF-7 further included NRG1, IGF1, and PIM (Eq 4). This modification was based on supporting evidence showing that the presence of NRG1, IGF1, or PIM could maintain the MTORC1 activity even when an inhibitor suppressed the AKT activity [9, 41].

In another instance, KMT2D, which is regulated positively and negatively by PAX7 and AKT, respectively, was described by a generic Boolean function (5)

The function was modified to (6) in MCF-7 to reduce the contribution of AKT in inactivating KMT2D. The modification was made to match the simulation result with the observed expression data of KMT2D in MCF-7, which has a scaled value of 0.78 (out of the maximum value of 1.0). In contrast, for MDA-MB-231, the generic function in Eq 5 was used as the scaled expression data of KMT2D in this cell line was 0.

The complete logical rules for each cell-line specific model are listed in S1 File. The 42 experiments whose information was used to modify and validate the Boolean rules will be discussed further in the Results section (Table 2 below).

4. Drug combination simulations

We selected 150 drug combinations from the DrugComb database (https://drugcomb.fimm.fi/) to investigate using our cell-line specific models (S1 File). To simulate the drug combination effects, we identified the target proteins of 33 unique drugs in the dataset from Drugbank (https://go.drugbank.com/) and Therapeutic Target Database (TTD; https://db.idrblab.net/ttd/) (Table 1). Each cell-line specific model was simulated without drug perturbation until the steady state of protein activities was reached (Eq 2). Then, the effects of the drug combinations were applied. Four doses (0, 0.25, 0.75, and 1) of each drug in combination were investigated for each drug pair. The doses represented the probability of the target-protein nodes being off or on (depending on whether the drug is an antagonist or agonist). For example, if fulvestrant (an ESR1 inhibitor) was used with a dose of 0.25, the ESR1 node had a probability of 0.25 to be off determined at the time the drug was applied to the model in each simulation repeat. Thirty simulation repeats were performed for each combination dose. The simulations continued until the steady state of protein activities was reached again.

thumbnail
Table 1. Drug target proteins identified from Drugbank and Therapeutic target database.

https://doi.org/10.1371/journal.pone.0298788.t001

Implementing drug effects on protein targets as a probability between 0 and 1 reflects the concentration ranges of monotherapy that cover the minimum and the maximum effect of drugs as implemented in the drug combination experiments [4244]. We assume that a drug concentration with the maximum effect (equivalent to a dose of 1.0 in our simulation) would fully inhibit the protein targets with a 100% chance. Similarly, a drug concentration at half the maximum effect (equivalent to a dose of 0.5 in our simulation) would inhibit the protein targets with a 50% chance. However, as simulating various dose-response profiles can be computationally intensive, we have selected to use only four doses for each drug. As a result, the simulations produced 4×4 protein activity profiles per drug pair per cell line.

To calculate drug synergy scores from our model and compare them to the values reported in the DrugComb database, we proposed the calculation of the PSS. First, we categorized proteins in the model into onco- or tumor-suppressor proteins. Then, we calculated the PSS using the protein activities at the steady state, according to Eqs 7 and 8: (7) (8) where XSS and YSS represent the steady-state protein activity of an oncoprotein and tumor-suppressor protein, respectively. The subscripts 1, 2, and 1+2 indicate the steady-state protein activity after being perturbed by drug 1, drug 2, and both, respectively. The calculation of PSSs for each dose followed the methodology for computing the highest single agent (HSA) synergy score. The value of the PSS is between −1 and +1, where the positive and negative values were interpreted as the drug pair having a synergistic and antagonistic effect on the protein level, respectively. The PSS of 0 indicated the additive effect. The PSS at the 75th percentile across all dose combinations was designated the consensus PSS for that protein. Finally, the PSSs of selected proteins were summed together to represent the predicted HSA synergy score for each drug pair.

To determine which proteins whose PSSs should be used to calculate the predicted HSA score, a genetic algorithm was implemented to find the best representatives from both onco- and tumor-suppressor proteins by the following steps:

  1. The initial population in the genetic algorithm comprised 60 chromosomes with a random length (the number of representative proteins) and randomly chosen proteins, except E2F1 and CASP3, which were selected by default.
  2. The predicted HSA scores of the 150 drug pairs were calculated from each chromosome by summing the PSSs of all proteins in the chromosome. Then, the chromosomes were evaluated by their resulting Pearson’s correlation coefficient between the predicted and observed HSA synergy scores. The top 30 chromosomes that yielded the highest Pearson’s correlation coefficients were selected.
  3. Genetic mutation and crossing-over were applied both with the probability of 0.9 to the selected chromosomes to create 30 offspring chromosomes. The high probability values (0.9) were assigned due to the large search space.
  4. New 30 chromosomes were created and combined with the offspring chromosomes to constitute the new population.
  5. The processes in steps 2–4 were repeated for 10,000 generations to ensure the optimally selected protein representatives.

5. Mechanism explanations

Finally, the potential synergistic and antagonistic mechanisms of drug combinations were investigated by clustering drug pairs into groups based on the similarity in the selected PSSs. Through an analysis of the PSS profile patterns, we identified the contributions of individual protein activities to the observed synergy of drug pairs, characterizing them as distinct modes of drug synergism.

Results

1. Protein signaling network reconstruction

Using the STRINGdb API (https://pypi.org/project/stringdb/), the majority of the original 76 signaling proteins from the KEGG database were significantly enriched in the PI3K-AKT (hsa04151: N = 51 FDR = 8.74×10−64) and MTOR (hsa04150: N = 29, FDR = 7.05×10−50) signaling pathways.

When 41 more nodes from the SIGNOR and SignaLink databases were included in the network, more pathways showed significant enrichment, e.g., apoptosis (hsa04210: N = 21 from KEGG and 8 from SIGNOR and SignaLink, FDR = 3.46×10−42), ERBB (hsa04012: N = 22 from KEGG and 8 from SIGNOR and SignaLink, FDR = 4.04×10−48), MAPK (hsa04010: N = 22 from KEGG and 9 from SIGNOR and SignaLink, FDR = 3.34×10−33), and FOXO (hsa04068: N = 29 from KEGG and 2 from SIGNOR and SignaLink, FDR = 7.49×10−46) signaling pathways.

Overall, the reconstructed protein signaling network comprehensively covered multiple breast cancer growth and proliferation pathways, including the PI3K-AKT, MTOR, apoptosis, ERBB, MAPK, and FOXO signaling pathways.

2. Model validation

The five cell-line specific models were simulated with and without drug perturbations to obtain protein activities at the steady state. The unperturbed protein activities from each cell line were compared to their respective cell lines’ actual gene expression data retrieved from the Cell Model Passports database. The drug perturbation simulations were compared to observed responses of the cell lines to drug treatments collected from 42 experiments. The results of these simulations demonstrated that the models could capture the unique characteristics of each cell line and their responses to drug perturbations.

2.1 Model simulation compared to gene expression profiles.

We conducted simulations for each cell-line specific Boolean model without drug perturbations to derive the steady-state protein activities. The simulated steady-state protein activities exhibited a reasonable level of correlation with the actual gene expression data, with Pearson’s correlation coefficients ranging between 0.67 and 0.85. MDA-MB-231 and T-47D were the cell lines from triple-negative and ER-positive subtypes with the highest correlation at 0.85 and 0.75, respectively. An identical Pearson’s correlation at 0.67 was observed for MDA-MB-468 and BT-549, whereas MCF-7 was achieved at 0.68. When we performed clustering analysis, combining these simulated steady-state protein activities with the actual gene expression data, the results accurately segregated into the ER-positive (T-47D and MCF-7) and triple-negative (BT-549, MDA-MB-231, and MDA-MB-468) sub-groups, as demonstrated in Fig 3. This suggests that our cell-line specific models behaved similarly to their unperturbed cell-line behaviors.

thumbnail
Fig 3. The simulated steady-state protein activities of each cell line.

The steady-state protein activities simulated from the five cell-line specific Boolean models were clustered with the actual gene expression profiles (labeled as ’exp’ in the figure) based on the correlation distance. The actual gene expression profiles were scaled with the Min-Max normalization method to yield a value between 0 and 1 before clustering.

https://doi.org/10.1371/journal.pone.0298788.g003

2.2 Model simulation compared to drug-treatment experiments.

The Boolean model of each cell line was further validated by drug-treatment experiments from the literature. Table 2 shows that the models captured many drug-resistance behaviors, mainly involving the release of negative regulation of an oncogenic pathway. We used the activities of E2F1 and CASP3 in representing cell proliferation and apoptosis in response to single-drug and combination treatments.

thumbnail
Table 2. Drug-treatment simulations from cell-line specific models.

https://doi.org/10.1371/journal.pone.0298788.t002

Some examples of simulations from Table 2 are illustrated in Fig 4. Fig 4A demonstrates the release of negative regulation on the MAPK (Ras-Raf-MEK-ERK) pathway from MTORC1 inhibition (rapamycin) in BT-549, leading to rapamycin resistance [8]. In Fig 4B, MTORC1 inhibition released the inhibition on IRS1, triggering the activation of AKT in MCF-7 [45]. In Fig 4C, a MEK inhibitor removed the inhibition on EGFR by ERK, leading to the activation of the PI3K/AKT pathway in MDA-MB-231 [47]. Fig 4D plots the simulation of a resistance mechanism in MDA-MB-468 in response to a PI3K inhibitor. The inhibition of PI3K released the negative regulation of the IRS1 protein by the AKT-MTORC1-S6K1 pathway, which subsequently activated the JAK2-STAT5 cascade, allowing the cells to escape cell death [10].

thumbnail
Fig 4. Simulations of resistance mechanisms in response to drug treatments.

(A) MTORC1 inhibition in BT-549. (B) MTORC1 inhibition in MCF-7. (C) MEK inhibition in MDA-MB-231. (D) PI3K inhibition in MDA-MB-468.

https://doi.org/10.1371/journal.pone.0298788.g004

3. Drug synergy predictions

Next, we used the validated cell-line specific models to calculate the PSSs for 150 drug pairs retrieved from the DrugComb database. The PSSs of 13 nodes (BCL2, CYCLIN D, MTORC1, ESR1, PGR, PAK1, STAT3, WNT1, BID, GATA3, FADD, P21, and CASP3) were selected by the genetic algorithm to calculate the predicted HSA synergy score (see Methods). The selected 13 nodes contribute to different breast cancer-related pathways: cell cycle (CYCLIN D and P21), apoptosis (BCL2, BID, CASP3, and FADD), estrogen signaling (ESR1), EGFR signaling (PAK1), progesterone signaling (PGR), Wnt signaling (WNT1), JAK/STAT signaling (GATA3 and STAT3), and PI3K/AKT signaling (MTORC1) pathways. The comparison between the predicted HSA synergy scores and the experimental HSA scores from MDA-MB-231 achieved the highest Pearson’s correlation coefficient of 0.58 with a great balance among the classification metrics (AUC = 0.74, sensitivity = 0.63, and specificity = 0.64). When considering all five cell-line specific models, the models achieved a Pearson correlation coefficient of 0.326 (Fig 5), an AUC of 0.62 with a sensitivity of 0.50, and a specificity of 0.67 (Fig 6).

thumbnail
Fig 5. Scatter plot between actual and predicted HSA synergy scores.

The actual and predicted HSA scores were scaled with the Min-Max normalization method to achieve a score between −1 and +1.

https://doi.org/10.1371/journal.pone.0298788.g005

thumbnail
Fig 6. Prediction performance of the models.

Pearson’s correlation coefficients, AUC, sensitivity, and specificity scores of the cell-line specific models are shown.

https://doi.org/10.1371/journal.pone.0298788.g006

4. Predicted mechanisms

Finally, the potential synergistic and antagonistic mechanisms of drug combinations were analyzed by clustering the drug pairs into groups based on the 13 selected PSSs. For the synergistic group, 42 drug pairs, correctly predicted as true positives, from five cell lines were categorized into five groups (Fig 7). For groups 1 and 2, the synergistic effect mainly relied on the induction of the apoptosis pathway (indicated by a positive PSS of CASP3). In group 1, the extrinsic apoptosis pathway was synergistically upregulated (positive PSSs of FADD and BID), while the intrinsic pathway was synergistically induced in group 2 (a positive PSS of BCL2). An example of a drug pair from group 1 was AZD5363 + MK-2206. MK-2206 (AKT3 inhibitor) alone can induce apoptosis in T-47D [52].Combining MK-2206 with AZD5363 (AKT1 and AKT3 inhibitor) resulted in synergistic activation of apoptosis mainly from the extrinsic pathway. However, the combination also activated ESR1 and CYCLIN D (the negative values of PSSs in Fig 7) due to the release of negative regulation on FOXO3 by AKT, which is consistent with experimental observations [7].

thumbnail
Fig 7. Clustering analysis of true-positive synergistic drug pairs based on the protein synergy scores of the genetic algorithm-selected proteins.

The positive scores (green) indicate increased or decreased protein activities due to the drug combination perturbation compared to the single-drug treatment for tumor-suppressor or oncoproteins, respectively. The negative scores (red) indicate vice versa.

https://doi.org/10.1371/journal.pone.0298788.g007

An example of a drug combination in group 2 included MK-2206 (AKT inhibitor) + selumetinib (MEK inhibitor), strongly inducing the intrinsic apoptosis pathway in MCF-7. It is possible that selumetinib increased the sensitivity of cells to MK-2206 in inducing intrinsic apoptosis because MK-2206 alone could not effectively cause apoptosis in MCF-7 when compared to unperturbed cells [53]. However, the negative effect on CYCLIN D was also observed due to the release of the inhibition of FOXO3 by AKT [7]. In MDA-MB-468, the combination of erlotinib hydrochloride (EGFR inhibitor) and imatinib (ABL1 inhibitor) in group 2 increased the intrinsic apoptosis and reduced the G1/S transition and protein translation (indicated by the positive PSSs of CYCLIN D and MTORC1, respectively). This predicted mechanism aligns with a previous report involving the combination of a different EGFR inhibitor (lapatinib) and imatinib (an ABL1 inhibitor), where it was observed that the level of MYC, a downstream target of MTORC1, exhibited a synergistic reduction in MDA-MB-468 cells [54].

For groups 3 to 5, the proliferation was mainly inhibited (indicated by the positive PSSs of CYCLIN D and/or MTORC1). An example of a drug pair in group 3 was MK-2206 + AZD5363 (both target the AKT proteins) treated in BT-549. Unlike the effect of the same drug combination in the ER-positive T-47D cell line, the drug pair synergistically reduced the CYCLIN D activity in BT-549. Drug pairs in group 4 were the combination of an estrogen receptor inhibitor (tamoxifen citrate or fulvestrant) and another drug. The combinations reduced the G1/S transition by synergistically inhibiting ESR1 and activating P21. The synergy of drug pairs in group 5 was contributed mainly by the down-regulation of CYCLIN D. An example was an EGFR inhibitor (gefitinib or lapatinib) in combination with antibiotic AY 22989 (MTORC1 inhibitor). The combination of EGFR and MTOR inhibitors was previously reported with a synergistic effect in MDA-MB-231 [49]. We predict the mechanism to be the synergistic reduction in G1/S transition and protein translation (due to the observed positive PSSs of CYCLIN D and MTORC1).

Table 3 lists seven drug pairs whose synergistic predictions were supported by in vitro assays from the literature. Here, we predict the synergistic contributions made by individual proteins and list them in Table 3.

thumbnail
Table 3. Synergistic drug pairs with the literature support.

https://doi.org/10.1371/journal.pone.0298788.t003

For the antagonistic group, Fig 8 shows three separate groups of 45 drug pairs, correctly predicted as true negatives. In group 1, many drug pairs contained megestrol acetate, which has an agonist effect on the hormone receptor PGR. Thus, drugs combined with megestrol acetate had a highly negative value of PSS in PGR (e.g., vandatinib + megestrol acetate in MCF-7 and ruxolitinib + megestrol acetate in MDA-MB-231).

thumbnail
Fig 8. Clustering analysis of true-negative antagonistic drug pairs based on the protein synergy scores of the genetic algorithm-selected proteins.

The positive scores (green) indicate increased or decreased protein activities due to the drug combination perturbation compared to the single-drug treatment for tumor-suppressor or oncoproteins, respectively. The negative scores (red) indicate vice versa.

https://doi.org/10.1371/journal.pone.0298788.g008

In group 2, most drug pairs contained an ESR1 agonist (emcyt, mitotane, or raloxifene). Therefore, drug pairs in this group exhibited a highly negative PSS value of ESR1 (e.g., gefitinib + emcyt in MDA-MB-468 and raloxifene + imatinib in MDA-MB-231 or nilotinib in BT-549). Finally, we observed in group 3 that most drug pairs failed to inhibit proliferation nor promote apoptosis (indicated by the negative PSSs of MTORC1 and CASP3). For example, trisenox, which has an agonistic effect on JUN, MAPK3, MAPK1, and AKT1, combined with either erlotinib hydrochloride or megestrol acetate in MDA-MB-468 exhibited the highly negative PSSs of MTORC1 and CASP3.

Discussion and conclusions

Finding drug pairs for effective breast cancer treatment is challenging due to a large combination space that requires resources and time to investigate. In the present study, we utilized a Boolean modeling approach to capture the drug-response behaviors of cancer cells and predict the synergistic effects between drug pairs. Our approach shares similarities with the methodology employed by Tsirvouli et al. (2020) [28] but extends the analysis to five breast cancer cell lines covering both ER-positive and triple-negative subtypes of breast cancers. In addition, a unique feature in the current study is the calculation of the protein synergy scores (PSSs) that determine the effect of drug combinations on individual proteins’ activities. The calculation allowed the synergistic effects of drug pairs to be decomposed into the contributions made by individual proteins, thereby enabling the interpretation of the mechanism of the synergy through the profile of protein activities.

Our results identified five patterns of synergy among drug pairs: 1) synergistic activation of FADD and BID (extrinsic apoptosis pathway), 2) synergistic inhibition of BCL2 (intrinsic apoptosis pathway), 3) synergistic inhibition of MTORC1, 4) synergistic inhibition of ESR1, and 5) synergistic inhibition of CYCLIN D. Interestingly, the combination of AZD5363 and MK-2206 exhibited distinct synergistic mechanisms in T-47D and BT-549. In T-47D, it synergistically activated the extrinsic apoptosis pathway, while in BT-549, it prompted a synergistic inhibition of the proliferation pathway. Furthermore, our models have indicated that antagonistic interactions among drug pairs predominantly arise from the agonistic properties of specific drugs, which activate proteins such as ESR1, PGR, AKT, and ERK. Consequently, it is important to account for these factors during the development and evaluation of drug combinations to ensure their efficacy.

Our current study focused on simulating combinations of only two drugs because most drug-perturbed experiments involved only drug pairs. In addition, the HSA synergy scores from the DrugComb database to which we compared our synergy predictions were also limited to two drugs. However, cancer patients can be treated with more than two drugs if they offer greater efficacy [12], help to overcome drug resistance [2], and avoid drug toxicity [61]. It is worth noting that our simulation framework allows the simulation of more than two drug combinations, thereby providing further valuable predictions, which are planned for future investigation of the study.

We acknowledge that there is room for significant improvements in the future development of our models. First, the models require the integration of additional drug resistance mechanisms associated with ESR1 and EGFR. These proteins were frequently present as drug-pair targets, resulting in false-positive predictions (S2 Fig). By incorporating the resistance mechanisms, cancer cells may enhance their ability to counteract the effects of drugs, shifting the false positive to true negative predictions. For the false negative predictions, we speculate that including more experimental observations for each specific drug could help rectify this issue. At present, the validation of our models relied on the effects of only a limited number of drugs per cell line (as indicated in Table 2), which differs from all simulated drugs used in synergy predictions (as listed in Table 1).

Second, we utilized the highest single agent (HSA) model among the available synergy scores to quantify the drug pair’s synergistic effect. HSA has been widely used for drug synergy prediction in many types of cancer due to the simplest assumption [2830]. It should be noted, however, that the synergy model selection depends on individual preference [62], and the conflict among them can be observed (see examples in S2 Table). To reduce the misclassification from a single synergy score metric, a method capable of prioritizing synergistic drug pairs based on the consensus of all synergy score metrics is needed.

Third, integrating phosphoproteomic data would significantly enhance the model’s performance. The genomic and transcriptomic data have been used for transforming a generic Boolean model into a cell-line specific one due to its abundant availability [63, 64]. However, protein signaling networks are predominately influenced by the enzymatic activity of signaling proteins. Thus, the phosphoproteomic data should be used more precisely in the model development.

Additionally, the use of asynchronous updates in the current study did not consider the varying time scales of different reactions (e.g., protein synthesis versus phosphorylation/dephosphorylation). Incorporating this aspect into future versions of the model may enhance the simulation’s accuracy. Furthermore, it should be mentioned that the five cell lines currently selected for the study are not representative of all subtypes of breast cancer cells. Among the five cell lines, three are triple-negative breast cancer (TNBC), and two are estrogen receptor-positive (ER-positive). This limited selection of cell lines only covers a small portion of breast cancer cell heterogeneity. Therefore, it is essential to investigate more cell lines using the proposed framework of this study to ensure the generalizability of the findings to a broader range of cell types.

Despite its limitations, our approach offers a mechanistic understanding of the efficacy of drug combinations. Additionally, it provides valuable direction for selecting promising drug pairs worthy of further laboratory investigation.

Supporting information

S1 File. Boolean rules of the generic and cell-line specific models.

https://doi.org/10.1371/journal.pone.0298788.s001

(XLSX)

S1 Table. The missense mutations of five cell lines from cell model passports and DepMap.

https://doi.org/10.1371/journal.pone.0298788.s002

(DOCX)

S2 Table. Examples of drug pairs with inconsistent synergy scores across five cell lines.

https://doi.org/10.1371/journal.pone.0298788.s003

(DOCX)

S1 Fig. Clustering analysis of gene expression profiles.

https://doi.org/10.1371/journal.pone.0298788.s004

(DOCX)

S2 Fig. The frequency of drug target proteins found in false-positive predictions.

https://doi.org/10.1371/journal.pone.0298788.s005

(DOCX)

Acknowledgments

The authors acknowledge computational resources from the NSTDA Supercomputer Center (ThaiSC). The authors appreciate the editor and reviewers’ valuable comments and suggestions, which have greatly improved the manuscript.

References

  1. 1. Ferlay J, Laversanne M, Ervik M, Lam F, Colombet M, Mery L, et al. Global cancer observatory: cancer tomorrow [Internet]. Lyon, France: International Agency for Research on Cancer. 2020 [cited 2023 Jun 1]. https://gco.iarc.fr/tomorrow
  2. 2. Fisusi FA, Akala EO. Drug combinations in breast cancer therapy. Pharm Nanotechnol. 2019;7(1):3–23. pmid:30666921
  3. 3. Chandarlapaty S. Negative feedback and adaptive resistance to the targeted therapy of cancer. Cancer Discov. 2012;2(4):311–9. pmid:22576208
  4. 4. Murugaesu N, Wilson GA, Birkbak NJ, Watkins TBK, McGranahan N, Kumar S, et al. Tracking the genomic evolution of esophageal adenocarcinoma through neoadjuvant chemotherapy. Cancer Discov. 2015;5(8):821–32. pmid:26003801
  5. 5. Glen CD, Dubrova YE. Exposure to anticancer drugs can result in transgenerational genomic instability in mice. Proc Natl Acad Sci U S A. 2012;109(8):2984–8. pmid:22308437
  6. 6. Jia J, Zhu F, Ma X, Cao ZW, Li YX, Chen YZ. Mechanisms of drug combinations: interaction and network perspectives. Nat Rev Drug Discov. 2009;8(2):111–28. pmid:19180105
  7. 7. 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 Med. 2015;7(283):283ra51. pmid:25877889
  8. 8. Carracedo A, Ma L, Teruya-Feldstein J, Rojo F, Salmena L, Alimonti A, et al. 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. pmid:18725988
  9. 9. 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. pmid:27604488
  10. 10. Britschgi A, Andraos R, Brinkhaus H, Klebba I, Romanet V, Müller U, et al. JAK2/STAT5 inhibition circumvents resistance to PI3K/mTOR blockade: a rationale for cotargeting these pathways in metastatic breast cancer. Cancer Cell. 2012;22(6):796–811. pmid:23238015
  11. 11. Mauri D, Polyzos NP, Salanti G, Pavlidis N, Ioannidis JPA. Multiple-treatments meta-analysis of chemotherapy and targeted therapies in advanced breast cancer. J Natl Cancer Inst. 2008;100(24):1780–91. pmid:19066278
  12. 12. Xu L, Wu X, Hu C, Zhang Z, Zhang L, Liang S, et al. A meta-analysis of combination therapy versus single-agent therapy in anthracycline- and taxane-pretreated metastatic breast cancer: results from nine randomized Phase III trials. Onco Targets Ther. 2016;9:4061–74. pmid:27445497
  13. 13. Foucquier J, Guedj M. Analysis of drug combinations: current methodological landscape. Pharmacol Res Perspect. 2015;3(3):e00149. pmid:26171228
  14. 14. Malyutina A, Majumder MM, Wang W, Pessia A, Heckman CA, Tang J. Drug combination sensitivity scoring facilitates the discovery of synergistic and efficacious drug combinations in cancer. PLoS Comput Biol. 2019;15(5):e1006752. pmid:31107860
  15. 15. Flobak Å, Niederdorfer B, Nakstad VT, Thommesen L, Klinkenberg G, Lægreid A. A high-throughput drug combination screen of targeted small molecule inhibitors in cancer cell lines. Sci Data. 2019;6(237). pmid:31664030
  16. 16. Gu Z, Shi C, Li J, Han Y, Sun B, Zhang W, et al. Palbociclib-based high-throughput combination drug screening identifies synergistic therapeutic options in HPV-negative head and neck squamous cell carcinoma. BMC Med. 2022;20(175). pmid:35546399
  17. 17. Perez DR, Edwards BS, Sklar LA, Chigaev A. High-throughput flow cytometry drug combination discovery with novel synergy analysis software, SynScreen. SLAS Discovery. 2018;23(7):751–60. pmid:29842834
  18. 18. He L, Kulesskiy E, Saarela J, Turunen L, Wennerberg K, Aittokallio T, et al. Methods for high-throughput drug combination screening and synergy scoring. Methods Mol Biol. 2018;1711:351–98. pmid:29344898
  19. 19. Chen G, Tsoi A, Xu H, Zheng WJ. Predict effective drug combination by deep belief network and ontology fingerprints. J Biomed Inform. 2018;85:149–54. pmid:30081101
  20. 20. Sidorov P, Naulaerts S, Ariey-Bonnet J, Pasquier E, Ballester PJ. Predicting synergism of cancer drug combinations using NCI-ALMANAC data. Front Chem. 2019;7:1–13.
  21. 21. Li X, Xu Y, Cui H, Huang T, Wang D, Lian B, et al. Prediction of synergistic anti-cancer drug combinations based on drug target network and drug induced gene expression profiles. Artif Intell Med. 2017;83:35–43. pmid:28583437
  22. 22. Jeon M, Kim S, Park S, Lee H, Kang J. In silico drug combination discovery for personalized cancer therapy. BMC Syst Biol. 2018;12:16. pmid:29560824
  23. 23. Preuer K, Lewis RPI, Hochreiter S, Bender A, Bulusu KC, Klambauer G. DeepSynergy: predicting anti-cancer drug synergy with deep learning. Bioinformatics. 2018;34(9):1538–46. pmid:29253077
  24. 24. Karlebach G, Shamir R. Modelling and analysis of gene regulatory networks. Nat Rev Mol Cell Biol. 2008;9(10):770–80. pmid:18797474
  25. 25. der Heyde S V., Bender C, Henjes F, Sonntag J, Korf U, Beißbarth T. Boolean ErbB network reconstructions and perturbation simulations reveal individual drug response in different breast cancer cell lines. BMC Syst Biol. 2014;8(1):75.
  26. 26. 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. 2017;1(1):5. pmid:29623959
  27. 27. Zañudo JGT, Mao P, Alcon C, Kowalski K, Johnson GN, Xu G, et al. Cell line-specific network models of ER+ breast cancer identify potential PI3KCA inhibitor resistance mechanisms and drug combinations. Cancer Res. 2021;81(17):4603–17.
  28. 28. Tsirvouli E, Touré V, Niederdorfer B, Vázquez M, Flobak Å, Kuiper M. A middle-out modeling strategy to extend a colon cancer logical model improves drug synergy predictions in epithelial-derived cancer cell lines. Front Mol Biosci. 2020;7:1–19.
  29. 29. Niederdorfer B, Touré V, Vazquez M, Thommesen L, Kuiper M, Lægreid A, et al. Strategies to enhance logic modeling-based cell line-specific drug synergy prediction. Front Physiol. 2020;11:1–16.
  30. 30. Flobak Å, Baudot A, Remy E, Thommesen L, Thieffry D, Kuiper M, et al. Discovery of drug synergies in gastric cancer cells predicted by logical modeling. PLoS Comput Biol. 2015;11(8):1–20.
  31. 31. Vitali F, Cohen L, Demartini A, Amato A, Eterno V, Zambelli A, et al. A network-based data integration approach to support drug repurposing and multi-target therapies in triple negative breast cancer. PLoS One. 2016;11(9):e0162407. pmid:27632168
  32. 32. Ma J, Motsinger-Reif A. Current methods for quantifying drug synergism. Proteom Bioinform. 2019;1(2):43–8. pmid:32043089
  33. 33. Zheng S, Aldahdooh J, Shadbahr T, Wang Y, Aldahdooh D, Bao J, et al. DrugComb update: a more comprehensive drug sensitivity data repository and analysis portal. Nucleic Acids Res. 2021;49:W174–84. pmid:34060634
  34. 34. Lim D, Cho JG, Yun E, Lee A, Ryu HY, Lee YJ, et al. MicroRNA 34a-AXL axis regulates vasculogenic mimicry formation in breast cancer cells. Genes (Basel). 2020;12(1):9. pmid:33374832
  35. 35. Sarmiento-Salinas FL, Delgado-Magallón A, Montes-Alvarado JB, Ramírez-Ramírez D, Flores-Alonso JC, Cortés-Hernández P, et al. Breast cancer subtypes present a differential production of reactive oxygen species (ROS) and susceptibility to antioxidant treatment. Front Oncol. 2019;9:480. pmid:31231612
  36. 36. Morotti M, Bridges E, Valli A, Choudhry H, Sheldon H, Wigfield S, et al. Hypoxia-induced switch in SNAT2/SLC38A2 regulation generates endocrine resistance in breast cancer. Proc Natl Acad Sci U S A. 2019;116(25):12452–61. pmid:31152137
  37. 37. Shi Y, Yu Y, Wang Z, Wang H, Bieerkehazhi S, Zhao Y, et al. Second-generation proteasome inhibitor carfilzomib enhances doxorubicin-induced cytotoxicity and apoptosis in breast cancer cells. Oncotarget. 2016;7(45):73697–710. pmid:27655642
  38. 38. Sporn MB, Roberts AB. Autocrine growth factors and cancer. Nature. 1985;313(6005):745–7. pmid:3883191
  39. 39. Manning BD, Toker A. AKT/PKB signaling: Navigating the network. Cell. 2017;169(3):381–405. pmid:28431241
  40. 40. Stoll G, Viara E, Barillot E, Calzone L. Continuous time boolean modeling for biological signaling: application of Gillespie algorithm. BMC Syst Biol. 2012;6:1–18.
  41. 41. Elkabets M, Vora S, Juric D, Morse N, Mino-Kenudson M, Muranen T, et al. mTORC1 inhibition is required for sensitivity to PI3K p110α inhibitors in PIK3CA-mutant breast cancer. Sci Transl Med. 2013;5(196):196ra99.
  42. 42. Houweling M, Giczewska A, Abdul K, Nieuwenhuis N, Küçükosmanoglu A, Pastuszak K, et al. Screening of predicted synergistic multi-target therapies in glioblastoma identifies new treatment strategies. Neurooncol Adv. 2023;5(1):vdad073. pmid:37455945
  43. 43. Jaaks P, Coker EA, Vis DJ, Edwards O, Carpenter EF, Leto SM, et al. Effective drug combinations in breast, colon and pancreatic cancer cells. Nature. 2022;603(7899):166–73. pmid:35197630
  44. 44. Mäkelä P, Zhang SM, Rudd SG. Drug synergy scoring using minimal dose response matrices. BMC Res Notes. 2021;14(1):1–6.
  45. 45. 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. pmid:16452206
  46. 46. 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. Science. 2017;355(6331):1324–30. pmid:28336670
  47. 47. Mirzoeva OK, Das D, Heiser LM, Bhattacharya S, Siwak D, Gendelman R, et al. Basal subtype and MAPK/ERK kinase (MEK)-phosphoinositide 3-kinase feedback signaling determine susceptibility of breast cancer cells to MEK inhibition. Cancer Res. 2009;69(2):565–72. pmid:19147570
  48. 48. Haga Y, Higashisaka K, Yang L, Sekine N, Lin Y, Tsujino H, et al. Inhibition of AKT/mTOR pathway overcomes intrinsic resistance to dasatinib in triple-negative breast cancer. Biochem Biophys Res Commun. 2020;533(4):672–8. pmid:33036754
  49. 49. You KS, Yi YW, Kwak SJ, Seong YS. Inhibition of RPTOR overcomes resistance to EGFR inhibition in triple-negative breast cancer cells. Int J Oncol. 2018;52(3):828–40. pmid:29344641
  50. 50. You KS, Yi YW, Cho J, Seong YS. Dual inhibition of AKT and MEK pathways potentiates the anti‐cancer effect of gefitinib in triple‐negative breast cancer cells. Cancers (Basel). 2021;13(6):1–25. pmid:33801977
  51. 51. Xu T, Liu P, Li Q, Shi C, Wang X. Inhibitory effects of everolimus in combination with paclitaxel on adriamycin-resistant breast cancer cell line MDA-MB-231. Taiwan J Obstet Gynecol. 2020;59(6):828–34. pmid:33218396
  52. 52. Ma CX, Sanchez C, Gao F, Crowder R, Naughton M, Pluard T, et al. A phase I study of the AKT inhibitor MK-2206 in combination with hormonal therapy in postmenopausal women with estrogen receptor-positive metastatic breast cancer. Clin Cancer Res. 2016;22(11):2650–8. pmid:26783290
  53. 53. Li YL, Weng HC, Hsu JL, Lin SW, Guh JH, Hsu LC. The combination of MK-2206 and WZB117 exerts a synergistic cytotoxic effect against breast cancer cells. Front Pharmacol. 2019;10(1311). pmid:31780937
  54. 54. Wang YL, Overstreet AM, Chen MS, Wang J, Zhao HJ, Ho PC, et al. Combined inhibition of EGFR and c-ABL suppresses the growth of triple-negative breast cancer growth through inhibition of HOTAIR. Oncotarget. 2015;6(13):11150–61. pmid:25883211
  55. 55. Drago JZ, Formisano L, Juric D, Niemierko A, Servetto A, Wander SA, et al. FGFR1 amplification mediates endocrine resistance but retains torc sensitivity in metastatic hormone receptor-positive (HR+) breast cancer. Clin Cancer Res. 2019;25(21):6443–51. pmid:31371343
  56. 56. Larsen SL, Laenkholm AV, Duun-Henriksen AK, Bak M, Lykkesfeldt AE, Kirkegaard T. SRC drives growth of antiestrogen resistant breast cancer cell lines and is a marker for reduced benefit of tamoxifen treatment. PLoS One. 2015;10(2):e0118346. pmid:25706943
  57. 57. Ribas R, Pancholi S, Guest SK, Marangoni E, Gao Q, Thuleau A, et al. AKT antagonist AZD5363 influences estrogen receptor function in endocrine-resistant breast cancer and synergizes with fulvestrant (ICI182780) in vivo. Mol Cancer Ther. 2015;14(9):2035–48. pmid:26116361
  58. 58. Hutcheson IR, Goddard L, Barrow D, McClelland RA, Francies HE, Knowlden JM, et al. Fulvestrant-induced expression of ErbB3 and ErbB4 receptors sensitizes oestrogen receptor-positive breast cancer cells to heregulin β1. Breast Cancer Res. 2011;13(R29).
  59. 59. Park BJ, Whichard ZL, Corey SJ. Dasatinib synergizes with both cytotoxic and signal transduction inhibitors in heterogeneous breast cancer cell lines—lessons for design of combination targeted therapy. Cancer Lett. 2012;320(1):104–10. pmid:22306341
  60. 60. Liu T, Yacoub R, Taliaferro-Smith LTD, Sun SY, Graham TR, Dolan R, et al. Combinatorial effects of lapatinib and rapamycin in triple-negative breast cancer cells. Mol Cancer Ther. 2011;10(8):1460–9. pmid:21690228
  61. 61. Nikanjam M, Liu S, Yang J, Kurzrock R. Dosing three-drug combinations that include targeted anti-cancer agents: analysis of 37,763 patients. Oncologist. 2017;22(5):576–84. pmid:28424323
  62. 62. Tang J, Wennerberg K, Aittokallio T. What is synergy? the Saariselkä agreement revisited. Front Pharmacol. 2015;6:1–5.
  63. 63. Beal J, Montagud A, Traynard P, Barillot E, Calzone L. Personalization of logical models with multi-omics data allows clinical stratification of patients. Front Physiol. 2019;9: 1965. pmid:30733688
  64. 64. Béal J, Pantolini L, Noël V, Barillot E, Calzone L. Personalized logical models to investigate cancer response to BRAF treatments in melanomas and colorectal cancers. PLoS Comput Biol. 2021;17(1):e1007900 pmid:33507915