[Frontiers in Bioscience S2, 454-467, January 1, 2010]
Characterization of a three-drug nonlinear mixture response model
Yseult Françoise Brun1,2 , William R. Greco1,2
1
TABLE OF CONTENTS
1. ABSTRACT
Our group recently developed a response-surface modeling paradigm (White et al: Curr Drug Metab 2, 399-409, 2003) and tested its application to both mixtures of anticancer agents and antifungals. This new model is a Hill-type equation, with the slope and potency parameters being functions of the normalized drug ratios, using polynomial expressions. Response surface methods allow one to model and interpret all of the information present in the full concentration-effect data set, to visualize local regions of synergy, additivity and antagonism, and even to quantify the degree of synergy or antagonism, both globally, and across local regions of the response surface. In the present article, we study the effect of changes in the different parameters of the polynomial expressions for two-drug and three-drug mixtures, on the geometrical shapes of several 2-dimensional representations of the 3-dimensional concentration-effect surface. A secondary goal of this report is to compare the mathematical characteristics of the rival White and (Minto et al: Anesthesiol 92, 1603-1616, 2000) modeling paradigms.
2. INTRODUCTION
In the past, many varied approaches to the assessment of synergy, additivity and antagonism have been developed and applied to real data (1). They include older numerical (e.g., 2) and graphical (e.g., 3) methods and older statistical response surface methods (e.g., 4). The older response surface methods are mostly limited to two-agent interactions (e.g., 4) or have other limitations, but some newer response surface methods can be used for three-agent and higher order combinations (e.g., 5-8).
Response surface methods allow one to model and interpret all of the information present in the full concentration-effect data set, to visualize local regions of synergy, additivity and antagonism, and even to quantify the degree of synergy or antagonism, both globally, and across local regions of the response surface.
Such information from in-vitro studies can be useful to pinpoint which drug combinations may be good candidates for further in-vivo studies, and subsequent clinical studies. In addition, this information may also help to suggest optimal proportions of the agents.
In this article we use Loewe additivity (defined in (1)) as the standard for defining "no interaction" among several drugs. The response surface model studied here reduces to Loewe additivity in the complete absence of synergy or antagonism. The model has already been presented in the literature for use with both combinations of anticancer drugs (5-6) and antifungals (7). This new model is a Hill-type equation, with the steepness and potency parameters being functions of the normalized drug ratios, using polynomial expressions. For examples of how our main model (equation 1 described in Section 3) has been fit to real laboratory data, and how we have dealt with the statistical complexities of experimental design, curve fitting, model building, and results interpretation, the reader is encouraged to peruse our three published applications (5-7). The current article showcases the mathematical characteristics of the model, in order to provide the potential user with a good geometrical understanding of the model structure and parameters. This type of geometrical understanding is critical for a deep biological understanding of results from fitting real data with the model.
Minto et al. (8) published an interaction model three years before we independently, without knowledge of the former work, published the White model (5). The two models have many similarities; we will compare and contrast the Minto and White models.
Both the Minto (8) and White (5) drug combination models have advantages over older drug combination response surface models (e.g., 1, 4). Both of these newer modeling paradigms allow: (a) local regions of synergism and antagonism throughout the multidimensional surface; (b) extreme antagonism; (c) more than 2 drugs to be modeled simultaneously; (d) complex patterns of other concentration-effect parameters (Econ, B%, m; explained below in equation 1). The disadvantages of the two newer models is that their increased complexities require larger experiments, more complex model-building routines, and more complex explanations of the results.
A secondary goal of this report is to compare the mathematical characteristics of the rival Minto (8) and White (5) modeling paradigms.
3. MAIN MATHEMATICAL MODEL
The overall response surface model that we will study is essentially a 4-parameter Hill model with polynomial expressions replacing the steepness m and potency CI50 parameters:
(equation 1)
With
E: the measured endpoint (or effect or response or output or dependent variable).
Econ: the measured endpoint at zero drug concentration (control).
B%: the measured endpoint at infinite concentration(s) (the background endpoint, expressed as a percentage of Econ).
m: the steepness function (expression).
CI50: combination index expression (normalized potency expression) at the 50% level.
X represents the vector of normalized fractions of each drug. For example, for drug A,
(equation 2)
in which IC50,A represents the median effective dose (or concentration) of drug A,
and T represents the total normalized amount of drug in the mixture
(equation 3)
(as an example for a three drug mixture).
The ms and CI50s are modeled as functions of drug fractions, using the following constrained polynomials. When dealing with two-drug mixtures:
For m (equation 4)
For CI50: (equation 5)
When dealing with three-drug mixtures:
For m (equation 6)
For CI50: (equation 7)
All of the , , and terms are polynomial coefficients (estimatible parameters).
We will particularly focus on the patterns that m and CI50 functions can take depending on the values of the polynomial coefficients , , and . For each specific group of drugs, , , and form the vector θ referenced in equation 1.
4. RESULTS AND DISCUSSION
In this section, we will discuss patterns obtained by changing the various polynomial coefficients , , and , first by varying each single polynomial coefficient by a systematic increment, while fixing the other coefficients at zero; and second by varying two or more coefficients simultaneously.
For two-drug mixtures, Figures 1, 3, 5, 7, 9, 11, 13, 15, 17, 19 and 21 include plots of CI50 versus fraction of drug A, and Figures 2, 4, 6, 8, 10, 12, 14, 16, 18, 20 and 22 are corresponding plots, which are traditional isobolograms at the 50% effect level. Figures 1 to 8 show the effect of varying each of the polynomial coefficients, one at a time, with typical small to moderate effects on CI50 and isobols. Figures 9 to 12 represent more extreme patterns in CI50 and isobols caused by changes in only one polynomial coefficient. Figures 13 to 22 show complex profiles for CI50 and isobols, using sets of polynomial coefficients listed in Table 1.
Figures 23 and 24 show the effects of various values of m12 and m12 on m. Figures 25 and 26 are examples of complex patterns for m, using sets of polynomial coefficients listed in Table 2.
Finally, we simulated examples for three-drug mixtures. Figure 27 includes several ternary plots for CI50 for various D123 values. Figure 28 has different examples of complex ternary plots for CI50, for the corresponding sets of polynomial coefficients listed in Table 3. And, Figure 29 includes several ternary plots for m for various m123 values.
4.1. Two-drug mixtures
4.1.1. Characterization of CI50
First we changed the values of D1, from -10 to 10 in increments of 1, with the other polynomial coefficients in Equation 5 fixed at zero. In Figure 1, one can see the resulting plot of CI50 (log scale) versus fraction of drug A, and in Figure 2, one can see the corresponding traditional isobol plot. The CI50 plots and isobols are purely synergistic for negative values of D1 (below the CI50 = 1 additivity line) and purely antagonistic for positive values. When D1 is at -10, -1, 1 and 10, the most extreme CI50 will be at 0.23, 0.86, 1.16 and 4.40, respectively. The most extreme CI50s all occur at a fraction of drug A of 0.67. The plots are asymmetrical (skewed toward higher proportions of drug A) relative to a vertical axis at a fraction of drug A of 0.5 for the top panel, and relative to the NE to SW diagonal line for the isobol. Note that the isobol plots are not symmetrical, and that the shape of the synergistic isobols (below the NW to SE diagonal line from (0,1) to (1,0)) differs from that of the antagonistic isobols (above the line of additivity).
If we change D2 (Figures 3-4), from -10 to 10 in increments of 1, with the other polynomial coefficients in Equation 5 fixed at zero, we get the mirror images of Figures 1-2.
When we change D12 (Figures 5-6) from -10 to 10 in increments of 1, with the other polynomial coefficients in Equation 5 fixed at zero, we obtain a symmetric figure (relative to a vertical axis at a fraction of drug A of 0.5) or a symmetric isobol (relative to the NE to SW diagonal). As with the previous polynomial coefficients, the CI50 plots and isobols are purely synergistic for negative values of D12 (below the CI50 = 1 additivity line) and purely antagonistic for positive values. When D12 is at -10, -1, 1 and 10, the most extreme CI50 will be at 0.54, 0.94, 1.06 and 1.87, respectively. The most extreme CI50s all occur at a fraction of drug A of 0.5. Note that the shape of the synergistic isobols (below the NW to SE diagonal line from (0,1) to (1,0)) is not the same as the antagonistic isobols (above the line of additivity).
Finally, we change D12 (Figures 7-8) from -10 to 10 in increments of 1, with the other polynomial coefficients in Equation 5 fixed at zero. The CI50 plots and isobols are first antagonistic (above the CI50 = 1 additivity line) then synergistic (below the CI50 = 1 additivity line) for negative values of D12; and first synergistic then antagonistic for positive values. The switching point between antagonism and synergy is always for a fraction of drug A of 0.5. The absolute value of the vertical distance between the highest point of the curve (upper panel) and the horizontal additivity line (CI50 = 1) is always the same as the absolute vertical distance between the lowest point and the horizontal additivity line. When D12 is at -10/10, and -1/1, the most extreme CI50 will be at 0.98/1.02 and 0.84/1.20, respectively. The most extreme CI50s all occur at fractions of drug A of 0.28 and 0.72. Note that the isobol plots are not symmetrical, and that the shape of the synergistic isobols (below the NW to SE diagonal line from (0,1) to (1,0)) differs from that of the antagonistic isobols (above the line of additivity).
Figures 9-10 and Figures 11-12 show patterns for extreme antagonism or extreme synergy respectively, using large values of each of the polynomial coefficients, one at a time, while the other polynomial coefficients are fixed at zero. We wanted to show CI50s around 100 or 0.01as representative examples of extreme antagonism or synergy. For extreme antagonism, we had to raise D1 to 31, or D2 to 31, or D12 to 73, or D12 to 255 (or lower to -255). For extreme synergy, we had to lower D1 to -31, or D2 to -31, or D12 to -73, or D12 to -255 (or raise to 255).
Figure 13 to Figure 22 show examples of complex profiles, the lists for the sets of polynomial coefficients used are in Table 1. Figures 13-14 (curves A1 to A6) show examples of synergy with two local peaks of synergy. To obtain such profiles, we assigned negative values to D1 and D2, and large values to D12 (the sign of D12 did not matter).
Figures 15-16 (curves B1 to B6) are the complements of Figures 14-15 for antagonism.
Figures 17-18 (curves C1 to C4) show composite profiles, with both asymmetrical local synergy and antagonism. The values of the parameters, D1, D2. D12 and D12, that generated curves C1 to C4 are listed in Table 1.
Figures 19-20 (curves D1 to D4) show seesaw patterns, with antagonism, then synergy, then antagonism (D1 and D4), or synergy, then antagonism, and then synergy (D2 and D3). The values of the parameters, D1, D2. D12 and D12, that generated curves D1 toD4 are listed in Table 1.
Figures 21-22 (curves E1 and E2) also show a seesaw pattern, but a symmetrical one (relative to a vertical axis at a fraction of drug A of 0.5 for the CI50 plot, and relative to the NE to SW diagonal line for the isobol), in contrast with Figures 19-20. D1 and D2 were always of the same sign, D12 was of the opposite sign, and D12 was fixed at zero.
We could not obtain patterns that resulted in the CI50 line crossing the additivity line three times or more. If such a pattern would be needed to model real data, then it would be necessary to modify or augment the polynomials.
4.1.2. Characterization of m
m1 and m2 are simply the m's for drug A alone and drug B alone.
Figure 23 shows m versus the fraction of drug A for m1 and m2 fixed at -2, m12 fixed at 0, and m12 ranging from -10 to 7 in increments of 1. Each resulting curve is symmetrical relative to a vertical line at a fraction of drug A of 0.5, and the curves with negative m12 values are symmetrical to the curves with positive m12 values, relative to a horizontal line at an m value of -2: curves with positive m12 values are above the m = -2 line, curves with negative m12 values are below the m=-2 line. If m12 is equal to -10, the extreme value for m is -4.5; if m12 is equal to 7, the extreme value for m is -0.25. The extreme values for m are observed at a fraction of drug A of 0.5. If m12 is higher than 7, the highest value of m would be positive, which is incompatible with the kind of real data that we want to model.
Figure 24 represents m versus the fraction of drug A for m1 and m2 fixed at 2, m12 fixed at 0, and m12 ranging from -10 to 10 in increments of 1. The plots are first below the m=-2 horizontal line, then above for positive values of m12; and first above then below for negative m12 values. The switching point between above and below the m=-2 line is always for a fraction of drug A of 0.5. The absolute vertical distance between the highest point and the horizontal m=-2 line is always the same as the absolute vertical distance between the lowest point and the horizontal m=-2 line. The extreme values of m are reached for fractions of drug A at 0.21 and 0.79. If m12 is equal to 10 or -10, the extreme values for m are -1.04 and -2.96.
Figure 25 and Figure 26 are examples of complex m profiles; the values of the polynomial coefficients used are listed in Table 2.
Figure 25 (curves A1 to A4) includes examples of profiles for m ranging to more extreme values than either m for the individual drugs alone. Drug A has a steepness parameter (m) of -0.5 and drug B has an m of -2 for profiles A1 and A2; A1 has as a minimum below -2 and A2 a maximum above -0.5. Drug A has a slope of -2 and drug B a slope of -0.5 for profiles A3 and A4; A3 has as a minimum below -2 and A4 a maximum above -0.5. m12 and m12 are always of opposite signs, and m12 is negative for the two curves with minimums below -2 (A1 and A3) but positive for the two curves with maximums above -0.5 (A2 and A4).
Figure 26 (curves B1 to B4) includes examples of curves where the m shifts rapidly from the slope of one drug alone to the slope of the other drug alone.
4.2. Three-drug mixtures
4.2.1. Characterization of CI50
To represent the different values of CI50 for various values of polynomial coefficients for 3-drug mixtures, we used ternary plots. In these plots, each axis represents a normalized fraction of the drug in the mixture. The colors represent the simulated values for the CI50s; the log-scaled color scheme is shown on the left side: red indicates synergy, yellow indicates additivity, and blue indicates antagonism.
Figure 27 shows an example of the effect of various D123 values on CI50. The other polynomial coefficients are fixed at the following values: D1 =-3, D2 =12, D3 =-8, D12 =3, D13 =2, D23 =-15, D12 =20, D13 =-34, D23 =8. D123 was varied from -80 to 80 in increments of 40. The fixed polynomial coefficients have been chosen so that, on the "edges" of the ternary plots, we would have the following patterns:
- For the mixture of drug A and drug B, mostly antagonism, with small synergy for high fractions of drug A.
- For the mixture of drug A and drug C, pure synergy, but with an irregular isobol.
- For the mixture of drug B and drug C, first synergy (for low proportions of drug B) then antagonism (for high proportions of drug B), crossing the additivity line at the fraction of drug B of 0.5.
When D123 is negative and decreases, the antagonistic area (blue region) increases while the synergistic area (red region) decreases; on the contrary, when D123 is positive and increases, the antagonistic area decreases while the synergistic area increases.
Figure 28 includes five examples of representative patterns for CI50. The values of the sets of various polynomial coefficients are summarized in Table 3. and coefficients were fixed at 0 for these examples.
Example A shows additivity for the three two-drug mixtures and synergy for the three-drug mixture. We can model this with parameters fixed at zero, and a negative D123.
Example B shows additivity for the three two-drug mixtures and antagonism for the three-drug mixture. We can model this with parameters fixed at zero, and a positive D123.
Example C shows synergy for the three two-drug mixtures and antagonism for the three-drug mixture. We can model this with negative parameters, and a positive D123, higher than in example B.
Example D shows synergy for the three two-drug mixtures and even more intense synergy for the three-drug mixture. We can model this with negative parameters, and a negative D123, lower than in example A.
Example E shows antagonism for the three two-drug mixtures and synergy for the three-drug mixture. We can model this with positive parameters, and a negative D123.
4.2.2. Characterization of m
Figure 29 shows examples of the effect of various m123 values on m. m123 was varied from -40 to 40 in increments of 20.The other polynomial coefficients were fixed at the following values: m1=-0.5,
m2=-2, m3=-5, m12=-2.5, m13=2.5, m23=-10, m12=-2, m13=6, m23=6. Yellow colors indicate slopes around -2; red colors are for slopes increasing toward -0.5, and blue colors indicate values of -5 and below -5.
The fixed polynomial coefficients have been chosen so that the m of drug A alone would be -0.5, m of drug B alone would be -2, m of drug C alone would be -5, and on the "edges" of the ternary plots, we would have the following patterns:
- For the mixture of drug A and drug B, m plateaus at -2 (m of drug B alone) for low and medium fractions of drug A, and changes to -0.5 (m of drug A alone) abruptly, only for high fractions of drug A.
- For the mixture of drug A and drug C, m changes progressively from -5 (m of drug C alone) to -0.5 (m of drug A alone) with increasing fraction of drug A, but also dips to around -0.3 before coming back to -0.5.
- For the mixture of drug B and drug C, m changes rapidly from -5 (m of drug C alone) to -6.5, and then dips slowly to -2 (m of drug B alone).
When m123 is negative and decreases, the region with extreme ms (blue, ms below -5) increases, while when m123 is positive and increases, the region with ms closer to zero (red) increases.
4.3. Comparison with the Minto model
For our comparison, we considered the following polynomial equations from Minto et al. (8):
For m (Equation 8)
For CI50: (Equation 9)
Here mA and mB are the slopes for drug A and drug B alone respectively, and 2,m, 3,m, 4,m, 2,D, 3,D, and 4,D are the polynomial parameters.
For our simulations, we fixed mA and mB to -2, and we changed the six polynomial parameters mentioned above from -10 to 10 in increments of 1 when doable (we had to have negative values for m, and positive values for CI50). The simulations are shown in Figure 30 (changing 2,m), Figure 31 (changing 3,m), Figure 32 (changing 4,m), Figure 33 (changing 2,D), Figure 34 (changing 3,D), and Figure 35 (changing 4,D) respectively.
First, it is to be noted that Minto et al. suggest using a polynomial function for the parameter characterizing the endpoint (effect) at infinite drug concentrations, that we called the background, B in this article. In our datasets, whether with antifungals (7) or anticancer drugs (6), we observed that a step function for B was more appropriate than a polynomial function.
For both the m and the CI50 simulations of the Minto model, we could observe that the curves were systematically biased toward the left (lower proportions of drug A). It is not a big problem for two-drug combinations as inverting A and B as needed could take care of this, but it will be more problematic for three-drug mixtures. We observe also that no single parameter allows immediately "seesaw" patterns, with curves going from one side of the additivity line or the mA to mB line to the other. In our model, this was made possible by the gamma parameters, D12 and.m12. Finally, in the Minto et al. article (8), no specific parameter was mentioned for three-drug interactions, that could have been compared with our delta parameter.
For the polynomials for the slope m, the first parameter 0,m had to be constrained to be equal to mA, like our parameter m1. The second parameter 1,m, however, had to be constrained to be equal to mB-mA-2,m-3,m-4,m, and no parameter reflects directly mB as our parameter m1 does, making it less intuitive than our polynomial equation for steepness (sigmoidicity).
For the polynomials for CI50, the first parameter 0,D had to be constrained to be equal to 1, and the second parameter 1,D had to be constrained to be equal to -2,D-3,D-4,D so that the profiles would be anchored at 1 for the drugs alone. Our use of the leading factor (1-XA)*(1-XB) and the power function provides this anchoring, without constraining any parameter. Also, it is to be noticed that Minto's polynomial for CI50 is not logarithmic, in contrast to ours. By its inherent mathematical and statistical nature, the CI50 parameter is most probably deserving of a logarithmic transformation, before expressing it as a polynomial function. Finally, since CI50 has to have only positive values, this limits the values that the different parameters may take in the Minto model; whereas, the polynomial coefficients (parameters) from the White model have no such constraints.
4.4. Conclusion
In conclusion, our model was reasonably flexible. Regarding the potency parameter CI50, we could model simultaneous synergy and antagonism for the same two-drug mixtures, and even more so for the three-drug mixtures, with different kind of patterns, with many local pockets of synergy and antagonism. For m, we could model many different relevant geometrical shapes.
However, some limitations were found. In particular, for two agents, we could not obtain patterns that resulted in the CI50 curve crossing the additivity line three times or more. And for m, we could only barely approach a sharply changing profile.
On the other hand, our group made a decision, based upon practical experience working with several complex 3-agent datasets (5, 6), that the maximum practical complexity for the two critical polynomials should be equations 4,5 for two agents and equations 6,7 for three agents.
The decision was based on problems in interpreting the meaning of each individual coefficient in a very large set of polynomial coefficients, and the lack of practical interest by biologists in an extremely complex, but mostly empirical model. Our compromise approach is to limit the maximum complexity of the two critical polynomials. Note that also when fitting the White model to real data, polynomial coefficients with 95% confidence intervals that initially include zero will be likely removed from the model for subsequent curve fittings. Thus, for example, the final model for a three-drug mixture, may end up with less than the 10 total parameters illustrated in this article.
Overall, equations 1-3,6-7 comprise a response surface concentration-effect model that has been shown to be adequate and useful for fitting three different examples of very large, complex 3-agent data sets (5-7). The general modeling paradigm has enough flexibility to follow complex, observed patterns of sigmoidicity (m) and normalized potency (CI50) versus drug fraction. Compared to the Minto model (8), the White model appears to include material improvements. The White model should be useful for a wide spectrum of future applications. Like all models, this response surface model has the potential to be improved.
5. ACKNOWLEDGMENTS
Supported in part by NIH RR10742 (WRG), and an educational grant from Novartis through the University at Buffalo, SUNY (YFB).
6. REFERENCES
1. Greco W.R., G. Bravo, J.C. Parsons: The search for synergy: a critical review from a response surface perspective. Pharm Rev 47, 221-285 (1995)
2. Berenbaum M.C.: What is Synergy? Pharm Rev 41, 93-141 (1989)
3. Loewe S. H. Muischnek: Effect of combinations: mathematical basis of problem. Arch Exp Pathol Pharmacol 11, 313-326 (1926)
doi:10.1007/BF01952257
4. Greco W.R., H.S. Park Y.M. Rustum: An application of a new approach for the quantitation of drug synergism to the combination of cis-diamminedichloroplatinum and 1-beta-D-arabinofuranosylcytosine. Cancer Res 50, 5318-5327 (1990)
5. White D.B., H.K. Slocum, Y.F. Brun, C. Wrzosek W.R. Greco: A new nonlinear mixture response surface paradigm for the study of synergism: a three drug example. Curr Drug Metab 2, 399-409 (2003)
doi:10.2174/1389200033489316
PMid:14529372
6. White D.B., H.M. Faessel, H.K. Slocum, L. Khinkis W.R. Greco: Nonlinear response surfaces and mixture experiment methodologies applied to the study of synergism. Biometrical J 46, 56-71 (2004)
doi:10.1002/bimj.200210002
7: Brun Y.F., C.G. Dennis, W.R. Greco, R.J. Bernacki, P.J. Pera, J.J. Bushey, R.C. Youn, D.B. White B.H. Segal: Modeling the combination of Amphotericin B, Micafungin, and Nikkomycin Z against Aspergillus fumigatus in vitro using a novel response surface paradigm. Antimicrob Agents Chemother 51, 1804-1812 (2007)
doi:10.1128/AAC.01007-06
PMid:17325217 PMCid:1855564
8: Minto C.F., T.W. Schnider, T.G. Short, K.M. Gregg, A. Gentilini S.L. Shafer: Response surface model for anesthetic drug interactions. Anesthesiol 92, 1603-1616 (2000)
doi:10.1097/00000542-200006000-00017
PMid:10839909
Key Words: Synergy Additivity Antagonism, Review
Send correspondence to: William Greco, 545 Hochstetter Hall, School of Pharmacy and Pharmaceutical Sciences, University at Buffalo, SUNY, Buffalo NY, 14260, Tel: 716-645-2842, ext 232, E-mail:rosgreco@hotmail.com