MuSyC is a consensus framework that unifies multi-drug synergy metrics for combinatorial drug discovery

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.

: Extension of MuSyC to combinations of three drugs. A) Following the cubic geometry of Figure 1B (main text), combinations of 3 drugs result in a cube. Numbered notation next to each edge corresponds to the ratio of the connected corners at equilibrium for the boundary conditions. For example, edge #10 annotation means A123 A12 → C3 α23d3 γ312h3 as doses d 1 → inf, d 2 → inf. For combinations of 4 drugs, the geometry is a tesseract. In general, MuSyC can describe combinations N -drug combinations by considering 2 N possible states with transitions defining the edges of an N -dimensional hypercube. Dose-response surfaces generalize to N -dimensional scalar functions. In the most general case, for N drugs there are 2 N − N − 1 distinct β parameters (one for each state characterized by the action of at least 2 drugs), and N · 2 N −1 − 1 distinct α and γ parameters (one for each edge, excluding edges connected to the undrugged state, which correspond only to single-drug potency and cooperativity). Thus, MuSyC can account for higher-order synergies (e.g., synergy that emerges from a combination of three drugs, but is not evident in any pairwise combination of those drugs), however the rapid growth of the number of synergy parameters with N suggests that significant quantities of data, or confident knowledge of pairwise synergies, would be needed to measure such higher-order synergies. The linear isoboles of Loewe result from asserting that for a given effect E-achievable either by dose d 1 of Drug 1 alone, or dose d 2 of Drug 2 alone-there is a constant ratio R = d1 d2 such that using ∆d 2 less of Drug 2 can always be compensated for with ∆d 1 = R∆d 2 more of Drug 1 to achieve the same effect. Synergy (Syn, red) occurs when less of one drug is required to compensate for a decrease in the other than expected. Likewise, antagonism (Ant, blue) occurs when more of one drug is required. B) MuSyC's isoboles become linear under the constraint that the drugs are mutually exclusive (α 12 = α 21 = 0) and h 1 = h 2 = 1. Additionally, the slope of the line (R) in MuSyC is equal to the ratio of the two drugs' EC50 (− C2 C1 ) (See eq. 27 in Supplemental Section Relationship between MuSyC and the MSP and DEP). C) When h < 1, the MuSyC isoboles bend inward. As a result, DEP-based frameworks (e.g. Loewe and CI), which assume linear isoboles, will quantify the red shaded region (between the straight, dotted, diagonal line, and the curved line below it) as synergistic thereby overestimating the synergy. D) Conversely, when h > 1, MuSyC isoboles bend outward causing DEP-based frameworks to overestimate antagonism. For clarity, panels C and D have been replicated in Figures 6A Figure S3: Fit algorithm robustness to different noise profiles. A) Varying conditions considered included percent noise, noise type, sampling design, data density, and synergy profile. Each condition was randomly sampled 10 times corresponding to 54,000 synthetic datasets. For sample design, the data density was distributed evenly such that a Full Matrix sample with 100 conditions has 10 single doses per drug while the CSS or the Edge sampling with 100 conditions has 25 single dose samples (100/4 lines). All surfaces have the following single drug parameters (h1 = h2 = 1,E0 = 1,E1 = E2 = 0.334,C1 = C2 = 1e − 5). Doses were sampled between (1e − 10,1) using a log sampling. B) Percent of conditions for which the fit converged out of 54,000. Only parameterized models of drug synergy were compared. Fits with R 2 < 0.3 were considered to not converge. MuSyC and ZIP converge most frequently (99%) while CI and Effective Dose (Eff. Dose) converge less frequently as they both assume the maximal and minimal effect of the drug is (0,1). C) Distributions of R 2 values for each method. ZIP and CI synergy are calculated as the deviation at each point from the null and therefore have no R 2 value. MuSyC has the highest mean R 2 value (0.95) followed by BRAID and GPDI(Loewe) at (0.94). D) Average Z-score between each 10 matched random samples is used as a measure of robustness. Lower average Z-score indicates less variation in synergy calculation as a result of noise. MuSyC fits are not significantly impacted by the density of the data (first column). CSS sampling increases the uncertainty of β (second column) which is expected because the point used to define β (max(d1),max(d2)) is absent. Overall, the Effective Dose (Eff. Dose) was the most robust. For all methods, heteroskedastic distribution of noise is better fit than homoskedastic, a difference that becomes more exaggerated at higher levels of noise.  Figure S4: Model selection prefers simpler MuSyC model. Difference in the Akaike Information Criterion (AIC) values for MuSyC models including fitting synergistic cooperativity (γ 12 and γ 21 ) or fixing γ to 1 thereby reducing the parameter count by 2. Models which minimize AIC are preferred. In most cases the simpler model is preferred. The mean AIC γ − AIC noγ for anti-cancer and anti-malarial datasets was 8 and 3312, respectively. The percent of combinations for which the model including γ had a lower AIC value was 18% and 5%, respectively.   Figure S6: Impact of bias in Bliss and Loewe. A) Schematic for interpreting bias figures. Synergy as calculated by Bliss or Loewe is plotted on the x-axis. Subtracting the MuSyC-estimated bias (calculated for each data point) results in a corrected value of Bliss or Loewe synergy, plotted on the y-axis. We consider efficacy-bias for Bliss (see Figure 4B, main text) and Hill-slope bias for Loewe (see Figure 5C, main text). Points along the line y = x have no bias. Points above the line indicate the uncorrected metric underestimated synergy. Points below the line indicate the uncorrected metric overestimated synergy. B) Bliss efficacy-bias plots. Bliss cannot be calculated for the Cokol and Holbeck datasets due to their effect metrics (Table 3, main text). There is substantial deviation from the red line indicating the magnitude of the efficacy-bias is comparable to the magnitude of uncorrected Bliss. C) Loewe Hill-bias plots. While, there is less deviation from the red line than Bliss, we observe whole datasets (e.g., Cokol Figure S7: Drugs with weak monotherapy maximum effect can cause problems for some synergy models. A) Dose-response surface for combination erlotinib and abt-888 in SW837 cells 1 . B) MuSyC parameter fits showing weak E max = E 2 = 0.9 for abt-888. C) Dose response surface according to Effective Dose model, which enforces E 0 = 1 and E 1 = E 2 = 0, which leads to poor fit to data.  Figure S8: Small Hill-slope biases Combination Index toward synergy. A) MuSyC predicts nonlinear isoboles (bottom) when h < 1. Here, CI would erroneously assess dose pairs in the red shaded region (between the straight, dotted diagonal line, and the solid curve beneath it) as synergistic. B) An example, synthetic two-drug dataset where this Hill dependent bias is apparent (h 1 = h 2 = 0.5,E 0 = 1, E 1 = E 2 = 0.0, C 1 = C 2 = 1, α 12 = α 21 = 0). This trend is most apparent near d 1 ≈ d 2 . C) MuSyC fit of the combination of mk-8776 (CHK1 inhibitor) and lapatinib (EGFR inhibitor) in ZR751 (luminal B breast cancer) 1 . D) Color indicates CI of the combination from (C), and shows the predicted CI bias (synergy along diagonal due to h < 1). Because E 0 = 1,E 1 = E 2 = 0, and (α 12 ,α 21 )<0.0001, the difference between MuSyC and CI is due to the Hill slope error in CI.  Figure S9: MuSyC distinguishes Non-efficacious and Efficacious clinical combinations based on in vitro screening data. A) Distributions of synergistic efficacy (β), synergistic potency (log(α 12/21 )), and Loewe additivity from in vitro measurements of clinical combinations in DCDB. Data were aggregated from 1-4 , see Figure S10 for details. Only combinations with two drugs were considered. p-values were calculated by one-sided t-test (null hypothesis green > purple) after outliers were removed. The white marks indicate the means, and black boxes indicate the 25 th and 75 th percentiles. Sample count for p-value calculations for β, log(α 12/21 ), and Loewe are (n=2538, 5266, 1198) efficacious combinations (purple) and (n=796, 1620, 315) non-efficacious combinations (green), respectively. B) Combination of gemcitabine (nucleoside analogue) and docetaxel (tubulin stabilizer) in prostate cancer PC3 cells (ALMANAC dataset 2 ). C) The combination in (B) is antagonistically efficacious (β=-0.31), but synergistic by Loewe due in part to the Hill bias ( √ h 1 · h 2 < 1, Figure 5, main text).   Figure S10: Clinical drug combinations with matching in vitro data. A) Five in vitro drug combination datasets (Table 3, main text) were cross-referenced with clinical drug combinations in DCDB. Total there were 126 matching combinations tested in 2,956 samples (n in chart). The ALMANAC screen (holbeck_anticancer) 2 contains the 77% of the 126 combinations. B) Of the clinically-tested drug combinations, 46% reached a Phase 3 or Phase 4 clinical trial and over 85% were at least Phase 1. C) Clinical combinations are annotated as efficacious according to DCDB. D) Distribution of Bliss, HSA, and CI for the efficacious and non-efficacious combinations. Bliss and CI cannot be calculated for any of the holbeck_anticancer or cokol_antifungal due to the measured drug effects not being a percentage in the range (0,1) resulting in only 3 and 1 combinations, respectively, which are non-efficacious. The classic HSA method depends on the assay scale. However, if HSA is normalized by dividing the effect metric range, similar to MuSyC's β, there is a statistical difference (middle panel, p-value calculated by one-sided t test, assuming normality; green>purple). Sample count for p-value calculations for Bliss, HSA, and CI are (n=403, 2497, 265) efficacious combinations (purple) and (n=3, 772, 1) non-efficacious combinations (green), respectively. The white mark shows the median. The thick portion of the vertical bar running through the middle of each violin plot shows the interquartile range, while the thin portion shows the minimum to maximum range.   Figure S12: Demo: Bliss shows antagonism for medium-efficacy drugs. Screenshot of the interactive Jupyter Notebook in File S1, and hosted online at https://mybinder.org/v2/gh/djwooten/natcomms-musyc2021/HEAD?filepath= demo.ipynb. On the left is the MuSyC dose-response surface corresponding to the parameters given by the sliders below. On the right is the Bliss excess synergy landscape. Because E 1 = 0.4 and E 2 = 0.5, Bliss expects the combination to achieve Because of this, Bliss rates this combination as antagonistic at high doses of both drugs. Figure S13: Demo: Synergistic efficacy may be reflected by Bliss at high doses of both drugs. Screenshot of the interactive Jupyter Notebook in File S1, and hosted online at https://mybinder.org/v2/gh/djwooten/ natcomms-musyc2021/HEAD?filepath=demo.ipynb. On the left is the MuSyC dose-response surface corresponding to the parameters given by the sliders below. On the right is the Bliss excess synergy landscape. Because E 1 = 0.4 and E 2 = 0.5, Bliss expects the combination to achieve E 3 = E 1 · E 2 = 0.2. However, because β = 0.5, this combination has E 3 = 0.1. Because of this, Bliss rates this combination as synergistic at high doses of both drugs. In general, adjusting the slider for beta affects both the dose response surface and the Bliss excess landscape at high doses of each drug. Figure S14: Demo: MuSyC can reproduce Bliss null model. Screenshot of the interactive Jupyter Notebook in File S1, and hosted online at https://mybinder.org/v2/gh/djwooten/natcomms-musyc2021/HEAD?filepath=demo.ipynb. On the left is the MuSyC dose-response surface corresponding to the parameters given by the sliders below. On the right is the Bliss excess synergy landscape. Clicking the "Bliss Null" button highlighted in red will adjust the MuSyC synergy sliders to have log(alpha) = log(gamma) = 0 and E 3 = E 1 · E 2 = 0.2 (note: E 3 is controlled by adjusting the beta slider, such that With these settings, as can be seen on the right, the MuSyC dose-response surface matches the Bliss null model, so that this combination has 0 synergy by Bliss at every dose. In general, adjusting the alpha sliders affects both the dose response surface and the Bliss synergy landscape near the EC50 of one drug, and at high doses of the other. Figure S15: Demo: Synergistic potency may be reflected near the EC50 of the drug whose potency is synergistically improved. Screenshot of the interactive Jupyter Notebook in File S1, and hosted online at https: //mybinder.org/v2/gh/djwooten/natcomms-musyc2021/HEAD?filepath=demo.ipynb. On the left is the MuSyC doseresponse surface corresponding to the parameters given by the sliders below. On the right is the Bliss excess synergy landscape. Starting from the Bliss null model (obtained by pressing the "Bliss Null" button), and adjusting the slider for alpha21 highlighted in red, we see that this parameter manifests as Bliss synergy at doses of Drug 1 near its EC50, and at high doses of Drug 2. This is because the parameter alpha21 describes how the presence of Drug 2 (meaning Drug 2 must be administered at high enough doses to be impactful) affects the potency of Drug 1 (meaning this is most apparent near the EC50 of Drug 1). Figure S16: Demo: Synergistic cooperativity may be reflected near the EC50 of the drug whose Hill-slope is synergistically altered. Screenshot of the interactive Jupyter Notebook in File S1, and hosted online at https: //mybinder.org/v2/gh/djwooten/natcomms-musyc2021/HEAD?filepath=demo.ipynb. On the left is the MuSyC doseresponse surface corresponding to the parameters given by the sliders below. On the right is the Bliss excess synergy landscape. Starting from the Bliss null model (obtained by pressing the "Bliss Null" button), and adjusting the slider for gamma21 highlighted in red, we see that this parameter manifests as both Bliss synergy and Bliss antagonism at doses of Drug 1 near its EC50, and at high doses of Drug 2. In general, adjusting the gamma sliders affects both the dose response surface and the Bliss synergy landscape near the EC50 of one drug, and at high doses of the other. MuSyC matches the Bliss null surface when there is no potency synergy (α 12 = α 21 = 1), no cooperativity synergy (γ 12 = γ 21 = 1), and E 3 = E 1 · E 2 (Figure 2A, main text). To show this, let each drug in isolation have a 1D hill response where U i reflects the portion of cells unaffected by drug i alone. For the 2D case, when α 12 = α 21 = 1 and γ 12 = γ 21 = 1, each edge in Figure 1B (main text) satisfies detailed balance and therefore the state occupancy is given by the MuSyC mass action model gives From this, it is easy to verify that is equivalent to the Bliss Independence null model.

35
Furthermore, We define U i , A i , and E i = U i + A i E i to be the fraction of unaffected cells, fraction of affected cells, and observed effect for treatment due to the single drug i, as described by equation 22. The overline distinguishes affects attributable to each drug, such that A 1 includes cells affected either by drug 1 alone, or by both drug 1 and drug 2, while A 1 only includes cells affected by drug 1, but not drug 2 (i.e., From 23, we know U = U 1 · U 2 , and A 1,2 = A 1 · A 2 , leading to Similarly, it is simple to show A 1 = U 2 − U , and similarly for A 2 If E 3 = E 1 · E 2 , then this is equivalent to equation (24). Therefore, given α 12 = α 21 = 1, γ 12 = γ 21 = 1, can always be compensated for by using ∆d 1 = R∆d 2 more of Drug 1 to achieve the same effect 6 . This 45 definition leads to the linear isoboles characteristic of the Loewe null model. 46 Chou and Talalay showed that linear isoboles emerge when the two drugs are mutually exclusive 7 , meaning that the double-drugged state (A 1,2 in Figure 1B, main text) is unreachable. In MuSyC, this requires setting α 12 = α 21 = 0, which reduces the 2D Hill equation (eq. 15, main text) to From this it is easy to see when h 1 = h 2 = 1, the equation describes a straight line line, equivalent to the canonical linear isoboles of Loewe Additivity and the CI null models. Further from equation (25), we find the slope of isoble is equal to − C2 C1 as shown by: Zimmer et. al. 9 introduced the effective dose model as a parameterized extension of Bliss, and to our knowledge were the first to account the asymmetric potency synergy, which is also present in MuSyC. The effective dose model is constructed by fitting the dose response of each single drug to a two-parameter 1D Hill equation in which E 0 and E max are fixed at 1 and 0, respectively To model combination synergy, the authors propose transforming the doses d i to "effective doses" via a 59 system of equations coupling effective doses to one another via a Michaelis-Menten term in the denominator 60 scaled by a synergy parameter a. Original variable name MuSyC variable name The parameter a 12 represents how drug 2 modifies the effective dose synergistically (a 12 < 0) or an-62 tagonistically (a 12 > 0) drug 1. Note that as a 12 → −   Figure S7). In contrast, MuSyC is able to fit dose response surfaces with arbitrary effect ranges.
In contrast to the Effective Dose Model, ZIP, accounts for changes in both the Hill slope and the potency across the dose-response surface. In ZIP, these changes are integrated into a single number (δ), given by ZIP is formulated for arbitrary E 0 and E max ; however, it assumes E max is the same for both drugs, as well 76 as the combination (E 1 = E 2 = E 3 ). To calculate δ, the ZIP method fixes the concentration of one drug, 77 then fits a Hill-equation dose response for the other drug. However, for combinations with efficacy synergy or 78 antagonism, slices of the dose-response surface can have non-Hill, and even non-monotonic shapes. In these 79 cases, ZIP parameter fits may not be meaningful. Because MuSyC accounts explicitly for efficacy synergy, 80 its surfaces are able to describe such complex drug combination surfaces where ZIP cannot.

81
Nevertheless, ZIP parameters µ and η are closely related to MuSyC parameters α and γ. In the absence of synergistic efficacy, slices of MuSyC dose-response surfaces are sigmoidal, though in general do not perfectly follow a Hill equation, and so the ZIP model is still not identical to MuSyC. However, at saturating concentrations of one or the other drug, MuSyC does reduce to the Hill equation. In these saturating cases, ZIP's δ can be related analytically to MuSyC's α and γ by

82
BRAID 12 is an extension of the DEP to effects exceeding the weaker drug and consequently reduces to 83 Loewe under particular conditions (Figure 2, main text). The authors propose three BRAID models with 84 increasing complexity, with eBRAID capable of describing the most general dose-interaction surfaces. We 85 focus our analysis on eBRAID, which assumes that each drug alone has a Hill-like response, and constructs 86 an Hill-like equation for the combination where The BRAID equation (eq. (32)) uses a dose parameter, which combines the doses of both individual drugs,   Original variable name MuSyC variable name

98
These effects can then be inserted into any model of drug effect additivity, such as simple additivity   Highest Single Agent (HSA) 6 is a parsimonious model that defines synergy as the net difference between the combination response and the stronger single-drug response This form assumes that drug decreases E, though it can also be defined for drugs that increase E. At

2D Hill PDE
The most recent framework is one by Schindler 14 which interpolates a null dose-response surface from the  It is generally known that Bliss, and more generally the MSP methods, are only applicable to percentage 147 data. Indeed, percent transformations are commonly applied to data in order to apply Bliss. In Bliss' 148 original study, drug effect was quantified as the percentage of eggs killed (the probability that each egg 149 would die at a given toxin dose), but in all cases the measurement was a discrete event (death of an insect 150 egg). Therefore, the metric of drug effect was the a percentage of affected eggs. However, to apply Bliss to 151 more general measures of drug effect for which discrete counts cannot be obtained, it has been ubiquitous 152 practice to normalize the drug effect as a percent relative to control (e.g., percent viability in cancer research).

153
Nevertheless, such normalization does not, in general, transform measures of efficacy into measures of percent 154 affect.

155
Suppose, as an example, analyzing a drug treatment which actually caused the cells to grow slightly 156 faster than control. By normalizing the drug effect to control, the percent viability is greater than 100% 157 which cannot mean that >100% of the cells were affected.

158
Alternatively, consider the case when a cytostatic drug causes all treated cells to halt both proliferation 159 and death. If the control population doubled twice over 72 hours, the percent viability would be 1 2 3 = 16.25%.

160
If the measure of percent viability was taken instead at 96 hours, the percent viability would be 1 2 4 = 6.25%.

161
At both time points the percent of affected cells was the same (100%); however, the percent of drug effect 162 changes due to normalization.

164
Most DEP frameworks can work with any effect metric or effect range. However, CI does expect data to measure percent affect. CI is defined as 8 where d 1 and d 2 are the doses of drugs 1 and 2, E is the effect (experimentally measured) achieved by the 165 combination of those doses, and f −1 1 (E) and f −1 2 (E) give the dose of Drug 1 or Drug 2 that, when treated 166 as a single-agent, achieve effect E. In CI, E = f (d) is determined by the mass-action derived median-effect 167 equation (Eq. 8, main text). As in Box 1, this leads to the two-parameter Hill equation (Eq. 9, main text)

171
But in drug assays, the measured effect is often not the fraction of (un)inhibited target molecules. Instead, 172 common metrics include population-level readouts such as percent viability, colony formation, or proliferation 173 rate. These data often do not approach 0 as d → ∞, but may instead approach some non-zero maximum 174 effect ( Figure 6A, main text). As shown in Figure 6A (main text), when a drug's maximum effect is non-175 zero, the CI method fits data poorly, which can introduce significant errors in f −1 (E). These errors can be 176 observed in systematic, dose-dependent biases in CI synergy ( Figure 6B, C, main text).  193 In reviewing the literature, we identified a two errors beyond the Hill slope bias ( Figure 5, main text) 194 pertaining to the sham experiment that merit addressing. Further we know This system can be solved to find demonstrating that m 1←2 = m 1 for the sham experiment in the ZIP model, contradicting their proof that 211 ZIP is sham-compliant.

212
To verify ZIP identifies synergy or antagonism for sham combinations, we generated a synthetic sham dose 213 response surface. Sham experiments can be generated exactly for any drug with a pre-defined dose-response 214 by asserting the condition in equation 18 (main text). We constructed a synthetic dataset describing a sham 215 dose response surface for a drug with h = 2, sampled at 2.5 orders of magnitude above and below the EC50.

216
One drug was sampled at 7 concentrations, the other at 12, defining a 7 x 12 sham dose response matrix. 217 We used the synergyfinder 5 R package to calculate synergy by both Loewe and ZIP ( Figure S11), and found 218 Loewe reported close to 0 synergy, as expected for a sham combination, but confirms ZIP detects synergy 219 and antagonism at several concentrations.

223
From the datasets in Table 3 (main text), we identified 126 combinations with clinical trial information 224 annotated in the Drug Combination Database (DCDB) 16 . These combinations were tested in a total of 225 2,996 in vitro disease models ( Figure S1A,B). Of these clinical combinations, 18% were annotated as "Non-226 efficacious" in DCDB, compared to 75% annotated "Efficacious" ( Figure S1C). We found clinically efficacious 227 combinations have higher synergistic efficacy than clinically non-efficacious combinations, based on data from 228 the in vitro screens ( Figure S9A left panel, one-sided t-test p-val=0.002). However, we found no significant 229 difference in synergistic potency between the two sets ( Figure S9A middle, p-val>0.05) pointing to the 230 importance of distinguishing different types of synergy. Loewe ( Figure S9A), Bliss, and CI ( Figure S10D) 231 were unable to distinguish clinically efficacious combinations. Only HSA, which is related to synergistic 232 efficacy ( Figure 2C, main text), is statistically higher in efficacious combinations ( Figure S10).

233
As an example, the combination of gemcitabine and docetaxel for the treatment of prostate cancer was 234 tested in two Phase II clinical trials, where it did not improve overall survival over docetaxel alone 17,18 .

235
Loewe rates this combination as synergistic ( Figure S9B), while MuSyC finds this combination to be an- When d 1 → ∞, then the 4 state reduces to 2 states transition model between A1 and A12.
At equilibrium When A 1 = A 1,2 this is the EC50 for drug 2 given saturating concentrations of drug 1 (C 2 ). This dose can be found by solving the above.
The boundary then reduces to a Hill equation of the form After logging in, click the "Create dataset" button. Enter a name for the dataset, and use the "browse" 273 button to select a CSV file from your computer. You can also select the orientation of the dose response 274 surface (whether Emax >E0, or otherwise), the name of the effect metric (e.g., percent effect, DIP rate), 275 apply any constraints to the E0 and Emax values (either unconstrained, fixed values, or range bounds).

276
Click the "Create dataset" button on this page, and the upload will begin. 277

278
The dataset is first uploaded to the server and validated. Each combination experiment is then sent to 279 a queue for processing. The processing step runs the MuSyC algorithm to fit the dose-response surface, and 280 returns the relevant fitting parameters. For small datasets, this process typically only takes a few minutes, 281 but this will vary depending on dataset size and server demand.

282
When the upload is complete, the web browser will redirect to a page showing the dataset name and a browser. For tasks marked as "FAILED", clicking on that word will show more details about the error (e.g., 290 if there was a data validation issue that the user should correct). There is also a link at the bottom of the 291 task result page to download that single combination's parameters as a CSV file. We created an interactive Jupyter Notebook that can be used to explore how MuSyC's 2D Hill-equation 294 parameters relate to the dose-response surface, as well as how they relate to Bliss and Loewe synergy.

295
This notebook, the backend code, and a conda environment file are included in File S1. The notebook 296 may be run locally, but we also provide a version online at https://mybinder.org/v2/gh/djwooten/ 297 natcomms-musyc2021/HEAD?filepath=demo.ipynb. To run the code locally, extract the contents of File 298 S1 to a directory. We recommend using Anaconda to manage installation of Python dependencies. On

299
MacOS or Linux, from a command line terminal in the directory of the contents of File S1, the following 300 commands will install the necessary Python dependencies and run the notebook: From Windows, we recommend using the Anaconda GUI application to set up the conda environment.

305
As noted in Table 2, Bliss and Loewe are not concentration independent synergy models, and thus synergy 306 can be calculated individually at each concentration. Doing so leads to a Bliss or Loewe "synergy landscape".

307
Using this interactive notebook, we explored how distinct types of synergy (e.g., potency, efficacy) may be 308 reflected in the Bliss or Loewe synergy landscapes. Figures S12-S16  a key advantage of MuSyC is that it is able to describe most real-world complex synergy landscapes in a 318 small number of easy to interpret synergy parameters.