Research ArticlesMEDICINE

Drug Synergy Screen and Network Modeling in Dedifferentiated Liposarcoma Identifies CDK4 and IGF1R as Synergistic Drug Targets

See allHide authors and affiliations

Science Signaling  24 Sep 2013:
Vol. 6, Issue 294, pp. ra85
DOI: 10.1126/scisignal.2004014


Dedifferentiated liposarcoma (DDLS) is a rare but aggressive cancer with high recurrence and low response rates to targeted therapies. Increasing treatment efficacy may require combinations of targeted agents that counteract the effects of multiple abnormalities. To identify a possible multicomponent therapy, we performed a combinatorial drug screen in a DDLS-derived cell line and identified cyclin-dependent kinase 4 (CDK4) and insulin-like growth factor 1 receptor (IGF1R) as synergistic drug targets. We measured the phosphorylation of multiple proteins and cell viability in response to systematic drug combinations and derived computational models of the signaling network. These models predict that the observed synergy in reducing cell viability with CDK4 and IGF1R inhibitors depends on the activity of the AKT pathway. Experiments confirmed that combined inhibition of CDK4 and IGF1R cooperatively suppresses the activation of proteins within the AKT pathway. Consistent with these findings, synergistic reductions in cell viability were also found when combining CDK4 inhibition with inhibition of either AKT or epidermal growth factor receptor (EGFR), another receptor similar to IGF1R that activates AKT. Thus, network models derived from context-specific proteomic measurements of systematically perturbed cancer cells may reveal cancer-specific signaling mechanisms and aid in the design of effective combination therapies.


Liposarcoma is the most common type of soft tissue sarcoma (1). Among the subtypes of liposarcoma, dedifferentiated liposarcoma (DDLS) is associated with the lowest survival rate (2) and often recurs or metastasizes despite treatment with surgery, radiation, or chemotherapy (3). Because response rates to classical chemotherapeutics are low (4), targeted agents have increasingly been under preclinical and clinical investigation for DDLS treatment (3). Unfortunately, drugs directed to kinases such as vascular endothelial growth factor receptor (VEGFR), platelet-derived growth factor receptor (PDGFR) (for example, sorafenib), and BCR-ABL (for example, imatinib) show limited response in phase 2 trials (5, 6). Profiling of DDLS reveals complex DNA copy number changes across the genome as well as recurrent focal alterations, including a high-frequency (~90%) 12q13-15 amplification that harbors the CDK4 and MDM2 oncogenes (7), which respectively encode cyclin-dependent kinase 4 (also known as cell division kinase 4) and mouse double minute 2 homolog, an E3 ubiquitin ligase. We previously evaluated the selective oral inhibitor of CDK4 and CDK6 (hereinafter referred to simply as CDK4), PD0332991, in a phase 2 clinical trial in patients with well-differentiated liposarcoma (WDLS) and DDLS and found prolongation of progression-free survival (8). However, the response rate remains low, suggesting that PD0332991 needs to be combined with other drugs to enhance its antitumor efficacy, exemplified by the advances in treatment of estrogen receptor–positive and human epidermal growth factor receptor 2–negative (ER+/HER2) advanced breast cancer with PD0332991 and the aromatase inhibitor letrozole (clinical trial number NCT01740427).

Targeted therapies are revolutionizing cancer treatment by acting on patient-specific genetic alterations with fewer side effects than conventional cytotoxic chemotherapy. Despite these advantages, single-target therapies for cancer often have limited clinical success because of resistance, such as with HER2-targeted therapies for breast cancer (9, 10). In the case of limited initial treatment response (referred to as primary resistance), combination therapy with two or more targeted drugs improves efficacy (11, 12), such as the combination of the HER2-specific antibody trastuzumab with the receptor tyrosine kinase inhibitor lapatinib in HER2-positive metastatic breast cancer (13). Ongoing clinical trials with combination therapies include BRAF and MEK [mitogen-activated protein kinase (MAPK) kinase] inhibition for BRAF mutant melanoma (clinical trial number NCT01072175) (14), phosphatidylinositol 3-kinase (PI3K) and MEK inhibition for PIK3CA/KRAS mutant colorectal cancer (clinical trial number NCT00996892), and AKT and MEK inhibition for advanced solid tumors (clinical trial number NCT01021748). The improvements of combination therapy are most likely because of cooperative inhibition of the multiple cellular signaling events typically altered in cancer (1517). In DDLS, for example, in addition to the amplicon with CDK4 and MDM2, there are complex genetic rearrangements, including partial deletion of chromosomes 11q and 19q, suggesting that several dysregulated pathways are likely involved in DDLS pathogenesis (18, 19). Targeting multiple pathways with combinations of drugs therefore represents a relevant strategy for DDLS treatment and could potentially lead to higher response rates and better clinical outcomes. Moreover, by targeting two different pathways converging on a common phenotype, such as cell proliferation, it is possible to obtain a more-than-additive (synergistic) response compared to the individual agents alone (15, 20), sometimes referred to as the parallel pathway inhibition model (17, 21, 22). Synergistic drugs have the potential to not only increase the rate of initial treatment response (16, 23) but also reduce the concentration of each needed to elicit a given effect and consequently improve the therapeutic index (24).

The issue of how best to optimize drug combinations becomes particularly acute in developing combinations of drugs that have distinct cell cycle effects. In a process termed cell cycle–mediated drug resistance, a cell cycle inhibitor that blocks cell growth in one phase of the cell cycle (such as G1) can then antagonize a second drug that exerts its cytotoxic effect within another phase of the cell cycle (such as mitosis) (25). For example, the pan-CDK inhibitor flavopiridol, which induces a G1 cell cycle arrest, essentially prevents mitotic spindle poisons such as paclitaxel or docetaxel from exerting their antitumor effects during mitosis (26, 27). A similar inability to enhance the effect of cell cycle–specific chemotherapies is shown for the CDK4 inhibitor PD0332991 (28). Thus, combinations of cell cycle inhibitors with conventional chemotherapy or with small-molecule inhibitors that could similarly affect the cell cycle are not intuitive.

With the goal of identifying effective, synergistic drug combinations for DDLS treatment, we performed a synergy screen with 14 targeted drugs in a cell line derived from a DDLS patient with the 12q13-15 amplicon. Of several identified synergistic drug pairs, we focused on an unexpected synergy observed with paired CDK4 and IGF1R (insulin-like growth factor 1 receptor) inhibition. To investigate the possible mechanisms underlying this synergy, we applied our recently developed hybrid experimental and computational approach for deriving context-specific signaling models (29, 30). In a three-step process, we (i) performed combinatorial drug perturbations and used high-throughput reverse-phase protein arrays (RPPAs) to measure the response of a panel of 13 proteins and phosphoproteins, (ii) derived network models based on the perturbation profiles and prior knowledge of signaling pathways, and (iii) quantitatively simulated the effects of perturbations on signaling events through our models. The network model was used to predict the key regulatory proteins mediating the synergy, which were then assessed experimentally. Although our cell line–specific models may not be generalizable to other contexts, our results illustrate the potential for computationally assisted design and analysis of systematic perturbation screens to efficiently explore therapeutically relevant drug combinations. We anticipate that such integrated approaches to combinatorial therapeutics may provide translational opportunities for further development of DDLS treatment.


Drug combination screen identifies synergistic drug targets in DDLS

To investigate synergistic drug interactions in DDLS, we selected a tumor-derived cell line, DDLS8817, which contains several genomic alterations characteristic of DDLS, including amplification of 12q13-15 and partial loss of segments in chromosomes 3p and 19q (19) (fig. S1). We then chose 14 small-molecule drugs (table S1 and Fig. 1) that target proteins within the canonical mTOR (mammalian target of rapamycin), PI3K/AKT, and MAPK pathways, which are altered across a wide range of cancers (31). In addition, some of these drugs are under preclinical and clinical investigation for treatment of various sarcoma subtypes, including drugs targeting CDK4 (8), PDGFR (platelet-derived growth factor receptor) (5), and IGF1R (32), although only partial therapeutic responses have been observed.

Fig. 1 Drug combination screen identifies synergistic and antagonistic drug targets in the dedifferentiated liposarcoma cell line DDLS8817.

(A) Example of CI score calculation of paired-drug perturbation with the MEK inhibitor SL327 (MEKi) and the AKT inhibitor AKT inhibitor VIII (AKTi). Cell viability was estimated after 72 hours of drug treatment by the resazurin assay, which measures cellular metabolic activity. Error bars represent SD of at least four biological replicates. (B) CI scores for a drug synergy screen performed in DDLS8817 cells using 14 targeted inhibitors (“Inh” or “i” in the labels). CI scores were derived as described in (A) and categorized as synergistic (<0.75, red), antagonistic (>1.5, blue), or additive (green). In some cases, CI scores could not be calculated (gray). Inset shows the distribution of CI scores. A complete list of targets and secondary targets of the drugs used is provided in table S1. Each CI score represents data from seven different drug doses of single- and paired-drug treatments with at least four biological replicates per condition.

To determine potential nonadditive effects (indicating the existence of common mechanisms, known as epistasis), the activity of combined agents is typically compared to single-agent activities and related to a null expectation model that assumes no interaction between the drugs (24). The most commonly used null expectation models are Loewe additivity and Bliss independence (15), whereby effects can be categorized as additive, synergistic, or antagonistic. For this study, we measured synergy by the Loewe additivity–based combination index (CI) score, which can handle cases where two drugs act on targets regulating a common pathway (33) (Fig. 1A). We then performed a systematic dose-response screen of single and combined agents at seven different concentrations administered for 72 hours and estimated the effect on cell viability, using metabolic activity–based resazurin assay. Using more than 10,000 cell viability measurements of the effect of all single- and dual-drug perturbations (table S2), we calculated CI scores for each drug combination on the basis of half-maximal effective concentration (EC50) values obtained from sigmoid-fitted dose-response curves (Fig. 1B and fig. S2). CI scores less than 0.75 were considered synergistic, scores larger than 1.5 were considered antagonistic, and the rest were considered additive (Fig. 1B and table S3). In some cases, CI scores could not be calculated; for example, the mTOR inhibitor rapamycin arrested cells in a dose-independent manner. Using these CI score cutoffs, we identified nine synergistic drug-drug interactions (of a possible 91), which corresponds with previous reports in which 4 to 10% of the combinations were synergistic (3437). The synergistic pairs included combined inhibition of epidermal growth factor receptor (EGFR) and IGF1R, MEK and either PI3K or AKT, among others [inclusively, AKT and MEK, ERK (extracellular signal–regulated kinase) and HDAC (histone deacetylase), ERK and MET, IGF1R and CDK4, IGF1R and EGFR, IGF1R and STAT3 (signal transducer and activator of transcription 3), MEK and MET, PDGFR and MEK, and PI3K and MEK].

CDK4 and IGF1R are synergistic drug targets in DDLS

We decided to further investigate the particular combination of CDK4 and IGF1R inhibition using the small molecules Ryuvidine and AG538 because of the recurrent CDK4 amplicon in DDLS and the clinical relevance of CDK4 inhibition for treatment of liposarcoma. In addition, the IGF1R target was a surprising combination candidate with CDK4 inhibition because the two molecules are generally presumed to control cell survival through separate pathways [the AKT/mTOR pathway by IGF1R, and the retinoblastoma (RB) pathway by CDK4].

We first tested whether the main drug targets were driving the identified synergy and whether it was exclusive to the DDLS8817 cell line. To investigate this, we used alternative inhibitors, PD0332991 (a CDK4 inhibitor) and R1507 (an IGF1R antibody), in two DDLS cell lines, DDLS8817 and LPS141 (fig. S1 and table S4). Mixing serial dilutions of R1507 and PD0332991 in all combinations (dose matrix) and measuring their effect on cell metabolic activity (as an indication of cell viability) after 6 days of drug treatment, we calculated the CI scores using CompuSyn (38), and found that these drugs were synergistic in both cell lines (average CI scores were 0.34 ± 0.19 for DDLS8817 and 0.63 ± 0.17 for LPS141) (Fig. 2). In addition, the IGF1R small-molecule inhibitor NVP-AEW541 appeared to be synergistic in combination with PD0332991, although this was more prominent for LPS141 than for DDLS8817 (fig. S3), most likely because of a higher baseline abundance of IGF1R in LPS141. Finally, in LPS141, PD0332991 but not IGF1R inhibition alone (neither by inhibitor NVPAE nor by antibody R1507) induced G1 cell cycle arrest, with greater effects from combined treatment (fig. S4). These results indicate that the synergy of CDK4 and IGF1R inhibitors is most likely mediated by inhibition of the main targets rather than off-target effects, and support the findings of our initial screen in an independent cell line and with different agents targeted to the same molecules.

Fig. 2 IGF1R and CDK4 are synergistic drug targets in DDLS.

(A) Dose-response effects of the CDK4 inhibitor PD0332991 and the IGF1R antibody R1507 on inhibition of cell viability in DDLS8817 cells estimated after 6 days of drug treatment with the cell counting kit-8 (CCK-8) assay, which measures cell metabolic activity (top). Each condition represents at least four biological replicates. CI scores were calculated at various effect levels (bottom). (B) Performed as in (A) in the DDLS cell line LPS141.

RPPA assay provides proteomic profiles for DDLS-specific network modeling

We next investigated possible mechanisms underlying the observed CDK4 and IGF1R inhibitor synergy by building network models of signaling pathways in DDLS. The approach consisted of profiling the cellular response (proteomic and phenotypic measurements) to a series of drug-pair perturbations, followed by computational transformation of the data into network models that quantitatively link proteins in signaling network circuits (29, 30).

For the first profiling step, we used the same cell line and the same set of drugs as in our initial drug synergy screen. In an effort to reveal interactions at relatively low doses, we selected drug concentrations that inhibited phosphorylation of the presumed target by 40% (IC40) at 24 hours (figs. S5 and S6 and table S1). Notably, the protein IC40 values used here were roughly one order of magnitude lower than the drugs’ IC40 values for cell viability.

We then treated DDLS8817 cells with all individual drugs and all possible pairs of drugs, resulting in a total of 105 different perturbations, not counting biological replicates and untreated controls (Fig. 3A). Biological replicates of cellular lysates were spotted in fivefold dilution series using the RPPA platform, and the abundance of more than 100 different proteins and phosphoproteins was assessed with a panel of antibodies (39). A final set of 13 antibodies was chosen for readouts based on quality control of slides, correlation between biological replicates, response to drug perturbation, and evidence of the protein being relevant for liposarcoma tumor biology (Fig. 3B). This final set of antibody readouts shows the overall pattern of increased (green) and decreased (blue) phosphoprotein abundance over the range of single- and dual-drug perturbations. Several consistent patterns emerge across similar drug conditions (horizontal “bars”), such as repression of phosphorylated AKT (at Ser473) in the majority of drug pairs with the AKT inhibitor (conditions 15 to 27) as well as repression of phosphorylated ERK (at Thr202/Tyr204) with pairs of drugs that include ERK inhibition (conditions 28 to 39) or EGFR inhibition (conditions 40 to 50). Although difficult to interpret without computational analyses, distinct vertical patterns also appear, such as repression of multiple phosphoproteins with combined MEK and PI3K inhibition (condition 74)—a drug combination that was also found to have a synergistic reduction in cell viability (Fig. 1B). We tested the selected set of RPPA antibodies using Western blotting and found no apparent cross-reactivity to nontarget antigens in the cell lines used (fig. S7). In addition, we profiled the two cell lines for mRNA expression using microarrays and found that the expression of the network genes ranked in the top half of all genes (fig. S8). Finally, we found multiple DNA copy number gains in the network genes, and we detected no homozygous deletions (fig. S9). Together, these observations on transcript, copy number, and protein abundance indicate that all selected genes were expressed at “steady state” in unperturbed conditions.

Fig. 3 Systematic drug treatments and large-scale proteomic profiling in DDLS8817.

(A) Experimental design of 14 individual and 91 pairwise drug perturbations with a set of targeted small-molecule drugs. (B) Response of 13 proteins or phosphoproteins [phosphorylated at the indicated residue(s)] after 24 hours of drug treatment as shown in (A), assessed by RPPA. The effect on cell viability was estimated with the resazurin assay (measures cell metabolic activity) after 72 hours of drug treatment. All readouts were z score–normalized. Each condition represents the mean of two biological replicates.

Proteomic data and prior knowledge interactions provide quantitative models of signaling pathways in DDLS

For inferring DDLS-specific network models, we used existing information on protein interactions from several databases (4042) and information on signaling pathways involved in sarcoma (18). With this information, we constructed a prior knowledge network of interactions (termed edges) in the selected set of proteins measured by RPPA (termed nodes) (Fig. 4A and table S5). For simplicity, we represented the combined effects of cell proliferation, survival, and death as a single node (cell viability).

Fig. 4 Network inference of proteomic data profiles and prior knowledge interactions provides signaling models specific to DDLS.

(A) A network of prior knowledge interactions (edges) in the selected set of protein readouts (nodes) measured by RPPA and their connection to cell viability, as measured by metabolic activity assays in Fig. 3B. (B) Inferred network models from perturbation-response data and prior knowledge information, with the linewidth reflecting the most probable interactions. The network represents the average network of the 100 lowest-error models. Predicted interactions (gray) and prior knowledge interactions that are retained (blue) and rejected are indicated. Nodes that are perturbed but not observed are termed “activity node” and preceded by an “a” (for example, aAKT) and represent the presumed, direct target of the drugs applied.

We next derived quantitative network models using three types of data: the perturbation-response profiles (Fig. 3B), the drug-target relationships (table S1), and the list of prior knowledge interactions (Fig. 4A). To maintain flexibility to infer interactions that are consistent with the experimental data, we softly enforced the prior knowledge interactions, such that inconsistent interactions could be rejected. We inferred interactions into a given “target” node one at a time by searching the data for weighted subsets of possible “source” nodes that collectively covary with the target node. Technically, this was done by minimizing a bipartite cost function that penalizes (i) discrepancies between model-predicted and measured values and (ii) a large number of interactions. In this way, edges are numerical parameters representing the relative influence of one node on the state of another, representing logical but not necessarily direct biological interactions. Explicitly searching all possible subsets of source nodes is computationally prohibitive because of combinatorial explosion of possibilities. Therefore, we used a statistical inference algorithm called Belief Propagation (BP), which is designed to find the most probable interactions in large systems (30) (described in fig. S10). The output of the BP method is a probability distribution for each of the N2 possible parameters, which we used to collect an ensemble of network models with the most probable interactions.

To visualize the inferred ensemble of network models, we plotted all of the highly probable interactions, with the edge width representing the inferred interaction strength (Fig. 4B). Qualitatively, we observed a fairly complex network of interactions with several phosphoproteins acting as hubs [for example, AKT at pSer473, p70S6 kinase (p70S6K) at pThr389, and glycogen synthase kinase α and β (GSK3α/β) at pSer21] and without any obvious modular or isolated pathways. The majority of the prior knowledge interactions (26 of 34) were consistent with the data (fig. S11) and were retained in the inferred network. Supporting this, models based on prior knowledge alone performed better than random models as estimated by the mean squared error between measured and model estimated data points (fig. S12). The rejected interactions were mostly those connecting the drugs (represented as activity nodes; for example, aPDGFR) to their presumed targets. Although some of this discrepancy may result from use of low drug concentrations, in some cases the measured phosphorylation site may not be affected by drug treatment despite the target being inhibited (for example, the abundance of phosphorylated MEK may not be affected by a MEK inhibitor). Although drug specificity is crucial for accurate network modeling, we found that off-target effects could be tolerated (fig. S13). Thus, using our network inference approach, we developed context-specific DDLS network models that can be used to simulate effects of perturbations.

Network models are predictive of cellular responses to drugs

To use the models for prediction, we generated a collection of quantitative models by sampling the probability distributions for each parameter and subsequently optimizing with a gradient descent method. Of a collection of 1000 distinct network models, we then used the 100 best-performing models (lowest error) to estimate how well the models quantitatively link nodes (measured proteins) in the network to the cell viability node. To assess this in our models, we simulated an inhibitory perturbation (single or paired) by applying a negative force on the target node(s), which in turn propagated changes to all other linked nodes sequentially, thus simulating a cascade of signaling events. To avoid overfitting, we performed this test using a cross-validation procedure in which we generated network models, where all measurements of paired-drug perturbations involving the drug in question were left out (using “leave-k-out” cross-validation). We found that the network modeling approach was reasonably accurate in predicting the effect of any pairwise node perturbation on cell viability (fig. S14). This indicated that central signaling connectivities were captured in a biologically relevant manner because the predicted cell viability outcome depends on how information is propagated in the network models.

Network models capture synergy and recapitulate the CDK4 and IGF1R inhibitor synergism

We next investigated if the network models were able to identify epistatic (synergistic and antagonistic) interactions in the data. Applying a similar leave-k-out approach, we determined possible epistatic effects on cell viability by inhibiting two nodes at various strengths in all permutations. We determined epistatic effects in the network models, synergy (S), as the difference between the effect of the paired simulated perturbations (Zsim) and the Loewe additivity surface (ZLoewe): S=(ZsimZLoewe), where ZLoewe was derived from the added effects of the single-node perturbations (43). For example, when applying this approach to investigate the effects of combined inhibition of the ERK activity (aERK) and 4EBP1 (eukaryotic translation initiation factor 4E–binding protein 1)–pSer65 nodes, we find that cell viability is inhibited more than expected over the additive independence model (Fig. 5, A and B). Although it was not experimentally tested, this epistatic interaction was derived nontrivially from the input data and showed an ability to model and predict new potential synergistic drug targets. We then calculated S for all possible node pairs and found that our network models were able to capture epistasis as numerous synergistic and antagonistic effects were identified (Fig. 5C). Although S is unitless and can be used only for comparisons within the same data set, the overall trend of the distribution is similar to that found experimentally (Fig. 1B). This indicates that the modeling approach is not over- or underestimating synergistic effects; however, further efforts are needed to confirm this observation. We then compared predicted versus experimentally tested drug combinations. Although there were several miscategorized predictions, many of the drug combinations with the strongest predicted synergistic and antagonistic effects were categorized in accordance with the experiments (fig. S15).

Fig. 5 Synergistic effects are captured in the network models, and the experimentally observed CDK4-IGF1R drug synergy is recapitulated.

(A) Example of calculation of model-based synergy scores (S) by in silico perturbation of 4EBP1-pSer65 (phosphorylated at serine 65) and ERK encoded as an activity node (aERK, external drug node). These nodes were inhibited with eight different perturbation strengths (u) in all possible combinations, and the effect was recorded as the response on the cell viability node (z score). (B) Nonadditive synergy effects (synergy over linear drug combination model) were determined as the difference between the effect of the paired inhibition and the added effects of the two single-node inhibitions (S = −0.24 for aERK and 4EBP1-pSer65). (C) Computed synergy scores for all node pairs where each synergy score was determined from 64 unique perturbations as in (A). Synergy scores were categorized into three classes: S < −0.20 was considered synergistic (red), >0.20 antagonistic (blue), and otherwise additive (green). The synergy score for IGF1R and CDK4 node inhibition is highlighted (S = −0.23).

Similar to the experimental results, the network models predicted a synergistic effect from combined inhibition of IGF1R and CDK4 (Fig. 5C and fig. S16). Although the predicted efficacy, unlike the experimental synergy, was relatively modest, it was significantly different from synergy scores determined from perturbing the same nodes in random but similarly parameterized models (fig. S16C). Thus, the models effectively recapitulated the phenotypic observations for combined CDK4 and IGF1R inhibition.

Models predict that the AKT pathway is involved in mediating CDK4-IGF1R drug synergy

In an effort to identify important network features mediating the synergy of IGF1R and CDK4 inhibition, we systematically altered the network models to find the interactions that were essential for the synergistic effect on the cell viability node. Each interaction observed in at least one of the top 100 models was removed in turn, then the synergy calculation was repeated. The interactions whose removal resulted in the most pronounced drop in the overall synergy score were ranked (Fig. 6). This analysis yielded several important insights. First, most of the model interactions did not contribute to the synergy. Second, in agreement with the parallel pathway hypothesis model for drug synergy (17, 21, 22), the most important interactions formed two seemingly parallel and nonoverlapping pathways controlling cell viability. The most essential interactions were those linking the drug targets (aCDK4 and aIGF1R) to cell viability either directly or through downstream effectors such as RB-pSer780 (for aCDK4) and S6-pSer235 (for aIGF1R). Although some of these interactions are likely indirect and did not occur in the average network, they do represent strong causal connection between the drug target and the predicted control over cell viability. In addition, this analysis suggested that control of the phosphorylation of AKT at Ser473 by activated EGFR (phosphorylated at Tyr992) also plays an important role in mediating the synergy. The connections between AKT-pSer473, p70S6K-pThr389, S6-pSer235, and cell viability represent the canonical AKT/mTOR pathway cascades regulating cell proliferation (44, 45) and are consistent with IGF1R-mediated cell proliferation regulation through this pathway. These predictions suggested that AKT was involved in mediating part of the observed synergy. Therefore, we decided to investigate this in further detail.

Fig. 6 Simulation of information flow in network models predicts several important interactions mediating the synergy of CDK4 and IGF1R inhibition, including AKT pathway members.

Network edges ranked by their contribution to the model-simulated synergy between CDK4 and IGF1R inhibitors. Each edge was removed in turn, and the effect on the cell viability synergy score was recalculated and expressed as the percent suppression of original synergy score. The leave-edge-out calculations were performed with the 100 lowest-error models.

AKT pathway–mediated control over cell viability is likely involved in the CDK4 and IGF1R inhibitor synergy

To begin interrogating the molecular pathways that might explain the CDK4 and IGF1R inhibitor synergy, we investigated the effects of their inhibition on IGF1R, PI3K/AKT, and mTOR signaling in the cell lines DDLS8817 and LPS141. To address specificity issues, we used two IGF1R inhibitors (the IGF1R antibody R1507 and the small-molecule IGF1R inhibitor NVP-AEW541) and the CDK4 inhibitor PD0332991 as well as CDK4-targeted small interfering RNA (siRNA). Under all of the conditions tested, the IGF1R inhibitor and antibody decreased the abundance of phosphorylated AKT (Fig. 7, A and B). With R1507, this decrease in phosphorylated AKT was coincident with decreased IGF1R abundance, more markedly in LPS141 cells, which have a high baseline abundance of IGF1R. PD0332991 or CDK4 siRNA alone resulted in variable suppression of mTOR signaling with a decrease in both phosphorylated S6 (at Ser235–236) and pS6K (at Thr389) under some of the conditions tested. However, with the combination treatment of PD0332991 and R1507, there was evidence of enhanced inhibition of mTOR signaling indicated by decreased S6 phosphorylation (pSer235–236) and S6K phosphorylation (pThr389) (fig. S17, B and C), as well as a similar but nonsignificant trend for suppression of phosphorylated AKT (fig. S17A). These data suggest that dual blockade of CDK4 and IGF1R signaling results in cooperativity such that two prosurvival pathways are inhibited. Despite the MAPK pathway being downstream of IGF1R, the same drug combination did not appear to cooperatively suppress the phosphorylation of ERK (Fig. 7A).

Fig. 7 The AKT pathway is likely involved in the synergy of CDK4 and IGF1R inhibitors.

(A) Western blot of DDLS8817 and LPS141 cells inhibited for 12 hours with R1507 (10 μg/ml) (IGF1R antibody), 1 μM PD0332991 (CDK4 inhibitor), and siRNA-mediated knockdown of CDK4. (B) Similar to (A), except cells were inhibited for 24 hours with 1 μM NVP-AEW541 (IGF1R inhibitor) and PD0332991. Western blots shown are representative data of at least two independent experiments. (C) Dose-response measurements of cell metabolic activity with the CCK-8 assay (correlates with cell viability) of DDLS8817 cells after drug treatment for 6 days with the AKT inhibitor MK2206 and the CDK4 inhibitor PD0332991. The CI score is indicated and was determined at EC50 levels indicated by dashed lines. Error bars represent SD of six biological replicates. (D) Similar to (C) but in LPS141 cells.

On the basis of these observations and our modeling results, we hypothesized that combining PD0332991 with inhibitors of the AKT pathway would also be synergistic. Indeed, measuring cell viability after 6 days of drug treatment, we found that combined AKT and CDK4 inhibition with MK2206 and PD0332991 resulted in a synergistic effect in both DDLS8817 and LPS141 (Fig. 7, C and D). Consistent with this, inhibiting EGFR, which is upstream of AKT, also synergized with PD0332991 (fig. S18A). In addition, in a triple-drug perturbation, we observed cooperative effects of combining CDK4 inhibitor, IGF1R inhibitor, and EGFR inhibitor (fig. S18, B and C). In contrast, neither the MEK inhibitor AZD6244 nor the ERK inhibitor FR180204 synergized with PD0332991 because the effects of these dual-drug perturbations were comparable to single PD0332991 treatment (fig. S19), arguing that the MAPK pathway is not a point of convergence in the CDK4 and IGF1R inhibitor synergy. Together, these observations indicate that part of the synergistic mechanism may result from a more-than-additive suppression of phosphorylated AKT (at Ser473) and that the added antiproliferative effect may involve inhibition of the AKT pathway rather than the MAPK pathway.


Here, we used a perturbation-based systems biology approach to analyze drug combination effects in DDLS. Performing a drug synergy screen in a patient-derived cell line, we identified CDK4 and IGF1R as synergistic drug targets. To investigate potential mechanisms of this synergy, we applied an integrated experimental-computational approach to infer network models from rich perturbation-response profiles. We used these network models to quantitatively describe signaling pathways in DDLS and model the observed CDK4-IGF1R drug synergy. Both predictions and experiments suggest that AKT or effectors of AKT are involved in mediating this more-than-additive effect. Although other mechanisms are likely involved and further methodological developments are needed, these results reveal the power of network pharmacology approaches for identifying and modeling drug synergy.

Several features of this work may contribute substantially to the discovery and analysis of effective combination therapies. We have extended the conventional drug synergy screening approach by combining cell viability measurements with high-throughput proteomic measurements of systematically perturbed cancer cells. This enables us to apply our powerful probability-based BP algorithm to derive network models that integrate both prior knowledge of pathways and direct measurements of signaling events in the system of interest. In this way, we infer network models of signaling connectivities in a context-specific manner. This integrated experimental-computational approach is predictive of cell viability outcome to network perturbations and allows for modeling and predicting mechanisms of synergistic drug combinations, which is typically not possible with other approaches. Although the generalizability of our models remains untested, our approach is generalizable to other contexts, and we are currently modeling oncogenic signaling in several additional cancer types.

The applicability of our approach is exemplified by the discovery and modeling of the CDK4 and IGF1R drug synergy. Increasing evidence points to CDK4 as a major oncogenic driver in DDLS (7, 18, 19); however, because IGF1R is not frequently altered in DDLS (7), this drug combination would not likely be identified with standard genomic screening methods. On the basis of these results, we are continuing preclinical investigations in DDLS with paired CDK4 and IGF1R inhibition as well as devising clinical trials with multicomponent therapy consisting of PD0332991 combined with either IGF1R or EGFR inhibition. Both experiments and model simulations indicate that CDK4 and IGF1R inhibitors operate through separate survival pathways (through RB and AKT/mTOR, respectively), supporting the parallel pathway hypothesis for synergism (17, 21, 22). In this light, we expect that the combination therapies suggested here may improve response rates in the initial treatment phase and thereby reduce primary resistance. Moreover, because CDK4 inhibition, but not IGF1R blockade, arrests cells in G1 (fig. S4), this particular combination may avoid problems of cell cycle–mediated drug resistance (25).

The network models predicted that regulation of AKT’s effect on cell viability is critical for the observed CDK4 and IGF1R drug synergy. The AKT pathway is one of the major downstream pathways of IGF1R (46), and indeed, we found that IGF1R inhibition with R1507 or NVP-AEW541 represses phosphorylation of AKT and downstream members of the AKT/mTOR pathway (p70S6K, S6). Consistent with the literature, CDK4 inhibition repressed phosphorylated RB (Fig. 7A) and had little to no repressive effect on the AKT/mTOR pathway by itself. However, the mTOR pathway was cooperatively repressed by paired inhibition of CDK4 and IGF1R compared to single treatment alone, and a similar trend was observed for the AKT pathway (fig. S17). Whereas AKT is an upstream regulator of cyclin D1 and CDK4 activity (47), CDK4 does not, to our knowledge, directly control AKT pathway activity. In this light, it is surprising that the paired drug perturbation has a more-than-additive suppressive effect on this pathway. This suggests the existence of an unknown connectivity or feedback mechanism, although further research is required to evaluate this idea. In accordance with the premise that the control of the AKT pathway contributes to the CDK4 and IGF1R inhibitor synergy, we found that combined inhibition of CDK4 and AKT or its upstream regulator EGFR resulted in a synergistic repression of cell proliferation. In contrast, combining CDK4 inhibition with either MEK inhibition or ERK inhibition showed no cooperative effects, eliminating the MAPK pathway by itself as a mechanism driving the synergy.

Designing the experimental setup for large-scale profiling, as performed here, requires a range of choices with respect to cell type, growth conditions, drug selection, drug concentrations, assays, and so on. We chose clinically relevant cell lines and drugs for this study. However, some of the drugs used in the screen were not identical to those under preclinical development. For example, we used the CDK4 inhibitor Ryuvidine in the synergy screen instead of the more specific and clinically relevant PD0332991. We found greater synergy with PD0332991 and R1507 than with Ryuvidine and AG538, possibly reflecting differences in specificities. Similarly, the CDK4-AKT synergy and the CDK4-EGFR synergy were not observed in the original screen in which the CDK4 inhibitor was Ryuvidine (Fig. 1B), but were observed in the follow-up experiments using PD0332991 (Fig. 7, C and D). Differences in the experimental setup, including different time points (3 versus 6 days) and different drug specificities, may explain some of this discrepancy. Another choice in the experimental design was in the determination of drug concentrations for the RPPA-based protein profiling, which was done on the basis of changes in the abundance of drug target proteins detected by Western blot. Ideally, the RPPA assay itself should be used for this purpose because targets of certain drugs elicited only moderate responses, possibly because of sensitivity differences between RPPA and Western blot.

Increasingly, there is a demand for developing data-driven network inference algorithms that link signaling events to phenotypic outcomes (17, 48). Our modeling approach falls in between fully parameterized mass action kinetic models (49, 50) and “black box” machine learning models (51), offering significant predictive power while maintaining biological interpretability (30). We use probability distributions for each possible model interaction, calculated from an iterative BP algorithm, and thereby identify likely interactions and parameterize hundreds of predictive models. The ability to quickly construct these models frees us from relying on a single model and enables us to attach probabilities to predicted outcomes and construct alternative hypotheses. For example, multiple models (20 of the 100 lowest-error models) predicted strong synergy when inhibiting the CDK4 and IGF1R nodes (fig. S20). Although these 20 models may differ in some predicted outcomes, they all predicted that tight control of AKT from IGF1R inhibition and other upstream regulators is essential for the observed synergism.

For deriving DDLS network models, we incorporated existing background knowledge of protein interactions because it (i) reduces the search space of all possible interactions, (ii) enables capturing cascades of signaling events with finer granularity, and (iii) facilitates assigning directionality between correlated events. However, heavy restrictions from prior knowledge come with disadvantages. For example, common protein interactions described in other contexts may not exist in the system of interest. To balance between the use of prior knowledge and the ability to infer new interactions, we use flexible constraints so that our method can accept or reject each prior knowledge interaction based on its fit with the experimental data and the inferred parameters. For example, an edge between phosphorylated RB (pSer807-RB) and cell viability was encoded as a prior knowledge interaction but was not present in the average of the 1000 network models. The lack of identification of this interaction could possibly result from technical issues with the RPPA assay, low drug concentration, or off-target effects of Ryuvidine. However, this edge was present in 49 of the top 100 models and thus narrowly missed being classified as an average interaction. This illustrates the importance of using an ensemble of models and not a single averaged model. Nevertheless, in our follow-up experiments, we recapitulated the CDK4-RB interaction because RB (phosphorylated at Ser780 and Ser807) was suppressed by PD0332991 and by CDK4-directed siRNA (Fig. 7A). We found that the selected set of prior knowledge interactions enhanced model performance (fig. S12), but at the same time, the models were not overfit to prior information (fig. S11). We realize that the final selection of prior knowledge interactions is qualitative, and additional known interactions could have been included for better performance. However, in our case, the performance of prior information alone was small compared to overall performance of BP-derived models trained on both data and prior knowledge (fig. S12 and table S6).

Drug specificity is a fundamental element of our modeling approach, and some of the drugs used in this work may have off-target effects. Whereas known off-target effects can directly be incorporated into the modeling approach, unknown off-targets may pose challenges in interpretation and general predictive power of the network models. We found that lack of drug specificity does not prohibitively confuse the inference of interactions because the BP method can infer interactions from the correlations in the data itself (fig. S13). Hereby, our modeling strategy can effectively predict on-target, off-target, and/or indirect drug effects because each drug is represented in the models as an activity node. Specificity is primarily a problem for predicting outcome of response to perturbation because this requires knowing which elements are being perturbed for any drug. Nevertheless, we were able to correctly predict the CDK4-IGF1R synergy outcome using models that were derived from perturbation data that included the CDK4 inhibitor Ryuvidine, which is known to have nonspecific effects. In the future, we expect that by looking at response to different drugs acting on the same primary target, we can disambiguate signals that are due to changes in the known target from those that are due to off-target effects.

Our network modeling approach recapitulated the experimentally observed CDK4-IGF1R drug synergy as well as several of the strongest predicted nonadditive combinations (fig. S15). However, several predictions were miscategorized. With the current modeling approach, we did not expect to accurately predict synergy. The major problem is that only cell viability measurements from one drug concentration (the one used in the RPPA screen) could be incorporated during model training, making it challenging to predict synergy from data obtained from two single-drug and one drug-pair perturbations. On the experimental side, seven different concentrations for both single and paired drugs were used. Furthermore, the experimental synergies were determined with high drug concentrations, whereas the models were trained to data collected from low drug concentrations (roughly 10-fold lower). These discrepancies were limitations in our experimental design and model implementation, which we expect to be resolved in future work. Other limiting factors for the miscategorizations include (i) different measurement time point between the RPPA and the cell viability screen (24 versus 72 hours), (ii) important regulators that were not measured and therefore absent from the derived networks, and/or (iii) possible technical limitations of the experimental or computational approaches. Future work will aim to address this by incorporating multiple doses and time points. With further experimental and algorithmic improvements, the number of measured proteins in the networks can, in principle, be expanded to handle hundreds of nodes. Furthermore, the approach can be extended to predict synergy and mechanisms beyond the set of tested drugs. This is possible because all measured proteins, whether they are perturbed or not, become part of the network models. In this way, the set of possible predictions is further expanded to include any measured protein, enabling efficient computational screening and identification of combinatorial interventions.

This study was an interdisciplinary study to find clinically relevant drug combinations for treatment of DDLS, explore the use of computational modeling to predict synergistic or efficacious drug combinations, and identify biological mechanisms driving synergistic drug combinations. As our data indicate, our integrated approach has contributed with several advancements in these areas. For future studies, we are planning to further test and optimize the CDK4 and IGF1R drug combination as a potential DDLS treatment while simultaneously developing the data-driven modeling strategy for predicting new drug targets and drug combinations.

Materials and Methods

Cell lines and drugs

The LPS141 and DDLS8817 cell lines were derived from high-grade retroperitoneal DDLS tumors (52). All cells were maintained in Dulbecco’s modified Eagle’s medium supplemented with antibiotics and 10% fetal bovine serum. A set of 14 small-molecule drugs was selected for the drug synergy screen and RPPA proteomic assay with DDLS8817 cells (table S1). For follow-up experiments with DDLS8817 and LPS141 cells, an additional seven drugs were used (table S4).

DNA copy number and mRNA expression profiling

RNA and DNA from asynchronously growing DDLS8817 and LPS141 cultures were isolated with RNeasy and DNeasy kits per the manufacturer’s specifications (Qiagen). DNA copy number profiling was performed with the 244K Agilent array-comparative genomic hybridization (aCGH) platform, and standard circular binary segmentation (R/Bioconductor DNAcopy library) was analyzed with the RAE algorithm (53). Transcript profiling was performed with the Illumina HT 12 V3 microarray platform, and data were processed with BeadStudio software version 3.3.7.

Resazurin cell viability assay

Although strictly a measure of cell metabolic activity, the resazurin assay is widely used as a measure of cell viability. DDLS8817 cells were seeded at 100-μl volume per well in 96-well plates (1000 cells per well) and grown for 24 hours. Cells were then inhibited with seven different concentrations (twofold dilution) of single and paired drugs in six replicates by adding 100 μl of drug solution per well. After 72 hours, resazurin (Sigma-Aldrich) was added to each well at a final concentration of 44 μM and incubated for 2 to 3 hours at 37°C. The fluorescent signals were measured with a SpectraMax microplate reader (Molecular Devices Corp.) with 530-nm excitation wavelength and 590-nm emission wavelength. Cell viability was normalized to control cells treated with drug vehicle [dimethyl sulfoxide (DMSO)]. The resazurin assay measures mitochondrial metabolic activity and correlates well with the number of cells (54).

CCK-8 colorimetric cell viability assay

Similar to the resazurin assay, the CCK-8 assay measures cell metabolic activity and is widely used as a measure of cell viability. The assay was done as per the manufacturer’s protocol (Dojindo Molecular Technologies Inc.). Briefly, 1700 cells were plated at 100-μl volume per well in a 96-well plate, and treatments were done 24 hours after plating. After 3 or 6 days of drug treatment, 10 μl of CCK-8 solution was added to each well and further incubated at 37°C for 1 to 4 hours. This assay quantifies the amount of formazan dye generated by the activity of cellular dehydrogenases, which is directly proportional to the number of living cells. Cell viability was determined by measuring the absorbance at 450 nm with SpectraMax 340 PC (Molecular Devices Corp.).

siRNA transfection

Cells (5 × 105) were plated on 60-mm plates, and transfections with Lipofectamine RNAiMAX (Invitrogen) were performed according to the manufacturer’s protocol. CDK4 siRNA and control siRNA were purchased from Santa Cruz Biotechnology. Combination treatments with drugs were done 24 hours after transfection.

Flow cytometry

For cell cycle analysis, cells were trypsinized, washed, and fixed in 75% ice-cold ethanol after 24 hours of drug treatment. Cells were stained with propidium iodide (50 μg/ml) containing ribonuclease (5 μg/ml) for the measurement of DNA content. Samples were analyzed on a FACScan (Becton Dickinson) for cell cycle distribution with CellQuest software. For this analysis, 10,000 events were examined per sample.

Protein extraction and immunoblotting

Cell lysates were prepared by lysing both floating and adherent cells in radioimmunoprecipitation assay buffer [50 mM tris (pH 7.4), 150 mM NaCl, 1% NP-40, 1 mM EDTA, 0.25% sodium deoxycholate, with protease inhibitor cocktail tablet (Roche)], allowed to lyse on ice for 10 min, syringed, and cleared by centrifugation in a microcentrifuge at 13,000 rpm for 10 min at 4°C. Protein (25 μg) was fractionated by SDS–polyacrylamide gel electrophoresis and transferred onto Immobilon membranes (Millipore). Equal protein loading was confirmed by amido black staining (Bio-Rad). After being blocked with 5% nonfat milk, the membranes were probed with primary antibodies. The following antibodies were used in this study: phosphorylated 4EBP1 (pSer65), phosphorylated AKT (pSer473), AKT, CDK4, GAPDH (glyceraldehyde-3-phosphate dehydrogenase), GSK3α/β (pSer9/pSer21), IGF1R, phosphorylated EGFR (pSer992), phosphorylated MEK (pSer217), PKCα (protein kinase Cα), phosphorylated RB (pSer780), phosphorylated RB (pSer807/pSer811), RB, phosphorylated S6 (pSer235–236), S6, phosphorylated S6K (pThr389), S6K, phosphorylated STAT3 (pTyr705), and phosphorylated Src (pTyr527) from Cell Signaling Technology; phosphorylated ERK (pTyr204) and ERK from Santa Cruz Biotechnology. Bound primary antibodies were detected with horseradish peroxidase–conjugated secondary antibodies (GE Healthcare UK Ltd.) and visualized by enhanced chemiluminescence reagent (GE Healthcare UK Ltd.).

RPPA assay

Drug concentrations for the RPPA screen were based on the IC40 of targets or downstream targets measured by Western blotting or based on information obtained from literature or manufacturer (table S1). IC40 was chosen to maximize the possibility of capturing synergistic and antagonistic effects of drug combinations while operating within the linear dynamic range of the antibody-based measurements. DDLS8817 cells were grown in six-well plates to around 60% confluence. Cells were inhibited with drugs and harvested after 24 hours by collecting and freezing the cell pellet. Nonperturbed control cells were treated with drug vehicle (DMSO) for 24 hours. Cells were thawed and lysed, and protein concentrations were determined by the Bradford assay (Bio-Rad). Protein concentrations were adjusted to 1 to 1.5 mg/ml and denatured in 2% SDS for 5 min at 95°C. Each condition was analyzed in duplicates from two independent biological samples. The RPPA assay was performed at the RPPA core facility at MD Anderson Cancer Center, where cell lysates were spotted on nitrocellulose-coated slides as described previously (39) and stained with about 100 different antibodies.

RPPA data analysis

The antibody staining intensities were quantified with the MicroVigene automated RPPA module (VigeneTech). Data handling consisted of correction of systematic spatial effects due to the location of the sample on the chip (55), conversion of dilution series to a single intensity readout (56), replicate averaging, and z-score scaling within each antibody readout. Z-score scaling was chosen as a normalization method because it robustly captures both activating and inhibiting effects over the entire series of perturbations (>600, including replicates, untreated controls, and additional drugs) and thereby minimizes the risk of systematic errors that may arise when normalizing to a subset of conditions, such as untreated controls. Supporting this, we found a better correlation between biological replicates when z-score normalizing compared to normalizing to untreated controls. After data handling, a subset of the total panel of antibody-stained chips was selected by first removing chips that were unevenly stained, saturated, or underexposed. Second, antibody readouts for which at least 10 conditions produced a z score with an absolute value greater than 1 were retained. Finally, antibodies with a Pearson correlation coefficient greater than 0.5 between biological replicates or antibodies that were present in the prior knowledge network were selected.

CI score

On the basis of cell viability measurements, we used the CI score to determine whether two drugs had synergistic, antagonistic, or additive effects. The CI score is well suited for estimating effects of drug combinations because it is based on the concept of dose substitution and can handle cases where the two drugs are the same, act on the same target, or act on targets converging in a common pathway (15). The CI score, CIX=C1ECX,1+C2ECX,2, measures the fractional shift between the combination doses (C1 and C2) and the single-agent inhibitory concentration for a given level of inhibition (ECX,1 and ECX,2). For this study, we used the EC50, which was the drug concentration that induced a response halfway between the maximum and minimum observed effects of the condition with the largest inhibitory effect (either single or dual). To obtain a confidence estimate, CI scores were also calculated at EC45 and EC55 levels, and SDs of the CI scores were reported (table S3). If a single agent did not reach the chosen effect level (percent inhibition level of the condition with the largest effect), we assumed no effect of this drug, and its contribution to the combination dose became negligible (C1 = 0 or C2 = 0). In cases where neither single nor dual perturbation resulted in an inhibitory effect less than 25% of control levels, CI scores were not determined and reported as “not a number” (NaN). For cell viability experiments where all mixtures of two serially diluted agents (dose matrix) were performed, CompuSyn was used for calculation of CI scores at various effect levels (38).

Network modeling

Cell viability in response to drug intervention is mediated through coordinated changes in the concentrations of proteins and phosphoproteins. We modeled these connected changes as a set of coupled nonlinear ordinary differential equations (ODEs) (fig. S10). Each model variable (or network node) represents the z score of a biological entity, where the z score is a measure of significant change from the average abundance measurement across a large set of drug conditions from which 105 conditions were kept for modeling. Positive and negative values correspond to concentrations above and below the average response, respectively. Our models capture epistasis-like effects and the tendency of a system to return to the baseline state, which we observe to be nearly identical to the untreated as well as average condition. The dynamics are driven by perturbations represented as external forces (μi) on one or more model variables (xi), which drive the system of equations to a unique steady state. The steady state is taken to be the model predicted outcome to the perturbation.

We have incorporated the use of so-called activity nodes, which represent the kinase activity of proteins (30). These activity nodes are introduced to distinguish between phosphorylated and active states of kinases. For example, MEK can be in a phosphorylated but not active state, where phosphorylated MEK does not propagate signals to its downstream effectors. In this example, the kinase activity of MEK (aMEK) is altered with drugs but not directly measureable with antibody. Without readouts of activity, there is no evidence to determine regulators of such activity nodes. Consequently, all activity nodes have only outgoing interactions and represent the points of entry of the various drugs on the network.

Network model parameterization

The models are fully parameterized by an N-by-N weight matrix (W), which contains signed values (wij) representing an interaction from model variable xj on to model variable xi. Given that our model equations are abstractions from the underlying biochemistry of signaling events, we expect that multiple models can fit the collected data set. We define a probabilistic description of model space that rewards fitness to experimental data and agreement with prior knowledge.P(W)=eβiNμM(xiμ*xiμ)2eλiNjiNδ(wij)δ(wij)={0if wij is consistent with prior knowledge1if wij is inconsistent with prior knowledge

The experimental and simulated outcomes for a perturbation condition μ are denoted xiμ* and xiμ. There are KN2 possible configurations of this weight matrix, where K is the number of possible parameter assignments for a single interaction. Explicit calculation of this probability distribution is computationally prohibitive, given the large number of possible model configurations.

Model construction

We adapted the BP inference method for approximating the probability distributions for each possible parameter (wij). A detailed description of our application of BP to perturbation data is provided elsewhere (30). We constructed 1000 high probability models by sampling from the BP-calculated probability distributions (fig. S10). We kept only the top 100 models after parameter refinement with gradient descent (57) and ranking by performance metrics on the training set, and only these 100 low-error models were used for direct simulation and prediction.

Average network models

We used an average network model for qualitative visual analysis of our inference results. An average model is simply the concatenation of all interaction parameters whose expected value from the BP probability distribution is nonzero. For this work, only average edges with interaction strength above or below 0.1 in the set of 1000 models were shown. Note that average network models do not capture mutually exclusive interactions because an average model may contain two interactions that never coexist in the individual models. The sparsity of the networks is tunable through the parameter λ, such that larger values of λ yield sparser networks. Here, we set λ and β equal to 2 and 3, respectively.

Leave-k-out cross-validation

To estimate the genuine power of our models to predict response to drug combinations, we partitioned the full data set into a training set and a test set. The test set consisted of all drug-pair conditions involving a drug of interest, thereby leaving the training set consisting of all other conditions including the drug of interest applied alone. We constructed a unique set of 100 top performing models for the training set, simulated the responses of those models under the conditions in the test set, and compared the results against the data from the test set.

Model prediction of cell viability outcome and synergy

We modeled single-target inhibition with a constant external force (ui as in the equation in fig. S10) on a single node in the interconnected system. Roughly, the value of the external perturbation is related to the expected deviation from average, such that the simulated effect of the cell viability node (ECV) in the presence of inhibition of node A (xCVA) is expressed relative to 0 (average). To obtain synergy scores for a combination of inhibition of nodes A and B, we subjected the models to eight increasing strengths of each perturbation and all 8-by-8 strength combinations. We estimated nonadditive effects by calculating the difference between the simulated effects (Zsim) and the Loewe additivity surface (ZLoewe), derived from adding the effects of single-agent perturbations (43). The synergy score was taken to be the sum of the differences across all 8-by-8 conditions: S=(ZsimZLoewe). The unitless score, S, is only a relative measure of synergy and does not map directly to CI values. Because our modeling approach does not explicitly incorporate drug concentrations, CI scores could not be used as a synergy measure for the computed perturbation effects.

Supplementary Materials

Fig. S1. WDLS/DDLS tumors and the two cell lines used in this study have complex DNA copy number alterations.

Fig. S2. EC50 determination with a cell metabolic activity assay shows good agreement between biological replicates.

Fig. S3. Effects of PD0332991 and NVP-AEW541 in DDLS8817 and LPS141 cells.

Fig. S4. CDK4 inhibition causes G1 cell cycle arrest in LPS141 cells.

Fig. S5. Phosphorylation of AKT in DDLS8817 cells is consistently suppressed for at least 24 hours after PI3K inhibition.

Fig. S6. Examples of determining the appropriate drug concentration for RPPA-based proteomic profiling in DDLS8817 cells.

Fig. S7. Western blot analysis of antibodies used in the RPPA assay shows no apparent cross-reactivity.

Fig. S8. mRNA expression of nodes in the network falls in the top half of all genes in both cell lines used.

Fig. S9. DNA copy number analysis of network genes shows copy gains in multiple genes and amplification of CDK4.

Fig. S10. Schematic illustration of the computational analysis on a fictional four-node system.

Fig. S11. Network models are not overfitted to prior knowledge interactions.

Fig. S12. BP improves the performance of prior knowledge interactions alone.

Fig. S13. BP inference is minimally sensitive to drug specificity.

Fig. S14. Liposarcoma-specific network models are predictive of cell response to drugs.

Fig. S15. Many of the drug combinations with the strongest predicted synergy scores are categorized in accordance with experiments.

Fig. S16. Combined inhibition of the CDK4 and IGF1R nodes is predicted to be synergistic by the network models.

Fig. S17. Combination treatment enhances repression of mTOR signaling compared to single-drug treatment.

Fig. S18. Inhibition of EGFR and CDK4 has synergistic effects on cell metabolic activity, and effects are enhanced in a triple perturbation adding an IGF1R inhibitor.

Fig. S19. Combining CDK4 inhibition with MEK or ERK inhibition does not result in synergistic effects on cell viability based on cell metabolic activity.

Fig. S20. Many of the top 100 models predict synergistic effects of combined CDK4 and IGF1R inhibition.

Table S1. Drugs used in the synergy screen (cell viability) and the proteomic screen (RPPA).

Table S2. Dose-response measurements of single- and paired-drug perturbations with the resazurin assay.

Table S3. Combination index scores.

Table S4. Drugs used in follow-up experiments.

Table S5. Interactions in the prior knowledge network.

Table S6. Bayesian-derived models have many interactions in common with BP-derived models.

References (5896)

References and Notes

Acknowledgments: We thank R. Sinha, M. B. Faura, and J. Saez-Rodriguez for technical support; N. P. Gauthier, P. Kaushik, W.-Q. Wang, S. Nelander, D. Bemis, D. S. Marks, and A. Koff for helpful discussions; and V. A. Pedicord for editing the manuscript. Funding: This work was funded in part by Center for Cancer Systems Biology grant U54 CA148967 (NIH), National Resource for Network Biology grant GM103504 (NIH), SPORE (Specialized Program of Research Excellence) Soft Tissue Sarcoma grant P50 CA140146 (NIH/NCI), and Physical Sciences-Oncology Center grant U54 CA143798 (NIH). Support for E.J.M. was provided by the Tri-Institutional Training Program in Computational Biology and Medicine (NIH training grant 1T32GM083937). Author contributions: Conception and design: M.L.M., E.J.M., S.S., G.K.S., and C.S. Development of methodology: M.L.M., E.J.M., A.K., G.K.S., and C.S. Acquisition of data: M.L.M., J.S.N., T.S., R.S., X.J., and Q.H. Analysis and interpretation of data: M.L.M. and E.J.M. Writing of the manuscript: M.L.M., E.J.M., J.S.N., A.M.C., S.S., G.K.S., and C.S. Administrative, technical, or material support: M.L.M., T.S., R.S., X.J., Q.H., A.M.C., and S.S. Study supervision: G.K.S. and C.S. Competing interests: S.S. is an advisory committee member of Pfizer, and G.K.S. has served on an Advisory Board for Pfizer for the development of PD0332991 in liposarcoma. Data and materials availability: Gene Expression Omnibus: GSE50749 (Illumina gene expression array experiments) and GSE50750 (Agilent 244k Comparative Genomic Hybridization array experiments). RPPA data:
View Abstract

Navigate This Article