A New Bliss Independence Model to Analyze Drug Combination Data

Abstract

The Bliss independence model is widely used to analyze drug combination data when screening for candidate drug combinations. The method compares the observed combination response (Y_O) with the predicted combination response (Y_P), which was obtained based on the assumption that there is no effect from drug-drug interactions. Typically, the combination effect is declared synergistic if Y_O is greater than Y_P. However, this method lacks statistical rigor because it does not take into account the variability of the response measures and can frequently cause false-positive claims. In this article, we introduce a two-stage response surface model to describe the drug interaction across all dose combinations tested. This new method enables robust statistical testing for synergism at any dose combination, thus reducing the risk of false positives. The use of the method is illustrated through an application describing statistically significant “synergy regions” for candidate drug combinations targeting epidermal growth factor receptor and the insulin-like growth factor 1 receptor.

Introduction

Various oncogenic cell signaling pathways are known to provide cross-talk and redundancy within tumors. ¹ Thus, inhibition of such pathways individually by a single targeted therapy has been shown to lead to compensation by other pathways. This, in turn, results in a loss of sensitivity to the original targeted therapeutic agent at the cellular level. In the clinic, this type of compensation leads to innate and/or acquired tumor resistance and relapse.

Because advanced tumors are often resistant to single agents, there is an increasing trend to combine drugs to achieve better treatment effect and reduce safety issues.^2,3 It is desirable that the combination drugs are synergistic; that is, better activity is achieved at lower dose levels when drugs are combined than individually observed at single drug doses. Depending on drug doses, drug combinations can yield activity that is synergistic, independent, or antagonistic.

Commonly used statistical models to evaluate drug combination efficacy are the Bliss independence and Loewe additivity models. Greco et al. ⁴ have discussed these two reference models in detail. Conceptually, the Loewe additivity model focuses on dose reduction and the Bliss independence model focuses on treatment effect enhancement. The two reference models handle the same question from two different perspectives. The Bliss independence model has been criticized for its potential to incorrectly claim synergy when two identical drugs are combined (i.e., the “sham” experiment). But in reality, ethical and economic constraints preclude “sham” testing of a combination of two identical drugs to evaluate synergy in a clinical setting. Furthermore, we find that Bliss independence is legitimately derived from the complete additivity of probability theory ⁵ and serves well as a reasonable reference model.

Suppose two drugs, A and B, both inhibit tumor growth: drug A at dose a inhibits Y_a percent of tumor growth and drug B at dose b inhibits Y_b percent of tumor growth. If two drugs work independently, the combined percentage inhibition Y_ab,P can be predicted using the complete additivity of probability theory as

$$Y_{ab,P}=Y_a+Y_b-Y_a{Y}_b.$$

The observed combined percentage inhibition Y_ab,O is then compared with Y_ab,P. Typically, if Y_ab,O > Y_ab,P, the combination treatment is thought to be more efficacious than expected; if Y_ab,O < Y_ab,P, the combination treatment is worse than expected; and if Y_ab,O =Y_ab,P, the combination is equal to a simple addition of two separate drugs. More concisely, the three scenarios are summarized as

$$Y_{ab,O}\{\begin{array}{cc}\begin{array}{l}>Y_{ab,P}\\ =Y_{ab,P}\end{array}& \begin{array}{l}Synergy\\ Independent\end{array}\\ <Y_{ab,P}& Antagonism\end{array}.$$

Typically, the comparison is examined at each dose combination. Average percentage inhibition is usually used when replicates are available and a conclusion can be drawn by strictly following equations (1) and (2). But performing tests repetitively at all combination doses becomes a multiple testing problem, which easily leads to false-positive claims.

Response surface methodology (RSM) is a well-established method to explore relationships between dose combinations and treatment effect in combination studies. Recently, a Loewe additivity–based two-stage response surface model for a drug interaction index was developed. ⁶ The method gives a global view of synergistic effects across all drug combinations. However, we find in practice that the Loewe additivity–based models work well only when monotherapy dose-response curves are well characterized through parametric models such as four-parameter logistic (4PL) nonlinear regression curves. When dose-response curves are nonstandard or not available, we find that the Bliss independence–based model provides a viable alternative. This is particularly important in high-throughput screening settings where large dose ranges are less practical. Even when large dose ranges are feasible for combination screening, nonstandard response curves compromise the reliability of the Loewe additivity model. In the following section, we adopt the two-stage framework to develop a response surface model based on Bliss independence.

Methods

We use a response surface model $M(\Gamma |a,b)$ to describe the difference between the observed percentage inhibition and the predicted percentage inhibition as

$$I=Y_{ab,O}-Y_{ab,P}=M(\Gamma |a,b)+\epsilon ,$$

where ε is a random error that is normally distributed with mean 0 and variance ${\sigma }^2$ , Г={γ₀, γ₁, γ₂, γ₃, γ₄, γ₅} is a set of response surface model parameters. I is called the interaction index, an indicator for synergism. In this article, the response surface model we use is a quadratic polynomial function,

$$M(\Gamma |a,b)={\gamma }_0+{\gamma }_{1}a+{\gamma }_{2}b+{\gamma }_{3}ab+{\gamma }_4{a}^2+{\gamma }_5{b}^2.$$

If I is known, Г can be easily solved using the multiple linear regression method. Moving $Y_{ab,P}$ to the right side of the equation (3), we have

$$Y_{ab,O}=Y_{ab,P}+M(\Gamma |a,b)+\epsilon .$$

Since $Y_{ab,P}$ is a random term itself, its distribution will be estimated in the first stage and then be incorporated in the equation (3) to estimate parameters of the response surface model. In this way, the variance of the response surface model parameters will include the variance of $Y_{ab,O}$ and the variance of $Y_{ab,P}$ .

First Stage Model of the Predicted Percentage Inhibition

We estimate the predicted percentage inhibition, $Y_{ab,P}$ , in the first stage. The monotherapy percentage inhibition models at dose (a, b) can be written separately using simple linear regression models as

$$\begin{array}{l}{Y}_a={\mu }_a+e_a,\\ Y_b={\mu }_b+e_b,\end{array}$$

where µ _i (i = a, b) is the true percentage inhibition for drug A and B at dose (a, b) and e_i is a normally distributed random error with mean 0 and variance ${\sigma }_i^2$ . Because negative inhibition values are often observed in real experiments, we model the percentage inhibition purposely using a Gaussian model instead of the more traditional logistic regression function to accommodate possible negative values. The similar modeling approach was adopted by other researchers as well. ⁷ The predicted percentage inhibition of the combination treatment can then be written as

$$Y_{ab,P}={\mu }_a+{\mu }_b-{\mu }_a{\mu }_b+e_a+e_b-{\mu }_a{e}_b-{\mu }_b{e}_a-e_a{e}_b.$$

Since the expectations of all the terms with random errors are zero, the expectation of Y_ab,P can be written as

$$E\left(Y_{ab,P}\right)={\mu }_a+{\mu }_b-{\mu }_a{\mu }_b,$$

and the variance of Y_ab,P and covariance between two combination doses, ab and $a\prime b\prime$ , can be written as

$$\mathrm{var}\left(Y_{ab,P}\right)={\left(1-{\mu }_b\right)}^2{\sigma }_a^2+{\left(1-{\mu }_a\right)}^2{\sigma }_b^2+{\sigma }_a^2{\sigma }_b^2,$$

$$\mathrm{cov}\left(Y_{ab,P},Y_{a^{\prime }{b}^{\prime },P}\right)=\{\begin{array}{cc}\begin{array}{l}\left(1-{\mu }_b\right)\left(1-{\mu }_{b^{\prime }}\right){\sigma }_a^2,\\ \left(1-{\mu }_a\right)\left(1-{\mu }_{a^{\prime }}\right){\sigma }_b^2,\end{array}& \begin{array}{l}\text{if}a=a^{\prime },b\ne b^{\prime }\\ \text{if}a\ne a^{\prime },b=b^{\prime }\end{array}\\ 0,& \text{if}a\ne a^{\prime },b\ne b^{\prime }\end{array}.$$

Second Stage Model to Calculate Parameter Variance

Conditional on the predicted percentage inhibition ${\widehat{Y}}_{ab,P}$ , the difference between $Y_{ab,O}$ and ${\widehat{Y}}_{ab,P}$ is modeled using a response surface model. All parameters with a ^ sign are estimated using statistical models.

Following the two-stage paradigm, ⁶ the estimated variance covariance matrix of model parameters $\widehat{\Gamma }$ can be written as

$$\mathrm{var}\left(\widehat{\Gamma }\right)=E\left(\mathrm{var}\left(\widehat{\Gamma }|{\widehat{Y}}_{ab,P}\right)\right)+\mathrm{var}\left(E\left(\widehat{\Gamma }|{\widehat{Y}}_{ab,P}\right)\right).$$

The first term on the right-hand side of equation (9) is the mean variance of $\widehat{\Gamma }$ assuming the predicted response is ${\widehat{Y}}_{ab,P}$ , and the second term is the variance of the expected $\widehat{\Gamma }$ . Properly calculating $\mathrm{var}\left(\widehat{\Gamma }\right)$ is of critical importance in constructing the confidence interval of the interaction index. The bootstrap method is commonly used for this type of calculation. First, a random sample of ${\widehat{Y}}_{ab,P}$ is simulated from a multivariate normal distribution with mean and covariance matrix parameters calculated using equations (6), (7), and (8). Conditional on ${\widehat{Y}}_{ab,P}$ , response surface parameters $\widehat{\Gamma }$ in equation (3) and $\mathrm{var}\left(\widehat{\Gamma }|{\widehat{Y}}_{ab,P}\right)$ in equation (9) are calculated using the linear model method. These steps are repeated for a number of times and then $\mathrm{var}\left(\widehat{\Gamma }\right)$ can be approximated by replacing expectation and variance with sample mean and sample variance in equation (9). $\mathrm{var}\left(\widehat{\Gamma }\right)$ roughly stabilizes after 50 iterations. Usually, the difference of $\mathrm{var}\left(\widehat{\Gamma }\right)$ between 50 and 10,000 iterations is within 15%, and the difference between 1000 and 10,000 repeats is within 5%. The 100(1 – α)% confidence interval for interaction index I at any combination dose (a, b) can then be calculated using the following formula:

$$CI=\left(\widehat{I}-z_{1-\alpha /2}\sqrt{x^{\prime }\mathrm{var}(\widehat{\Gamma })x},\widehat{I}+z_{1-\alpha /2}\sqrt{x^{\prime }\mathrm{var}(\widehat{\Gamma })x}\right),$$

where $x={\left(1,a,b,ab,a^2,b^2\right)}^{\prime }$ and $z_{\alpha /2}$ is the 100(α/2)th percentile of the standard normal distribution. α is the type I error rate, which is the probability to falsely claim the combination is synergistic when it is not. α = 0.05 is commonly used as a benchmark error rate.

Results

Data for combinations studies reported here arise from in vitro combination experiments testing antibodies against two well-characterized and clinically validated oncology targets: the epidermal growth factor receptor (EGFR) and the insulin-like growth factor 1 receptor (IGF1R). Each of these targets has been shown to play a role in tumor biology. ⁸ Both are frequently coexpressed in tumor cells ⁹ and have been shown to be a part of the same cellular complexes using coimmunoprecipitation experiments. ¹⁰ Furthermore, cross-talk between the IGF1R and EGFR receptor has been shown to result in compensation and resistance during targeted inhibition of either receptor. ¹¹ Finally, various preclinical studies have shown improved activity when both IGF1R and EGFR are inhibited concomitantly.^12,13

Given that anti-IGF1R and anti-EGFR antibodies have shown toxicity clinically, it is important to study the combined effects of anti-IGF1R and anti-EGFR treatment quantitatively. Of particular interest is the hypothesis that subtoxic doses of each agent might yield efficacy comparable to that seen with toxic doses separately because of synergism when tested in combination. These data are evaluated using the proposed model with a view toward predicting superior efficacy in combination than that observed when the same antibodies are dosed separately.

In this experiment, both antibodies were tested at four doses (0.5, 5, 50, and 500 nM) in quadruplicate. Including monotherapy data, there are 96 data points in total. Figure 1 clearly shows that drug A does not follow a classic 4PL dose-response curve, which precludes the use of Loewe additivity–based models. Table 1 gives the estimated parameters of the response surface model. The 95% lower bound of the estimated interaction index is shown in Figure 2 , and the contour plot with log₁₀-transformed doses is shown in Figure 3 . Strong synergistic interaction is seen at lower doses of both antibodies. Subsequent animal studies using the same antibodies in a xenograft model revealed synergistic activity arising from this target combination. ¹⁴

Figure 1.

Monotherapy dose-response curves for the anti–insulin-like growth factor 1 receptor (IGF1R) antibody and anti–epidermal growth factor receptor (EGFR) antibody for the non–small cell lung cancer (NSCLC) data.

Table 1.

Estimated Parameters for the Non–Small Cell Lung Cancer Data.

${\gamma }_0$	${\gamma }_1$	${\gamma }_2$	${\gamma }_3$	${\gamma }_4$	${\gamma }_5$
0.077	0.018	0.026	−0.001	−0.014	−0.012

Figure 2.

Estimated 95% lower confidence bound of interaction index based on Bliss Independence at all dose combinations. EGFR, epidermal growth factor receptor; IGF1R, insulin-like growth factor 1 receptor.

Figure 3.

Contour plot of interaction index for the anti–insulin-like growth factor 1 receptor (IGF1R) antibody and anti–epidermal growth factor receptor (EGFR) antibody combination for the non–small cell lung cancer (NSCLC) data.

Discussion

The Bliss independence–based model has fewer restrictions than Loewe additivity–based models. First, Loewe additivity models rely on accurately estimated dose-response curves to support the calculation of the effective dose for a given response. When a 4PL model is used for estimating a dose response, it is mandatory that the response has to fall between the estimated E_min and E_max parameters—a result that is often not possible. When the data point is not between E_min and E_max, it has to be manually removed from the analysis, which is undesirable for statistical analysis. Second, the Loewe additivity model becomes unusable when a dose-response curve is not available or difficult to model. Third, the Loewe additivity model is far more computationally challenging than the Bliss independence model. When both methods have solutions, their results are very similar (comparisons not shown). The Bliss independence model presented in this article requires only a linear model technique, which is readily available in most software.

Screening for effective combination drugs among hundreds of possible candidates is a challenging task. Those analytical methods without proper statistical modeling can easily lead to false decisions. The methodology introduced in this article rigorously considers the variances in both monotherapy and combination drug experiments and provides a practical way to integrate those variances into one comprehensive model. The confidence intervals constructed using this method have successfully helped us to identify synergistic regions of the test antibodies. Finally, it should be emphasized that most other methods provide only a simple score to quantify synergism. In contrast, the method reported here provides both a synergy score as well as a statistical confidence interval to qualify that synergy score.

The analysis was conducted using Windows R 2.14.1. Code is available upon request.