Response Surface Model for Anesthetic Drug Interactions

Click on the links below to access all the ArticlePlus for this article.

Please note that ArticlePlus files may launch a viewer application outside of your web browser.

• https://links.lww.com/ALN/A47

DRUG interactions are the basis of anesthetic practice. For example, induction of anesthesia may consist of intravenous administration of a benzodiazepine before induction, a hypnotic to achieve loss of consciousness, and an opioid to blunt the response to noxious stimulation. Similarly, anesthesia often is maintained with a combination of a hypnotic (e.g. , propofol, isoflurane) and an analgesic (e.g. , fentanyl, nitrous oxide). Anesthetic drugs are often combined because they interact synergistically to create the anesthetized state.

Pharmacodynamic drug interactions are typically described using mathematical models. The basic model is that of an isobole. Isoboles are iso-effect curves, curves that show dose combinations that result in equal effect. 1The combination of two doses (d¹and d²) can be represented by a point on a graph, the axes of which are the dose axes of the individual drugs (fig. 1). The isobole connects isoeffective doses of the two drugs when administered alone, D¹and D². If the isobole is straight (fig. 1A), then the relation is additive. If the isobole bows toward the origin (fig. 1B), then smaller amounts of both drugs are needed to produce the drug effect when administered together, so the relation is supraadditive or synergistic. If the isobole bows away from the origin (fig. 1C), then greater amounts of both drugs are needed to produce the drug effect when administered together, so the relation is infraadditive.

In table 1we propose a set of criteria that pharmacodynamic models of drug interactions should meet. In this article we propose an interaction model that meets these criteria, based on response-surface methodology. Response surfaces are a powerful statistical methodology for estimating and interpreting the response of a dependent variable to multiple inputs. 2Response-surface methodology is used for two principal purposes; to provide a description of the response pattern in the region of the observations studied and to assist in finding the region in which the optimal response occurs. Our model is a straightforward extension of the sigmoidal concentration–response relation for individual drugs. We test the proposed model using data from a study of the interaction of midazolam, propofol, and alfentanil with loss of consciousness. 3 

This article only considers pharmacodynamic interactions, the type of interaction most relevant to the practice of anesthesia. Pharmacokinetic interactions are entirely different and will not be considered. Appendix 1 (which can be found on the Anesthesiology Web site at www.anesthesiology.com) reviews several commonly used pharmacodynamic models of drug interactions and shows areas in which existing models fail to meet the criteria in table 1.

The effects of individual drugs are often modeled by relating drug effect (E) to drug concentration (C) using a sigmoid model:

where E⁰is the baseline effect when no drug is present, E^maxis the peak drug effect, C⁵⁰is the concentration associated with 50% drug effect, and γ is a “sigmoidicity factor” that determines the steepness of the relation.

This relation is shown graphically in figure 2. The concentration term often is defined as the concentration at the site of drug effect, but the model can be generalized to any measure of exposure (e.g. , dose, plasma concentration, or area under the curve). For models of probability, such as the probability of moving in response to surgical incision, E⁰is 0 and E^maxis the maximal probability (usually assumed to be 1).

Dividing the numerator and denominator of equation 1by C^50γ, we obtain an alternate form:

In this model, concentration has been normalized to the concentration that results in 50% of maximal drug effect. This is a natural way to think about drug concentration—as a fraction of some measure of potency. For example, anesthesiologists are accustomed to thinking about volatile anesthetics in terms of minimum alveolar concentration (MAC), rather than in absolute concentration terms. This is precisely the concept of normalizing drug concentration to potency.

The basic concept of our proposed interaction model is simple. Consider two drugs, each of which has a sigmoidal concentration–response relation. We will think of any given ratio (i.e. , B/(A + B), called θ herein) of the two drugs as behaving as a new drug. This new drug, which is actually a fixed ratio of the two drugs, has its own sigmoidal concentration–response relation, as shown in figure 3. This is the basic premise of our interaction model. The mathematics are simply an extension of the model for a single drug to a model that considers each ratio of two drugs as a drug in its own right.

We will express the concentrations of drugs A and B as [A] and [B]. As suggested by equation 2, we must first normalize each drug to its potency, C⁵⁰, and express the results in units (U) of potency.

where U^Ais the normalized concentration of drug A, and U^Bis the normalized concentration of drug B. We can define a family of “drugs,” each being a unique ratio of U^Aand U^B. Each drug will be defined in terms of θ, where θ is defined as

By definition, θ ranges from 0 (drug A only) to 1 (drug B only). The “drug concentration” is simply U^A+ U^B. We can extend equation 2to describe the concentration–response relation for any ratio, θ, of the two drugs in combination:

where θ is the ratio of the two drugs, the drug concentration is U^A+ U^B, γ(θ) is the steepness of the concentration–response relation at ratio θ, U⁵⁰(θ) is the number of units (U) associated with 50% of maximum effect at ratio θ, and E^max(θ) is the maximum possible drug effect at ratio θ. Because E^max, C⁵⁰, and γ in equation 2have been replaced by functions of θ, each ratio has the potential to have its own E^max, C⁵⁰, and γ. This allows each ratio of drug A and drug B to behave as its own drug, with its own sigmoidal concentration–response relation, which is the basic premise of the model.

The term “U⁵⁰(θ)” is the potency of the drug combination at ratio θ relative to the normalized potency of each drug by itself. This requires careful explanation. Let us assume that only drug A is present, in a concentration of C^50,A. In this case, the drug effect is half of the maximal effect, U^A= 1, U^B= 0, θ= 0, and the drug concentration is U^A+ U^B= 1. Because we have 50% of the maximum drug effect, and 1 unit of drug, then the number of units associated with 50% drug effect when only drug A is present, U⁵⁰(0), must be 1. Similarly, let us assume that only drug B is present and the concentration of drug B is C^50,B. In this case, the drug effect is half of the maximal effect, U^A= 0, U^B= 1, θ= 1, and the drug concentration is U^A+ U^B= 1. Because we have 50% of the maximum drug effect, and 1 unit of drug, then the number of units associated with 50% drug effect when only drug B is present, U⁵⁰(1), must again be 1. By definition, if only drug A or drug B is present, U⁵⁰(θ) = 1.

Now, let us assume that drug A and drug B both are present, each in exactly half of the concentration that would cause 50% of the drug effect when administered alone. In this case, U^A= 0.5, U^B= 0.5, θ= 0.5, and the drug concentration is U^A+ U^B= 1. If this causes 50% of maximum effect, then the drugs are simply additive at θ= 0.5, and U⁵⁰(0.5) = 1. However, if this combination produces more than a half-maximal effect, then 1 unit of this combination, at θ= 0.5, is more potent than 1 unit of either drug alone (i.e. , synergistic). In this case, U⁵⁰(0.5) < 1. Conversely, if this combination produces less than a half-maximal effect, then 1 unit of this combination, at θ= 0.5, is less potent than either drug alone (i.e. , infraadditive). In this case, U⁵⁰(0.5) > 1. Thus, U⁵⁰(θ) is the potency of the combination compared with the potency of either drug alone, which is 1 by definition.

Thus, the units of U⁵⁰(θ) are not concentration units, but rather the number of units, at ratio θ, associated with 50% of maximal drug effect. U⁵⁰(θ) is 1 for θ= 0 and θ= 1. For all values of θ between 0 and 1 (i.e. , all possible ratios of the two drugs), U⁵⁰(θ) assumes a value determined by the data. If this value is 1, then the interaction is additive at θ. If the value is less than 1, then the drug effect is synergistic at θ. If the value is greater than 1, then the interaction is antagonistic at θ.

Figure 4shows the relation between a three-dimensional response surface and a conventional two-dimensional isobolographic analysis. The two-dimensional isobologram is a cut through the three-dimensional surface, generally taken at the 50% response level. In this particular example, synergy is evident in the three-dimensional model as a bowing of the surface toward the reader. This bowing causes the conventional isobologram to deviate toward the origin from the straight line of additivity.

Much pharmacodynamic literature supports the sigmoid relation in equation 1, equation 2, and equation 5. There is only modest information specifying the functions E^max(θ), U⁵⁰(θ), and γ(θ). Our choice is to use functions that are capable of taking a variety of shapes, so that good approximations to the true relations can be determined empirically. To provide these flexible functions we chose fourth-order polynomials of the form

where f(θ) is E^max(θ), U⁵⁰(θ), or γ(θ). The coefficients (β⁰, β¹, β², β³, β⁴) are model parameters that are either constrained by the model or estimated from the data. Fortunately, two of these terms, β⁰and β¹, can be replaced by other terms already defined.

We already defined the values E^max(θ), U⁵⁰(θ), and γ(θ) when only drug A is present, E^max,A, U^50,A, and γ^A, respectively. Note in equation 6that when θ= 0 (only drug A is present), f(0) =β⁰. Therefore, when f(θ) is E^max(θ), U⁵⁰(θ), or γ(θ), β⁰must be E^max,A, U^50,A, and γ^A, respectively.

Similarly, we also defined the values E^max(θ), U⁵⁰(θ), and γ(θ) when only drug B is present, E^max,B, U^50,B, and γ^B, respectively. Referring again to equation 6, when θ= 1 (only drug B is present), f(1) =β⁰+β¹+β²+β³+β⁴. We can rearrange this as β¹= f(1) −β⁰−β²−β³−β⁴. Thus, when f(θ) is E^max(θ), U⁵⁰(θ), or γ(θ), β¹must be E^max,B− E^max,A−β^2,Emax−β^3,Emax−β^4,Emax, U^50,B− U^50,A−β^{2, U50}−β^3,U50−β^4,U50, or γ^B−γ^A−β^2,γ−β^4,γ, respectively.

This permits us to develop models that incorporate the individual drug parameters for E^max(θ), U⁵⁰(θ), and γ(θ) as functions of θ. The equation for E^max(θ), using the substitutions previously mentioned for β⁰and β¹, is

U^50,Aand U^50,B, [equivalent to U⁵⁰(θ) and U⁵⁰(1)], are both 1 by definition. Thus, when f(θ) = U⁵⁰(θ), the values of β⁰and β¹in equation 6are 1 and −β²−β³−β⁴, respectively. Therefore, the equation for potency as a function of θ can be simplified to

Many isobolograms have a simple inward or outward curvature, which can be readily encompassed with a simple quadratic form of equation 8with just one coefficient:

If β^2,U50is 0, then the value of U⁵⁰(θ) will be 1 for all values of θ. This means that the interaction will be additive. If β^2,U50is a positive number, then U⁵⁰(θ) will be less than 1 for all values of θ between 0 and 1. The effect is to magnify the term

in equation 5, making it appear that there is more drug present. This will produce a greater than additive effect, i.e. , synergy. If β^2,U50is a negative number, then U⁵⁰(θ) will be greater than 1 for all values of θ between 0 and 1. This reduces the term

in equation 5, making it appear that there is less drug present. This will produce a less than additive effect. This assumes that drugs A and B have the same maximal effect. It is possible for some approaches to synergy analysis to show apparent synergy if the maximal effects of drugs A and B are not identical, even if U⁵⁰(θ) = 1 for all values of θ.

The model for the steepness term, γ(θ), can similarly be described from equation 6, with appropriate substitutions for β⁰and β¹. The resulting equation is

Equations 6–10describe straight lines (simple additivity) when the coefficients (i.e. , β², β³, β⁴) are 0. They are the equations for parabolas if the respective β²coefficient is nonzero, and β³and β⁴are 0. More complex shapes are generated when β³and β⁴are nonzero.

Figure 5shows E^max, U⁵⁰, and γ as functions of θ for the synergistic interaction seen in figure 4. E,^maxand γ are constant, and thus have no interaction. U⁵⁰is necessarily 1 at θ= 0 and θ= 1, but is less than one between these extremes. This increases the potency of the drugs when administered in combination, resulting in the synergy seen in figure 4.

The model can be readily expanded to show the interaction of more than two drugs. In the case of three drugs (A, B, and C) the proportion of each drug present can be expressed by θ^A, θ^B, and θ^C, where

We can define the ratio of three drugs from just two of these ratios because θ^A+θ^B+θ^C=1. For our purposes here, we will use θ^Band θ^C. We again assume that for any fixed value of θ^Band θ^C, there is a sigmoidal relation between concentration and response. Therefore, if the three drugs could be administered to the effect site in an exactly fixed proportion, they would show a sigmoidal total concentration–response relation, where the “concentration” was the sum of the three normalized concentrations. This is precisely the notion that underlies the two-drug model. The equation for the three-drug model follows:

In the three-drug model the parameters of the sigmoidal relation, E^max, γ, and U⁵⁰, are functions of θ^Band θ^C. The functions E^max(θ^B,θ^C), U⁵⁰(θ^B,θ^C), and γ(θ^B,θ^C) are described in Appendix 2 (which can be found on the Anesthesiology Web site at www.anesthesiology.com). The important point is that instead of describing a curve, as in equations 6–10, when three drugs are present, the parameters of the sigmoidal relation (equation 12) are themselves surfaces. Figure 6shows representative surfaces for U⁵⁰as functions of θ^Band θ^C.

The response-surface model was written in the C programming language by the authors and compiled as a dynamic link library for Excel (Excel 7.0 for Windows; Microsoft Corporation, Redmond, CA).

We also implemented the model for NONMEM (NONMEM Project Group, University of California at San Francisco, CA). 3This code is available on request to the authors. **Graphs were prepared using Excel and Mathematica (Mathematica for Windows, Version 4.0; Wolfram Research, Champaign, IL).

To demonstrate the use of our response-surface model, we analyzed data previously published by Short et al.  4These data are also available via  the Anesthesiology Web site (www.anesthesiology.org). Dose–response relations were established after intravenous bolus doses of midazolam, propofol, and alfentanil administered individually and in combination in 400 women patients undergoing elective gynecologic surgery. Patients were assessed for hypnosis (defined as failure to open the eyes to verbal command) 4 min after midazolam administration and 2 min after propofol or alfentanil injection (the approximate times to peak effect after an intravenous bolus). When the combination being tested included midazolam, it was administered 2 min before the other drugs. The doses of midazolam, propofol, and alfentanil used and the proportion of patients achieving hypnosis for each dose category are shown in table 2. This data set was selected because it provided three two-drug combinations that could be used to demonstrate response–surface relations. In addition, the data are notable for the relatively large number of subjects and the uniformity of the experimental conditions. Lastly, it provides the ability to include a three-drug interaction model, showing the flexibility of the proposed response-surface model.

We assumed in the model that all subjects were awake when no drug was present; therefore, E⁰was 0 by definition. In addition, we assumed that each drug was capable of causing complete hypnosis if administered in a high enough concentration.

These assumptions permitted a simplification of equation 12. E⁰(no hypnosis) and E^max(θ^B, θ^C) (complete hypnosis) can be constrained to 0 and 1 respectively, with no interaction across the E^max(θ^B, θ^C) parameter surface. Thus, the probability of hypnosis for any combination is

where U^A, U^B, and U^Care the doses of midazolam, propofol, and alfentanil, respectively. The units are fractions of each drug’s dose necessary to cause hypnosis in 50% of the population (D⁵⁰). Four separate analyses based on equation 13were performed: (1) the data for the single and paired combination of midazolam and propofol were modeled; (2) the data for the single and paired combination of midazolam and alfentanil were modeled; (3) the data for the single and paired combination of propofol and alfentanil were modeled; and (4) the complete data set for the single, paired, and three-drug combinations (table 2) were modeled simultaneously. Model parameters and standard errors were estimated using NONMEM, by maximizing the log likelihood (LL) for all 400 observations:

where Rⁱis the observed response of the ith individual, either 0 (respond to voice) or 1 (not respond to voice), and P  is the probability of response to voice for each dose combination. Equation 14can be expressed in words as the sum of the natural logarithm of the probabilities of “no response in the nonresponders” and “response in the responders.” The contribution of the polynomial coefficients to the model was evaluated by excluding the coefficients one at a time, by determining whether the model deteriorated significantly (likelihood ratio test), and by assessment of residuals (observed vs.  predicted probability of hypnosis for each dose combination). Parameter estimates and standard errors were tabulated. The response surfaces for the paired interactions and the U⁵⁰(θ^B,θ^C) parameter surface were graphed.

Excel Solver was used to find the maximally synergistic combinations of the paired and the three-drug interactions based on the parameter estimates obtained from the simultaneous analysis of the entire data set. The intravenous bolus doses necessary to achieve 95% probability of hypnosis in this population were tabulated for each drug alone, for each paired combination, and for the three-drug combination. To illustrate the application of the response-surface model with concentrations, the pharmacokinetic parameters of Bührer et al. , 5Schnider et al. , 6,7and Scott and Stanski 8were used to calculate the C⁵⁰for hypnosis for midazolam, propofol, and alfentanil, respectively. This was calculated as the predicted effect-site concentration at the time of assessment after the respective D⁵⁰bolus. Equation 13was then used to calculate the predicted time course of effect (probability of hypnosis in the population) after these maximally synergistic doses of midazolam, propofol, and alfentanil, administered alone and in combination. The pharmacokinetic parameters, the calculated C⁵⁰values and the time from hypnosis in 95% of the population to “no hypnosis” in 95% of the population were tabulated. The predicted time course of effect was graphed. All simulations were performed using PKPD Tools for Excel (PKPD Tools for Excel; written by C. Minto and T. Schnider, Stanford University, Palo Alto, CA).

The shapes of the response surfaces generated by the aforementioned equations are not readily visualized. We therefore used three-dimensional graphics to illustrate the response surface for a variety of interactions between two drugs, by varying the model parameters of equation 5. These interactions include additive, synergistic, and antagonistic interactions between two agonists and interactions between full agonist–partial agonists, agonist–competitive antagonists, and agonist–inverse agonists.

We were able to include the data for all 400 patients in the analysis, in contrast to the method used by Short et al. , 4which necessitated the exclusion of the data for 120 patients. Parameter estimates and standard errors for the analysis of the three paired drug interactions are shown in table 3. Parameter estimates and standard errors for the simultaneous analysis of the entire data set are shown in table 4. There were no significant differences in the slope parameter of the three drugs, nor were there significant drug interactions across the γ(θ^B, θ^C) surface. The D⁵⁰values for each drug tended to be the same, regardless of the combination being modeled (tables 3 and 4). This is an expected result. Because the data for single administration was included in the interaction models, the D⁵⁰values were almost entirely determined from subjects receiving each drug alone. Figure 7shows the response surfaces for each of the paired interactions. All paired drug interactions showed significant synergy. The triple synergy parameter in the three-drug model was not statistically significant (see Appendix 2 on the Anesthesiology Web site at www.anesthesiology.com). This indicates that, when all three drugs are present, there is not synergy beyond that expected from the paired interactions of all three drugs.

The U⁵⁰(θ^B, θ^C) surface is shown in figure 8. The maximum decrease in U⁵⁰(θ^B, θ^C) values for the paired combinations are represented by the midpoint of the three edges of the triangular surface. Values follow (expressed as a percentage decrease): midazolam–propofol = 35%; midazolam–alfentanil = 44%; propofol–alfentanil = 16%. The maximum decrease (45%) in U⁵⁰(θ^B, θ^C) for the three-drug combination is found at the lowest point of the U⁵⁰(θ^B, θ^C) surface, which is located at θ^B= 0.144 and θ^C= 0.396. This corresponds in figure 8to θ^Mid–Prop= 0.24, θ^Mid–Alf= 0.47, and θ^Alf–Prop= 0.73, where the point appears as a dot on the axis plane.

Table 5shows the doses for the single, paired, and three-drug combinations optimized for maximum synergy associated with 95% probability of hypnosis (in a population of 50-kg women) based on the parameters shown in figure 6. The time between 95% probability of hypnosis to 95% probability of no hypnosis for each combination was calculated based on the simulations shown in figure 9. Although the maximally synergistic three-drug combination enables administration of one tenth the propofol dose, this results in a threefold increase in the time for 95% to awaken. Propofol alone is the drug of choice when the desired end point is a brief period of hypnosis. A knowledge of the occurrence and time course of other effects, such as anxiolysis, amnesia, analgesia, bradycardia, hypotension, muscular rigidity, and respiratory depression is necessary to determine the optimum drug combinations for other clinical end points.

TABLE 6 

Figure 10shows our two-drug, response-surface model (equation 5) for an additive interaction (A), a synergistic interaction (B), and an antagonist interaction (C). It also shows the interaction between a full agonist and a partial agonist (D), a full agonist and a competitive antagonist (E), and a full agonist and an inverse agonist (F).

One way to describe our model for two drugs, A and B, is that along the drug A axis and the drug B axis there are two one-drug sigmoid E^maxequations (equation 1). Our proposed model connects these two sigmoid curves via  polynomial functions of θ (equation 4). Drug interactions are then sought as coefficients of the polynomials that connect the “one-drug” sigmoids to each other. Although each value of θ can have its own E^max, C⁵⁰, and γ, the model assumes that the underlying sigmoidal shape is preserved for all values of θ. This concept of a fixed ratio (θ) of two drugs having their own properties is not new. Carter et al.  9considered each drug ratio to be essentially similar to a single agent, whose response could be described by a single two-dimensional concentration–response curve. Short et al.  4proposed that when using combinations of sedatives “the combination should be regarded as a new ‘drug’ with individual properties, rather than merely reflecting the known properties of the individual agents.”

The use of functions (equations 7–10) to connect the parameters of drug A to drug B assumes that the response surface is smooth and continuous. Although the basic model (equation 5) contains the usual parameters estimated in sigmoid E^maxmodels, the polynomials that underlie U⁵⁰(θ), E^max(θ), and γ(θ) are more complex. The use of polynomial response functions to approximate complex response surfaces is common in many experimental situations. However, polynomial models with the independent variable present at powers higher than the second power are not often used because it becomes difficult to interpret the coefficients. 2We proposed a very flexible model. Potentially, poor trial design, incomplete sampling of the response surface, or careless application of the curve-fitting procedures could result in estimation of parameters that provide a “good fit,” but when evaluated more closely generate an unrealistic shape for the surface. These concerns are not unique to this model and can be alleviated by proper model-selection techniques.

When evaluating the polynomial functions, one is particularly interested in whether the higher order terms can be dropped from the model. The statistical significance of the higher order terms should be tested against the simpler models using the likelihood ratio test. In addition to standard methods of model development, such as evaluation of standard errors and the pattern of residuals, the model should also be evaluated by graphing the entire response surface, and by graphing the individual model parameters as functions of θ (e.g. , equations 7–10) to ensure that the pharmacodynamic parameters do not extrapolate to unreasonable values. For instance, it must be ensured that U⁵⁰(θ) and γ(θ) are always positive in the range of 0 ≤θ≤ 1. In the case of three drugs, the parameter surface should be graphed as shown in figure 8.

Berenbaum 1developed a mathematical proof showing that zero-interaction holds if and only if the isobolograms are straight lines. For this to be true, E^maxmust be constant with respect to θ. The case of γ(θ) is more complicated. Ifγ is constant with respect to θ, there is zero-interaction if E^max(θ) and U⁵⁰(θ) are also constant. But ifγ^Ais not equal to γ^B, there is no easy way to use our interaction model to test for zero-interaction (as defined by Berenbaum 1). Although we agree with Berenbaum’s 1treatment of the zero-interactive surface, the description of an interaction as zero-interactive, synergistic, or antagonistic may be too simplistic. For example, a drug combination can be synergistic in certain regions and antagonistic in others. 10To our thinking, the focus on whether drug interactions can be reduced to simple descriptors such as “additive,”“synergistic,” or “antagonistic” completely misses the point. The interaction has the potential to be highly dimensional and complex. Rather than worry about which descriptor applies to the relation, the goal should be to characterize the response surface. From the surface one can identify the best combination to produce the specified therapeutic effect.

Our model is empirical. It makes no assumptions about the mechanism of interaction between the drugs. However, we assume that the concentration–response relation for each of the interacting drugs is described by a direct pharmacodynamic model. We have not attempted to describe interactions between drugs that act directly and those that act indirectly, 11although we are not aware of any reason that our approach could not be combined with the more kinetic approach of indirect interactions.

The fundamental model for each drug does not have to be the sigmoid E^maxmodel (equation 1). For example, the model could be a linear, polynomial, or exponential response function. The model could also be a bimodal response function, as seen with some hypnotics. 6In fact, the model can be any function, so long as it has parameters that permit connection among the individual single-drug models. Thus, the only constraint is that the interaction model reduces to the model for drug A when θ= 0 and to the model for drug B when θ= 1.

Despite the flexibility of response-surface approaches, difficulties can still arise. For example, the occurrence of adverse effects might prevent the collection of data in a sufficient range to identify the model parameters for one of the drugs. An alternative to specifying a administered function for the concentration–response relation is to use a semiparametric function, as described by Troconiz et al.  12These authors discuss the use of flexible nonparametric functions (splines) that are forced to obey certain constraints (for example, a spline can be constrained to resemble an E^maxmodel). Although they illustrate their approach in the context of antagonistic interactions, it can be extended to model additive and synergistic drug interactions.

We used our response-surface model to analyze the drug interactions for the hypnotic end point between midazolam, propofol, and alfentanil. Based on this end point, a maximum effect of 100% for all three drugs described the data well. Although the high dose of propofol (alone) and alfentanil (alone) abolished the response to voice in 100% of patients, the high dose of midazolam (alone) abolished the response to voice in only 80% of patients (table 2). These results do not suggest that alfentanil will prevent purposeful movement or awareness in response to a surgical incision; nor do they prove that higher doses of midazolam will abolish response to voice in 100% of patients.

Our computer simulations showed the ability of our response-surface model to take into account potential interactions of two or three drugs at the effect site over time. We emphasize that these simulated effect-site concentrations are based on published information regarding pharmacokinetic parameters and do not take into account possible pharmacokinetic interactions between the three drugs. Notwithstanding this reservation, we were able to demonstrate that the maximally synergistic dose ratio is not necessarily advantageous if a brief hypnotic effect is the desired end point.

Although the three-dimensional view may aid in understanding and comparing different methods used to study drug interactions, it does not contain information that cannot be presented by a sequence of two-dimensional figures such as the isobologram (horizontal “slices” of the response surface, representing typical contour representations of three-dimensional graphs), the response to increasing concentration of one drug in the presence of fixed concentration of the other drug (vertical slices parallel to one axis), and the response to fixed concentration ratios of the two agents (radial slices). With any of these methods, a three-dimensional picture of the response surface can be constructed if sufficient slices are studied. In the case of three drugs, it is no longer possible to graph the response surface because it is a four-dimensional surface, although the model parameters themselves can be graphed in three dimensions (fig. 6). Development of techniques using parallel coordinate axes might prove to be useful aids in the visual interpretation of the interactions between three or more drugs. 13 

In conclusion, the study of drug interactions in anesthesia traditionally has used isobolographic analysis or multiple logistic regression. Both approaches have significant limitations. In particular, the multiple logistic regression approach is so beset with flaws (as described in Appendix 1 in the web supplement), including fundamental incompatibility with accepted pharmacodynamic principles, that it should be abandoned. The application of response-surface methodology to the study of drug interactions has the potential to overcome the limitations of these models. We proposed a flexible model for drug interactions, which describes the full relation between the concentrations of two or three drugs and drug effect. We illustrated our new model using previously published data and showed that this model can also describe differing types of interaction between an agonist, a partial agonist, a competitive antagonist, and an inverse agonist. Application of response-surface methodology permits characterization of the full concentration–response relation and therefore can be used to develop practical guidelines for optimal drug dosing.