Olivier Mélanie D'almeida

Salim Sire

Hélène Lardjane

Duval

Mélanie D'almeida

Olivier Sire

Salim Lardjane

Hélène Duval
email: helene.duval@univ-ubs.fr.

Development of a new approach using mathematical modeling to predict cocktail effects of micropollutants of diverse origins

Keywords: micropollutant, surface water, cocktail effect, orthogonal array composite design, response surface model micropollutant, EC: Effective Concentration, EDC: endocrine disrupting compound, HUVEC: Human Umbilical Vein Endothelial Cells, OACD: Orthogonal Array Composite Design, REACH: Registration, Evaluation, Authorisation and restriction

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Complex contamination of water by organic chemicals threatens ecosystems and human health as micropollutants (MP) exhibit potential harmful effects on non-target organisms. The numerous chemical components are generally found in all matrices and also accumulate in wild-life organisms. [START_REF] Pico | Contaminants of emerging concern in freshwater fish from four Spanish Rivers[END_REF] recently showed that contaminants of emerging concern are present as mixtures in most freshwater fish species at low levels in European rivers. MP are predominantly aromatic compounds used as pesticides, pharmaceuticals or industrial chemicals (41%, 30% and 16% of total contaminating organic compounds, respectively; [START_REF] Busch | Micropollutants in European rivers: A mode of action survey to support the development of effect-based tools for water monitoring[END_REF]. Pesticides widely used for agriculture, are constantly present in the environment and consequently, included in environmental monitoring programs. Poorly reduced by municipal wastewater treatment processes [START_REF] Stamatis | Occurrence and removal of fungicides in municipal sewage treatment plant[END_REF], significant levels enter river streams. Co-occurrence of multiple pesticides in a given water sample is often observed (2 or more pesticides in 90% urban and agricultural streams). Cyproconazole, one of the most frequent detected fungicide in influent and effluent samples with 93% and 79% detection percentages and also present in river samples [START_REF] Stamatis | Occurrence and removal of fungicides in municipal sewage treatment plant[END_REF], was chosen in our study, associated with diuron and terbutryn.

Pharmaceuticals are a class of emerging contaminants found in the aquatic environment with poor degradation and growing utilization rates. Like pesticides, pharmaceuticals are not completely removed from wastewater treatment plants, leading to their potential presence in drinking water [START_REF] Jones | Human pharmaceuticals in the aquatic environment a review[END_REF]. Among them, diclofenac, a prevalent non-steroidal anti-inflammatory drug (veterinary and human use) is often detected. Anti-inflammatories and analgesic drugs such as tramadol belong to the therapeutic group the most commonly detected in soils and sediments [START_REF] Biel-Maeso | Determining the distribution of pharmaceutically active compounds (PhACs) in soils and sediments by pressurized hot water extraction (PHWE)[END_REF]. A recent study [START_REF] Schmitz | Environmentally-relevant mixture of pharmaceutical drugs stimulates sex-steroid hormone production and modulates the expression of candidate genes in the ovary of juvenile female rainbow trout[END_REF] evidenced that combined occurrence of common pharmaceutical drugs at concentrations present in surface water environments may act as endocrine-disrupting compounds (EDC).

EDC are defined as exogenous substances (individual or mixture) belonging to a broad family of compounds that alter hormone secretion, potentially affecting different functions of the organism (individual, offspring or subpopulations): metabolism, reproductive functions, nervous system...The sources of exposure are numerous and difficult to control. Among them, bisphenol A, classified as category 1B reprotoxic substance from REACH legislation (Registration, Evaluation, Authorisation and restriction of CHemicals), was included in our study. Bisphenol A is ubiquitous in the environment as a result of its high production, consumption, and subsequent environmental introduction [START_REF] Corrales | Global Assessment of Bisphenol A in the Environment: Review and Analysis of Its Occurrence and Bioaccumulation[END_REF]. It is well documented that bisphenol A affects aquatic organism reproduction [START_REF] Luo | Differential effects of bisphenol A toxicity on oyster (Crassostrea angulata) gonads as revealed by label-free quantitative proteomics[END_REF][START_REF] Wang | Toxic effects of bisphenol A on goldfish gonad development and the possible pathway of BPA disturbance in female and male fish reproduction[END_REF].

Current risk assessment studies mostly focus on exposure to single chemicals. Studies involving mixtures of MP mainly deal with chemicals of the same group and are used as binary mixtures.

Pesticides and antifouling chemicals are used in more than 90 % and 75 % of the publications, respectively (reviewed by [START_REF] Cedergreen | Quantifying synergy: a systematic review of mixture toxicity studies within environmental toxicology[END_REF], without considering the large variety of chemicals present in a same location at a same time. Different mathematical models have been developed to predict the effects of mixtures. Hence, using Concentration Addition (CA) or Independent Action (IA) models, the overall toxicity of a mixture of known composition can be determined to some extent [START_REF] Kienzler | Regulatory assessment of chemical mixtures: Requirements, current approaches and future perspectives[END_REF]. Synergistic or antagonistic effects are also evidenced when compounds interact together [START_REF] Laetz | The synergistic toxicity of pesticide mixtures: implications for risk assessment and the conservation of endangered Pacific salmon[END_REF]. The effect resulting from the interaction between chemicals is called cocktail effect. Other studies use environmental samples where the mixture is tested as a whole. A screening approach consisting of the analysis of a maximum of chemicals is then compared with the individual chemical known effects. The monitoring of water quality is mostly based on chemical analyses of surface water using liquid chromatography coupled to tandem mass spectrometry. This approach gives accurate information regarding the structure and the number of numerous contaminants but is complex and cost-effective.

In this context, the objective of our study is to propose a new methodology using mathematical modeling to easily determine cocktail effects within complex mixtures containing diverse categories of MP, whatever their mode of action. We use a statistical design termed Orthogonal Array Composite Design (OACD) based on a two-level factorial design (two concentrations of MP) and a three-level orthogonal array (three concentrations of MP). This new class of composite design recently developed by [START_REF] Ding | Parliament and of the Council of 23 October 2000 establishing a framework for Community action in the field of water policy, Off[END_REF] and [START_REF] Xu | Combining two-level and three-level orthogonal arrays for factor screening and response surface exploration[END_REF] is effective for response surface modeling. A second-degree polynomial model is then applied to estimate the potential interactions between MP. Briefly, the experimental procedures were carried out with Human Umbilical Vein Endothelial Cells (HUVEC).

These cells are directly exposed to chemicals circulating in the blood and may consequently be a site for adverse effects of xenobiotics, metabolically activated by cytochrome P450 enzymes and/or other Ah receptor-regulated enzymes [START_REF] Annas | Differential response of cultured human umbilical vein and artery endothelial cells to Ah receptor agonist treatment: CYP-dependent activation of food and environmental mutagens[END_REF][START_REF] Dauchy | Expression and transcriptional regulation of ABC transporters and cytochromes P450 in hCMEC/D3 human cerebral microvascular endothelial cells[END_REF]. Six MP were selected: bisphenol A used in the chemical industry, diclofenac and tramadol from the pharmaceutical industry and cyproconazole, diuron and terbutryn, three pesticides from agriculture activities. The cytotoxicity of the six MP was first evaluated individually. Then different mixtures were tested at the concentration at which 10 % (Effective Concentration EC10) or 25% (EC25) of their maximal effect on human cytotoxicity was observed. The OACD allowed to reduce the number of the subsequent in vitro experiments using MP mixtures. A response surface model was then applied to look for cocktail effects.

Materials and Methods

Chemicals

Six chemicals corresponding to micropollutants of different categories and origins were used for the in vitro tests. 75 mM stock solutions were aliquoted in 100% methanol and stored at -20°C.

100X solutions were prepared in 100% methanol for the treatments. Bisphenol A (#15824) was purchased from Acros Organics. Other chemicals: cyproconazole (#46068), diclofenac (#93484), diuron (#45463), terbutryn (#45677) and tramadol (#42965) were purchased from Sigma.

Cell culture and single compound treatment

HUVEC obtained from Lonza, were maintained in Endothelial Growth medium (EGM-2, Lonza) at 37°C and 5% CO2. The medium was changed every 2 days, and the cells were passaged when they reached 80% confluence. They were used at passage 6 for all the treatments. HUVEC cells were seeded into 96-well plates (10 4 cells/well) and were allowed to adhere overnight in phenol red-free M200 medium (Gibco), supplemented with LSGS (Low Serum Growth Supplement kit), 10% (v/v) charcoalstripped-fetal bovine serum and 1% (v/v) penicillin-streptomycin before starting the experiments. The primary cell line was first treated with increasing concentrations of individual MP (0, 100, 200, 300, 400 and 500 µM) to determine their cytotoxicity and in particular the effective concentration for which 10% (EC10) and 25% (EC25) between the minimum and the maximum cytotoxic effect was obtained.

The untreated condition (0 µM) consisting of M200 contained 1% methanol (reference) as this concentration of methanol was used when performing all the treatments. Each treatment was performed in triplicate and three independent experiments were used for EC determination. Mixed pollutants were then tested as described in the OACD design section.

Crystal violet staining

HUVEC cells seeded into 96-well plates were treated with individual or mixed MP for 24 h. The medium with 1% methanol was used as a blank. At the end of the treatment, the medium in each well was discarded, wells were rinsed with Phosphate Buffer Saline, then 100 µl of formaldehyde 3.7% was added to fix the cells (10 min, room temperature). The plate was washed three times with deionized water before adding 0.1 % crystal violet (30 min incubation.) After 3 additional rinses, acetic acid 10% was added for 30 min prior to absorbance reading at 590 nm using FLUOstar® Omega Microplate Reader (BMG LABTECH, Ortenberg, Germany). % cell survival evaluated by the percentage of adherent cells was calculated as follows:

% = - - × 100 (1)
where Abs is the measured absorbance in each well; blank, the mean absorbance in blank (no cell in the well); reference, the mean absorbance for the reference (untreated cells). % mortality was then deduced (100 -% survival cells) to evaluate the cytotoxic effect.

Single micropollutant EC determination

Dose response curves were obtained for the six micropollutants: bisphenol A (BIS), cyproconazole (CYP), diclofenac (DIC), diuron (DIU), terbutryn (TER) and tramadol (TRA). Each dose response curve was used to determine the cytotoxic effective concentration EC50 using the following formula: = ( ) + ( 2)

with

( ) = + - 1 + ( !50 # ) $ (3)
where y is 100 times the percentage of surviving cells, x the concentration of micropollutant, a the minimal survival rate, b the maximal survival rate, EC50 the micropollutant concentration for which y=(a + b)/2. ε is the Gaussian observation error. An adjustment was done using the maximum likelihood method i.e. by minimizing the least square criterion:

! ( , , !50, ℎ) = (( ) -( ) )) * (4) , )-.
with a > 0, b > a, EC50 > 0, h > 0 and where n is the number of observations collected. yi is the observed percentage of surviving cells at concentration xi.

In order to obtain the EC10 (10% cytotoxicity) and EC25 (25% cytotoxicity) values, the adjusted function was reversed, that is to say ECy was calculated using the following formula:

! / = !50 × ( -1 + 1 -- ) . $ ⁄ (5)
The R software was used for all calculations.

OACD design

Since there were six micropollutants, it was not feasible to test in vitro all the combinations using three conditions corresponding to 0, EC10 and EC25 values (respectively coded by -1, 0 and +1).

Mathematical approach for cocktail effect analysis

The biological data obtained with MP mixtures were collected to perform the mathematical analysis, which was based on the adjustment of a second-order Response Surface Model (RSM) to the experimental data. The RSM is a second-order polynomial approximation for the unknown true functional relationship between the percentage of surviving cells and the concentration of micropollutants. This methodology was previously described in [START_REF] Cheng | Factor screening and response surface exploration[END_REF], [START_REF] Ding | Parliament and of the Council of 23 October 2000 establishing a framework for Community action in the field of water policy, Off[END_REF], [START_REF] Xu | Combining two-level and three-level orthogonal arrays for factor screening and response surface exploration[END_REF] and [START_REF] Weiss | A streamlined search technology for identification of synergistic drug combinations[END_REF]. The second-order RSM of degree two is given by the following equation:

= 1 2 + ( 1 ) ) + 3 )-.
( 1 ))

3 )-. )

* + ( ( 1 )4 ) 3 4-)5.

6

)-.

+ 7 (6)

where y is 100 times the percentage of surviving cells, xi (i = 1 to 6) represents the concentration of the micropollutants, ε is the observation error term. The adjustment was made using the least square method and the β's are the parameters of the model.

Physicochemical characteristics of micropollutants

SMILES (Simplified Molecular-Input Line-Entry System) strings using short ASCII (American Standard Code for Information Interchange) format were imported from PubChem database [START_REF] Weininger | SMILES, a chemical language and information system. 1. Introduction to methodology and encoding rules[END_REF]. The six SMILES strings were compared using Hierarchical Clustering Analysis (HCA) with Ward's method in the R core environment. n-octanol/water partition coefficient or Kow is an important parameter to estimate environmental exposure risk as it is assumed that the uptake of an organic substance depends on its hydrophobicity. Attempts were made to correlate Log(Kow) values with the observed EC10 and EC25 values (linear regression).

Results

HUVEC cells were first treated with increasing concentrations of each micropollutant for 24h to establish dose response curves. Figure 1 shows the dose response curves we obtained for the six individual MP. Three independent experiments testing each dose in triplicate were performed to consolidate the concentration response pattern. The Hill equation was selected for the best-fitting regression model. The survival data collected using crystal violet staining were used to calculate cytotoxic effects. 100% mortality was obtained at higher doses for three MP: bisphenol A, cyproconazole and terbutryn. Tramadol used for human pain management was the less cytotoxic. The effective concentration (EC) values obtained for the different micropollutants are indicated in Table 1.

Distinct EC10 and EC25 values were observed among the different MP, which highlights the need of concentration-effect analysis instead of using fixed and identical concentrations. Additionally, the Q-Q plot analyzing the sample quantile versus theoretical quantile appears linear, which means that the residuals come from the expected normal distribution (Figure 3).

There is a slight deviation at each extremity, suggesting slightly more extreme values than expected.

This analysis of residues confirmed that the fit was correct. The proposed model contains regression parameters that include coefficients for main effects regarding one MP, coefficients for quadratic main effects regarding one MP and coefficients for two MP interaction effects. A first exploratory adjustment of the response surface, including a screening of the highest cell survival values, highlighted that relatively few quadratic terms had significant coefficients different from 0, at the level of 10 % statistical significance. The model was readjusted by keeping only the terms whose coefficient was significantly different from 0. We get an R 2 of 0.7037 and an adjusted R 2 of 0.6907. The coefficients obtained are shown in Table 2. The Fisher statistic value is 54.38 with 10 and 229 degrees of freedom, with a p-value lower than 2.2e-16, which confirms that the response surface model is globally significant.

Table 2: Estimates of the parameters for the response surface model. All coefficients are significantly different from 0 to 10% level. "." corresponds to a p-value less than 0.1 and "*" to a p-value less than 0.05. "**" and "***" correspond to lower levels.

Regarding the linear effects, all micropollutants except diclofenac are identified as significant.

Concerning the other parameters, the quadratic term CYP 2 and the bilinear effects CYP × DIC and BIS × TRA are the most significant terms (p-value < 0.05).

Next, the response surfaces corresponding to individual micropollutants were deduced from the adjusted response surface, replacing the level of all others by -1. Each response surface is noted as a function f1(x1), ..., f6(x6) where x1 is the concentration of BIS, ..., x6 the concentration of TRA. When assuming independent effects, the cell survival rate is given by the formula:

100 × 1 100 6 9 ) (: ; ) 3 )-.

(7) Thus, we can compare the survival rate obtained under the assumption of independence with the survival rate obtained from the global response surface, for all combinations of levels (-1, 0, +1) of the micropollutants. This is done graphically in Figure 4. Under the assumption of independence of the effects, these must be placed along the first bisector of the Cartesian plane. As shown in Figure 4, these are well above. Hence the hypothesis of independence seems unjustified from the graph. The values obtained using the response surface being mostly much higher than those obtained under the hypothesis of independent effects, it was deduced that there is a strong global cocktail effect between the micropollutants. This effect can be investigated by considering micropollutant response surfaces taken in pairs, the concentrations of the others being set at -1 in the overall response area. For example (Figure 5), for BIS and TER, there is no interaction as the contour diagram reflects it and as it is reflected in the response surface equation, where there is no significant interaction term corresponding to BIS and TER. It is different for BIS and TRA: the interaction term in the response surface having a significantly positive coefficient, it can be concluded that there is a neutralization between these two micropollutants. Eventually, Log(Kow) values of the six MP, used to measure their hydrophobicity, were compared to the observed EC10 and EC25 values (Figure 6A). The absence of correlation suggests the absence of a link between the uptake of each MP and its biological effect. We then used SMILES strings [START_REF] Weininger | SMILES, a chemical language and information system. 1. Introduction to methodology and encoding rules[END_REF] to classify the six MP according to their physicochemical characteristics (Figure 6B).

The dendrogram calculated with the Ward's Hierarchical Clustering Method shows two clusters: one containing bisphenol A and another one containing tramadol (neutralization with bisphenol A) and terbutryn (no interaction with bisphenol A). These results suggest that the cocktail effect can't be predicted using structural analysis of MP.

Discussion

We present a study focusing on MP cocktail effects where a selection of micropollutants was made to reflect the heterogeneity of water contamination. Hence, six organic compounds, namely cyproconazole, diuron, terbutryn, bisphenol A, diclofenac and tramadol were tested. Single MP doseresponse curves were first generated for experimental data collection. Three doses per MP were then retained for complex mixture analysis: EC25, the highest tested concentration, the dose where 25 % of the maximal response was observed, EC10, where 10 % of the maximal response was observed, and dose 0, where no MP was present. An OACD plan lead to test 55 combinations in 4 replicates, instead of the complete factorial design plan made of 3 6 combinations. Then a response surface model allowed to generate significant regression coefficients corresponding to single-chemical linear and quadratic effects as well as two-chemical pair-wise interaction effects. A neutralizing effect between micropollutants was observed, in particular, between bisphenol A and tramadol.

The proposed methodology was inspired by clinical methodologies developed to test the effect of drug combinations [START_REF] Weiss | A streamlined search technology for identification of synergistic drug combinations[END_REF]. The advantage of using drug mixtures is that they often have higher efficacy and lower drug dosages than individual drugs. From an ecotoxicological point of view, with the increase of MP present in the aquatic environment, important questions are emerging: is the predicted toxicity of single compounds realistic in presence of other micropollutants, with same or different modes of action?

The mixture risks are a major concern as illustrated with the recent growing number of studies dealing with mixtures of MP and mathematical modeling for combined effect prediction [START_REF] Altenburger | What contributes to the combined effect of a complex mixture?[END_REF][START_REF] Altenburger | Algal toxicity of nitrobenzenes: combined effect analysis as a pharmacological probe for similar modes of interaction[END_REF][START_REF] Ermler | Genotoxic mixtures and dissimilar action: concepts for prediction and assessment[END_REF]. [START_REF] Altenburger | Mixture effects in samples of multiple contaminants -An inter-laboratory study with manifold bioassays[END_REF] compared predicted with observed combined effects using twelve micropollutants at the estimated low effect levels (10 or 20% maximal effect). The use of 21 bioassays show variable results: response most of the time in agreement with concentration additive model, or statistically higher or lower activity than expected. The independent action model was the best way to predict the observed mixture effects when looking at apical effects (observable outcome in a whole organism). Usually, concentration addition model is retained for the similarly acting components, and subsequently, the dissimilarly acting compound groups are analyzed using independent action model [START_REF] Altenburger | What contributes to the combined effect of a complex mixture?[END_REF]. In environmental mixtures, the information on modes of action for all compounds is, however, not always easily available.

Very few studies investigate the robustness and validity of mathematical models to assess the additive effect of micropollutants with the same mode of action in complex and more realistic mixtures. Using human and fish reporter gene assays, [START_REF] Serra | Combined effects of environmental xeno-estrogens within multi-component mixtures: Comparison of in vitro humanand zebrafish-based estrogenicity bioassays[END_REF] tested the combined effects of environmental xenoestrogens, among them the binary mixture bisphenol A and genistein, within multicomponent mixtures, comprising diuron and diclofenac. They evidenced that non-estrogenic chemicals can interfere with the response to a mixture of xenoestrogens; for example, the triazole fungicide propiconazole decreased estradiol-induced estrogenic activity in a non-specific manner.

Using the response surface model, we show for the first time that bisphenol A interacts with the painkiller tramadol. It is largely described that bisphenol A perturbs normal endocrine functions. Surprisingly, it can also act on the nervous system [START_REF] Ogiue-Ikeda | Rapid modulation of synaptic plasticity by estrogens as well as endocrine disrupters in hippocampal neurons[END_REF], modulate synaptic plasticity in the hippocampus, interfere with anxiety-related behavior in rats [START_REF] Farabollini | Perinatal exposure to the estrogenic pollutant bisphenol A affects behavior in male and female rats[END_REF] or disrupt physiologic neurotransmitter function [START_REF] Alyea | Differential regulation of dopamine transporter function and location by low concentrations of environmental estrogens and 17beta-estradiol[END_REF]. More recently, [START_REF] Mizuta | Bisphenol A inhibits compound action potentials in the frog sciatic nerve in a manner independent of estrogen receptors[END_REF] showed that bisphenol A could inhibit compound action potentials (nerve conduction block) in the frog sciatic nerve, like tramadol [START_REF] Katsuki | Tramadol, but not its major metabolite (mono-O-demethyl tramadol) depresses compound action potentials in frog sciatic nerves[END_REF]. Our results show that these two chemicals can interact within the cell, whatever their mode of action. This cocktail effect can't be predicted using structural analysis of MP.

In conclusion, there is an urgent need to prioritize the risk assessment of chemical mixtures by testing MP under REACH legislation. We propose a method that can be transposed to in situ studies and monitoring programs. Hence, once the most proeminent micropollutants in a given site have been identified, our mathematical approach can be applied to predict situations at risk. Other bioassays such as Reactive Oxygen Species detection, nuclear receptor-based responses… can be used to evaluate

  Hence, the complete factorial design plan 3 k (three levels and 6 factors, making 3 6 = 729 tests without any replicate) was replaced by a reduced factorial design using FrF2 and DoE.base packages (R software). The proposed OACD was composed of a subset of the experimental conditions carefully chosen in order to maintain the desired level of statistical power. This fractional factorial design of resolution V combining a two-level fractional design and a three-level orthogonal array comprised a total of 55 conditions (supplementary data). These 55 tests, for k=6 factors, denoted by x1, …, xk, consisted of 34 cube points with all xi = -1 or +1 (2 levels), 6 center points with all xi = 0 and 15 additional points with xi = -1, 0 or +1 (3 levels). We performed four replicates, to which 45 central tests (center points) were added. Eventually, 265 tests were performed in vitro. Experiments were carried out in a randomized order (shuffling program).

Figure 1 :

 1 Figure 1: Dose response curves following single MP treatment. After 24h treatment of HUVECs with increasing concentrations of each micropollutant, a dose-response curve was obtained (BIS: bisphenol A; CYP: cyproconazole; DIC: diclofenac; DIU: diuron; TER: terbutryn; TRA: tramadol).

Figure 2 :

 2 Figure 2: Comparison of residuals and predicted adjusted values. The plot of residuals against predicted adjusted values allowed to assess the response surface model.

Figure 3 :

 3 Figure 3: Analysis of the residuals' distribution. Data distribution was checked by Q-Q Plot analysis.

Figure 4 :

 4 Figure 4: Comparison of cytotoxic effects estimated by the response surface model versus the independent action model.

Figure 5 :

 5 Figure 5: Contour plot of response surface for BIS and TER (left panel) and for BIS and TRA (right panel) when the concentration of the other MP is set at -1 in the overall response surface.

Figure 6 :

 6 Figure 6: Physicochemical characteristics of MP. (A) Log(Kow) values were compared with biological activity of MP. (B) SMILES strings were used for Hierarchical Clustering Analysis of MP with Ward's method (Numbers indicate Log(Kow) values).

environmentally relevant MP concentrations. Hopefully, new legislation dealing with priority cocktail monitoring will appear to reinforce micropollutant monitoring policy.

Conflicts of Interest:

The authors declare no conflict of interest