Theory of synergistic effects: Hill-type response surfaces as ‘null-interaction’ models for mixtures

Background The classification of effects caused by mixtures of agents as synergistic, antagonistic or additive depends critically on the reference model of ‘null interaction’. Two main approaches are currently in use, the Additive Dose (ADM) or concentration addition (CA) and the Multiplicative Survival (MSM) or independent action (IA) models. We compare several response surface models to a newly developed Hill response surface, obtained by solving a logistic partial differential equation (PDE). Assuming that a mixture of chemicals with individual Hill-type dose-response curves can be described by an n-dimensional logistic function, Hill’s differential equation for pure agents is replaced by a PDE for mixtures whose solution provides Hill surfaces as ’null-interaction’ models and relies neither on Bliss independence or Loewe additivity nor uses Chou’s unified general theory. Methods An n-dimensional logistic PDE decribing the Hill-type response of n-component mixtures is solved. Appropriate boundary conditions ensure the correct asymptotic behaviour. Mathematica 11 (Wolfram, Mathematica Version 11.0, 2016) is used for the mathematics and graphics presented in this article. Results The Hill response surface ansatz can be applied to mixtures of compounds with arbitrary Hill parameters. Restrictions which are required when deriving analytical expressions for response surfaces from other principles, are unnecessary. Many approaches based on Loewe additivity turn out be special cases of the Hill approach whose increased flexibility permits a better description of ‘null-effect’ responses. Missing sham-compliance of Bliss IA, known as Colby’s model in agrochemistry, leads to incompatibility with the Hill surface ansatz. Examples of binary and ternary mixtures illustrate the differences between the approaches. For Hill-slopes close to one and doses below the half-maximum effect doses MSM (Colby, Bliss, Finney, Abbott) predicts synergistic effects where the Hill model indicates ‘null-interaction’. These differences increase considerably with increasing steepness of the individual dose-response curves. Conclusion The Hill response surface ansatz contains the Loewe additivity concept as a special case and is incompatible with Bliss independent action. Hence, when synergistic effects are claimed, those dose combinations deserve special attention where the differences between independent action approaches and Hill estimations are large.


Background
Quantifying the effect of an active substance upon application to its target organism is a central topic in life sciences. The description and prediction of dose-response curves of pure compounds and estimating the combined effect of simultaneous administration of several active ingredients (a.i.s) is a field of active research in pharmacology, anesthesiology, toxicology, environmental science, and agrochemistry. The detection of virtually any kind of interaction between drugs is, e.g., of utmost importance for risk assessment in the registration process of chemicals in general. While synergistic adverse effects present a challenge to say the least, synergistic effects in drug administration are desirable as they might be used to reach the same goal with less effort. Disregarding the steadily growing area of possible interactions of chemicals with their environment, this article concentrates on interactions between a.i.s in the life-sciences. The terms 'agent' , 'drug' , 'chemical' , and 'a.i. ' are used interchangeably.
Generally, biological effects of an a.i. follow a dose-response relation, starting from a no-effect level (NOEL) and ending at the maximum effect corresponding to a saturation dose. In order to quantitatively describe interactions between substances, it is necessary to define a reference behaviour, namely the response of a system of compounds acting independently. Deviations from this ideal reference can then be classified as synergistic or antagonistic, being aware of the fact that these deviations can be caused by a multitude of physico-chemical and biochemical effects, especially in whole organisms (for a recent study see e.g. [9]).
Whenever the effect observed after applying a mixture exceeds the expectation, the action of the agents is called synergistic, and if it is smaller than expected it is antagonistic. This seemingly simple definition is not simple at all, as there is an ongoing debate since the beginning of the last century on how to correctly define this reference of 'no interaction' . In the literature, synergy is defined either phenomenologically or based on assumptions on the modes of action of the a.i.s involved. The observed effects like zerointeraction, synergism or antagonism, are often quantified by calculating interaction- [10] or combination-indices [11], however, often based on differing definitions of additivity [2,4].
Mainly two types of reference models are currently in use [12,13]. One class of models, the Additive Dose Model (ADM) pioneered by Loewe [14], assumes an additive behaviour of doses, similar modes of action and consequently parallel dose-response curves with identical slopes. Its basic concept is that of equipotent doses, which means that an arbitrary dose of a.i. A can be replaced by an isoeffective dose of a.i. B. The terms 'effect summation' , 'dose addition' , and 'concentration addition (CA)' are used in this context. As the joint action of the a.i.s is not independent in this model, the term 'mutually exclusive' [11] is more appropriate here. Implicitly the model assumes that the maximum possible effect achievable by the mixture is that of its most potent partner. A generalized concentration addition (GCA) model [6] extended the original ADM approach [10,[14][15][16] to mixtures with partially overlapping agonists.
The other class of models, called Multiplicative Survival Models (MSM), assumes that the effects caused by the two a.i.s are mutually non-exclusive and originate from different modes of action, and hence the asymptotically achievable effect is the sum of the individually possible ones. In agrochemical research, it is associated with the names of Colby [17] and Limpel [18], while in toxicology these models of independent action (IA) have been described by Bliss [19] and Finney [20] and in entomology by Abbott [21].
Chou's widely used combination indices [11] cover both classes, his formula for the mutually exclusive combination is of ADM type while his description of mutually nonexclusive mixtures is MSM like and reduces to Colby's formula for reaction orders of one. The approach presented here relies neither on Bliss independence or Loewe additivity nor uses Chou's unified general theory. It focuses on logistic (Hill) response surfaces as 'null-effect' models. For mixtures of n partners they result from solving an n-dimensional partial differential equation (PDE). Appropriate boundary conditions guarantee that the expressions describing the Hill-surfaces are asymptotically correct, i.e., lead to Hill's wellknown dose-response curves in the one-dimensional case of pure a.i.s. In addition, the solutions are sham-compliant, meaning that a mixture consisting of combinations of a drug with itself shows no synergistic effect.
Sigmoid dose-response curves like Hill's equation [22,23] are solutions of a class of ordinary differential equations (ODEs) which was originally used to describe population dynamics [24]. Examples of the phenomena described by this type of functions are titration curves in chemistry, dose-response curves in enzyme kinetics, and predator-prey models in biology.
Following a short recapitulation of the ODE leading to Hill-type dose-response curves for pure compounds, Hill response surfaces for binary and n-component mixtures are introduced as solutions of the corresponding PDEs. In the subsequent sections the properties of Hill surfaces are compared to other response surface models and applied to some literature examples, namely binary and ternary mixtures from various fields of life-science. A summary of our findings will be given in the last section.

Methods
While Hill's equation can be obtained by solving a first-order ODE, the logistic differential equation, its n-dimensional generalization results from solving a semilinear PDE with the appropriate boundary conditions. These are the requirements that in the limit of one dimension the original Hill equation results, and that the solution of the PDE is shamcompliant, meaning that an artificial partitioning of the one-dimensional problem into an n-dimensional one does not change the results.
Here an analytical expression of the n-dimensional extension of Hill's equation is provided. Verification is done by substituting the solution into the PDE by using Mathematica 11 [25]. Checking the fulfillment of the boundary conditions is achieved by performing the corresponding limiting processes.

Logistic functions and the Hill response surface
Let A and B be active ingredients with dose-response curves a(x) and b(y) depending on the variables x and y and let U be a combination of A and B with the dose-response surface u(x, y). Let us further assume that we can describe the effect z (being a(x), b(y) or u(x, y)), by a logistic function. Then z is characterized by 4 parameters, the minimum and maximum responses z min and z max , the position of the inflection point and a slope parameter, i.e., the slope at the inflection point of the curve. All effects are limited by the no-effect-and the full-effect levels, 0 ≤ z min ≤ z ≤ z max ≤ 1.
The differential equations for the one-and two-dimensional cases, describing the variation of the effects a(x), b(y) and u(x, y) with the variation of x and y, are where u x = ∂u(x, y)/∂x and u y = ∂u(x, y)/∂y denote the partial derivatives of u, and α, β and a max , b max are constants. In Eq. 2, γ and u max are functions of x and y.
Denoting by x 50 and y 50 the positions of the half-maximum effects a(x 50 ) = a max /2 and b(y 50 ) = b max /2 the final forms of a(x) and b(y) are the logistic functions with x = x − x 50 , y = y − y 50 and α and β being the slopes at the inflection points x 50 and y 50 . They show an exponentially increasing effect at low doses, becoming linear close to x = x 50 and y = y 50 and finally an exponentially decreasing growth until the limiting effect z max is reached. a(x) and b(y) are intimately connected to Hill's equation, providing a relation between the effect E and the dose C.
Actually, x and y are the natural logarithms of doses with −∞ ≤ x, y ≤ ∞, whereas the doses themselves (i.e., e x and e y ) are ≥ 0. Hence, we can identify the effects E 0 , E max and the shape parameter α with z min , z max and α or β, and the doses C and EC 50 with e x and e x 50 or e y and e y 50 of Eq. 3, meaning that Hill's equation is the solution of a logistic ODE, subject to the appropriate boundary conditions. Our approach to handle mixtures is completely analogous to that used for the pure compounds. The solution of the PDE describing the effect u(x, y) of a binary mixture (Eq. 2) is the logistic Hill-surface as the response surface of 'null-interaction' .
The boundary condition that in the limit of one vanishing agent the response due to the second one has to result, and the so-called sham-compliance requirement have to be fulfilled: If a dose d of one a.i. is artificially split into two contributions nd and (1 − n)d with 0 ≤ n ≤ 1, the response due to this 'mixture' has to be identical to the response of the pure compound, i.e., a(d) = u(n × d, (1 − n) × d), irrespective of any possible interaction between different a.i.s. In our formalism a dose d is equal to e x . Hence, sham partitioning means setting e x = ne x + (1 − n)e x , leading to the requirement u(x + ln n, x + ln(1 − n)) = a(x).
These conditions are satisfied by the ansatz γ (x, y) = αe x + βe y e x + e y u max (x, y) = a max e x + b max e y e x + e y (6) which provides smooth transitions between α and β and between a max and b max .
To facilitate a comparison with literature expressions for response surfaces, u Hill can be re-written by substituting the doses d a and d b for e x and e y and using the doses scaled by their median effects m a = d a /d a 50 = e x and m b = d b /d b 50 = e y . Then the 'nullinteraction' reference response surface is with As the ansatz presented here uses the known dose-response curves of the pure mixture partners to unambiguously define the 'null-interaction' surface for the binary mixture, there is no obvious way to introduce interaction parameters (and functions) which account for deviations from this reference [26][27][28]. At present the only way to detect synergism or antagonism is the direct comparison of experimental data with the Hill response surface, either visually or by statistical means like the root-mean squared error of prediction (RMSE).
A way out of this dilemma might be to obtain the necessary parameters directly from fitting the experimental mixture data to a logistic surface which could then be compared to the 'null-interaction' model. However, the functional form of u Hill is not flexible enough for this purpose and some sort of perturbation theory might be more appropriate to handle deviations from the Hill surface.

The logistic PDE for mixtures of n components
The extension of the Hill formalism to mixtures of n agents A i is straightforward. The corresponding semilinear logistic PDE [29] is with describes an n-dimensional logistic surface. It satisfies the boundary conditions, and the sham compliance condition, that any partition of the dose of a pure drug into various artificial new drugs must cause the same effect as the pure drug alone.
As an example, u Hill for a ternary mixture of A, B and C satisfying the sham condition where in a self-explaining notation

Discussion
From a theoretical point of view it is an advantage of the Hill surface approach that it does not rest on assumptions on maximum effects or restrictions on specific parameter combinations of the mixture components. This distinguishes our approach from other response surface models as will be outlined in the next section. In practical applications, however, these differences become significant only if steep dose-response curves and/or strongly differing maximum-effect parameters of the individual agents are involved.

A Comparison of empirical response surface models
Many of the response surface models are based on two synergy approaches, the Bliss independence [19,20] and Loewe additivity [14] models. The characteristics of some of them are listed in Table 1. There their functional forms describing the mixture effects are classified by their assumptions on maximum effects and slopes and by analyzing whether they are asymptotically correct and sham compliant. Apparently only a few authors [26,28,30] use the freedom to define the shape function γ ( x) and the maximum effect function u max ( x) (Eq. 6). Bliss independent joint action (IA) [3] can be formulated in terms of fractions of possible response unaffected fu or fractions affected fa, with fa + fu = 1. For a binary mixture Fu ab = fu a × fu b and Fa ab = 1 − Fu ab . Hence It is related to the expression for the probability P(A ∪ B) of an event A or B if the basic events with probabilities P(A) and P(B) are independent. In crop science it is known as MSM [13] or Colby [17] model and is widely used to classify mixture effects [31]. Recently Colby's formula has been extended to multi-compound mixtures [32]. As the IA ansatz is not sham-compliant it is not compatible with the Hill approach.
The ADM [13] or CA model of mutually exclusive action [1,14] for two noninteracting isoactive drugs A and B is where f x −1 is the dose or concentration of compound x that causes the specified effect. If the agents are acting according to Hill dose-response functions with slopes γ a and γ b , it is given by Adapting Berenbaum's approach, Greco derived a model for two-agent combined action by adding an interaction term, parameterized by a factor α. Assuming that the Hill-type dose-response curves of A and B differ only in the slope parameters, he gets [33] Although analytical expressions for u can be obtained from Eqs. 15 and 16 only under the restrictions of either a fixed maximum effect a max = b max and identical slope parameters γ a = γ b or of different maximum effects and identical slopes of unity γ a = γ b = 1, they are the starting points for several response surface models, e.g., Greco's model from Eq. 16 or the GCA expression [6] from Eq. 15. It permits different maximum effects but is limited to γ = 1.
Hence, u Greco (for α = 0) and u GCA are special cases of u Hill . The same holds true for Chou and Talalay's mutually exclusive model [11]. It was derived from the the median effect principle, assuming both a constant u max and γ Their mutually non-exclusive model [11] is an ad hoc extension of Eq. 19 Although it has been criticized by several authors [1,5,33] because of its questionable validity, it is one of the most often used models in the literature [34]. For γ = 1 it becomes the Bliss IA expression (Eq. 13, with α = β = 1 and u max = a max = b max ). Chou's models are related to Greco's expression for identical slopes γ = γ a = γ b and identical u max , i.e., u Greco = u Chou ex for α = 0 and u Greco = u Chou nex for α = 1.
Minto [30] proposed a model that solved the problem of the different denominators in Eq. 15 by expanding u max and γ in polynomials in a parameter . Fidler [28] extended Minto's approach by adding an interaction term. Their model is where α indicates the type of interaction. Minto's model corresponds to α = 0, α > 0 means synergism and α < 0 antagonism. f (s, w, p ) resembles a generalizeddistribution, and u max ( p ) and γ ( p ) are functions of the potency fraction p and f (s, w, p ). p ranges from 0 (drug A only) to 1 (drug B only).
Minto's model differs from the logistic Hill surface only in the functional forms of u max and γ . By truncating their polynomial ansatz for u max and γ after the linear terms in p , we have (Hill) Thus by making γ ( p ) and u max ( p ) symmetric with respect to m a and m b , u Minto becomes identical to u Hill . However, inclusion of higher powers of p leads to violations of the boundary conditions and the sham compliance requirement.

Examples
As identical theoretical concepts handling mixture effects have been developed by scientists from different disciplines, there is some confusion in the nomenclature used in the literature. We shall use the term 'Bliss independent action' in the subsequent discussions when referring to the probabilistic independent action ansatz, the exception being examples from agrochemistry, where we use the term 'Colby's formula' . This is justified because of its extensive use in the agrochemical literature.
In general, however, one has to keep in mind that a numerical evaluation of a model by comparison with experimental data is difficult because often the error bars of the experiments are large or unknown.

Crop protection agents
To demonstrate the applicability of the Hill model old mixture data from the crop protection area were chosen. All pesticides involved achieve an u max of 100%. For the pairs of atrazine/alachlor [12] (herbicides), aldrin/dieldrin [35] (insecticides), and oxadixyl/mancozeb [36] (fungicides) the root mean square errors of prediction are shown in Table 2. The mixtures are equally well described by using variable slopes or slopes of 1. All pesticides are assumed to achieve individual maximum effects of 100%. The magnitude of the deviations of mixture predictions from the experimental data is an indication of synergistic effects a LC 50 (in kg/ha) for herbicides, LD 50 (in g a.i./100 ml kerosene solution) for insecticides, EC 50 (in mg/l) for fungicides; b root-mean-square error of efficacy prediction for mixture For all pairs of pesticides the original conclusions are confirmed: neither the herbicidenor the insecticide-mixtures show deviations from the 'expected' response, whereas in the case of the fungicide mixtures the differences between expection and experiment indicate the presence of synergism.

Pyrethroids
Mixtures of pyrethroids can act both as synergists and antagonists. While permethrin/etofenprox and permethrin/cypermethrin show antagonism, cypermethrin/etofenprox act synergistically. These effects are only insufficiently described by percentages of 41, 58 and 332% [37] measuring the deviation from additivity. The full picture is given in Fig. 1, where Hill-, Bliss-and experimental dose-responses curves are overlaid with the Hill response surfaces. In this example the differences between the theoretical models are extraordinarily large. The contour plots of u Hill − u Colby reveal that especially for x ≤ x 50 and y ≤ y 50 Colby's independent action formula predicts much lower effects than the Hill model. As compared with the latter, the IA approach underestimates antagonism and overestimates synergism. These large differences are caused by the extremely steep dose response curves of permethrin (λ = 25), etofenprox (λ = 11) and cypermethrin (λ = 106), meaning that the effects vary strongly within very narrow dose-ranges.

Dioxin-like chemicals
Aryl hydrogen receptor (AhR) ligands were used to compare the toxic equivalence factor (TEF) approach and the more general GCA ansatz by [38] in predicting the expected effect of mixtures containing partial agonists or competitive antagonists. From their supplementary material the Hill curves with variable γ were derived and used to predict u Hill . GCA surfaces are Hill-surfaces with slopes of γ = 1. Hence, one might expect that Hill-surfaces with variable γ exhibit slightly improved Mann-Whitney (MW) statistics. However, as the slope parameters from the fits are only marginally different from 1 (c.f. Table 3), the differences between the GCA and Hill in the E max -, EC 50 -values and consequently in their response surfaces (not shown) are much smaller than the error bars of the experimental data, while the surface predictions are of comparable quality.

Response surfaces and isoboles
As the classification of mixture effects in terms of interaction indexes is not an adequate means to describe this complex phenomenon, an analysis of responses by looking at the form of iso-effect levels, isoboles, is certainly more appropriate [10,14,16]. Deviations from straight lines were used to classify the effects as (Loewe and Bliss) agonistic or antagonistic. From a response surface view [3,33,39] isoboles are cuts of the response surface at defined effect levels and the corresponding contour plots are the easiest way to get a full picture for the whole range of dose-combinations. This facilitates the understanding of the fact that linear isoboles are the exception and not the rule and are not a general means of detecting synergism. For critical discussions on the interpretation of the shape of isoboles see, e.g., [5,[40][41][42].
The comparison will be mainly between u Hill and u Bliss = u Colby , as the GCA model u GCA and Chou's mutually exclusive model u Chou ex are special cases of the logistic surface u Hill , and under special assumptions Chou's mutually non-exclusive model u Chou nex reduces to Bliss IA.
In Fig. 2 differences between the two surfaces show up at high doses of both mixture partners. There IA would find large synergistic effects of up to 30% where the Hill model would indicate 'null-interaction' .
Some trends can be observed: For identical maximum effects and slopes, u Hill − u Bliss is positive at doses below d 50 and increases with increasing u max and λ, and it is negative else, almost independent of the dose range. This means that Bliss independent action tends to overestimate synergism at doses smaller than their median effect doses and to underestimate synergism at doses above. At low doses the size of the difference increases with u max and λ. If the maximum effects differ, the effect differences at higher doses increase with increasing a max − b max . Some of these findings are illustrated in Fig. 3 for binary mixtures of agents with identical parameters: maximum effects (u max = 1 and u max = 0.4) were combined with three different slope parameters (λ = 0.5, 1.0, 5.0). Using scaled doses m = d/d 50 for the axes simplifies the picture without obscuring the essentials.

Iso-surfaces
For ternary mixtures the pendants of isoboles are iso-surfaces, i.e., projections of the four-dimensional response surfaces to three-dimensional ones defined by triples of doses leading to a specified response. A visualization of effects resulting from more than three agents in one graphical object is hardly possible. As shown for some fictitious ternary mixtures in Fig. 4, planarity is achieved only if all mixture components have identical slope parameters. In this case different maximum effects of the mixture partners do not affect planarity. Iso-surfaces of mixtures of agents having different slopes are non-planar. This is shown for ternary anesthetic mixtures [43] as analyzed in great detail by [28] and [30]. For the present purpose the data for midazolam, propofol, and alfentanil were fitted to Hill dose-response curves assuming effect ranges from 0 (no hypnosis) to 1 (full hypnosis). The characteristics of fits and of predictions for the mixtures are summarized in Table 4. The corresponding iso-surfaces in Fig. 5 show slight deviations from planarity which result from different slopes (4.8-11.1) of the pure agents. As the Hill model provides only null-interaction surfaces, the size of the deviations of experimental data from the Hill prediction may give some hints on the presence of synergy. While the RMSEs of fit for the pure compounds are of the order of 5%, those for mixture prediction are of the order of 40%, indicating the presence synergistic effects of increasing magnitude for the binary mixtures of propofol/alfentanil, midazolam/propofol and midazolam/alfentanil, while the ternary mixture is comparable to the best binary one.

Conclusions
Starting from a logistic PDE, analytical expressions for the response surfaces of ncomponent mixtures have been derived under the sole provision that each a.i. is described by a sigmoid dose-response curve. No further assumptions are required. The resulting 'null-interaction' surfaces, i.e., Hill-surfaces in the absence of of synergistic or antagonistic effects, provide the 'expected' response for each dose combination. Deviations from this reference response in order to quantify synergism or antagonism should possibly be handled by some sort of perturbation theory.
The Hill approach provides a framework to classify several models describing mixtures like ADM (Loewe CA, GCA), MSM (Colby, Bliss, Finney, Abbott) or Chou's 'unified general theory' and various response surface models. It can be applied to mixtures of compounds having different maximum effects and differing slope-parameters. Many Loewe-additivity based approaches are found to be special cases of the Hill surface while Bliss IA is incompatible with the logistic ansatz.
The independent action model is frequently used, e.g., in patent applications to quantify synergistic effects [44]. Its outcome should be checked by comparison with the Hill  For γ 1 the differences between Hill-and MSM-surfaces become extremely large, but they are restricted to those small dose-ranges where the effects change rapidly.
As a consequence for claiming synergistic effects at doses below the respective d 50 -values, special caution is required whenever the predictions from MSM and Hill approaches differ considerably.
An in-depth analysis of real mixtures would have to take into account not only errors in the models but also other sources of errors like the error statistics of experimental data.