[Frontiers in Bioscience E2, 582-601, January 1, 2010]
Emax model and interaction index for assessing drug interaction in combination studies
J. Jack Lee,1 Heather Y. Lin,1 Diane D. Liu,1, Maiying Kong2
1
TABLE OF CONTENTS
1. ABSTRACT
Applying the Emax model in a Lowe additivity model context, we analyze data from a combination study of trimetrexate (TMQ) and AG2034 (AG) in media of low and high concentrations of folic acid (FA). The Emax model provides a sufficient fit to the data. TMQ is more potent than AG in both low and high FA media. At low TMQ:AG ratios, when a smaller amount of the more potent drug (TMQ) is added to a larger amount of the less potent drug (AG), synergy results. When the TMQ:AG ratio reaches 0.4 or larger in low FA medium, or when the TMQ:AG ratio reaches 1 or larger in high FA medium, synergy is weakened and drug interaction becomes additive. In general, synergistic effect in a dilution series is stronger at higher doses that produce stronger effects (closer to 1-Emax) than at lower dose levels that produce weaker effects (closer to 1). The two drugs are more potent in the low compared to the high FA medium. Drug synergy, however, is stronger in the high FA medium.
2. INTRODUCTION
Due to complex disease pathways, combination treatments can be more effective and less toxic than treatments with a single drug regimen. Successful applications of combination therapy have improved treatment effectiveness for many diseases. For example, the combination of a non-nucleoside reverse transcriptase inhibitor or protease inhibitor with two nucleosides is considered a standard front-line therapy in the treatment of AIDS. Typically, a combination of three to four drugs is required to provide a durable response and reconstitution of the immune system (1). Another example is platinum-based doublet chemotherapy regimens as the standard of care for patients with advanced stage non-small-cell lung cancer (2). Combination treatments have also been shown to prevent and to overcome drug resistance in infectious diseases such as malaria and in complex diseases such as cancer (3, 4). Emerging developments in cancer therapy involve combining multiple targeted agents with or without chemotherapy, or combining multiple treatment modalities such as drugs, surgical procedures, and/or radiation therapy (5, 6).
How do we assess the effect of a combination therapy? It is a simple question, yet it requires a complex answer. A superficial way to answer the question is to determine that a combination therapy is working if its effect is greater than that produced by each single component given alone. The notion of classifying drug interaction as additive, synergistic, or antagonistic is logical and easily understood in a general sense, but can be confusing in specific application without consensus on a standard definition. Excellent reviews of drug synergism have been written by Berenbaum (7), Greco et al (8), Suhnel (9), Chou (10), and Tallarida (11), to name a few. In essence, to quantify the effect of combination therapy, we must first define drug synergy in terms of "additivity." An effect produced by a combination of agents that is more (or less) than the additive effect of the single agents is considered synergistic (or antagonistic). Then, we must further assess drug interaction in a statistical sense. Under a more rigorous definition, synergy occurs when the combined drug effect is statistically significantly higher than the additive effect. Conversely, antagonism occurs when the combination effect is statistically significantly lower than the additive effect.
Despite the controversy arising from multiple definitions of additivity or no drug interaction, the Loewe additivity model is commonly accepted as the gold standard for quantifying drug interaction (7-11). The Loewe additivity model is defined as
(E 1)
Here y is the predicted additive effect at the combination dose (d1, d2) when the two drugs do not interact. Dy,1 and Dy,2 are the respective doses of drug 1 and drug 2 required to produce the same effect y when used alone. Note that the Loewe additivity can be easily demonstrated in a "sham combination" (i.e., a drug combined with itself or its diluted form). For example, suppose drug 2 is a 50% diluted form of drug 1. The combination of one unit of drug 1 and one unit of drug 2 will produce the same effect as 1.5 units of drug 1 or 3 units of drug 2. Plugging the respective values in equation (E 1), we have 1/1.5 + 1/3 = 1. Given the dose-effect relationship for each single agent, say Ei(d)= fi(d) for agent i (i=1,2), Dy,i can be obtained by using the inverse function of fi, say, fi-1(y). Replacing Dy,1 and Dy,2 in equation (E 1) with f1-1(y) and f2-1(y), respectively, we can rewrite equation (E 1) as
(E 2)
Note that (E 2) involves an unknown variable y. By solving equation (E 2), the predicted additive effect yadd can be obtained under the Loewe additivity model. Denote that the observed mean effect is yobs at the combination dose (d1, d2). The drug combination at that dose is considered synergistic, additive, or antagonistic when the effect yobs is greater than, equal to, or less than yadd, respectively. When the dose-effect curve is decreasing (or increasing), a synergistic effect corresponds to a smaller (or larger) value than the predicted quantity.
Alternatively, to measure and quantify the magnitude of drug interaction, the interaction index (II) can be defined as
(E 3)
Note that II < 1, II =1, and II >1 correspond to the drug interaction being synergistic, additive, and antagonistic, respectively. Chou and Talalay (12) proposed the following median effect equation (E 4) to characterize the dose-effect relationship in combination studies:
(E 4)
where ED50 is the dose required to produce 50% of the maximum effect. Although the median effect equation can be applied in many settings, it assumes that when m is positive, E(d)=0 for d=0 and E(d)=1 for d=¥ . On the other hand, when m is negative, E(d)=1 for d=0 and E(d)=0 for d=¥ . If we assume that the data follow the median effect equation, a linear relationship can be found by plotting the logit transformation of the effect versus the logarithm transformed dose. A more detailed account of the interpretation and use of the interaction index can be found in a number of references (13-16). Several methods for constructing the confidence interval estimation of the interaction index were proposed by Lee and Kong (17).
To help advance research developing and comparing methods for analyzing data for combination studies, Dr. William R. Greco at the Roswell Park Cancer Institute has organized an effort and invited several groups to participate in an exercise to compare rival modern approaches to model data from two-agent concentration-effect studies. We describe the data and statistical methods, including the Emax model, and the calculation of the interaction index under the Emax model in Section 3. We describe an exploratory data analysis in Section 4, and data preprocessing for outlier rejection and standardization in Section 5. We present the main results of the data analysis in Section 6 and summarize our findings in Section 7. We close with a discussion in Section 8.
3. STATISTICAL METHODS
3.1. Data sets
Two data sets provided by Dr. Greco are used to examine the effect of the combination treatment of trimetrexate (TMQ) and AG2034 (AG) in HCT-8 human ileocecal adenocarcinoma cells. The cells were grown in a medium with two concentrations of folic acid: 2.3 m M (the first data set, called low FA) and 78 m M (the second data set, called high FA). Trimetrexate is a lipophilic inhibitor of the enzyme dihydrofolate reductase; and AG2034 is an inhibitor of the enzyme glycinamide ribonucleotide formyltransferase. The experiment was conducted on 96-well plates. The endpoint was cell growth measured by an absorbance value (ranging from 0 to 2) and recorded in an automated 96-well plate reader. Each 96-well plate included 8 wells as instrumental blanks (no cells); the remaining 88 wells received drug applications. The experiments were performed using the "ray design," which maintains a fixed dose ratio between TMQ and AG in a series of 11 dose dilutions. With 88 wells in each plate, each 5-plate stack allowed for an assessment of the combination doses at 7 curves (i.e., design rays) plus a "curve" with all controls. Two stacks were used for studying 14 design rays: TMQ only, AG only, and twelve other design rays with a fixed dose ratio (TMQ:AG) for each ray. The fixed dose ratios in the low FA experiment were 1:250, 1:125, 1:50, 1:20, 1:10, 1:5 (2 sets), 2:5, 4:5, 2:1, 5:1, and 10:1. The fixed dose ratios in the high FA experiment were 1:2500, 1:1250, 1:500, 1:200, 1:100, 1:50 (2 sets), 1:25, 2:25, 1:5, 1:2, and 1:1. Data from each of the 16 curves (2 for controls, 2 for single agents, and 12 for combinations) are grouped together. Curves 1-8 were performed on the first stack with curve 8 serving as the "control" experiment while curves 9-16 were performed on the second stack with curve 16 serving as the "control" experiment. The assignment of different drug combinations to the cells in the wells was randomized across the plates. Five replicate plates were used for each set of two stacks, resulting in a total of 10 plates for each of the two medium conditions (low FA and high FA). The maximum number of treated wells per medium condition is 880 (16 curves x 11 dilutions x 5 replicates). Complete experimental details and mechanistic implications were reported by Faessel et al (18).
3.2. Emax model
Due to a plateau of the measure of cell growth such that it does not reach zero at the maximum dose levels used in the experiments, the median effect equation (E 4) does not fit the data. Instead, we take the Emax model (19) to fit the data at hand.
(E 5)
where E0 is the base effect, corresponding to the measurement of cell growth when no drug is applied; Emax is the maximum effect attributable to the drug; ED50 is the dose level producing half of Emax; d is the dose level that produces the effect E(d), and m is a slope factor (Hill coefficient) that measures the sensitivity of the effect within a dose range of the drug. Thus, E0 - Emax is the asymptotic effect when a very large dose of the drug is applied. Figure 1 shows a few examples of the Emax model where E0 is assumed to be 1. The parameter m governs how quickly the curve drops. For the three cases in the first row in Figure 1, ED50 is fixed at 2 and Emax is at 0.8, while the slope varies. When m=1 (Figure 1.A), the dose-response curve drops slowly; when m=5 (Figure 1.B and E), a sigmoidal curve is formed, and when m=20 (Figure 1.C and F), the drop of the sigmoidal curve becomes very steep. In the three curves in the first row, as the dose increases, the curves drop, and the effect asymptotes to 1 - Emax = 0.2. In the second row, the three plots are set at Emax = 1, which means that as the dose increases, the treatment will reach the theoretical full effect. For example, if the measure of the treatment effect is cell count, all the cells will be killed at very high doses of the treatment when Emax = 1. The figures also show that, as ED50 increases, the curves shift to the right, indicating that the treatment is less potent. In all cases when m increases, the effect drops more rapidly. We apply the nonlinear weighted least squares method to estimate the parameters in the Emax model. Due to the heteroscedascity observed in the data, which means that the variance increases as the observed response increases, we use the reciprocal of the fitted response as the weight function (20). We use S-PLUS, R (21), and SAS (22) to carry out the estimation.
3.3. Interaction index under the Emax model
As when using the median effect model, the Emax model can be applied to fit the single-drug and combination drug dose-response curves, and then the interaction index can be calculated accordingly. Although equation (E 5) allows for different values of E0 and Emax for different curves, when calculating the interaction index, we need to assume all curves have the same E0 so that the "base measure" of no drug effect is the same in all curves. This can be achieved by dividing all of the effect measures with the mean of the controls. Note that Emax can vary in different curves to signify different drug potencies. However, the calculation of the interaction index will be a little more complicated when different drugs or combinations produce different values of Emax.
Hereafter, we assume the dose-response curve follows the Emax model given in (E 6):
(E 6)
The experiments we analyzed studied the ability of the combination treatments to inhibit the growth of cancer cells. The measure of the treatment effect was cell growth corresponding to the number of cells observed. Hence, the height of the dose-effect curve decreases when the dose increases. In this case, we have m > 0. In addition, as d goes to infinity, the effect plateaus at 1-Emax. Hence, Emax must be between 0 and 1.
In a study of two-drug combinations, we need to fit three curves using the Emax model: curve 1 for drug 1 alone, curve 2 for drug 2 alone, and curve c for the drug combinations. Denote Emax, i , ED50, i, and mi as the three parameters for drug i (i=1,2, c). Given an effect e (e>1-Emax), the corresponding dose d(e) can be calculated as
(E 7)
Note that the dose for the combination treatment is simply the sum of the doses of the single agents. This approach works well for a ray design with constant or varying relative potency between the two drugs (12, 17). Without loss of generality, we can assume that Emax, 1 > Emax, 2. In addition, we assume that the dose ratio for the two drugs in the combination treatment (dc=d1+d2) is fixed with d1/ d2 =p. Upon fitting the three dose-response curves, the interaction index at a fixed effect e where e Î
(1 - 
, 1) can be calculated as follows:
(E 8)
, the interaction index cannot be calculated. However, the combination effect in this range is more than additive because it reaches an effect level that no single agent can achieve alone. If Emax, 1 = Emax, 2, the interaction index can be calculated using the first formula in (E 8).
3.4. Confidence interval for the interaction index
We can apply the delta method to calculate the (large sample) variance of the interaction index (23). From our previous work (17), we found that better estimation of the confidence interval for the interaction index can be achieved by working on the logarithmic transformation of the interaction index. 
where
for i=1, 2, c.
Upon calculation of the variance for log(
), the point-wise (1-a
)100% confidence interval for the interaction index (II) for a specified effect can be constructed as
, (E 11)
where 
is the upper a
/2 upper percentile of the standard normal distribution. We also construct the simultaneous confidence band for the interaction index over the range of estimated responses. Because the estimation process involves estimating nine parameters from three curves, to construct a Scheffe-type simultaneous confidence band, we simply replace in equation (E 11) by (chi2p(a
))1/2 where p=9 (24).
3.5. Data analysis plan
The overall objective of the data analysis is to assess the synergistic effect of the combination of TMQ and AG in both low and high FA media. We apply the exploratory data analysis first, and then estimate the dose-response relationship using the Emax model. We evaluate the drug interaction by calculating the interaction index under the Loewe additivity model. We perform an exploratory data analysis in order to understand the data structure and patterns and to determine whether preprocessing of the data in terms of outlier rejection and standardization would be required prior to data modeling. We analyze the low FA and high FA experiments separately then compare the results. For each experiment, we apply the Emax model to fit the two marginal and twelve combination dose-response curves. We compute the interaction index and its 95% confidence intervals for each of the twelve combinations, and assess the overall pattern of drug interaction by examining the interaction index from the 12 fixed-ratio combinations together. We apply a one-dimensional distribution plot via the BLiP plot (25) to display the data. We use a two-dimensional scatter plot, a contour plot, and an image plot as well as a three-dimensional perspective plot to show the dose-response relationship. We also apply a trellis plot (26) to assemble the individual plots together into consecutive panels conditioning on different values of fixed dose ratios.
4. EXPLORATORY DATA ANALYSIS
As in all data analyses, we begin with an exploratory data analysis. For the low and high FA experiments, there are 871 and 879 readings, respectively. Only 9 and 1 observations, respectively, are missing out of a maximum of 880 readings in each experiment. The data include designated curve numbers ranging from 1 to 16 and data point numbers ranging from 1 to 176. Each curve number indicates a specific dose combination. We re-label the curves as A-P where A and B correspond to the control (no drug) curves; C and D correspond to the curves of TMQ and AG administered alone, and curves E through P correspond to the combination curves with fixed dose ratios in ascending order. Each point number indicates the readings at each specific dilution of each curve. Because five duplicated experiments were performed, there are up to five readings for each specific point number. There is, however, no designation of the plate number in the data received. Figure 2 shows the variable percentile plot of the distribution of the effect from the low FA and high FA experiments using the BLiP plot, with each segment corresponding to a five percent increment (25). The plot gives an overall assessment of the distribution of the outcome variable of cell growth without conditioning on experimental settings. The middle 20% of the data (40th to 60th percentiles) are shaded in a light orange color. This figure indicates that the data have a bimodal distribution with most data clustered around either a low value of 0.2 or a high value of 1.2. For the low FA experiment, the distribution of the effect ranges from 0.072 to 1.506 with the lower, middle, and upper quartiles being 0.149, 0.449, and 1.150, respectively. Similarly, for the high FA experiment, the effect ranges between 0.070 and 1.545. The three respective quartiles are 0.213, 0.990, and 1.1495. The median of the data from the low FA experiment is smaller than the median of the data from the high FA experiment. The bimodal distributions could result from steep dose-response curves. As a consequence, the slope may not be estimated well in certain cases.
To help understand the pattern of the fixed ratio dose assignment in a ray design and the relationship between the fixed ratio doses and curve numbers, we plot the logarithm transformed dose of TMQ and AG in Figure 3 for both the low FA and high FA experiments. As can be seen, curves A and B are the controls with no drugs. Curves C and D correspond to the single drug study of TMQ and AG, respectively. Curves E through P are the various fixed ratio combination doses of TMQ and AG. Note that curves J and K have the same dose ratios. Within each curve, the 11 dilutions are marked by 11 circles. For the combination studies, the curves for different dose ratios are parallel to each other on the log dose scale. If the same plot is shown in the original scale, these lines will form "rays," radiating out from the origin like sun rays. Hence, the term "ray design" is an appropriate name for this type of experiment. The corresponding dose ranges used for each drug alone are 5.47
10-6 to 0.56 mM for TMQ in both the low FA and high FA experiments, and 2.7110-5 to 2.78 mM for AG2034 in the low FA experiment and 2.7110-4 to 27.78 mM in the high FA experiment.
Figures 4 and 5 show the raw data of the effect versus dose level by curve for the low FA and high FA experiments, respectively. Instead of using the actual dose, we plot the data using a sequentially assigned dose level to indicate each dilution within each curve such that the data can be shown clearly. In addition, the data points at each dilution for each curve are coded from 1 to 5 according to the order of the appearance in the data set. We assume that these numbers correspond to the replicate number for each design point (the well position in the stack of 5 plates). Because the plate number was not listed in the data, we are not certain that this is the case. From the plot, we can see that there are outliers in several dilution series. Of note, in Figure 4, the effects from plate (replicate) #1 in curves B, E, F, and K tend to be lower than all other replicates. There are also some unusually large values, for example, in replicate 2 in curve A, dose level (dilution series) 6; replicate 3 in curve L, dose level 4; and replicate 2 in curve M, dose level 1. Similar observations can be made for the high FA experiment: plate #1 seems to have some low values in curves B, C, H, I, and J, and plate #4 seems to have some low values in curves E, K, N, O, and P. These findings indicate that certain procedures need to be performed to remove the obvious outliers in order to improve the data quality before the data analysis.
Figure 6 shows the perspective plot, contour plot, and image plot for the low FA experiment. From the perspective plots in Figure 6.A (back view), B (front view), and C (side view), we can see that the effect starts at a high plane plateau at an effect level of about 1.2 when the doses of TMQ are AG are small. As the dose of each drug increases, the effect remains approximately constant for a while and then a sudden drop occurs. This steep downward slope can be found by taking the trajectory of any combination of the TMQ and AG doses; it is also evident in the dose-response curves shown in Figures 4 and 5. The steep drop of the effect can also be found in the contour plot and the image plot. Similar patterns in the dose-response relationship are shown in Figure 7 for the high FA experiment. The steep drop of the effect occurs at smaller doses in the low FA experiment and at larger doses in the high FA experiment.
5. DATA PREPROCESSING: OUTLIER REJECTION AND DATA STANDARDIZATION
5.1. Outlier rejection
To address the concern that outliers may adversely affect the analysis outcome, we devise the following simple plan. For each of the 176 point numbers (16 curves x 11 dilutions), the five effect readings should be close to each other because they are from replicated experiments. However, because the plate number is not in the data set, we cannot assess the plate effect. Neither can we reject a certain replicate plate entirely should there be a plate with outlying data, nor apply a mixed effect model treating the plate effect as a random effect. For the four or five effect readings in each point number (only 9 point numbers in the low FA and 1 in the high FA experiments have 4 readings), we compute the median and the interquartile range. An effect reading is considered an outlier if the value is beyond the median ± 1.4529 times the interquartile range. If the data are normally distributed (i.e., follow a Gaussian distribution), the range expands to cover the middle 95% of the data. Hence, only about 5% of the data points (2.5% at each extreme) are considered outliers. The number 1.4529 is obtained by qnorm(.975)/( qnorm(.75) - qnorm(.25)) where qnorm(x) is a quantile function which returns the xth percentiles from a normal distribution. Upon applying the above rule, 129 out of 871 (14.8%) effect readings in the low FA experiment and 126 out of 879 (14.3%) in the high FA experiment are considered outliers and are removed before proceeding to further analysis. The numbers of outliers in replicates 1 to 5 are 60, 28, 19, 14, and 8 for the low FA experiment and 35, 18, 21, 34, and 18 for the high FA experiment, indicating a non-random pattern of outliers that could be attributed to experimental conditions. Note that the outlier rejection algorithm is only applied "locally." In other words, it only applies to the replicated readings up to five replicates in each of the 176 experimental conditions.
5.2. Data standardization
After outliers are removed from the data, we compute the mean of the control curves. The means for curves 8 and 16 are 1.1668 and 1.1534 for the low FA experiment and 1.1483 and 1.1477 for the high FA experiment, respectively. To apply the Emax model in equation (E 6) with E0 = 1, we standardize the data by dividing the effect readings of respective curves 1-7 by the mean of curve 8 and the effect readings of respective curves 9-15 by the mean of curve 16.
6. RESULTS
6.1. Results for the low folic acid experiment
The Emax model in equation (E 6) is applied to fit all of the dose-response curves. For the low FA experiment, the parameter estimates, their corresponding standard errors, and the residual sum of squares are given in Table 1. The dose-response relationships showing the data and the fitted curves are displayed in Figure 8. Note that although model fitting is performed on the original dose scale, the dose is plotted on the logarithmically transformed scale to better show the dose-response relationship. The fitted marginal dose-response curves for TMQ (curve C) and AG (curve D) are shown in a blue dashed line and a red dotted line, respectively. From Table 1, we see that 
is 0.00133 for TMQ and 0.00621 for AG, indicating that TMQ is about 4.7 times more potent than AG at the 
level. For curves E through P, the fitted dose-response curve for the combination treatment is shown as a solid black line superimposed on the marginal dose-response curves. The proposed Emax model fits all curves well except for curves G, H and K. For curve G, although the model estimates converge in an initial attempt, the parameter m is estimated with a standard error of 30.3. The large standard error essentially indicates that the estimate 
is not reliable. For curve K, the model does not converge on the original dose scale but converges on the logarithmically transformed dose scale. However, the standard error of the estimate 
is still very large, which leads us to believe that the model is not very stable. For curve H, as can be seen in Figure 8, there are no observed effects between 0.3 and 1 from the second to the fifth dilutions. The parameter m cannot be estimated and the model fails to converge on both the original scale and the logarithmic scale. To address these problems, we conclude that the data do not provide us sufficient information to yield a reasonable estimate of the parameter m. Therefore, we take a remedial approach by fixing m, and then proceed to estimate the other two parameters. Upon checking the data, we set the parameter m as 5, 4.5, and 5 for curves G, H, and K, respectively. The choice of m is somewhat arbitrary with a goal of yielding a good fit to the data and producing a small residual sum of squares. The resulting "reduced" models fit the data reasonably well but with a consequence that there is no standard error estimate for
, which affects the variance estimation of the interaction index (to be shown later). Based on limited sensitivity analysis, the estimation of the interaction index remains reasonably robust.
In all dose-response curves, the standardized effect level starts to drop between dose levels (dilutions) 3 to 6. Once the effect starts to drop, it drops quickly and plateaus at the 1 - 
level. There are ample data points at the effect levels around 1 (dose levels 1-4) and 1 - 
(dose levels 8-11). However, due to the sharp drop in the dose-response curves, fewer data points can be found in the middle of the effect range. When the data points become too few or do not spread out to cover enough range, it becomes harder for the model to converge, as seen in curves G, H, and K. The overall results for the curve fitting of the low FA experiment are that the values of 
range from 0.863 to 0.890; the values of
range from 0.00133 to 0.00621; and the values of 
range from 1.971 to 5.473. The residual sum of squares ranges from 0.0599 to 0.1025 without large values, suggesting that the model fits the data reasonably well.
Based on the fitted dose-response curve, the interaction index (II) can be calculated over the entire effect range and at specific dose combinations. Table 2 gives a detailed result of the estimated interaction index and its 95% point-wise confidence interval at each dose combination for each combination curve. The II is calculated at the predicted effect level from the combination curve and not at the observed effect level. The results are shown in a trellis plot in Figure 9 where the red lines represent the 95% point-wise confidence intervals at each specific effect level and the black dashed lines indicate the 95% simultaneous confidence bands of the II for the entire range. From the figure we find that the interaction index can be estimated with very good precision in all curves except at the two extremes when the effect is close to 1 or 1 - 
. The trend and the pattern of the interaction index are clearly shown in these figures. For curves E through K, i.e., with a TMQ:AG dose ratio ranging from 0.004 to 0.2, synergy is observed in the effect range between 0.2 to 0.9. For curves L and M, which have TMQ:AG ratios of 0.4 and 0.8, we see that synergy is observed at the low effect level from 0.2 to about 0.5. Beyond 0.5 the combinations are generally additive. For curves N, O, and P, with TMQ:AG ratios of 2, 5, and 10, the synergistic effect is lost and we see additivity in all dose ranges.
6.2. Results for the high folic acid experiment
Table 3 lists the parameter estimate, corresponding standard error, and sum of squares for all the curves in the high FA experiment. Unlike in the experiment using low FA media, the model fitting for all curves in the high FA experiment converge when using the Emax model. The estimated 
ranges from 0.831 to 0.893; 
ranges from 0.0137 to 0.1943 except for curve D (AG alone with
= 0.5224); and 
ranges between 1.468 and 3.625. The residuals sum of squares ranges from 0.0615 to 0.1134. Compared to the low FA experiment, the values of
are greater in the high FA experiment, indicating that the drugs are less potent when applied to a high FA medium. Note that the doses of TMQ are the same between the two experiments but the doses of AG are 10 times higher in the high FA experiment. In addition, 
= 0.0137 and 0.00133 for TMQ alone in the high and low FA experiments, respectively, which indicates that the drug is 10 times less potent in the high FA medium compared to the low FA medium. The potency of AG is even more dramatically reduced. In Figure 10 we see that the Emax model provides an excellent fit to all the curves. Table 4 gives a detailed account of the interaction index in all dilutions for all of the combination curves. The results are summarized in a trellis plot in Figure 11. Again, the red lines represent the 95% point-wise confidence intervals at each specific effect level and the black dashed lines correspond to the simultaneous confidence bands of the II for the whole range. Using the high FA medium, synergy can be achieved for most of the drug combinations in all the effect ranges, with the exception of the very low or very high effect ranges. The confidence intervals are still very tight although they are a little wider compared to their counterparts from the low FA experiment. As the TMQ:AG ratio increases from 0.0004 to 0.5, synergy is observed across all dilution series. In addition, higher synergy is observed at the lower effect levels, particularly when the TMQ:AG is at 0.01 or lower (curves E, F, G, H, and I). In the middle effect levels (effects between 0.2 and 0.8), the II ranges from about 0.1 in curves J and K, to 0.12 in curve L, 0.15 in curve M, 0.25 in curve N, and 0.35 in curve O. The higher the TMQ:AG ratio, the less synergy it achieves. In curve P, for example, when the TMQ:AG ratio reaches 1, synergy is lost.
7. SUMMARY
In both the low FA and high FA experiments, TMQ is more potent than AG. At low TMQ:AG ratios, i.e., when a small amount of the more potent drug (TMQ) is added to a larger amount of the less potent drug (AG), synergy is achieved. However, when the TMQ:AG ratio reaches 0.4 or larger for the low FA medium, or when the TMQ:AG ratio reaches 1 or larger for the high FA medium, synergy decreases, or the interaction becomes additive. In general, a synergistic effect in a drug combination dilution series is stronger at higher doses that produce stronger effects (effects closer to 1-Emax) than at lower dose levels that produce weaker effects (effects closer to 1). The two drugs in this study are more potent in the low FA medium compared to the high FA medium. The drug synergy, however, is stronger in the high FA medium.
8. DISCUSSION AND PERSPECTIVE
The data supplied by Dr. Greco provide an excellent opportunity to apply and compare various approaches for studying the effects of combination drug treatments. For the median effect model, a linear relationship between the logit transformed effect and the log-dose makes the model fitting straightforward and easy. However, when measuring cell growth, as in the experiments we analyzed, if the maximum drug effect reaches a plateau and does not kill all the cancer cells, even at the highest experimental doses, the median effect model (12) does not apply. We used the Emax model (19), which provides an adequate fit for most data. Parameter estimation under the Emax model requires the use of iterative procedures such as the nonlinear weighted least squares method, which can address the heteroscedascity problem. Model convergence is not guaranteed; whether or not the model converges depends on the data and the choice of the initial values. We find that PROC NLIN in SAS provides a more comprehensive and robust environment for estimating parameters with nonlinear regression compared to the nls() function in S-PLUS/R. It can be useful to apply SAS first to estimate the parameters and then feed the results into S-PLUS/R for further data analysis and production of graphics. Unlike fitting the linearly-transformed median effect model via linear regression, for which a solution can always be found, fitting the Emax model via nonlinear regression may result in nonconvergence of the model in some cases. This nonconvergence may indicate aberrant conditions in the data such that the data do not provide adequate information for model fitting. We had convergence problems with the curves G, H, and K in the low FA experiment. In these cases, there were insufficient data in the middle of the effect range; hence, the parameters could not be estimated reliably. We had to fix the m parameter before we could estimate the other two parameters. From the dose-response curves, we found that TMQ was more potent than AG, and that the drug combination was more potent in the low FA medium than in the high FA medium.
Upon construction of the marginal and combination dose-response curves, we applied the Loewe additivity model to compute the interaction index. We note that a definition of drug interaction such as the interaction index is model dependent. Additionally, no matter which model is used, based on the definition of the interaction index (7,8), the dose levels used in calculating the interaction index must be translated back to the original units of dose measurement. Under the given model, we found that the drug interaction between TMQ and AG was largely synergistic. Synergy was more clear and evident in the high FA experiment than in the low FA experiment. In addition, synergy was more likely to be observed when a small dose of the more potent drug (TMQ) was added to a large dose of the less potent drug (AG). When a large amount of a more potent drug is present, adding the less potent drug does not show synergy because the effect is already largely achieved by the more potent drug. In addition, the interval estimation showed that the 95% confidence intervals were wider at the two extremes of the effect, which were closer to 1 or to 1-Emax. This result is consistent with that of many regression settings in which estimation achieves higher precision in the center of the data distribution but lower precision at the extremes.
We have provided a simple, yet useful approach for analyzing drug interaction for combination studies. The interaction index for each fixed dose ratio is computed and then displayed together using a trellis plot. This method works well for the ray design. Other methods have been proposed to model the entire response surface using the parametric approach (27) or the semiparametric approach (28). The results from applying the semiparametric model are reported in a companion article (29).
9. ACKNOWLEDGEMENT
We thank Dr. William R. Greco at the Roswell Park Cancer Institute for organizing this project comparing rival modern approaches to analyzing combination studies, for supplying the data sets, and for the invitation to present this manuscript. The authors also thank Lee Ann Chastain for her editorial assistance. This work is supported in part by grants W81XWH-05-2-0027 and W81XWH-07-1-0306 from the Department of Defense, and grant CA16672 from the National Cancer Institute. J. Jack Lee's research was supported in part by the John G. & Marie Stella Kenedy Foundation Chair in Cancer Research.
10. REFERENCES
1. R. Chris Rathbun, Staci M. Lockhart, Johnny R. Stephens: Current HIV treatment guidelines - an overview. Curr Pharm Des 12, 1045-1063 (2006)
doi:10.2174/138161206776055840
2. Julian R. Molina, Alex A. Adjei, James R. Jett: Advances in chemotherapy of non-small cell lung cancer. Chest 130, 1211-1219 (2006)
doi:10.1378/chest.130.4.1211
3. Myaing M. Nyunt, Christopher V. Plowe: Pharmacologic advances in the global control and treatment of malaria: combination therapy and resistance. Clin Pharmacol and Ther 82, 601-605 (2007)
doi:10.1038/sj.clpt.6100361
4. Kamalesh K. Sankhala, Kyriakos P. Papadopoulos: Future options for imatinib mesilate-resistant tumors. Expert Opinion Investig Drugs 16, 1549-1560 (2007)
doi:10.1517/13543784.16.10.1549
5. Roberto Bianco, Vincenzo Damiano, Teresa Gelardi, Gennaro Daniele, Fortunato Ciardiello, Giampaolo Tortora: Rational combination of targeted therapies as a strategy to overcome the mechanisms of resistance to inhibitors of EGFR signaling. Curr Pharm Des 13, 3358-3367 (2007)
doi:10.2174/138161207782360537
6. Sith Sathornsumetee, David A. Reardon, Annick Desjardins, Jennifer A. Quinn, James J. Vredenburgh, Jeremy N. Rich: Molecularly targeted therapy for malignant glioma. Cancer 110, 12-24 (2007)
doi:10.1002/cncr.22741
7. Morris C. Berenbaum: What is synergy? Pharmacol Rev 41, 93-141 (1989)
8. 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)
9. Jurgen Suhnel: Parallel dose-response curves in combination experiments. Bull Math Biol 60, 197-213 (1998)
doi:10.1006/bulm.1997.0009
10. 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
11. Ronald J. Tallarida: An overview of drug combination analysis with isobolograms. J Pharmacol Exp Ther 319, 1-7 (2006)
doi:10.1124/jpet.106.104117
12. Ting-Chao Chou, Paul Talalay: Quantitative analysis of dose effect relationships: the combined effects of multiple drugs or enzyme inhibitors. Adv Enzyme Reg 22, 27-55 (1984)
doi:10.1016/0065-2571(84)90007-4
13. Stella G. Machado, Grant A. Robinson: A direct, general approach based on isobolograms for assessing the joint action of drugs in pre-clinical experiments. Stat Med 13, 2289-2309 (1994)
doi:10.1002/sim.4780132202
14. Chris Gennings: On testing for drug/chemical interactions: definitions and inference. J Biopharm Stat 10, 457-467 (2000)
doi:10.1081/BIP-100101978
15. Kathryn S. Dawson, Walter H. Carter, Jr., Chris Gennings: A statistical test for detecting and characterizing departures from additivity in drug/chemical combinations. J Agric Biol Environ Stat 5, 342-359 (2000)
doi:10.2307/1400458
16. J. Jack Lee, Maiying Kong, Gregory D. Ayers, Reuben Lotan: Interaction index and different methods for determining drug interaction in combination therapy. J Biopharm Stats 17, 461-480 (2007)
17. 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
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. Naitee Ting: Dose finding in drug development. Springer, NY, 127-145 (2006)
20. John M. Chambers, Trevor J. Hastie: Statistical models in S. Chapman & Hall/CRC Press, Boca Raton FL, 450 - 452 (1992)
21. Sylvie Huet, Anne Bouvier, Marie-Anne Poursat, Emmanuel Jolivet: Statistical tools for nonlinear regression: a practical cuide with S-PLUS and R examples, 2nd ed. Springer, NY (2003)
22. SAS/STAT 9.1 User's guide. SAS Institute, Cary, NC (2007)
23. Peter J. Bickel, Kjell A. Doksum: Mathematical statistics: basic ideas and selected topics, Vol I, 2nd ed. Prentice Hall, New Jersey pp. 306-314 (2000)
24. Christopher Cox, Guangpeng Ma: Asymptotic confidence bands for generalized nonlinear regression models. Biometrics 51, 142-150 (1995)
doi:10.2307/2533321
25. J. Jack Lee, Z. Nora Tu: A versatile one-dimensional distribution plot: the BLiP plot. Amer Statistician 51(4), 353-358 (1997)
doi:10.2307/2685905
26. Michael J. Crawley: The R book. Wiley, NY (2007)
doi:10.1002/9780470515075
27. 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
28. 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
29. Maiying Kong, J. Jack Lee: Applying Emax model and bivariate thin plate splines to assess drug interactions. Front Biosci E2, 279-292, (2010)
Abbreviations: AG: AG2034, an inhibitor of the enzyme glycinamide ribonucleotide formyltransferase; ED50: dose required to produce 50% of the maximum effect; Emax: maximum effect attributed to the drug; FA: folic acid; II: interaction index; TMQ: Trimetrexate, a lipophilic inhibitor of the enzyme dihydrofolate reductase;
Key Words: additivity, antagonism, confidence interval estimation, Emax model, Loewe additivity model, nonlinear regression, synergy, trellis plot.
Send correspondence to: J. Jack Lee, Department of Biostatistics, Unit 1411, The University of Texas M. D. Anderson Cancer Center, P. O. Box 301402, Houston, Texas 77230-1402, Tel: 713-794-4158, Fax: 713-563-4242, E-mail:jjlee@mdanderson.org