- Open Access
Systems biology approaches for advancing the discovery of effective drug combinations
Journal of Cheminformaticsvolume 7, Article number: 7 (2015)
Complex diseases like cancer are regulated by large, interconnected networks with many pathways affecting cell proliferation, invasion, and drug resistance. However, current cancer therapy predominantly relies on the reductionist approach of one gene-one disease. Combinations of drugs may overcome drug resistance by limiting mutations and induction of escape pathways, but given the enormous number of possible drug combinations, strategies to reduce the search space and prioritize experiments are needed. In this review, we focus on the use of computational modeling, bioinformatics and high-throughput experimental methods for discovery of drug combinations. We highlight cutting-edge systems approaches, including large-scale modeling of cell signaling networks, network motif analysis, statistical association-based models, identifying correlations in gene signatures, functional genomics, and high-throughput combination screens. We also present a list of publicly available data and resources to aid in discovery of drug combinations. Integration of these systems approaches will enable faster discovery and translation of clinically relevant drug combinations.
Despite increasing investments in pharmaceutical research and development, the rate of introduction of successfully translated drugs has decreased . Reasons for increased attrition rates in drug development include toxicity and inadequate efficacy due to individual variation in therapeutic response and development of drug resistance . This increased attrition rate has coincided with increased interest in seeking highly specific ligands affecting single targets for treatment of disease . Pharmaceutical research has increasingly relied on reductionist approaches, even though systemic diseases such as cancer and heart disease are managed by large, interconnected networks with many pathways affecting pathological signaling [4,5]. The redundancy and feedback in these networks allows for robustness of phenotype and maintenance of homeostasis [6,7]. This network complexity has hindered development of new therapies and indicates a need for more integrative systems approaches to make better predictions of drug responses [8,9].
The failure of single targets to successfully translate into clinical practice and the problem of development of drug resistance with single target cancer therapies has increased interest in discovery of effective drug combinations. Administering drug combinations has been effective in overcoming resistance to anti-microbial therapies for treatment of infectious diseases such as HIV and tuberculosis . In cancer, drug resistance can occur through mutation of the drug target , amplification of an alternate pathway , or intrinsic resistance of a subset of the cancer cells . Combinations of drugs could potentially overcome these resistance strategies by limiting the potential of escape mutations and pathways .
While combination therapies may dramatically improve efficacy of cancer therapies, the discovery of effective combinations is a challenging endeavor. With over 1,500 FDA approved compounds, experimentally testing every possible combination of these drugs would be unfeasible, even with high-throughput experimental methods . Therefore, new systems approaches are needed to reduce the search space and prioritize combinations for experimental testing (Figure 1) .
Here, we review computational and experimental methods for accelerating the discovery of effective drug combinations for complex diseases, with special focus on cancer. In addition, we include a list of publicly available resources as a reference for future drug combination studies (Table 1).
Quantifying synergistic drug combinations
When presenting the results of drug combination studies, it is important to have a standard to statistically define synergistic drug pairs. Two commonly used methods for quantifying synergy between drug combinations are Loewe additivity and Bliss independence. Loewe additivity is based on the assumption that the two inhibitors act through a similar mechanism while Bliss independence assumes independent mechanisms .
Using Loewe additivity, the concentration of two inhibitors (A and B) which alone results in X% inhibition of the target ([I A]X%, [I B]X%) can be used to calculate the theoretical concentrations of each inhibitor needed to achieve the same X% inhibition when combined ([C A]X%, [C B]X%).
The Loewe additivity applies the isobologram analysis to evaluate the combination effects of two drugs at a given effect. For example, in a Cartesian coordinate plot where x and y-axes represent concentrations of drugs A and B to achieve a defined effect X% (e.g., X = 50% for half maximal inhibitory concentration (IC50) of [I A]50% and [I B]50%), respectively. The coordinates ([I A]50%,0) and (0, [I B]50%) represent the concentration for drugs A and B, respectively. The line of additivity is constructed by connecting these two points for a 50% effect isobologram plot. The concentrations of the two drugs used in combination to provide the same effect X% (e.g. X = 50%) will be denoted by point ([C A]50%,[C B]50%) and are placed in the same plot. Synergy, additivity, or antagonism will be determined when this point ([C A]50%,[C B]50%) is located below, on, or above the line, respectively. More generally, linear, concave, and convex isoboles represent non-interacting, synergy, and antagonistic drug combination, respectively (Figure 2A).
This approach led to the development of the combination index (CI) popularized by Chou and Talalay . Here, the CI provides a quantitative measure of the extent of drug interaction at a given effect. It measures the combination concentrations of drugs A and B to produce a effect X%, [C A] and [C B], normalized by their corresponding concentrations that produces the same effect as a single agent, [I A] and [I B], respectively. CI value is calculated by:
where CI value <1, =1, and >1 represent synergy, additivity, and antagonism, respectively.
Bliss independence is based on probability theory and assumes the two inhibitors are working through independent mechanisms . The inhibitors do not interfere with each other, but contribute to a common result. Unlike Loewe additivity, calculating Bliss independence does not require determination of dose–response curves for the individual compounds to determine the theoretical results, making it easier to compute . Bliss independence models the combined effect (E T) as the product of the individual effects with drugs A (E A) and B (E B). The predicted combined effect (E T) is computed by:
where each effect (E) is expressed as fractional activity compared to control between 0 (100% inhibition) and 1 (0% inhibition). For example, if drug A and drug B each result in 40% tumor growth compared to control, then the predicted tumor growth when combined would be (0.4)*(0.4)*(100%) = 16% of control according to Bliss Independence. The predicted combined inhibition level would therefore be 100%-16% = 84% inhibition of tumor growth. If the actual tumor growth when drug A and B are combined is less than 16% of control (greater than 84% growth inhibition), then the compounds would be synergistic by Bliss Independence. If the tumor growth level is greater than 16% of control (less than 84% growth inhibition), then the compounds would be defined as antagonistic (Figure 2B).
These two methods produce different results, and it is uncertain which method performs better with uncertainty of mechanism and noisy data [17,20]. Drugs inhibiting parts of the same linear pathway may act according to Loewe additivity . Drugs nonexclusively affecting parallel pathways may act according to Bliss independence. Experimental characterization of drug combinations typically involves generating dose response curves with the inhibitors separately and combined. The experimental dose response curve data can then be compared to the predictions of Loewe additivity or Bliss independence to determine if the drugs are acting synergistically.
Computational models of signaling networks
Given the complexity of the signaling networks controlling systemic diseases such as cancer, computational models of cell signaling pathways are important tools for increasing understanding of pathological signaling and prioritizing targets to test experimentally . Models can be used to quantify systems properties that are often not apparent in individual experiments. Through model simulations, one can predict the relative importance of various proteins in the network, the presence of signal amplification, and the role of feedback and cross-talk . These features will be important in the prediction of viable drug combinations. While model predictions require experimental validation, they are useful tools for prioritizing targets for experimental planning.
Mass action and enzyme kinetics-based models
Three predominant signaling network modeling approaches are mass action and enzyme kinetics-based, logic-based, and statistical association-based models. Mass action models are biochemically detailed kinetic models that typically represent interactions between molecular species in the signaling network as ordinary differential equations (ODEs) and require selection of parameter values for concentrations of species in the network and rate constants controlling protein-protein associations . Many of these parameters may be unavailable in the literature and can either be measured experimentally or fit to the data by minimizing an objective function such as the sum of squared error.
When training parameters to data, it is important to determine the importance of parameter selection on model predictions. For example, Chen et al. measured parameter sensitivity for several independent fits and saw that the rank order of the most sensitive parameters was nearly the same across the fits for a given output, therefore parameter uncertainty did not affect major model predictions . In another approach, Iadevaia et al. developed a mass-action model of the IGF-1 signaling network in a breast cancer cell line with 161 unknown parameters and fit the model to the time courses of six proteins measured with reverse-phase protein array . Given the uncertainty in parameter estimation with so many unknown parameters, they identified ten sets of parameters using particle swarm optimization that equally fit the experimental data. Model predictions were averaged from three randomly sampled sets of the ten parameter sets. The trained model was then used to identify beneficial drug combinations in a breast cancer cell line.
Mass-action network models have been used to predict new beneficial drug combinations for cancer. As an example, Faratian et al. used a mass-action model of heregulin-induced HER2/3 signaling through MAPK and PI3K to study the role of PIK3CA activation in Receptor Tyrosine Kinase (RTK) inhibitor resistance . Model results demonstrated that the ratio of PTEN to activated PIK3CA predicted resistance to RTK inhibitors. This finding could therefore be used to predict patient response to anti-HER2 therapies based on clinical measurements of PTEN. It predicts that PIK3CA inhibition should be paired with RTK inhibitors in patients with tumors with low PTEN, a negative regulator of PI3K signaling. Another group developed a mass action kinetics model of PI3K signaling by ERBB receptors including ligand binding, dimerization, internalization, recycling, and degradation . Sensitivity analysis of this model predicted an important role of ERBB3 in AKT activation, which was then validated in mice xenografts. Sensitivity analysis could be used in future work to find drug combinations that may work synergistically with ERBB3 inhibition.
A limitation of mass-action modeling approaches is the amount of data required to generate specific values for the abundance and rate constant parameters, which can be prohibitive for large scale network reconstructions (Figure 1) . Logic-based models use network topology without the need for specific parameter values. Network interactions are modeled with OR, AND, and NOT Boolean logic gates. Each species in the network takes a value of 0 (inactive) or 1 (active) based on the state of its effectors . As an example, Sahin et al. developed a Boolean model of ERBB signaling of G1/S cell cycle transition . The group used computational knockouts of network proteins, validation experiments with RNAi, and model revision based on proteomic data, to predict the effects of combined inhibition of ERBB2 and c-MYC or EGFR. A combination therapy targeting c-MYC and ERBB2 was predicted to improve treatment for breast cancer that is de novo resistant to ERBB2 inhibition. Another group developed a Boolean logic model of apoptosis signaling in Leukemic T-Cell large granular lymphocytes . The authors used the model to determine species that controlled apoptosis and experimentally validated two of these species, sphingosine kinase 1 and NFκB. Given the limitations in representing species as either on or off, this modeling approach has been extended to accommodate intermediate activity states using fuzzy logic .
Normalized hill differential equation modeling approach
While logic-based modeling approaches benefit from simple construction using network topology, results can be difficult to interpret due to assignment of discrete values to continuous variable such as concentration of active species, sensitivity to temporal node-updating schemes, and incompatibility with many systems analysis tools such as quantitative sensitivity analysis . To address the limitations of mass-action and logic-based models, Kraeutler et al. developed the normalized Hill differential equation modeling approach, which uses logic-based differential equations to represent activation or inhibition by molecular species in the network . Cross-talk is represented with AND and OR gates and species activation is continuous over time and in units of fractional activation instead of concentration. Therefore protein abundance parameters are not required like with mass-action models. Interactions between species in the network are modeled with normalized Hill equations with 3 parameters: reaction weight, half maximal effective concentration (EC50), and Hill coefficient. While these parameters can be fit to data, using default values generated highly similar quantitative predictions as a previously constructed detailed biochemical model of the same pathway which used 88 parameters from literature [33,34]. Therefore, this approach allows for straightforward model construction of a known network topology even if kinetic and abundance parameters are unknown, like with logic-based modeling, while also allowing for prediction of dynamics and systems analysis tools such as quantitative sensitivity analysis.
The normalized-Hill modeling approach is a valuable tool for model construction of larger networks with more unknown parameters. As an example, Ryall et al. used this approach to model the cardiac hypertrophy signaling network, which contained 106 species and 193 reactions . Since cardiac myocytes have minimal capacity for proliferation, many of these pathways also regulate proliferation in cancer cells . Quantitative systems analysis revealed the most prevalent species involved in growth of cardiac myocytes, prioritizing future experimental targets . While Ras, the largest signaling hub, was the highest influencer on cell size, the correlation between the number of connections a species has and its influence was low. Moreover, highly influential species were at many levels in the network, not just close to the output level. These findings demonstrate the need for model reconstructions to predict important drug targets in cell signaling networks. Highly influential species are not obvious from intuition alone or data from gain or loss of function studies of single genes .
Ryall et al.’s analysis of the hypertrophy signaling network also looked at the presence of different signaling motifs such as bi-fan and feed-forward loops. Motifs can affect network properties such as signal filtering, acceleration, pulse generation, ultra-sensitivity, stability, and robustness [38-40]. Yin et al. modeled three-node enzymatic networks with many different topologies to study the effect of topology on drug combinations . Model simulations were conducted to identify motifs that could result in synergy. Most of the combinations were not dependent on parameter selection, demonstrating that network topology can be used to predict synergistic combinations. Moreover, synergistic drug combinations were found in both parallel and series drug combinations. In a similar study, Zhang et al. made reduced models of the convergence of two signaling pathways on a target and observed synergy in only a subset of the motifs . Synergy had a greater likelihood in motifs with negative feedback between the target and an upstream effector or mutual inhibition between parallel signaling pathways. These findings suggest that searching for synergistic motifs within a cancer signaling network topology can be a useful strategy in prioritizing drug combinations to test experimentally. Networks exported into Cytoscape , a open source software platform for network visualization, can use the Netmatch plug-in  to quickly search for motifs of interest.
Statistical association-based modeling approach
Network modeling approaches are useful when network topology is known, but these approaches can be biased towards established pathways and may miss novel interactions. Statistical association-based models do not depend on prior knowledge of pathways and instead use correlations and patterns in experimental data to predict network structure. As an example, Ryall et al. used data of correlations among cell shape features and expression of 12 genes relevant to cardiac hypertrophy to identify a network map linking input modules to output modules . Drug combinations could then be prioritized by selecting targets that enabled adaptive module signaling and prevented maladaptive module signaling. Molinelli et al. developed a network inference algorithm based on Belief Propagation  to construct networks from phenotypic screen data . They applied their method to screen data from a melanoma cell line and identified both new and established pathway interactions and then used the network to predict efficacious drug targets. Another useful approach for network inference is Bayesian network computational methods .
Signature-based approaches for predicting drug combinations
Many effective drug combinations have been discovered using correlations in gene expression signatures. One useful tool for this is the Connectivity Map (CMap) database . The first-generation CMap contained gene expression profiles from three cancer cell lines perturbed by 164 distinct small-molecule compounds. The second generation of CMap (CMap 2.0) includes gene expression profiles from 1,309 small molecules including FDA approved drugs tested in five human cancer cell lines [49,50]. This method assumes gene expression changes can be used as a “universal language” to connect distinct biological states (e.g. diseases), allowing for the successful repurposing of compounds [51,52]. In short, drugs known to be effective in one disease can serve as candidates for use in other diseases marked by similar gene expression changes. Users query the database to compute similarity metrics between a test gene expression signature and each reference set. Similarity metrics are scaled between −1 and +1, where a positive score indicates positive correlation and a negative score indicates negative correlation. An advantage of this approach is the ability to query CMap with publicly available gene expression data from sources such as Gene Expression Omnibus (GEO) , therefore facilitating rapid drug combinations prediction for experimental validation .
As an example, Riedel et al. applied CMap to predict drugs that would prevent resistance to chemotherapy agents in lung cancer cell lines . Genes with the highest changes after treatment with docetaxel were analyzed using CMap to identify drugs with negative connectivity scores, indicating these drugs had antagonistic effects on the genes associated with docetaxel resistance. PI3K inhibitor LY294002, which was highly ranked among these antagonistic compounds, was tested in vitro with docetaxel and found to synergistically increase cytotoxicity. Wei et al. used a similar approach to predict drugs to overcome resistance to glucocorticoid treatment in acute lymphoblastic leukemia . Microarrays of pre-treated cell lines either sensitive or resistant to glucocorticoid in vitro were used to define a sensitive/resistant gene signature. Using CMap, mTOR inhibitor rapamycin was found to induce a highly similar signature, leading to the hypothesis that mTOR inhibition could induce glucocorticoid sensitivity. Follow-up experiments supported this hypothesis and showed that rapamycin conferred sensitivity through down-regulation of MCL1. Therefore, CMap is a useful tool for using gene signatures to predict drug combinations that may overcome drug resistance.
Inspired by the CMap concept, Kim et al. recently developed K-Map (Kinase Inhibitor Connectivity Map) that systematically connects a set of query kinases to kinase inhibitors based on quantitative profiles of the kinase inhibitor activities . Instead of gene expression signatures, Kim et al. used the kinase activity profiles as the “language” for connecting kinases and small molecules in K-Map to reveal the complex interactions of kinases and inhibitors. By querying K-Map with the essential kinases mediating resistance to EGFR-inhibitor gefitinib in an EGFR mutant non-small cell lung cancer (NSCLC) cell line, bosutinib was predicted to be a more effective drug for killing EGFR mutant cancer cells. Follow up in vitro experiments confirmed that bosutinib alone is a more effective agent than gefitinib, and that the combination of bosutinib and gefitinib had synergistic effects in EGFR mutant NSCLC cells . This demonstrates the utility of K-Map in connecting kinases with kinase inhibitors and suggesting candidates for drug combinations.
Network-based approaches for predicting drug combinations
Other computational approaches have been developed to predict drug combinations using data from high-throughput screens and drug databases. Pal and Berlow developed an algorithm based on set theory that uses tumor drug sensitivities and kinase inhibition profiles for a set of individual drugs to predict the tumor sensitivity to new drugs or drug combinations . The algorithm applies the following rules to generate circuit representations of tumor pathways: 1) drugs that inhibit a superset of an effective set of inhibited kinases will also be successful in inhibiting tumor growth and 2) drugs inhibiting subsets of ineffective sets of inhibited kinases will also be unsuccessful. These circuits reveal a set of kinases that are most predictive of drug sensitivity and depict combinations of kinases that need to be inhibited to prevent tumor growth. This analysis is helpful for identifying drug combinations that inhibit a minimal set of kinases with as few off target effects as possible to minimize negative side effects. This approach was validated using data from four canine cancer cell lines given 60 different drugs at four different concentrations to generate IC50 values . Tang et al. expanded on this algorithm to improve the computational cost and accuracy with drug screen data with little overlap between drug target profiles . While this approach requires a lot of experimental data from drug screens to be useful (Figure 1), technological advancements are enabling larger drug screens at lower costs.
Given the high cost of exhaustive drug screens, Gujral et al. exploited the polypharmacology of kinase inhibitors by developing an approach to select the most predictive kinase targets from a smaller scale drug screen of multi-target drugs . They performed a phenotypic screen to identify kinases regulating cell migration using an optimal set of 32 kinase inhibitors. Elastic net regularization was then used to deconvolute the polypharmacology of the kinase inhibitors, identifying kinases with the greatest explanatory power for the phenotype. Elastic net regularization regresses an output variable against a set of predictor variables (kinase activity) and invokes a penalty on the number of variables in order to eliminate kinases with insignificant contributions. This is a useful approach for reducing experimental time and cost by extracting more information from a smaller set of compounds in drug screens.
Huang et al. used drug genomic profiles from the Connectivity Map database to construct a drug functional network and then grouped drugs into modules with similar transcriptional responses . They then built disease-signaling networks highlighting defective signaling modules based on patient genomic profiles and protein interaction data. The DrugComboRanker algorithm ranked potential drug combinations by selecting drugs with high overlap in the disease network, affecting multiple key signaling modules. In a similar approach, Pang et al. developed an algorithm to identify combination therapies by building a network of drug-target interactions from the DrugBank database . Given an input disease gene set, the algorithm selected drugs that maximized on target effects and minimized off target effects. This algorithm also identifies drug combinations of more than two drugs, which would be unfeasible to predict using a high-throughput screen. These approaches allow researchers to take advantage of publicly available drug data to prioritize combinations for experimental validation.
Moreover, Liu et al. have developed a database of both successful and unsuccessful drug combinations (DCDB) . The current version of the database contains 1,363 drug combinations involving 904 individual drugs and 805 targets. The database provides information about the potential mechanism, drug interactions, indication, published study, development status, and targets. The ability to analyze patterns in successful and unsuccessful combinations with this database will be useful for systems analysis of drug combinations and rational experimental design. As an example, Xu et al. constructed a network of successful drug interactions using data from the Drug Combination Database . Analysis revealed that drugs with similar therapeutic effects tended to cluster together in the network and targets of hub drugs were often membrane or membrane-associated proteins. They used these observations to develop a statistical approach to predict new drug combinations.
Cheng et al. constructed a global human kinome interaction map by integrating kinase-substrate interactions, kinase-drug interactions, protein-protein interactions, and atomic resolution three-dimensional structural protein-protein interactions . Their analysis of the topological features of these networks revealed an enrichment of hubs as drug targets. While targeting hubs can be beneficial due to cascading effects, their analysis revealed that targeting hubs also increases risk of adverse drug reactions and drug resistance through feedback and crosstalk.
Network based approaches have also been used to study drug-drug interactions. As an example, Cheng and Zhao developed a comprehensive drug-drug interaction network incorporating 6946 interactions of 721 approved drugs using data from DrugBank . They then calculated drug-drug pair similarities using four features: phenotypic similarity, therapeutic similarity, chemical structure similarity, and genomic similarity. They applied five machine learning-based models to the dataset to predict drug-drug interaction, with the overall hypothesis that drugs with similar chemical structure, target proteins, adverse drug reactions, and therapeutic purposes have high probability of drug-drug interaction. They tested the model on antipsychotic drug-drug interactions and found literature support for predictions of drug-drug interactions involving weight gain and P450 inhibition. This approach demonstrates the power of harnessing network-based drug-drug interactions to reveal new information on adverse drug effects and provide additional filtering rules for drug combination studies.
Integrating functional genomics and computational methods for identifying drug combinations
Large scale knockdown screens using RNA interference (RNAi) can also be used to identify potential drug combinations . RNAi screens can identify genes that lead to sensitivity or resistance to a drug of interest. As an example, an RNAi screen conducted by Berns et al. showed that knockdown of PTEN decreased sensitivity to trastuzumab in BT-474 breast cancer cells . Follow-up studies showed that assessment of both the loss of PTEN expression and activating mutations in PIK3CA could predict the risk for HER2 amplified tumor progression. Drugs reducing PI3K signaling may therefore increase response to trastuzumab. In another study, Prahallad et al. used an RNAi screen to identify kinases whose knockdown synergized with BRAF(V600E) inhibition in colon tumors . Follow-up experiments demonstrated synergy between cetuximab (EGFR inhibitor) and vemurafenib (BRAF inhibitor). The rational combination of cetuximab and vemurafenib is currently being evaluated in clinical trials.
Pritchard et al. used RNAi signatures of eight cell death genes to determine the mechanism of drug combination effects in lymphoma cells . Single drugs were classified based on their similarity to the RNAi signatures of well-characterized compounds with known mechanisms. They then generated signatures for drug combinations to see if the signature was more similar to results from one of the drugs alone, an average of the two, or a unique signature. Results showed that the combination signature was usually a weighted composite of single drug effects where one drug potentiated the mechanism of the other or the two drugs produced an additive effect. Interestingly, they observed that applying larger pools of drugs to tumors reduced the genetic heterogeneity, which could be prohibitive in selection of personalized drug treatments for patients based on biomarkers.
Spreafico et al. identified non-canonical Wnt pathway mediated resistance to MEK1/2 inhibitor Selumetinib in colorectal cancer cells by integrating gene set enrichment analysis and synthetic lethality screens . Using cyclosporine A (CsA) as a drug to inhibit non-canonical Wnt pathway, they validated that the combination of Selumetinib and CsA has synergistic anti-proliferative effects both in vitro and in vivo patient-derived xenografts. This rational combination is now being translated into a Phase I clinical study (ClinicalTrials.gov ID: NCT02188264). This illustrates the utility of integrating functional genomics screens with bioinformatics to identify and translate drug combinations into clinical study.
High-throughput drug combination screens
While the number of FDA approved drugs makes exhaustive drug screens unfeasible (Figure 1), there have been efforts to reduce the search space in an unbiased manner. Tan et al. used pools of ten drugs in 384-well plates to study all possible pairs of 1,000 compounds in the minimum number of wells possible in order to find drugs combinations that synergistically prevent HIV replication . Synergistic wells from the primary screen are then deconvolved into possible synergistic pairs for a secondary screen. Results revealed an enrichment of anti-inflammatory drugs in the synergistic pairs.
Roller et al. conducted a functional chemical genetic screen of 300 drug combinations in nine melanoma cell lines and identified pairs of compounds that synergistic increase cytotoxicity . Interestingly, the synergistic cytotoxicities identified did not correlate with the known oncogene RAS and RAF mutational status of the melanoma cell lines. From this screen, they identified sorafenib (a multi-kinase inhibitor) and diclofenac (a non-steroidal anti-inflammatory drug) to be the most robust drug combination that has synergistic effects across the melanoma cell lines. By using this functional chemical genetic screen, the authors uncovered novel interactions between signaling inhibitors that would not be predicted based on current understanding of the signaling networks. Their results suggest that the underlying signaling networks controlling drug responses can vary substantially based on unidentified elements of cell genotype. In another study, Griner et al. conducted a large scale screen of multiple concentrations of 500 compounds with ibrutinib in activated B-cell-like diffuse large lymphoma cells . They discovered many compounds that interacted favorably with ibrutinib, including inhibitors of PI3K signaling, the Bcl2-family, and the B-cell receptor pathway.
One of the major ongoing initiatives at the Developmental Therapeutics Program of the National Cancer Institute, U.S. National Institute of Health, is the large-scale high-throughput drug combination screening of 100 FDA approved drugs. The first set of screening results generated 5,000 possible drug combinations in the NCI-60 cancer cell lines panel . The goal of this project is to identify novel drug combinations that are more active than the single agents alone. As all the drugs tested have been FDA approved, any drug combinations identified from this screen may rapidly translate into the clinic. As the NCI-60 cell lines panel have been fully characterized by various “omics” technologies, the release of this drug combination matrix to the public could facilitate the development of novel computational methods to integrate, predict, and mine the interactions between molecular markers and drug combinations.
While substantial intratumoral heterogeneity has been detected in cancer patients using next generation sequencing technologies , current drug combination prediction methods have primarily focused on targeting the predominant tumor subpopulation. To study the effect of different tumor subpopulations on treatment efficacy, Zhao et al. developed a multi-objective linear optimization algorithm to select optimal drug combinations for heterogeneous tumors by maximizing efficacy and minimizing toxicity . Their goal was to determine the best drug combinations to minimize all subpopulations. They experimentally validated the algorithm’s prediction of two-drug combinations with three-component heterogeneous tumors created using RNA interference . They then expanded the model to simulate more complex tumors and greater numbers of drugs . Their results revealed that intratumor heterogeneity influences the prediction of effective drug combinations. Different predictions are made depending on if all tumor subpopulations are considered or just the predominant subpopulation. This approach represents a step forward of predicting drug combinations to tackle tumor heterogeneity in the era of precision oncology.
Given the experimental costs of exhaustively testing drug combinations, computational models of signaling networks will be especially useful in pre-clinical screening of combinations of compounds. Model simulations reveal non-intuitive effects of drug combinations . Due to the size and complexity of cancer signaling, modeling strategies accommodating reconstruction of larger networks while still being compatible with quantitative systems analysis tools will be especially useful . Global network models allow for more comprehensive and unbiased discovery of therapeutic targets than experimental approaches based on a priori selection of important pathways . One of these methods is the recently described normalized-Hill modeling approach, which even with minimal biochemical data from literature enables a global view of quantitative functional relationships between every node in a network . This method was used to identify the most important nodes in a integrative network model of cardiac hypertrophy with 106 nodes and 193 reactions and default parameter values . Many of these cardiac hypertrophy signaling pathways also play important roles in tumor growth. While model predictions need to be experimentally validated, models can substantially improve hypothesis generation and experimental planning.
In addition to larger network reconstructions, future modeling efforts will benefit from tighter integration of high-throughput sequencing, proteomics, and phenotypic screen data (Figure 3A). This will enable tuning of a model to an individual patient’s tumor, which would be beneficial for use in personalized medicine. Comprehensive double knockdown model simulations would enable personalized prediction of drug combinations for patients. The ability to readily adapt models to different situations is important in cancer research since the molecular networks are not fixed within a particular cancer type . Patients that share the same mutation and tumor type can have different responses to a drug . Genetic background, cell lineage, and exogenous signals can influence the network behavior . Efficacy data identified from in vitro and in vivo experiments would then be used for model refinement so more informed predictions of drug combinations can be made in future studies.
Predicted drug combinations should be validated in cancer cell lines and in relevant in vivo human disease models such as patient-derived tumor xenografts . These models, however, typically overestimate the clinical benefit due to factors such as tumor heterogeneity, differences in tumor microenvironment, and inaccurate estimates of drug exposure . Therefore, it is important to have a high threshold when choosing effective combinations, ignoring modest inhibitions of tumor growth in favor of combinations promoting cancer cell death and tumor regression.
Another opportunity for improved design of combination therapies is through quantitative systems pharmacology approaches integrating cell signaling network models with pharmacokinetic-pharmacodynamic (PK/PD) models . Quantitative systems pharmacology uses multi-scale data to better understand and ultimately predict how drugs affect cellular networks and human pathophysiology . Mechanistic models of cell signaling networks are linked to PK/PD models of physiological processes at the level of tissues and organisms. These models will enable patient-specific prediction of therapeutic and toxic drug responses and drug resistance mechanisms, improve translation of in vitro discoveries to patients, and enhance discovery of pharmacodynamic biomarkers.
Additionally, future work will benefit from parallel integration of computational modeling, preclinical testing, and clinical trials, where data from each approach can be used for refinement of the other. Computational models tuned to specific cancer cell lines using bioinformatics and experimental data could be perturbed to make predictions of effective drug combinations to validate in preclinical models (Figure 3A). For clinical application, similarity scores between patients and previously modeled cell lines could be calculated using statistical clustering (Figure 3B). Drug combinations predicted using the most similar model could then be applied in the clinic.
It is becoming increasingly apparent that drug combinations will be essential for improving therapies for complex diseases such as cancer . The signaling pathways controlling these systemic diseases are highly interconnected, with cross-talk, redundancy, and feedback, making single-target therapies much less effective . While combination therapies have the potential to prevent the development of resistance seen in many single drug therapies, it is prohibitively expensive to experimentally test every potential combination, especially when considering combinations of more than two drugs. Here, we highlighted a variety of systems biology applications for advancing the prediction of effective drug combinations, as summarized in Table 2. These methods include computational modeling, gene signature analysis, functional genomics, and high-throughput drug combination screening. Utilization and integration of these systems biology approaches hold great promise in speeding up the development of clinically relevant drug combinations.
Pammolli F, Magazzini L, Riccaboni M. The productivity crisis in pharmaceutical R&D. Nat Rev Drug Discov. 2011;10(6):428–38.
Hopkins AL. Network pharmacology: the next paradigm in drug discovery. Nat Chem Biol. 2008;4(11):682–90.
Sams-Dodd F. Target-based drug discovery: is something wrong? Drug Discov Today. 2005;10(2):139–47.
Heineke J, Molkentin JD. Regulation of cardiac hypertrophy by intracellular signalling pathways. Nat Rev Mol Cell Biol. 2006;7(8):589–600.
Breitkreutz D, Hlatky L, Rietman E, Tuszynski JA. Molecular signaling network complexity is correlated with cancer patient survivability. Proc Natl Acad Sci U S A. 2012;109(23):9209–12.
Kitano H. Cancer as a robust system: implications for anticancer therapy. Nat Rev Cancer. 2004;4(3):227–35.
Albert R, Jeong H, Barabasi AL. Error and attack tolerance of complex networks. Nature. 2000;406(6794):378–82.
Lusis AJ, Weiss JN. Cardiovascular networks: systems-based approaches to cardiovascular disease. Circulation. 2010;121(1):157–70.
Berger SI, Iyengar R. Network analyses in systems pharmacology. Bioinformatics. 2009;25(19):2466–72.
Glickman MS, Sawyers CL. Converting cancer therapies into cures: lessons from infectious diseases. Cell. 2012;148(6):1089–98.
Gorre ME, Mohammed M, Ellwood K, Hsu N, Paquette R, Rao N, et al. Clinical resistance to STI-571 cancer therapy caused by BCR-ABL gene mutation or amplification. Science. 2001;293(5531):876–80.
Engelman JA, Zejnullahu K, Mitsudomi T, Song Y, Hyland C, Park JO, et al. MET amplification leads to gefitinib resistance in lung cancer by activating ERBB3 signaling. Science. 2007;316(5827):1039–43.
Sharma SV, Lee DY, Li B, Quinlan MP, Takahashi F, Maheswaran S, et al. A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations. Cell. 2010;141(1):69–80.
Weinstein IB, Joe A. Oncogene addiction. Cancer Res. 2008;68(9):3077–80.
Sun X, Vilar S, Tatonetti NP. High-throughput methods for combinatorial drug discovery. Sci Transl Med. 2013;5(205):205rv1.
Tang J, Aittokallio T. Network pharmacology strategies toward multi-target anticancer therapies: from computational models to experimental design principles. Curr Pharm Des. 2014;20(1):23–36.
Fitzgerald JB, Schoeberl B, Nielsen UB, Sorger PK. Systems biology and combination therapy in the quest for clinical efficacy. Nat Chem Biol. 2006;2(9):458–66.
Chou T, Talalay P. Analysis of combined drug effects: a new look at a very old problem. Trends Pharmacol Sci. 1983;4:450–4.
Keith CT, Borisy AA, Stockwell BR. Multicomponent therapeutics for networked systems. Nat Rev Drug Discov. 2005;4(1):71–8.
Goldoni M, Johansson C. A mathematical approach to study combined effects of toxicants in vitro: evaluation of the Bliss independence criterion and the Loewe additivity model. Toxicol In Vitro. 2007;21(5):759–69.
Kreeger PK, Lauffenburger DA. Cancer systems biology: a network modeling perspective. Carcinogenesis. 2010;31(1):2–8.
Klipp E, Liebermeister W. Mathematical modeling of intracellular signaling pathways. BMC Neurosci. 2006;7 Suppl 1:S10.
Schoeberl B, Eichler-Jonsson C, Gilles ED, Müller G. Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors. Nat Biotechnol. 2002;20(4):370–5.
Chen WW, Schoeberl B, Jasper PJ, Niepel M, Nielsen UB, Lauffenburger DA, et al. Input–output behavior of ErbB signaling pathways as revealed by a mass action model trained against dynamic data. Mol Syst Biol. 2009;5(239):239.
Iadevaia S, Lu Y, Morales FC, Mills GB, Ram PT. Identification of optimal drug combinations targeting cellular networks: integrating phospho-proteomics and computational network analysis. Cancer Res. 2010;70(17):6704–14.
Faratian D, Goltsov A, Lebedeva G, Sorokin A, Moodie S, Mullen P, et al. Systems biology reveals new strategies for personalizing cancer medicine and confirms the role of PTEN in resistance to trastuzumab. Cancer Res. 2009;69(16):6713–20.
Schoeberl B, Pace EA, Fitzgerald JB, Harms BD, Xu L, Nie L, et al. Therapeutically targeting ErbB3: a key node in ligand-induced activation of the ErbB receptor-PI3K axis. Sci Signal. 2009;2(77):ra31.
Papin JA, Palsson BO. Topological analysis of mass-balanced signaling networks: a framework to obtain network properties including crosstalk. J Theor Biol. 2004;227(2):283–97.
Morris MK, Saez-Rodriguez J, Sorger PK, Lauffenburger DA. Logic-based models for the analysis of cell signaling networks. Biochemistry. 2010;49(15):3216–24.
Sahin O, Fröhlich H, Löbke C, Korf U, Burmester S, Majety M, et al. Modeling ERBB receptor-regulated G1/S transition to find novel targets for de novo trastuzumab resistance. BMC Syst Biol. 2009;3:1.
Zhang R, Shah MV, Yang J, Shah MV, Yang J, Nyland SB, et al. Network model of survival signaling in large granular lymphocyte leukemia. Proc Natl Acad Sci U S A. 2008;105(42):16308–13.
Aldridge BB, Saez-Rodriguez J, Muhlich JL, Sorger PK, Lauffenburger DA. Fuzzy logic analysis of kinase pathway crosstalk in TNF/EGF/insulin-induced signaling. PLoS Comput Biol. 2009;5(4):e1000340.
Kraeutler MJ, Soltis AR, Saucerman JJ. Modeling cardiac beta-adrenergic signaling with normalized-Hill differential equations: comparison with a biochemical model. BMC Syst Biol. 2010;4(1):157.
Saucerman JJ, Brunton LL, Michailova AP, McCulloch AD. Modeling beta-adrenergic control of cardiac myocyte contractility in silico. J Biol Chem. 2003;278(48):47997–8003.
Ryall KA, Holland DO, Delaney KA, Kraeutler MJ, Parker AJ, Saucerman JJ. Network reconstruction and systems analysis of cardiac myocyte hypertrophy signaling. J Biol Chem. 2012;287(50):42259–68.
Molkentin JD, Dorn GW. Cytoplasmic signaling pathways that regulate cardiac hypertrophy. Annu Rev Physiol. 2001;63:391–426.
Kestler HA, Wawra C, Kracher B, Kühl M. Network modeling of signal transduction: establishing the global view. Bioessays. 2008;30(11–12):1110–25.
Alon U. Biological networks: the tinkerer as an engineer. Science. 2003;301(5641):1866–7.
Alon U. Network motifs: theory and experimental approaches. Nat Rev Genet. 2007;8(6):450–61.
Pelaez N, Carthew RW. Biological robustness and the role of microRNAs: a network perspective. Curr Top Dev Biol. 2012;99:237–55.
Yin N, Ma W, Pei J, Ouyang Q, Tang C, Lai L. Synergistic and antagonistic drug combinations depend on network topology. PLoS One. 2014;9(4):e93960.
Zhang Y, Smolen P, Baxter DA, Byrne JH. Computational analyses of synergism in small molecular network motifs. PLoS Comput Biol. 2014;10(3):e1003524.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.
Ferro A, Giugno R, Pigola G, Pulvirenti A, Skripin D, Bader GD, et al. NetMatch: a Cytoscape plugin for searching biological networks. Bioinformatics. 2007;23(7):910–2.
Ryall KA, Bezzerides VJ, Rosenzweig A, Saucerman JJ. Phenotypic screen quantifying differential regulation of cardiac myocyte hypertrophy identifies CITED4 regulation of myocyte elongation. J Mol Cell Cardiol. 2014;72:74–84.
Yedidia JS, Freeman WT, Weiss Y. Understanding belief propagation and its generalizations. In: Exploring Artificial Intelligence in the New Millennium. 2002. p. 239–69.
Molinelli EJ, Korkut A, Wang W, Miller ML, Gauthier NP, Jing X, et al. Perturbation biology: inferring signaling networks in cellular systems. PLoS Comput Biol. 2013;9(12):e1003290.
Sachs K, Perez O, Pe’er D, Lauffenburger DA, Nolan GP. Causal protein-signaling networks derived from multiparameter single-cell data. Science. 2005;308(5721):523–9.
Lamb J, Crawford ED, Peck D, Modell JW, Blat IC, Wrobel MJ, et al. The connectivity map: using gene-expression signatures to connect small molecules, genes, and disease. Science. 2006;313(5795):1929–35.
Qu XA, Rajpal DK. Applications of connectivity map in drug discovery and development. Drug Discov Today. 2012;17(23–24):1289–98.
Hieronymus H, Lamb J, Ross KN, Peng XP, Clement C, Rodina A, et al. Gene expression signature-based chemical genomic prediction identifies a novel class of HSP90 pathway modulators. Cancer Cell. 2006;10(4):321–30.
Wei G, Twomey D, Lamb J, Schlis K, Agarwal J, Stam RW, et al. Gene expression-based chemical genomics identifies rapamycin as a modulator of MCL1 and glucocorticoid resistance. Cancer Cell. 2006;10(4):331–42.
Edgar R, Domrachev M, Lash AE. Gene Expression Omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res. 2002;30(1):207–10.
Sirota M, Dudley JT, Kim J, Chiang AP, Morgan A, Sweet-Cordero A, et al. Discovery and preclinical validation of drug indications using compendia of public gene expression data. Sci Trans Med. 2011;3(96):96ra77.
Riedel RF, Porrello A, Pontzer E, Chenette EJ, Hsu DS, Balakumaran B, et al. A genomic approach to identify molecular pathways associated with chemotherapy resistance. Mol Cancer Ther. 2012;11(5):1214–5.
Kim J, Yoo M, Kang J, Tan AC. K-Map: connecting kinases with therapeutics for drug repurposing and development. Hum Genomics. 2013;7(1):20.
Kim J, Vasu VT, Mishra R, Singleton KR, Yoo M, Leach SM, et al. Bioinformatics-driven discovery of rational combination for overcoming EGFR-mutant lung cancer resistance to EGFR therapy. Bioinformatics. 2014;30(17):2393–8.
Pal R, Berlow N. A kinase inhibition map approach for tumor sensitivity prediction and combination therapy design for targeted drugs. Pac Symp Biocomput. 2012;351–62.
Berlow N, Davis LE, Cantor EL, Séguin B, Keller C, Pal R. A new approach for prediction of tumor sensitivity to targeted drugs based on functional data. BMC Bioinformatics. 2013;14(1):239.
Tang J, Karhinen L, Xu T, Szwajda A, Yadav B, Wennerberg K, et al. Target inhibition networks: predicting selective combinations of druggable targets to block cancer survival pathways. PLoS Comput Biol. 2013;9(9):e1003226.
Gujral TS, Peshkin L, Kirschner MW. Exploiting polypharmacology for drug target deconvolution. Proc Natl Acad Sci U S A. 2014;111(13):5048–53.
Huang L, Li F, Sheng J, Xia X, Ma J, Zhan M, et al. DrugComboRanker: drug combination discovery based on target network analysis. Bioinformatics. 2014;30(12):i228–36.
Pang K, Wan Y, Choi WT, Donehower LA, Sun J, Pant D, et al. Combinatorial therapy discovery using mixed integer linear programming. Bioinformatics. 2014;30(10):1456–63.
Liu Y, Hu B, Fu C, Chen X. DCDB: drug combination database. Bioinformatics. 2010;26(4):587–8.
Xu KJ, Song J, Zhao XM. The drug cocktail network. BMC Syst Biol. 2012;6 Suppl 1:S5.
Cheng F, Zhao Z. Machine learning-based prediction of drug-drug interactions by integrating drug phenotypic, therapeutic, chemical, and genomic properties. J Am Med Inf Assoc. 2014;21(e2):e278–86.
Iorns E, Lord CJ, Turner N, Ashworth A. Utilizing RNA interference to enhance cancer drug discovery. Nat Rev Drug Discov. 2007;6(7):556–68.
Berns K, Horlings HM, Hennessy BT, Madiredjo M, Hijmans EM, Beelen K, et al. A functional genetic approach identifies the PI3K pathway as a major determinant of trastuzumab resistance in breast cancer. Cancer Cell. 2007;12(4):395–402.
Prahallad A, Sun C, Huang S, Nicolantonio FD, Salazar R, Zecchin D, et al. Unresponsiveness of colon cancer to BRAF(V600E) inhibition through feedback activation of EGFR. Nature. 2012;483(7387):100–3.
Pritchard JR, Bruno PM, Gilbert LA, Capron KL, Lauffenburger DA, Hemann MT. Defining principles of combination drug mechanisms of action. Proc Natl Acad Sci U S A. 2012;110(2):E170–9.
Spreafico A, Tentler JJ, Pitts TM, Tan AC, Gregory MA, Arcaroli JJ, et al. Rational combination of a MEK inhibitor, selumetinib, and the Wnt/calcium pathway modulator, cyclosporin A, in preclinical models of colorectal cancer. Clin Cancer Res. 2013;19(15):4149–62.
Tan X, Hu L, Luquette LJ, Gao G, Liu Y, Qu H, et al. Systematic identification of synergistic drug pairs targeting HIV. Nat Biotechnol. 2012;30(11):1125–30.
Roller D, Axelrod M, Capaldo B, Jensen K, Mackey A, Michael J. Synthetic lethal screening with small molecule inhibitors provides a pathway to rational combination therapies for melanoma. Mol Cancer Ther. 2012;11(11):2505–15.
Matthews Griner LA, Guha R, Shinn P, Young RM, Keller JM, Liu D, et al. High-throughput combinatorial screening identifies drugs that cooperate with ibrutinib to kill activated B-cell-like diffuse large B-cell lymphoma cells. Proc Natl Acad Sci U S A. 2014;111(6):2349–54.
Holbeck S, Collins JM, Doroshow JH. NCI-60 combination screening matrix of approved anticancer drugs. Eur J Cancer. 2012;48(Suppl:6):11.
Marusyk A, Almendro V, Polyak K. Intra-tumour heterogeneity: a looking glass for cancer? Nat Rev Cancer. 2012;12(5):323–34.
Zhao B, Hemann MT, Lauffenburger DA. Intratumor heterogeneity alters most effective drugs in designed combinations. Proc Natl Acad Sci U S A. 2014;111(29):10773–8.
Zhao B, Pritchard JR, Lauffenburger DA, Hemann MT. Addressing genetic tumor heterogeneity through computationally predictive combination therapy. Cancer Discov. 2014;4(2):166–74.
Pe’er D, Hacohen N. Principles and strategies for developing network models in cancer. Cell. 2011;144(6):864–73.
Sharma SV, Haber DA, Settleman J. Cell line-based platforms to evaluate the therapeutic efficacy of candidate anticancer agents. Nat Rev Cancer. 2010;10(4):241–53.
Tentler JJ, Tan AC, Weekes CD, Jimeno A, Leong S, Pitts TM, et al. Patient-derived tumour xenografts as models for oncology drug development. Nat Rev Clin Oncol. 2012;9(6):338–50.
Lieu CH, Tan AC, Leong S, Diamond JR, Eckhardt SG. From bench to bedside: lessons learned in translating preclinical studies in cancer drug development. J Natl Cancer Inst. 2013;105(19):1441–56.
Iyengar R, Zhao S, Chung SW, Mager DE, Gallo JM. Merging systems biology with pharmacodynamics. Sci Transl Med. 2012;4(126):126ps7.
Sorger PK, Allerheiligen SRB. Quantitative and systems pharmacology in the post-genomic Era: New approaches to discovering drugs and understanding therapeutic mechanisms. An NIH white paper by the QSP workshop group. 2011. p. 0–47.
Feala JD, Cortes J, Duxbury PM, Piermarocchi C, McCulloch AD, Paternostro G. Systems approaches and algorithms for discovery of combinatorial therapies. Wiley Interdiscip Rev Syst Biol Med. 2010;2(2):181–93.
We thank the Tan lab members for useful comments on this manuscript.
This work is partly supported by the National Institutes of Health under Ruth L. Kirschstein National Research Service Award T32CA17468 (K.A.R.), the National Institutes of Health P50CA058187, Cancer League of Colorado, the Department of Defense Award W81XWH-11-1-0527 and the David F. and Margaret T. Grohne Family Foundation. Its contents are solely the responsibility of the authors and do not necessarily represent the official views of the funders.
The authors declare that they have no competing interests.
KAR and ACT wrote the manuscript. Both authors read and approved the final manuscript.
Karen Ryall is a post-doctoral research fellow at the University of Colorado School of Medicine. Her research interests include systems biology, computational modeling, and bioinformatics.
Aik Choon Tan is an Associate Professor of Bioinformatics at the University of Colorado School of Medicine. His research interests include translational bioinformatics and cancer systems biology.