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 massaction 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 Hillslope and maximal effect of singledrugs. 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 NCIDREAM consortia^{17,18}, as well as combinatorial CRISPR screens^{19}. Yet, the paucity of successful clinical combinations explicable by true pharmacological interaction, rather than patienttopatient 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 Multidimensional 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) MSPbased frameworks are biased toward antagonism for drugs with intermediate efficacy; and (3) DEPbased frameworks contain a Hillslope dependent bias. The Hillslope bias results from satisfying the famous “sham” combination thought experiment, thus arguing against the merit of shamcompliance as a measure of validity for synergy frameworks. Using five large combination datasets^{25,26,27,28,29}, MuSyC identifies realworld 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 Actionbased, 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 statetransition model to measure multidrug synergistic effects
The 4parameter Hill equation is commonly used to fit doseresponse 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 twostate 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 doseresponse curve (Fig. 1A, right). To generalize this onedrug formalism to two concurrent drugs, we developed a fourstate transition model of combination pharmacology (Fig. 1B, left)^{21}. From this model, we derive a twodimensional (2D) Hill equation for two drugs (Box 1, Eq. (15)) defining a doseresponse surface (Fig. 1B, middle). The 2D Hill equation contains five additional parameters, not present in the singledrug 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 growthrate data as the sum of growthrate 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}. DEPbased 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) doseresponse 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 MSPderived 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 Hillslope 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 Reexamining the sham experiment: Sham compliance introduces Hilldependent 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 drugresponse data, we generated synthetic doseresponse 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 anticancer 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 (BRCA1mutant ovarian carcinoma) cells is antagonistically efficacious and potent by MuSyC, but synergistic by Loewe (Fig. 3C).
Some assays, such as the antimalarial combination screen from Mott et al.^{25} are designed so that the drug effect spans the entire range from E_{0} = 1 to E1 = E2 = E3 = 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 antimalarial dataset^{25}, Bliss is consistently synergistic when log(α_{12}, α_{21}) > 0, and antagonistic if log(α_{12}, α_{21}) < 0; however, when log(α_{12}) < 0 < log(α_{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 multidrug 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 doseresponse surface (compare \(({{{\mathrm{log}}}}\,({\alpha }_{12}),{{{\mathrm{log}}}}\,({\alpha }_{21}),\beta )=(,,+)\) and ( − , − , − ) quadrants of Fig. 3G and Video S2). In the anticancer 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 mk8669 (PI3K/mTOR dual inhibitor) in a colorectal adenocarcinoma cellline 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, doseresponse 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 anticancer dataset is biased toward antagonism for moderately efficacious combinations 0.35 < (E1, E2) < 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 pancancer trends in the combination of drugs which have, on average, intermediate effect over the entire cellline panel (Fig. 4D). As a particular example, the synergistic efficacy of paclitaxel (targets microtubule stability) and mk2206 (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 MuSyCpredicted bias was sufficient to obfuscate many of the strongest synergies and antagonisms according to Bliss (Fig. S6A, B).
Other MSPbased methods, such the Effective Dose model also assume data measures percent affect and fit a simplified 2parameter 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.
Reexamining the sham experiment: Sham compliance introduces Hilldependent 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 mixedinhibition (black circles, Fig. 5A, B). In sham combinations, these mixedinhibition states are equivalent to single drug completeinhibition states (Fig. 5A, cyan circles), while for real combinations, these are not equivalent (Fig. 5B). When h = 1, these intermediate mixedinhibition 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 anticancer dataset, the average trend of Loewe synergy closely follows the Hill slope bias predicted by MuSyC (Fig. 5C). Further, subtracting the MuSyCpredicted 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 mk4827 (PARP inhibitor) have small Hill slopes when applied to MBAMB436 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 MuSyCpredicted bias revealed a screenwide 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 Hillslope and efficacy range bias, like those we reported above for DEP and MSPbased frameworks, respectively. Like MuSyC, CI is based on the Law of MassAction, 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 Hillslope dependent bias we show in Fig. 5. We show the consequence of this bias in illustrative twodrug 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 singledrug doseresponse curves using the twoparameter medianeffect equation^{30} (Eq. (8)), see also Supplemental Section Percent Affect vs Percent Effect—Combination Index). In contrast, MuSyC is based on a fourparameter Hill equation (Box 1, Eq. (10)), and thus can describe metrics of drug effects with arbitrary ranges. It is common to observe in many doseresponse assays, such as percent viability, drug effects that do not reach 0% (i.e. E_{1} > 0)^{38}. For such assays, the twoparameter medianeffect equation fits poorly (Fig. 6D). These poor fits lead to an effectdependent error in CI quantification of synergy as observed in synthetic (Fig. 6E) and experimental combinations (Fig. 6F). We note, this effectdependent error is not the same as the effectrange bias we report for MSPbased 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 statetransition model results in a single framework unifying both principles. By mapping all frameworks on a common landscape, MuSyC facilitates rigorous investigation of oftcited, 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 doseresponse 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 inbetween, so no choice is required. Further, while MuSyC’s massaction formalism applies most directly to molecular inhibition, its synergy parameters describe geometric transformations of efficacy, potency, and Hill slope, which are wellestablished 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 10fold 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, C2), (d_{1,max}, C2) 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 onoff 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 nonefficacious combinations (see Supplemental Section MuSyC statistically distinguishes efficacious and nonefficacious 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 batcheffects 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 Hilldependent bias.
CI has previously been criticized for its procedure to fit the medianeffect equation. Specifically, CI depends on a logtransformation in order to linearize the twoparameter median effect equation “in similar logic to the defunct Scatchard analysis in pharmacology, which has been replaced by nonlinear modeling”^{36}. That is, this logtransformation 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 nonlinear ways, the expected error when applying CI is unique to each particular combination and assay design.
The prospects of higherorder 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 higherorder 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 higherorder interactions than combinations all sharing a single synergy type. However, the number of synergy parameters in MuSyC scales as 2^{N}(N + 1) − 3N − 1 (including γ) where N is the number of drugs (Fig. S1), and the commensurate data necessary to fully constrain MuSyC hypersurfaces 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 stepwise from the boundaries (see Supplemental Section Proof of boundary behavior of the 2D Hill equation).
MuSyC expects singledrug doseresponse curves to be sigmoids well fit by the 1D Hill equation (Eq. (10)), and doseresponse 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 anticancer and antimalarial datasets had R^{2} > 0.7, respectively). However, it is possible for drugs to have multiphasic responses due to polypharmacology 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 nonparametric frameworks such as Bliss, Loewe, or HSA. Nevertheless we note that without fitting doseresponse curves to a parametric model, these metrics are sensitive to noise in individual datapoints. 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 massaction, 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 commaseparated value (CSV) format; data format details and usage instructions are in the supplemental materials. The application fits doseresponse surfaces using MuSyC and offers the results both as a CSV download of fit parameters, and interactive plots of the doseresponse 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 taskworker 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., antiproliferative 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 nonlinear 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 multicollinearity 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 4parameter 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 colinear. 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 multicollinearities described above. The Monte Carlo algorithm for fitting the 2D Hill equation is as follows. First, the 4parameter 1D Hill equation (Eq. (10)) is fit to the doseresponse of each drug in isolation. The fit uses the Trust Region Reflective (TRF) algorithm implemented in the curve_fit() 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}, h1, h2), 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 bestfit values of the 2D Hill equation for all drug doses. The data plus noise is then fit again initializing the optimization from the bestfit 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 nonGaussian distribution of uncertainty for many fits (e.g. the distribution of h is lognormal) as well as being robust to the colinear parameters in the 2D Hill equation. The asymmetric confidence interval is particularly salient when the doserange is insufficient to observe the lower plateau of the doseresponse. 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 Zscore 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 doseresponse 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. anticancer screen
The anticancer 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. antimalarial screen
The antimalarial drug combination data were downloaded from https://tripod.nih.gov/matrixclient/ 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. antiHIV screen
The antiHIV 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. antifungal screen
The antifungal drug combination data were downloaded from supplemental dataset one in Cokol et al.^{29}. The raw cell growth measurements for all 200 drugdrug 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. anticancer screen
The ALMANAC anticancer drug combination data were downloaded from https://wiki.nci.nih.gov/display/NCIDTPdata/NCIALMANAC file ComboDrugGrowth_Nov2017. zip. The matching single doseresponse data were downloaded from https://wiki.nci.nih.gov/display/NCIDTPdata/NCI60+Growth+Inhibition+Data, June 2016 release DOSE_RESPONSE 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 Ed(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 Hilldependent 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 Hx_{inv} was was evaluated at the MuSyCpredicted Ed(d1, d2) for the combination retaining all the single drug parameters (E_{0}, E_{1}, E_{2}, C_{1}, C_{2}, h1, h2) 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.rproject.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 doseresponse surface is derived the single doseresponse 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 medianeffect equation, a 2parameter, loglinearized 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 yintercept of the loglinearized 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 medianeffect 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 (C1, C2, h1, h2, E1, E2) were based on the single doseresponse data alone (See Table S5 for parameter translation). Fitting for the parameters (INT_{21}, INT_{12}, C_{INT,12}, C_{INT,21}) was done using curve_fit 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 matching_drug_names−11−29−2019_final.csv 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 (Q1−Q3) 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/15357163.MCT150843, Figs. 3–4 and S4–S9), Holbeck (10.1158/00085472.CAN170489, 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/natcommsmusyc2021/.
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 combinedaction 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 doseresponse landscapes using an interaction potency model. Comput. Struct. Biotechnol. J. 13, 504–513 (2015).
 8.
Schindler, M. Theory of synergistic effects: Hilltype response surfaces as ‘nullinteraction’ 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 doseresponse matrix data. Bioinformatics 33, 2413–2415 (2017).
 14.
Flobak, A., Vazquez, M., Lægreid, A. & Valencia, A. CImbinator: a webbased tool for drug synergy analysis in small and largescale datasets. Bioinformatics 33, 2410–2412 (2017).
 15.
He, L. et al. Methods for highthroughput 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 patienttopatient variability without drug additivity or synergy theory combination cancer therapy can confer benefit via patienttopatient 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. Highthroughput 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 choutalalay 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., UlrichMerzenich, 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 multihit models. Philos. Trans. R. Soc. B: Biol. Sci. 371, 20150294 (2016).
 34.
Chou, T. C. & Talalay, P. Quantitative analysis of doseeffect 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 choutalalay 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.
FallahiSichani, 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 higherorder 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 doseresponse 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/WBIVietnamOctober2011/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 MichaelisMenten and higherorder 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) [U54CA217450 and U01CA215845]; and V.Q. was supported by NIH [R01186193].
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 multidrug synergy metrics for combinatorial drug discovery. Nat Commun 12, 4607 (2021). https://doi.org/10.1038/s4146702124789z
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s4146702124789z
Further reading

Nonparametric synergy modeling of chemical compounds with Gaussian processes
BMC Bioinformatics (2022)
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.