Applying E_max model and bivariate thin plate splines to assess drug interactions

Maiying Kong¹ , J. Jack Lee²

1Department of Bioinformatics and Biostatistics, School of Public Health and Information Sciences, University of Louisville, Louisville, Kentucky 40292, U.S.A.,²Department of Biostatistics, The University of Texas M. D. Anderson Cancer Center, Unit 1411, 1515 Holcombe Boulevard, Houston, Texas 77030, U.S.A.

TABLE OF CONTENTS

1. Abstract

2. Introduction

2. Introduction

3. Statistical methods

3.1. Estimating dose-effect curves

3.2. Predicting additive effects

3.3. Assessing drug interactions using bivariate thin plate splines

4. Case studies

4.1. Case study 1: cancer cells grown in a medium with 2.3 .μM folic acid (low FA experiment)

4.2. Case study 2: cancer cells grown in a medium with 78 μM folic acid (high FA experiment)

5. Summary and perspective

6. Acknowledgement

References

1. ABSTRACT

We review the semiparametric approach previously proposed by Kong and Lee and extend it to a case in which the dose-effect curves follow the E_max model instead of the median effect equation. When the maximum effects for the investigated drugs are different, we provide a procedure to obtain the additive effect based on the Loewe additivity model. Then, we apply a bivariate thin plate spline approach to estimate the effect beyond additivity along with its 95% point-wise confidence interval as well as its 95% simultaneous confidence interval for any combination dose. Thus, synergy, additivity, and antagonism can be identified. The advantages of the method are that it provides an overall assessment of the combination effect on the entire two-dimensional dose space spanned by the experimental doses, and it enables us to identify complex patterns of drug interaction in combination studies. In addition, this approach is robust to outliers. To illustrate this procedure, we analyzed data from two case studies.

2. INTRODUCTION

Studies of interactions among biologically active agents, such as drugs, carcinogens, or environmental pollutants have become increasingly important in many branches of biomedical research. For example, in cancer chemotherapy, the therapeutic effect of many anticancer drugs is limited when they are used as single drugs. Finding combination therapies with increased treatment effect and decreased toxicity is an active and promising research area (1). An effective and accurate evaluation of drug interaction for in vitro and/or in vivo studies can help to determine whether a combination therapy should be further investigated.

The literature supports the Loewe additivity model as the gold standard for defining drug interactions (2-5). The Loewe additivity model defines an additive effect based on the following equation,

<  (E 1)

Here y is the predicted additive effect, which is produced by the combination dose (d₁, d₂) when the two drugs do not interact; and D_y,1 and D_y,2 are the respective doses of drug 1 and drug 2 required to produce the same effect y when applied alone. If we know the dose-effect relationship for each single agent, say E(d)=f_i(d) for agent i (i=1,2), we are able to obtain the dose D_y,i by using the inverse function of f_i, denoted as f_i^-1(y). By replacing D_y,1 and D_y,2 in equation (E 1) with f₁^-1(y) and f₂^-1(y), respectively, we can obtain an equation that includes the single variable y, i.e.,

<  (E 2)

By solving equation (E 2), we can obtain the predicted additive effect y. If the observed effect at (d₁, d₂) is more than (equal to, or less than) the predicted effect, we say that the combination dose (d₁, d₂) is synergistic (additive, or antagonistic).

In our previous studies (6-8), we found that Chou and Talalay's (9) median effect equation was appropriate to describe the dose-effect relationships. Chou and Talalay's median effect equation, in its nonlinear form, can be written as follows:

<  (E 3)

where ED₅₀ is the dose required to produce 50% of the maximum effect, and m is the slope factor (Hill coefficient), measuring the sensitivity of the effect to the dose range of the drug. For the data in the case studies (see Section 4 for details), we found that the median effect equation (E 3) could not adequately describe the marginal dose-effect relationship because the plateau of the effect does not go to zero when a large dose level of a drug is applied. Instead, the following E_max model (E 4) presented by Ting (10) describes the dose-effect relationship very well:

< (E 4)

In the E_max model (E 4), E₀ is the base effect, corresponding to the measurement of response when no drug is applied; E_max is the maximum effect attributable to the drug; ED₅₀ is the dose level producing half of E_max, i.e., ED₅₀ is the dose level required to produce the effect at a value of E₀-0.5E_max (Figure 1.A); d is the dose level, which produces the effect E. Thus, E₀-E_max will be the asymptotic net effect when a large dose of the drug is applied. Different maximum effects for agents may reflect different mechanisms of action for the drugs (11). For in vitro studies, cell growth is commonly used as an endpoint to measure the effect of inhibitors. When no drugs (or, no inhibitors) are applied, the cell proliferation obtains its largest value. In this case, the dose-effect curve is similar to the one shown in Figure 1.A, where E_max>0. The effect range determined by the dose-effect curve lies between (E₀-E_max, E₀), and the asymptotic measurement for the maximum drug effect is E₀-E_max.

In the investigation of drug interactions, theoretically, we expect the measurements for the endpoints to be similar when no drug is applied. We use the measurements that are made without any drugs applied as controls. However, we realize that environmental factors other than the experimental conditions may lead to different measurements for the controls under different environments. Thus, we may need to standardize the observed effects by the mean of the control for each environmental condition (1,6), and then take E₀=1. In this paper, we consider the following dose-effect curve for each drug:

<  (E 5)

which assumes an effect at value 1 when no drug is applied. Once we obtain the dose-effect curve for each single drug, we can use the Loewe additivity model (E 1) to obtain the additive effect for any combination dose, particularly, for the combination dose with observed effects. Thus, we may obtain the differences in observed effects and the predicted additive effect at each observed combination dose. We use the bivariate thin plate splines approach (12) to estimate the relationship between these differences and the combination doses. Consequently, we obtain a response surface of the differences over the combination doses, and can construct 95% confidence surfaces of the response surface. When the dose-response curves decrease with increasing dose, an observed effect that is smaller in magnitude than the prediction of Loewe additivity implies that the observed effect is stronger than the predicted effect, indicating that the combination dose is synergistic. Conversely, an observed effect that is larger in magnitude than the prediction of Loewe additivity implies that the observed effect is weaker than the predicted effect, indicating that the combination dose is antagonistic. However, these inferences should be made based on sound statistical considerations. Based on the fitted response surface and its upper and lower confidence surfaces, we can judge whether the difference is significantly less than zero, not different from zero, or greater than zero, and we can determine the patterns of drug interaction in terms of synergy, additivity, and antagonism. This paper is organized as follows. In Section 3.1, we describe the underlying stochastic assumption for the dose-effect curve and the procedure to estimate the parameters in each marginal dose-effect curve. In Section 3.2, we explain how we obtained the additive response surface based on the Loewe additivity model, in particular for studies in which the maximum effects of the drugs are different. In Section 3.3, we explain how we assessed the response surface beyond the additivity surface and how we constructed its 95% confidence surfaces. These procedures allow us to identify drug interactions in terms of synergy, additivity, or antagonism for all of the combination doses in the region containing the combination design points. In Section 4, we illustrate the procedure introduced in Section 3 by analyzing real data in two case studies. The last section is devoted to a short summary and perspective.

3. STATISTICAL METHODS

Assume that the observed data are (d_1i, d_2i, E_i) for i=1, ... , n. For each i, (d_1i, d_2i) is the observed combination dose and E_i is the corresponding observed effect. When only a single drug is applied (drug 1 or drug 2), we refer to the observations as marginal observations. That is, the marginal observations for drug 1 are the observations (d_1i, d_2i, E_i) with d_2i=0 (i=1, ... , n), and the marginal observations for drug 2 are the observations (d_1i, d_2i, E_i) with d_1i=0 (i=1, ... , n). The marginal dose-effect curves are estimated based on the marginal observations, which we present in Section 3.1. It is commonly accepted that the additive effect should be obtained based on the dose-effect relationships for each individual drug. In Section 3.2, we explain how we obtained the predicted effect at combination dose (d₁, d₂) based on the Loewe additivity model (E 1) and the marginal dose-effect curves (E 5). We denote the predicted effect as<. By definition, there is no drug interaction when only a single drug is applied. Therefore, the term for drug interaction is meaningful only for the combination dose (d₁ ,d₂) with nonzero d₁ and d₂. In Section 3.3, we develop a procedure to estimate the effect beyond additivity for any combination dose (d₁ ,d₂) with nonzero d₁ and d₂, denoted by <

3.1. Estimating dose-effect curves

Chou and Talalay (9), Chou (4), and Kong and Lee (6) estimated the parameters in the median effect equation (E 3) by using the transformation < and applying the least squares method in the linear regression setting, where α=-m log( ED₅₀). The case studies we evaluated (see Section 4) included many low doses, the effects of which are larger than 1 after setting the effect at the control level to be 1. Thus, a similar transformation for models (E 3) and (E 5) cannot be carried out. Because the measurements are continuous, we propose applying nonlinear least squares regression to estimate the parameters in models (E 3) and (E 5) with the assumption that a stochastic error with N(0, σ²) exists on the right-hand side of the two models. We note that estimating the dose-effect curve for drug i requires only the marginal observations for drug i with i=1, 2. We apply nonlinear least squares regression to estimate the parameters in the marginal dose-effect curves in the two case studies (shown in Section 4).

3.2. Predicting additive effects

We obtain the predicted effect based on the Loewe additivity model (E 1) when model (E 5) is applied as the marginal dose-effect curve for each drug. When model (E 5) is applied, the dose required to produce effect E is given by

< 

However, the maximum effects for the two drugs may be different. Without a loss of generality, we assume that the maximum effect of drug 1 is larger than the maximum effect of drug 2, i.e., E_max,1>E_max,2. For this case, when the dose-effect curves are decreasing, neither drug applied alone can produce an effect in (0, 1- E_max,1) (Figure 1.B). Based on the Loewe additivity model (E 1), we can see that the predicted effect will be in the interval of (1- E_max,1, 1) for any combination dose (d₁, d₂).

Recall that the Loewe additivity model (E 1) can be rewritten as <, and the ratio < (denoted as ρ(y)), is often called the relative potency of drug 2 versus drug 1 at effect level y, which means that the effect of 1 unit of drug 2 produces the same effect as ρ(y) units of drug 1. Generally speaking, the relative potency ρ(y) is dose-dependent (7). When there is no drug interaction, the effect of the combination dose (d₁, d₂) produces the same effect as drug 1 alone at dose level D_y,1, which equals d₁+ρ(y)d₂, or drug 2 alone at dose D_y,2 , which equals ρ(y)^-1d₁+d₂ (Figure 2.A). All the combination doses (d₁, d₂) on the line < have the predicted effect y, where <is the line connecting the points P=( D_y,1, 0) and Q=(0, D_y,2 ) (Figure 2.A). This line < is often called an additive isobole (3, 4).

When E_max,1>E_max,2, as illustrated in Figure 1.B, we can calculate the dose of drug 1 required to produce the maximum effect of drug 2, i.e., < Note that the range of the effect for drug 2 is (1- E_max,2 ,1), which could be produced by drug 1 alone at a dose level between 0 and <. Based on the Loewe additivity model, for any level of effect y in (1-E_max,2 ,1), the associated additive isobole is the line connecting (D_y,1, 0) and (0, D_y,2). When y varies from 1-E_max,2 to 1, the dose of drug 1 required to produce effect y varies from<to 0, while the dose of drug 2 required to produce effect y varies from infinitely large to 0. In particular, when y is close to 1- E_max,2, the dose of drug 1 required to produce effect y is close to<, and the dose of drug 2 required to produce effect y goes to infinity. Figure 2.B shows four typical additive isoboles (dashed lines), which connect equally effective doses of drug 1 and drug 2 at different effect levels. From left to right, the effect level decreases in magnitude. The additive isoboles may not be parallel because the relative potency may not be constant. When y varies in (1-E_max,2 ,1), all the additive isoboles cover the region between the two solid vertical lines (Figure 2.B). Meanwhile, any combination dose (d₁, d₂) with <must lie on one of these isoboles. Therefore, for any combination dose (d₁, d₂) with<, the predicted additive effect, say y, can be obtained by solving the following nonlinear equation for E:

<

Now we examine the predicted effect for the combination dose (d₁, d₂) with<. When <, drug 1 alone at dose d₁ produces an effect<, an effect beyond 1- E_max,2, which cannot be produced by drug 2 alone at any dose level. In this case, if the effect of the combination dose is more than the effect produced by drug 1 alone, then drug 2 potentiates the effect of drug 1, and synergy occurs because the predicted additive effect is the effect produced by drug 1 alone at dose level d₁. Alternatively, because drug 2 alone cannot produce such an effect, we could consider D_y,2 to be infinitely large. Thus, the Loewe additivity model is reduced to d₁/D_y,1=1, and the predicted additive effect is the effect produced by drug 1 alone at D_y,1. No matter which approach we take, the predicted effect y is the same. We can determine the predicted effect y from the following equation:

< 

Thus, we can obtain the predicted effect for any combination dose (d₁, d₂). Using notation similar to that from the previous study of Kong and Lee (6), we denote the predicted effect at the combination dose (d₁, d₂) as <. In the following subsection, we develop a procedure to estimate the effect beyond additivity, denoted by< and to construct its 95% point-wise confidence interval and simultaneous confidence interval. We assess the drug-drug interaction based on the estimated< and its confidence intervals.

3.3. Assessing drug interactions using bivariate thin plate splines

By definition, there is no drug interaction when a single drug is used alone. Therefore, we set the differences between the observed and predicted effects at zero for the marginal observations, that is, for the combination doses (d₁, d₂) with only one nonzero component. We apply a bivariate thin plate spline to estimate the differences as a function of the combination dose, say, f(d₁, d₂). When the dose-effect curves are decreasing, f(d₁, d₂)<0 indicates that the observed effect is more than the predicted effect at (d₁, d₂), thus the combination dose (d₁, d₂) is synergistic. Inversely, f(d₁, d₂)>0 indicates that the combination dose (d₁, d₂) is antagonistic. Kong and Lee (6) used the different observed combination doses as the knots for the bivariate thin plate splines (12). The choice of knots works well when the number of combination doses is not large and the combination doses are not close, such as those from factorial designs or uniform designs (13). However, when ray designs are applied and the doses are low, the combination doses are very close and some columns of the design matrix (i.e., Ω and Z₁ in the following notations) may be highly correlated, which results in a nearly singular matrix for estimating the parameters in the function f . If that happens, a low rank smoothing thin plate spline (14) should be applied to avoid the singularity of the involved matrix due to the low rank of the design matrix. An example of such a low rank smoothing thin plate spline is the knots formed by selecting the observed combination doses with distances larger than some pre-specified small number. Alternatively, one may use an appropriate transformation of the dose, such as the log-transformation, to evenly distribute the experimental combination doses in certain regions in order to improve the ability to estimate the effect beyond additivity using bivariate thin plate splines.

Suppose the selected knots are < (k=1,...,K) , then the bivariate thin plate spline can be expressed by the following form:

<

where < and < are the parameters in the thin plate spline function f, and

< The distance in the expression is the Euclidean distance. Let us denote

<

The coefficient γ and <can be obtained by minimizing the following penalized residual sum of squares:

< subject to <. (E 6)

Following the notation by Kong and Lee (6) and Green and Silverman (12), consider a QR decomposition of T^T, say T^T=FG, where F is a <orthogonal matrix and G is a <upper triangular matrix. Let F₁ be the first three columns of F, and F₂ be the last K-3 columns of F. We can show that <if and only if <can be expressed as F₂ξ for some ξ. Thus, the penalized residual sum of squares can be expressed as

<.

Set < and < where < is the matrix square root of <. The penalized residual sum of squares can be expressed as

<  (E 7)

Based on the approach proposed by Ruppert, Wand, and Carroll (15) and Wang (16), and detailed by Kong and Lee (6) in this setting, the parameters in terms of < and u can be obtained by solving the following mixed effect model:

< with <. (E 8)

Thus, the parameters can be estimated by < with <, <, and D=diag(0, 0, 0, 1, ..., 1), where the number of zeros in the matrix D corresponds to the number of < (i=0,1,2) and the number of ones corresponds to the number of < Under these notations, for any combination dose (d₁, d₂), f(d₁, d₂) can be predicted by < with <, and an approximate 100(1-α)% point-wise confidence interval for f(d₁, d₂) can be constructed by

< (E 9)

where C_d=(1, d₁, d₂, Z₀) and< is the upper < percentile of the standard normal distribution. Thus, we can construct a 95% point-wise confidence band for f = 0 by taking the intercept lines of the lower and upper confidence surfaces with the dose plane. We then claim that the combination doses in the area outside the confidence bound with f < 0 are synergistic, the combination doses inside the bound are additive, and the combination doses in the area outside the bound with f > 0 are antagonistic.

Based on the 95% point-wise confidence intervals constructed from (E 9), some combination doses that are additive may be claimed as synergistic or antagonistic. To be conservative and to control the family-wise error rate, we also construct the 95% lower and upper simultaneous confidence surfaces, which are based on a format similar to that of equation (E 9) except that < is replaced by < (17), where EDF is the effective degrees of freedom from the resulting bivariate smoothing splines (12) and is defined as the trace of the matrix <, and < is the upper < percentile of the F distribution with EDF and n-EDF degrees of freedom. Here n is the total number of observations except controls. A 95% simultaneous confidence band for f = 0 can be formed by taking the intercept lines of the 95% lower and upper simultaneous confidence surfaces with the dose plane. For the two case studies presented in the next section, we will present the plots of different patterns of drug interaction based on the respective 95% point-wise confidence intervals and the 95% simultaneous confidence band.

4. CASE STUDIES

The following two data sets were provided by Dr. William R. Greco. The data were collected during a study of the joint effect of trimetrexate (TMQ) and AG2034 on cells grown in medium with different concentrations of folic acid (FA): 2.3 m M in the first experiment (the low FA experiment), and 78 m M in the second experiment (the high FA experiment). TMQ is a lipophilic inhibitor of the enzyme dihydrofolate reductase, and AG2034 is an inhibitor of the enzyme glycinamide ribonucleotide formyltransferase. The unit of drug concentration is the micromole (m M) for all data analyzed in the case studies. The endpoint was the growth of human ileocecal adenocarcinoma (HCT-8) cells in 96-well assay plates as measured by the sulforhodamine B (SRB) protein stain. The drug treatments were randomly assigned to the cells in the assay wells. Each 96-well plate included 8 wells as instrumental blanks (no cells); thus 88 wells were used for drug treatments. Five replicate plates were used for each set of 88 treated wells. Each of the two large data sets were obtained from two 5-plate stacks with a maximum of 880 treated wells per experiment. Each experiment included 110 control wells, in which no drugs were applied to the cells. Ray designs were used, with the experimental doses being distributed in 14 rays, including two rays for TMQ and AG2034 when used alone. The complete details and mechanistic implications of the study were reported by Faessel et al (18). Assuming that the first observation recorded in each dose or combination dose from the first 5-plate stack was from the same plate, say the 1^st plate, the second observation from the 2^nd plate, and so on, and also assuming that the first observation recorded in each dose or combination dose from the second 5-plate stack was from the same plate, say the 6^th plate, the second observation from the 7^th plate, and so on, we have a total of 10 plates for each of the two data sets.

To examine whether there is a significant difference among the plates, we applied one-way analysis of variance (ANOVA) to the controls in each data set. The p-values were 0.001 for the low FA experimental data and 0.005 for the high FA experimental data. The results indicate a significant plate effect among the 10 plates for each experiment, that is, the inter-plate variability is high. To attenuate the effect from the inter-plate variability, we applied a standardization procedure to each data set, dividing the effect readings by the mean of the controls in each associated plate. Thus, the mean for the controls within each plate was standardized to 1, and the effect for the controls was treated as 1. In addition to 110 controls for each experiment, the data included 761 observations for the low FA experiment and 769 observations for the high FA experiment. We applied the statistical method described in Section 3 to each of the two standardized data sets. We present the results for each experiment in the following two subsections.

Lee et al (19) performed extensive exploratory analyses on the same data sets and identified 129 outliers out of 871 (14.8%) effect readings in the low FA experiment and 126 outliers out of 879 (14.3%) effect readings in the high FA experiment. To compare our findings with the results previously obtained (19), we removed the outliers from the data and then again applied the statistical method described in Section 3. For each experiment, we report the detailed analyses of the original data set and of the modified data set that excluded the outliers.

4.1. Case study 1: cancer cells grown in a medium with 2.3 m M folic acid (low FA experiment)

In this experiment, the cells were grown in a medium with 2.3 m M folic acid. We fitted marginal dose-effect curves for TMQ and AG2034 by using both the median effect equation (E 3) and the E_max model (E 5). The dose levels for TMQ when applied alone were 5.47<10^-6, 4.38<10^-5, 1.38<10^-4, 4.38<10^-4, 8.75<10^-4, 1.75 <10^-3, 3.5<10^-3, 7<10^-3, 2.21<10^-2, 7 <10^-2, and 0.56 m M, and the dose levels for AG2034 when applied alone were 2.71<10^-5, 2.71<10^-4, 6.87<10^-4, 2.17<10^-3, 4.3<10^-3, 8.7 <10^-3, 1.74<10^-2, 3.48<10^-2, 0.11, 0.3475, and 2.78 m M. Note that some effect readings at low doses or combination doses are greater than 1, thus, the logit transformation could not be carried out. We applied nonlinear least squares regression to estimate the parameters in models (E 3) and (E 5). Figures 3.A and 3.B show the respective fitted marginal dose-effect curves for TMQ and AG2034 with the dose levels shown on a log scale. The dotted-dashed lines are the curves based on the median effect model (E 3), and the solid lines are the dose-effect curves based on the E_max model (E 5). From the fitted dose-effect curves, we found that the E_max model (E 5) provided a much better fit than the median effect equation for the marginal data. Therefore, we chose the E_max model (E 5) to describe the dose-effect relationship in this case study. The parameters estimated for TMQ and AG2034 are shown in the three columns under the title "Low FA" in Table 1. Here, the estimate of E_max,TMQ is slightly larger than the estimate of E_max,AG2034. We plotted the distribution of the combination doses using the original scale (not shown) and found that most of the combination doses were crowded in the region of the low doses, which could cause a singularity of the involved matrices due to the low rank of Ω and Z₁ used for estimating the effect beyond additivity when using bivariate thin plate splines (see Section 3.3). Hence, we applied a log transformation of the form log(dose+δ) for each dose level, where δ is a small number, say 2.74<10^-6, half of the smallest dose level for the two drugs when applied alone. We plotted the distribution of the combination doses on the log(dose+δ) scale, as shown in Figure 3.C. The points on the horizontal line in Figure 3.C are the doses of TMQ on the log(dose+δ) scale; the points on the vertical line are the doses of AG2034 on the log(dose+δ) scale; and the points on each of the remaining 12 design rays are the combination doses at each ray with each dose component on the log(dose+δ) scale. The 12 design rays for combination doses appearing left to right in Figure 3.C correspond to the combination doses at 12 ratios of TMQ to AG2034, i.e., 1:250, 1:125, 1:50, 1:20, 1:10, 1:5, 1:5, 2:5, 4:5, 2:1, 5:1, 10:1. The 12 rays are denoted by the letters E, F, G, H, I, J, K, L, M , N, O, P, which represent the respective curves 15, 13, 11, 7, 5, 3, 9, 4, 6, 10, 12, 14 in the original data set for the low FA experiment. Note that rays 3 and 9, denoted by J and K, indeed have the same fixed dose ratio. To obtain the predicted additive effects, we applied the procedure described in Section 3.2, keeping the dose levels on the original scale. The contour plot of the predicted additive effect is shown in Figure 3.D. Note that the effect levels for TMQ applied alone are obtained from (1-E_max,TMQ, 1), which is (0.1190, 1), and the effect levels for AG2034 applied alone are obtained from (1-E_max,AG2034, 1), which is (0.1312, 1). The vertical line with contour level 0.13 is the predicted effect produced by TMQ alone. The plot of the differences of the observed effects and predicted effects versus the dose levels of AG2034 on log(dose+δ) scale is shown in Figure 3.E. That plot shows that the differences are not distributed around zero, rather, for some observations of AG2034, the differences are significantly less than zero, in the range of (-7, -4) on the log (dose+ δ) scale with δ=2.74<10^-6, i.e., in the range of 0.001 μM to 0.018 μM on the original dose scale. Therefore, the pure additive effect model could not describe the data well. We used bivariate thin plate splines to fit the differences versus the transformed doses, with the knots at all the distinct transformed dose levels. The transformation is taken as log (dose+ δ) for all single and combination doses. By convention, there is no drug interaction when a single drug is applied. Therefore the differences were set to zero for the marginal doses. Applying the bivariate thin plate splines (Section 3.3), we obtained < < Next, we constructed 95% point-wise upper and lower confidence surfaces for the fitted bivariate spline function f(d₁, d₂) based on equation (E 9). Figure 3.F shows the contour plot of the fitted spline function f(d₁, d₂) at the levels of -0.1, 0, and 0.1 as thin solid lines; the intercept lines of the corresponding 95% point-wise upper confidence surface with the dose plane as thick dashed lines; and the intercept lines of the corresponding 95% point-wise lower confidence surface with the dose plane as thick solid lines. The combination doses inside the thick dashed curves, shown in light blue, are synergistic because the effects beyond additivity at these combination doses are significantly smaller than zero. The combination doses inside the thick solid curves, shown in light pink, are antagonistic because the effects beyond additivity at these combinations were significantly larger than zero. The combination doses in the uncolored region, which lie between the thick solid curves and the thick dashed curves, are additive because the effects beyond additivity are not significantly different from zero. In particular, the combination doses with AG2034 in the transformed scale in the range of (-7, -4) inside the thick dashed lines are synergistic, which is consistent with the residual plot in Figure 3.E. The fitted response surface was obtained by adding the fitted spline function f (i.e., the effect beyond additivity) to the predicted additive surface. The contour plot of the fitted response surface at the contour levels of 0.2, 0.5, and 0.9 is shown in Figure 3.I. The final residuals were obtained by subtracting the fitted effects from the observed effects. The plots of final residuals versus the dose levels of TMQ and AG2034 on the log (dose+ δ) scale are shown in Figures 3.G and 3.H, respectively. From these two plots, we see that the residuals are centered around zero along the experimental dose range. We conclude that the model fits the data reasonably well.

To examine the patterns of drug interactions in different rays and different experimental combination doses, we combined Figures 3.F and 3.I, that is, we plotted the contour curves of the fitted response surface at the levels of 0.2, 0.5, and 0.9 in Panel F. We also plotted the representative design rays, with the experimental combination doses shown as dots on these rays (Figure 4.A), As seen in that figure, the combination doses on the rays E through K (curves 15, 13, 11, 7, 5, 3, 9 in the original data set) are synergistic when the effect levels are between 0.9 and a number smaller than 0.2. The combination doses on these rays are additive when the effect level is less than this small number, and the combination doses at low levels on these lines are either additive or antagonistic. The combination doses on the rays N, O, and P (curves 10, 12, and 14 in the original data set) are additive when the effects are less than 0.9, and the combination doses at low dose levels are antagonistic.

In addition to the 95% point-wise confidence interval, we constructed the 95% simultaneous confidence band for f=0 with <, and < = 12.20, where n=761, EDF=119, and <. The resulting patterns of drug interactions are shown in Figure 4.B, in which the thick dashed line is the intercept line of the 95% upper simultaneous confidence surface with the dose plane. Based on Figure 4.B, we conclude that the combination doses inside the thick dashed curves, shown in light blue, are synergistic. The combination doses outside the thick dashed curves are additive. As seen in the figure, compared to the point-wise confidence interval approach, the simultaneous confidence band method shrinks the synergistic area and results in the disappearance of the antagonistic area. A point-wise confidence interval is appropriate for making inferences for each observed combination dose. The simultaneous confidence band is suitable for making a global assessment; however, it can be overly conservative.

Following those analyses, we evaluated the low FA experimental data set from which we had removed the outliers (19). The results of our assessment of drug interactions for this data set are presented in Figures 4.C and 4.D. Figures 4.A and 4.C contain parallel information, as do Figures 4.B and 4.D. Comparing the plots across the panels, we conclude that the results from fitting the original data set and those from fitting the data set excluding outliers are very similar. Therefore, the semiparametric method presented in Section 3 is robust to outliers in this example.

We recommend using caution when considering extrapolations based on spline estimations. The fitted response surface for the differences between the observed effects and predicted effects gives an overall picture of the drug interactions (see Figures 4.A and 4.B). However, the fitted results on the two larger areas outside experiment rays E and P should not be over-interpreted because (i) there are no experimental data in those areas and (ii) we forced the differences of the observed effects and predicted additive effects to be zero at the marginal observed dose levels.

4.2. Case study 2: cancer cells grown in a medium with 78 m M folic acid (high FA experiment)

In the high FA experiment, the dose levels for TMQ when applied alone were 5.47<10^-6, 4.38<10^-5, 1.38<10^-4, 4.38<10^-4, 8.75<10^-4, 1.75 <10^-3, 3.5<10^-3, 7<10^-3, 2.21<10^-2, 7 <10^-2, and 0.56 m M, and the dose levels for AG2034 when applied alone were 2.71<10^-4, 2.17<10^-3, 6.87<10^-3, 2.17<10^-2, 4.34<10^-2, 8.68 <10^-2, 1.74<10^-1, 3.47<10^-1, 1.1, 3.47, and 27.8 m M. The procedure we used to analyze this data set was the same as that used in case study 1. By applying nonlinear least squares regression, we estimated the marginal dose-effect curves using the median effect equation (E 3) (dotted-dashed lines) and the E_max model (E 5) (solid lines), (shown in Figures 5.A and 5.B). It is clear that the E_max model (E 5) fits the data better than the median effect equation, thus, we chose the E_max model (E 5) as the dose-effect curve for this data set. The estimated parameters for the marginal dose-effect curves for the E_max model are shown in the three columns under the title "High FA" in Table 1. The combination doses on the original scale (not shown) are crowded in the region of the low doses, thus we applied the transformation in the form of log(dose+δ) to each dose level, where δ is a small number, say 2.74<10^-6, one half of the lowest dose level for TMQ and AG2034 when applied alone. The distribution of the experimental dose levels on the log(dose+δ) scale is shown in Figure 5.C. The 12 design rays for the combination doses correspond to the 12 dose ratios of TMQ versus AG2034 at 1:2500, 1:1250, 1:500, 1:200, 1:100, 1:50, 1:50, 1:25, 1:12.5, 1:5, 1:2, and 1:1, which are denoted by the letters E, F, G, H, I, J, K, L, M , N, O, and P, representing the curves 15, 13, 11, 7, 5, 3, 9, 4, 6, 10, 12, 14 in the original data set for the high FA experiment. Applying the procedure described in Section 3.2, we obtained the contour plot of the predicted additive effect that is shown in Figure 5.D. We see that the contour line at level 0.15 is the predicted effect produced by TMQ alone because AG2034 could not produce such an effect when applied alone; the effect levels for AG2034 applied alone range from 0.1816 to 1. Figure 5.E shows the differences of the observed effects and predicted effects versus the dose levels of AG2034 on the log(dose+δ) scale. That plot shows that the differences are not centered around zero, rather for some observations of AG2034, the differences are significantly less than zero, in the range of (-5, 0) on the log(dose+δ) scale, i.e., in the range of 6.7<10^-3 μM to 1.0 μM on the original dose scale. These findings indicate that some combination doses were synergistic and that the pure additive effect model could not describe the data well. We used bivariate thin plate splines to fit these differences versus the transformed doses or combination dose, with the knots at all the distinct transformed dose levels. The transformation is taken as log (dose+ δ) for all single doses and combination doses. We constructed its 95% point-wise confidence surfaces based on equation (E 9). The estimated <, < Figure 5.F shows the contour plot of the fitted spline function f at the levels of -0.1, 0, and 0.1 as thin solid lines; the intercept lines of its corresponding 95% point-wise upper confidence surface with the dose plane as thick dashed lines; and the intercept lines of its corresponding 95% point-wise lower confidence surface with the dose plane as thick solid lines. The combination doses inside the thick dashed curves, shown in light blue, are synergistic; the combination doses inside the thick solid curves, shown in light pink, are antagonistic; and the combination doses in the uncolored area are additive. We obtained the fitted response surface by adding the fitted spline function f to the predicted additive surface (shown in Figure 5.I). The plots of the final residuals versus the dose levels of TMQ and AG2034 on the log(dose+δ) scale are shown in Figures 5.G and 5.H, respectively. From these two plots, we see that the residuals are centered around zero along the experimental dose range, indicating that the model describes the data reasonably well.

To examine the patterns of drug interactions in different rays and different experimental combination doses, we combined the plots in Figures 5.F and 5.I to form Figure 6.A (as we did when analyzing the data from the low FA experiment). From Figure 6.A, we see that the combination doses on all 12 rays are synergistic when the effect levels are between 0.9 and 0.15. The combination doses at high dose levels are additive, and most of the combination doses at low dose levels are additive. In addition, we constructed a 95% simultaneous confidence band based on equation (E 9) with <replaced by <. Here EDF=91, n=769, and <=10.77. The results are presented in Figure 6.B, in which the thick dashed line is the intercept line of the upper 95% simultaneous confidence surface with the dose plane. Based on Figure 6.B, we conclude that the combination doses inside the thick dashed curves, shown in light blue, are synergistic. The combination doses outside the thick dashed curves are additive. As in our analysis of the data from case study 1, in this analysis, we found the simultaneous confidence band to yield more conservative results and to be more suitable for a global assessment. For this case study, we also assessed the data set from which we had removed the outliers. The results for assessing drug interactions are presented in Figures 6.C and 6.D. Comparing Figures 6.A to C and 6.B to D, we conclude that the results from fitting the original data set and those from fitting the data set excluding outliers are very similar. Thus, the results indicate that the semiparametric method is robust to outliers.

5. SUMMARY AND PERSPECTIVE

We extended the approach proposed by Kong and Lee (6) to a situation for which the E_max model is more appropriate to describe the marginal dose-effect relationship. We considered the possibility that some effect readings at low doses may fall beyond the mean of the controls. Under such circumstances, the standardized effect is greater than 1 and a logit transformation to a linear model (4, 8, 9) cannot be carried out. Hence, other models such as the E_max model are needed and nonlinear least squares regression methods can be applied for estimating parameters for the dose-effect curves. We applied nonlinear least squares regression in the case studies to estimate the parameters for the dose-effect curves specified by the median effect equation and the E_max model.

Additionally, we extended the approach of Kong and Lee (6) as a solution to the problem arising when the experimental points are very close and the low rank of the design matrix may cause computational problems in matrix inversion. For this situation, we considered a low-rank thin plate spline (14) to estimate the surface beyond additivity, and, alternatively, we applied an appropriate transformation to the doses so that the combination doses on the transformed scale were more evenly distributed. In the case studies, we first applied the transformation log(dose+δ) to each component of the combination doses and then applied bivariate thin plate splines with knots representing all the different observed doses on the log(dose+δ) scale. In both case studies, we chose δ as half of the smallest non-zero dose among TMQ and AG2034 when applied alone, that is, δ=2.74<10^-6 for both experiments. The value of δ selected should not be too small or too large compared with the magnitude of the dose levels. An extremely small δ results in a relatively large distance between the marginal doses and combination doses. Conversely, a large δ dominates in the transformation log(dose+δ) when the dose levels are low. From the final residual plots, it is evident that the current transformation works well.

It is well known that the smoothing parameter λ governs the trade-off between the goodness-of-fit and the smoothness of the function f. When λ becomes larger, the fitted function f tends to be smoother and the residuals tend to be larger. The selection of the smoothing parameter plays a key role in the fitted results. In the case studies, the smoothing parameter, <, was selected as <, which is almost identical to the selected smoothing parameter based on the generalized cross validation criterion and "leave-out-one" cross validation criterion. For example, for the low FA experimental data, the selected parameters based on the mixed model approach used in this paper, cross validation criterion, and generalized cross validation criterion were 0.0178, 0.0112, and 0.0071, respectively. For the high FA experimental data, the corresponding selected parameters were 0.0842, 0.0842, and 0.0531, respectively. Indeed, Kohn, Ansley, and Tharm (20) showed that the estimation of the smoothing parameter based on a mixed model approach is comparable with the standard method of the generalized cross validation criterion. By applying a mixed effects model, the smoothing parameter can be automatically determined by <. This method has been implemented in S-PLUS by Ruppert et al (15) using the lme function (21). In our previous study (1), based on extensive simulations, we showed that the selection of the smoothing parameter provides a good estimate to the underlying function in general.

In the two case studies, we also performed the same analyses for the two reduced data sets analyzed by Lee et al (19), and achieved almost identical results, which indicates that this semiparametric method is robust to outliers. The semiparametric method we have developed can also be used to assess drug interactions for the combination doses not on the design rays, and to identify complex patterns of drug interaction in combination studies. In addition, the method gives an overall assessment of the combination effect in the entire two-dimensional dose space spanned by the experimental doses with a caveat that extrapolation beyond data points can be risky.

We also note that the estimated function f(d₁,d₂) and its 95% confidence surfaces can guide the exploration of whether some parametric models are sufficient to describe the data. Many parametric models have been proposed in the literature. Greco, Bravo, and Parsons (3) provided an excellent review of the response surface approach. However, without prior knowledge of the response surface model, or adequate representation of the data by a parametric model, most parametric approaches will fail. Blindly using any parametric model can be dangerous and may lead to wrong conclusions of drug interactions. In our proposed approach, there is no need to assume any parametric models for f(d₁,d₂). We provide a promising approach by modeling the mixture effect data with spline techniques via a mixed-effect model. We advocate the use of the semiparametric method for model building because the true patterns of drug interactions are typically not known. The conclusions regarding drug interactions are based on the estimated f and its confidence surfaces, which are determined by the underlying data. The S-PLUS code used to evaluate the case studies can be obtained from the first author.

6. ACKNOWLEDGEMENT

We thank Dr. William R. Greco at Roswell Park Cancer Institute for supplying the data sets and for his invitation to submit this paper. We thank the two referees for their constructive comments. MK is supported in part by a National Institutes of Health grant P20RR024489. JJL is supported in part by the Department of Defense grants W81XWH-05-2-0027 and W81XWH-07-1-0306 and a National Cancer Institute grant CA16672.

7. REFERENCES

1. Fumihiko Kanzawa, Kazuto Nishio, Kazuya Fukuoka, Toshihiko Sunami, and Nagahiro Saijo: In vitro interactions of a new derivative of spicamycin, KRN5500, and other anticancer drugs using a three-dimensional model. Cancer Chemother Pharmacol 43, 353-363 (1999)
doi:10.1007/s002800050907
PMid:10100589

2. Morris C. Berenbaum: What is synergy? Pharmacol Rev 41, 93-141 (1989)

3. William R. Greco, Gregory Bravo, John C. Parsons: The search of synergy: A critical review from a response surface perspective. Pharmacol Rev 47(2), 331-385 (1995)

4. Ting-Chao Chou: Theoretical basis, experimental design, and computerized simulation of synergism and antagonism in drug combination studies. Pharmacol Rev 58, 621-681 (2006)
doi:10.1124/pr.58.3.10
PMid:16968952

5. J. Jack Lee, Maiying Kong, Gregory D. Ayers, Reuben Lotan: Interaction index and different methods for determining drug interaction in Combination Therapy. J Biopharm Stat 17, 461-480 (2007)

6. Maiying Kong, J. Jack Lee: A semiparametric response surface model for assessing drug interactions. Biometrics 64, 396-405 (2008)
doi:10.1111/j.1541-0420.2007.00882.x
PMid:17900314

7. Maiying Kong, J. Jack Lee: A general response surface model with varying relative potency for assessing drug interactions. Biometrics 62 (4), 986-995 (2006)
doi:10.1111/j.1541-0420.2006.00579.x
PMid:17156272

8. J. Jack Lee, Maiying Kong: A confidence interval for interaction index for assessing multiple drug interaction. Stat Biopharm Res 1, 4-17 (2009)
doi:10.1198/sbr.2009.0001

9. Ting-Chao Chou, Paul Talalay: Quantitative analysis of dose effect relationships: the combined effects of multiple drugs or enzyme inhibitors. Adv Enzyme Regul 22, 27-55 (1984)
doi:10.1016/0065-2571(84)90007-4
PMid:6382953

10. Naitee Ting: Dose Finding in Drug Development. Springer, New York, USA, 127-145 (2006)
doi:10.1007/0-387-33706-7
PMid:AMBIGUOUS (161 citations)    PMCid:AMBIGUOUS (22 citations)

11. Daniel M. Jonker, Sandra A. G. Visser, Piet H. van der Graaf, Rob A. Voskuyl, and Meindert Danhof: Towards a mechanism based analysis of pharmacodynamic drug-drug interaction in vivo. Pharmacol and Ther 106, 1-18 (2005)
doi:10.1016/j.pharmthera.2004.10.014

12. Peter J. Green, Bernard W. Silverman: Nonparametric Regression and Generalized Linear Models. Chapman & Hall, London (1994)

13. Ming Tan, Hong-Bin Fang, Guo-Liang Tian, and Peter J. Houghton: Experimental design and sample size determination for testing synergism in drug combination studies based on uniform measures. Stat Med 22, 2091-2100 (2003)
doi:10.1002/sim.1467
PMid:12820275

14. Haonan Wang, and M. Giovanna Ranalli: Low-rank smoothing splines on complicated domains. Biometrics 63, 209-217 (2007)
doi:10.1111/j.1541-0420.2006.00674.x
PMid:17447947

15. David Ruppert, Matthew P. Wand, and Raymond J. Carroll: Semiparametric Regression. Cambridge University Press, UK (2003)

16. Yuedong Wang: Mixed effect smoothing spline analysis of variance. J Roy Stat Soc B 60, 159-174 (1998).
doi:10.1111/1467-9868.00115

17. Henry Scheffe: The Analysis of Variance. New York: John Wiley & Sons (1959)

18. Helene M. Faessel, Harry K. Slocum, Robert C. Jackson, Theodore J. Boritzki, Youcef M. Rustum, M. G. Nair, and William R. Greco: Super in vitro synergy between inhibitors of dihydrofolate reductase and inhibitors of other folate-requiring enzymes: The critical role of polyglutamylation. Cancer Res 58, 3036-3050 (1998)

19. J. Jack Lee, Heather Y. Lin, Diane D. Liu, and Maiying Kong: Applying Emax model and interaction index for assessing drug interaction in combination studies. Front Biosci (2009) (in press)

20. Robert Kohn, Craig F. Ansley, David Tharm: The performance of cross validation and maximum likelihood estimators of spline smoothing parameters. J Am Stat Assoc 86, 1042-1050 (1991)
doi:10.2307/2290523

21. Jose C. Pinheiro, Douglas M. Bates: Mixed-Effects Models in S and S-PLUS. New York: Springer-Verlag (2000)

Abbreviations: ANOVA = analysis of varianceFA = folic acidHCT-8 cells = human illeocecal adenocarcinoma cellslme = linear mixed effects�M = micromoleTMQ = trimetrexateSRB = sulforhodamine B

Key Words: Additivity, Antagonism, Bivariate splines, Combination therapy, E_max model, the Loewe additivity model, Review, Synergy, Review