## Abstract

Drug combination discovery depends on reliable synergy metrics but no consensus exists on the correct synergy criterion to characterize combined interactions. The fragmented state of the field confounds analysis, impedes reproducibility, and delays clinical translation of potential combination treatments. Here we present a mass-action based formalism to quantify synergy. With this formalism, we clarify the relationship between the dominant drug synergy principles, and present a mapping of commonly used frameworks onto a unified synergy landscape. From this, we show how biases emerge due to intrinsic assumptions which hinder their broad applicability and impact the interpretation of synergy in discovery efforts. Specifically, we describe how traditional metrics mask consequential synergistic interactions, and contain biases dependent on the Hill-slope and maximal effect of single-drugs. We show how these biases systematically impact synergy classification in large combination screens, potentially misleading discovery efforts. Thus the proposed formalism can provide a consistent, unbiased interpretation of drug synergy, and accelerate the translatability of synergy studies.

## Introduction

Throughout the preceding century, two principles have been used to quantify synergy of drug combinations: the Dose Equivalence Principle (DEP), introduced by Loewe^{1,2}, and the Multiplicative Survival Principle (MSP), introduced by Bliss^{3}. In 1992, a committee was convened in Saariselkä, Finland seeking to find a consensus between these principles and unify the field^{4,5}. Unable to reconcile their differences, the committee’s conclusion (The Saariselkä Agreement) did not reach a consensus and simply recommended that drug combination studies explicitly state how synergy was calculated^{4,5}. Multiple synergy models have since emerged, often derived as extensions of either the DEP or MSP, further splintering the field^{6,7,8,9,10,11,12}. In the absence of a consensus framework for drug synergy, discovery efforts for combinations often calculate all available synergy metrics^{13,14,15}, as first recommended by Greco and colleagues following Saariselkä^{16}. However, there remains no basis for choosing one metric over another, which becomes particularly problematic when synergy metrics conflict. This “calculate everything” paradigm thus hampers reproducibility between studies, delays progress in the discovery of synergistic drug combinations, and negatively impacts the translatability of combination discovery efforts.

Despite the lack of consensus on how to quantify synergy, drug combination screens remain essential to both pharmaceutical and academic discovery efforts, as shown by recent studies by AstraZeneca and the NCI-DREAM consortia^{17,18}, as well as combinatorial CRISPR screens^{19}. Yet, the paucity of successful clinical combinations explicable by true pharmacological interaction, rather than patient-to-patient variability^{20}, is symptomatic of the challenges facing the field. Therefore, the need identified at Saariselkä still exists: a consensus framework to interpret drug combination pharmacology.

We recently introduced a framework to quantify synergy based on the Law of Mass Action, named Multi-dimensional Synergy of Combinations (MuSyC)^{21}, that distinguishes between different synergy types (e.g., potency, efficacy). In the present work, we build upon our previous findings to show how MuSyC generalizes the DEP and MSP, thereby unifying the field of drug synergy, as sought at Saariselkä. Further, we map the landscape of current synergy metrics, including: Bliss Independence^{3}, Loewe Additivity^{1}, Combination Index (CI)^{22}, Highest Single Agent (HSA)^{23}, Effective Dose model^{6}, ZIP^{7}, a partial differential equation (PDE) Hill model by Schindler^{8}, BRAID^{11}, and the General Pharmacodynamic Interaction (GPDI) model^{24}. In mapping relationships between these various frameworks, we identified systematic differences impacting the interpretation of synergy in drug combination experiments. Specifically, we found: (1) the conflation of synergistic potency and efficacy masks synergistic interactions; (2) MSP-based frameworks are biased toward antagonism for drugs with intermediate efficacy; and (3) DEP-based frameworks contain a Hill-slope dependent bias. The Hill-slope bias results from satisfying the famous “sham” combination thought experiment, thus arguing against the merit of sham-compliance as a measure of validity for synergy frameworks. Using five large combination datasets^{25,26,27,28,29}, MuSyC identifies real-world examples where the conflicting assumptions of previous drug synergy frameworks misleads or impedes drug discovery efforts through these pervasive and predictable biases. Additionally, we show that MuSyC uncovers two consequential errors in the highly cited CI^{22,30} which has been proposed as the standard Mass Action-based, synergy framework^{31}. We therefore propose MuSyC as a consensus framework to interpret combination pharmacology and signify its broad applicability to the study of drug mixtures.

## Results

### A state-transition model to measure multi-drug synergistic effects

The 4-parameter Hill equation is commonly used to fit dose-response data from in vitro and in vivo assays (see Box 1 Eq. (10) and Table 1 for parameter annotation). Here we derive this equation from the equilibrium of a two-state transition model of drug effect based on the Law of Mass Action (Fig. 1A, left). Traditionally, the parameters of the Hill equation are interpreted as a drug’s efficacy (*E*_{0} − *E*_{1}), potency (*C*), and cooperativity (*h*), also known as the Hill slope. These parameters correspond to three possible geometric transformations of a dose-response curve (Fig. 1A, right). To generalize this one-drug formalism to two concurrent drugs, we developed a four-state transition model of combination pharmacology (Fig. 1B, left)^{21}. From this model, we derive a two-dimensional (2D) Hill equation for two drugs (Box 1, Eq. (15)) defining a dose-response surface (Fig. 1B, middle). The 2D Hill equation contains five additional parameters, not present in the single-drug Hill equation, which measure different types of drug interactions. These additional parameters measure changes in a drug’s efficacy (*β*), potency (*α*_{12} and *α*_{21}), and cooperativity (*γ*_{12} and *γ*_{21}) in a combination—corresponding to three distinct types of synergy (Fig. 1B, right, Table 1). See Supplemental Code 1 and Supplemental Section Interactive MuSyC Jupyter Notebook for an interactive demonstration of the 2D Hill equation parameters. As we show below, these parameters are conflated in traditional synergy metrics (e.g. Loewe, Bliss, and HSA), as well as in recently proposed ones obscuring the true origin and magnitude of drug synergy or antagonism.

### Mapping the landscape of prominent synergy models within a consensus framework

Multiple alternative synergy models have been proposed, most broadly derived from one of two guiding principles: the Multiplicative Survival Principle (MSP) or the Drug Equivalence Principle (DEP) (Table 2). Prior work has shown contradictory results when comparing between MSP and DEP frameworks^{4,12,32}, and a lack of consensus remains on the commonality between the two principles^{4,7,9,11}. Here we show MuSyC satisfies both the DEP and MSP under specific parametric constraints (Fig. 2A, B), thereby unifying the foundational principles of drug synergy.

The MSP was first described by Bliss^{3} and is the foundation of the Bliss Independence framework. MSP assumes the probability of a cell surviving treatment by drug 1 (*U*_{1}) is independent of the probability of the same cell surviving treatment by drug 2 (*U*_{2}). Therefore, the probability of surviving both Drug 1 and Drug 2 is *U* = *U*_{1} ⋅ *U*_{2}^{3}. Synergy or antagonism occur when *U* ≠ *U*_{1} ⋅ *U*_{2}. A method to define an alternative Bliss null model has been reported for growth-rate data as the sum of growth-rate inhibition^{33}, but this formulation is uncommon, and not classified as MSP. MuSyC satisfies the MSP under the following conditions: (1) the effect metric is expressed as a percent (*E*_{0} = 1, and *E*_{3} = *E*_{1}*E*_{2}), (2) there is no potency synergy (*α*_{12} = *α*_{21} = 1), and (3) there is and no cooperativity synergy (*γ*_{12} = *γ*_{21} = 1) (Fig. 2A, see Supplemental Section Multiplicative Survival Principle for details).

The DEP was first established by Loewe^{1,2}. DEP-based methods are characterized by linear isoboles (contours of equal effect) (Fig. S2A). A combination of doses *d*_{1} and *d*_{2}, achieving effect *E*, satisfies the DEP when \(\frac{{f}_{1}({d}_{1})}{{f}_{1}^{-1}(E)}+\frac{{f}_{2}({d}_{2})}{{f}_{2}^{-1}(E)}=1\), where *f*_{i}(*d*_{i}) represents the monotherapy response of drug *i*. MuSyC satisfies the DEP under the following conditions: (1) the drugs’ actions are mutually exclusive (*α*_{12} = *α*_{21} = 0) and (2) *h*_{1} = *h*_{2} = 1 (Figs. 2B and S2B, see Supplemental Section Dose Equivalence Principle for details).

From the literature, we identified several prominent synergy models beyond Bliss and Loewe including: CI^{34}, HSA^{23}, Effective Dose model^{6}, ZIP^{7}, Hill PDE^{8}, and GPDI^{24}. Table 2 compares key features and assumption of the different synergy models. Each of these methods, as well as MuSyC, defines synergy based on the experimental deviation from a null (additive) dose-response surface. Because almost all synergy frameworks are founded on either the DEP or MSP, we standardized relationships between the various models, mapping the global landscape of drug synergy (Fig. 2C, see Supplemental Section Relationships between different synergy frameworks for details).

In deriving this map, we uncovered potential sources of error when using MSP or DEP methods which impact interpretation of synergy studies. Specifically, we identified three recurrent considerations meriting attention from the field. (1) Previous synergy metrics conflate different synergy types (i.e. potency, efficacy, cooperativity) in ways that can mask synergistic and antagonistic interactions (Fig. 3). (2) The connection between MuSyC and the MSP-derived frameworks depend on the single drugs’ efficacy (*E*_{1}, *E*_{2}), and as a result, MSP frameworks are biased against the combination of moderately efficacious single agents (Fig. 4). (3) The connection between the DEP and MuSyC is constrained by single drugs’ Hill slopes (*h*), and therefore DEP frameworks impose a Hill-slope dependent bias, artificially inflating the synergy for drugs with low Hill slopes (Fig. 5). This bias is a consequence of satisfying the sham experiment as described in Section Re-examining the sham experiment: Sham compliance introduces Hill-dependent bias in DEP models. While there may be different approaches to quantify synergy relative to a given null model, the biases we discuss are intrinsic to the null models themselves, while the magnitude of bias depends on the precise synergy quantification. To assess impact of these considerations on synergy calculations in different fields, we analyzed five large publicly available datasets (Table 3) using MuSyC and other synergy frameworks (see Methods section for description of fitting methods). We find the MuSyC algorithm to be robust to different types and magnitude of noise as well as sampling designs facilitating the analysis of many types of data (Fig. S3). Note, while synergistic cooperativity (*γ*) is theoretically plausible (as initially postulated by^{7}), including it did not increase the fit quality (Fig. S4) as measured by AIC and therefore we do not explore synergistic cooperativity in subsequent analyses.

### Conflating synergistic potency and efficacy masks synergistic interactions

To determine how conflation of distinct synergy types impacts the interpretation of drug-response data, we generated synthetic dose-response surfaces using MuSyC (Eq. (15)) across a range of *α* and *β* values and calculated the synergy according to Loewe, Bliss, and Highest Single Agent (HSA) at the EC50 of both drugs (Fig. 3A, D, G and Video S2). In each case, many distinct sets of (*α*_{12}, *α*_{21}, *β*) are indistinguishable (e.g. the black contour line on the spheres).

Figure 3A shows that near the boundary between synergism and antagonism, Loewe is insensitive to changes in synergistic potency, tracking instead with synergistic efficacy. Consequently, in the anti-cancer dataset from O’Neil et al.^{26}, Loewe misses potency antagonism in combinations with synergistic efficacy (Fig. 3B middle distribution, see Fig. S5 for an example surface). This reflects Loewe’s assertion of infinite potency antagonism (*α*_{12} = *α*_{21} = 0, Fig. 2A) in its null model. Therefore, combinations that are antagonistically potent (*α* < 1) are all synergistic by Loewe in the absence of sufficient antagonistic efficacy (values above black contour in Fig. 3A). Indeed, Loewe is frequently synergistic even in cases of antagonistic potency and efficacy (Fig. 3B bottom distribution). As an example, the combination of methotrexate (targets folate synthesis) and erlotinib (EGFR inhibitor) in UWB1289 (BRCA1-mutant ovarian carcinoma) cells is antagonistically efficacious and potent by MuSyC, but synergistic by Loewe (Fig. 3C).

Some assays, such as the anti-malarial combination screen from Mott et al.^{25} are designed so that the drug effect spans the entire range from *E*_{0} = 1 to *E*1 = *E*2 = *E*3 = 0. In such cases there is no synergistic efficacy, because each drug alone already achieves the greatest efficacy measurable. Nevertheless, even if all synergy detected by an assay is synergistic potency, traditional synergy metrics can still conflate asymmetric synergistic potencies. This can occur when *α*_{12} is synergistic while *α*_{21} is antagonistic, or vice versa (Fig. 3D, Bliss synergy, black contour line through *β* = 0 plane). In the anti-malarial dataset^{25}, Bliss is consistently synergistic when *l**o**g*(*α*_{12}, *α*_{21}) > 0, and antagonistic if *l**o**g*(*α*_{12}, *α*_{21}) < 0; however, when *l**o**g*(*α*_{12}) < 0 < *l**o**g*(*α*_{21}), Bliss will strictly classify a combination as either synergistic or antagonistic (Fig. 3E bottom distribution) despite the asymmetric interactions. As an example, Bliss conceals that halofantrine (inhibits polymerization of heme molecules) reduces the potency of mefloquine (targets phospholipids) against the multi-drug resistant malaria strain HB3 (Fig. 3F).

HSA is commonly thought to quantify synergistic efficacy. However, for antagonistically potent combinations, HSA cannot distinguish synergistic and antagonistic efficacy because it does not account for the topology of the dose-response surface (compare \(({{{\mathrm{log}}}}\,({\alpha }_{12}),{{{\mathrm{log}}}}\,({\alpha }_{21}),\beta )=(-,-,+)\) and ( − , − , − ) quadrants of Fig. 3G and Video S2). In the anti-cancer combination dataset^{26}, we observe this trend (Fig. 3H middle vs bottom distributions). As an example, the synergistically efficacious combination of dexamethasone (agonist of the glucocorticoid receptor) and mk-8669 (PI3K/mTOR dual inhibitor) in a colorectal adenocarcinoma cell-line is masked by HSA due to antagonistic potency (Fig. 3I). Repressing glucocorticoid signaling has previously been shown to repress mTOR signaling^{35} providing a potential molecular mechanism explaining the synergy.

### MSP is biased against combinations of drugs with intermediate efficacy

MSP frameworks, such as Bliss, explicitly expect drug effects to measure the “percentage of cells affected”, which is by definition bounded within the closed interval *E* ∈ [0, 1]. Nevertheless, dose-response data is usually not a measure of percent affect, but rather of relative percent effect (see Supplemental Section Percent Affect vs Percent Effect for an example). This distinction, maintained by MuSyC (Box 1), is critical because percent effect data commonly saturates (i.e., percent affect is near 100%) at intermediate effect (i.e., relative percent effect is near 50%). For combinations of these moderately efficacious drugs, Bliss expects a large increase in effect over the single agents, even when each drug is administered at saturating concentrations (Fig. 4A, middle panel). In contrast, if combining drugs with high or low efficacy, Bliss expects a more modest increase (Fig. 4A, left and right panels).

Based on this expectation that *E*_{3} = *E*_{1} ⋅ *E*_{2} (Fig. 2A), MuSyC predicts Bliss would be biased toward antagonism in combinations of moderately efficacious drugs (Fig. 4B yellow shading around *E*_{1} = *E*_{2} ≈ 0.5). As expected, the median Bliss score in the anti-cancer dataset is biased toward antagonism for moderately efficacious combinations 0.35 < (*E*1, *E*2) < 0.65 (Fig. 4C, cyan square), though the magnitude of bias is less than predicted in Fig. 4B. This bias persists even when looking at pan-cancer trends in the combination of drugs which have, on average, intermediate effect over the entire cell-line panel (Fig. 4D). As a particular example, the synergistic efficacy of paclitaxel (targets microtubule stability) and mk-2206 (AKT inhibitor) in KPL1 cells is masked by Bliss’s high expectation for moderately efficacious drugs (Fig. 4E, gray plane). In several datasets, the magnitude of this MuSyC-predicted bias was sufficient to obfuscate many of the strongest synergies and antagonisms according to Bliss (Fig. S6A, B).

Other MSP-based methods, such the Effective Dose model also assume data measures percent affect and fit a simplified 2-parameter Hill equation enforcing *E*_{0} = 1 and *E*_{1} = 0. This assumption can lead to poor fits of percent effect data for moderately efficacious drugs, and thus invalid synergy scores (Fig. S7). Therefore, the distinction between percent affect and percent effect is a critical component of MuSyC.

### Re-examining the sham experiment: Sham compliance introduces Hill-dependent bias in DEP models

A new synergy model’s consistency is traditionally tested with the “sham” combination thought experiment. In a sham experiment, a single drug is considered as though it were a combination, with the expectation that the drug should be neither synergistic nor antagonistic with itself. DEP frameworks, characterized by linear isoboles, are known to satisfy the sham experiment, while MSP frameworks famously do not^{9}.

In Box 2, we show MuSyC only satisfies the sham experiment when *h* = 1, which makes sense as MuSyC produces linear isoboles only in this condition (Fig. S2B and Eq. 25). Further, our analysis in Box 2 revealed that sham combinations exhibit unique biochemistry, only equivalent to true combinations in the case *h* = 1. When *h* ≠ 1, combinations contain intermediate states representing mixed-inhibition (black circles, Fig. 5A, B). In sham combinations, these mixed-inhibition states are equivalent to single drug complete-inhibition states (Fig. 5A, cyan circles), while for real combinations, these are *not* equivalent (Fig. 5B). When *h* = 1, these intermediate mixed-inhibition states do not exist, explaining the concordance between sham and true combinations in this case. We expect when *h* ≠ 1, enforcing sham compliance leads to predictable systematic biases. We expect when *h* < 1 DEP frameworks will overestimate synergy (Fig. S2C), and when *h* > 1, DEP frameworks will overestimate antagonism (Fig. S2D).

In combinations from the anti-cancer dataset, the average trend of Loewe synergy closely follows the Hill slope bias predicted by MuSyC (Fig. 5C). Further, subtracting the MuSyC-predicted bias from Loewe values for each combination results in a distribution independent of Hill slope (bottom panel). The bias toward synergy is particularly large for drugs with low Hill slopes. As an example, both doxorubicin (DNA damaging agent) and mk-4827 (PARP inhibitor) have small Hill slopes when applied to MBA-MB-436 cells, and their combination is synergistic by Loewe. However, using MuSyC, we see this combination is both antagonistically efficacious and antagonistically potent (Fig. 5D). In one dataset (Cokol et al.), this MuSyC-predicted bias revealed a screen-wide underestimation of synergy by Loewe (Fig. S6A,C).

Therefore, satisfying sham compliance biases models toward synergy for drugs with low Hill slopes, regardless of with what these drugs are combined. This bias—which stems from enforcing a biochemical reaction scheme only appropriate for sham combinations—should be sufficient grounds for dismissing the sham experiment as a measure of a new synergy framework’s validity.

### MuSyC reveals errors in the derivation and application of the Combination Index

Recent reports have identified potential flaws with the use of CI^{11,36}, yet it remains the most highly cited synergy metric^{30,34,37}. CI has recently been proposed to the Food and Drug Administration (FDA) and National Institutes of Health (NIH) as the de facto definition of drug synergy^{31}. Due of its prominence in the field, here we specifically examine its behavior with respect to the biases discussed above. We find CI and MuSyC have the same null model when *h*_{1} = *h*_{2} = 1, *α*_{12} = *α*_{21} = 0, *E*_{0} = 1, *E*_{1} = *E*_{2} = 0 (Supplemental Section Relationship between different synergy frameworks). The presence of constraints on *h* and *E* indicates CI could be impacted by both a Hill-slope and efficacy range bias, like those we reported above for DEP and MSP-based frameworks, respectively. Like MuSyC, CI is based on the Law of Mass-Action, facilitating a direct comparison of their formulations. In doing so, we found two errors in CI with significant consequences: (1) a fundamental math error in its derivation impacting combinations with *h* ≠ 1 (Box 3, Fig. 6A–C), and (2) a fitting error that arises when applied to drugs with partial efficacy (*E*_{1} or *E*_{2} > 0) (Fig. 6D–F).

Details of the error in the derivation of CI are in Box 3. Because of this error, the CI equation is only valid for combinations of drugs with Hill slopes equal to one (*h*_{1} = *h*_{2} = 1). In this regime, MuSyC also results in linear isoboles (Fig. S2B). When *h* ≠ 1 (Figs. 6A and S8A), the CI equation incorrectly factors the exponent outside the sum, introducing the exact same cross terms that we show (Box 2) are only valid for sham combinations, and not valid for all combinations when *h*_{1} ≠ *h*_{2} ≠ 1. Therefore, when *h* ≠ 1, CI suffers from the same Hill-slope dependent bias we show in Fig. 5. We show the consequence of this bias in illustrative two-drug examples using synthetic (Figs. 6B and S8B) and experimental combinations (Figs. 6C and S8C).

Even when *h* = 1, CI requires that the maximal and minimal effects (*E*_{0}, *E*_{1}, Fig. 1) be fixed at 1.0 and 0. Subsequently, CI fits the single-drug dose-response curves using the two-parameter median-effect equation^{30} (Eq. (8)), see also Supplemental Section Percent Affect vs Percent Effect—Combination Index). In contrast, MuSyC is based on a four-parameter Hill equation (Box 1, Eq. (10)), and thus can describe metrics of drug effects with arbitrary ranges. It is common to observe in many dose-response assays, such as percent viability, drug effects that do not reach 0% (i.e. *E*_{1} > 0)^{38}. For such assays, the two-parameter median-effect equation fits poorly (Fig. 6D). These poor fits lead to an effect-dependent error in CI quantification of synergy as observed in synthetic (Fig. 6E) and experimental combinations (Fig. 6F). We note, this effect-dependent error is not the same as the effect-range bias we report for MSP-based frameworks in Section “MSP is biased against combinations of drugs with intermediate efficacy”.

## Discussion

Herein, we have demonstrated four key advances of MuSyC^{21} germane to the study of combination pharmacology: (1) the unification of the DEP and MSP; (2) the decoupling of three distinct types of synergy; (3) the revelation of biases emerging from constraints on the single drug pharmacological profile inherent in the DEP and MSP; and (4) reporting on flaws of highly cited CI including a mathematical, derivational error which impacts its reliability for synergy quantification.

The DEP and MSP have formed the foundational principles of most synergy frameworks over the last century; however, the connection between these principles has remained unknown^{4,5}. Here, approaching combination pharmacology using the Law of Mass Action applied to a state-transition model results in a single framework unifying both principles. By mapping all frameworks on a common landscape, MuSyC facilitates rigorous investigation of oft-cited, contradictory conclusions between existing frameworks^{16}—contradictions that preclude reproducibility between synergy studies. Specifically, as is seen in Fig. 2C, there is no combination which can simultaneously satisfy the conditions required by both DEP and MSP synergy frameworks. Previous works advocate prioritizing combinations that are synergistic by all methods^{5}, or choosing a synergy model carefully based on the data measured and shape of dose-response curves^{32}. In contrast, MuSyC’s unification of DEP and MSP means it is applicable regardless of whether the data are best described by the DEP or MSP, or something in-between, so no choice is required. Further, while MuSyC’s mass-action formalism applies most directly to molecular inhibition, its synergy parameters describe geometric transformations of efficacy, potency, and Hill slope, which are well-established quantities used to describe sigmoidal responses, regardless of mechanism.

One key advance of MuSyC, facilitating this unification, was the decoupling of *α*, *β*, and *γ*. These synergy parameters correspond directly to classic, pharmacological measures of a drug’s potency, efficacy, and cooperativity. By calculating synergy in this way, interpretation of synergy does not depend on arbitrary expectations or thresholds. Rather, an *α* of 10 corresponds to a 10-fold increase in a compound’s potency, as a result of the other drug, regardless of whether we define *α* = 1 or *α* = 10 as the “threshold” for synergy. As practical advice for accurately fitting all synergy parameters, we recommend sampling (*d*_{1}, *d*_{2}) around the four corners (0, 0), (*d*_{1,max}, 0), (0, *d*_{2,max}), and (*d*_{1,max}, *d*_{2,max}) to best constrain synergistic efficacy, and around the four edges (*C*_{1}, 0), (*C*_{1}, *d*_{2,max}), (0, *C*2), (*d*_{1,max}, *C*2) to best constrain synergistic potencies and cooperativities, where *d*_{i,max} is an asymptotically high dose of drug *i*. We refer interested readers to Supplemental Section “Interactive MuSyC Jupter Notebook”, Figs. S12–S16, and Supplemental Code 1, which provide an interactive demonstration that shows how each synergy parameter results in different outputs across multiple concentrations. We envision distinguishing synergies of potency, efficacy, and cooperativity will be of differential consequence in alternate contexts. For example, in cancer synergistic efficacy may be most important, while for neurological disorders, synergistic cooperativity—i.e. sharp on-off drug response —may be preferred. In an analysis of clinical trials of combination therapies, we find synergistic efficacy is statistically higher in clinically efficacious combinations than clinically non-efficacious combinations (see Supplemental Section MuSyC statistically distinguishes efficacious and non-efficacious drug combinations in clinical trials based on combination screens, Figs. S9 and S10).

The relationship between MuSyC and the MSP and DEP frameworks (Fig. 2) is constrained by monotherapy parameters (*E*_{1}, *E*_{2} for MSP, *h* for DEP). These constraints suggested systematic biases in MSP and DEP frameworks contingent on a single drug’s efficacy (MSP, Fig. 4) and Hill slope (DEP, Fig. 5). These systematic biases merit consideration when using these frameworks for drug discovery in large screens, or when accounting for batch-effects across different datasets (Fig. S6). Such systematic biases can confound machine learning models to predict synergy, decreasing their utility. Additionally, the constraint on *h* highlighted a discrepancy between the biochemistry of true sham experiments and real combinations. The centrality of the sham experiment in the drug synergy literature cannot be overstated; however, we argue enforcing sham compliance comes at the cost of improperly modeling real combinations, leading to a predictable Hill-dependent bias.

CI has previously been criticized for its procedure to fit the median-effect equation. Specifically, CI depends on a log-transformation in order to linearize the two-parameter median effect equation “in similar logic to the defunct Scatchard analysis in pharmacology, which has been replaced by non-linear modeling”^{36}. That is, this log-transformation alters the noise profile such that small deviations for effects at low concentrations result in large deviations in the synergy calculations. Additionally, this transformation results in information loss, since undefined values for effects outside the range of (0,1) are forcibly removed, and as a consequence, CI synergy estimates become “statistically unstable” for noisy experimental data^{11}. Beyond these valid points regarding CI’s practical application, here we uncover a mathematical error in the derivation of CI (Box 3) causing a systematic bias depending on the Hill slope. Because these flaws compound in non-linear ways, the expected error when applying CI is unique to each particular combination and assay design.

The prospects of higher-order synergies (i.e., interactions beyond pairwise) and scaling laws for drug mixtures, while provocative, have remained contentious^{6,12,39,40}. MuSyC’s cubic geometry allows it to be easily extended to three or more drugs (Fig. S1), and we expect MuSyC will enable a more refined search for higher-order interactions. For instance, combinations that mix different synergy profiles (e.g., drugs 1 and 2 are synergistically potent, and drugs 2 and 3 are synergistically efficacious) may exhibit different higher-order interactions than combinations all sharing a single synergy type. However, the number of synergy parameters in MuSyC scales as 2^{N}(*N* + 1) − 3*N* − 1 (including *γ*) where *N* is the number of drugs (Fig. S1), and the commensurate data necessary to fully constrain MuSyC hyper-surfaces invokes a parameter identifiability problem (“the curse of dimensionality”). Nonetheless, MuSyC’s geometry could be leveraged to guide sampling schemes to constrain the boundaries, allowing the solution to be built up step-wise from the boundaries (see Supplemental Section Proof of boundary behavior of the 2D Hill equation).

MuSyC expects single-drug dose-response curves to be sigmoids well fit by the 1D Hill equation (Eq. (10)), and dose-response surfaces to be well fit by the 2D Hill equation (Eq. (15)). In our experience, these expectations are met by real data, as most single drugs have monotonic, sigmoidal responses, and even complex drug interactions can be modeled using various mixtures of *α*, *β*, and *γ* (96% and 88% of combinations in anti-cancer and anti-malarial datasets had *R*^{2} > 0.7, respectively). However, it is possible for drugs to have multiphasic responses due to poly-pharmacology which are not well fit by a Hill curve. It may be possible to extend MuSyC to encompass such drugs—for instance by including a multiphasic Hill model^{41} or modeling effects of “partially affected” states (Fig. 5A, B and Box 3). In extreme cases, it may only be possible to apply non-parametric frameworks such as Bliss, Loewe, or HSA. Nevertheless we note that without fitting dose-response curves to a parametric model, these metrics are sensitive to noise in individual data-points. Replicate measurements may be able to reduce this sensitivity. Additionally, MuSyC assumes all drugs are administered concurrently, whereas patient treatments are often staggered. New theory and experimental methods are needed to address the synergy of combinations which are staggered temporally, bridging the synergy of pharmacodynamics with the synergy of pharmacokinetics. Finally, in the datasets we analyzed here, we did not find a role for synergistic cooperativity (*γ*). Future studies in other systems are needed to better understand situations when synergistic cooperativity is expected.

By viewing the landscape of drug synergy through the lens of mass-action, we have demonstrated the underlying assumptions, limitations, and biases of commonly applied synergy methods. We have shown how MuSyC unifies the DEP and MSP thus providing a consensus framework for the study of combination pharmacology. These findings provide much needed clarity to the field and should promote the reproducibility of drug synergy studies across drug combination discovery efforts. Such a rigorous approach to the discovery and application of drug combinations will open the door to the discovery of new and previously discarded avenues for therapeutic mixtures.

## Methods

We note the synergy calculations conducted for the different published datasets were not necessarily the same as those used in the original paper. Indeed the limitations of the current frameworks forced customized analysis for each publication highlighting the need for a consensus framework. However, in order to compare between datasets, we have calculated Bliss, Loewe, HSA, and other synergy frameworks, as described below, from the raw data.

### Software

#### Implementation and website

A web application to calculate MuSyC synergy parameters for users’ data is available at https://musyc.lolab.xyz/. Experimental data are uploaded in comma-separated value (CSV) format; data format details and usage instructions are in the supplemental materials. The application fits dose-response surfaces using MuSyC and offers the results both as a CSV download of fit parameters, and interactive plots of the dose-response surface.

The web application uses the Django web application framework (djangoproject.com) and Python 3.7. Fitting tasks are processed asynchronously using a message queue (RabbitMQ; rabbitmq.com) and task-worker framework (Celery; celeryproject.org). Data are organized in a Postgres relational database (postgres.org). The following packages were used for fitting, data analysis, or visualization: SciPy v1.1.0^{42}, Numpy v1.14.3^{43}, Pandas v0.23.0^{44}, Matplotlib v2.2.2^{45}, uncertainties v3.0.2^{46}.

#### Fitting 2D Hill equation

Here we describe fitting protocol for drug metrics where the metric of drug effect decreases as dose increases (*E*_{0} > *E*_{3}) (e.g., anti-proliferative drugs); however, the framework is equally valid if increasing the drug corresponds to increases the effect (*E*_{0} < *E*_{3}) (e.g., percent effect).

Previously, we found it necessary to use a Metropolis Hastings Monte Carlo (MCMC) seeded with a particle swarm optimization (PSO) to fit the 2D Hill equation^{21}. This was prompted by the inconsistent performance of standard non-linear least squares (NLLS) regression. In particular, we observed instances of calculated uncertainties in NLLS which were several orders of magnitude greater than the parameter value. This, we have discovered, is because the multi-collinearity between the Hill slope and the EC50 (*C*) inherent in the structure of the Hill equation—collinearities which are amplified when extending the Hill equation to 2D. The correlation between variable *h* and *C* is easiest to observe in a linearized 1D Hill equation (Eq. (1)).

In Eq. (1), the intercept of the line (\(h* {{{\mathrm{log}}}}\,(C)\)) depends on the slope of the line (*h*). This correlation is problematic when trying to estimate the parameter uncertainty (*σ*) from NLLS regression because *σ* is estimated as the square root of the inverse Hessian, approximated to be **J**^{T}**J** (where **J** is the Jacobian at the solution). **J** of the 4-parameter Hill equation is

When the Hill slope is large (e.g., *h* > 5), the second two terms of the **J** cause the numerical approximation of the inverse of *J* to be undefined. This problem is compounded in the 2D Hill equation where, in addition to *h* and *C*, the parameters *α* and *γ* are co-linear. However, this does not affect the accuracy of the fitted parameter values from the NLLS regression—only the parameter uncertainty^{47}.

For the fitting the 2D Hill equation in this study, we adopted a Monte Carlo sampling approach to calculate the fit parameter uncertainty. This is significantly faster than our previous method (PSO + MCMC)^{21} while maintaining reasonable calculations of the parameter uncertainties accounting for the multi-collinearities described above. The Monte Carlo algorithm for fitting the 2D Hill equation is as follows. First, the 4-parameter 1D Hill equation (Eq. (10)) is fit to the dose-response of each drug in isolation. The fit uses the Trust Region Reflective (TRF) algorithm implemented in the *c**u**r**v**e*_*f**i**t*() module of the scipy.optimization package. *h* and *C* were unconstrained while *E*_{max} and *E*_{0} are constrained for each dataset as annotated in the Methods section, data acquisition, preparation, and analysis. In general, adjusting the parameter bounds to closely match what is feasible for the given dataset will lead to better parameter estimates, helping the curve fitting algorithm to avoid becoming stuck in a suboptimal local minima. The initial 1D Hill fits provide estimates for (*E*_{0}, *E*_{1}, *E*_{2}, *C*_{1}, *C*_{2}, *h*1, *h*2), because the 2D Hill equation becomes equivalent to the 1D Hill equation in the limit as *d*_{i} → 0. In practice, best fits of these parameters in the 2D Hill equation which have counterparts from the 1D Hill equation tend to be similar to their monotherapy fits which were used as initial guesses. However we note that it is possible for these values to differ significantly from their monotherapy best guesses when the monotherapy data are noisy, and thus can have wide uncertainties. Next the 2D Hill equation (Eq. (15)) is fit using the TRF algorithm with initial values based on the 1D Hill equation fits and with bounds based on the parameter uncertainty calculated for the 1D Hill fits. The initial values for parameters unique to the 2D Hill equation, *E*_{3}, *α*_{21}, *α*_{12}, *γ*_{12}, *γ*_{21} are (\(\min ({E}_{1},{E}_{2})\),1,1,1,1). For all combinations *r*_{1} = *r*_{2} = 100. The bounds for \({{{\mathrm{log}}}}\,({\alpha }_{21}),{{{\mathrm{log}}}}\,({\alpha }_{12})\) are set to [−4,4]. From this initial fit, 100 Monte Carlo samples are used to calculate the parameter uncertainty as described by Motulsky and Christopoulos^{47}, (Chapter 17: Generating confidence intervals by Monte Carlo, pg. 104). Specifically, noise, with a distribution N(0,*σ*), where *σ* is equal to the root mean square (RMS) of the best fit, is added to best-fit values of the 2D Hill equation for all drug doses. The data plus noise is then fit again initializing the optimization from the best-fit parameters of the original data. This is done 100 times. From this ensemble of fits, the 95% confidence interval of each parameter can be calculated. This Monte Carlo approach results in asymmetric confidence intervals which better captures the non-Gaussian distribution of uncertainty for many fits (e.g. the distribution of *h* is log-normal) as well as being robust to the co-linear parameters in the 2D Hill equation. The asymmetric confidence interval is particularly salient when the dose-range is insufficient to observe the lower plateau of the dose-response. Only combinations for which *R*^{2} > 0.7 and the fitted EC50s of both drugs was less than maximum tested dose for each (\({C}_{1} < \max (d1),{C}_{2} < \max (d2)\)) were included for subsequent analysis.

#### Comparing fitting algorithm robustness between different synergy frameworks

We additionally examined the performance of the MuSyC fitting when the raw data is subject to different types of noise (Fig. S3A). We synthetically generated 10 random samples of 5400 noise and synergy profiles and compared between all parameterized models of synergy the percent of convergence (Fig. S3B), fit quality assessed by *R*^{2} (Fig. S3C), and variation in synergy parameters assessed by Z-score between the random samples (Fig. S3D). Fitting algorithms for the different models is described in Section Calculation of other synergy metrics. Overall, we find MuSyC performs as well as or better than comparable parameterized models of synergy on the tested synergy profiles. We note this analysis is hampered by the lack of a “true” standard for synergy necessitating dose-response surfaces to be generated based on a defined model—in this case the MuSyC model. Bliss, Loewe, and HSA in general do not require fitting an equation to the data, and were thus excluded from this analysis. Nevertheless we note such synergy metrics may be more sensitive to noise in the data, as noise in individual datapoints is not smoothed out via a curve fit.

### Data acquisition, preparation, and analysis

#### ONeil et al. anti-cancer screen

The anti-cancer drug combination data were downloaded from the supplemental materials of ONeil et al.^{26}. Single agent and combination datasets were merged. Drug effect was the mean normalized percent viability (*X*/*X*_{0} column) calculated as detailed in ONeil et al.^{26}. The minimum and maximum bounds for [*E*_{0}, *E*_{1}, *E*_{2}, *E*_{3}] during 2D Hill equation fits were [0.99,0.0,0.0,0.0,0.0] and [1.01,2.5,2.5,2.5], respectively.

#### Mott et al. anti-malarial screen

The anti-malarial drug combination data were downloaded from https://tripod.nih.gov/matrix-client/ from the Malaria Matrix project. Blocks downloaded for analysis were: [1601,1602,1603,1701,1702,1703,1761,1764]. Only blocks with a mqcConfidence of >0.9 were included. The drug effect was calculated as described in Mott et al.^{25}. Effects <−20% and >120% were removed. The minimum and maximum bounds for [*E*_{0}, *E*_{1}, *E*_{2}, *E*_{3}] during 2D Hill equation fits were [90.,0.0,0.0,0.0] and [110,200,200,200], respectively.

#### Tan et al. anti-HIV screen

The anti-HIV drug combination data were downloaded from the supplemental table four of Tan et al.^{28}. Drug effect was one minus the normalized infection rates as detailed in Tan et al.^{28}. The minimum and maximum bounds for [*E*_{0}, *E*_{1}, *E*_{2}, *E*_{3}] during 2D Hill equation fits were [0.99,0.0,0.0,0.0,0.0] and [1.01,1.5,1.5,1.5], respectively.

#### Cokol et al. anti-fungal screen

The anti-fungal drug combination data were downloaded from supplemental dataset one in Cokol et al.^{29}. The raw cell growth measurements for all 200 drug-drug interaction assays were stored as a 96 × 64 matrix of numbers. Rows were time points with 15 min intervals and columns are the indices of an 8 × 8 drug matrix. Drug dilutions were linear between the maximum reported in Table 1 of Cokol et al.^{29} and 0. The drug effect was quantified using the area under the growth curve (AUGC), calculated using Simpson’s integration, after the first 10 time points (150 min). The background unique to each experiment was removed by subtracting the minimum observed growth rate for each pair individually. The minimum and maximum bounds for [*E*_{0}, *E*_{1}, *E*_{2}, *E*_{3}] during 2D Hill equation fits were [0, 0, 0, 0] and [*∞*, *∞*, *∞*, *∞*], respectively.

#### Holbeck et al. anti-cancer screen

The ALMANAC anti-cancer drug combination data were downloaded from https://wiki.nci.nih.gov/display/NCIDTPdata/NCI-ALMANAC file *C**o**m**b**o**D**r**u**g**G**r**o**w**t**h*_*N**o**v*2017. *z**i**p*. The matching single dose-response data were downloaded from https://wiki.nci.nih.gov/display/NCIDTPdata/NCI-60+Growth+Inhibition+Data, June 2016 release *D**O**S**E*_*R**E**S**P**O**N**S**E* link. Single agent and combination datasets were merged using pandas dataframe operations. Drug effect was the percent growth inhibition calculated as detailed in^{27}. The minimum and maximum bounds for [*E*_{0}, *E*_{1}, *E*_{2}, *E*_{3}] during 2D Hill equation fits were [99, −100, −100, −100, −100] and [101, 350, 350, 350], respectively.

### Calculation of other synergy metrics

#### Bliss, Loewe, and HSA

Bliss, Loewe, and HSA depend on the concentration of drugs so a combination can be synergistic at one dose, but antagonistic at another dose. Several methods have been proposed for extracting a single synergy metric per combination from these frameworks to enable comparisons between combinations^{13,14,15,26}. For our analysis, we calculate the synergy score at the combination of each drug’s EC50 (Figs. 3 and 5) as proposed by Malyutina et al.^{48} or at the maximum tested concentration of each drug (Fig. 4). The EC50 of each drug was calculated from the fits to the 2D Hill Eq. (15) which we have observed to be more robust to noise when estimating the single drug pharmacologic profile. Assuming the notation for the 1D Hill equation and inverse Hill equation—which calculate effect (*E* given a dose (*d*) and a dose given an effect, respectively—are given by

where *E*_{x} < *E*_{0}, then equations for Bliss, Loewe, and HSA at the EC50 are:

where *E**d*(*C*_{1}, *C*_{2}) is the measured effect of combining *C*_{1} of *d*_{1} and *C*_{2} of *d*_{2}. And equations for Bliss at the max of each drug is:

Thus, Loewe synergy is calculated using an equation similar to CI^{34}, while Bliss and HSA are calculated using an “excess over” approach, which calculates the raw difference between the expected and observed responses. While the reference models are always the same, we note alternative equations may be used to quantify synergy of a combination^{33,49}, though the biases we report are intrinsic to the reference models, not the synergy quantification approach. These equations assume the metric of drug effect decreases as the dose increases. Because many single agents did not reach 0% maximum efficacy, the EC50s (*C*_{1}, *C*_{2}) were not necessarily 50% (Fig. S7). If *E*(*C*_{1}, *C*_{2}) < *E*_{1}, *E*_{2} then Loewe was undefined. We apply a \(-{{{{\mathrm{log}}}}\,}_{10}\) transformation the scale Loewe to match the ranges Bliss and HSA are synergistic; therefore, f \(-{{{{\mathrm{log}}}}\,}_{10}(Loewe)\, > \, 0\) the combination is synergistic and if \(-{{{{\mathrm{log}}}}\,}_{10}(Loewe)\, < \, 0\) the combination is antagonistic. Additionally, for Figs. 3 and 5 we had to calculate the Hill-dependent bias in Loewe. For Fig. 3, we subtracted the Hill slope bias to only study the impact of conflating synergistic potency and efficacy. To calculate the bias, Loewe was calculated as in Eq. (5) where *H**x*_{inv} was was evaluated at the MuSyC-predicted *E**d*(*d*1, *d*2) for the combination retaining all the single drug parameters (*E*_{0}, *E*_{1}, *E*_{2}, *C*_{1}, *C*_{2}, *h*1, *h*2) and assuming (*α*_{12} = *α*_{21} = 0). This resulted in an estimate of the bias purely due to the Hill slope in the Loewe calculation.

#### ZIP and BRAID

Both ZIP and BRAID were calculated for each combination using the R packages available for each method: (ZIP’s R code is in the supplemental file 1 of the manuscript^{7} and BRAID’s package is available from: https://cran.r-project.org/web/packages/braidReports/braidReports.pdf.

#### Effective dose model

To fit Zimmer et al.’s Effective Dose Model we used the scipy.optimization.curve_fit module in Python 2.7. Specifically, the Effective Dose Model, Eq. 2 in Zimmer et al.^{50} (Eq. 30 in Supplement), contains parameters (*C*_{1}, *C*_{2}, *a*_{12}, *a*_{21}, *h*_{1}, *h*_{2}) where the *a* parameters are the synergy values. 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 (Fig. S7A, B).

#### Schindler’s Hill PDE model

The Hill PDE model has no parameters to fit as the dose-response surface is derived the single dose-response curves. 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^{8}. Therefore, to calculate the synergy of this model, we defined the sum of residuals between the null surface and the experimental data as the metric of synergy.

#### Combination Index

Following the CI fitting algorithm in the CompuSyn software, we fit the median-effect equation, a 2-parameter, log-linearized form of the Hill equation to each drug alone obtained by assuming *E*_{0} = 1 and *E*_{1} = 0. *C* and *h* were then obtained from the slope and y-intercept of the log-linearized data using the scipy.stats.linregress module in Python 2.7. While CI assumes the drug effect is scaled between (0, 1), when this is not the case, a poor fit results (Figs. 6 and S7C, D). All data points with percent viability >1 were excluded from the CI fit because the median-effect equation becomes complex in those cases. For some drugs, this left too few points to fit a line, such that CI was undefined. In other cases, the fitted value for *h* was <0, a physically unrealistic result; therefore, those combinations were also considered undefined. After that, CI was calculated using Eq. (5).

As with Loewe, we apply a \(-{{{{\mathrm{log}}}}\,}_{10}\) transformation to scale CI synergy such that \(-{{{{\mathrm{log}}}}\,}_{10}({{{\mathrm{CI}}}})\; > \; 0\) the combination is synergistic and if \(-{{{{\mathrm{log}}}}\,}_{10}({{{\mathrm{CI}}}})\, < \, 0\) the combination is antagonistic.

#### GPDI model

The GPDI model fitting algorithm was developed in Python by the authors based on description from^{24}. Fits for the single drug parameters (*C*1, *C*2, *h*1, *h*2, *E*1, *E*2) were based on the single dose-response data alone (See Table S5 for parameter translation). Fitting for the parameters (*I**N**T*_{21}, *I**N**T*_{12}, *C*_{INT,12}, *C*_{INT,21}) was done using *c**u**r**v**e*_*f**i**t* function in python in either the Loewe or Bliss version of GPDI (see Supplemental Section 4.5). For all conditions (*h*_{INT,12} = *h*_{INT,21} = 1) as was asserted in the original paper.

### Drug combination database analysis

Initial possible matching drug names between the in vitro experiments and the Drug Combination Database (DCDB) were determined using fuzzy string matching in Python (https://github.com/seatgeek/fuzzywuzzy v0.17.0). Drugs which had a sorted token ratio of 85 were initially included. Of the 427 drugs, there were 172 single drug matches. Matches included structural analogs. See *m**a**t**c**h**i**n**g*_*d**r**u**g*_*n**a**m**e**s*−11−29−2019_*f**i**n**a**l*.*c**s**v* for complete matching list. Of these matching drugs, there were 126 tested combinations in clinical trials according to DCDB. Outliers in the synergy calculation were considered values 1.5 times the interquartile range (*Q*1−*Q*3) above or below Q1 or Q3 respectively.

### Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

## Data availability

The datasets analyzed in this study were obtained from publicly available sources with DOIs Mott (10.1038/srep13891, Figs. 3, 6, S4, S6, and S9), O’Neil (10.1158/1535-7163.MCT-15-0843, Figs. 3–4 and S4–S9), Holbeck (10.1158/0008-5472.CAN-17-0489, Figs. 7, S6, and S9), Tan (10.1038/nbt.2391, Fig. S6), and Cokol (10.1038/msb.2011.71, Figs. S6 and S9). Clinical trial data was collected from the Drug Combination Database (DCDB) (10.1093/database/bau124, http://public.synergylab.cn/dcdb/, Figs. S9–S10). The synergy datasets generated in this study are available in a repository at https://bitbucket.org/meyerct1/musyc_theory/.

## Code availability

All code required for recreating manuscript analyses from the MuSyC fits are available for review in the repository https://bitbucket.org/meyerct1/musyc_theory/. A web application to calculate MuSyC parameters is available at https://musyc.lolab.xyz/. The code for the interactive Jupyter notebook demonstration of MuSyC is available at https://github.com/djwooten/natcomms-musyc2021/.

## References

- 1.
Loewe, S. über Kombination swirkungen.

*Arch. Exp. Pathol.***114**, 313–326 (1926). - 2.
Loewe, S. Versuch einer allgemeinen Pharmakologie der Arznei- kombinationen.

*Klin. Wochenschr.***6**, 1078–1085 (1927). - 3.
Bliss, C. I. The toxicity of poisons applied jointly.

*Ann. Appl. Biol.***26**, 585–615 (1939). - 4.
Greco, W. et al. Consensus on concepts and terminology for combined-action assessment: the saariselka agreement.

*ACES***4**, 65–69 (1992). - 5.
Tang, J., Wennerberg, K. & Aittokallio, T. What is synergy? The Saariselkä agreement revisited.

*Front. Pharmacol.***6**, 181 (2015). - 6.
Zimmer, A., Katzir, I., Dekel, E., Mayo, A. E. & Alon, U. Prediction of multidimensional drug dose responses based on measurements of drug pairs.

*Proc. Natl Acad. Sci. USA***113**, 10442–7 (2016). - 7.
Yadav, B., Wennerberg, K., Aittokallio, T. & Tang, J. Searching for drug synergy in complex dose-response landscapes using an interaction potency model.

*Comput. Struct. Biotechnol. J.***13**, 504–513 (2015). - 8.
Schindler, M. Theory of synergistic effects: Hill-type response surfaces as ‘null-interaction’ models for mixtures.

*Theor. Biol. Med. Model.***14**, 15 (2017). - 9.
Foucquier, J. & Guedj, M. Analysis of drug combinations: current methodological landscape.

*Pharmacol. Res. Perspect.***3**, e00149 (2015). - 10.
Geary, N. Understanding synergy.

*AJP: Endocrinol. Metab.***304**, E237–E253 (2013). - 11.
Twarog, N. R., Stewart, E., Hammill, C. V. & A. Shelat, A. BRAID: a unifying paradigm for the analysis of combined drug action.

*Sci. Rep.***6**, 25523 (2016). - 12.
Meyer, C. T., Wooten, D. J., Lopez, C. F. & Quaranta, V. Charting the fragmented landscape of drug synergy.

*Trends Pharmacol. Sci.***41**, 266–280 (2020). - 13.
Ianevski, A., He, L., Aittokallio, T. & Tang, J. SynergyFinder: a web application for analyzing drug combination dose-response matrix data.

*Bioinformatics***33**, 2413–2415 (2017). - 14.
Flobak, A., Vazquez, M., Lægreid, A. & Valencia, A. CImbinator: a web-based tool for drug synergy analysis in small- and large-scale datasets.

*Bioinformatics***33**, 2410–2412 (2017). - 15.
He, L. et al. Methods for high-throughput drug combination screening and synergy scoring.

*Methods Mol. Biol.***1711**, 351–398 (2018). - 16.
Greco, W. R., Bravo, G. & Parsons, J. C. The search for synergy: a critical review from a response surface perspective.

*Pharmacol. Rev.***47**, 331–85 (1995). - 17.
Bansal, M. et al. A community computational challenge to predict the activity of pairs of compounds.

*Nat. Biotechnol.***32**, 1213–22 (2014). - 18.
Menden, M. P. et al. Community assessment of cancer drug combination screens identifies strategies for synergy prediction.

*bioRxiv*https://doi.org/10.1101/200451 (2017). - 19.
Han, K. et al. Synergistic drug combinations for cancer identified in a CRISPR screen for pairwise genetic interactions.

*Nat. Biotechnol.***35**, 463–474 (2017). - 20.
Palmer, A. C. & Sorger, P. K. Combination cancer therapy can confer benefit via patient-to-patient variability without drug additivity or synergy theory combination cancer therapy can confer benefit via patient-to-patient variability without drug additivity or synergy.

*Cell***171**, 1678–1691.e13 (2017). - 21.
Meyer, C. T. et al. Quantifying drug combination synergy along potency and efficacy axes.

*Cell Syst.***8**, 97–108 (2019). - 22.
Chou, T.-C. et al. Analysis of combined drug effects: a new look at a very old problem.

*Trends Pharmacol. Sci.***4**, 450–454 (1983). - 23.
Gaddum, J.

*Pharmacology*. (Oxford University Press, London, 1940). - 24.
Wicha, S. G., Chen, C., Clewe, O. & Simonsson, U. S. A general pharmacodynamic interaction model identifies perpetrators and victims in drug interactions.

*Nat. Commun.***8**, 2129 (2017). - 25.
Mott, B. T. et al. High-throughput matrix screening identifies synergistic and antagonistic antimalarial drug combinations.

*Sci. Rep.***5**, 13891 (2015). - 26.
O’Neil, J. et al. An unbiased oncology compound screen to identify novel combination strategies.

*Mol. Cancer Ther.***15**, 1155–1162 (2016). - 27.
Holbeck, S. L. et al. The National Cancer Institute ALMANAC: a comprehensive screening resource for the detection of anticancer drug pairs with enhanced therapeutic activity.

*Cancer Res.***77**, 3564–3576 (2017). - 28.
Tan, X. et al. Systematic identification of synergistic drug pairs targeting HIV.

*Nat. Biotechnol.***30**, 1125–1130 (2012). - 29.
Cokol, M. et al. Systematic exploration of synergistic drug pairs.

*Mol. Syst. Biol.***7**, 544 (2011). - 30.
Chou, T.-C. Drug combination studies and their synergy quantification using the chou-talalay method.

*Cancer Res.***70**, 440–446 (2010). - 31.
Chou, T.-C. The combination index (CI < 1) as the definition of synergism and of synergy claims.

*Synergy***7**, 49–50 (2018). - 32.
Vlot, A. H., Aniceto, N., Menden, M. P., Ulrich-Merzenich, G. & Bender, A. Applying synergy metrics to combination screening data: agreements, disagreements and pitfalls.

*Drug Discov. Today***24**, 2286–2298 (2019). - 33.
Baeder, D. Y., Yu, G., Hoze, N., Rolff, J. & Regoes, R. R. Antimicrobial combinations: bliss independence and loewe additivity derived from mechanistic multi-hit models.

*Philos. Trans. R. Soc. B: Biol. Sci.***371**, 20150294 (2016). - 34.
Chou, T. C. & Talalay, P. Quantitative analysis of dose-effect relationships: the combined effects of multiple drugs or enzyme inhibitors.

*Adv. Enzyme Regul.***22**, 27–55 (1984). - 35.
Wang, H., Kubica, N., Ellisen, L. W., Jefferson, L. S. & Kimball, S. R. Dexamethasone represses signaling through the mammalian target of rapamycin in muscle cells by enhancing expression of REDD1.

*J. Biol. Chem.***281**, 39128–34 (2006). - 36.
Ashton, J. C. Drug combination studies and their synergy quantification using the chou-talalay method–letter.

*Cancer Res.***75**, 2400–2400 (2015). - 37.
Chou, T.-C. Theoretical basis, experimental design, and computerized simulation of synergism and antagonism in drug combination studies.

*Pharmacol. Rev.***58**, 621–681 (2006). - 38.
Fallahi-Sichani, M., Honarnejad, S., Heiser, L. M., Gray, J. W. & Sorger, P. K. Metrics other than potency reveal systematic variation in responses to cancer drugs.

*Nat. Chem. Biol.***9**, 708–714 (2013). - 39.
Wood, K., Wood, K., Nishida, S. & Cluzel, P. Uncovering scaling laws to infer multidrug response of resistant microbes and cancer cells.

*Cell Rep.***6**, 1073–1084 (2014). - 40.
Tekin, E. et al. Prevalence and patterns of higher-order drug interactions in Escherichia coli.

*npj Syst. Biol. Appl.***4**, 31 (2018). - 41.
Di Veroli, G. Y. et al. An automated fitting procedure and software for dose-response curves with multiphasic features.

*Sci. Rep.***5**, 14701 (2015). - 42.
Jones, E., Oliphant, T. & Peterson, P.

*SciPy: Open Source Scientific Tools for Python*(2001). - 43.
Oliphant, T. E.

*Guide to NumPy*(2006). http://web.mit.edu/dvp/Public/numpybook.pdf. - 44.
McKinney, W.

*Data Structures for Statistical Computing in Python*(2010). http://conference.scipy.org/proceedings/scipy2010/mckinney.html. - 45.
Hunter, J. D. Matplotlib: a 2D graphics environment.

*Comput. Sci. Eng.***9**, 90–95 (2007). - 46.
Lebigot, E. O.

*Uncertainties: a Python package for calculations with uncertainties*(2011). http://pythonhosted.org/uncertainties/. - 47.
Motulsky, H. & Christopoulos, A.

*Fitting Models to Biological Data using Linear and Nonlinear Regression A practical guide to curve fitting Contents at a Glance*(GraphPad Software Inc., San Diego, 2003). http://www.facm.ucl.ac.be/cooperation/Vietnam/WBI-Vietnam-October-2011/Modelling/RegressionBook.pdf. - 48.
Malyutina, A. et al. Drug combination sensitivity scoring facilitates the discovery of synergistic and efficacious drug combinations in cancer.

*PLoS Comput. Biol.***15**, 1–19 (2019). - 49.
Demidenko, E. & Miller, T. W. Statistical determination of synergy based on bliss definition of drugs independence.

*PLoS ONE***14**, 1–22 (2019). - 50.
Zimmer, A., Katzir, I., Dekel, E., Mayo, A. E. & Alon, U. Prediction of multidimensional drug dose responses based on measurements of drug pairs.

*Proc. Natl Acad. Sci. USA***113**, 10442–7 (2016). - 51.
Elstrodt, F. et al. Brca1 mutation analysis of 41 human breast cancer cell lines reveals three new deleterious mutants.

*Cancer Res.***66**, 41–45 (2006). - 52.
Chou, T. C. Relationships between inhibition constants and fractional inhibition in enzyme catalyzed reactions with different numbers of reactants, different reaction mechanisms, and different types and mechanisms of inhibition.

*Mol. Pharmacol.***10**, 235–47 (1974). - 53.
Harris, L. A. et al. An unbiased metric of antiproliferative drug effect in vitro.

*Nat. Methods***13**, 497–500 (2016). - 54.
Weiss, J. N. The Hill equation revisited: uses and misuses.

*FASEB J.***11**, 835–841 (1997). - 55.
Chou, T. C. & Talalay, P. Generalized equations for the analysis of inhibitions of Michaelis-Menten and higher-order kinetic systems with two or more mutually exclusive and nonexclusive inhibitors.

*Eur. J. Biochem.***115**, 207–16 (1981).

## Acknowledgements

The authors would like to thank Corey Hayford, Sarah Groves, Darren Tyson, Leonard Harris, Joshua Bauer, James Pino, Ken Lau, and Chris Wright for insightful conversations and critical feedback. The authors would also like to thank Monica Del Valle for improvements to the software. This work was supported by the following funding sources: CTM was supported by National Science Foundation (NSF) Graduate Student Fellowship Program (GRFP) [Award #1445197]; C.F.L. was supported by the National Science Foundation [MCB 1411482 and CAREER 1942255]; C.F.L. and V.Q. were supported by the National Institutes of Health (NIH) [U54-CA217450 and U01-CA215845]; and V.Q. was supported by NIH [R01-186193].

## Author information

### Affiliations

### Contributions

D.J.W., C.T.M., C.F.L., and V.Q. conceived the study. C.T.M. and D.J.W. designed and executed the analyses. All authors reviewed the final manuscript; C.T.M. and D.J.W. performed the visualizations; A.L.R.L. implemented the webserver for MuSyC calculations. C.F.L., C.T.M., and D.J.W. wrote the manuscript.

### Corresponding authors

## Ethics declarations

### Competing interests

C.T.M. and V.Q. are academic founders and part equity holders in Parthenon Therapeutics Inc. All other authors declare no competing interests.

## Additional information

**Peer review information** *Nature Communications* thanks Åsmund Flobak and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Publisher’s note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

## About this article

### Cite this article

Wooten, D.J., Meyer, C.T., Lubbock, A.L.R. *et al.* MuSyC is a consensus framework that unifies multi-drug synergy metrics for combinatorial drug discovery.
*Nat Commun* **12, **4607 (2021). https://doi.org/10.1038/s41467-021-24789-z

Received:

Accepted:

Published:

## Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.