In anesthesia, a synergistic effect of the combined treatments gives patients the following favorable outcomes: increasing the efficacy of the anesthetic effect, reducing doses while increasing or maintaining the same efficacy to avoid toxicity of anesthetic drugs, minimizing or slowing the development of drug resistance, and providing selective synergism against the target (or efficacy synergism) versus the host (or toxicity antagonism) [1]. Many such studies can be found in related literature [2-8].

There are 3 classes of drug interaction: additive interaction, in which the response of a drug combination is just what is expected from the dose-response relationships of drugs; synergistic interaction, in which the response is greater than expected; and antagonistic interaction, in which the response is less than expected. There is no uniform agreement on the terminology of drug interaction. Synonyms for what is termed additive interaction are zero, or null interaction, independence (Bliss), and indifference. Synonyms for synergistic interaction are positive or supra-additive interaction, potentiation, and augmentation. Synonyms for antagonistic interaction are negative, infra-additive or subadditive interaction, and negative synergy. These definitions of drug interaction suggest that a zero interaction model should exist to be used as a reference in classifying the interaction of drug combinations. A pharmacodynamic model of most anesthetic drugs is the sigmoid E_max model (Hill equation); let's discuss drug interactions and response surface focusing on this model.

Because the classification of drug interactions is determined as a greater or less pharmacological effect of a two-drug combination than what would be predicted for zero interaction from the knowledge of the effects of each drug individually, its determination absolutely depends on the reference model of zero interaction. A review by Berenbaum [9] has listed 8 approaches, and another review by Greco et al. [10] categorized 13 approaches and methods for determining the class of drug interaction. To define the class of drug interactions, first, an adequate reference model of additive (or zero) interaction should be established as a universal standard reference. Historically, 3 zero interaction models are commonly used [9,10]: Bliss independence [11], Loewe additivity [12], and the median-effect method from the law of mass action [13].

Bliss independence implies that 2 drugs do not pharmacologically or physiologically cooperate with each other for each drug behaves independently of the other [10]. The general equations of Bliss independence (Eq. (1)) and a specific one (Eq. (2)), which assumes that a sigmoid E_max relationship [14,15] is an appropriate pharmacodynamic model for each drug, are as follows:

In Eq. (1), fu₁, fu₂, and fu₁₂ are the fractions of response for drug 1, drug 2, and combination unaffected [16]. In Eq. (2), is a fraction of the effect resulting from the mixture of d₁ and d₂; d₁ and d₂ are doses of each drug in the mixture that yield the effect E; ED_{50,1 or 2} is a dose of drugs to produce 50% of the maximal effect (E_max) for each drug acting alone; and γ₁ and γ₂ are Hill coefficients (slope of dose-response curve). When Eq. (1) and (2) are recast in terms of the fraction of effect affected (fa_{1 or 2}) (if 1 - fa exchanges fu in Eq. (1), because fa + fu = 1), then Eq. (3) and (4) are the results.

Eq. (3) is similar to the equation of probabilities of independent events [17]. The curvature of the Bliss independence line depends on the Hill coefficient of each drug (Fig. 1).

Lines connecting points representing equi-effective dose combinations are termed isoboles (Fig. 1). They were introduced by Fraser [18,19], who claimed that, in using this device for evaluating the antagonism between the actions of physostigmine and atropine, "the results of the experiment will be rendered apparent by a mere glance at the diagram." Loewe extended this method to synergism [12,20]. Loewe used a straight line isobole for zero interaction when the combined drugs had similar mechanisms of action and similar dose-response relationship, and asserted that the isobole would be curved when combined drugs had different mechanisms of action and dissimilar dose-response relationship. Recently, Tallarida [21] supported Loewe's assertion. When a full agonist and a partial agonist have different maximum values, or when Hill coefficients (γ of sigmoid E_max model) of two full agonist drugs are different, each case has a varying potency ratio. For a full agonist and a partial agonist, the isobole of no interaction is no longer straight but curved [22]. For two full agonists with a varying potency, termed "heterodynamic" by Loewe [12], the application of dose equivalence leads to not just one, but to two nonlinear additive isoboles [21,23]. However, Berenbaum [9] maintained that a straight-line isobole is appropriate for a zero interactive combination, for it is not derived from the knowledge of the shapes of the dose-response curves or of their mechanisms of action, and he concluded that it is valid irrespective of the shapes of these curves, similar or dissimilar, and also, irrespective of their mechanisms of action. Although the shape of the isobole for Loewe additivity remains to be further studied, a straight line of additivity is commonly employed to distinguish synergistic and antagonistic from additive interactions. Although there are many reasons to prefer the Berenbaum method, the most predominate reason to use it may be because of easy calculations.

The general equation of the Loewe additivity (Eq. (5)) and a specific form (Eq. (6)), which assumes that a sigmoid E_max relationship is appropriate for each drug, is defined as

where d₁ and d₂ are doses of each drug in the mixture that yield an effect, E, while D_{1 or 2} is the dose of each drug to produce the same effect of E when given alone. The term is called the interaction index at the combination dose (d₁, d₂). The dose (D₁, or D₂) for each single drug producing the same effect of E is expressed as

If the interaction index at (d₁, d₂) is equal to, less than, or greater than 1, the combination dose is asserted to be additive, synergistic, or antagonistic, respectively.

The median-effect principle (Eq. (7)) was originated by Chou [24,25]. The median-effect equation describes dose-response relationship in the simplest terms, which is given by

where d is the dose of a drug; f_a is the fraction affected by d; f_u is the fraction unaffected by d (i.e., f_u = 1 - f_a, and f_a / f_u = odds ratio); D_m is the median effect dose (e.g. EC₅₀, or ED₅₀); and m is the slope coefficient of dose-response curve. When m is greater than 1, the dose-response curve becomes S-shaped, so m is equivalent to the Hill coefficient of sigmoid E_max equation.

Chou and Talalay [16] extended the median-effect equation for a single drug to multiple drugs. Thus, for a two-drug interaction, they got the same equation of combination index (CI) (Eq. (8)) as that of the interaction index by Berenbaum [26]. This equation is based on the assumption that two drugs bind the same receptors and are mutually exclusive. This gives

where d₁ and d₂ are doses of each drug in the mixture that yield an effect, f_a, while D_{1 or 2} is the dose of each drug to produce the same effect of f_a when given alone. D_{m,1 or 2} is the median effect dose of drug 1 or drug 2, and m_{1 or 2} is the slope coefficient of dose-response curve of drug 1 or drug 2. In the view of interaction, CI < 1, = 1, and > 1 indicate synergy, additivity, and antagonism, respectively.

If it is assumed that two drugs do not share receptors, do not interfere with each other, and are mutually nonexclusive, the resulting equation (Eq. (9)) should include a third term, the product of the first two terms, thus,

As shown in Fig. 1, the isobole of CI for the mutually nonexclusive model is concave upward, and is identical to that of the Bliss independence model (Eq. (4)) in which the Hill coefficient (γ) of sigmoid E_max model is 1. Therefore, this model under-estimates synergistic effects.

Because partially exclusive or partially nonexclusive cases may exist, and the equation with the third term may underestimate synergistic effects, Chou [1] claimed that the equation without the third term should be used as the "base equation" and that mutual non-exclusivity should be used as a contributing factor for the intrinsic synergistic effect under the assumption of mutual exclusivity. To be consistent with the classic isobologram and its equation (interaction index of Berenbaum), he decided that mutually exclusive assumption will be used for the analysis of drug interaction.

The Bliss independence and Loewe additivity models are the two most cited reference models for determining drug interaction [9,10,27,28]. According to Berenbaum, only when the dose-response relationship of each drug follows an exponential decline with the dose, the two models result in the same equation and are consistent with each other [10]. When dose-response curves are steeper than the exponential model, Bliss independence will result in Loewe antagonism; whereas, when dose-response curves are less steep than the exponential model, Bliss independence will result in Loewe synergism (Fig. 1). Greco et al. [10] revealed that they prefer Loewe additivity to Bliss independence for the lack of synergism or antagonism. When curves are steep (γ > 2), the Bliss independence criterion may overestimate synergism, whereas the Loewe additivity model can overemphasize antagonism. That is to say, a class of drug interaction is dependent on a reference model because a given effect seems to be either synergistic with the Loewe additivity or antagonistic with the Bliss independence. However, the use of both models is not practical because one of the two models could be more plausible than the other [27]. Goldoni and Johansson [27] said that although the Loewe additivity model is slightly more plausible and is preferred in toxicology, it could be incorrect under certain conditions.

The interpretation of an assessment of drug interactions by the Loewe additivity reference model is free of mechanistic restrictions and implications. It is possible to determine that the combination of two or more inhibitors is more effective than their individual use by means of isobolographic analysis, even when no information about their mechanism of action is available [29]. Accordingly, the Loewe additivity reference model is consistent with the graphical isobologram approach [29]. The sigmoidicity of a dose-effect curve greatly magnifies the differences among the different methods or theories (Fig. 1). Thus, the main controversies in drug combination analysis in the past century can be readily resolved [1]. Drug additivity is substantiated under the Loewe additivity model but not the Bliss independence model [30]. Therefore, the Loewe additivity model is a universal reference model for classifying drug interaction.

There have been several methods for evaluating drug inter-action: diagonal array [17], isoboles [9,17,31], isobolar surface for three-drug combination [17], interaction index [9,17], and response surface [10]. Isobolographic analysis and interaction index method are not universally applicable, and they rely on the level of effect; one isobole relates only to a single effect, and one interaction index does only to a single dose pair. Therefore, separate analyses should be conducted at numerous levels (ED₅, ED₂₅, ED₅₀, ED₇₅, etc.) in order to classify all possible interactions. To depict the entire set of levels of effect (from ED₁ to ED₉₉), one may build all the isoboles in three dimensions, thus, creating a surface of isoboles, the so-called response surface. It illustrates the drug effect (Z axis) versus two-drug doses (X and Y axes) [10] and presents an entire drug interaction at all dose pairs.

Over the last decade, drug interaction studies have been increasingly evaluated on the basis of a response surface model. The approaches to choosing a response surface model are empirical or functional [32]. Many examples of response surface methods employed a single interaction parameter to catch synergy, additivity, or antagonism [33-37]. These models are sound if only synergy, additivity, or antagonism exists throughout the entire surface; they are incompetent to describe the interspersion of regional synergy or regional antagonism in different areas of drug combinations [38].

Isoboles are not necessarily consistent or similar because they may cross the additivity line so that some combinations with a specific effect are synergistic, and others antagonistic. Berenbaum maintained that the conclusion as to whether a combination has a specific class of drug interaction applies to that combination but not to untested combinations [9]. A model with a single interaction parameter never reflects an inconsistency in the patterns of drug interactions. To solve this problem, some researchers proposed a few parametric models with a polynomial interaction function for two drugs [38-40]. The interaction functions capture synergy, additivity, and antagonism scattering all over the surface of drug combinations.

The analyses of the following data were performed using the 7 methods reviewed here. The aim of the study was to investigate the combinations of vecuronium and rocuronium to have synergistic interaction. Left phrenic nerve-hemidiaphragms of 48 male Sprague-Dawley rats (150-250 g in weight) were mounted in Krebs solution. Supramaximal electrical stimulation (0.2 ms, rectangular) of 50 Hz for 1.9 s was applied to the phrenic nerve, and the evoked tetanic contractions were measured with a force transducer. Each preparation was exposed to one of four vecuronium concentrations (0.0, 1.5, 2.5, 3.0 µM), or one of four rocuronium concentrations (0.0, 3.0, 4.5, 5.5 µM). Subsequently the adequate amount of rocuronium to a vecuronium bath or vecuronium to a rocuronium bath was cumulatively added until an 80-90% increase of tetanic fade was reached.

In this data set, E₀'s of two drugs are 0, and E_max's are 1 (or 100%). For interaction parameter and interaction index of the Greco, and Machado models, the values of Hill slope and potency of each drug acting alone were calculated (vecuronium: 6.5125, 2.98 µM, rocuronium: 7.8975, 5.34 µM, respectively) and then substituted for ED₅₀'s and γ's of the two models. Interaction indices of the Berenbaum method (Eq. (5)) were directly calculated using the sigmoid E_max equation substituted with the above calculated values. Plummer model (Eq. (14)) and Carter Model (Eq. (17)) were fitted with the data by using Tablecurve3D®.

The iso-effective contours (Fig. 2) and the plot of interaction indices of actual data points (Fig. 3) were drawn. All the isoboles of the 4 models are upward concave and consistently indicate synergistic interaction (Fig. 2). The interaction index of the Machado model cannot be calculated due to the limitation of the model. Table 1 shows the results of five interaction analysis methods. Four models, except for the Berenbaum method, consistently show interaction indices of less than 1, which mean synergistic interaction (Table 1 and Fig. 3). In Fig. 3, the Berenbaum method reveals that the interaction indices are less than 1 at the high effect region and are greater than 1 at the low effect region. According to the Berenbaum method, combinations of the two drugs have all types of interaction. Therefore, the response surface models with a single interaction parameter are inadequate to analyze the data which has complicated scatter of different types of drug interactions (Fig. 3 and 4, Panel A).

Next, three aforementioned flexible response surface models with an interaction function were fitted to the same data using Tablecurve3D®, and the parameters, the response surfaces, isoboles, and contours of the interaction index of each model were determined, and plots were made for seeking the differences of their results (Fig. 4 and 5). The functions for the parameters U₅₀(θ) and γ(θ) of the Minto model were calculated and illustrated (Fig. 6). The full model (10 parameters) for the Minto model was analyzed; the Fidler model (9 parameters) was analyzed when w = 1; and the Kong model (8 parameters) was analyzed when the relative potency ratio was constant. The isobole plot of raw data is rendered via Renka I (nonparametric algorithm) procedure, and it is shown in Fig. 4, panel A. Pharmacodynamic parameters of response surface models are listed at Table 2. The values are different depending on the model. In the Kong and Plummer models, because relative potency ratio is constant, the Hill slopes of each model are identical.

The isobole plots of the Minto model, Fidler model, and Kong model were relatively similar to the isobole plot of raw data. Minto et al. determined the pattern of drug interaction based on the functions of E_max, U₅₀(θ), and γ(θ), especially the value of U₅₀(θ) (Fig. 6). The isobole plot of the Minto model is shown in Fig. 4, panel B. In this article, with the equation of the interaction index, the contour plot of the interaction index (Fig. 5, panel A) can be illustrated and help to make the quantitative evaluation. Chou and Hayball [45] proposed that an interaction index from 0.9 to 1.1 is designated as being additive. In the Minto model, the lower extreme effect region has greater synergistic interaction. Chou [13] said that for anticancer or antiviral agents, synergy at high effect levels (fa > 0.8) is more relevant to therapy than at low effect levels (fa < 0.2). On the contrary, in the anesthetic area, synergy at low effect levels may be more crucial than at high levels because the residual effects of anesthetic agents could threaten the safety of patients.

In the Fidler model, the pattern of drug interaction is determined based on contour of iso-effective combinations, contour of interaction index, a direct transformation of the α term (% contribution of α on a 50% effect), and the best ratio for synergistic effect. The isobologram (Fig. 4, panel C) shows the deepest upward concavity around the combination of the unit ratio of 0.78 : 0.22, vecuronium : rocuronium (Table 2). The interaction indices range from 1 to 0.78. According to Chou and Hayball [45], the contour plot of interaction index (Fig. 5, panel B) reveals synergistic interaction, except for the area in which the combination of two drugs has a higher concentration of rocuronium than of vecuronium. A direct transformation of the α term (% contribution of α on a 50% effect) is 5% in this data set. In other words, the interaction of vecuronium and rocuronium requires 95% of the additive drug combination to produce the 50% effect at the most effect ratio. The minimum interaction index (0.78) is in about a 99% effective zone. The way in which to make this discrepancy has not been determined within the course of this research.

In the Kong model, it is assumed that relative potency is constant (that is, two drug dose-effect curves are parallel) like the analysis of the Plummer model. The shape of isobole in this model is most similar to that in the raw data among seven models (Fig. 4, panel A and D). The pattern of drug interaction is determined based on isobole plot, contour of interaction index, and polynomial function (not seen here). According to the Chou and Hayball [45], the isobologram (Fig. 4, panel D) shows synergistic, additive, and antagonistic interactions, and the interaction indices range from more than 1 to less than 0.6 (Fig. 5, panel C). The extreme effect regions get greater synergistic interaction.

The distribution of patterns of drug interaction depends on a response surface model. Which model best approximates the information in data should be decided using the model selection.

The information scientists want to reveal is in the observed data. It has been expressed as a model. Making valid inferences from scientific data depends on a model of the information in the data [46]. For selecting the best model, it is assumed that there are data and a priori set of models and that statistical inference is to be model-based.

Researchers should conceive many models (i.e., multiple working hypotheses) and at the same time, cope with how to choose among the possibilities. The approaches to model selection use classical hypothesis testing [47], the Akaike information criterion (AIC) [48], the Bayesian information criterion (BIC) [49], and cross validation, etc. [46]. Generally, hypothesis testing is a very poor basis for model selection [48], especially, tests between models that are not nested [46]. In cross validation, the data are divided into two partitions, the first partition is used for model fitting; and the second one is used for model validation. A new partition is continuously selected, and this whole process is repeated hundreds or thousands of times. This method may be an alternative method for model selection but is time-consuming.

The target models of AIC and BIC are different [50]. The target model of AIC is specific to the sample size; it is the best model to minimize the estimated Kullback-Leibler information loss [51,52] when models are used to approximate the full reality (or truth). This target model depends on the change of sample size. Even the priori set of models may change when the sample size changes greatly. In contrast with AIC, the derivative assumption of BIC is that there is a true model in a model set. A true model is one that generates the data, independent of data size. The target model of BIC is the true model. The goal of model selection is to obtain a perfect model so that no information is lost in fitting data to a model of the information in data. It is impossible to achieve this ideal purpose because models are only approximations. However, AIC helps us choose the best model, which loses as little information as possible and allows the efficient and objective filtration of information from noise [46].

The principle of parsimony is that model selection is a bias versus variance trade-off. A parsimonious model must be as simple as possible in terms of variables and model structure. Inference from a model with too few variables may be biased, and inference from a model with too many variables may not be precise and can result in spurious conclusions. Parsimony lies between underfitting and overfitting [46]. For linear models, such as multiple regression, a minimum of 10 to 15 observations (or events) per predictor variable will generally allow good estimates [53-56]. Many researchers do not notice the problem of overfitting in regression-type modeling. Overfitting results in overly-optimistic findings that do not present in the real population and will not reproduce in a future data set. Although overfitting is not a bad beginning to seek a good model, spurious conclusions generated by overfitted model may lead doctors the wrong way in clinical decision-making.

In this article, 7 models are introduced, and the last 3 models have several nested submodels. A model should be selected to best fit a dose-effect data set using AIC, not hypothesis testing.

In anesthetic practice, anesthesiologists administer multiple drugs to patients simultaneously. They have been long aware of the importance of drug interaction and have performed many investigations for the efficient combination of anesthetic agents. The study of drug interaction is aiming at choosing the optimal drug combination.

Determining what pattern of drug interaction is acting is dependent on a reference model [10]. Although further study should be continued, Loewe additivity is a universal reference model for the classification of drug interaction. The reference model has been established, and then, many analytic models for drug interaction were proposed to be validated for their performances.

A lot of models have prevailed. In the anesthetic area, models with a single interaction parameter or flexible models with an interaction function are adequate to be applied to analyzing the data of anesthetic drug interaction. Due to detecting the best approximating model among the aforementioned models, AIC is the most popular approach for choosing the best model.

However good a chosen model may be, it never finds out the truth which generates the data set, and, of course, it is uncertain whether a chosen model is the best model for approximating information in the data. Instead of making an inference from only a chosen model estimated as the best, a more robust inference can be extracted from averaging several models that are considered relevant. This procedure is called multimodel inference [57], and this concept needs further study.