Article| Volume 8, ISSUE 2, P97-108.e16, February 27, 2019
• Top

# Quantifying Drug Combination Synergy along Potency and Efficacy Axes

• Author Footnotes
9 These authors contributed equally
Open ArchivePublished:February 20, 2019

## Highlights

• MuSyC is a synergy framework applicable to any metric of drug combination effect
• Unlike other methods, MuSyC decouples synergy of potency and efficacy
• It subsumes traditional synergy methods, resolving ambiguities and biases in the field
• MuSyC reveals optimal co-targeting strategies in NCSLC and melanoma

## Summary

Two goals motivate treating diseases with drug combinations: reduce off-target toxicity by minimizing doses (synergistic potency) and improve outcomes by escalating effect (synergistic efficacy). Established drug synergy frameworks obscure such distinction, failing to harness the potential of modern chemical libraries. We therefore developed multi-dimensional synergy of combinations (MuSyC), a formalism based on a generalized, multi-dimensional Hill equation, which decouples synergistic potency and efficacy. In mutant-EGFR-driven lung cancer, MuSyC reveals that combining a mutant-EGFR inhibitor with inhibitors of other kinases may result only in synergistic potency, whereas synergistic efficacy can be achieved by co-targeting mutant-EGFR and epigenetic regulation or microtubule polymerization. In mutant-BRAF melanoma, MuSyC determines whether a molecular correlate of BRAFi insensitivity alters a BRAF inhibitor’s potency, efficacy, or both. These findings showcase MuSyC’s potential to transform the enterprise of drug-combination screens by precisely guiding translation of combinations toward dose reduction, improved efficacy, or both.

## Introduction

Recent decades have witnessed an exponential expansion of available drugs for the treatment of diseases (
• Gong Z.
• Hu G.
• Li Q.
• Liu Z.
• Wang F.
• Zhang X.
• Xiong J.
• Li P.
• Xu Y.
• Ma R.
• et al.
Compound libraries: recent advances and their applications in drug discovery.
). This expansion has been concomitant with an evolving understanding of disease complexity—complexity commonly necessitating combination therapy (
• He B.
• Lu C.
• Zheng G.
• He X.
• Wang M.
• Chen G.
• Zhang G.
• Lu A.
Combination therapeutics in complex diseases.
). However, clinical applications of combination therapy are often limited by tolerable dose ranges, and, therefore, it is desirable to identify combinations that enable dose reduction (
• Tallarida R.J.
Quantitative methods for assessing drug synergism.
), i.e., synergistic potency. Additionally, combining drugs does not guarantee a priori an increase in efficacy over the single agents, and, therefore, it is desirable to identify combinations with effects greater than what is achievable with either drug alone (
• Foucquier J.
• Guedj M.
Analysis of drug combinations: current methodological landscape.
), i.e., synergistic efficacy. To assess a combination’s performance toward these goals, several drug synergy metrics have been proposed (
• Foucquier J.
• Guedj M.
Analysis of drug combinations: current methodological landscape.
). The roots of current synergy metrics can be traced back to either Loewe, who advanced the dose additivity principle (
• Loewe S.
• Muischnek H.
über Kombination swirkungen.
) or Bliss who first described the multiplicative survival principle (
• Bliss C.I.
The toxicity of poisons applied jointly.
). Nearly a century later, methods to quantify drug synergy continue to appear (
• Chou T.-C.
• Talalay P.
Analysis of combined drug effects: a new look at a very old problem.
,
• Wennerberg K.
• Aittokallio T.
• Tang J.
Searching for drug synergy in complex dose–response landscapes using an interaction potency model.
,
• Twarog N.R.
• Stewart E.
• Hammill C.V.
• Shelat A.A.
BRAID: A unifying paradigm for the analysis of combined drug action.
,
• Zimmer A.
• Katzir I.
• Dekel E.
• Mayo A.E.
• Alon U.
Prediction of multidimensional drug dose responses based on measurements of drug pairs.
,
• Schindler M.
Theory of synergistic effects: hill-type response surfaces as “null-interaction” models for mixtures.
) based on these two principles. However, none of these methods distinguish between synergistic potency and synergistic efficacy. Instead, they either make no distinction or tacitly assume the only form of synergism is through potency.
Nevertheless, this distinction is essential to arrive at an unambiguous definition of synergy and properly rationalize the deployment of drug combinations, e.g., in personalized medicine. Indeed, conflating them may mislead drug combination discovery efforts. For instance, a search for improved efficacy based on traditional synergy frameworks may be confounded by an inability to sort out synergistically potent combinations.
To address this critical shortcoming and resolve these two independent types of synergy, herein we propose a synergy framework termed multi-dimensional synergy of combinations (MuSyC), which is based on a two-dimensional (2D) extension of the Hill equation derived from mass action kinetics. The 2D Hill equation extends dose-response curves to dose-response surfaces.
MuSyC distinguishes between synergistic potency and synergistic efficacy based on parameters in the 2D Hill equation. These synergy parameters are extensions of standard pharmacologic measures of potency and efficacy and define a dose-response surface onto which changes in potency and efficacy are orthogonal. We visualize synergy of potency and efficacy on drug synergy diagrams (DSDs), which globally stratify drug combinations along orthogonal axes of synergy facilitating comparisons between the synergistic profiles of many combinations.
To demonstrate the value of MuSyC, we investigate a panel of anti-cancer compounds in combination with a third-generation mutant-EGFR inhibitor, osimertinib, in EGFR-mutant non-small-cell lung cancer (NSCLC). We find that drugs targeting epigenetic regulators or microtubule polymerization are synergistically efficacious with osimertinib. In contrast, drugs co-targeting kinases in the MAPK pathway affect potency, not efficacy of osimertinib. These conclusions have implications for drug combination deployment in NSCLC where increasing the efficacy of EGFR inhibitors has historically relied on trial and error, with no overarching principles to guide development (
• Schiffmann I.
• Greve G.
• Jung M.
• Lübbert M.
Epigenetic therapy approaches in non-small cell lung cancer: update and perspectives.
).
We also apply MuSyC to study the well-established, clinically relevant combination targeting RAF and MEK in BRAF-mutant melanoma (
• Long G.V.
• Stroyakovskiy D.
• Gogas H.
• Levchenko E.
• de Braud F.
• Larkin J.
• Garbe C.
• Jouary T.
• Hauschild A.
• Grob J.J.
• et al.
Combined BRAF and MEK inhibition versus BRAF inhibition alone in melanoma.
). We find this combination to be synergistically efficacious, though, in several cases, at the cost of potency. We then identify NADPH oxidase 5 (NOX5) as a previously unsuspected molecular determinant of sensitivity to BRAF inhibition (BRAFi) in BRAF-mutant melanoma. Applying MuSyC, we find that NOX5 expression levels affect BRAFi efficacy but not potency.
In direct comparisons, we found that traditional synergy frameworks are biased and ambiguous even for the most synergistically efficacious of the NSCLC and melanoma combination studies, leading to misclassifications of combination synergy. We further show how MuSyC addresses and corrects these problems by generalizing the traditional models.

## Results

### 2D Hill Equation Decouples Synergy of Efficacy from Synergy of Potency

The dose-effect relationship of a single drug is traditionally quantified by the Hill equation, which contains parameters describing efficacy (Emax) and potency (EC50) of a dose-response curve (see STAR Methods for equation derivation; Table 1 for definitions) (Figure 1A). The Hill equation is derived from a phenomenological 2-state model of drug effect (Figure S1A). Therefore, to characterize the dose-effect relationship for drug combinations, we extended this model to a 4-state model (Figure S1B) to derive a 2D generalization of the Hill equation, using principles of mass action kinetics (see STAR Methods). The 2D Hill equation parameterizes a dose-response surface (Figure 1B; Table S1 for parameter descriptions) (
• Greco W.R.
• Bravo G.
• Parsons J.C.
The search for synergy: a critical review from a response surface perspective.
), a 2D extension of 1D dose-response curves (Figure 1A). In this equation, the changes in the efficacy and potency resulting from the combination are quantified by parameters for synergistic efficacy, denoted by β, and synergistic potency, denoted by α (Table S1). These parameters govern the shape of the dose-response surface and can capture complex patterns in experimental data.
Table 1Key Definitions
 Potency The amount of drug required to produce a specified effect. A highly potent drug is active at low concentrations. Classically quantified as the required concentration to achieve half the maximal effect (EC50). Efficacy The degree to which a drug can produce a beneficial effect. Classically quantified as the maximal effect (Emax). Synergistic potency The magnitude of the change in the drug potency, owing to the presence of another drug. Synergistic efficacy The percent change in the maximal efficacy of the combination compared to the most efficacious single agent.
The parameter β is defined as the percent increase in a drug combination’s effect beyond the most efficacious single drug. For instance, in the case of synergistic efficacy (β > 0), the effect at the maximum concentration of both drugs (E3) exceeds the maximum effect of either drug alone (E1 or E2) (Figure 1C, quadrants I and II). For antagonistic efficacy (β < 0) (Figure 1C, quadrants III and IV), at least one or both drugs are more efficacious as single agents than in combination (see Video S1 for an animated example of how the dose-response surface changes as a function of β).
The parameter α quantifies how the effective dose of one drug is altered by the presence of the other. In the case of synergistic potency (α > 1), the EC50 (denoted C in Figure 1B) decreases because of the addition of the other drug (Figure 1C, quadrants I and IV), corresponding to an increase in potency. In the case of antagonistic potency (0 ≤ α < 1), the EC50 of the drug increases as a result of the other drug (Figure 1C, quadrants II and III), corresponding to a decrease in potency (see Video S1 for an animated example of how the dose-response surface changes as a function of α). Since each drug can modulate the effective dose of the other independently (
• Zimmer A.
• Katzir I.
• Dekel E.
• Mayo A.E.
• Alon U.
Prediction of multidimensional drug dose responses based on measurements of drug pairs.
), the 2D Hill equation contains two α values (α1 and α2) (Figure S1B, bottom and right edges of surface). This separation of α values in the 2D Hill equation makes it possible for a given drug combination to have synergism of potency in one direction (α1 > 1), and antagonism of potency in the other (α2 < 1), or vice versa (see Figure S1C for example surfaces).
Both MuSyC parameters for synergy of efficacy (β) and synergy of potency (α) correspond to geometric transformations of the dose-response surface, analogous to the parameters for efficacy (Emax) and potency (EC50) that transform the single-drug dose-response curve in classic pharmacology. We surveyed eight synergy methods to understand how they might account for these distinct types of synergy, including traditional methods of Bliss, Loewe, and highest single agent (HSA) (
Pharmacology.
), as well as more recent frameworks including the combination index (CI) (
• Chou T.-C.
• Talalay P.
Analysis of combined drug effects: a new look at a very old problem.
), Zimmer et al.’s equivalent dose model (
• Zimmer A.
• Katzir I.
• Dekel E.
• Mayo A.E.
• Alon U.
Prediction of multidimensional drug dose responses based on measurements of drug pairs.
), Schindler’s PDE-Hill model (
• Schindler M.
Theory of synergistic effects: hill-type response surfaces as “null-interaction” models for mixtures.
), ZIP (
• Wennerberg K.
• Aittokallio T.
• Tang J.
Searching for drug synergy in complex dose–response landscapes using an interaction potency model.
), and BRAID (
• Twarog N.R.
• Stewart E.
• Hammill C.V.
• Shelat A.A.
BRAID: A unifying paradigm for the analysis of combined drug action.
). We find Bliss, Loewe, HSA, PDE-Hill, ZIP, and BRAID conflate synergy of efficacy and potency (Figures S2A–S2F) so that a drug combination with high synergistic potency scores identical to a combination with high synergistic efficacy (Figure S2A). This conflation, even in methods classically regarded as quantifying exclusively changes in efficacy, such as HSA, underscores the necessity of considering the entire topology of the dose-response surface in order to decouple synergistic efficacy from synergistic potency.
In other methods (equivalent dose and CI), only synergistic potency is tacitly assumed by asserting the maximal effect of each drug and of the combination is equal to zero (Figures S2G–S2J). (See STAR Methods, section Methods Details, subsection Comparison to alternative synergy models, for a case-by-case comparison of MuSyC with other synergy frameworks.)
By using the Hill equation as the basis for MuSyC, the metric of drug effect is not bounded to range between 0 and 1, as is the case for Bliss, CI, and the equivalent dose model, providing a greater versatility for application to other systems. Indeed, the challenges in applying prior synergy frameworks to our recently proposed metric of drug effect, the drug-induced proliferation (DIP) rate (
• Harris L.A.
• Frick P.L.
• Garbett S.P.
• Hardeman K.N.
• Paudel B.B.
• Lopez C.F.
• Quaranta V.
• Tyson D.R.
An unbiased metric of antiproliferative drug effect in vitro.
), provided the initial impetus for developing this framework.
In summary, the 2D Hill equation enables a formalism, termed MuSyC, in which synergistic efficacy and synergistic potency are orthogonal and quantified by the parameters β and α, respectively. We have provided an interactive MuSyC demo (see Data and Software Availability section in STAR Methods,) to facilitate an intuitive understanding of the relationship between different parameter values and the shape of the dose-response surface.

### MuSyC Quantifies Synergy of Potency and Efficacy in a Drug Combination Screen

We applied MuSyC to evaluate the synergistic potency and efficacy of a 64-drug panel (see Table S2 for drugs, drug classes, nominal targets, and tested concentration ranges) in combination with osimertinib, a mutant EGFR-tyrosine kinase inhibitor recently approved for first-line treatment of EGFR-mutant NSCLC (
• Soria J.C.
• Ohe Y.
• Vansteenkiste J.
• Reungwetwattana T.
• Lee K.H.
• Dechaphunkul A.
• Imamura F.
• Nogami N.
• Kurata T.
• et al.
Osimertinib in untreated EGFR -mutated advanced non–small-cell lung cancer.
). The selected drugs span a diverse array of cellular targets that can be broadly grouped into four categories: kinases, receptors and channels, epigenetic regulators, and mitotic checkpoints (Figure 2D), each with several sub-categories. The combinations were tested in PC9 cells, a canonical model of EGFR-mutant NSCLC (
• Jia P.
• Jin H.
• Xia J.
• Ohashi K.
• Liu L.
• Pirazzoli V.
• Dahlman K.B.
• Politi K.
• Michor F.
• et al.
Next-generation sequencing of paired tyrosine kinase inhibitor-sensitive and -resistant EGFR mutant lung cancer cell lines identifies spectrum of DNA changes associated with drug resistance.
) using a high-throughput, in vitro, drug-screening assay (Figure 2A). We quantified drug effect using the DIP rate metric (
• Harris L.A.
• Frick P.L.
• Garbett S.P.
• Hardeman K.N.
• Paudel B.B.
• Lopez C.F.
• Quaranta V.
• Tyson D.R.
An unbiased metric of antiproliferative drug effect in vitro.
), a metric that avoids temporal biases characteristic of traditional endpoint assays (see Quantification and Statistical Analysis section in STAR Methods, section).
To fit the resulting dose-response surfaces, we developed a Bayesian fitting algorithm, using a particle swarm optimizer (PSO) to seed priors for a Markov chain Monte Carlo (MCMC) optimization (Figures S3A and S3B; Quantification and Statistical Analysis section in STAR Methods). The algorithm also accounts for non-optimal drug dosage selection, since dose ranges that are insufficient to observe saturating effects—owing to limited solubility or potency of the drug—result in a commensurate increase in the uncertainty of MuSyC’s synergy parameters (Figures S3C–S3E).
Applying this algorithm, we extracted synergy parameters (α1, α2, and βobs) from fitted surfaces for all osimertinib combinations (βobs is the observed synergistic efficacy at the maximum tested dose range) (see STAR Methods, section Methods Details).
The drug panel displays wide ranges of efficacy (E2) and potency (C) for single agents (Figure S4A). The efficacy and potency of the single agents have no relationship with the synergistic efficacy and synergistic potency when combined with osimertinib (p value > 0.2) (Figure S4B), confirming MuSyC’s synergy parameters are independent of single agents’ dose-response curve and therefore, as expected, cannot be predicted from the single-agent, pharmacologic profiles.
Inspection of dose-response surfaces from this combination screen highlights the significance of resolving synergistic potency and efficacy. For instance, the dose-response surface for the osimertinib combination with M344 (a histone deacetylase [HDAC] inhibitor) exhibits synergistic efficacy (βobs = 1.25 ± 0.03, reflecting a 125% increase in efficacy over osimertinib alone) (Figures 2B and 2E). However, this improved efficacy comes at the cost of potency (log(α2) = −0.90 ± 0.01) as observed in the shift in the EC50 of osimertinib in the presence of 1 uM M344 (Figure 2B; red to purple dotted line). In contrast, ceritinib, an ALK inhibitor with off-target effects on IGF1R (
• Shaw A.T.
• Kim D.W.
• Mehra R.
• Tan D.S.
• Felip E.
• Chow L.Q.
• Camidge D.R.
• Vansteenkiste J.
• Sharma S.
• De Pas T.
• et al.
Ceritinib in ALK -Rearranged non–small-cell lung cancer.
), increases osimertinib’s potency (log(α2) = 6.25 ± 0.50) (Figure 2C; green to orange dotted line) at 4 uM (maximal tested concentration) but with inconsequential improvement of efficacy (βobs = 0.28 ± 0.003).
To visualize synergy globally, we plotted drug combinations on DSDs, with observed synergistic efficacy (βobs) and potency on the vertical and horizontal axes, respectively (Figure 2E). These DSDs reveal distinguishing trends between the four drug categories tested.
Within the mitotic checkpoint drugs, tubulin destabilizers (including vindesine and vinorelbine) showed an upward shift along the axis of synergistic efficacy (Figure 2E). The marginal distribution confirmed this trend in comparison to all the drugs (Figure 2F, blue versus black vertical distributions). Similar results were obtained for the histone deacetylase inhibitor (HDACi) subgroup within the epigenetic regulators (Figures 2E and 2F). As expected, we observed limited synergistic or antagonistic efficacy for drugs targeting G-protein-coupled receptors (GPCRs) (Figures 2E and 2F; red versus black distributions). We also observed limited synergistic efficacy in directly co-targeting kinases in the MAPK pathway, suggesting this may be an unproductive avenue in EGFR-mutant NSCLC (Figures 2E and 2F; purple to black comparison along vertical axis).
In summary, by quantifying synergy of potency separate from synergy of efficacy, MuSyC reveals drug-class trends, which can be used to guide subsequent screens and drug combination deployment in NSCLC.

### MuSyC Validates Co-targeting RAF and MEK in BRAF-Mutant Melanoma

The NSCLC drug screen (Figure 2) suggests combinations targeting molecules within the same signaling pathway may not be productive avenues for increasing efficacy. However, a combination used clinically in BRAF-mutant melanoma co-targets kinases BRAF and MEK in the MAPK pathway (
• Long G.V.
• Stroyakovskiy D.
• Gogas H.
• Levchenko E.
• de Braud F.
• Larkin J.
• Garbe C.
• Jouary T.
• Hauschild A.
• Grob J.J.
• et al.
Combined BRAF and MEK inhibition versus BRAF inhibition alone in melanoma.
,
• Eroglu Z.
• Ribas A.
Combination therapy with BRAF and MEK inhibitors for melanoma: latest evidence and place in therapy.
). To investigate this combination in more detail, we screened a panel of 8 BRAFV600-mutant melanoma cell lines (see
• Paudel B.P.
• Harris L.A.
• Hardeman K.N.
• Abugable A.A.
• Hayford C.E.
• Tyson D.R.
• Quaranta V.
A nonquiescent “idling” population state in drug-treated, BRAF-mutated melanoma.
for cell-line information) against 16 BRAFi/MEKi combinations (see Table S2 for drug information and tested dose ranges).
Based on the mean βobs across cell lines, all 16 combinations were synergistically efficacious (Figures 3A and S5C), indicating MuSyC would have identified this treatment strategy prospectively. In contrast, conventional methods produce ambiguous results (Figure S6, top 3 panels in each cell line group), such that this combination strategy could have not been identified. Furthermore, MuSyC detected variations in synergistic efficacy between cell lines (Figures 3A and S5C), underscoring its sensitivity and pointing to heterogeneous, cell-intrinsic mechanisms modulating the efficacy of combined BRAF and MEK inhibition. In particular, A2058 displayed low average synergistic efficacy, suggesting that its canonical insensitivity to BRAFi does not depend on MEK reactivation but rather on altered metabolic phenotype (
• Parmenter T.J.
• Kleinschmidt M.
• Kinross K.M.
• Bond S.T.
• Li J.
• Rao A.
• Sheppard K.E.
• Hugo W.
• Pupo G.M.
• et al.
Response of BRAF-mutant melanoma to BRAF inhibition is mediated by a network of transcriptional regulators of glycolysis.
,
• Hardeman K.N.
• Peng C.
• Paudel B.B.
• Meyer C.T.
• Luong T.
• Tyson D.R.
• Young J.D.
• Quaranta V.
• Fessel J.P.
Dependence on glycolysis sensitizes BRAF-mutated melanomas for increased response to targeted BRAF inhibition.
).
MuSyC also provides information on synergistic potency for these combinations. A clinically deployed combination (dabrafenib and trametinib) is synergistically efficacious but antagonistically potent in all cell lines except one (Figure S5), a trade-off that may be relevant in the clinic.
Together, MuSyC analyses of NSCLC and of melanoma combination screens indicate that the magnitude of a drug combination’s synergistic efficacy depend upon the oncogenetic context, i.e., co-targeting within the MAPK pathway may work for mutant-BRAF melanoma but not for mutant-EGFR NSCLC.

### MuSyC Reveals Whether Molecular Correlates of Insensitivity Alter Synergistic Efficacy or Potency

While drug combinations are commonly identified from top-down approaches, e.g., high-throughput drug screens, others, including BRAFi/MEKi, were discovered from a bottom-up approach by investigating molecular correlates of insensitivity. However, these molecular correlates may alter either the potency or the efficacy of the primary drug (or both). MuSyC can distinguish among these possibilities, enabling an informed choice between improving either efficacy or potency. As an example, we looked for molecular correlates of BRAFi insensitivity between subclones of a BRAF-mutant melanoma cell line (SKMEL5) with differential sensitivity to BRAFi (Figure 4A). Specifically, we quantified gene expression using RNA sequencing (RNA-seq) and identified the top 200 differentially expressed genes (DEGs) (FDR < 0.001; see STAR Methods section Quantification and Statistical Analysis). This gene set was significantly enriched in processes, cellular components, and molecular functions relating to metabolism (Figure 4B), aligning with previous reports on the relationship between altered metabolism and resistance to BRAFi (
• Parmenter T.J.
• Kleinschmidt M.
• Kinross K.M.
• Bond S.T.
• Li J.
• Rao A.
• Sheppard K.E.
• Hugo W.
• Pupo G.M.
• et al.
Response of BRAF-mutant melanoma to BRAF inhibition is mediated by a network of transcriptional regulators of glycolysis.
,
• Hardeman K.N.
• Peng C.
• Paudel B.B.
• Meyer C.T.
• Luong T.
• Tyson D.R.
• Young J.D.
• Quaranta V.
• Fessel J.P.
Dependence on glycolysis sensitizes BRAF-mutated melanomas for increased response to targeted BRAF inhibition.
). We computed the correlation of the 200 DEGs’ expression to BRAFi sensitivity across a 10-cell-line panel (see STAR Methods) using expression data from
• Subramanian A.
• Narayan R.
• Corsello S.M.
• Peck D.D.
• Natoli T.E.
• Lu X.
• Gould J.
• Davis J.F.
• Tubelli A.A.
• Asiedu J.K.
• et al.
A next generation connectivity map: L1000 platform and the first 1,000,000 profiles.
. NOX5 stood out as one of five genes with a significant, positive correlation with BRAFi insensitivity (Pearson r = 0.65; p value = 0.042) (Figures 4C and 4D; Table S3 for quantification of BRAFi insensitivity and Table S4 for genes correlated with BRAFi insensitivity) and was significantly up-regulated in the BRAFi-insensitive subclone (SC10) compared with the sensitive subclone (SC01) (Figure 4E). Previously unconsidered, NOX5 is an interesting target because of its convergent regulation on metabolic and redox signaling at mitochondria (
• Lu W.
• Hu Y.
• Chen G.
• Chen Z.
• Zhang H.
• Wang F.
• Feng L.
• Pelicano H.
• Wang H.
• Keating M.J.
• et al.
Novel role of NOX in supporting aerobic glycolysis in cancer cells with mitochondrial dysfunction and as a potential target for cancer therapy.
), processes significantly enriched in the DEGs (Figure 4B).
To study NOX5’s contribution to the potency or efficacy of BRAFi, we tested PLX4720 in combination with a NOX5 inhibitor, DPI (
• Jaquet V.
• Marcoux J.
• Forest E.
• Leidal K.G.
• McCormick S.
• Westermaier Y.
• Perozzo R.
• Plastre O.
• Fioraso-Cartier L.
• Diebold B.
• et al.
NADPH oxidase (NOX) isoforms are inhibited by celastrol with a dual mode of action.
), in a panel of 7 melanoma cell lines selected based on differential NOX5 expression. We found synergistic efficacy correlated with NOX5 expression (Pearson r = 0.77; p value = 0.043) (Figures 4G and 4H); however, synergistic potency did not (Pearson r = 0.01; p value = 0.96) (Figures 4G and 4I). Of note, A2058, well known for its resistance to BRAFi, exhibited the highest NOX5 expression among the cell lines and the highest synergistic efficacy (βobs = 1.42 ± 0.05) (Figure 4F), which was more synergistically efficacious than all tested MEKi/BRAFi combinations (Figure 3A).
Taken together, these results suggest co-targeting NOX5 in BRAF-mutant melanoma could lead to improved outcomes for BRAF-mutant melanoma patients with a unique metabolic program for which NOX5 is a biomarker. Furthermore, this study demonstrates the utility of MuSyC for distinguishing a molecular constituent’s role in modulating the potency or efficacy of a drug.

### MuSyC Generalizes Traditional Synergy Metrics and Removes Biases and Ambiguities

To investigate how results from MuSyC compare with the most frequently used synergy metrics, we calculated synergy using Loewe additivity, CI, and Bliss on data from the NSCLC (Figure 2) and the melanoma (Figure 3A) screens. Loewe synergy was calculated directly from the DIP rate data, while CI and Bliss, which require percent metrics, were calculated from 72-h percent viability (
• Barretina J.
• Caponigro G.
• Stransky N.
• Venkatesan K.
• Margolin A.A.
• Kim S.
• Wilson C.J.
• Lehár J.
• Kryukov G.V.
• Sonkin D.
• et al.
The Cancer Cell Line Encyclopedia enables predictive modelling of anticancer drug sensitivity.
) imputed from the growth curves (see STAR Methods section Quantification and Statistical Analysis). Unlike MuSyC, these metrics are evaluated at every concentration, resulting in dose-dependent distributions of synergy (Figures 5A and S6) commonly resulting in an ambiguous classification of a combination. By the median of each distribution, none of the metrics can statistically discriminate between the MuSyC DSD quadrants (Figures 5A and S6; Kruskal-Wallis p value > 0.05).
Examining the models underlying these metrics revealed several limitations and biases accounting for their ambiguity. For Loewe additivity, synergy is undefinable for many tested concentrations as Loewe cannot be calculated at combination conditions with effects exceeding the maximum effect of the weaker drug (
• Foucquier J.
• Guedj M.
Analysis of drug combinations: current methodological landscape.
). This is particularly limiting for synergistically efficacious combinations, which, by definition, achieve greater effect than either drug alone. In the NSCLC screen, because osimertinib alone was not sufficient to achieve a negative DIP rate (i.e., regressing population), Loewe is undefinable for all conditions where the DIP rate was less than zero (Figure 5B). For conditions where Loewe is defined, Loewe additivity has been reported to be most appropriate for combinations of mutually exclusive inhibitors (
• Chou T.C.
• Talalay P.
Quantitative analysis of dose-effect relationships: the combined effects of multiple drugs or enzyme inhibitors.
). Accordingly, we found Loewe emerges from MuSyC as a special case under the conditions of both α1 = α2 = 0 (i.e., the drugs are mutually exclusive) and h1 = h2 = 1 (see STAR Methods, section Methods Details, subsection 2.1). If the condition h1 = h2 = 1 is not satisfied (Figure 5C), MuSyC predicts that when the geometric mean of the hill slopes is less than 1 $(h1∗h2<1)$, the linear model of Loewe will overestimate synergy and when $h1∗h2>1$, Loewe will underestimate synergy (Figure 5C). Correspondingly, we found the median value of Loewe synergy was negatively correlated with the geometric mean of the hill coefficients in both the NSCLC and melanoma screens (Figure 5D, Spearman r = −0.51 and −0.41; p value = 1e−3 and 8e−4, respectively); that is, the synergy of a combination according to Loewe additivity could be estimated based on the hill slope of a single drug alone in contrast to MuSyC where synergistic potency and efficacy are decoupled from the single drug’s pharmacologic profile (Figure S3B).
CI is a special case of Loewe additivity that adds the additional condition that E0 = 1, E1 = E2 = E3 = 0, such that the drug effect is equated with percent inhibition (
• Chou T.-C.
• Talalay P.
Analysis of combined drug effects: a new look at a very old problem.
). The condition on effect range assumes all drugs achieve the same maximum effect, and thus, unlike Loewe additivity, CI range is not limited by the weaker drug. However, in percent viability data, many drugs do not achieve 0% viability (e.g., methotrexate, which reaches a maximum effect of 52% viability) (Figure 5E). In these cases, fits for the single-drug dose-response curves used to calculate CI are poor (Figure 5E). CI is thus inappropriate for cell-based assays of drug effect where the correspondence between percent inhibition and cell viability is not one-to-one.
Bliss, similar to CI, can only be applied to percent metrics with the condition E0 = 1, E1 = E2 = E3 = 0. As with CI, because most drugs in combination do not satisfy this condition, Bliss is also an inappropriate model to use. However, if this condition is satisfied, Bliss emerges as a special case of MuSyC under the conditions α1 = α2 = 1 (see STAR Methods section Methods Details, subsection 2.2).
In summary, MuSyC subsumes Loewe (and therefore CI) and Bliss into a single framework satisfying both the dose additivity and the multiplicative survival principles under certain conditions. For combinations that do not satisfy these conditions, we show the traditional metrics lead to biased and ambiguous results, while MuSyC’s generality resolves these limitations. Specifically these limitations are as follows traditional methods cannot distinguish synergy of potency from synergy of efficacy (Figures 5A and S2); Loewe is undefined for combinations with synergistic efficacy (Figure 5B); Loewe (and by extension CI) contain an artificial bias toward synergy for drugs with hill slopes much less than 1 (Figures 5C and 5D); and CI leads to poor fits because it disregards synergistic efficacy by assuming that the maximal effect of a drug reaches 0%, even when this is not the case (Figure 5E).

## Discussion

The goal of using synergistic drugs is to achieve more with less. It is therefore intuitive that two types of synergy exist: one corresponding to how much more is achievable (synergistic efficacy) and the other to how much less is required (synergistic potency). Finding such combinations is vital for optimizing therapeutic windows, as there exists a fundamental trade-off between clinical efficacy and tolerable doses. Diseases for which single-drug efficacy is sufficient would benefit from synergistically potent combinations to drive down toxicity and/or side effects. Diseases with treatments of insufficient efficacy are in pressing need of synergistically efficacious combinations in order to improve the depth and durability of response. By stratifying synergy along distinct axes of potency and efficacy using MuSyC, informed choices can be made about this trade-off. The distinction facilitates identifying drug-class trends that can be iteratively expanded in future screens to optimize synergistic efficacy or synergistic potency, whichever is desirable for a particular disease.
In this respect, MuSyC provides a global view of the synergistic behavior of whole classes of drugs, e.g., from a high-throughput drug screen, via DSDs. In this work, MuSyC revealed a subclass of epigenetic regulators as potentially interesting targets for combination therapy in an EGFR-oncogene-addicted background. Epigenetic regulators have previously been suggested to prime NSCLC for sensitivity to EGFRi (
• Schiffmann I.
• Greve G.
• Jung M.
• Lübbert M.
Epigenetic therapy approaches in non-small cell lung cancer: update and perspectives.
), and the HDACi entinostat in combination with erlotinib (first generation EGFR-TKI) has been shown to increase overall survival in EGFR-mutant NSCLC cases with high expression of E-cadherin (
• Witta S.E.
• Gemmill R.M.
• Hirsch F.R.
• Coldren C.D.
• Hedman K.
• Ravdel L.
• Helfrich B.
• Chan D.C.
• Sugita M.
• et al.
Restoring E-cadherin expression increases sensitivity to epidermal growth factor receptor inhibitors in lung cancer cell lines.
,
• Witta S.E.
• Jotte R.M.
• Konduri K.
• Neubauer M.A.
• Spira A.I.
• Ruxer R.L.
• Varella-Garcia M.
• Bunn P.A.
• Hirsch F.R.
Randomized phase II trial of erlotinib with and without entinostat in patients with advanced non-small-cell lung cancer who progressed on prior chemotherapy.
). Consistent with this, we also observe entinostat was synergistically efficacious with osimertinib (βobs = 0.84 ± 0.027) in PC9 cells, an E-cadherin high-expressing cell line (
• Shimoyama Y.
• Nagafuchi A.
• Fujita S.
• Gotoh M.
• Takeichi M.
• Tsukita S.
• Hirohashi S.
Cadherin dysfunction in a human cancer cell line: possible involvement of loss of alpha-catenin expression in reduced cell-cell adhesiveness.
). As is typical of high-throughput screens, there were results of undetermined significance, including dronedarone (an anti-arrhythmic sodium channel inhibitor) and GW694590a (an anti-angiogenesis compound targeting the TIE2 receptor), which were the most antagonistic and synergistically efficacious compounds out of the receptors and channels drug classes, respectively. Further studies are needed to verify these results. Nonetheless, MuSyC provides a quantitative foundation to further investigate unsuspected combinations.
The global views provided by the MuSyC DSDs also reveal synergistic trends that vary according to disease context. For example, co-targeting the MAPK pathway in NSCLC or BRAF-mutant melanoma yields different outcomes: in the former, only synergistic potency is observed, while in the latter, synergistic efficacy, and sometimes potency, are registered. The disparity emphasizes that synergistic trends require data-driven metrics that distinguish between synergy of efficacy and potency.
MuSyC dose-response surfaces facilitate evaluating the significance that combination synergy should be assigned; that is, MuSyC’s synergy parameters quantify the relative increase in efficacy or potency of the combination, with respect to single agents, and therefore, the improvements should be interpreted in the context of the absolute potency and efficacy. This information is directly conveyed in the topology of the dose-response surface. As an example, in the NSCLC screen, the combination of osimertinib with quisinostat exhibited the greatest total efficacy. However, since quisinostat is already significantly efficacious on its own, that combination ranks lower than the M344-osimertinib combination along the axis of synergistic efficacy on a DSD. Thus, DSDs are useful to rank relative increase in potency or efficacy, whereas surfaces convey the absolute efficacy and potency achieved by a combination.
MuSyC is also useful for investigating a molecular species’ contribution to the potency and efficacy of a compound. Here, we demonstrated NOX5 activity modulates the efficacy but not the potency, of BRAFi. However, the NOX5i used, DPI, is known to have off-target effects (
• Altenhöfer S.
• Kleikers P.W.
• Wingler K.
• Schmidt H.H.
Evolution of NADPH oxidase inhibitors: selectivity and mechanisms for target engagement.
); therefore, further evidence for the role of NOX5 in BRAFi efficacy will require extending MuSyC to studies combining drugs and gene silencing technology (e.g., RNAi or CRISPR).
To fit the dose response surface and extract synergy parameters, MuSyC utilizes a Bayesian approach combining PSO and a multi-tier MCMC walk in order to track uncertainty in the values for synergistic potency and efficacy. The sources for this uncertainty include noise, partial dose-response curves, and data density. A similar Bayesian approach was previously implemented for Loewe (
• Hennessey V.G.
• Rosner G.L.
• Bast Jr., R.C.
• Chen M.
A Bayesian approach to dose-response assessment and synergy and its application to in vitro dose-response studies.
).
Loewe additivity and Bliss independence have maintained dominance in the field, along with the related work of Chou and Talalay. Yet there is no consensus regarding the appropriate use of these methods because they are based on distinct foundational principles, often leading to incompatible results (
• Greco W.
• Unkelbach H.-D.
• Pöch G.
• Sühnel J.
• Kundi M.
• Bödeker W.
Consensus on concepts and terminology for combined-action assessment: the Saariselka Agreement.
). MuSyC removes these sources of confusion by unifying these methods into a consensus framework, within which Loewe and Bliss emerge as special cases.
There has been much critical analysis over the past 25 years on the term “synergy” (
• Greco W.
• Unkelbach H.-D.
• Pöch G.
• Sühnel J.
• Kundi M.
• Bödeker W.
Consensus on concepts and terminology for combined-action assessment: the Saariselka Agreement.
), arguably rooted in the practice of defining synergy with respect to arbitrary expectations of drug additivity implicitly codified in previous methods’ foundational principles. In contrast, ambiguity about the meaning of synergy disappears in MuSyC because its synergy parameters relate directly to the textbook pharmacology concepts of efficacy and potency. Indeed, a major advance of MuSyC is the decisive shift toward synergy calculations directly related to an observable change in efficacy and/or potency. Thus, ambiguous questions such as “Is there synergy?” can be recast into more precise questions, such as “How much does efficacy or potency of drug X change when drug Y is added?” Such precise language should promote a move away from arbitrary cutoffs for “significant synergy,” which are context dependent.
While we focused on the DIP rate as our metric of effect, MuSyC may be applied to any quantifiable phenotype whose dose response is suitable to be fit by a Hill equation. In contrast, all other synergy models we surveyed impose strict constraints on the type and/or magnitude of the drug effect metric. Thus, MuSyC opens up the potential to study synergy of drug effects previously impossible to address by existing methods. Examples of metrics include immune activation, growth in 3D culture, or second messenger efflux. The flexibility is particularly critical in translating drug combinations to the clinic by using models of increasing complexity, such as organoids, which better represent the drug sensitivity of a patient (
• Jabs J.
• Zickgraf F.M.
• Park J.
• Wagner S.
• Jiang X.
• Jechow K.
• Kleinheinz K.
• Toprak U.H.
• Schneider M.A.
• Meister M.
• et al.
Screening drug effects in patient-derived cancer cells links organoid responses to genome alterations.
). Indeed, that most clinical combinations can be explained by patient-to-patient variability (
• Palmer A.C.
• Sorger P.K.
Combination.
) is a strong rationale for translating combination screens to more complex, pre-clinical models. Subsequent work will be devoted to scaling the combination drug screening pipeline developed here to pre-clinical experimental models of increasing complexity, such as organoids.
In conclusion, we have presented MuSyC, a drug synergy framework that maintains a distinction between two intuitive types of pharmacological synergy and that may be applied to any drug effect metric. We showed this framework allows for a richer understanding of drug interactions, with practical, translational consequences. We foresee this approach will streamline drug discovery pipelines and facilitate the deployment of precision approaches to therapeutic combinations.

## STAR★Methods

### Key Resources Table

Tabled 1
REAGENT or RESOURCESOURCEIDENTIFIER
Chemicals, Peptides, and Recombinant Proteins
TrizolInvitrogen15596026
FBSGibco10437-028
PBSCorning21-040-CV
DMEMGibco11965-092
RPMICorning10-040-CV
TryplEGibco12604-013
DMEM/F12Gibco11330-032
SytoxGreenThermoFisherS7020
5-IodotubericidinENZOEI-293
Abexinostat (PCI-24781)SelleckChemS1090
AcetylcysteineSelleckChemS1623
Afatinib (BIBW2992)LC LaboratoriesA8644
AG-879ENZOEI-258
Alisertib (MLN8237)MedChemExpressHY-10971
Amiodarone HClSelleckChemS1979
AprepitantSelleckChemS1189
Bazedoxifene HClSelleckChemS2128
Beclomethasone dipropionateLight Biologicals (NIH Clinical Collection II)MZ-3012
BendroflumethiazideLight Biologicals (NIH Clinical Collection II)B-8008
BML-259ENZOEI-344
Bosutinib (SKI-606)LC LaboratoriesB-1788
Brigatinib (AP26113)SelleckChemS8229
Buparlisib (BKM120, NVP-BKM120)SelleckChemS2247
CabozantinibLC LaboratoriesC-8901
CarfilzomibLC LaboratoriesC-3022
CarmustineNCI Chemotherapeutic Agents Repository409962
CephalomannineSelleckChemS2408
Ceritinib (LDK378)SelleckChemS7083
CisplatinSigma470306
CobimetinibMedChemExpressHY-13064
CrizotinibLC LaboratoriesC-7900
DabrafenibLC LaboratoriesD-5678
DactolisibLC LaboratoriesN-4288
DasatinibLC LaboratoriesD-3307
DocetaxelSelleckChemS1148
Dronedarone HCl (Multaq)SelleckChemS2114
Ensartinib (X-396)SelleckChemS8230
Entinostat (MS-275)SelleckChemS1053
ErlotinibLC LaboratoriesE-4007
Foretinib (GSK1363089)SelleckChemS1111
Gefitinib (ZD1839)LC LaboratoriesG-4408
Givinostat (ITF2357)SelleckChemS2170
GSK1751853AGSK PKISN/A
GSK994854AGSK PKISN/A
GW458787AGSK PKISN/A
GW644007XGSK PKISN/A
GW694590AGSK PKISN/A
GW770249X (GW770249A)GSK PKISN/A
Homoharringtonine (Omacetaxine mepesuccinate)Sequoia Research Products Ltd. (NIH Clinical Collection II)SRP02125h
Ivacaftor (VX-770)SelleckChemS1144
(+)-JQ1SelleckChemS7110
Linsitinib (OSI-906)SelleckChemS1091
LY294002ENZOST-420
M344SelleckChemS2779
MethotrexateMedChemExpressOL-14377
MG-132SelleckChemS2619
ML-9-HClENZOEI-153
Mocetinostat (MGCD0103)SelleckChemS1122
NaftopidilSelleckChemS2126
NateglinideSelleckChemS2489
Nebivolol HClSelleckChemS1549
Olaparib (AZD2281, Ku-0059436)LC LaboratoriesO-9201
Osimertinib (AZD9291)SelleckChemS7297
PaclitaxelSigma17191
PanobinostatNCI Chemotherapeutic Agents Repository761190
2'-Amino-3'-methoxyflavoneLC LaboratoriesP-4313
PimobendanSelleckChemS1550
PLX-4720SelleckChemS1152
Ponatinib (AP24534)LC LaboratoriesP-7022
PP2ENZOEI-297
Pracinostat (SB939)SelleckChemS1515
Primaquine DiphosphateSelleckChemS4237
QuercetinENZOAC-1142
Quisinostat (JNJ-26481585)SelleckChemS1096
RAF-265MedChemExpressHY-10248
Rapamycin (Sirolimus)SelleckChemS1039
SB-253226GSK PKISN/A
Selumetinib (AZD-6244)LC LaboratoriesS-4490
SP 600125ENZOEI-305
Sunitinib MalateSelleckChemS1042
TAK-632SelleckChemS7291
Tanespimycin (17-AAG)SelleckChemS1141
Thioridazine hydrochlorideSelleckChemS5563
TrametinibLC LaboratoriesT8123
AG-370ENZOEI-229
U-0126ENZOEI-282
Ulixertinib (BVD-523, VRT752271)SelleckChemS7854
Vemurafenib (PLX4032)SelleckChemS1267
VerteporfinSelleckChemS1787
VindesineSequoia Research Products Ltd. (NIH Clinical Collection)SRP01038v
Vinorelbine TartrateSelleckChemS4269
ZM 447439SelleckChemS1103
Critical Commercial Assays
Tru-Seq stranded mRNA sample prep kitIlluminaCat # RS-122-2101
Reverse Transcription KitQuantiTectCat # 205311
IQTM SYBR Green SupermixBioRadCat # 170
Deposited Data
Fitted combination surface plotsThis Paperhttps://github.com/QuLab-VU/MuSyC_Cell.git; In folder(s): Code_Paper_Figures/Fig2(3)/html
Code for Generating Paper PlotsThis Paperhttps://github.com/QuLab-VU/MuSyC_Cell.git; In folder(s): Code_Paper_Figures/
Table of fitted parameters for all experimentsThis Paperhttps://github.com/QuLab-VU/MuSyC_Cell.git; In folder: Data; Files: MasterResults.csv and MasterResults_plx_dpi_melPanel
RT-qPCR quantification of NOX5 expressionThis Paperhttps://github.com/QuLab-VU/MuSyC_Cell.git; In file: Data/nox5Expr.csv
DIP Rate CalculationsThis Paperhttps://github.com/QuLab-VU/MuSyC_Cell.git; In folder(s): Data; Files: HTS018_rates, HTS022_timeavg_rates_sub2.csv, -03-27-2018-dpi+plx-cm_bp_timeavg_preCalcDIP_timSub.csv, dasatinib_osimertinib_cellavista_cm_8-24-17.csv, linsitinib_osimertinib_cellavista_cm_8-24-17.csv, HTS015_017_Combined.csv
cFP Raw DataThis Paperhttps://github.com/QuLab-VU/MuSyC_Cell.git; In folder(s): Data/SKMEL5_cFP
List of DEGSThis Paperhttps://github.com/QuLab-VU/MuSyC_Cell.git; In folder(s): Data/DEGs_GO_Analysis
Raw RNAseq data for subclonesGEOGEO: GSE122041
This data is also available from Mendeley Data at the following doiThis Paper[https://doi.org/10.17632/n8bp8db5ff.1]
Experimental Models: Cell Lines
PC9-H2B.RFP
• Tyson D.R.
• Garbett S.P.
• Frick P.L.
• Quaranta V.
Fractional proliferation: a method to deconvolve cell population dynamics from single-cell data.
(W. Pao at UPenn)
N/A
SKMEL5-H2B.RFP
• Paudel B.P.
• Harris L.A.
• Hardeman K.N.
• Abugable A.A.
• Hayford C.E.
• Tyson D.R.
• Quaranta V.
A nonquiescent “idling” population state in drug-treated, BRAF-mutated melanoma.
(ATCC)
HTB-70
WM1799-H2B.RFP
• Paudel B.P.
• Harris L.A.
• Hardeman K.N.
• Abugable A.A.
• Hayford C.E.
• Tyson D.R.
• Quaranta V.
A nonquiescent “idling” population state in drug-treated, BRAF-mutated melanoma.
(M. Herlyn at Wistar Institute)
N/A
WM983B-H2B.RFP
• Paudel B.P.
• Harris L.A.
• Hardeman K.N.
• Abugable A.A.
• Hayford C.E.
• Tyson D.R.
• Quaranta V.
A nonquiescent “idling” population state in drug-treated, BRAF-mutated melanoma.
(M. Herlyn at Wistar Institute)
N/A
A375-H2B.RFP-FUCCI
• Paudel B.P.
• Harris L.A.
• Hardeman K.N.
• Abugable A.A.
• Hayford C.E.
• Tyson D.R.
• Quaranta V.
A nonquiescent “idling” population state in drug-treated, BRAF-mutated melanoma.
(ATCC)
CRL-1619
SKMEL28-H2B.RFP-FUCCI
• Paudel B.P.
• Harris L.A.
• Hardeman K.N.
• Abugable A.A.
• Hayford C.E.
• Tyson D.R.
• Quaranta V.
A nonquiescent “idling” population state in drug-treated, BRAF-mutated melanoma.
(ATCC)
HTB-72
WM2664-H2B.RFP
• Paudel B.P.
• Harris L.A.
• Hardeman K.N.
• Abugable A.A.
• Hayford C.E.
• Tyson D.R.
• Quaranta V.
A nonquiescent “idling” population state in drug-treated, BRAF-mutated melanoma.
(M. Herlyn at Wistar Institute)
N/A
A2058-H2B.RFP
• Paudel B.P.
• Harris L.A.
• Hardeman K.N.
• Abugable A.A.
• Hayford C.E.
• Tyson D.R.
• Quaranta V.
A nonquiescent “idling” population state in drug-treated, BRAF-mutated melanoma.
(ATCC)
CRL-11147
SKMEL5.SC10-H2B.RFP
• Paudel B.P.
• Harris L.A.
• Hardeman K.N.
• Abugable A.A.
• Hayford C.E.
• Tyson D.R.
• Quaranta V.
A nonquiescent “idling” population state in drug-treated, BRAF-mutated melanoma.
Derived Subclone from SKMEL5-H2B.RFP
N/A
SKMEL5.SC07-H2B.RFP
• Paudel B.P.
• Harris L.A.
• Hardeman K.N.
• Abugable A.A.
• Hayford C.E.
• Tyson D.R.
• Quaranta V.
A nonquiescent “idling” population state in drug-treated, BRAF-mutated melanoma.
Derived Subclone from SKMEL5-H2B.RFP
N/A
SKMEL5.SC01-H2B.RFP
• Paudel B.P.
• Harris L.A.
• Hardeman K.N.
• Abugable A.A.
• Hayford C.E.
• Tyson D.R.
• Quaranta V.
A nonquiescent “idling” population state in drug-treated, BRAF-mutated melanoma.
Derived Subclone from SKMEL5-H2B.RFP
N/A
Oligonucleotides
NOX5_Forward Primer: GGCTCAAGTCCTACCACTGGAThis paperN/A
NOX5_Reverse Primer: GAACCGTGTACCCAGCCAATThis paperN/A
HPRT_Forward Primer: TGCTCGAGATGTGATGAAGGAGThis paperN/A
HPRT_Reverse Primer: TGATGTAATCCAGCAGGTCAGCThis paperN/A
36B4_Forward Primer: CATGTTGCTGGCCAATAAGGThis paperN/A
36B4_Reverse Primer: TGGTGATACCTAAAGCCTGGAAThis paperN/A
PGC1a_Forward Primer: TGCCCTGGATTGTTGACATGAThis paperN/A
PGC1a_Reverse Primer: TTTGTCAGGCTGGGGGTAGGThis paperN/A
Recombinant DNA
• Welm B.E.
• Dijkgraaf G.J.
• Bledau A.S.
• Welm A.L.
• Werb Z.
Lentiviral transduction of mammary stem cells for analysis of gene function during development and cancer.
Plasmid #:18982
Software and Algorithms
Scikit-learn
• van der Walt S.
• Schönberger J.L.
• Nunez-Iglesias J.
• Boulogne F.
• Warner J.D.
• Yager N.
• Gouillart E.
• Yu T.
scikit-image contributors
Scikit-image: image processing in Python.
N/A
RabbitMQ/Celerywww.celeryproject.orgN/A
GNU parallel
• Tange O.
GNU parallel: the command-line power tool.
N/A
Scipy
• Jones E.
• Oliphant T.
• Peterson P.
SciPy: open source scientific tools for python.
N/A
Matplotlib
• Hunter J.D.
Matplotlib: a 2D graphics environment.
N/A
Pandas

McKinney, W. (2010) Data structures for statistical computing in python, Proceedings of the 9th Python in Science Conference, pp. 51–56.

N/A
Numpy
• Oliphant T.E.
Guide to NumPy.
N/A
Pymc3
• Salvatier J.
• Wiecki T.V.
• Fonnesbeck C.
Probabilistic programming in Python using PyMC3.
N/A
HISAT2
• Kim D.
• Salzberg S.L.
HISAT: a fast spliced aligner with low memory requirements.
N/A
featureCounts
• Liao Y.
• Smyth G.K.
• Shi W.
featureCounts: an efficient general purpose program for assigning sequence reads to genomic features.
N/A
Bioconductor (R)www.bioconductor.orgN/A
ENRICHR (R)
• Kuleshov M.V.
• Jones M.R.
• Rouillard A.D.
• Fernandez N.F.
• Duan Q.
• Wang Z.
• Koplev S.
• Jenkins S.L.
• Jagodnik K.M.
• Lachmann A.
• et al.
Enrichr: a comprehensive gene set enrichment analysis web server 2016 update.
N/A
DESeq2
• Love M.I.
• Huber W.
• Anders S.
Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.
.
N/A

### Contact for Reagent and Resource Sharing

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Vito Quaranta ([email protected]).

### Experimental Model and Subject Details

PC9 (previously PC-14, gender unknown) cells were obtained from W. Pao (U Penn.) and were cultured in RPMI 1640 medium containing 10% FBS at 37C and 5% CO2. Cells were engineered to express histone 2B-mRFP via lenti-viral transfection using the pHIV-H2B-mRFP plasmid (
• Welm B.E.
• Dijkgraaf G.J.
• Bledau A.S.
• Welm A.L.
• Werb Z.
Lentiviral transduction of mammary stem cells for analysis of gene function during development and cancer.
) as previously described (
• Tyson D.R.
• Garbett S.P.
• Frick P.L.
• Quaranta V.
Fractional proliferation: a method to deconvolve cell population dynamics from single-cell data.
). A single-cell derived clonal population demonstrated to exhibit proliferation characteristic of the parental population was then selected by limiting dilution. BRAFV600-mutant melanoma lines cells (A2058 (M), WM1799 (U), A375 (F), WM983B (M), SKMEL5 (F), SKMEL28 (M), WM2664 (F), M=Male, F=Female, U=Unknown) were obtained from ATCC or M. Herlyn (Wistar Institute)(see Key Resources Table) and were cultured in DMEM media containing 2 mM glutamine, 4.5 g/L glucose, 10% FBS and no sodium pyruvate (catalog 11965-092) as previously described (
• Hardeman K.N.
• Peng C.
• Paudel B.B.
• Meyer C.T.
• Luong T.
• Tyson D.R.
• Young J.D.
• Quaranta V.
• Fessel J.P.
Dependence on glycolysis sensitizes BRAF-mutated melanomas for increased response to targeted BRAF inhibition.
). SKMEL5.SC10, SKMEL5.SC07, and SKMEL5.SC01 are single cell-derived subclones from SKMEL5 (
• Paudel B.P.
• Harris L.A.
• Hardeman K.N.
• Abugable A.A.
• Hayford C.E.
• Tyson D.R.
• Quaranta V.
A nonquiescent “idling” population state in drug-treated, BRAF-mutated melanoma.
). Cell lines were tested for mycoplasma before each experiment.

### Methods Details

#### Key Equations

For full derivation of these equations, see subsection 4 Derivation of generalized 2-dimensional hill equation. This section is meant to serve as a quick reference guide for the main equations used in the paper.
If the behavior of the drugs in the model formulation of Figure S1B obey detailed balance, then the effect of the combination (i.e., the height of combination surface) is described by
$Ed=C1h1C2h2E0+d1h1C2h2E1+C1h1d2h2E2+(α2d1)h1d2h2E3C1h1C2h2+d1h1C2h2+C1h1d2h2+(α2d1)h1d2h2,$
(Equation 1)

where Ed represents the expected effect for a given dose pair d1,d2 and is specified with 9 parameters defined in Table S1. In addition, detailed balance enforces the constraint that
$α1h2=α2h1.$
(Equation 2)

α is a unitless scalar transforming dose d into an effective dose αd and is used to quantify synergistic potency in MuSyC.
Synergistic efficacy (β) is calculated from E0,E1,E2,E3. β is defined in Equation 3 and is interpreted as the percent increase in maximal efficacy of the combination over the most efficacious single agent. The observed β at the maximum of tested concentrations is defined in Equation 4.
$β=min(E1,E2)−E3E0−min(E1,E2)$
(Equation 3)

$βobs=min[E1(d1max),E2(d2max)]−E3(d1max,d2max)E0−min[E1(d1mx),E2(d2max)]$
(Equation 4)

Equation 1 can be re-written to include β explicitly by replacing E3 with min(E1,E2)−β(E0min(E1,E2)) resulting in the following equation.
$Ed=C1h1C2h2E0+d1h1C2h2E1+C1h1d2h2E2+(α2d1)h1d2h2(min(E1,E2)−β∗(E0−min(E1,E2))C1h1C2h2+d1h1C2h2+C1h1d2h2+(α2d1)h1d2h2$
(Equation 5)

For drugs that do not follow detailed balance, we have derived a more general formulation with 12 parameters:
$Ed=[E0E1E2E3]⋅[−(r1d1h1+r2d2h2)r−1r−20r1d1h1−(r−1+r2(α1d2)h2)0r−2r2d2h20−(r1(α2d1)h1+r−2)r−11111]−1⋅[0001]$
(Equation 6)

where again E3 can be replaced to include β explicitly.
Because we do not know a priori whether combinations will follow detailed balance, we use an information theoretic approach to pick the best model for the data. We have defined six tiers of model complexity, and the best model is selected based on minimizing the deviance information criterion. (See section Quantification and Statistical Analysis, subsection 1 Fitting Dose-response Surfaces for description of fitting algorithm and Table S5 for description of model tiers).

#### Comparison to Alternative Synergy Models

Several other methods for calculating synergy exist, including long-standing traditional methods Loewe (
• Loewe S.
• Muischnek H.
über Kombination swirkungen.
,
• Loewe S.
Versuch einer allgemeinen Pharmakologie der Arznei- Kombinationen.
), Bliss (
• Bliss C.I.
The toxicity of poisons applied jointly.
), HSA (
Pharmacology.
,
• Greco W.R.
• Bravo G.
• Parsons J.C.
The search for synergy: a critical review from a response surface perspective.
), and CI (
• Chou T.C.
• Talalay P.
Quantitative analysis of dose-effect relationships: the combined effects of multiple drugs or enzyme inhibitors.
), as well as more recent methods such as ZIP (
• Wennerberg K.
• Aittokallio T.
• Tang J.
Searching for drug synergy in complex dose–response landscapes using an interaction potency model.
), BRAID (
• Twarog N.R.
• Stewart E.
• Hammill C.V.
• Shelat A.A.
BRAID: A unifying paradigm for the analysis of combined drug action.
), the effective dose model (
• Zimmer A.
• Katzir I.
• Dekel E.
• Mayo A.E.
• Alon U.
Prediction of multidimensional drug dose responses based on measurements of drug pairs.
), and Schindler’s Hill-PDE model (
• Schindler M.
Theory of synergistic effects: hill-type response surfaces as “null-interaction” models for mixtures.
). All of these methods, as well as our own, define a null surface. Combinations with effects greater than or less than expected based on the null surface are deemed synergistic or antagonistic respectively. These methods broadly use one of two approaches to quantify synergy. Loewe, Bliss, CI, HSA, Schindler’s Hill-PDE, and ZIP quantify synergy at every concentration based on how the experimentally measured response deviates from the null surface. BRAID, the effective dose model, and MuSyC provide equations with synergy parameters describing the entire surface which is fit to experimental data using non-linear curve-fitting techniques.
Here, we briefly compare our model to each of these others and show that our model (1) describes distinct combination surfaces, (2) results in synergy parameters which are straight forward to interpret, (3) is not restricted to a special class of effects with bounded scales, and (4) reduces to many of these other approaches in special cases thereby unifying and generalizing seemingly disparate synergy principles.

#### The Dose Equivalence Principle: Loewe and CI

The first prevalent foundational principle, established by Loewe (
• Loewe S.
• Muischnek H.
über Kombination swirkungen.
) and subsequently expanded on by CI (
• Chou T.C.
• Talalay P.
Quantitative analysis of dose-effect relationships: the combined effects of multiple drugs or enzyme inhibitors.
), is the Dose Equivalence Principle. This principle states that for a given effect magnitude E achieved by dose x of drug X alone or dose y of drug Y alone, there exists a constant ratio R = x/y such that using $Δx$ less of drug X can always be compensated for by using $Δy=RΔx$ more of drug Y. Therefore, the null surface is only defined for combinations whose magnitude of effect is less than the weaker drug’s maximal effect. This is because for combination effects greather than the effect of the weakest drug, no amount of the weaker drug can compensate for reducing the dose of the stronger drug.
The resulting null surfaces have linear isoboles. Our model recovers this under the constraint that the two drugs are maximally antagonistic. This can be seen by setting α=0, and reducing Equation 6 to
$(E−E0)+(E−E1)(d1Φ1)h1+(E−E2)(d2Φ2)h2=0.$

By this it is easy to see when h1 = h2 = 1, iso-effect lines $(∂∂E=0)$ are represented by the linear isoboles characteristic of Loewe Additivity and the CI null models. However, even in this case MuSyC is not limited by the weaker drug, and can therefore extend Loewe’s isoboles to any combination doses.
The requirement that α = 0 means the Loewe and CI null models assumes infinite potency antagonism (α1 = α2 = 0). Therefore, combinations with (0 <α <1) may be deemed synergistic by Loewe or CI. However, these values directly reflect a decrease in potency, and our formulation accurately identifies this as antagonistic. Finally, their null model also ignores the possible effect of hill slopes not equal to 1. For drugs with h < 1, they will tend to overestimate synergy, while drugs with h > 1 will lead to underestimated synergy (Figures 5C and 5D). Because their null model relies on such specific assumptions, which are not true for many drugs, it is generally impossible to know whether their results reflect true underlying synergy/antagonism, or simply stem from an inappropriate null surface.

#### The Multiplicative Survival Principle: Bliss and Effective Dose Model

The other prevalent foundational synergy principle is multiplicative survival, described by Bliss (
• Bliss C.I.
The toxicity of poisons applied jointly.
). Bliss’ null model assumes the probability of a cell being unaffected by drug 1 (U1) is independent of the probability of a cell being unaffected by drug 2 (U2). From this, the null surface states the probability of being unaffected by both drug 1 and drug 2 in combination is U1,2 = U1U2. When there is no potency synergy or antagonism, MuSyC reproduces this behavior in the following manner.
Setting α1 = α2 = 1, consider the fraction of unaffected cells, U, for each drug in isolation:
$Ui=11+(diΦi)hi.$

And for the two drugs in combination, solving Equation 67 for U we get
$U1,2=11+(d1Φ1)h1+(d2Φ2)h2+(d1Φ1)h1(d2Φ2)h2.$

From this, it is easy to verify that U1,2=U1U2, which is equivalent to Bliss Independence. However, the Bliss method explicitly requires the effect being measured in the combination surface is "percent affected", such as percent of cells killed vs. percent of cells remaining. For drugs which induce different maximum effects, Bliss is unable to account for the difference between being affected by drug 1 (E1), drug 2 (E2), and or both (E1,2), and may give unreliable results. Our model addresses this by decoupling the effect of a drug (E0,E1,E2,E3) and the "percent affected" by a drug (U,A1,A2,A1,2). If the effect itself is measuring percent (un)affected, that corresponds to the case where E0 = 1, E1 = E2 = E3 = 0, in which case MuSyC’s null model is identical to Bliss’.
Zimmer et. al. introduced the effective dose model (
• Zimmer A.
• Katzir I.
• Dekel E.
• Mayo A.E.
• Alon U.
Prediction of multidimensional drug dose responses based on measurements of drug pairs.
) as a parameterized version of Bliss, and shares the same null surface. However, while Bliss defines synergy at every concentration independently, the effective dose model introduces a parameter ai,j to quantify synergy, similar to MuSyC’s potency synergy (α). The ai,j parameter reflects how the presence of drug i modulates the potency of drug j. However, like Bliss, the effective dose model can only be applied to drug responses where the measured drug effect is "percent affected" thereby implicitly requiring the maximum effect of both drugs and the combination is 100% affected which is commonly not observed in dose-response studies (
• Fallahi-Sichani M.
• Heiser L.M.
• Gray J.W.
• Sorger P.K.
Metrics other than potency reveal systematic variation in responses to cancer drugs.
).

#### ZIP

Like the equivalent dose model (
• Zimmer A.
• Katzir I.
• Dekel E.
• Mayo A.E.
• Alon U.
Prediction of multidimensional drug dose responses based on measurements of drug pairs.
), as well as our potency synergy (α), ZIP (
• Wennerberg K.
• Aittokallio T.
• Tang J.
Searching for drug synergy in complex dose–response landscapes using an interaction potency model.
) works by quantifying how one drug shifts the potency of the other. ZIP is formulated for arbitrary E0 and Emax; however, it assumes Emax is the same for both drugs, as well as the combination (explicitly E1 = E2 = E3). To identify potency shifts, the ZIP method fixes the concentration of one drug, then fits a Hill-equation dose response for the other drug. However, for combinations with efficacy synergy or antagonism, dose responses can have non-Hill, and even non-monotonic shapes. In our data, several drugs displayed this behavior. Because our method accounts explicitly for efficacy synergy, our surfaces are able to describe such complex drug combination surfaces where ZIP fails.
Furthermore, ZIP calculates synergy at every concentration. This is similar to the approach taken by Bliss, Loewe, and CI, and can be used to find doses which “maximize” the observed synergy. However, quantifying synergy on a dose-by-dose basis confounds synergy of potency and efficacy which emerge only on inspection of the global dose-response surface. Additionally, this dose dependent synergy often leads to ambiguous results about whether a given combination is synergistic or not, as it synergizes at some concentrations, and antagonizes at others (Figure 5A).

#### BRAID

Like ZIP, BRAID (
• Twarog N.R.
• Stewart E.
• Hammill C.V.
• Shelat A.A.
BRAID: A unifying paradigm for the analysis of combined drug action.
) assumes that each drug alone has a sigmoidal dose-response, and constructs a Hill-like equation for the combination. This equation uses a single dose parameter κ which combines the doses of both individual drug. To uniquely solve for κ, this formalism, like Loewe additivity, adds the constraint that a drug in combination with itself must be neither synergistic nor antagonistic. By adjusting κ, BRAID is able to fit complex drug combination surfaces, including non-monotonic responses. Because BRAID fits the whole combination surface using a single parameter, it can be used to make unambiguous statements about whether the combination is synergistic or antagonistic. Nevertheless, BRAID does not account for differences in synergy due to efficacy vs. potency, whereas we find many combinations that are synergistic with respect to one, but antagonistic with respect to the other. Further, the biochemical interpretation of κ is not straightforward. And finally, the BRAID model is unable to fit combination surfaces with synergistic efficacy, as it assumes that the maximum effect of the combination is equal to the maximum effect of the stronger single-drug.

#### Highest Single Agent (HSA)

HSA, originally proposed by
Pharmacology.
and then revived later by Greco (
• Greco W.R.
• Bravo G.
• Parsons J.C.
The search for synergy: a critical review from a response surface perspective.
), is a simple heuristic which argues synergy is any combination effect which exceeds the effect of either single agent. While β is conceptually similar to HSA, β provides a global view of the possible increase in effect rather than a point-by-point dose comparison as done in HSA. Because HSA is calculated at every dose it cannot distinguish between synergistically efficacious combinations and synergistically potent combinations as both will increase the effect at intermediate doses (Figure S2). Additionally, as HSA is only defined on a dose-by-dose basis with no model fit, it is sensitive to the dose range selected.

#### Schindler 2D-Partial Differential Equation (PDE) Model

Schindler’s Hill PDE was derived to impute the dose-response surface from the single dose-response curves alone (
• Schindler M.
Theory of synergistic effects: hill-type response surfaces as “null-interaction” models for mixtures.
). Therefore, it does not contain any fit parameters, but rather defines a null surface for which synergy results in deviations in the surface. While Schindler did not specify how to account for these devations, he postulates some implementation of perturbation theory would be sufficient. Like CI and the Equivalent Dose Model, Schindler’s framework requires effects in a range between 0 and 1, based on the assumption that the metric is a percent. Therefore, Schindler cannot be applied to data collected with other types of metrics (e.g., DIP Rates). Additionally, Schindler’s maximum effect of the combination (E3) is set equal to the average of the single drug maximal effect. This allows for smooth transitions between the two single dose-response curves but results in some non-intuitive solutions. For example, if drug 1 has a maximal effect of 50% and drug 2 has a maximal effect of 70% the expected additive effect of the combination in the null model is 60% which is less than the maximal effect of drug 1. Therefore, an effect of 65% in combination, though less than achievable with one drug, is designated synergistic by Schindler.

#### Sham Experiment

It is common for synergy metrics to examine the special case in which the two drugs being combined are actually the same drug in a so called sham experiment first postulated by Chou (
• Chou T.C.
• Talalay P.
Quantitative analysis of dose-effect relationships: the combined effects of multiple drugs or enzyme inhibitors.
). Famously, Loewe, Combination Index, and other methods based on the Dose Equivalence Principle are sham compliant while Bliss and other methods based on the Multiplicative Survival Principle are not. Because our method distinguishes between two types of synergy, we tested sham compliance for each independently. It is immediately apparent synergistic efficacy is sham complaint in all conditions. This can be observed by substituting E1 = E2 = E3, as the maximum effect of the drug remains constant, into the definition for β in Equation 3
$β=min(E1,E2)−E3E0−min(E1,E2)=0.$
(Equation 7)

To test the sham compliance of synergistic potency, we can write the full dose response surface as a direct 2D extension of the 1D dose-response curve in Equation 12 by replacing d with d1 + d2.
$Ed=Em(d1+d2)h+E0Ch(d1+d2)h+Ch$
(Equation 8)

Our 2D generalization of Equation 12, given by Equation 1 can be rewritten for the case of 2 identical drugs by observing that C1 = C2 = C, h1 = h2 = h, and E1 = E2 = E3 = Em, resulting in
$Ed=C2hE0+d1hChEm+Chd2hEm+(α2d1d2)hEmC2h+d1hCh+Chd2h+(α2d1d2)h.$
(Equation 9)

Setting Equations 8 and 9 equal to one another, we find
$α2=Ch(d1+d2)h−d1h−d2h(d1d2)h.$
(Equation 10)

This equality is true when α2 = α1 = 0 and h = 1. This makes sense as our model reduces to Loewe additivity under those conditions, and Loewe additivity was developed to explicitly address the sham-combination case. In conclusion, MuSyC satisfies the sham experiment in all conditions where Loewe is the appropriate model.

#### One-Dimensional Sigmoidal Dose-Response Curve

In pharmacology, the effect of a drug is usually described by the Hill equation, which arises from the equilibrium of a reversible transformation between an unaffected population (U) and an affected population (A)
(Equation 11)

Here, [d] is the concentration of the drug, h is the Hill slope, and r1 and r−1 are constants corresponding to its rate of action. Solving for the equilibrium results in
$∂U∂t=A⋅r−1−U⋅r1dh≡0$

$AU=r1dhr−1.$

When $dh=r−1r1$, then half the population is affected, and half is unaffected (A = U). This dose is the EC50, denoted as $Ch=r−1r1$. Adding the constraint that U + A = 1, which states that 100% of the population is either unaffected or affected, we find the classic Hill equation:
$U=ChCh+dh.$

If the unaffected and affected populations differ phenotypically by some arbitrary effect (e.g., proliferation rate), the observed effect over the whole population at dose d of some drug will be a weighted average of the two effects by the percent affected and unaffected. Namely,
$Ed=U⋅E0+A⋅Em,$

where E0 is the effect characteristic of the unaffected population, and Em is the effect characteristic of the affected population. From this we find the final form of a 4-parameter sigmoidal equation describing dose-response due to Hill-kinetics:
$Ed−EmE0−Em=ChCh+dh$
(Equation 12)

#### Extending the Mass Action Paradigm to Simple Four-State Model Assuming Detailed Balance

Consider a cell type U that can transition into a “drugged” state A1 in the presence of drug d1 and into a different drugged state A2 in the presence of drug d2 (Figure S1B). We can write these transitions as
$U↔r−1r1∗[d1]A1$
(Equation 13)

$U↔r−2r2∗[d2]A2,$
(Equation 14)

where [di] denotes concentration of drug di. At equilibrium, the forward and reverse rates of these processes are equal, i.e.,
$r1[d1][U]=r−1[A1]$
(Equation 15)

$r2[d2][U]=r−2[A2],$
(Equation 16)

where [Ai] is the population of cell state Ai. Defining Θx as the ratio of forward and reverse rates ($Θx≡r−xrx$) and assuming the system obeys detail balance, we find
$[U][A1]=Θ1[d1]$
(Equation 17)

$[U][A2]=Θ2[d2].$
(Equation 18)

Now assume that a fourth state exists, A1,2, corresponding to a “doubly” drugged state (Figure S1B). A1 cells can transition into this state in the presence of drug d2 and A2 cells can transition into this state in the presence of drug d1. We can write these processes as
$A1↔r−2r2∗α1[d2]A1,2$
(Equation 19)

$A2↔r−1r1∗α2[d1]A1,2.$
(Equation 20)

Note that without loss of generality, we set the forward rate constant for (19) equal to the same value in (14) multiplied by a factor α1>0. Similarly, the rate constant for (20) is the same as in (13) multiplied by a factor α2>0. Here α represents how each drug potentiates the action of the other and can be interpreted as a change in the “effective” dose of one drug given the presence of the other. When α=1 the effective dose of the first drug is the same given the presence of the second drug. When α<1, more of the first drug is required to observe the same effective concentrations due to the presence of the second drug. Finally, when α>1 the same concentration of the first drug is effectively increased by the second drug.
Again asserting the system obeys detailed balance at equilibrium, we have
$[A1][A1,2]=1α1Θ2[d2]$
(Equation 21)

$[A2][A1,2]=1α2Θ1[d1]$
(Equation 22)

We can derive the relationship between the multiplicative factors α1 and α2 by rearranging Equation 17 as
$[U]=Θ1[d1][A1]$
(Equation 23)

Substituting for [A1] from Equation 21 gives
$[U]=1α1Θ1[d1]Θ2[d2][A1,2]$
(Equation 24)

Substituting for [A1,2] from Equation 22 gives
$[U]=α2α1Θ2[d2][A2]$
(Equation 25)

Finally, substituting for [A2] from Equation 18 gives
$[U]=α2α1[U]$
(Equation 26)

i.e., α1=α2 = α. Note this equality only holds for systems obeying detailed balance. In general, we do not assume this (See Section 4.4 ’Generalized derviation without assuming detailed balance’) and α1 and α2 are independent (Figure S4). However, assuming detailed balance facilitates in deriving a more intuitive form of the 2D Hill equation (Equation 1) compared to the full form (Equation 6).
Now, we define the total cell count
$CT≡[U]+[A1]+[A2]+[A1,2].$
(Equation 27)

Substituting for [A1], [A2], and [A1,2] from Equations 17, 18, and 24, respectively, gives
$CT=[U]+[d1]Θ1[U]+[d2]Θ2[U]+α[d1]Θ1[d2]Θ2[U].$
(Equation 28)

Solving for [U] gives
$[U]=Θ1Θ2CTΘ1Θ2+[d1]Θ2+Θ1[d2]+α[d1][d2].$
(Equation 29)

Substituting Equation 29 into Equation 17 and rearranging gives
$[A1]=[d1]Θ2CTΘ1Θ2+[d1]Θ2+Θ1[d2]+α[d1][d2].$
(Equation 30)

Similarly, from Equation 18 we get
$[A2]=Θ1[d2]CTΘ1Θ2+[d1]Θ2+Θ1[d2]+α[d1][d2],$
(Equation 31)

and from Equation 24
$[A1,2]=α[d1][d2]CTΘ1Θ2+[d1]Θ2+Θ1[d2]+α[d1][d2].$
(Equation 32)

As in the derivation of the 1D Hill equation, the measured effect ($Ed$) is then the relative proportion of cells in each state multiplied by the effect characteristic of that state as in
$Ed=E0∗U+E1∗A1+E2∗A2+E3∗A1,2.$
(Equation 33)

Here we define the effect of each state (E0,E1,E2,E3) to be proliferation rate in the following way.
We assume that cells in each state can divide and die at rates characteristic of the state, i.e.,
$Ci→kidivCi+Ci$
(Equation 34)

$Ci→kidie∅$
(Equation 35)

where Ci is specific state of the cell.
We define the drug-induced proliferation (DIP) rate for each state as the difference between the division and death rate constants, i.e.,
$kidip≡kidiv−kidie.$
(Equation 36)

Using Equation 27, the rate of change of the total cell population is
$dCTdt=d[U]dt+d[A1]dt+d[A2]dt+d[A1,2]dt.$
(Equation 37)

From (13), (14), (19), (20), (34)–(36), we get
$dCTdt=k0dip[U]+k1dip[A1]+k2dip[A2]+k3dip[A1,2].$
(Equation 38)

Substituting Equations 29, 30, 31, and 32 and rearranging, we get
$dCTdt=kTdipCT,$
(Equation 39)

with
$kTdip≡Θ1Θ2k0dip+[d1]Θ2k1dip+Θ1[d2]k2dip+α[d1][d2]k3dipΘ1Θ2+[d1]Θ2+Θ1[d2]+α[d1][d2].$
(Equation 40)

Note that with a slight modification, Equation 40 can be written as
$kTdip=Θ1k0dip+[d1]k1dip+Θ1[d2]Θ2k2dip+α[d1][d2]Θ2k3dipΘ1+[d1]+Θ1[d2]Θ2+α[d1][d2]Θ2.$
(Equation 41)

Therefore, if [d2]=0 (i.e., single-drug treatment) we get
$kTdip=Θ1k0dip+[d1]k1dipΘ1+[d1]$

$=Θ1k0dip+[d1]k1dip+(Θ1k1dip−Θ1k1dip)Θ1+[d1]$

$=(Θ1+[d1])k1dip+Θ1(k0dip−k1dip)Θ1+[d1]$

$=k1dip+Θ1Θ1+[d1](k0dip−k1dip).$
(Equation 42)

Rearranging gives
$kTdip−k1dipk0dip−k1dip=Θ1Θ1+[d1].$
(Equation 43)

Comparing to Equation 12, we see that Equation 43 is a one-dimensional sigmoidal dose-response curve with $Ed=kTdip$, $E0=k0dip$, $Em=k1dip$, $C=Θ1$, and h=1. By analogy, we surmise that Equation 40 is a 9 parameter, two-dimensional generalization of Equation 12, i.e.,
$Ed=C1h1C2h2E0+d1h1C2h2E1+C1h1d2h2E2+(α2d1)h1d2h2E3C1h1C2h2+d1h1C2h2+C1h1d2h2+(α2d1)h1d2h2,$
(Equation 44)

with $Ed=kTdip$, $E0=k0dip$, $E1=k1dip$, $E2=k2dip$, $E3=k3dip$, $C1=Θ1$, $C2=Θ2$, h1=1, h2=1, and the additional parameter α2. Note that under the assumption of detailed balance we found α1=α2 for the case when h=1. Therefore, in the general case when h≠1, $α1h2=α2h1$. By fitting the Equation 44, α1 is uniquely determined.

#### Four-state model with multiple steps between states

Let us assume instead of occurring in a single step, the cell state transitions are h step processes, i.e.,
$Ci↔r−x,1rx,1∗αy[dx]Cij1↔r−x,2rx,2∗αy[dx]Cij2⋯↔r−x,h−1rx,h−1∗αy[dx]Cijh−1↔r−x,hrx,h∗αy[dx]Cj.$
(Equation 45)

Assuming that all steps are in rapid equilibrium, it is straightforward to show that
$[Ci][Cj]=∏m=1hΘx,m[αydx]h,$
(Equation 46)

where $Θx,m≡r−x,m/rx,m$. Defining $Φx≡∏m=1hΘx,mh$, Equation 46 can be written as
$[Ci][Cj]=Φxh[αydx]h,$
(Equation 47)

which is the well-known Median-Effect Equation from Chou (
• Chou T.-C.
• Talalay P.
Analysis of combined drug effects: a new look at a very old problem.
,
• Chou T.C.
• Talalay P.
Quantitative analysis of dose-effect relationships: the combined effects of multiple drugs or enzyme inhibitors.
,
• Chou T.C.
Drug combination studies and their synergy quantification using the Chou-Talalay method.
). Replacing reactions (13) and (14) with multi-step processes of the form (45), gives us
$[U][A1]=Φ1h1[d1]h1$
(Equation 48)

$[U][A2]=Φ2h2[d2]h2.$
(Equation 49)

Similarly, we replace reactions (19) and (20) with the same multi-step process except with the rate constant for the $Ci→Ci1$ transition (entry into the cascade) equal to $αyrx,1[dx]$, giving
$[A1][A1,2]=Φ2h2[α1d2]h2$
(Equation 50)

$[A2][A1,2]=Φ1h1[α2d1]h1.$
(Equation 51)

Note, we assume that the number of steps in the cascade (45) is dependent on the drug type (i.e., $U→A1$ and $A2→A1,2$, both driven by d1, take h1 steps, while $U→A2$ and $A1→A1,2$, both driven by d2, take h2 steps). Using Equations 48, 49, 50, and 51 and again defining the total cell count CT as in Equation (27), we derive
$[U]=Φ1h1Φ2h2CTΦ1h1Φ2h2+[d1]h1Φ2h2+Φ1h1[d2]h2+[α2d1]h1[d2]h2$
(Equation 52)

$[A1]=[d1]h1Φ2h2CTΦ1h1Φ2h2+[d1]h1Φ2h2+Φ1h1[d2]h2+[α2d1]h1[d2]h2$
(Equation 53)

$[A2]=Φ1h1[d2]h2CTΦ1h1Φ2h2+[d1]h1Φ2h2+Φ1h1[d2]h2+[α2d1]h1[d2]h2$
(Equation 54)

$[A1,2]=α[d1]h1[d2]h2CTΦ1h1Φ2h2+[d1]h1Φ2h2+Φ1h1[d2]h2+[α2d1]h1[d2]h2.$
(Equation 55)

Therefore, in the same way that we arrived at Equation 40, we can derive
$kTdip≡Φ1h1Φ2h2k0dip+[d1]h1Φ2h2k1dip+Φ1h1[d2]h2k2dip+[α2d1]h1[d2]h2k3dipΦ1h1Φ2h2+[d1]h1Φ2h2+Φ1h1[d2]h2+[α2d1]h1[d2]h2,$
(Equation 56)

which is of the form Equation 44 with $Ed=kTdip$, $E0=k0dip$, $E1=k1dip$, $E2=k2dip$, $E3=k3dip$, $C1=Φ1$, and $C2=Φ2$. From this it is clear the hill coefficent (h) is related to the number of intermediate steps in the system.
The derivation of Equation 56 assumes that the populations of all intermediate cell states $Cijm$ ($m∈{1…h−1}$) in (45) are small ($≈0$). (This is most evident in our use of Equation 27 for the total cell population, where we only consider the end states. However, it is also implicit in our use of Equation 46, which is used to derive Equations 52, 53, 54, and 55 that lead to Equation 56 via Equation 38. In other words, we are assuming that the intermediate states do not significantly contribute to the dynamics of the total cell population. Since it is not reasonable to assume that cells in these states do not divide and die, we must assume the percent occupancy of these states is near zero.) We can satisfy this assumption by requiring that all $rx,m$, $r−x,m≫1$ $(m∈{1…h})$ and $Θx,1≫Θx,2≈…≈Θx,h−1≫Θx,h$. To see this, consider cell state U and all of its intermediate states between states A1 and A2. Let us define
$UT≡[U]+∑m=1h1−1[C01m]+∑m'=1h2−1[C02m'].$
(Equation 57)

From (45), we see that $Θx,1[dx]=[Ci][Cij1],Θx,2[dx]=[Cij1][Cij2],$ etc. Therefore,
$UT=[U](1+[d1]Θ1,1+[d1]2Θ1,1Θ1,2+…+[d1]h1−1∏m=1h1−1Θ1,m+[d2]Θ2,1+[d2]2Θ2,1Θ2,2+…+[d2]h2−1∏m'=1h2−1Θ2,m').$
(Equation 58)

If $Θ1,1≫1$, $(m∈{2…h1−1})$ and $Θ2,1≫1$, $(m'∈{2…h2−1})$, we get $UT≈[U]$, i.e., the populations of all intermediate states are ≈0. Now consider cell state A1 and all of its intermediate states between cell state A1,2. Similar to above, we have
$A1T≡[A1]+∑m=1h2−1[C13m]$
(Equation 59)

and
$A1T=[A1](1+α[d2]Θ2,1+[α1d2]2Θ2,1Θ2,2+…+[α1d2]h2−1∏m=1h2−1Θ2,m).$
(Equation 60)

Thus, as before, if $Θ2,1≫1$ and $(m∈{2…h2−1})$ we have $A1T≈[A1]$. However, from (45) we also have
$[A1]=[d1]Θ1,h1[C01h1−1]$

$=[d1]h1Θ1,1Θ1,2…Θ1,h1$
(Equation 61)

Since, from above, $Θ1,1≫1$ and $(m∈{2…h1−1})$, we must require that $Θ1,h1≪1$ and $(m∈{2…h1−1})$ in order to offset the large value of $Θ1,1$ and to ensure that , and to ensure. The latter condition means that $Θ1,m≈1$ $(m∈{2…h1−1})$. Therefore, we have the condition that $Θ1,1≫Θ1,2≈…≈Θ1,h1−1≫Θ1,h1$. Similarly, we can derive that $Θ2,1≫Θ2,2≈…≈Θ2,h2−1≫Θ2,h2$ by considering cell state A2 and all of its intermediate states between cell state A1,2 (not shown).

#### Generalized Derivation without Assuming Detailed Balance

More generally if we do not assume detailed balance, the state occupancy of U,A1,A2,A1,2 are defined by the partial equilibrium equations
$∂U∂t=−U⋅(r1d1+r2d2)+A1⋅r−1+A2⋅r−2$
(Equation 62)

$∂A1∂t=−A1⋅(r−1+α1r2d2)+U⋅r1d1+A1,2⋅r−2$
(Equation 63)

$∂A2∂t=−A2⋅(α2r1d1+r−2)+U⋅r2d2+A1,2⋅r−1$
(Equation 64)

$∂A1,2∂t=−A1,2⋅(r−1+r−2)+A1⋅α1r2d2+A2⋅α2r1d1.$
(Equation 65)

A final constraint is
$U+A1+A2+A1,2=CT.$
(Equation 66)

At equilibrium, the Equations 62, 63, 64, and 65 must be equal to zero; however, the system only defines a rank 3 matrix, necessitating Equation 66. Thus we find
$[−(r1d1+r2d2)r−1r−20r1d1−(r−1+r2(α1d2))0r−2r2d20−(r1(α2d1)+r−2)r−11111]⋅[UA1A2A1,2]=[000CT].$
(Equation 67)

Equations of the form
$Y⋅x→=b→,$

can be solved as
$x→=Y−1⋅b→.$

Thus we find the expected effect, as in the 1D case, is the weighted average of the characteristic effect of each state weighted by the state occupancy as governed by the 2D Hill Equation.
$E=[E0E1E2E3]⋅Y−1⋅[0001].$
(Equation 68)

This is derived assuming mass action from the reaction rules $U+d→A1+d$, $A1→U$. If instead we assume a multi-step transition as in section 4.3, we can simply replace the following in 68
$d1→d1h1$

$d2→d2h2$

$α2d1→(α2d1)h1$

$α1d2→(α1d2)h2,$

resulting in
$Ed=[E0E1E2E3]⋅[−(r1d1h1+r2d2h2)r−1r−20r1d1h1−(r−1+r2(α1d2)h2)0r−2r2d2h20−(r1(α2d1)h1+r−2)r−11111]−1⋅[0001].$
(Equation 69)

Equation 69 has the following twelve explicit parameters: r1, r−1, r2, r−2, E0, E1, E2, E3, h1, h2, α1 and α2. There is a relationship defined between a drug’s EC50 (C in our derivation), the transition rates (ri, r−i), and the hill slope (hi), given by $Cihi=r−iri$ | {i = 1 or 2}.

#### Combination Experiments Protocol

Experiments were conducted in the Vanderbilt High Throughput Screening Facility. Cells were seeded at approximately 300 cells per well in 384-well plates and allowed to adhere overnight. A preliminary image of each plate was taken approximately 8 hours after seeding to verify sufficient numbers of cells for each experiment. Images were taken on either the ImageXpress Micro XL (Molecular Devices) or CellaVista. The matrix of drug concentrations was prepared using a row-wise and column-wise serial 2X or 4X dilution in 384 well plates using a Bravo Liquid Handling System (Agilent) or manually in 96-well plates. See Table S2 for dose ranges tested. After allowing to adhere overnight, medium containing drugs and 5 nM Sytox Green (to detect dead cells) was added (time = 0) and replaced after 72 hours. Images were obtained at intervals ranging from every 4 to 8 h, depending on the experiment, for >120 hours. Cell counts were determined using custom-image segmentation software developed in Python using scikit-image package (
• van der Walt S.
• Schönberger J.L.
• Nunez-Iglesias J.
• Boulogne F.
• Warner J.D.
• Yager N.
• Gouillart E.
• Yu T.
scikit-image contributors
Scikit-image: image processing in Python.
) and run in parallel using RabbitMQ/Celery (http://www.celeryproject.org/).

#### RNA-seq of Melanoma Cell Lines

Total RNA was isolated from untreated SKMEL5 single-cell derived sublines, each in triplicate, using Trizol isolation method (Invitrogen) according to the manufacturer’s instructions. RNA samples were submitted to Vanderbilt VANTAGE Core services for quality check, where mRNA enrichment and cDNA library preparation were done with Illumina Tru-Seq stranded mRNA sample prep kit. Sequencing was done at Paired-End 75 bp on the Illumina HiSeq 3000. Reads were aligned to the GRCh38 human reference genome using HISAT2 (
• Kim D.
• Salzberg S.L.
HISAT: a fast spliced aligner with low memory requirements.
) and gene counts were obtained using featureCounts (
• Liao Y.
• Smyth G.K.
• Shi W.
featureCounts: an efficient general purpose program for assigning sequence reads to genomic features.
). All downstream analyses were performed in R (https://www.r-project.org) using the Bioconductor framework (https://www.bioconductor.org)

#### RT-qPCR Quantification of NOX5 Expression

Total RNA was extracted using Trizol isolation method (Invitrogen) according to the manufacturer’s instructions. cDNA synthesis was performed with QuantiTect Reverse Transcription Kit (Cat# 205311) from Qiagen. RT-qPCR was performed using the IQTM SYBR Green Supermix from BioRad (Cat# 1708880). Amplifications were performed in BioRad CFX96 TouchTM Real-Time PCR Detection System. All experiments were done at least in 3+ technical replicates. Log2 of the transcript expressions were normalized to SKMEL5 subline SC01. HPRT or 36B4 were used as housekeeping gene in all the experiments. Primers used in RT-qPCR are listed in Key Resources Table.

### Quantification and Statistical Analysis

#### Fitting Dose-Response Surfaces

We developed a fitting algorithm, implemented in Python, to fit the combinations experiments to the 2D Hill equation. The fitting is done in three steps, first estimates of the single dose-response parameters (C1,C2,h1,h2,E0,E1,E2) are extracted from fits to the single dose-response curves using the Pythonic implementation of a Levenburg-Marquart (LM) least squares optimization (scipy.optimize.curve_fit). The fit uncertainty (σ) is then the square root of the covariance matrix which is approximated as the inverse of the Hessian matrix (equal to JTJ in LM where J is the Jacobian) at the solution. In the second step, a Particle Swarm Optimizer (10,000 particles, 100 iterations) fits the full 2D Hill equation using the single parameter fits and uncertainties as initial values and bounds (±2σ). In the last step, the PSO optimized values are used to construct priors for a Metropolis-Hastings Monte-Carlo Markov Chain (MCMC) Optimization (Metropolis Hastings 10,000 iterations). Convergence is tested by checking all parameters' Geweke Z-score. If the Z-score range is (-2,2) over the sampling time frame, the optimization is considered to converge (Figures S3D and S3E). We found it necessary to use both the PSO and MCMC in order to fit a wide range of dose-response surfaces (Figure S3). To test the sensitivity of our fitting algorithm, we generated in silico data for 125 different dose-response surfaces at different data densities. The density of data tested were square matrixes of rank 5, 7, 10, 15, and 25. At each density 25 different dose-response surfaces were sampled across a 5X5 grid of log(α) and β values ranging from [-2,2] and from [-0.5,0.5], respectively. The parameters for E0, the single drug hill slope, EC50, and maximal effects were held constant at (0.3, 1, 10e-5, and 0.0), respectively. Random noise equal to the average uncertainty in the DIP Rate fits from the NSCLC screen was added to the data (0.001). In all conditions we observed a PSO particle count of 10,000 converged to a minimum in <60 iterations (Figure S3A). However, this minimum was not the optimal solution. The addition of an MCMC walk further improved fits (Figure S3B) (Pymc3 package). The MCMC walk calculates the posterior distribution for each parameter from which each parameter’s value (mean of trace) and uncertainty (standard deviation of trace) is calculated. Uncertainty in (E0, E1, E2, and E3) was propagated when calculating β using the equation

where Ex and σEx are min(E1,E2) and σ(min(E1,E2)) respectively. All other uncertainty propagations were handled with python package uncertainties (
• Lebigot E.O.
Uncertainties: a Python package for calculations with uncertainties.
). By calculating the uncertainty in the synergy parameters from the posterior distributions, the significance of synergy can be assessed in an unbiased way. Multiple factors contribute to increasing uncertainty in the fitted parameters. Dose-selection, an important consideration in all drug response profiling, changes the certainty of the fits (Figures S3C–S3E). While we are unable to observe saturating effects implicit in the model for some of our drug combinations – due to limited solubility or potency of the drug – by keeping careful account of the uncertainty in our synergy calculations we can still interpret the synergy of non-optimal dose-regimes. To demonstrate this, we generated the same 25 dose-response surfaces with varying log(α) and β values ranging from [-2,2] and from [-0.5,0.5] respectively but at different coverage of the dose-response curve. The uncertainty in the synergy parameters increases for decreased dose range (Figures S3C and S3D). It is important to note that in general the uncertainty is a function of many different aspects other than data density including the hill slope of the single curves (high hill slopes can result in higher uncertainty in log(α)), noise of experimental data, and quality of priors resulting from the single-drug fits. We posit the rigorous approach taken here accounts for all these sources resulting in a true estimate of confidence in a particular synergy value.
To prevent over fitting the data, we have defined six different model tiers which have increasing degrees of freedom (Table S5). To select the correct model tier, we penalize models with higher degrees of freedom by selecting the model based on minimizing the deviance information criterion (DIC) (
• Berg A.
• Meyer R.
• Yu J.
Deviance information criterion for comparing stochastic volatility models.
). Fits for each nest are used to inform priors for subsequent nests. Only drug combinations which converged to the full model (tier 5 with fits for all 12 parameters – Equation 69 in section 4.4 of Methods Details) were used for subsequent analysis. The MCMC optimization additionally allows for quantification of parameter confidence given the data. The following packages were used for fitting, data analysis, or visualization: GNU parallel (
• Tange O.
GNU parallel: the command-line power tool.
), SciPy (
• Jones E.
• Oliphant T.
• Peterson P.
SciPy: open source scientific tools for python.
), Numpy (
• Oliphant T.E.
Guide to NumPy.
), Pandas (

McKinney, W. (2010) Data structures for statistical computing in python, Proceedings of the 9th Python in Science Conference, pp. 51–56.

), Matplotlib (
• Hunter J.D.
Matplotlib: a 2D graphics environment.
). Pymc3 (
• Salvatier J.
• Wiecki T.V.
• Fonnesbeck C.
Probabilistic programming in Python using PyMC3.
).

#### Calculating the DIP Rate

Traditionally, the efficacy of an anti-proliferative compound is measured as the percent of viable cells (relative to control) after a treatment interval (
• Fallahi-Sichani M.
• Heiser L.M.
• Gray J.W.
• Sorger P.K.
Metrics other than potency reveal systematic variation in responses to cancer drugs.
); however, it has been recently shown this metric is subject to temporal biases (
• Hafner M.
• Niepel M.
• Chung M.
• Sorger P.K.
Growth rate inhibition metrics correct for confounders in measuring sensitivity to cancer drugs.
;
• Harris L.A.
• Frick P.L.
• Garbett S.P.
• Hardeman K.N.
• Paudel B.B.
• Lopez C.F.
• Quaranta V.
• Tyson D.R.
An unbiased metric of antiproliferative drug effect in vitro.
). To address these biases, we previously developed an unbiased metric of drug effect termed the drug-induced proliferation (DIP) rate (
• Harris L.A.
• Frick P.L.
• Garbett S.P.
• Hardeman K.N.
• Paudel B.B.
• Lopez C.F.
• Quaranta V.
• Tyson D.R.
An unbiased metric of antiproliferative drug effect in vitro.
). The DIP rate is defined as the steady state proliferation rate after drug equilibration. A positive DIP rate indicates an exponentially growing population, while a negative DIP rate indicates a regressing one. A rate of zero indicates a cytostatic effect, which may result from cells entering a non-dividing state or from balanced death and division (
• Paudel B.P.
• Harris L.A.
• Hardeman K.N.
• Abugable A.A.
• Hayford C.E.
• Tyson D.R.
• Quaranta V.
A nonquiescent “idling” population state in drug-treated, BRAF-mutated melanoma.
). We used the available findDIP R package for calculating DIP rates from growth curves which automatically selects the interval after drug equilibration (https://github.com/QuLab-VU/DIP_rate_NatMeth2016.git).

#### Calculating Loewe, CI, Bliss, and HSA

To compare our method to the prevailing methods for computing synergy, we calculated Loewe, CI, Bliss for the data from the osimertinib screen in Figure 2 and melanoma BRAF/MEK data of Figure 3. Loewe is agnostic to effect metric, and so we applied it directly to the DIP rate. To calculate CI and Bliss, we imputed the percent viability at 72 hours from the DIP rate for each condition. Percent viability is defined as in equation 1.
(Equation 70)

Estimates of percent viability are sensitive to even small differences between initial cell counts in the control and treated wells due to exponential amplification (
• Harris L.A.
• Frick P.L.
• Garbett S.P.
• Hardeman K.N.
• Paudel B.B.
• Lopez C.F.
• Quaranta V.
• Tyson D.R.
An unbiased metric of antiproliferative drug effect in vitro.
). To correct for this the bias, a 'matching' control cell count at 72-hours for each treated condition was calculated using equation
(Equation 71)

where Control Growth Rate is the median of the fitted growth rates for all control wells. Because the automated microscope did not image all the conditions at exactly zero or seventy two hours, we extrapolate and interpolate respectively the cell count at these times from the measured time series.
The Bliss metric only requires marginal data. For each experiment, individually, we calculated a Bliss score as
$Bliss=PV1|d1∗PV2|d2−PV1,2|d1,d2$
(Equation 72)

where $PVi|di$ is the %-viability measured for treatment with drug i alone at dose di, and $PV1,2|d1,d2$ is the %-viability measured for the combination treatment. The first term corresponds to the expected viability, assuming independence, while the second term is the measured viability. By this definition, Bliss>0 is synergistic, and Bliss<0 is antagonistic.
Loewe and CI require parameterization of a 1D Hill equation for each drug alone.
$E=Emax+E0−Emax(dC)h+1.$
(Equation 73)

CI, as per standard calculations (
• Chou T.C.
• Talalay P.
Quantitative analysis of dose-effect relationships: the combined effects of multiple drugs or enzyme inhibitors.
), further requires that E0 = 1 and Emax = 0 and is fit to a linearized, log-transformed version of the hill equation(
• Chou T.C.
Drug combination studies and their synergy quantification using the Chou-Talalay method.
) which has been previously critiqued for artificial compression of uncertainty in experimental data leading to poor model fits compared with nonlinear regression (
• Ashton J.C.
Drug combination studies and their synergy quantification using the Chou–Talalay method—letter.
). CI dose-response curves were fit using the scipy.stats.linregress module. All data points with percent viability greater than 1 were excluded from the CI fit, as log(1−pervia/pervia) becomes complex. For some drugs, this left too few points to fit a line, such that CI was undefined for combinations with those drugs. In other drugs, the fit hill coefficient was negative, and likewise all CI values were undefined for those drugs. For Loewe, we used the single-drug parameters fit by MuSyC.
From these parameterized hill equations, Loewe and CI were calculated using
(Equation 74)

where Di is the amount of drug i which, alone, achieves an effect equal to the combination effect, and is calculated from the Hill equation fit for that drug. We take the negative log to transform the synergy values to match Bliss, such that S>0 is synergistic, while S<0 is antagonistic. Because Loewe allows the two drugs to have different Emax, Loewe synergy cannot be calculated for measurements which exceed the weaker drug's Emax because no amount of the weaker drug alone would be sufficient to achieve that effect; therefore, those conditions are undefined.
For calculating HSA (
Pharmacology.
), we calculate the difference between the observed effect at each combination concentration and the most efficacious single agent effect at those doses. This difference is integrated across the surface to yield a single value for a particular combination.

#### Fitting ZIP, BRAID, Schindler’s Hill PDE, and Equivalent Dose Models

Theoretical dose-response surfaces with different combinations of synergistic potency and efficacy were generated then fit to estimate the synergy according to these methods (Figure S2). Both ZIP and BRAID were calculated using the R packages available for each method (ZIP’s R code is in the supplemental file 1 of the manuscript (
• Wennerberg K.
• Aittokallio T.
• Tang J.
Searching for drug synergy in complex dose–response landscapes using an interaction potency model.
) and BRAID’s package is available from:https://cran.r-project.org/web/packages/braidReports/braidReports.pdf). Schindler’s Hill PDE model contains no fitting parameters as the dose-response surface is derived purely from the marginal data. In fact, Schindler does not propose a method to estimate synergy from experimental data, but postulates some implementation of perturbation theory could be used to fit experimental data (
• Schindler M.
Theory of synergistic effects: hill-type response surfaces as “null-interaction” models for mixtures.
). Therefore, to calculate the synergy of this model, we defined the sum of residuals between the null surface and the experimental data to the metric of synergy. Finally, to fit Zimmer et al.’s Equivalent Dose Model we used the curve_fit() module of the scipy.optimization package in python. Specifically, the Equivalent Dose Model, equation 2 in (
• Zimmer A.
• Katzir I.
• Dekel E.
• Mayo A.E.
• Alon U.
Prediction of multidimensional drug dose responses based on measurements of drug pairs.
), contains parameters for C1,C2, a12, a21, h1, and h2 where the C parameters are the EC50 of the single agents, the ai,j parameters are the synergy values corresponding to a change in potency, and the h parameters are the hill slopes of the single agents. In the model, there are no parameters for efficacy because it is assumed the drug effect ranges between zero and one. When this is not true, the Effective Dose Model results in poor fits to the data (Figure S2) similar to CI.

#### Identifying DEGs for GO Enrichment Analysis

Differentially Expressed Genes (DEGs) were selected by ANOVA on baseline gene expression data on three clones based on a statistical cutoff of Likelihood Ratio Test (LRT) (p-values < 0.001). Functional enrichment analyses, including GO Term Enrichment and Pathway Enrichment Analysis were done using CRAN Package “Enrichr” (https://cran.r-project.org/web/packages/enrichR/index.html), based on a web-based tool for analyzing gene sets and enrichment of common annotated biological functions (
• Kuleshov M.V.
• Jones M.R.
• Rouillard A.D.
• Fernandez N.F.
• Duan Q.
• Wang Z.
• Koplev S.
• Jenkins S.L.
• Jagodnik K.M.
• Lachmann A.
• et al.
Enrichr: a comprehensive gene set enrichment analysis web server 2016 update.
). The enriched GO terms and enriched KEGG pathways were restricted to those with p-values corrected for multiple testing less than 0.001. The top GO Biological Processes included generation of precursor metabolites and energy, electron transport chain, inorganic cation transmembrane transport, and metabolic process. The top GO Molecular Function terms included inorganic cation transmembrane transporter activity, cofactor binding, NAD binding, and ATPase activity. The top GO Cellular Component term was the mitochondria membrane. Top KEGG pathways enriched in the DEGs included metabolic pathways, oxidative phosphorylation, carbon metabolism and TCA cycle (Figure 3B). Overall, these enriched GO terms and pathways point toward differences in the regulators of metabolic function in the three subclones. This is consistent with previous reports that suggest altered metabolism is implicated in drug sensitivity and melanoma resistance to BRAFi (
• Parmenter T.J.
• Kleinschmidt M.
• Kinross K.M.
• Bond S.T.
• Li J.
• Rao A.
• Sheppard K.E.
• Hugo W.
• Pupo G.M.
• et al.
Response of BRAF-mutant melanoma to BRAF inhibition is mediated by a network of transcriptional regulators of glycolysis.
,
• Hardeman K.N.
• Peng C.
• Paudel B.B.
• Meyer C.T.
• Luong T.
• Tyson D.R.
• Young J.D.
• Quaranta V.
• Fessel J.P.
Dependence on glycolysis sensitizes BRAF-mutated melanomas for increased response to targeted BRAF inhibition.
).
Correlation of BRAFi insensitivity was computed for each identified DEG according to DIP Rate at 8uM PLX-4720 for a 10 cell line panel (Table S3) Pair-wise comparisons of DEGs was performed on genes (after low count genes were removed) using DESeq2 pipeline (
• Love M.I.
• Huber W.
• Anders S.
Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.
).

### Data and Software Availability

All raw cell counts, calculated DIP rates, DEGs between subclones, and expression data are available in the github repo: https://github.com/QuLab-VU/MuSyC_Cell.git in the folder Data. Additionally, the repo contains all the code required to reproduce all the figures and supplemental figures from the data and is found in the Code_Paper_Figures folder. The subfolders Fig2 and Fig3 contain html folders with interactive plots of all the screened combinations. Open the .html files using a browser. The raw RNAseq is available from GEO at the accession number GSE122041. The software for interactive manipulation of the different parameters to study their contribution to the contours of the dose-response surface is also available in the github repo in the folder MuSyC_App. This folder contains both the matlab source code and a compiled application for the different operating systems. A copy of the github repo at the time of publication is also available from Mendeley Data via the following https://doi.org/10.17632/n8bp8db5ff.

## Acknowledgments

The authors would like to thank Corey Hayford, Sarah Maddox, Carlos Lopez, and Chris Wright for insightful conversations; Jing Hao for the help with experimental preparation and reagent procurement; and James Pino and Alexander Lubbock for the pythonic implementation of PSO and findDIP libraries, respectively.
Experiments were performed in the Vanderbilt HTS Core Facility, an institutionally supported core with assistance provided by Debbie Mi, Page Vinson, and Joshua Bower.
This work was supported by the following funding sources: C.T.M. was supported by the National Science Foundation (NSF) Graduate Student Fellowship Program (GRFP) [Award #1445197]. V.Q. was supported by the National Institutes of Health (NIH) (U54-CA217450, R01-186193, and U01-CA215845); J.B. was supported by the National Cancer Institute (NCI) (R50-CA211206). D.W. was supported by Ruth L. Kirschstein National Research Service Award (2T32HL094296-06A1). C.M.L. was supported in part by the NIH (P30-CA086485, R01-CA121210, and P01-CA129243). C.M.L. additionally was supported by a V Foundation Scholar-in-Training Award, an AACR-Genentech Career Development Award, a Damon Runyon Clinical Investigator Award, a LUNGevity Career Development Award, and a Lung Cancer Foundation of America/International Association for the Study of Lung Cancer Lori Monroe Scholarship.

### Author Contributions

Conceptualization, C.T.M., D.J.W., L.A.H., D.R.T., and V.Q.; Methodology, C.T.M., D.R.T., and J.B.; Software, C.T.M. and D.J.W.; Investigation, C.T.M., D.W., B.B.P., K.N.H., D.R.T., L.A.H., and J.B.; Formal Analysis, C.T.M., D.J.W., B.B.P., and D.R.T.; Resources, D.R.T., J.B., and D.W.; Data Curation, C.T.M. and D.R.T.; Writing – Original Draft, C.T.M., D.J.W., L.A.H., and V.Q.; Writing – Review & Editing, C.T.M., D.J.W., B.B.P., D.R.T., C.M.L., D.W., J.B., L.A.H., and V.Q.; Visualization, C.T.M., D.J.W, and B.B.P.; Supervision, L.A.H., D.R.T., and V.Q.; Funding Acquisition, D.R.T. and V.Q.

### Declaration of Interests

The authors declare no competing interests.

## Supplemental Information

• Document S1. Figures S1–S6 and Tables S1–S5

## References

• Altenhöfer S.
• Kleikers P.W.
• Wingler K.
• Schmidt H.H.
Evolution of NADPH oxidase inhibitors: selectivity and mechanisms for target engagement.
Antioxid. Redox Signal. 2015; 23: 406-427
• Ashton J.C.
Drug combination studies and their synergy quantification using the Chou–Talalay method—letter.
Cancer Res. 2015; 75: 2400
• Barretina J.
• Caponigro G.
• Stransky N.
• Venkatesan K.
• Margolin A.A.
• Kim S.
• Wilson C.J.
• Lehár J.
• Kryukov G.V.
• Sonkin D.
• et al.
The Cancer Cell Line Encyclopedia enables predictive modelling of anticancer drug sensitivity.
Nature. 2012; 483: 603-607
• Berg A.
• Meyer R.
• Yu J.
Deviance information criterion for comparing stochastic volatility models.
J. Bus. Econ. Stat. 2004; 22: 107-120
• Bliss C.I.
The toxicity of poisons applied jointly.
Ann. Appl. Biol. 1939; 26: 585-615
• Chou T.C.
Drug combination studies and their synergy quantification using the Chou-Talalay method.
Cancer Res. 2010; 70: 440-446
• Chou T.-C.
• Talalay P.
Analysis of combined drug effects: a new look at a very old problem.
Trends Pharmacol. Sci. 1983; 4: 450-454
• Chou T.C.
• Talalay P.
Quantitative analysis of dose-effect relationships: the combined effects of multiple drugs or enzyme inhibitors.
Adv. Enzyme Regul. 1984; 22: 27-55
• Eroglu Z.
• Ribas A.
Combination therapy with BRAF and MEK inhibitors for melanoma: latest evidence and place in therapy.
Ther. Adv. Med. Oncol. 2016; 8: 48-56
• Fallahi-Sichani M.
• Heiser L.M.
• Gray J.W.
• Sorger P.K.
Metrics other than potency reveal systematic variation in responses to cancer drugs.
Nat. Chem. Biol. 2013; 9: 708-714
• Foucquier J.
• Guedj M.
Analysis of drug combinations: current methodological landscape.
Pharmacol. Res. Perspect. 2015; 3: e00149
Pharmacology.
Oxford University Press, 1940
• Gong Z.
• Hu G.
• Li Q.
• Liu Z.
• Wang F.
• Zhang X.
• Xiong J.
• Li P.
• Xu Y.
• Ma R.
• et al.
Compound libraries: recent advances and their applications in drug discovery.
Curr. Drug Discov. Technol. 2017; 14: 216-228
• Greco W.
• Unkelbach H.-D.
• Pöch G.
• Sühnel J.
• Kundi M.
• Bödeker W.
Consensus on concepts and terminology for combined-action assessment: the Saariselka Agreement.
Arch. Complex Environ. Stud. 1992; 4: 65-69
• Greco W.R.
• Bravo G.
• Parsons J.C.
The search for synergy: a critical review from a response surface perspective.
Pharmacol. Rev. 1995; 47: 331-385
• Hafner M.
• Niepel M.
• Chung M.
• Sorger P.K.
Growth rate inhibition metrics correct for confounders in measuring sensitivity to cancer drugs.
Nat. Methods. 2016; 13: 521-527
• Hardeman K.N.
• Peng C.
• Paudel B.B.
• Meyer C.T.
• Luong T.
• Tyson D.R.
• Young J.D.
• Quaranta V.
• Fessel J.P.
Dependence on glycolysis sensitizes BRAF-mutated melanomas for increased response to targeted BRAF inhibition.
Sci. Rep. 2017; 7: 42604
• Harris L.A.
• Frick P.L.
• Garbett S.P.
• Hardeman K.N.
• Paudel B.B.
• Lopez C.F.
• Quaranta V.
• Tyson D.R.
An unbiased metric of antiproliferative drug effect in vitro.
Nat. Methods. 2016; 13: 497-500
• He B.
• Lu C.
• Zheng G.
• He X.
• Wang M.
• Chen G.
• Zhang G.
• Lu A.
Combination therapeutics in complex diseases.
J. Cell. Mol. Med. 2016; 20: 2231-2240
• Hennessey V.G.
• Rosner G.L.
• Bast Jr., R.C.
• Chen M.
A Bayesian approach to dose-response assessment and synergy and its application to in vitro dose-response studies.
Biometrics. 2010; 66: 1275-1283
• Hunter J.D.
Matplotlib: a 2D graphics environment.
Comput. Sci. Eng. 2007; 9: 90-95
• Jabs J.
• Zickgraf F.M.
• Park J.
• Wagner S.
• Jiang X.
• Jechow K.
• Kleinheinz K.
• Toprak U.H.
• Schneider M.A.
• Meister M.
• et al.
Screening drug effects in patient-derived cancer cells links organoid responses to genome alterations.
Mol. Syst. Biol. 2017; 13: 955
• Jaquet V.
• Marcoux J.
• Forest E.
• Leidal K.G.
• McCormick S.
• Westermaier Y.
• Perozzo R.
• Plastre O.
• Fioraso-Cartier L.
• Diebold B.
• et al.
NADPH oxidase (NOX) isoforms are inhibited by celastrol with a dual mode of action.
Br. J. Pharmacol. 2011; 164: 507-520
• Jia P.
• Jin H.
• Xia J.
• Ohashi K.
• Liu L.
• Pirazzoli V.
• Dahlman K.B.
• Politi K.
• Michor F.
• et al.
Next-generation sequencing of paired tyrosine kinase inhibitor-sensitive and -resistant EGFR mutant lung cancer cell lines identifies spectrum of DNA changes associated with drug resistance.
Genome Res. 2013; 23: 1434-1445
• Jones E.
• Oliphant T.
• Peterson P.
SciPy: open source scientific tools for python.
• Kim D.
• Salzberg S.L.
HISAT: a fast spliced aligner with low memory requirements.
Nat. Methods. 2015; 12: 357-360
• Kuleshov M.V.
• Jones M.R.
• Rouillard A.D.
• Fernandez N.F.
• Duan Q.
• Wang Z.
• Koplev S.
• Jenkins S.L.
• Jagodnik K.M.
• Lachmann A.
• et al.
Enrichr: a comprehensive gene set enrichment analysis web server 2016 update.
Nucleic Acids Res. 2016; 44: W90-W97
• Lebigot E.O.
Uncertainties: a Python package for calculations with uncertainties.
• Liao Y.
• Smyth G.K.
• Shi W.
featureCounts: an efficient general purpose program for assigning sequence reads to genomic features.
Bioinformatics. 2014; 30: 923-930
• Loewe S.
• Muischnek H.
über Kombination swirkungen.
Archiv. Experiment. Pathol. Pharmakol. 1926; 114: 313-326
• Loewe S.
Versuch einer allgemeinen Pharmakologie der Arznei- Kombinationen.
Klin. Wochenschr. 1927; 6: 1078-1085
• Long G.V.
• Stroyakovskiy D.
• Gogas H.
• Levchenko E.
• de Braud F.
• Larkin J.
• Garbe C.
• Jouary T.
• Hauschild A.
• Grob J.J.
• et al.
Combined BRAF and MEK inhibition versus BRAF inhibition alone in melanoma.
N. Engl. J. Med. 2014; 371: 1877-1888
• Love M.I.
• Huber W.
• Anders S.
Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.
Genome Biol. 2014; 15: 550
• Lu W.
• Hu Y.
• Chen G.
• Chen Z.
• Zhang H.
• Wang F.
• Feng L.
• Pelicano H.
• Wang H.
• Keating M.J.
• et al.
Novel role of NOX in supporting aerobic glycolysis in cancer cells with mitochondrial dysfunction and as a potential target for cancer therapy.
PLoS Biol. 2012; 10: e1001326
1. McKinney, W. (2010) Data structures for statistical computing in python, Proceedings of the 9th Python in Science Conference, pp. 51–56.

• Oliphant T.E.
Guide to NumPy.
CreateSpace Independent Publishing Platform, 2006
• Palmer A.C.
• Sorger P.K.
Combination.
Cell. 2017; 171: 1678-1691.e13
• Parmenter T.J.
• Kleinschmidt M.
• Kinross K.M.
• Bond S.T.
• Li J.
• Rao A.
• Sheppard K.E.
• Hugo W.
• Pupo G.M.
• et al.
Response of BRAF-mutant melanoma to BRAF inhibition is mediated by a network of transcriptional regulators of glycolysis.
Cancer Discov. 2014; 4: 423-433
• Paudel B.P.
• Harris L.A.
• Hardeman K.N.
• Abugable A.A.
• Hayford C.E.
• Tyson D.R.
• Quaranta V.
A nonquiescent “idling” population state in drug-treated, BRAF-mutated melanoma.
Biophys. J. 2018; 114: 1499-1511
• Salvatier J.
• Wiecki T.V.
• Fonnesbeck C.
Probabilistic programming in Python using PyMC3.
PeerJ Comput. Sci. 2016; 2: e55
• Schiffmann I.
• Greve G.
• Jung M.
• Lübbert M.
Epigenetic therapy approaches in non-small cell lung cancer: update and perspectives.
Epigenetics. 2016; 11: 858-870
• Schindler M.
Theory of synergistic effects: hill-type response surfaces as “null-interaction” models for mixtures.
Theor. Biol. Med. Modell. 2017; 14: 15
• Shaw A.T.
• Kim D.W.
• Mehra R.
• Tan D.S.
• Felip E.
• Chow L.Q.
• Camidge D.R.
• Vansteenkiste J.
• Sharma S.
• De Pas T.
• et al.
Ceritinib in ALK -Rearranged non–small-cell lung cancer.
N. Engl. J. Med. 2014; 370: 1189-1197
• Shimoyama Y.
• Nagafuchi A.
• Fujita S.
• Gotoh M.
• Takeichi M.
• Tsukita S.
• Hirohashi S.
Cadherin dysfunction in a human cancer cell line: possible involvement of loss of alpha-catenin expression in reduced cell-cell adhesiveness.
Cancer Res. 1992; 52: 5770-5774
• Soria J.C.
• Ohe Y.
• Vansteenkiste J.
• Reungwetwattana T.
• Lee K.H.
• Dechaphunkul A.
• Imamura F.
• Nogami N.
• Kurata T.
• et al.
Osimertinib in untreated EGFR -mutated advanced non–small-cell lung cancer.
N. Engl. J. Med. 2018; 378: 113-125
• Subramanian A.
• Narayan R.
• Corsello S.M.
• Peck D.D.
• Natoli T.E.
• Lu X.
• Gould J.
• Davis J.F.
• Tubelli A.A.
• Asiedu J.K.
• et al.
A next generation connectivity map: L1000 platform and the first 1,000,000 profiles.
Cell. 2017; 171: 1437-1452.e17
• Tallarida R.J.
Quantitative methods for assessing drug synergism.
Genes Cancer. 2011; 2: 1003-1008
• Tange O.
GNU parallel: the command-line power tool.
• Twarog N.R.
• Stewart E.
• Hammill C.V.
• Shelat A.A.
BRAID: A unifying paradigm for the analysis of combined drug action.
Sci. Rep. 2016; 6: 25523
• Tyson D.R.
• Garbett S.P.
• Frick P.L.
• Quaranta V.
Fractional proliferation: a method to deconvolve cell population dynamics from single-cell data.
Nat. Methods. 2012; 9: 923-928
• van der Walt S.
• Schönberger J.L.
• Nunez-Iglesias J.
• Boulogne F.
• Warner J.D.
• Yager N.
• Gouillart E.
• Yu T.
• scikit-image contributors
Scikit-image: image processing in Python.
PeerJ. 2014; 2: e453
• Welm B.E.
• Dijkgraaf G.J.
• Bledau A.S.
• Welm A.L.
• Werb Z.
Lentiviral transduction of mammary stem cells for analysis of gene function during development and cancer.
Cell Stem Cell. 2008; 2: 90-102
• Witta S.E.
• Gemmill R.M.
• Hirsch F.R.
• Coldren C.D.
• Hedman K.
• Ravdel L.
• Helfrich B.
• Chan D.C.
• Sugita M.
• et al.
Restoring E-cadherin expression increases sensitivity to epidermal growth factor receptor inhibitors in lung cancer cell lines.
Cancer Res. 2006; 66: 944-950
• Witta S.E.
• Jotte R.M.
• Konduri K.
• Neubauer M.A.
• Spira A.I.
• Ruxer R.L.
• Varella-Garcia M.
• Bunn P.A.
• Hirsch F.R.
Randomized phase II trial of erlotinib with and without entinostat in patients with advanced non-small-cell lung cancer who progressed on prior chemotherapy.
J. Clin. Oncol. 2012; 30: 2248-2255