Prediction of multidimensional drug dose responses based on measurements of drug pairs
See allHide authors and affiliations
Edited by Ruslan Medzhitov, Yale University School of Medicine, New Haven, CT, and approved July 19, 2016 (received for review April 22, 2016)

Significance
We present a mechanism-free formula that predicts effects of multiple drugs at all doses based on measurements of drug pairs at a few doses. The formula bypasses the combinatorial explosion problem by greatly reducing the number of measurements needed to design optimal cocktails for cancer and infection.
Abstract
Finding potent multidrug combinations against cancer and infections is a pressing therapeutic challenge; however, screening all combinations is difficult because the number of experiments grows exponentially with the number of drugs and doses. To address this, we present a mathematical model that predicts the effects of three or more antibiotics or anticancer drugs at all doses based only on measurements of drug pairs at a few doses, without need for mechanistic information. The model provides accurate predictions on available data for antibiotic combinations, and on experiments presented here on the response matrix of three cancer drugs at eight doses per drug. This approach offers a way to search for effective multidrug combinations using a small number of experiments.
To kill cancer cells or bacteria, combination therapy can be more effective than individual drugs (1⇓⇓⇓⇓–6). Combination therapy is thought to allow increased efficacy at low doses, thus reducing side effects and toxicity; it is also believed to minimize the chances of resistance (7⇓–9), a pressing problem in treating cancer and infectious diseases.
Much work has been devoted to classifying how pairs of drugs interact (10⇓⇓⇓–14). Across systems, a good first approximation is the Bliss independence model (15, 16), in which the pair effect is the product of the individual drug effects: If the effect of the drugs are g1 and g2, the effect of the combination is g12 = g1 · g2. The Bliss model ignores interactions in which drugs enhance each other effects—synergism—or inhibit each other’s effects—antagonism. Some drugs even inhibit each other so much that the combined effect is lower than either drug alone, an effect called hyperantagonism (17⇓⇓⇓–21).
Going beyond drug pairs has been difficult. Experimentally testing high-order combinations beyond pairs in a systematic way is challenging because it requires an exponentially large number of experiments (22⇓⇓–25): For N drugs at D doses, one needs DN experiments. For N = 10 drugs and D = 8 doses, this means ∼109 measurements. The combinatorial explosion makes exhaustive testing of drug dose combinations unfeasible. This problem is especially acute in cases where material is scarce, such as testing of patient-derived samples (26⇓⇓–29). Hence, models for predicting high-order effects are essential.
Apart from detailed simulations of particular systems (22, 30, 31), there has been little study of general mechanism-independent models for multiple drugs. An exception is the elegant study by Wood et al. (32) that showed that combinations of antibiotics can be predicted by an Iserliss-like formula that uses pair effects to predict the effects of higher-order mixtures. For example, the effect of a drug triplet is modeled as
Another line of research uses machine learning algorithms to make predictions based on a small number of measurements (23, 33⇓⇓⇓–37). Such approaches use regression to extrapolate from data, and often require iteration cycles of prediction and experiment to home in on solutions. Such machine learning approaches can benefit from a predictive formula that can make the search less prone to experimental error and can help avoid local maxima.
Here we study how multiple cancer drugs interact to affect cells by measuring the response matrix of three cancer drugs at eight doses per drug. We find that existing approaches do not accurately predict the response matrix. We therefore introduce a relatively simple model for drug mixtures. This model predicts high-order combination effects at all doses, based on measurements on single drugs and drug pairs at a small number of doses. The model provides accurate predictions for the three anticancer drug dose–response matrix measured in this study, and available data on antibiotic effects on bacterial growth rate (32), even in cases where synergy and antagonism is strong. This model is especially noise-resistant because it interpolates between dose measurements. The present model may be used to make accurate predictions of multiple drugs, based on only a few measurements of single and pair effects, thus bypassing the need for an exponentially large number of measurements.
Results
Current Approaches Do Not Describe the Response Matrix of Three Cancer Drugs.
We experimentally measured the dose–response matrix for eight doses of three commonly used chemotherapy drugs—doxorubicin, taxol, and cisplatin—on survival of A549 human lung cancer cells. Survival was measured using a cell viability assay (Methods).
The single-drug dose curves are given by the response matrix when the dose of the two other drugs is zero. The dose–response curves were well fit by Hill curves (10, 38, 39), characterized by a steepness parameter (Hill coefficient) n and a halfway point
Response matrix of three chemotherapy drugs at eight doses shows nonmonotonic behavior that is well captured by the dose model but not by other models. (A) Survival of A549 lung cancer cells when treated with the three pairs of drugs (taxol−doxorubicin, cisplatin−doxorubicin, cisplatin−taxol) at eight doses each: (Left) the measured response, (Middle) the predicted response using the present dose model, and (Right) the predicted response by Bliss independence. (B) Slices of the three-drug dose–response matrix. The first column is the measured response, followed by the prediction of the three drugs’ interactions using different models. Note that the Isserlis and regression models apply only to triplets and above, not to pairs.
To evaluate synergy and antagonism, we computed the Bliss approximation, given by the product of the single-drug response curves
Response of three cancer drug mixtures at multiple doses is best captured by the present dose model. Experimental measurements of the mixture of the three anticancer drugs taxol, doxorubicin, and cisplatin at all dose combinations compared with (A) Bliss independence, (B) Isserlis predictions, (C) regression model, and (D) present dose model predictions. R2 = 1 − ∑(Model − Experiment)/var(Experiment) can yield negative values when the mean model prediction is different from the mean experimental measurement, indicating a very poor fit.
We next compared the data to the Isserlis formula
We also performed regression similar to that used in machine learning approaches (36, 37, 43), which results in the formula
A Simple Dose Model for Drug Combinations.
We present a simple model for drug combinations. This is a model for three or more drugs, based on measurements of the response of single drugs and drugs pairs at a few dose combinations.
The model extends the Bliss formula; it is based on the product of the effects of all drugs in the mixture, not at their true doses but rather at effective doses that differ from the true doses due to interactions with the other drugs in the mixture. This interaction is modeled by introducing interaction terms between drug pairs (14). In the model, the effective dose of each drug, Dieff, is its true dose Di times a product of the Michaelis−Menten-like terms for all of the other drugs
To test this model, we begin by considering two-drug response surfaces (3, 45, 46). Thus, we consider Eq. 2 with only two drugs. The response surface in the model is a product of Hill functions gi of the effective doses D1eff and D2eff
The model captures all three pairs of cancer drugs studied here with high accuracy (Fig. 1A). It shows R2 = 0.87, 0.91, and 0.86 for the pairs doxorubicin−taxol, doxorubicin−cisplatin, and taxol−cisplatin, respectively. For comparison, Bliss independence gives R2 = 0.16, 0.83, and 0.7. All three pairs of drugs show antagonism (aij > 0), with the antagonism of taxol and doxorubicin being most significant.
For drug pairs, there exists an additional wealth of published data to test the model. The model fits well all interactions of 20 antibiotic pairs and two pairs of anticancer drugs measured in ref. 14 (Fig. 3 and Fig. S1A). The model captures almost all of the variation in the data (<R2> = 0.95). For comparison, the Bliss independence model gives <R2> = 0.58. The model is even capable of capturing hyperantagonism (47), where the pair has a lower effect than either of the two single drugs [such as the pair trimethoprim−lincomycin (Fig. 3B)].
The present dose model accurately describes response surfaces of antibiotic and chemotherapy drug pairs. Examples of response surfaces for antibiotic pairs that are (A) antagonistic (B) hyperantagonistic, and (C) synergistic. (D) Interactions of anticancer drug pairs also well described by the model. (E) R2 values of the 20 drug pairs considered here for Bliss independence (o), and the present dose model of Eqs. 3 and 4 (x). Experimental data from Wood et al. (14). Cip, ciprofloxacin; Cm, chloramphenicol; Dox, doxycycline; Ery, erythromycin; Linc, lincomycin; Ofl, ofloxacin; Sal, salicylate; Tet, tetracycline; Tmp, trimethoprim (see also Fig. S1A).
The errors of the dose model seem to result mainly from experimental noise. (A) RMSE of the dose model for 20 drug pairs (o) as a function of the maximal RMSE of the Hill model for the single drugs. The correlation between the two shows that much of the error in the dose model errors can be attributed to noise. Similar results for three anticancer drugs are also presented (□). (B) RMSE of the dose model for three- and four-drug combinations as a function of the maximal RMSE of the dose model for the pairs of drugs in the combination (o). The correlation between the two suggests that much of the dose model error can be attributed to noise. Similar results for three anticancer drugs are also presented (□).
Importantly, the present dose model requires only a few measurements to accurately fit the entire dose surface of a drug pair, because it has only two (or one) free parameters. In all cases considered here, 10 measurements on pairs of doses suffice to provide an excellent fit to the entire experimental dose surface of a drug pair (Fig. S2).
Accuracy analysis suggests that a single interaction parameter suffices for the pair model in most cases, and that only about 10 measurements suffice to specify the pair parameters. (A) The loss of accuracy (ΔRMSE) when using only one parameter (a12 or a21) to describe pair interactions. The best results (lower RMSE) out of the two options are presented. In most cases, the loss of accuracy is small (<0.5%). In the case of superantagonism (Tmp-Ofl), the loss of accuracy is large, because both interaction terms are needed to describe the effect. (B) The loss of accuracy as a function of the number of pairs of doses measured to estimate the parameters of a given pair of drugs. Accuracy is relative to the measurement of all dose pairs in the dataset (∼100). For 10 measurements, the loss of accuracy is about 0.5%.
The model can be made even simpler in many cases: For many drug pairs, setting one of the two interaction parameters a12 or a21 to zero obtains a model with a single fitting parameter. This model reduction gives a loss of accuracy of less than 0.5% for 80% of the drug pairs (Fig. S2).
The finding that a single interaction parameter captures many of the pairs provides a picture of a hierarchy between the drugs: Drug 1 changes the effective dose of drug 2 but not vice versa (1→2). A hierarchy is obtained because we find that that these relations are transitive: If drug 1 affects drug 2 (1→2) and drug 2 affects drug 3 (2→3), then drug 1 will affect drug 3 (1→3) but not vice versa. The hierarchy between antibiotics is shown in Fig. 4. For example, salicylate is high on the hierarchy because it activates a multidrug efflux system that transports other drugs out of the bacterial cell (48) (Fig. S3).
A hierarchy between antibiotics. For most of the drug pairs, one of the two interaction parameters a12 or a21 can be set to zero without significant loss of fit quality, resulting in a single parameter for the interaction. In these cases, we can say that drug 1 changes the effective dose of drug 2 but not vice versa (1→2). We find that these relations are transitive: If drug 1 affects drug 2 (1→2) and drug 2 affects drug 3 (2→3), then drug 1 will affect drug 3 (1→3) but not vice versa. The direction of the interaction is transitive, but the absolute strength of the interaction is not always transitive. To quantify the strength of the hierarchy we use the goodness of fit (RMSE) difference in the case of a12 = 0 and a21 = 0. This quantity is indicated by the arrow thickness. The strength of the interaction is indicated by the color of the arrow (red, strong antagonistic; blue, strong synergistic). We plotted only interaction with nonnegligible hierarchy (∆RMSE > 0.005%). Note that 8 out of the 28 possible interactions were not present in the published dataset.
A mechanistic motivation for the two-drug interaction term based on drug efflux pump induction. We use the well-understood drug interaction mechanism in the mar system. In this system, salicylate and other drugs interact through the induction of the mar efflux system by salicylate. The efflux system exports the other drug and effectively reduces its dose (14). The figure shows data from ref. 14 on the promoter activity of the mar system as a function of salicylate dose. We find that this induction curve is fit very well by the present Michaelis–Menten-like interaction term between the two drugs salicylate and tetracycline. Additional mechanisms are probably at play in other cases.
The Dose Model Describes Combinations of Three and Four Antibiotics Using Only Pair Data.
We now turn to the anticancer drug triplet measured here. We used the aij interaction parameters for the three pairs in the triplet. We find that the dose model for three drugs (Eq. 2) describes the full three-drug interaction matrix very well (R2 = 0.82; Figs. 1 and 2).
We further compared the model for triplets and quadruplets of antibiotics published previously (32). The dataset includes six triplets and two quadruplets at a total of 1,384 dose combinations. The drug pairs that make up these combinations have synergy and antagonism, and several of the triplet and quadruplets are markedly non-bliss-like, so that prediction is challenging. We find that the model of Eq. 2, with aij parameters fitted from drug pair data, works very well (R2 ∈ [0.85, 0.93]) for all of these combinations (Fig. 5 and Fig. S1B).
Combinations of three or four drugs are well described by the present dose model (x). Bliss independence (o) results are displayed for comparison. The three-drug combination Cm-Ofl-Sal shows strong antagonism (B) The three-drug combination Cm-Ofl-Tmp show synergism. (C) The three-drug combination Cm-Ery-Tmp shows synergism. (D) The three-drug combination Sal-Ery-Cm shows mild antagonism (E) The four-drug combination Linc-Cm-Ofl-Tmp shows complex interactions captured by the model. (F) The four-drug combination Dox-Ery-Linc-Sal shows antagonism. (G) R2 values for three- and four-drug combinations: Bliss model (o) and the present dose model (x). Antibiotic data are from Wood et al. (32); results for the present experiments on three anticancer drugs on A549 cells (Figs. 1 and 2) are also presented. Cm/C, chloramphenicol; Dox/D, doxycycline; Ery/E, erythromycin; Linc/L, lincomycin; Ofl/O, ofloxacin; Sal/S, salicylate; Tmp/T, trimethoprim. Isserlis model for quadruplets is g1234 = g12g34 + g13g24 + g14g23 − 2g1g2g3g4 (see also Fig. S1B).
For example, consider the antibacterial triplet chloramphenicol, ofloxacin, and salicylate. Salicylate is highest in the hierarchy and antagonistic to the other two antibiotics (48, 49). In the model, this antagonism is represented by the parameters a23 = 7.2 ± 0.4 and a13 = 2.8 ± 0.3. The remaining pair has mild antagonism (a21 = 0.5 ± 0.02). The triplet shows very strong antagonism predicted accurately by the model [R2 = 0.88; compare with Bliss R2 ≈ 0 (Fig. 5A)].
A second example that includes both synergy and antagonism is the triplet chloramphenicol, ofloxacin, and trimethoprim (Fig. 5B). Here ofloxacin and trimethoprim are synergistic (a23 = −0.8 ± 0.05) as are chloramphenicol and trimethoprim (a31 = −0.57 ± 0.08), whereas chloramphenicol and ofloxacin are antagonistic (a21 = 0.99 ± 0.07). Thus, a priori, it is difficult to predict the synergy/antagonism of the triplet. The model predicts overall synergy, and fits the experiment very well (R2 = 0.93 compared with Bliss R2 = 0.86).
The Dose Model May Be Used to Optimize Drug Combinations for Minimal Side Effects.
The dose model described here can be used to navigate the space of multidrug doses to find combinations of interest. As a schematic example, we consider the problem of finding a drug combination with a given efficacy that has minimal side effects. This example is meant to demonstrate a general approach, which can be made more realistic with accurate models for side effects or other features that need to be optimized.
To illustrate this, we consider a triplet of antibiotics—chloramphenicol, erythromycin, and trimethoprim. The model allows calculating all dose combinations with a given effect, say a 10-fold growth reduction of the bacteria, g = 0.1. The dose combinations that give g = 0.1 form a surface in the 3D space of doses (Fig. 6).
Search for a mixture with high efficacy at low doses. We seek a combination of three antibiotics with a given large effect (g = 0.1). The doses that satisfy this condition form a surface in the 3D dose space for chloramphenicol (Cm), erythromycin (Ery), and trimethoprim (Tmp). On this surface, we seek the combination with minimal side effects. Assuming that side effects increase with drug doses, we can seek to minimize, for example, the maximal dose in the combination M = max(Di/D0i). The colors indicate the decrease in M relative to the lowest Di/D0i of the individual drugs. The optimal drug combination indicated by an arrow on the g = 0.1 surface has a fourfold-lower M than the best individual drug (see also Fig. S4).
On this surface, we seek the combination with minimal side effects. Assuming that side effects increase with drug doses, we can seek to minimize, for example, the sum of normalized doses S =
Search for high efficacy at low dose combinations. We seek a combination with a given large effect (g = 0.1). The doses that satisfy this condition form a surface in the 3D dose space, for example, the g = 0.1 surface for the three antibiotics: Cm, Ery, and Tmp. On this surface, we seek the combination with minimal side effects. Assuming that side effects increase with drug doses, we can seek to minimize, for example, the sum of normalized doses S =
Discussion
We presented a model for the effects of multidrug mixtures for anticancer drugs and antibiotics. The model does not require knowledge of the mechanism of drug action, and is based on relatively few measurements on drug pairs. The model accurately describes published data on antibiotic triplets and quadruplets, as well as experimental data presented here on dose–response of combinations of three cancer drugs.
The model has two interaction parameters for each pair, a12 and a21. These parameters describe how the effective dose of a drug is affected by the presence of the other drug. The parameters can be accurately estimated based on about 10 dose combinations for each pair. The use of several dose combination measurements to estimate the parameters helps reduce the effects of experimental noise, partly explaining why the present model is more accurate than previous mechanism-free models such as the Isserlis model and regression-based models (Fig. S5).
The Isserlis model agrees with the present model in the limit of weak interactions and low noise. (A) We generated synthetic data using the present model and fit the data using the Isserlis model. We used doses D/D0 = 0 to Dmax, where g(Dmax) = 0.2, Hill coefficients for the two simulated drugs were in the range n = 1–7, and interaction parameters aij in the range −1 to 6 (note that the smallest possible value for aij is −1). The models agree well when the interactions are weak (n < 2, |aij| < 3) (R2 > 0.9). Agreement breaks down when the interactions are stronger or Hill coefficients for at least one drug exceed 2. (B) The agreement breaks down even for weak interactions in the presence of noise. Here sigma is the SD of Gaussian noise added to the simulated log effect data. Hill coefficient for the two drugs in this example was n = 2. (C and D) We also tested a smoothed version of the Isserlis model on the present cancer drug dose–response matrix, in which the single and pair drug inputs to the Isserlis model are taken from the present model. Because the present model uses multiple dose measurements, this input has reduced noise. This smoothed Isserlis model fits the present cancer drug data better than the unsmoothed Isserlis model [R2 = 0.7 (C) vs. R2 = 0.6 (D)] but not as well as the present model (R2 = 0.82) (Fig. 2D). This indicates that part of the improvement in R2 is due to noise insensitivity (∼0.1) and part of it is systematic (∼0.1).
The fitted interaction parameters reveal a hierarchy between the drugs. Some drugs impact the effective dose of the other drug but not vice versa. This relationship is transitive and can offer clues about the mechanisms of synergy and antagonism. For example, a drug high in the hierarchy is likely to set off general drug resistance systems, as in the case of salicylate, which activates the mar system in bacteria (Fig. S3).
The present approach can overcome the combinatorial explosion problem in which the number of experiments rises exponentially with the number of drugs and doses, because it requires few measurements. Instead of DN measurements for N drugs and D doses, we can use the single-drug dose–response curves plus a sampling of the pair dose combinations to fit the interaction parameters aij. In total, this model requires only a quadratic number of measurements, N · D + cN(N − 1)/2, where c is ∼10. This means that for N = 6 drugs and D = 8 doses, we need only 198 measurements instead of DN ≈ 3 · 105. For N = 10 drugs and D = 8 doses, we need only 530 measurements instead of ∼109. This approach therefore potentially allows navigating the space of high-order drug combinations to search for mixtures with high efficacy at low dose.
Methods
Regression Model.
Machine learning approaches often use regression models. A commonly used regression formula employs the variables xi = 0 or 1 that denote the absence or presence of drug i in the mixture. The effect of the mixture g is described in the model as
Estimation of Model Parameters.
We determine the values and confidence intervals of the model parameters (ni, D0i, and aij) using the MATLAB function “fit.” We use the single-drug dose–response curve to estimate the values and uncertainty of ni and D0i for each drug. The interaction parameters aij are determined by the two-drug response matrix. In the process of fitting aij, we allow ni and D0i to change within their uncertainty estimated from the single-drug measurements.
Computation of Effective Doses in the Model.
In the general case aij ≠ 0, we determine the effective dose by numerically solving Eq. 2 using the MATLAB function “fmincon.” When there exists a hierarchy between the drugs in the combination and {aij = 0|aji = 0}, there is an analytic solution to Eq. 2 that one can use instead of a numerical solution
Cell Lines and Culture Conditions.
Measurements were performed by BioDuro (email bioduroinfo{at}bioduro.com). A549 (non-small cell lung cancer) cells were grown in 384-well plates (#3707; Corning) plated at 104 cells per well. Plates were incubated at 37°C and 5% (vol/vol) CO2 for 24 h in membrane-filtered growth media composed of RPMI (Roswell Park Memorial Institute) 1640 with l-glutamine (#11875093; Invitrogen) supplemented with 10% FCS [certified FBS, FBS (#SH30084.03; HyClone)], and 0.05% Penicillin-Streptomycin (#15140-122; Invitrogen). After 24 h, the drugs (taxol, doxorubicin, and cisplatin) were added (2.5 μL per well) to final concentrations of 20 μM, 10 μM, 3.3 μM, 1.1 μM, 370 nM, 123 nM, 41 nM, and 13.7 nM and two controls: 0.1% DMSO and 2.5 μL of regular media (untreated). Cells were incubated with drugs for 48 h, followed by a cell survival assay [3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide].
Cell Survival Assay.
The viability assay was done using an MTT-based in vitro toxicology assay (#C0009; Beyotime). Cells were incubated with reconstituted MTT for 4 h, and absorbance was measured at 570 nm.
Acknowledgments
U.A. is the incumbent of the Abisch-Frenkel Professorial Chair. Support was provided by the Clore-Katz-David Foundation and Fela Shapell Family Foundation, and the Israel National Center for Personalized Medicine (INCPM) Fund for Preclinical Studies.
Footnotes
↵1A.Z. and I.K. contributed equally to this work.
- ↵2To whom correspondence should be addressed. Email: uri.alon{at}weizmann.ac.il.
Author contributions: U.A. designed research; A.Z., I.K., E.D., and A.E.M. performed research; A.Z. performed experiments; I.K. contributed new reagents/analytic tools; I.K. designed the model and analyzed data; and A.Z., I.K., and U.A. wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission.
See Commentary on page 10231.
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1606301113/-/DCSupplemental.
Freely available online through the PNAS open access option.
References
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Held MA, et al.
- ↵
- ↵.
- Stylianou M, et al.
- ↵
- ↵
- ↵.
- Bliss CI
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Blom K,
- Nygren P,
- Alvarsson J,
- Larsson R,
- Andersson CR
- ↵.
- Pemovska T, et al.
- ↵.
- Sequist LV, et al.
- ↵
- ↵.
- Janes KA, et al.
- ↵.
- Wood K,
- Nishida S,
- Sontag ED,
- Cluzel P
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Chou T-C
- ↵.
- Chou T-C
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Lehár J,
- Krueger A,
- Zimmermann G,
- Borisy A
- ↵
- ↵
- ↵.
- Cohen SP,
- Levy SB,
- Foulds J,
- Rosner JL
- ↵
Citation Manager Formats
Article Classifications
- Biological Sciences
- Systems Biology
See related content: