Assessing the Effect of Mycotoxin Combinations: Which Mathematical Model Is (the Most) Appropriate?

In the past decades, many studies have examined the nature of the interaction between mycotoxins in biological models classifying interaction effects as antagonisms, additive effects, or synergisms based on a comparison of the observed effect with the expected effect of combination. Among several described mathematical models, the arithmetic definition of additivity and factorial analysis of variance were the most commonly used in mycotoxicology. These models are incorrectly based on the assumption that mycotoxin dose-effect curves are linear. More appropriate mathematical models for assessing mycotoxin interactions include Bliss independence, Loewe’s additivity law, combination index, and isobologram analysis, Chou-Talalays median-effect approach, response surface, code for the identification of synergism numerically efficient (CISNE) and MixLow method. However, it seems that neither model is ideal. This review discusses the advantages and disadvantages of these mathematical models.


Introduction
Mycotoxins are secondary metabolites mainly produced by fungi belonging to the genera of Aspergillus, Penicillium, or Fusarium [1]. Although the role of mycotoxins is not yet fully understood, it has been shown that mycotoxins form an integral part of microbial interactions in ecological niches where they protect fungi from competing or invading microbes (e.g., by antimicrobial activity and/or quorum sensing disruption) [2,3]. Throughout history, these fungal toxic metabolites have been recognized as harmful contaminants in crops, causing acute toxic, carcinogenic, mutagenic, teratogenic, immunotoxic, and oestrogenic effects in humans and animals [1,4]. From the public health point of view, the most important foodborne mycotoxins are aflatoxins (AFs), fumonisins (FBs), trichothecenes (including deoxynivalenol (DON) and T-2 and HT-2 toxins), ochratoxin A (OTA), patulin (PAT) and zearalenone (ZEN) and maximum levels have been set in European Union legislation to control these mycotoxin levels in food and feed [4,5]. Analytical methods based on the liquid chromatography tandem mass spectrometry (LC-MS/MS) have been developed for the simultaneous detection of multiple mycotoxins in foods which facilitated and enabled survey of their co-occurrence in various food matrices [6,7]. This methodology enabled the simultaneous detection of more than one hundred fungal metabolites including major mycotoxins as well as masked (e.g., DON-3-glucoside The most comprehensive review on mycotoxin interactions in cell cultures of human and animal origin was given by Alassane-Kpembi et al. [18]; the majority of conducted studies used the arithmetic definition of additivity. In the studies conducted in the last four years (Tables 1 and 2) the interactions between mycotoxins in vitro were evaluated using more appropriate mathematical models than the arithmetic definition of additivity.

Mathematical Models for Assessing Mycotoxin Interactions
In this paper, E will serve as an abbreviation for "effect" in equations. It is also assumed that effect is relative to maximal effect, i.e., percentage of cell viability suppression, where suppression is equal to difference between negative control (100% viability) and treated cells (100%-effect viability).

Simple Addition of Effects
The simplest method for estimating interactions between mycotoxins is the assumption of effect additivity known as arithmetic definition of additivity or response additivity (Equation (1)): (1) where E exp is the expected effect of combination of mycotoxin M 1 in dose D 1 and mycotoxin M 2 in dose D 2 , while E M1 and E M2 are the effects of single tested mycotoxins M 1 and M 2 in doses D 1 and D 2 , respectively. That simple addition of effect was applied by Šegvić Klarić et al. [34] for assessing the combination effect of beauvericin (BEA) and OTA using Equation (1) and observed synergistic effect for two combinations. Mathematically, this approach would be incorrect most of the time because the dose-effect curve is not linear. Using the data on cytotoxicity of OTA alone of the mentioned paper, it is easy to see that using this method we can prove that OTA applied in combination with itself at concentrations of 5 µM and 5 µM revealed an antagonistic effect; the expected cell viability would be around 20%, while the observed value for cell viability after treatment with 10 µM ochratoxin A was around 50% (Figure 1). Interestingly, despite an inaccurate estimation of expected effects, this model was widely applied; Alassane-Kpembi et al. [18] in their review cited 52 studies out of 83 that used this method.  [34]; arithmetic additivity calculation shows that upon treatment with 5 + 5 M of OTA expected viability is much lower than observed viability indicating antagonism (no copyright permission needed as we created this figure).

Factorial Analysis of Variance
This model uses simple 2-way ANOVA for modelling the detection of interactions between two mycotoxins (Equation 2): (2) Figure 1. Cytotoxicity of OTA (5 µM and 10 µM observed) on PK15 cells after 24 h of exposure [34]; arithmetic additivity calculation shows that upon treatment with 5 + 5 µM of OTA expected viability is much lower than observed viability indicating antagonism (no copyright permission needed as we created this figure).
Some studies presented in Table 1 [35,36] used simple addition of effects according to Weber et al. [24] who modified Equation (1) by subtracting the 100% (or 1) from the sum of the mean effects. Needless to say, the unexplained subtraction of 100% did not account for the non-linearity of the dose response curves.

Factorial Analysis of Variance
This model uses simple 2-way ANOVA for modelling the detection of interactions between two mycotoxins (Equation (2)): where E is the estimated effect, β 0 is the part of the effect achieved by negative control, β 1 /β 2 is the coefficient that increases effect for each increase in one unit of dose D 1 /D 2 of mycotoxin M 1 /M 2 and β 3 is the interaction term. Eight studies that have used this approach to define mycotoxin interactions were reviewed in detail by Alassane-Kpembi et al. [18]. If the interaction term was significantly (in a statistical manner) different than zero, it was concluded that an interaction between mycotoxins occurred. The main problem with this method is that ANOVA can be very misleading, similarly to the simple addition of effects method because ANOVA is based on linear modelling which is not useful for modelling nonlinear dose-effect curves [25]. This method was recently applied in only one study for testing the dual combination effects of ZEN and OTA or α-ZEL in HepG2 cells [37], as summarized in Table 1.

Bliss Independence Criterion
Bliss introduced this model in 1939 for predicting the proportion of animals that will die after combining two toxins under the assumption that there is no interaction between the toxins (i.e., they have completely different mechanisms of action or act in different compartments): where E exp is the expected effect of a combination of mycotoxin M 1 in dose D 1 and mycotoxin M 2 in dose D 2 , while E M1 and E M2 are the effects of single tested mycotoxins M 1 and M 2 in doses D 1 and D 2 , respectively [26], all effects need to be expressed as proportions ranging from 0 to 1 (Equation (3)).
Similarly to the simple addition of effects, Bliss can result in a detection of an interaction of some mycotoxin with itself but that is not possible in model validation since this would a priori violate the assumption of two toxins acting independently.
Several of the recent studies listed in Table 1 simultaneously used different mathematical models, e.g., response additivity and Bliss independence criterion [38,39] or Bliss independence and Loewe additivity [40] or Chou-Talalay method [39,41]. As expected, these studies obtained different conclusions on mycotoxin interactions depending on the mathematical models that have been applied. For example, Smit et al. [39] obtained a synergism of DON + ZEN at low and medium concentrations by both response additivity and Bliss independence model; while at high concentrations in combinations, an additive effect was obtained with Bliss independence model and antagonism by response additivity.

Loewe's Additivity Law
Loewe's additivity law (also called isobolografic method, concentration additivity or dose additivity) assumes that mycotoxins act within the same compartment on the same biological size by the same mechanism. The only difference is in their potency. This model is based on the dose equivalence principle and the sham combination principle; in short, every dose D 1 of mycotoxin M 1 gives an equal effect as D 2(1) of mycotoxin M 2 , and vice versa, and any D 2 (1) can be added to any other dose of D 2 to show the additive effect [27] as presented by Equation (4): where E is the effect, D 1 is the dose of mycotoxin M 1 , D 2 is the dose of mycotoxin M 2 , D 1(2) dose of mycotoxin M 1 that provokes same effect as D 2 dose of mycotoxin M 2 , D 2(1) dose of mycotoxin M 2 that provokes same effect as D 1 of mycotoxin M 1 . For additive effects, the following Equation (5) is valid: where D 1 and D 2 are the doses of mycotoxins M 1 and M 2 applied in combination, and D E1 and D E2 are the dose of mycotoxin M 1 and M 2 applied alone. All doses (D 1 +D 2 , D E1 or D E2 ) result with the same effect E. Additionally, Loewe's additivity law makes a larger number of assumptions; each mycotoxin in a mixture must have an equal maximum effect and all log(dose)-effect curves must be parallel and have constant relative potency [42,43], according to Equation (6): Finding two mycotoxins in a combination that fulfils all of these assumptions seems somewhat impossible. For example, apart from the Bliss independence criterion, Li et al. [44] also used this method (as a concentration addition model) to assess the nature of interaction between OTA and ZEN. Since their dose-effect curves did not meet all of the assumptions, it is easy to see that Equation (4), on which Loewe's additivity law is based, does not hold true when we assign the values EC 10 (OTA) = 0.8 µM and EC 10 (ZEN) = 11.84 µM [44], and try to apply the main principles of dose equivalence and sham combination of this model (Equations (7) and (8) This does not seem to be correct according to the dose-response curves for OTA (E (1.60 µM of OTA) ≈ 30%) and ZEN (E (23.68 µM of ZEN) ≈ 50%) presented in aforementioned article [44], which raises the question: can the observed synergies be trusted at all? Even though this model is mathematically valid, due to the excessive number of assumptions that need to be fulfilled, this model probably remains inapplicable for assessing combinations of mycotoxins [43].

Response Surface
Some authors expanded the Loewe's additivity law and Bliss independence criterion to the whole surface defined by all predicted additive concentration combinations (in all ratios, for all effects) [45,46] as presented in Table 1. In mycotoxicology, Assunção et al. [46] implemented model generalization built by Jonker et al. [28]. They estimated the deviation from Loewe's additivity law by Equation (9): where G is the deviation function defined separately for 4 models. If G = 0, then Equation (9) collapses to Equation (5), suggesting an additive effect. To test for synergism or antagonism G is substituted with (Equation (10)): G (z 1 , z 2 ) = a × z 1 × z 2 (10) where parameter a is less than zero for synergisms and greater than zero for antagonisms, z 1 and z 2 are relative contribution to toxicity, i.e., for z 1 as presented by Equation (11): Jonker et al. [28] also define more complicated interaction patterns between two toxins and with the inclusion of parameters b 1 for detection of dose ratio-dependent deviation (Equation (12)), and parameters b DL for the detection of dose level-dependent deviations (Equation (13)) from a non-interacting additive model: The procedure by Jonker et al. [28] suggests fitting all four models (defined by four deviation functions) and then choosing the best one to make conclusions about the nature of the interaction at different dose ratios or dose levels based on parameters a, b 1 , and b DL according to Table 1 of Jonker et al. [28].
This method provides more information than the other methods mentioned in this article, but it comes with a greater cost of the experiment since a checkerboard experimental design is needed, with dense concentration ranges in all combinations.

Highest Single Agent (HSA) Model
This model is also referred to as the Gaddums non-interaction [32], it defines the expected effect as the maximum of single mycotoxin effects (Equation (14)): where E exp is the expected effect of a combination of mycotoxin M 1 in dose D 1 and mycotoxin M 2 in dose D 2 , while E M1 and E M2 are the effects of single tested mycotoxins M 1 and M 2 in doses D 1 and D 2 , respectively. Because of underestimations of the expected combination effect, this model is not appropriate for detection of synergistic effects, except in cases: (i) where one compound is completely inactive at any concentration for the measured effect (which is rare in the field of mycotoxins); (ii) where a mycotoxin with maximal effect does not reach full effect (i.e., never suppresses viability to 0%). On the other hand, this method is useful for detecting antagonistic effects since observing a combination effect less than the maximal effect of a mycotoxin alone clearly demonstrates an interaction of antagonistic nature. However, underestimations of the expected combination effect can hide milder antagonistic effects. The great advantage of this model is the financial cost of the experiment: to prove an antagonistic effect, it is sufficient to test three concentrations, each mycotoxin alone and a combination of the mycotoxins. Another advantage is that this method is also independent of the mechanism of action, and it does not make any assumptions on the dose-effect curve. However, this simple approach has never been applied in mycotoxicology. Table 1. Interactions between mycotoxin combinations in vitro assessed by simple addition of effects, full factorial analysis, Bliss independence criterion, Loewe additivity law and response surface. Response additivity, CI RA ) and Bliss independence criterion (independent joint action, CI IjA ); IC 10 (1:1) and IC 30

Combination Index and Isobologram Analysis
Applying Loewe's additivity law or similar methods can allow researchers to use the Interaction/combination index which is based on Equation (5) for describing the nature of the combination effect (Equation (15)): where CI is the interaction/combination index: CI < 1 indicates synergism, CI = 1 indicates an additive effect and CI > 1 indicates an antagonism [29]. Isobologram analysis is just a "fancy" name for the graphical representation of the combination index for the same effect in different ratios of two mycotoxins. It is a simple plot with the dose/concentration of mycotoxin 1 on the x axis and the dose/concentration of mycotoxin 2 on the y axis. The characteristic line, isobole, connects the y intercept and x intercept which represents the doses needed for achieving a defined effect (i.e., 50%) for single acting mycotoxins. Plotting the dot with coordinates of doses in combination that achieve the same defined effect gives us a clue about the nature of the combination effect. All of the dots below the isobole indicate synergy, the dots above the isobole indicate antagonism, while the dots on the isobole indicate a possible additive effect [49]. The combination index and isobologram method were applied in 15 studies reviewed in Alassane-Kpembi et al. [18] and was the second most used method for assessing mycotoxin interactions and much more appropriate than the arithmetic definition of additivity or factorial design. The problems of not meeting the assumptions of Loewe's additivity law affect the combination index and isobologram. For example, if the two dose-response curves are not parallel, instead of one linear isobole, there will be two curvilinear isoboles around the former, linear one. The area between the two new curvilinear isoboles is not an area of synergy, nor is it an area of antagonism [43]. A recent study by Anastasiadi et al. [50] generalized the Loewe's model accounting for nonparallel dose-response curves. As a result, Equation (15) was expanded to Equation (16): where m 1 and m 2 are the slopes of the dose-response curves for mycotoxin 1 and mycotoxin 2.

Chou and Talalay's Median Effect Approach
Chou and Talalay developed a unified general theory for the Michaelis-Menten, Hill, Henderson-Hasselbalch, and Scatchard equations, mathematically presented by Equation (17): where E is the effect (between 0 and 1), D is the dose, D M is the median effective dose (i.e., EC 50 ) and m is a parameter for shape definition (if m < 1 dose-effect curve is hyperbolic, and if m ≥ 1 dose-effect curve is sigmoidal) [30]. Using Equation (17), it is possible to estimate the doses needed to achieve a particular effect which can be used in Equation (15) for the estimation of CI, which is then used for assessing the nature of the combination effect. Similarly to Loewe's additivity model, the isobologram can be constructed. The Chou-Talalay model combined with an isobologram has been applied in the majority of the recently published studies [39,[51][52][53][54][55][56][57][58][59][60][61][62][63][64][65] listed in Table 2. Its great advantage is the recent development of a method for the estimation of confidence intervals for the combination index which enables the application of statistics [66]. This method can easily be implemented using the web-based CalcuSyn software which automatically calculates dose-effect curves and combination indices.

MixLow Method
Compared to the Chou-Talalay method, the MixLow method (Table 2) used by Lin et al. [67] improves model fitting and removes bias by fitting the log-logistic curve without prior linearization, similarly to the CISNE method (discussed in Section 2.10). However, another improvement of the MixLow method is the inclusion of random effects in a model that can account for different batches (trays) in the experiment and fit the model for both toxins and combination simultaneously [31]. Mixed modelling enables a more precise estimation of the combination index's (CI, here called Loewe's index) and more reliable confidence intervals or standard errors by accounting for both the error of single applied mycotoxins and combinations.
The MixLow method comes with the mixlow R package, which also includes functions for straigthforward data import and minimal data preprocessing, especially if the pattern suggested on the tray is followed during experimental design [68].

CISNE (Code for the Identification of Synergism Numerically Efficient)
Even though Chou-Talalay's method exceeded two and a half thousand citations in relevant article databases, it does possess some technical problems in model fitting leading to bias inclusion in parameter estimation. By Chou-Talalay's protocol Equation (17) is rearranged and transformed to linear form (Equation (18)): where y is log[E/(1 − E)], the intercept is -m × log(D M ), the slope is m, and x is the log(D) of the linear equation form. Estimating slope and intercept by least squares fit, and calculating D m as presented by Equation (19): This leads to bias, along with the exclusion of data points with effects smaller than 0% or larger than 100% (i.e., stimulation) which could not be used in the logarithm on the left side of Equation (18). García-Fuente et al. [33] showed that these biases can lead to significant false positive or false negative errors, depending on the slope of the dose-response curve. They also found that fitting the same equation as a non-linear regression model estimates model parameters better and reduces the rate of false positives or negatives, especially when the slope (m) deviates from 1. This non-linear regression can be easily applied using the free CISNE software [69]. In contrast, it has not yet been applied in mycotoxicology combination testing.

Other methods
Most of the recent studies used mathematical modelling according to Bliss or/and Loewe (or some modified Loewe's method) for assessing the nature of the effect of combination of mycotoxins. However several in vitro studies assed mycotoxin combined effects comparing the effect of combination to the effect of single mycotoxin [70][71][72] or only to negative controls [73][74][75][76] without estimating the theoretical (expected) effect of the combination (Table 3). Conclusions based on those studies are unreliable because the question of the nature of interaction of combination has not even been asked in a scientific manner to get a clear and exact answer. For example, Smith et al. [75] did not define the nature DON + ZEN interaction in HepRG cells; since the cytotoxic effect of a single DON was similar to the effect of DON + ZEN, it was concluded that a combined effect could not be classified as antagonistic nor synergistic. Any conclusion about an antagonistic or synergistic effect should include the effect of ZEN too, since it is a part of the mycotoxin combination.

Conclusions
Some of the methods found in studies assessing the effects of mycotoxins combination have been incorrectly based on the assumption that mycotoxin dose-effect curves are linear (simple addition of effects, factorial analysis of variance). For that reason, many conclusions have been derived incorrectly in published articles or review articles based on published data. There are many articles reviewing methods and discussing the problem of the misuse of some method, but it seems that the problem persists. The only appropriate approach to assess the nature of an interaction is to correctly estimate the dose-effect curves of each mycotoxin and combination and apply a well-defined model (based on Bliss or Loewe's theory) with respecting the model's assumptions and fitting the model by a direct estimation of all model parameters from a nonlinear least squares fitting. Results should be presented in a simple and clearly defined way (i.e., isobologram or combination index) with some of the most expected (mean) values accompanied by uncertainty bounds, where a 95% confidence interval should have priority over the standard error due to asymmetrical distributions.
Improvements to the presented methods are continuously being made but are not readily applied in the field of mycotoxicology.