# 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 *g*_{1} and *g*_{2}, the effect of the combination is *g*_{12} = *g*_{1} · *g*_{2}. 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 *D*^{N} experiments. For *N* = 10 drugs and *D* = 8 doses, this means ∼10^{9} 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

To evaluate synergy and antagonism, we computed the Bliss approximation, given by the product of the single-drug response curves *A*). The drug combinations show reduced effect compared with Bliss, which means that they show antagonism (Fig. 2*A*). Antagonism is partly due to cell cycle effects where taxol tends to prevent entry of cells into the S phase in which the other two drugs are most effective (40⇓–42).

We next compared the data to the Isserlis formula *R*^{2} = 0.6 (Fig. 2*B*). One reason for this relative lack of fit is experimental noise: The Isserlis model depends sensitively on errors in the three input pairs *g*_{12}, *g*_{13}, and *g*_{23} and the three single-effect data points *g*_{1}, *g*_{2}, and *g*_{3}.

We also performed regression similar to that used in machine learning approaches (36, 37, 43), which results in the formula *Methods*). This regression formula shows poor performance (Fig. 2*C*). We conclude that an improved formula is needed to describe the present three-drug response matrix.

### 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, *D*_{ieff}, is its true dose *D*_{i} times a product of the Michaelis−Menten-like terms for all of the other drugs*g*_{i} are Hill functions appropriate to each single drug (Eq. **1**) with steepness *n _{i}* and halfway point

*a*

_{ij}interaction parameters are computed from two-drug data (

*Methods*). This model assumes that the effects of drugs on each other’s effective doses are multiplicative, and that one can neglect third-order and higher interactions (44).

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 *g*_{i} of the effective doses *D*_{1eff} and *D*_{2eff}*a*_{12} and *a*_{21}, which are fit from the pair data (*Methods*).

The model captures all three pairs of cancer drugs studied here with high accuracy (Fig. 1*A*). It shows *R*^{2} = 0.87, 0.91, and 0.86 for the pairs doxorubicin−taxol, doxorubicin−cisplatin, and taxol−cisplatin, respectively. For comparison, Bliss independence gives *R*^{2} = 0.16, 0.83, and 0.7. All three pairs of drugs show antagonism (*a*_{ij} > 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. S1*A*). The model captures almost all of the variation in the data (<*R*^{2}> = 0.95). For comparison, the Bliss independence model gives <*R*^{2}> = 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. 3*B*)].

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).

The model can be made even simpler in many cases: For many drug pairs, setting one of the two interaction parameters *a*_{12} or *a*_{21} 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).

### 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 *a*_{ij} 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 (*R*^{2} = 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 *a*_{ij} parameters fitted from drug pair data, works very well (R^{2} ∈ [0.85, 0.93]) for all of these combinations (Fig. 5 and Fig. S1*B*).

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 *a*_{23} = 7.2 ± 0.4 and *a*_{13} = 2.8 ± 0.3. The remaining pair has mild antagonism (*a*_{21} = 0.5 ± 0.02). The triplet shows very strong antagonism predicted accurately by the model [*R*^{2} = 0.88; compare with Bliss *R*^{2} ≈ 0 (Fig. 5*A*)].

A second example that includes both synergy and antagonism is the triplet chloramphenicol, ofloxacin, and trimethoprim (Fig. 5*B*). Here ofloxacin and trimethoprim are synergistic (*a*_{23} = −0.8 ± 0.05) as are chloramphenicol and trimethoprim (*a*_{31} = −0.57 ± 0.08), whereas chloramphenicol and ofloxacin are antagonistic (*a*_{21} = 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 (*R*^{2} = 0.93 compared with Bliss *R*^{2} = 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).

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* = *D*_{i}/*D*_{0i} (Fig. S4). Another possibility is to minimize the maximal dose in the combination *M* = max(*D*_{i}/*D*_{0i}). In both cases, an optimal drug combination is found that uses lower doses than the best individual drug (indicated by an arrow on the *g* = 0.1 surface in Fig. 6). The combination that minimizes summed doses has a twofold lower total dose than the best individual drug (Fig. S4). The combination that minimizes maximal dose has a fourfold lower dose than the best individual drug (Fig. 6).

## 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, *a*_{12} and *a*_{21}. 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 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 *D*^{N} 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 *a*_{ij}. 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 *D*^{N} ≈ 3 · 10^{5}. For *N* = 10 drugs and *D* = 8 doses, we need only 530 measurements instead of ∼10^{9}. 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 *x*_{i} = 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 *b*_{i} and *c*_{ij} are fit from single and pair data. For triplets, the equivalent formula is

### Estimation of Model Parameters.

We determine the values and confidence intervals of the model parameters (*n*_{i}, *D*_{0i}, and *a*_{ij}) using the MATLAB function “fit.” We use the single-drug dose–response curve to estimate the values and uncertainty of *n*_{i} and *D*_{0i} for each drug. The interaction parameters *a*_{ij} are determined by the two-drug response matrix. In the process of fitting *a*_{ij}, we allow *n*_{i} and *D*_{0i} to change within their uncertainty estimated from the single-drug measurements.

### Computation of Effective Doses in the Model.

In the general case *a*_{ij} ≠ 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 {*a*_{ij} = 0|*a*_{ji} = 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 10^{4} cells per well. Plates were incubated at 37°C and 5% (vol/vol) CO_{2} 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

↵

^{1}A.Z. and I.K. contributed equally to this work.- ↵
^{2}To 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

*Escherichia coli*: Activation of the mar operon and a mar-independent pathway. J Bacteriol 175(24):7856–7862. - ↵

## Citation Manager Formats

## Article Classifications

- Biological Sciences
- Systems Biology

## See related content: