1. Numerical methods

Hereafter we will call effector any agent able to cause a (typically sigmoidal) response in a population of biological entities, and perturbator to any agent that can alter the response to an effector, itself being unable to cause it. Receptor is any biological structure acting as ligand of effectors or perturbators. The term dose is reserved to the concentration of an effector.

In the simulation procedures that are described later, the Weibull's random numbers w:(θ;α) were obtained, from the uniform random numbers u:[0,1] provided by the spreadsheet, through [4], [5]:where mean (μ_w) and standard deviation (σ_w) of the corresponding distribution are:

To facilitate comparisons, doses were coded into the [0-(0.1)-1] interval, and responses were calculated by considering Y as a total number of biological entities, S the survive population at a given dose, therefore the surviving population response R can be expressed in a useful coded range [0,1] as R = 1−(S/Y). The simulated and experimental results were adjusted to the proposed models by non-linear least squares methods (quasi-Newton), in Microsoft Excel spreadsheet, using Solver complement for parametric estimates, and Solver Aid macro [6], [7] for confidence intervals and model consistency (Student's t and Fisher's F tests, respectively, with α = 0.05 in both cases).

In the most complex cases (models with interactions), the fitting process involved always a progressive hypotheses of contrast, as it is usual in any stepwise multiple regression method, to select the interactive mode providing the most statistically consistent interpretation. Parametric confidence intervals, coefficient of multiple determination, residual bias and sensitivity analysis were applied as selection criteria. An efficient way of proceeding is by following the next steps: 1) calculation of the sigmoidal parameters from the individual responses; 2) use of these estimates as provisional fixed values of the model, and assay of different interaction hypotheses, rejecting those that lead to not statistically significant coefficients; 3) refinement of the model by recalculation, now allowing the variation of all the accepted parameters.

Although the initial number of parameters is high, it means only a high number of potential alternatives, some of which are mutually exclusive, and others easily rejected in the course of the fitting to a concrete data set. Trying of initial values is facilitated by follow up of the variations produced by the Solver results in graphic displays of response surface and residuals, and convergence is usually immediate. Details about experimental design are provided as Supporting Information (see Figure S1).

2. Null interaction, synergy and antagonism. IA and CA hypotheses

In any system (as defined in the Bertalanffy's sense: a set of interacting elements), an important and characteristic problem is to know whether the joint effect of two or more elements on the system behaviour is deducible from their individual effects. This issue, with a long history of controversy whose first known attempt goes back to Aristotle, is often stated by replacing «deducible» with «the sum», what leads to define the notions of synergy and antagonism as those interactions by virtue of which the joint effect of two (or more) effectors is greater (synergy) or lesser (antagonism) than the sum of the individual effects.

For examining the meaning of this non-acceptable statement, let us to consider a system in which two effector elements (E₁ and E₂) can act, and let us to denote the possible behaviours or responses of such a system as R₁ and R₂ (if only E₁ or E₂ are present), R_1,2 (if E₁ and E₂ do not interact), and R_1&2 (if E₁ and E₂ interact). Under these conditions, to obtain any initial analysis on the possible interactions, it is an essential requirement to compare the responses R_1,2 against R_1&2, defining synergy and antagonism respectively, as those interactions in which R_1,2<R_1&2 and R_1,2>R_1&2.

It is obvious that the key aspect to assess any interaction is the response of the system in the absence of interactions (the null interaction). Thus, to define the null interaction, any mechanism must be postulated as underlying any specific behaviour of the system. The addition is the simplest mechanism, and its most immediate options suppose that the added magnitudes can be the effectors (acting as one alone), or the effects on the system of their independent actions (the responses of the system to such actions).

Now, if we imagine that the response has a superior limit, as it happens, for example, in the case of the death-survival alternative in a microbial population under increasing doses of two toxics. If the added magnitudes are the concentrations of the toxics, the response obtained, would be identical to the resulting response values of a single toxic dose –although it can be expressed as a function of two independent variables–. The additive response is more problematic, since it is obvious that if one cell dies when one of the doses reaches a given level, it will die independently on the level of the other dose, simply because it cannot die twice.

These two types of sum are the foundations of the two basic accepted modes for describing the joint action of two effectors under null interaction conditions: concentration addition and independent action. Although both can be applied –at least in theory– to any number of effectors, here will be discussed in their simplest forms, for two effectors (their generalizations are, anyway, immediate).

2.2. Concentration addition hypothesis.

In its classical formulation [2], [3], null interaction is not defined as a relation between the individual responses, but through the following criterion: since the concentration (c) of an effector whose action obeys the equation R = f(c) can be considered as a fictitious combination of c₁ and c₂ concentrations (c = c₁+c₂), it is obvious that the response to c will be described by the equation R = f(c), with c = c₁+c₂. If the response to a mixed dose of two effectors behaves as the response to the “mixed” dose of the same effector, it is accepted that the interaction between them is null, implying that any effector concentration can be substituted by the equieffective concentration of the other one.

The conventional practice avoids an explicit algebraic formulation when evaluating the joint response of effectors, using the isobolographic analysis [9] (or lines on the plane of the independent variables that represent the dose combinations that produce an equal response) to solve the lack of a clear mathematical process. Thus, in the isobolographic analysis, if D₁ and D₂ are the doses of two effectors that produce the individual response R_a, and d₁ and d₂ any dose combination that produces the same joint response R_a (Figure 1), under null interaction conditions the isobole of the response R_a will be necessarily described by the following lineal equation:(3)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Isobole of a response R_a.

In 1, 2 and 3 the geometric logic underlying the analysis of the CA hypothesis through the equation (4) is shown. Type 4 isoboles arise in many real responses corresponding to the IA hypothesis with null interaction, and illustrate the limitations of the relation between factual and formal aspects of isobole analysis beyond a particular case of the CA mode of action (see results, section 5).

https://doi.org/10.1371/journal.pone.0061391.g001

Consequently, if the individual DR models are R_i = f_i(D_i) and their reciprocal functions D_i = g_i(R_i) exist, it can be established that:(4)

In other words: straight isoboles indicate null interaction, and concave and convex up isoboles indicate synergy and antagonism, respectively. The dimensionless quotients d_i/g_i(R_a) are called toxic units and represent the relative contribution of each effector to the joint response R_a.

3. Some problems associated with the AI-AC approach

A first unsatisfactory aspect of this approach is the difference between the formal criteria applied to each mode of action. IA hypothesis proposes an explicit response surface model, but a general agreement does not exist about the ways in which synergy and antagonism must be formulated. CA hypothesis avoids the explicit model, but equation (4) is accepted as a criterion for detecting synergy and antagonism. This criterion, however, is not transferable to the IA framework, whose mathematical form prevents straight isoboles in null interaction if –as it occurs in most DR relationships– the individual responses R₁ and R₂ are not proportionally constant. As a general rule, IA isoboles in null interaction are convex up at low doses, concave up at high doses, and with two branches of opposite curvature in a transitional zone when R₁≠R₂ (Figure 1). But the factual meaning of this formal property cannot be attributed to synergy or antagonism, only implying that the probability that at least one dose is lethal is low at low doses and high at high doses (another outcome of the statistical independence, as the above mentioned impossibility of dying twice).

In CA hypothesis, as already pointed out, the absence of an explicit mathematical model is another disadvantage. A measure of the sign and degree in which an isobole deviates from linearity is provided by (4), and isoboles with a variable degree of curvature and asymmetry can be defined by using alternative expressions as equation (5) described by [10] or equation (6) described by [11].(5)(6)where the additional parameters β or γ are more precise indexes of synergy and antagonism. Although these equations would enable the construction of isobole maps by solving them numerically, the procedure is laborious and in practice is usually applied only to the isobole at the half-maximum response. Furthermore, in this approach a homogeneous isobolic curvature is postulated for the whole considered domain, a property that as we shall see later, is not necessarily true.

Some joint responses to chemically similar/dissimilar effectors were suitably described with CA/IA models, respectively [12], [13]. However, there are evidences, as well, that these modes of action are not always obeyed by the reality. Jonker et. al [14] have proposed in each case, besides synergy and antagonism, other deviations from null interaction which were defined a priori as effects depending on the absolute dose levels or dose ratios, thus enabling synergistic and antagonistic displays on different regions of the same response surface. But even so, the reality seems to be richer: in a revision of 158 cases [15] was found that 20% of the responses can be adequately predicted by IA, 10% by CA, another 20% admitted both models and half of the cases were not correctly described with either of them. Moreover, neither of the models was significantly better than the other on assessing synergy and antagonism at the 50% effective doses.

These results are not really surprising. Figure 2 represents a simple hypothetical metabolic pathway that can lead to a set of different situations depending on the considered inhibition mechanisms (competitive, non-competitive, acompetitive), the kinetic constants that are null, the values of the rest of such constants and the definition of the response (individual or simultaneous drops of the products P₁, P₂ or P₃). But even without including an experimental error, none of these possible situations obeys unambiguously IA or CA modes of action. In our laboratory, these ambiguities have been detected not only in the joint action of hydrocarbons and dispersants on the larval growth of sea urchin [16], but also, as we will see, in the simpler context of the oxidation inhibition of a substrate by the joint action of two antioxidants.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Ambiguity of the IA-CA dualism.

Rates of the enzymatic reactions yielding products P₁ to P₃ from substrates S₁ and S₂ are affected by the inhibiting effectors E₁ and E₂ with the specified inhibition constants k_ij. Under these conditions, responses measured as drops of the levels of any of the products are dependent of the nature (competitive, non-competitive, acompetitive) of the inhibition and the values of the inhibition constants. But any result can be unambiguously attributed to IA or CA modes of action.

https://doi.org/10.1371/journal.pone.0061391.g002

4. DR model for a single effector

The natural form of a DR model is a cumulative (mass) probability function, and it translates the response of a population with a given sensitivity distribution to an effector. Four additional conditions seem reasonable as well: 1) the model should have an explicit algebraic form; 2) it should be lacking of an intercept (null response at null dose); 3) an asymptote equal or lesser than 1 should be enabled; and 4) the parameters with important factual meaning should be explicitly included, to facilitate the trial of initial values and the calculation of confidence intervals when non-linear fitting methods are applied.

Although normal and log-normal distributions have been the basis of the classical DR analysis, they have the disadvantage of lacking of an explicit form for their mass functions. Logistic-type equations are more useful and they can be expressed in forms easily modifiable to comply with the above mentioned conditions [17], [18]. However, their derivatives (their density functions) show only right bias, which can be a restriction scarcely realistic. Another option is the Gompertz [19] equation, but its use is prolix, especially with the modifications that are required to apply in the DR context. The mass function of the Weibull distribution [20] can be expressed in a suitable reparametrized form [21]–[23], providing the adequate DR analytical model, whose expression can be write as: (7)where D is the dose, R the response (with K as asymptotic maximum, not necessarily 1), m the dose producing half-maximum response, and a a shape parameter related to the maximum slope of the response (r_m) by:(8)

It should be noted that m is the abscissa of the inflection point, which represents the accumulated modal response. The basic parameter of the DR analysis is ED₅₀ or EC₅₀, defined as the effective dose for the 50% of the assay population. Thus, m and EC₅₀ are coincident values if the maximum response implies the whole population (K = 1). The reciprocal function of equation (7), which is a requirement when performing the isobole analysis in the CA hypothesis, it can be written as:(9)

The use of the equation (7) as DR model is interesting for several reasons. Its density function (the sensitivity distribution of the population) can be symmetrical or asymmetrical with right or left bias, which makes it very versatile. It produced the best fittings, among the above mentioned alternatives, when it was applied to the simulations that will be described in the next sections, and this result was repeatedly confirmed before by experimental data [21], [24]–[27]. Moreover, the Weibull distribution is the conventional model for the failure of complex devices, making the equation more attractive because it unifies phenomena in which underlies a profound analogy.

5. Simulation of the response to a single effector

Since the basic sigmoidal profile of the DR relationships translates a macroscopic, statistical phenomenon, it should be possible to simulate it as a result of the microscopic behaviour of a population of elementary biological entities, which we will call cells. Such a simulation can be carried out on the following basis:

B1.One cell is defined by means of three aleatory magnitudes: ρ, or number of receptors of an effector, α, or number of active receptors –ready for linking the effector– in a given instant, and λ, threshold or minimum number of active receptors that must be linked to an effector to produce a response.

B2.The dose D is defined as the number of effector units per cell, accepting that every unit is capable of linking to one receptor.

B3.The cell response r is limited to two modalities: death (r = 0) and survival (r = 1), obeying the following logical rule:

(10)or, as a Boolean proposition (: and, : or, : no, 1: true, 0: false):

Notice that α<λ implies r = 1 at any dose. It translates possible limitations in the effector bioavailability, resistant cells or other conditions which, relatively frequent in DR assays, produce lesser than 1 asymptotes.

Although ρ (highest limit of α) is interesting in some cases, hereafter it will be omitted without loss of generality, and a cell will be defined through the pair α, λ. Thus, if we assign to α and λ probability distributions defined by their mean and variance values, we can create, in a spreadsheet, a virtual cell population whose sensitivity distribution depends on the parameters α and λ. The population response R to an increasing dose series is simulated by applying the rule (10) to each cell and defining R = 1−(S/Y) when S cells, of a total number Y, survive at a given dose. The typical shape of the response arises even at moderate population sizes (Y∼100) and becomes highly stable when Y≥1,000. At low population sizes, the variability of the result represents a simulation of the natural variability, the experimental error, or both.

These premises define the minimum complexity of a system able to generate a sigmoidal response and, despite their schematic character, they produce a great variety of profiles, depending on the α and λ parameters. Although the microscopic solution of the system is determined by such parameters, the obtained profiles can be macroscopically described by the conventional equations of the DR analysis.

The use of the Weibull distribution for defining α and λ is not essential. With normal, log-normal, Poisson, binomial, or even uniform distributions, the best descriptions of the responses are obtained, as it has been mentioned. However, distributions with domain [0,∞) are obviously preferable, since low means and high variances in (−∞,∞) distributions lead to a finite probability of negative values of α and λ, which lacks of physical meaning.

6. Simulation of perturbations of the response to an effector

Using the elements above defined, we will admit that a perturbator does not produce a response by itself, but it can modify the response to an effector by altering the number of active receptors (α), the threshold (λ) or the effective dose corresponding to the nominal dose (D) (hereafter α, λ and D-perturbations, respectively). α and D-perturbations can be illustrated in molecular terms through the well-known key-lock analogies (Figure 3). λ-perturbations require to suppose an intermediate fast process modifying the cell sensitivity.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Some key-lock analogies.

Possible response modifications involving alterations of the effective dose or the number of active receptors are illustrated (see Table 1). Notice that the alterations of the effector-receptor affinity (for simplicity reasons only perturbations are illustrated) do not modify the response to a given dose, but the response to a given time.

https://doi.org/10.1371/journal.pone.0061391.g003

All of them can be exemplified by common physiological mechanisms as those that take place in trans-membrane proteins, second messengers or enzymatic systems, in which any DR assay allows to distinguish at least between α and D-perturbations. Indeed, in the presence of an excess of effector, a moderate increment of a perturbator modifies or not the response depending on whether the perturbation acts over the receptors or the effector.

The effect of these perturbations on the response can be simulated by using the rule R0 (10) and adding a vector that represents increasing perturbator concentrations, as well as a criterion for modifying the values of α, λ or D as a function of the perturbator concentration. Since direct or inverse ratios are the simplest criteria, a perturbation term can be formulated as:(11)where P is the perturbator concentration, p_ε the proportionality coefficient, and the ε subscript indicates the element affected by the perturbation. The term U_ε multiplies or divides the values of D, α or λ depending on the effect that we are trying to achieve, and it can be replaced by any other algebraic expression able to describe any other type of alteration.

Now, it can be pointed out that if a DR curve is sigmoidal because the most sensitive elements of a population die at lower doses than the most resistant ones, therefore a time-response curve will be sigmoidal if the sensitivity distribution of the population is translated into responses at shorter or longer times. Although this problem will not be considered here, the time-course of the response could be treated also in terms of molecular interactions, which in this case, would accelerate or delay the progress of the effect (Figure 3). In some fields –i.e. marine toxins–, the time of death of a target organism often replaces the response in the framework of the empirical dose-time death models [28]–[31]. We do not allude here to this approach, but to the time-response relationships, which can be described by using the same models as DR analysis [26].

7. Simulation of the response to two effectors

In this case, the perturbator must be replaced by a second effector, but the essential issue is the need to include in the simulation algorithm, even in the absence of any interaction, any rule about the joint action. The simplest hypothesis in this regard leads to admit that: 1) the receptors are specific of the effectors; 2) the cell response is r = 0 if any of the doses exceeds the corresponding threshold and there is a sufficient number of receptors. Thus, the rule is:(12)

These conditions are reminiscent of IA hypothesis and in fact, as we shall see later, they produce the typical IA response surfaces. Now then, both conditions can be denied, what seems quite reasonable for the second one [32], since it prevents to accept that two biological subsystems –e.g. glycolysis and β-oxidation– can be affected at individually sub-lethal, but jointly lethal levels. By denoting the first condition (specificity) as S⁺ and the second one (independence) as I⁺, three additional rules, besides S⁺I⁺, can be considered:

S⁺I⁻. It admits that 1) the effect Gi of a dose Di below threshold λi is Gi = Di/λi; 2) Gi values are additive; 3) r = 0 if G1+G2≥1 and there enough receptors:(13)

S⁻I⁺. It admits that any effector has access to the whole of the receptors (α₁+α₂), what produces competence if they are insufficient. Competence depends on factors as the relative doses of the effectors, their diffusivity or their affinity for the receptors, but to simplify, only relative doses will be considered here: C₁ = D₁/(D₁+D₂) and C₂ = D₂/(D₁+D₂). Thus, the number of receptors linked to D_i will be C_i(α₁+α₂), what leads to the rule:(14)

S⁻I⁻. It admits simultaneously competence (if the receptors are insufficient) and additivity of below-threshold effects. The rule is:(15)

None of these rules implies the sum of the doses required by the CA hypothesis. This condition plays down the receptor specificity, but again the competence arises if the receptors are insufficient. If both effectors have the same threshold, the response is not modified by the competence. If the thresholds are different, the number of occupied receptors will depend on the relative doses and will be C_i(α₁+α₂), turning the dose addition effect in below-threshold addition effect: G_i = C_i(α₁+α₂)/λ_i. Therefore, two rules are possible, any of which produces the response surfaces with straight isoboles that characterize the null interaction in the CA hypothesis:

R5a: Without competence (equal thresholds):(16)

R5b: With competence (different thresholds):(17)

8. Inherent and accidental mechanisms

Concrete assumptions (about specificities, thresholds, competence and dose or effect addition) as those included in the rules R1 to R5 are indispensable for any joint action hypothesis, and they represent mechanisms that are inherent to a given null interaction model. Over these minimal inherent modes of action, accidental mechanisms (interactions) can be superimposed, in which each effector can perturb the response to each other by modifying the number of active receptors, the threshold or the effective dose corresponding to the nominal dose. Thus, when the effector E₁ perturbs the response to E₂, we can write an interaction term like (11):(18)which is included into the rules R1 to R5, as a factor or divisor of D₂, α₂ or λ₂.

The interactions can be unidirectional (E₁ alters some factor related with E₂, but E₂ does not alter E₁ accordingly) or reciprocal (mutual alterations), and these last ones can or cannot be symmetrical (the same or different strength in both directions). Whereas the reference of a perturbation is the response in the absence of perturbator, the reference of an interaction is the inherent mechanism or null interaction. This implies that a given interaction has different consequences depending on the inherent mechanism and the frame of application. A systematics in this regard is proposed in Table 1, which attempts to preserve the main senses of the usual terminology [33], and whose categories allow simulations and specific formal descriptions. Although this terminology can be used without ambiguity, it seems preferable to simplify it by defining stimulation/inhibition as perturbations that increase/decrease the response to an effector, and sinergy/antagonism as interactions that increase/decrease the joint response to two effectors with respect to the response promoted by the null interaction.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. A possible systematics on the modifications of the response to an effector.

https://doi.org/10.1371/journal.pone.0061391.t001

1. Perturbations of the response to a single effector

In Figure 4 the three conditions that depress the response simulations to an effector in the presence of increasing concentrations (P) of a perturbator are shown (decrease of the effective dose, decrease of the number of active receptors and increase of the threshold). The individual fittings to the resulting profiles with the model (7) proved that, in each series, the estimates of K and m parameters varied as a function of P in specific ways (Table 2), thus enabling to identify the underlying perturbation. Although the maximum slope varied as a consequence of the variations of K and m, in agreement with (8), the parameter a remained constant in all cases.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Effect of a perturbator on the response (R) to a same dose series (D) of an effector and the parameters of the model (22).

The three cases in which the perturbation depresses the response are illustrated: reduction of the effective dose corresponding to the nominal one (left), reduction of the number of active receptors (canter), and increase of the threshold (right). Dots are simulated results in the absence (○) and presence (•) of increasing concentrations of the perturbator, and lines the respective fittings to model (22). See also tables 2 and 3.

https://doi.org/10.1371/journal.pone.0061391.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Variations (+: increase; −: decrease; 0: no change) in the parameters of the response to an effector, as described by the equation (7), due to the presence of an agent which produces the specified perturbations.

https://doi.org/10.1371/journal.pone.0061391.t002

Such a constancy of a is not surprising, since the simulations with the rule R0 prove that the variation of this parameter is related with the variations of the variances of α and λ, a condition that was not considered in any case. Such a restriction simplifies the analysis and is not arbitrary. From the microscopical point of view adopted for the effector and perturbator actions, variations in the number of receptors and threshold are clear possibilities, but it is more difficult to justify an action on a populational property such as the variance.

1.1. The perturbation function.

To obtain a simultaneous solution for every series of profiles, an auxiliary function (a perturbation function π_θ) is required for describing the possible variations of any parameter θ as a function of P (Figure 5). This can be achieved by means of different biparametric (exponential, polynomial or hyperbolic) expressions:where the θ subscript represents the modified parameter (K, m), θ₀ is the parametric value when P = 0, and the pairs b_θ, c_θ are fitting coefficients. It should be noted that, in the absence of perturbation, the first function requires b_θ = 0 and c_θ = 1, whereas the other two require b_θ = 0 and c_θ = 0. Moreover, the condition of the third function a singularity is obtained if c_θP = −1. In order to avoid it, when P is coded in the interval [0,1] it is advisable to include the restriction c_θ>−0.999 in the fitting algorithm. Thus, the model (7) turns into:(22)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 5. Possible variations of a parameter (θ) of the response to an effector, as a function of the concentration of a perturbator (P).

Any of the functions (19), (20) and (21) can produce all the profiles. Parametric value can be increased (+) or decreased (−), with constant (L), decreasing (A) or increasing (C) slope.

https://doi.org/10.1371/journal.pone.0061391.g005

With any of the mentioned forms of π_θ, the (22) led to excellent simultaneous fittings, which confirmed the specificity of the parametric variations found in the individual descriptions. Figure 4 and Table 3 show the results obtained with the hyperbolic function (21), which will be the form used now on. The description of the inhibition of the haemolytic activity of palytoxin by ouabain [22], illustrates the application of equation (22) for perturbations involving a time-response profile.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Simulation conditions of the responses to an effector as perturbed according to the three modalities that cause response drop, and respective fittings to the model (22).

https://doi.org/10.1371/journal.pone.0061391.t003

When the model (22) is satisfied with uniparametric π_θ (b_θ = 0 or c_θ = 0), the relationships between confidence intervals (ci) and parametric values are approximately of the same numerical size in all the cases. When b_θ≠0 and c_θ≠0 are required, both ci are penalized by the linear correlation between both coefficients, and therefore, with high experimental error and low number of observations, some b_θ, c_θ pair can become not statistically significant, even in highly predictive models. This problem is solved if the model is recalculated fixing the value of one of the coefficients and excluding it of the Student test, or –with a small loss of predictive capability– making zero the coefficient of the pair whose suppression does not alter the increasing or decreasing trend of the parametric variation due to the perturbation function.

Finally, in D-perturbations is interesting to observe that two equivalent solutions are possible. One of them is that provided by the model (22), describing the variation of m (the only affected parameter) due to the presence of the perturbator. The other one describes directly the variation of the effective dose through the equation:(23)

Later on we will see that this dual solution is not possible in interactions between two effectors.

2. The response surfaces to two effectors

Among the five basic types of response produced by the rules R1 to R5, only two of them were in accordance with those corresponding to IA and CA hypotheses, and the diversity increases as interactions are included. Two possible reasons to avoid describing this diversity as unrealistic are because it involves options that are implicit in the same rules that produce IA and CA responses and it translates common physiological mechanisms. In fact, the situation seems to be just the opposite one. As Rovati et. al [34] pointed out regarding the cases of partial agonism, or effectors able to interact with receptors promoting opposite effects, there are responses «that were often disregarded by the experimentalists, or considered as artefacts, in the absence of a biological and/or mathematical theory to justify them».

The parametric variations due to interactions showed the same specific increasing and decreasing trends as those due to perturbations (Table 2), in linear or asymptotic forms (Figure 5) depending on quantitative factors, and such fact is hardly surprising.

3. Modelling of interactions under the independent action hypothesis

As expected, the response surfaces produced by the rule R1 were consistent with the IA hypothesis (Figure 6, Table 4) and could be accurately described by transferring model (7) to equation (2):(24)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 6. Joint response to two effectors in the four suppositions resulting from combining the two implicit key conditions of the IA hypothesis (see text).

In the first column, dots are the result of simulations and surfaces the respective fittings to the model (26). Isobolograms, correlations between observations and predictions and parametric variations (K_i: ○, m_i: •) of the response to an effector as a function of the dose of the another are added. D₁ and D₂: doses; R: response. Numerical data in Table 4.

https://doi.org/10.1371/journal.pone.0061391.g006

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Simulation conditions and respective fittings (α = 0.05) to the generalized IA model in the specified examples.

https://doi.org/10.1371/journal.pone.0061391.t004

Since each effector can act as perturbator of the response to the other, by altering effective doses, active receptors or thresholds, auxiliary functions π_θi can be defined in terms as those applied to the perturbations. Thus, the (24) will be written, in its most complex form as: (25)

When this equation was used to describe the corresponding simulations, the specific variations in the parameters K_i and m_i led to discriminate all the modalities of interaction (Table 2), whose main types are summarized in Figure 6, Figure 7 and Table 4. For those interactions affecting the effective dose, the dual solution that is possible in the homologous case of the perturbations, here, it would be reduced to parametric variations affecting only to m_i.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 7. Some examples of joint responses to two effectors under IA mode of action.

Concrete types of interactions are specified and adjusted to the generalized model (26). Keys and graphic criteria as in Figure 6. See also Table 4.

https://doi.org/10.1371/journal.pone.0061391.g007

4. Modelling of interactions under the concentration addition hypothesis

As it has been already mentioned, the application of the CA hypothesis is usually carried out through the isobole analysis. However, the definition of null interaction according to Berenbaum provides the key for establishing an explicit model. Indeed, if the response to a mixed dose of two effectors should behave as a fictitious mixed dose of a single effector, the model is necessarily:(27)

In fact, the simulations obtained with any of the rules R5 were accurately described by this equation, which produced straight isoboles with equal intersection points on the doses axis (Figure 8 and Table 5). Although the competence affects the parametric values, it does not alter the functional form, making equivalent the two alternatives of the rule R5.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 8. Joint response to two effectors under CA mode of action.

Concrete types of interactions are specified and adjusted to the generalized model (32). Keys and graphic criteria as in Figure 6. See also Table 5.

https://doi.org/10.1371/journal.pone.0061391.g008

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Simulation conditions and respective fittings (α = 0.05) to the generalized CA model in the specified examples.

https://doi.org/10.1371/journal.pone.0061391.t005

As a consequence of the equation (27), the CA hypothesis can be accepted when the individual responses differ in their m parameters and –given the relation (8)– in their maximum slopes, but should be rejected when such responses differ in their K parameters. This is consistent with the isobole approach: if the asymptotes differ, expressions as (4), (5) or (6) could only be applied to lower responses than the lowest asymptote, since the inverse function (9) only exists if R<K.

When interactions are included, the key notion of concentration addition should be preserved in any modification of the (27), what means that the doses should act as an additive block in a function with a single set of sigmoidal parameters (K, m, a). These conditions enable the cases that are described next.

4.2. Synergy and antagonism.

To obtain surfaces with isoboles like those associated with synergy and antagonism, the rule R5a must include the condition that an effector alters, unidirectional or reciprocally, the effective dose of the other one. By using π_Di terms like those π_θi defined by equation (22), the corresponding simulations (Figure 8 and Figure 9) were described by means of:(29)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 9. More examples of joint responses to two effectors under CA mode of action.

Concrete types of interactions are specified and adjusted to the generalized model (32). Keys and graphic criteria as in Figure 6. See also Table 5.

https://doi.org/10.1371/journal.pone.0061391.g009

Contrary to the case of the perturbations, the equivalent dual solution based in the variation of the m parameter is not possible here.

5. Broad and strict sense of the notions of synergy and antagonism

The frequent and not always enlightening discussions about synergy and antagonism have led some authors to consider synergy as an «ineffective mixture risk definition» [10]. However, we believe that both concepts are important, and that their confused sides can be debugged by suppressing some common, not justified assumptions.

Firstly, their formal aspects should be distinguished from the factual ones. From a factual perspective, as it has been said, synergy and antagonism are the result of interactions that increase or decrease the response with respect to that expected in null interaction. This is a broad and unambiguous definition, but it requires taking into account the following issues:

1. Synergistic and antagonistic responses can be generated by any of the interactions (D, α or λ) considered in preceding sections. This implies that such responses can have diverse origins, that these origins require different formal descriptions, and that these descriptions are dependent not only on the elements involved in a given interaction, but also on the inherent mechanism of the system under null interaction conditions.
2. Association between concave/convex isoboles and synergy/antagonism is limited and misleading. In IA mode, as it has been seen, it lacks sense, and in CA mode it lacks general validity. As an example, the reciprocal synergistic α-interaction (Figure 9) increases the response in the entire domain with respect to null interaction, in spite of which its isoboles are straight. In fact, that association is only applicable to D-interactions in CA mode, where the corresponding regular series of concave/convex isoboles enable clear contrasts with the straight series which are typical –although no exclusive– of the null interaction.
3. If the usual isobolographic convention in CA framework is followed, synergy/antagonism could be defined in a strict sense as those effector-effector interactions that increase/decrease the response. Such a definition, however, would exclude arbitrarily other interactions (as α or λ ones) with similar net effects.
4. Unfortunately, satisfactory expedients for typifying synergy/antagonism by means of a single value or a particular isobole do not exist, since the differences between these conditions and null interaction vary throughout the corresponding response surfaces (Figure 10C2).
5. Furthermore, there are not theoretical reasons which prevent interactions with opposite regional effects on the response surface (e.g. E₁ increases the effective dose of E₂, and E₂ increases the parameter m of the response to E₁). Consequently, synergy and antagonism can be simultaneously detected in different regions throughout a given response surface.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 10. An example of ambiguous interpretation of a response through IA and CA hypothesis.

A: IA response surface (defined by K₁ = K₂ = 1.00; m₁ = m₂ = 0.30; a₁ = a₂ = 2.00). B: result obtained when this surface is interpreted, through CA hypothesis, as a case of symmetrical antagonism (K = 1.00; m = 0.316; a = 2.348; c_D2 = c_D1 = 1.311. All coefficients statistically significant, with α = 0.05). Notice the lack of residual randomness. C: differences between CA minus IA (C1) and CA null interaction minus CA antagonism (C2) responses.

https://doi.org/10.1371/journal.pone.0061391.g010

6. More insufficiencies and ambiguities of the IA-CA dualism

The frequent inconclusive character of IA and CA hypothesis is a well-documented experimental fact [39], whose justification was suggested in theoretical section 3. A different justification is provided by the simulations described Simulation of perturbations of the response to an effector and Partial forms of independent action sections. Indeed, when the conditions that define the IA hypothesis are altered in biologically plausible forms, response surfaces with cooperative or competitive effects are obtained, which cannot be acceptably described by any of the possibilities of the IA-CA scheme.

An additional cause of ambiguity –difficult to detect in the absence of explicit models– derives from a more formal issue. The IA and CA response surfaces are in general clearly distinguishable when the asymptotes of the individual responses are less than 1, as in the examples of Figure 6, Figure 7, Figure 8 and Figure 9. But the distinction turns more problematic as these asymptotes move closer to 1, since in such a case the IA surface losses its peculiar top region, in which the joint response surpasses the asymptotes of the individual responses. Figure 10 shows an IA surface obtained by assigning arbitrary parametric values (with K₁ = K₂ = 1) to the model (1), which could be significantly typified as a case of antagonism in CA mode by any assessment method, especially if the results are “blurred” by the experimental error. Similarly, a reciprocal asymmetrical synergy in IA mode is practically indistinguishable from synergy in CA mode. In such cases, the false hypothesis could only be detected by the lack of randomness of the residuals, if the number of observations is sufficient and the experimental error is reasonable.

All these reasons lead to doubt about the generalization to more than two effectors of any IA-CA discrimination method, since the probability of all kind of ambiguities increases with the number of agents considered. The main justification of this generalization is the experimental economy in the research of the joint effect of many effectors at moderate levels, an important issue in environmental toxicology, which seems practically unapproachable through the assay of binary combinations. However, it seems as well that the simplification of the experimental arrays only will produce even more ambiguous results.

Experimental examples

Some elements of the approach proposed here are exemplified by two above mentioned cases of study: the larval growth inhibition in sea urchin by the joint action of hydrocarbons and dispersants [16], and the inhibitory perturbation by ouabain of the haemolytic action of palytoxin [27]. In fact, the need to clarify problems as those arose in these cases was the origin of the simulation-based systematics we have attempted in this work.

We present now an experimental example in connection with the antioxidant activity, a field in which the possible interactions between both natural and synthetic products are formally equivalent to the toxicological ones, and they raise similar discussions. The example refers to the inhibition of crocin oxidation by the joint action of two well-known antioxidants: ascorbic acid and trolox. Triplicate analytical results were obtained by monitoring of the oxidation kinetics in 8×8 concentration arrays, using a microplate method [40], [41]. Inhibitory responses were quantified through the ratio between areas under kinetic profiles at the end point of the control [40], [42]–[44].

Antioxidants can compete with the oxidizable substrate for oxygen or the source of radicals (primary antioxidants), or for radicalized products that are formed in more advanced oxidation phases (secondary antioxidants). Therefore, their modes of action are in agreement with the diagrams of Figure 2 and they can be analysed using the general equations [31] or [37].

In this regard, the example of ascorbic acid and trolox is interesting, because it shows some of the features in which a clear decision is difficult (Figure 11 and Table 6). Although null interaction cannot be accepted under both IA and CA hypotheses, when a synergistic effect is admitted, the difference between the two (statistically significant) options becomes small. Residual distribution inclines the decision towards IA, but a more accurate characterization is a predominantly independent action, with synergy and a cooperative unspecific effect. On the other hand, it should be pointed out that a conventional analysis (use of a model (1) for IA hypothesis and a contrast on the 50% isobole for CA hypothesis) would lead to decide a CA

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 11. Joint effect of ascorbic acid (A₁) and trolox (A₂) on crocin oxidation under different hypotheses.

Keys and graphic criteria as in Figure 6, with parametrtic variations replaced by residuals, which are more informative in this case. See details in text, and numerical results in.

https://doi.org/10.1371/journal.pone.0061391.g011

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 6. Antioxidant joint action of ascorbic acid (A₁) and trolox (A₂) on crocin oxidation.

https://doi.org/10.1371/journal.pone.0061391.t006