How to analyse and account for interactions in mixture toxicity with toxicokinetic-toxicodynamic models

and antagonism complexify the risk assessment of chemical


Introduction
Ecosystems are subject to many chemical inputs from agricultural, industrial and domestic sources, resulting in a wide range of mixture exposure scenarios for non-target organisms (Ballabio et al., 2018;Boxall, 2010).Chemical mixtures appear as the dominant scenario in ecosystems, and understanding the consequence of this is one of the major challenges in ecotoxicology (Escher et al., 2020).Over the past decades, the joint effects of chemicals have been studied with dose-response models.Typically, two basic models are used: the concentration addition (CA) model for chemicals with the same mode of action, and the independent action (IA) model for chemicals with different mode of action (Cassee et al., 1998;Eggen et al., 2004;Van Gestel et al., 2010).These two reference models have shown their ability to predict simple and complex mixture effects (Altenburger et al., 2000;Faust et al., 2001).However, deviations from predicted effects by those models were highlighted for many chemical mixtures (Cedergreen, 2014), and in different species (Gottardi and Cedergreen, 2019;Meled et al., 1998).This over-or under-estimation of the observed effects, commonly named antagonism and synergism respectively, may be due to one chemical affecting the effects of another or indeed through a more complex mechanism.These interactions add complexity to the risk assessment of mixture hazards, but it needs to be pointed out that they are identified as deviations from a null model: additivity.The assumptions underlying additivity can differ and so can the resulting additivity model (Schäfer and Piggott, 2018).The CA and IA model, with the addition of a new parameter and a statistical test, are commonly used to screen mixture effects and identify synergistic or antagonistic interactions of chemical mixtures when observations deviate from expected additivity (Spurgeon et al., 2010;Cedergreen, 2014).However, this framework presents two main weaknesses: first, these models are derived from single time point dose-response curves, and so ignore the time dimension while toxicity is a process in time (Baas et al., 2010).Second, when joint effects deviate from either CA and IA (or both), it is concluded interaction occurs (Cedergreen, 2014), but within the analyses, neither model can provide an explanation for the causation of any interaction, nor predict such joint effects for different scenarios (e.g., time variable concentration) than those used to calibrate the model.
A promising alternative to the classic ad hoc approach is the use of toxicokinetic-toxicodynamic (TKTD) models which simulate the timecourse of processes leading to toxicity (Ashauer and Escher, 2010;Jager et al., 2011).In contrast to dose-response models, TKTD models have meaningful biological parameters and can deal with concentrations that vary over time (Ashauer and Escher, 2010;Péry et al., 2001), or when chemicals are applied sequentially in differing order.The latter is important because different exposure sequences of the same chemicals can lead to differences in effects (Ashauer et al., 2017).Among TKTD models, the General Unified Threshold model for Survival (GUTS) framework has been recently extended to mixture toxicity based on the concept of damage dynamics which links the internal concentration to the effects (Bart et al., 2021).In the same way as the classic mixture models, the GUTS mixture models may find deviations from the predicted mixture toxicity and interpret these as synergistic or antagonistic effects (Bart et al., 2021).However, when such interactions are observed, GUTS models (and in general all TKTD models) have the potential to provide explanations of the synergism or antagonism, and predictions thereafter.The interaction between chemicals leading to synergistic or antagonistic effects may involve toxicokinetic and/or toxicodynamic processes (Spurgeon et al., 2010), and when the mechanisms of these interactions are understood, they can be incorporated into GUTS models.Such next generation mixture effect models might accurately predict mixture effects with interactions.In this way, interactions will be recognised as artefacts of underdeveloped null models (Schäfer and Piggott, 2018).In a recent cutting-edge study, Cedergreen et al. (2017) used a GUTS model to simulate the synergistic effect between azole fungicides and cypermethrin on Daphnia magna.These authors extended GUTS by incorporating the ability of the azole fungicides to decrease the biotransformation rate of cypermethrin, leading to an increase of the internal concentration of cypermethrin and therefore a synergistic effect of those two pesticides.The authors used a full GUTS model including body residues to demonstrate how TKTD models can account for interactions.Unfortunately, body residues data are scarce and the use of GUTS reduced (RED) models, which link external concentration to the effect, is more reasonable as a first approach for screening of many mixtures for potential interactions and their predictions thereafter.Baas et al. (2007) used a GUTS-RED model with an interaction factor, a statistical interaction term that allows comparison of the model with and without an interaction term using a likelihood-ratio test (models are nested).The disadvantage of this approach was that the interaction factor was a descriptive term without biological meaning (so cannot be used as a metric for risk assessment nor used for interpretation of toxicological processes) and was applied on the hazard rate as calculated in the stochastic death (SD) approach, and so excluded the potential use of the individual tolerance (IT) approach of the GUTS framework.A trade-off between a full GUTS model and a GUTS-RED model with an additional ad hoc parameter is necessary to improve the risk assessment of mixture hazards and improve our understanding of synergistic and antagonistic interactions.Up to now a GUTS model for mixture interactions was missing.As a result, the process-based identification and interpretation of mixture interactions was not possible, and we did not have a method to predict effects of interaction mixtures for time-variable exposures with changing mixture compositions.However, this is necessary for realistic environmental risk assessments as such exposure scenarios are the reality.
We here aimed at extending the approach proposed by Baas et al., 2007 through the use of a more mechanistic module that is i) able to account for synergism and antagonism, and ii) working for both, the SD and the IT approach.In a previous paper (Bart et al., 2021), we fully described the extension of GUTS models to assess mixture effects with dedicated experiments and data extracted from the literature.We here used this framework to implement an interaction module allowing simulation and prediction of binary synergistic or antagonistic mixture effects.We tested the interaction module with dedicated experiments with two different species, and on data extracted from literature with two other species exposed to mixtures reported to be synergistic or antagonistic.

Species, chemicals and experimental design
Enchytraeus crypticus (Enchytraeidae, Oligochaeta, Annelida) were originally sourced from the laboratory of the Department of Ecological Science, Vrije Universiteit, Amsterdam, The Netherlands, and were maintained in culture at the UK Centre for Ecology and Hydrology, Wallingford (UK).The cultures were kept on 1.5 % agar, made up with artificial fresh water with the following concentrations: CaSO 4 0.12 g/L, MgSO 4 0.246 g/l, NaHCO 3 0.192 g/L and KCl 0.008 g/L.These were maintained in a climate room at a constant temperature (20 ± 1 °C) in darkness.Diet was prepared using oatmeal (150 g), dried baker's yeast (50 g), baby formula milk powder (50 g), and fish oil (50 ml), which was heated with 400 ml milk until boiling.Cultures were fed ad libitum.For all experiments, adult individuals were exposed for 96 h at 15 °C in darkness, in 1 ml of artificial fresh water (Roembke and Knacker, 1989) in 24 well plates, with one individual per well, so 24 individuals monitored per treatment.Survival was monitored at 3, 6, 24, 48, 72 and 96 h.Individuals were classified as dead if they did not respond to repeated touching with a pin.
Mamestra brassicae larvae were from a culture maintained at the UK Centre for Ecology and Hydrology.The cultures were kept as larvae in rearing boxes containing artificial diet (see recipe in the supplementary material); and as adults in flight cages on a 10 % honey water solution.Both the adult and larvae culture were maintained in a controlled climate room at a constant temperature of 20 ± 1 °C and a 16:8 light:dark pattern.For all experiments, 4th instar larvae were exposed to 1 ml of pesticide stock solution spiked at the concentration required to give the appropriate exposure level within the artificial diet (see supplementary material) and kept at 20 ± 1 °C in darkness for 96 h.Experiments used 12 well plates, with one individual larvae per well.Survival was recorded at 24, 48, 72 and 96 h.Larvae were classified as dead if they did not respond after 5 s of stimulation on the second segment behind the head with a fine paintbrush.
For all tested binary mixtures, we first exposed the individuals to chemical A separately (single), chemical B separately (single), and then a mixture of chemical A and B. The single exposures were designed to cover the toxic effect, from no effect to 100 % mortality (based on a pilot study), to calibrate the GUTS mixture models.The mixture exposures were design to cover the potential interactions between the two chemicals, and estimate the parameters zI and bI as robustly as possible, and their uncertainties (if interactions occur).For that purpose, we selected different ratios of the two chemicals (based on a pilot study).In addition, we selected only concentrations that were in the single exposures to reduce variability in the experimental set-up and in the model calibration.A summary of the concentrations applied is presented in Table S1.
In addition to the bespoke experimental studies, we also extracted previously published data from the literature to assess if the new interaction module is able to account for interactions occurring within further species and substances.We used data from Robinson et al. (2017) who exposed Bombus terrestris to single and mixture treatments of clothianidin (CLO) and propiconazole (PRP), and CLO and dimethoate (DIM), and reported synergism and antagonism interactions for those mixtures, respectively.We also used data from Loureiro et al. (2010) who exposed Daphnia magna to single and mixture treatments of imidacloprid (IMI) and thiacloprid (THI) and reported a synergistic effect.Finally, we tested our modelling approach with data from Cedergreen et al. (2017) who modelled the synergistic effects between CYP and propiconazole (PRP) on Daphnia magna with a GUTS model including body residues.We only used the external concentrations and survival data over time from the latter study.

The GUTS framework and the extension for mixture effects
The GUTS models simulate the time-course of processes leading to the death of an organism (Jager et al., 2011;Jager and Ashauer, 2018).They account for the accumulation of, and recovery from, damage (the toxicodynamics, TD), which results due to the bioaccumulation, distribution, biotransformation and elimination of the chemicals in the organisms (the toxicokinetics, TK).In the absence of body residues data, the TK and TD are combined.The models using this approach are commonly named the GUTS reduced (RED) models.In GUTS-RED models, the dominant rate constant 'k d ' describes the dynamics of the 'scaled' damage and will represent the one-compartment approximation of the "true" twocompartment behaviour (TK and damage dynamics).To describe the death mechanism related to the damage, two causations of the process affecting survival are formalised: the stochastic death (SD) and individual tolerance (IT) approach.The SD approach assumes that individuals are identical and have a probability to die upon chemical stress, which increases with increasing damage once some threshold damage has been exceeded.The IT approach assumes that individuals have differences in their sensitivity to chemical stress, and when the damage exceeds an individual's threshold it dies instantly.
The GUTS framework was extended for mixture effects, and is fully presented in Bart et al. (2021).Two models were developed: the GUTS-RED damage addition (DA) to account for the mixture effects of chemicals with the same mode of action, and the GUTS-RED independent action (IA) to account for chemicals with different mode of action.To summarize, in both models, each chemical keeps its one dominant rate constant.Then, the DA model considers chemicals leading to the same form of damage and therefore it is based on the addition of their scaled damage with a weight factor accounting for the difference in toxic potency.Therefore, the chemicals share the related mortality parameter.In the IA model, each chemical keeps its own GUTS parameters and the overall survival probability is obtained by multiplying the survival probability due to each chemical in the mixture.The DA model can be considered only when chemicals can share their related mortality parameter.This allows the simultaneous fits of the single chemical exposures to predict the mixture effect thereafter.Both mixture models can be used with the SD, and IT approach.The equations of the GUTS mixture model are presented in the supplementary material.

The interaction module, hypothesis and formulation
In situations where interaction occurs between chemicals, the GUTS mixture model will over-or under-estimate the survival probability over time.Therefore, the model needs an extension to be able to account for interactions and provide a more accurate prediction of the survival over time.
An extension of the model to account for interaction between chemicals must i) be able to change the final effect on survival (increase or decrease), and ii) reflect the chemical/biological/biochemical processes of the interaction.From a modelling perspective, there are two possibilities to change the output value of the GUTS-RED mixture model (i.e., the survival probability): a change in parameter values, which leads to a change in model variables, or a direct change in variables.
Change in parameter values (model parameters remain constant through a simulation).There are five parameters in the GUTS-RED models: the background hazard rate h b that we did not consider changing because it is not related to a toxicological process, the dominant rate constant k d , the median of the threshold distribution m w , the killing rate b w for the SD approach and the spread factor of the threshold distribution F s for the IT approach.
The dominant rate constant (k d ) governs the time needed for the damage density to reach steady state with the external concentration.A change in the k d value would assume that the interaction leads to a change in the toxicokinetics and/or toxicodynamics.An increase in k d would reflect a synergistic interaction, which will lead to steady state between the external concentration and the scaled damage being reached more rapidly.However, increasing this parameter would not be able to increase the damages beyond this steady state.In other words, the effects could appear faster in the simulation, but could not scale beyond the default model without interaction.Therefore, we discarded a change in k d .
The median of the threshold accounting for the effect (m w ) is part of the death mechanism and has the units of an external concentration in the GUTS-RED model.A lower threshold would lead to stronger effects, but the concept of a threshold for effect is rather vague mechanistically and it is hard to link this parameter to any biochemical process or measurement to support this hypothesis.In addition, m w is different between the SD and IT approach, so the change would be different between the two approaches.Therefore, we also discarded a change in m w as a solution.
The F s and b w parameters are part of the death mechanism of the IT and SD approach respectively, which would mean that any changes in these parameters will be specific to the chosen approach.In addition, a change in F s would mean that some individuals will become more sensitive, but others more resistant at the same time.This possibility is not logical in terms of observed synergistic and antagonistic effect.Therefore, we also discard a change in F s or b w .
Change in variables (variable values change through the course of a model simulation).There are two variables in the GUTS-RED model: the scaled damage D w (t), and the Survival probability S(t).Any change in S would be mechanistically meaningless because we would deal only with the final effect, while the survival probability is intended to be the result of toxicological mechanisms, as described in the GUTS framework.Therefore, we assume D w as the best candidate to account for interactions in GUTS-RED models.The interactions between chemicals can be driven by effects on seven processes (or groups of processes) that can drive toxicity of a chemical: the bioavailability, uptake, internal transportation, metabolization, excretion, binding at the target site and the Adverse Outcome Pathway (AOP) at the individual level (i.e., excluding the population scale and higher tier processes) (Groh et al., 2015;Spurgeon et al., 2020) (Fig. 1).We assume that any change in those processes, will lead to a change in the scaled damage (Fig. 1), and therefore to the simulated effect as described by the death mechanism in the GUTS-RED model (SD or IT).Therefore, we choose to change the scaled damage D w variable as the proxy for the interaction.
2.2.3.How to characterize the modification of the scaled damage dynamics due to the interaction?
The GUTS reduced model combines the toxicokinetics and the damage dynamics (i.e., the toxicodynamics), with an unknown partition coefficient K Dw between the external concentration and the damage level.Any change in toxicokinetics and/or toxicodynamics due to the changes in processes driving the toxicity should lead to a change in K Dw accordingly.For an antagonistic interaction, K Dw will decrease leading to a decrease in the scaled damage, and for a synergistic interaction, K Dw will increase leading to an increase in the scaled damage.In a GUTS-RED model, K Dw is unknown (and effectively scaled out), therefore, a factor must be applied directly to the scaled damage.We propose the introduction of an interaction factor 'I' which will increase or decrease the scaled damage, reflecting this change in K Dw .Considering a mixture of two chemicals A and B, and that the toxicokinetic and/or the toxicodynamic of the chemical A is affected by the chemical B, it leads to the following equation: where A and B refer to two chemicals, D w (t) is the scaled damage [with an external concentration unit, e.g., mg L −1 or μM], C w (t) is the external concentration [e.g.mg L −1 or μM], k d is the dominant constant rate [d −1 ] and I is the interaction factor [unitless].We assume that I B will depend on the characteristic of the chemical B which creates the interaction.It has been suggested that interaction is concentration dependent and in most situations a threshold above which the interaction occurs (Cedergreen, 2014).Therefore, we assumed that the compound creating the interaction (B) is characterized by i) a concentration threshold from which the interaction starts (i.e., a certain quantity of chemicals is required to start creating an interaction), and ii) an interaction strength which is multiplied with the concentration above the threshold.We have translated this into Eq.( 2) with two new parameters that describe I B as follows: where b I_B is the strength of the interaction created by the compound B [1/ external concentration, e.g., 1/mg L −1 ], Cw B (t) is the external concentration of the interacting compound (B), and z I_B is the threshold from which the interaction starts [with an external concentration unit, e.g., mg/L −1 ].
In a situation where C wB (t) < z I_B , I B (t) = 1, and will not change the scaled damage of chemical A. In the case of an interaction where C wB (t) > z I_B ; if b I_B > 0 ➔ I B (t) > 1 ➔ I B (t) will increase the scaled damage, reflecting the situation currently named synergism.And, if b I_B < 0 ➔ 0 < I B (t) < 1 ➔ I B (t) will decrease the scaled damage reflecting the situation currently named antagonism.The value of the b I_B parameter will reflect the strength of the interaction (for synergism, with b I_B positive, or antagonism, with b I_B negative).The Fig. 2 presents how this parameter will change the scaled damage and how it can bring it beyond the non-interaction scaled damage steady state for a synergistic mixture.
The interaction in a binary mixture can come from chemical A acting on chemical B, or B on A, or a double interaction.For each tested binary mixture and species, we tested these three situations.In the situation of a double interaction, the b I_A and b I_B parameters must be constrained to be both positive or both negative to avoid a situation where the two could cancel each other out.

Workflow and assessment of the interactions
The first step was to predict the mixture effect based on the calibration on single exposures with the same workflow presented in Bart et al. (2021).We first assessed if the GUTS-RED DA model is suitable by testing if the two chemicals can share the mortality related parameters and assessing the simultaneous fits.Then, we simulated predictions of the mixture effects with all models (GUTS-RED-SD IA, GUTS-RED-IT IA, and, if applicable, GUTS-RED-SD DA and GUTS-RED-IT DA).
The second step was the fit of the interaction module with the three different configurations (A modifies B, B modifies A, or double interaction) on the mixture data with all models.In this step, we kept fixed those GUTS parameters that were previously fitted on the single exposures (used to provide the predictions of the mixture effects without interaction).This procedure ensures that the model does not compensate for a poor representation of the mixture treatment data by changing the modelled effects of chemicals in single exposures.The final step was to assess whether the interaction module improves the simulation of the mixture data, and so if an interaction occurs.The model without the interaction is a simplification of the model with the interaction module, they are thus nested.This allows us to statistically test if the interaction module significantly improves the fit of the mixture data using a likelihood-ratio test.In other words, the test answers the question "Does the addition of the interaction module provide a better simulation of the data than the model without interaction?".If the interaction module does provide a significantly better fit of the mixture data, it is concluded that an interaction occurs, and if the b I parameter is positive, it is concluded synergism occurs (or potentiation if one of the chemicals does not have effects in single exposure); if the b I parameter is negative, it is concluded antagonism occurs.In the same line, the model with simple interactions (A on B, or B on A) is a simplification of the model with double interaction and can also be compared, answering the question 'Does the double interaction provide a better simulation of the data than the simple interaction?'All together these steps allow to statistically test if i) an interaction occurs in the considered mixture and ii) which chemical is most likely to cause the interaction (based on the result of the likelihood-ratio test), and finally the positive or negative value of the b I parameter indicates whether the data is best modelled with a synergistic or antagonistic interaction.

Model calibration and statistics
The Wilson score interval (Brown et al., 2001) was used to express the uncertainty in the survival data and allows plotting confidence intervals on each data point.All calculations were performed in Matlab 2020 using the BYOM modelling platform (www.debtox.info/byom.html).The optimisation of the parameter values was performed with the parameter space explorer (Jager, 2020).This algorithm is optimised for GUTS analyses, and combines grid search, a genetic algorithm, and likelihood profiling to calculate the confidence intervals of the parameter values.To produce confidence intervals on the model curve, a sample from the parameter space explorer was used.This ensures parameter covariance is taken into account when propagating parametric uncertainty to predictions with confidence intervals.For all calibrations, the background mortality (h b ) was fitted to the survival in the control treatment and kept fixed while fitting the toxicity parameters.For all simulations, we provided the minus log-likelihood (MLL) and r2 values to assess the goodness of fit.Fig. 2. Illustration of the change in scaled damage, and survival probability thereafter, due to the interaction between chemical A and B. We simulated a theoretical exposure of a hypothetical species to chemical A at 15 mg/L and B at 10 mg/L, and in which only the chemical A had an effect in single exposure, and the chemical B was the chemical creating the interaction.The left side is for the stochastic death approach (SD) with k dA = 1.5 (d −1 ), m wA = 5 (e.g., mg L −1 ) and b wA = 0.1 (e.g., 1/mg L −1 ), and the right side with the individual tolerance approach (IT) with k dA = 1 (d −1 ), m wA = 7.5 (e.g., mg L −1 , F sA = 2 (−).For both approaches h b = 0 (d −1 ).The threshold for interaction z I_B was fixed at 0 (e.g., mg L −1 ) and b I_B (e.g., 1/mg L −1 ) was fixed at 0 (no interaction, black line), or 0.05 (cyan line) and 0.2 (red line) for synergism, and − 0.05 (green line) and −0.2 (blue line) for antagonism.

Ability of the interaction module to detect and simulate synergistic and antagonist mixtures
Our mechanistic modelling approach was able to detect interactions that were previously identified with the classic approach (i.e., based on dose-response model).For instance, Loureiro et al. (2010), concluded on a synergistic interaction for D. magna exposed to imidacloprid and thiacloprid in mixture.The interaction module coupled with the GUTS mixture model significantly improved the fits, allowing to conclude that there was a synergistic interaction in this case (Table 1, S15).Similarly, our modelling approach identified an antagonistic interaction for B. terrestris exposed to clothianidin and dimethoate in mixture (Table 1, S15), and a synergistic interaction for the same species exposed to clothianidin and propiconazole in mixture (Tables 1, S15), both previously identified with the classic approach in Robinson et al. (2017).For the latter studies, the authors did not conclude on the synergistic interaction when using a reduced GUTS model with the addition of an ad hoc interaction factor (as presented in Baas et al. (2007), while our model did.These results showed that our mechanistic interaction module coupled with the GUTS mixture model is able to detect interaction, like the classic approach, whilst additionally accounting for the entire exposure time and all observations at intermediate observation times, thereby, making the conclusion more robust.Moreover, contrary to the classic approach, our toxicokinetics-toxicodynamics model was able to simulate the time-course of the mixture effects with interactions (i.e., for the above discussed mixtures and species, Fig. 3, Fig. S30, Figs.S36-S39, Figs.S40-43).Therefore, the model can be used for the prediction of synergistic or antagonistic effects resulting from untested exposure scenarios such as longer exposures, time-variable concentrations, or sequential exposures.
Generating data with bespoke experimental designs for model and hypothesis testing enabled us to test temporal aspects of mixture effects on E crypticus and M. brassicae, and to assess if interactions occur.The prediction of survival of E. crypticus exposed to the mixture azoxystrobin and prochloraz with the GUTS mixture model greatly underestimated the observed effects (Figs. 4 left panel).The addition of the interaction module, calibrated on the mixture data, significantly improved the simulation of the mixture effects (whatever the model used; i.e., SD or IT approach with the IA and DA model, Table S15, Fig. 4; S11-S13 right panel), highlighting a synergistic interaction for this mixture and species combination (Table 1, Table S15).This result is in line with previous studies which found synergistic effects with this mixture in aquatic invertebrates (Fu et al., 2018).For M. brassicae exposed to cypermethrin and chlorothalonil in mixture, we identified for the first time a synergistic interaction (Figs. 5 and S28, Table 1, Table S15).For the mixture imidacloprid and prochloraz in E. crypticus, and cypermethrin and prochloraz in M. brassicae, the modelling approach did not converge on the same conclusion for the stochastic death or the individual tolerance approach (Table 1, S15).Therefore, we cannot firmly conclude on a synergistic interaction for these two cases, and further investigation would be required.
For M. brassicae exposed to cypermethrin and chlorpyrifos in mixture (Table 1, Table S15, Fig. S4) an antagonistic interaction was identified.Antagonism between these two insecticides was previously

Table 1
Summary of the nature of interactions identified, with toxicokinetic-toxicodynamic modelling, in the joint effects of binary combinations of chemicals in four species.The full table with all model configurations, and log-likelihood ratio test results, is presented in the supplementary material (Table S15).identified in the whitefly species Bemisia tabaci (Ahmad, 2007).Finally, for M. brassicae exposed to cypermethrin and azoxystrobin in mixture, where only cypermethrin had an effect in single exposure, the interaction module did not improve the fit.The GUTS-RED mixture model alone already provided an accurate prediction of the effects (Fig. 6), leading to the conclusion that no interaction occurs between those two chemicals in this species, so the mixture is likely additive.

New insights in detecting synergistic and antagonistic mixtures
The interaction module offers the possibility for hypothesis testing by comparing model fits with different configurations (A on B, B on A, or double interaction).It provides not only the positive or negative value of b I , indicating the orientation of the interaction (synergy or antagonism), but also indication about which of the chemicals has the dominant responsibility in the interaction.This information can help in identifying causal relationships and improve understanding of interaction mechanisms and pathways, either from literature or through experimental studies.Below, we discuss the plausibility of the interactions identified by the model, as a further step toward testable mechanistic hypotheses.
Considering the experiment with E crypticus and the mixture azoxystrobin and prochloraz, the likelihood ratio test highlighted that a simple interaction with prochloraz as the synergist (b I positive, Table 2, i.e., responsible for the increase in azoxystrobin effects) is the best configuration to describe the data (Table S15).This result is in line with what is known about prochloraz.Indeed, prochloraz, like other azole fungicides, has been shown to act as synergist in a range of studies, enhancing the toxicity of pyrethroid insecticides (Kretschmann et al., 2015), strobilurin fungicides (e.g., azoxystrobin, Fu et al., 2018), or neonicotinoids (Robinson et al., 2017), and in a range of organisms (Dalhoff et al., 2016;Kuhlmann et al., 2019).The mixture prochloraz and azoxystrobin in E. crypticus can be added to this list of cases.Potentially the mixture of cypermethrin and prochloraz in M. brassicae follows the same pattern, but this would need to be confirmed, because not 0.5 mg kg -1 CYP + 1 mg kg -1 AZX 2.5 mg kg -1 CYP + 20 mg kg -1 AZX 10 mg kg -1 CYP + 20 mg kg -1 AZX 10 mg kg -1 CYP + 50 mg kg -1 AZX 50 mg kg -1 CYP + 20 mg kg -1 AZX 50 mg kg -1 CYP + 50 mg kg -1 AZX 0.5 mg kg -1 CYP + 1 mg kg -1 AZX 2.5 mg kg -1 CYP + 20 mg kg -1 AZX 10 mg kg -1 CYP + 20 mg kg -1 AZX 10 mg kg -1 CYP + 50 mg kg -1 AZX 50 mg kg -1 CYP + 20 mg kg -1 AZX 50 mg kg -1 CYP + 50 mg kg -1 AZX PredicƟon (without interacƟon) GUTS-RED Independent acƟon (IA), stochasƟc death model GUTS-RED Independent acƟon (IA), individual tolerance model Fig. 6.Observed and predicted survival over time of Mamestra brassicae exposed to cypermethrin and azoxystrobin in mixture.The predictions are simulated with the GUTS-RED-SD independent action model calibrated on single exposures, with the stochastic death approach (first row of plots) or the individual tolerance approach (second row of plots).The confidence intervals of the prediction (green shadows) propagate the uncertainty from the model calibration on the single exposures to the prediction of the mixture effects.The dashed lines show the background mortality, fitted to the control treatment.

Table 2
Parameter values of the interaction module for each species and mixture tested.z I is the threshold for interaction and b I is the strength of the interaction, positive for synergism, and negative for antagonism.The full table including all model configurations is presented in the supplementary material (Table S14).all models reach that conclusion (Table S15).In M. brassicae, chlorothalonil was identified as a synergist in mixture with cypermethrin, by all models (Tables 1, S15).The scientific literature provides some mechanistic support for this finding.For instance, cypermethrin detoxification appears, at least somewhat, dependent on glutathione conjugation, and it has been shown that Glutathione S-Transferase (GST) activities are up-regulated in pyrethroid-exposed insects (Nielsen et al., 2000;Yan et al., 2013;Zhang et al., 2013).Critically, chlorothalonil is known to rapidly form substituted chlorothalonil-glutathione derivatives, thereby depleting glutathione reserves (Arvanites and Boerth, 2001;Tillman et al., 1973).Such depletion of this potential conjugate may inhibit cypermethrin detoxification and explain the synergism observed in M. brassicae.Another hypothesis to explain cypermethrin-chlorothalonil synergism in M. brassicae is that chlorothalonil induces behavioral and/or physiological changes that increase cypermethrin uptake.Indeed, cypermethrin has been associated with a significant repellency in adult Apis mellifera leading to a decrease in the amount of sucrose consumed, but this change in feeding was dramatically reduced with chlorothalonil co-exposure (Thompson and Wilkins, 2003).A similar mechanism in M. brassicae, exposed through feeding (and their body were also in contact with the diet), could account for the synergy.
Another interesting insight provided by the modelling, is for the antagonistic mixture clothianidin and dimethoate in B terrestris, for which it is suggested the interaction is due to clothianidin whatever the GUTS model used (damage addition and independent action model with the stochastic death or individual tolerance approach, Table S15).Several studies have reported that neonicotinoid exposure leads to increased expression of cytochrome P450 (CYP450) genes (Yang et al., 2018), and increased CYP450 expression is often linked to neonicotinoid resistance (Puinean et al., 2010).Both, the xenobiotic metabolism regulating receptors and the CYP450s they regulate are promiscuous, such that the same receptor and CYP450 can have multiple ligands and substrates respectively (Fujita, 2004;Li and Wang, 2010).The need for this redundancy means two toxicants metabolised by the same CYP450 may not have identical binding interactions with the same xenobiotic receptor (Kojima et al., 2011), leading to differing induction of the CYP450 (or CYP450s) responsible for their metabolism.Therefore, it is plausible that clothianidin induces expression of a CYP450 that, while capable of metabolising dimethoate, is not induced (or is only weakly induced) by dimethoate itself, leading to notably increased breakdown of dimethoate in cases of clothianidin-dimethoate coexposure.Alternatively, it is possible that clothianidin reduces expression of a CYP450 that acts to activate dimethoate, leading to reduced dimethoate toxicity and antagonism.A further hypothesis that may explain this antagonistic interaction relates to the disruption of microbial content, a component of an organism's biology known to have an important role in toxicological response (Jin et al., 2015;Duperron et al., 2020).Clothianidin is known to alter the gut microbiome of non-target species (Khoury et al., 2021), and the toxicity of dimethoate is largely dependent on its metabolic desulfuration into an active dimethoate-oxon form (Buratti and Testai, 2007).Given that dimethoate can be biodegraded by bacteria (Ishag et al., 2016), it is plausible that species present in the gut microbiome also metabolise dimethoate into its oxon form, as is the case with other organophosphates (Daisley et al., 2018).The clothianidin-induced changes to the microbiome composition of the B. terrestris gut could reduce microbialinduced dimethoate activation, thus lowering its overall toxicity, leading to clothianidin-dimethoate antagonism.
The cypermethrin-chlorpyrifos antagonism we observed in M. brassicae has been also reported in the lepidopteran species Helicoverpa armigera (Ahmad, 2004), in the silverleaf whitefly Bemisia tabaci (Ahmad, 2007), and for some mixtures of these compounds in the killifish Jenynsia multidentate (Bonansea et al., 2016).To our knowledge, this antagonistic interaction has not been mechanistically explained but the existing literature does offer plausible hypotheses.Our modelling approach indicated cypermethrin or chlorpyrifos could be responsible for the antagonistic interaction in M. brassicae.Like for clothianidin-dimethoate antagonism in B. terrestris, it is plausible that cypermethrin induces expression of a CYP450 that, while capable of metabolising chlorpyrifos, is not induced (or is only weakly induced) by chlorpyrifos itself.Indeed, both cypermethrin (Baek et al., 2010) and chlorpyrifos (Roh et al., 2014) are associated with the up-regulation of specific CYP450s.Again, as for the combination of clothianidin and dimethoate discussed above, an alternative explanation is that cypermethrin inhibits the activation of chlorpyrifos to chlorpyrifosoxon.Likewise, cypermethrin and/or chlorpyrifos exposure may alter the gut microbiome, leading to altered metabolism and/or uptake of the second compound.Indeed, both pyrethroids and organophosphates have been found to impact microbiome compositions (Giambò et al., 2021;Krishnaswamy et al., 2021), and the differential microbiome modulation by pesticides does raise plausible hypotheses for chemical interactions (Jin et al., 2015;Liu et al., 2022).For instance, the toxicity of chlorpyrifos is known to be largely dependent on its metabolic desulfuration into an active chlorpyrifos-oxon form (Eyer et al., 2009), which is carried out, at least in part, by microbiota inhabiting the gut (Daisley et al., 2018).Finally, there is evidence that esterases play important roles breaking down pyrethroids (Bhatt et al., 2020), including an isomer of cypermethrin in lepidopterans (Bai et al., 2019).Along with sodium channel mutations, constitutively increased esterase expression is associated with pyrethroid resistance (Bhatt et al., 2020) and pyrethroid exposure is associated with increased esterase activity and up-regulation of esterase gene expression (Bilal et al., 2018;Feng et al., 2018;Leong et al., 2018).These esterases bind to, and have their activity inhibited by, the chlorpyrifos-oxon (Chanda et al., 1997;Sanchez-Hernandez et al., 2018;Vejares et al., 2010), irreversibly binding with a 1:1 ratio to organophosphates (Maxwell, 1992).This has led to the suggestion that binding of the chlorpyrifos-oxon to non-neural esterases leads to stoichiometric detoxification, as it prevents a given oxon from interacting with the neural AChE target (Vejares et al., 2010).Therefore, a cypermethrin-chlorpyrifos antagonistic interaction may occur in M. brassicae if cypermethrin exposure induces notable up-regulation of esterase expression, reducing chlorpyrifos toxicity by increasing the protective capacity of non-neural esterases.

Advantages and limitations
The development of a full GUTS mixture model including body residues measurement (of parent compounds and metabolites) is the most sophisticated option to accurately predict mixture effects over time with synergistic or antagonistic interaction.This strategy was followed by Cedergreen et al. (2017), where a full GUTS model was used to simulate the synergistic interaction between cypermethrin and propiconazole on D. magna.However, in a regulatory risk assessment context, the first step is the detection of such a mixture interaction, and a detailed study including body residues is frequently not feasible for the thousands of chemicals in circulation for practical reasons.Therefore, risk assessors need a quick and efficient approach for assessment of mixture toxicity.We tested our model using data presented in Cedergreen et al. (2017), and our simplified model was able to identify the interaction (Table S15), and suggested propiconazole was responsible for the interaction.This is in line with what Cedergreen et al. (2017) described mechanistically in their model, highlighting that our simpler approach, without explicit consideration of body residues, can be sufficient to identify such an interaction.In addition, our model was able to simulate the synergistic effects on survival over time (Fig. 7), although the model of Cedergreen et al. (2017) was slightly more accurate than our own (Fig. 3B in Cedergreen et al. (2017)).However, this better accuracy of the Cedergreen et al. (2017) approach was achieved at a significant cost in terms of the extra information required for parameterisation.Indeed, these authors needed for input, the external concentration and internal concentration of cypermethrin and propiconazole, the cypermethrin metabolization rate, and the survival over time, to allow the 14 parameters of their model to be estimated.In contrast, our model only used the external concentration and survival data over time, and is composed of 9 parameters.The difference in the accuracy of the simulation is probably due to the nature of the interaction.For this mixture, propiconazole reduces the biotransformation rate of cypermethrin by inhibiting cytochrome P450 activity, and this enzymatic phenomenon is better described with a sigmoid curve based on internal concentrations of propiconazole.Our interaction module uses a simple linear relationship based on the external concentration of propiconazole (above the threshold for interaction, z I ).We made this simplification, not only to reduce the number of parameters and the data required for calibration, but also to offer the opportunity to identify and model other types of interactions, and not only the specific case of decrease of biotransformation of a co-occurring chemical.Our model is, therefore, both simpler and more generic, whereas the model of Cedergreen et al. (2017) is more mechanistic and specific.In this respect, the two approaches are highly complementary.
This focus of the current modelling approach on simplicity and flexibility, is made to offer the possibility to screen many binary mixtures relatively quickly to identify interactions and to allow as many interaction types as possible to be indicated.Our modelling approach only requires the exposure concentrations and the observed effects over time as input.However, more complex situations can occur, for example mixtures which can be synergistic at low concentration, and antagonistic at higher concentration (Cedergreen, 2014).Such complex situations will not be identified by the current interaction module, which is only able to identify solely synergistic or antagonistic cases working in one direction.M. brassicae exposed to cypermethrin and prochloraz is a good example of such a more complex scenario (Figs.S24, S25).This mixture was suspected to be synergistic, for reasons discussed previously regarding the chemicals involved.The prediction which underestimated the effects could be seen as evidence for synergism, but this was not observed in all combinations of concentrations.The interaction module improved the fit only with the GUTS-RED-IT IA model (Table S15, Fig. S25), and the best simulation of this dataset was from the GUTS-RED-SD IA model without interaction (MLL = 123.8,Table S15, Fig. S24).The synergistic power of prochloraz might be weaker for this species.Also, the interaction appears to be non-linear with respect to the concentration.This dataset showed one of the limits of the interaction module for identifying and simulating interactions.However, it is important to note that, for this species, the exposure was through feeding.Effects on feeding behaviour or repellency (e.g., food avoidance, or reduced assimilation rates) cannot be excluded and may have interfered with the interpretation of the toxic effects of the mixture.
Our modelling approach provided interesting insights for ecological risk assessment, such as the parameter values of the interaction module given by the model calibration on mixture data, which provide useful information (Table 2).The positive or negative value of b I indicates the orientation of the interaction (synergy or antagonism), and the threshold for interaction provides valuable information for the risk assessment.For example, it is interesting to note that for the mixture azoxystrobin and prochloraz, the threshold for interaction of prochloraz is nearly zero, meaning that any co-occurrence of those two fungicides will lead to synergism in E. crypticus.On the other hand, for the mixture cypermethrin and prochloraz in M. brassicae, the threshold for interaction is 14.96 (mg/kg) suggesting that a minimum concentration of prochloraz is required to trigger the synergistic interaction.
Finally, it is important to stress that the detection of interaction with our modelling approach is binary, interaction or no interaction, and we cannot tell how strong the interaction is, like in the classic approach with doseresponse models (EFSA, 2013).However, we consider the entire time course of the mixture effects.For instance, the interaction (e.g., a synergistic effect) can be strong during the first days of exposure and disappear afterward if the concentration of the chemical responsible for the interaction drops below the threshold for interaction.Therefore, the strength of the interaction will depend on 'when' we look at the effects.Model predictions for effects in untested exposure scenarios can be made for mixture effects with or without the interaction module at the desired time point

Fig. 1 .
Fig. 1.High-level conceptualisation of the processes involved in toxicity that can be affected by an interacting chemical translated into a reduced GUTS model (General Unified Threshold Model of Survival).A change in any of these processes will lead to a change in the scaled damage.

Fig. 3 .
Fig. 3. Observed and simulated (or predicted) survival over time of Bombus terrestris exposed to clothianidin and propiconazole in mixture.The lines in the top row are the predictions with the GUTS-RED-IT independent action model, and the confidence intervals (green shadows) represent the uncertainty of the model propagated from the calibration performed on the single exposures to the simulation of the mixture effects.The simulations in the bottom row are from the same model coupled with the interaction module fitted to the mixture data, and the confidence intervals (green shadows) represent the parametric uncertainty of the model with interaction module, calibrated on the mixture data.The dashed lines show the background mortality, fitted to the control treatment.Data from Robinson et al. (2017).

Fig. 4 .Fig. 5 .
Fig. 4. Observed and simulated (or predicted) survival over time of Enchytraeus crypticus exposed to azoxystrobin and prochloraz in mixture.The simulations in the left panel are the predictions with the GUTS-RED-SD damage addition model, and the confidence intervals (green shadows) represent the uncertainty of the model propagated from the calibration performed on the single exposures to the simulation of the mixture effects.The simulations in the right panel are from the same model coupled with the interaction module fitted to the mixture data (with prochloraz as the synergist), and the confidence intervals (green shadows) represent the parametric uncertainty of the model with interaction module, calibrated on the mixture data.The dashed lines show the background mortality, fitted to the control treatment.

Fig. 7 .
Fig. 7. Observed and simulated (or predicted) survival of Daphnia magna over time when exposed to cypermethrin alone (A), and with the addition of 1.4 μM of propiconazole (B and C).A) is the calibration of the GUTS-RED-SD model (comparable to Fig. 3A in Cedergreen et al., 2017), B) is the prediction based on A) coupled with the interaction module fitted (comparable to Fig. 3B in Cedergreen et al., 2017).And C) is the prediction based on A) without interaction.This figure is presented in the same manner as Fig. 3 in Cedergreen et al., 2017 for easier comparison.