Nathanael Dougier

Pierre Garambois
email: pierre.garambois@ensam.eu

Julien Gomand

Lionel Roucoules

Efficient energy system modelling for multi-objective optimisation

Keywords: decision support, energy systems modelling, multi-objectives optimisation

Energy systems, on which our modern society rely, are in constant transformation. Technological evolution, climate change or the finitude of fossil fuels are some reasons to rethink the centralized, carbon-based energy networks. This way, the design of future energy systems have to take into account multiple concerns, such as local resi lience, in addition to technical and economic ones. This paper presents a decision-support tool for the conception of energy systems focusing on the electric vector. The tool was designed using an energy system model implemented in an optimisation algorithm. It takes into account several constraints simultaneouslyequilibrium between production and consumption as well as resources availabilityand assess the influence of technical parameters on the global performances of the system. An energy system is considered as a combination of production, storage and transport technologies with their operating strategies. The tool's modularity allows to choose the models adapted to a quick optimisation of energy systems or to an analysis of technical parameters. The second part of the paper presents the optimisation of a local energy system. Search space is composed of production and storage technologies' number and their operating strategies. Main goals are to find trade-offs between different economic and technical objective-functionssuch as levelized cost of energy or local autonomy. Therefore, a genetic algorithm method was used to perform a multi-objective optimisation based on the model. The impact of the operating strategy adopted is underlined.

INTRODUCTION

Energy consumption has increased for more than a century to sustain our changing way of life. Among energy vectors, electricity presents the biggest development of the past fifty years and may highlights advantages in terms of environmental concerns. Huge networks were built around main power plants, nuclear and hydroelectric ones in France. Economic challenges of the twentieth century have evolved and society's concerns are more related to social and environmental impacts. Same set of technologies cannot be designed to balance the increase in consumption and the end-of-life of current power plants. Therefore, new energy systems (ES)i.e. sets of production, storage and transport technologiesneed to be created. The development of distributed energy sources may represent an alternative to the current technologies. However, as they are dependent on intermittent resources, they can be supplemented by storage technologies. These configurations allow a local design for energy systems, microgrids, partially autonomous from the main grid if necessary. Economic, technical and environmental indicators have to be taken into account to assess the relevance of these new systems. This paper focuses on the conception of ES. The purpose of the work described below is to create a decisionsupport tool able to find, for a local context, various optimal energy systems that are compromises between economic and technical objectives. The tool uses a physical modelling that simulates the operation of an energy system. Then it runs this model into an optimisation algorithm to find optimal energy systems.

Energy system modelling can be performed following two main methods according to the purpose. Prospective studies that assess huge networks and determine the best choice to make between a range of technologies will use a simplified modelling. Some of them use bottom-up "technological" methods [START_REF] Assoumou | Modélisation MARKAL pour la planification énergétique long terme dans le contexte français[END_REF] with some physical background but others prefer top-down "economical" modelling with aggregated variables [START_REF] Nordhaus | Evolution of modeling of the economics of global warming: changes in the DICE model, 1992-2017[END_REF]. The time step of their simulation can range from hours to years [START_REF] Assoumou | Modélisation MARKAL pour la planification énergétique long terme dans le contexte français[END_REF]. On the other hand, studies analyse precisely the operation of one specific energy system [START_REF] Darras | PV output power fluctuations smoothing: The MYRTE platform experience[END_REF]. Their models are very detailed according to the technical parameters they want to assess and the time step used is small enough to capture the physical phenomena analysed, influence of temperature on the production of PV panels [START_REF] Xu | Optimized sizing of a standalone PV-wind-hydropower station with pumped-storage installation hybrid energy system[END_REF]. Whereas prospective studies consider often only a few variables like the technologies' number and use a long time step, detailed studies take also into account management parameters and short time steps. The study of local technological parameters requires detailed models. On the other hand, a model designed to minimise the computation-time (CPU-time) of the optimisation needs to be simplified. This conflict can be solved by finding compromises between accuracy and CPU-time. The present work tries to solve the conflict by using a systemic approach, i.e. several models with different levels of detail.

Once the energy system modelled and its operation simulated, an optimisation process can be performed to find the best solutions. Although systems need to be compared following many various indicators, optimums are often found using single-objective (SO) optimisation with an economic indicator [START_REF] Wasilewski | Optimisation of multicarrier microgrid layout using selected metaheuristics[END_REF]. Optimisation is sometimes performed with both economic and environmental objectives but using a prior weighting to be adapted to SO algorithm [START_REF] Kayal | Optimal mix of solar and wind distributed generations considering performance improvement of electrical distribution network[END_REF]. More rarely, multi-objectives (MO) optimisations without prior weighting are performed [START_REF] Kamjoo | Multi-objective design under uncertainties of hybrid renewable energy system using NSGA-II and chance constrained programming[END_REF]. To tackle such complex problems in the most exhaustive manner, it seems important to take into account several impacts of energy systems, especially as they can be antagonists. Moreover, finding compromises does not overly restrict the tool user's choice and allows the analysis of interesting systems that would not have been considered initially. The latter can apply further discrimininating criteria to selected the best solution for its situation.

In the case of energy systems modelling and sequential simulation, indicators cannot be formulated explicitely, hence classical Lagrangian optimisation methods cannot be used. The use of metaheuristic optimisation algorithms allowing studying non-explicit and non-weighted objective-functions is required. Numerous methods exist such as gradient-method, Particle Swarm Optimization or MO Evolutionary Algorithms (MOEA). The advantage of evolutionary, also called genetic, algorithms is their ability to deal with 0-order objective functions and to explore the entire design space. In the literature, several methods have been used like Strength Pareto Evolutionary Algorithms (SPEA [START_REF] Zitzler | Multiobjective evolutionary algorithms: a comparative case study and the strength Pareto approach[END_REF] and SPEA-II [START_REF] Zitzler | SPEA2 : Improving the Strength Pareto Evolutionary Algorithm[END_REF]), Non-dominated Sorting Genetic Algorithms (NSGA and NSGA-II [START_REF] Deb | A Fast Elitist Non-dominated Sorting Genetic Algorithm for Multiobjective Optimization: NSGA-II[END_REF]), Pareto Archived Evolution Strategy (PAES [START_REF] Knowles | The Pareto Archived Evolution Strategy: A New Baseline Algorithm for Pareto Multiobjective Optimisation[END_REF]) or Adaptive Pareto Algorithm (APA [START_REF] Dumitrescu | A New Evolutionary Adaptive Representation Paradigm[END_REF]). They allow finding Pareto-optimal solutions between the different objectives, i.e. ES that dominate the others on at least one objective, but they are time-consuming due to the repetition of the criteria's evaluation for each individual. The ES modeling choices and performance-to-CPU-time ratio of the simulation is all the more important.

The paper is structured as follow. Section 1 describes the modelling and the simulation of energy systems. Optimisation method, objectives and parameters are detailed in Section 2. Section 3 presents and discusses the results of the casy study. Eventually a conclusion underlines the main results and presents some perspectives.

MODELLING AND SIMULATION

The challenge of energy system modelling is to represent accurately the behaviour of the system with a computation-time (CPU-time) adapted to the use of a genetic optimisation algorithm.

Models

In order to find efficient and innovative energy mixes, this work aims at analyzing the influence of local technological parameters on the global performances of the system. Therefore, detailed models need to be developped. However, to remain compatible with the optimisation needs, several layers of modelisation have been implemented, from global pattern to very detailed models, all represented by interconnected blocs. Simulation and optimisation may then be able to use a layer or another according to the precision-to-CPUtime ratio. The advantage of this approach is its modularity and its ability to analyse both the global behaviour and the physical and technological frontiers of energy systems. So far, only low level models have been implemented, allowing access to the main parameters of each technology, surface and efficiency of photovoltaic (PV) panel or power curve of the wind turbine for example. Uncontrollable renewable energy sources are modeled thanks to basic proven models [START_REF] Ogunjuyigbe | Optimal allocation and sizing of PV/Wind/Split-diesel/Battery hybrid energy system for minimizing life cycle cost, carbon emission and dump energy of remote residential building[END_REF]. For storage technologies, a common model is used with different parameters according to the technology [START_REF] Multon | Systèmes de stockage d'énergie électrique[END_REF]. All technologies are defined by their capacity, their charging and discharging rates, their loss coefficient during charge and discharge and their coefficient and reference period for auto-discharge. The evolution of their energy level is measured with a common criterion, the State of Energy defined as the ratio between the contained energy in the storage at a time over the installed capacity [START_REF] Multon | Systèmes de stockage d'énergie électrique[END_REF]. Controllable technologies, such as gas power plants, are characterised by their installed power and the ramp to increase the output electrical power.

Management parameters

Models briefly presented above concern the design of technologies, not the way they work together. The set of rules describing the operation of an energy sub-system is called here a management strategy. Studies on microgrids focusing on conception choices usually consider fixed strategies [START_REF] Wang | Multicriteria Design of Hybrid Power Generation Systems Based on a Modified Particle Swarm Optimization Algorithm[END_REF]. When consumption needs to be balanced, the developed behaviour model considers that technologies have to deliver the correct amount of electricity in a sequential manner, one after the other. It relies on the production priority order defined between technologies, usually considered fixed in prospective models [START_REF] Assoumou | Modélisation MARKAL pour la planification énergétique long terme dans le contexte français[END_REF].

Many performances indicators used in this study are based on the actual production of each technology, like the greenhouse gases emissions or economic indicators. The energy production depends on the management parameters described above, strategy and production priority order. Thus, it seems important to take into account these parameters among other technological ones. Depending on the technology's categories (controllable / uncontrollable / storage technologies), strategies may be applicable or not (see Table 1).

Table 1: Description of the management strategies and their potential application to production controllable (CTRL), uncontrollable (UNCTRL) and storage technologies (STOR)

Strategy number

Strategy description CTRL UNCTRL STOR

1 Only balance consumption in the limit of 𝑃 𝑛𝑜𝑚 * X X 2
Balance consumption and store in the limit of 𝑃 𝑛𝑜𝑚 X X 3

Balance consumption in the limit of the 𝑃 𝑛𝑜𝑚 and store in the limit of 𝑃 𝑜𝑝𝑡𝑖 * X 4

Only balance consumption in the limit of 𝑃 𝑜𝑝𝑡𝑖 X X 5

Balance consumption and store in the limit of 𝑃 𝑜𝑝𝑡𝑖 X 6

Produce at 𝑃 𝑜𝑝𝑡𝑖 in the limit of the consumption and storage capacities X 7

Store only in the limit of 𝑃 𝑜𝑝𝑡𝑖 X *𝑃 𝑛𝑜𝑚 is the nominal power of the technology and 𝑃 𝑜𝑝𝑡𝑖 its optimal power, for which efficiency is the best

Simulation

A sequential simulation is performed using the two management parameters: strategy for each technology and a general priority order. The time step may be chosen in a range from minutes up to hours, depending on the phenomena we want to take into account (management, wind speed, river water flow), the accuracy needed and the computation-time requirements. According to the management strategy and the power balance at each time step, energy system model allocates sequentially the power of each power plant to balance consumption and even store electricity in some cases. Simulation's output are the powers produced, stored or lost by each technology for each time step, along with the state of energy of the storage devices and various performance indicators detailed below. An example of output power and state of energy is shown on Figure 1.

Multi-objectives optimisation method

Among the possible metaheuristic optimisation methods able to consider non-explicit and non-weighted objective-functions evolutionary algorithms have been chosen for their ability to explore the entire design space. Among the different evolutionary algorithms, choice has been made to implement the developed model in the NSGA-II algorithm which is efficient, freely available and already mastered in the laboratory [START_REF] Deb | A Fast Elitist Non-dominated Sorting Genetic Algorithm for Multiobjective Optimization: NSGA-II[END_REF].

Genetic algorithms work like the evolution of species: an initial population with individuals characterized by their genes evolves over generations and the most adapted ones are more likely to survive. Here, the initial population is composed of energy systems characterized by the decision variables of the optimisation (number of each technology). Adaptation to the environment is measured through objective-functions and a generation is an iteration of the systems' selection and evolution processes. The final population presents different combinations of genes than the original one due to the need to adapt to its environment. The process followed by such an algorithm is detailed in Figure 2.

The evolution of genescombinations of technologies with their management strategiesover generation happens thanks to two processes, crossover and mutations. The first method consists in the random combination of two individuals' genes to create a new individual, i.e. the number of wind turbines of a system applied to another. Mutation is the random mutation of an individual's genes. It allows exploring the entire conception space and avoiding local optimums. NSGA-II algorithm is usually used in the literature performing one of the two processes with a crossover probability of 90% and thus a mutation probability of 10%. Same configuration was used in this study. The best energy systems, i.e. the best compromises between the different optimisation's objectives, between the initial generation and the new systems created are selected to form a new generation the same size as the initial one. The last generation is composed of Pareto-optimal solutions.

Objective-functions

The purpose of the multi-objectives optimisation is to find energy systems that are trade-offs between various objectives. According to the situation comparison criteria can change but the three main categories of criteria taken into account in this work are economic and technical ones. In order to underline the interest of the multiobjectives approach without prior weighting, one objective of each category has been chosen for this study.

Economic objective

The main objective used to compare energy systems is the cost of the system. It can be the investment cost, the operation and maintenance (O&M) cost, the lifecycle cost or the payback period for example. An indicator often used in the energy industry is the levelized cost of energy (LCOE, in €/kWh produced) which represents the total amount of money spent over the lifetime of a power plant divided by the energy produced during the same period, as described by Equation 1 [START_REF] Ouyang | Levelized cost of electricity (LCOE) of renewable energies and required subsidies in China[END_REF]. Its popularity can be explained by the ease to understand it and the fact that it covers the overall lifetime costs and takes into account currency discount over time. LCOE can only compare technologies for a same situation (location and time). Unfortunately, its does not take into account the provided services like the frequency regulation or the reliability of the production. In a multiobjectives approach, other technical indicators can balance this lack.

Another major issue with LCOE is that elements considered to compute the lifecycle costland lease, insurance costs, taxes or carbon emissions for examplevary according to the source. Therefore, it is important to ensure a consistency between data for all technologies. In this work, LCOE formulas were taken from the Fifth report of the IPCC [17, p. 1333]. Nevertheless, storage technologies are not evaluated in this report and their LCOE values were taken from another study [START_REF] Schmidt | Projecting the Future Levelized Cost of Electricity Storage Technologies[END_REF]. In a future work, a unique methodology to compute LCOE could increase the accuracy. With 𝑛 the year, 𝑁 the expected lifetime, 𝐼 𝑛 the investment costs during year 𝑛, 𝑂 𝑛 the operation costs during year 𝑛, 𝑀 𝑛 the maintenance costs during year 𝑛, 𝐷 𝑛 the residual value and 𝑑 the discount rate.

In order to discriminate more efficiently energy systems, the economic objective used is the ratio of the lifecycle cost (LCCost, in €) of the system correspondig to the simulated period, deduced from Equation 2.

𝐿𝐶𝐶𝑜𝑠𝑡 𝑠𝑦𝑠𝑡𝑒𝑚 = ∑ 𝐿𝐶𝑂𝐸 𝑡𝑒𝑐ℎ𝑛𝑜 * 𝑃𝑟𝑜𝑑𝑢𝑐𝑡𝑖𝑜𝑛 𝑡𝑒𝑐ℎ𝑛𝑜 𝑡𝑒𝑐ℎ𝑛𝑜 (2)

Technical objective

The main purpose of an energy system is to provide enough electricity to the consumers at the right time.

Technical indicators measure the reliability of power supply. Among existing indicators in the literature, we can encounter the energy not supplied [START_REF] Sanajaoba | Maiden application of Cuckoo Search algorithm for optimal sizing of a remote hybrid renewable energy System[END_REF], loss of load expectation (LOLE) [START_REF] Sanajaoba | Maiden application of Cuckoo Search algorithm for optimal sizing of a remote hybrid renewable energy System[END_REF], loss of load probability (LOLP) [START_REF] Khatib | Optimal Sizing of Hybrid PV/Wind Systems for Malaysia Using Loss of Load Probability[END_REF], wasted renewable energy [START_REF] Perera | A hybrid tool to combine optimization and multi-criterion decision making in designing standalone hybrid energy systems[END_REF] or autonomy level [START_REF] Celik | Techno-economic analysis of autonomous PV-wind hybrid energy systems using different sizing methods[END_REF]. This work focuses on local energy systems, so the designed system can still rely on the main grid if needed and does not require a full autonomy. Therefore, the autonomy level (LA) for the main grid seems to be adapted to measure the reliability of local supply.

Two kind of autonomies can be analysed. Autonomy in terms of time measures the percentage of time when the energy demand is balanced by enough production. Energy autonomy, as defined here, measures the total energy produced to balance consumption and losses but also to store additional energy in the storage devices. Energy indicator ranges from zero, when there is no production, to the sum of the energy demand plus the additional energy that can be stored over the energy demand. When the energy autonomy is above one, it that the total energy added to the system during the operation period represent more than the consumption and could have been used some other time before.

𝐿𝐴 𝑒𝑛𝑒𝑟𝑔𝑦 = 𝐸 𝑝𝑟𝑜𝑑𝑢𝑐𝑒𝑑 𝐸 𝑑𝑒𝑚𝑎𝑛𝑑 +𝐸 𝑙𝑜𝑠𝑠𝑒𝑠 (3)
With 𝐸 𝑝𝑟𝑜𝑑𝑢𝑐𝑒𝑑 the total energy produced during the period, 𝐸 𝑑𝑒𝑚𝑎𝑛𝑑 the cumulated energy consumed during the period and 𝐸 𝑙𝑜𝑠𝑠𝑒𝑠 the total losses during the period (transport and storage losses).

CASE STUDY

The following case study aims at underlining the ability of the designed tool to find optimal energy systems and the interest of the multi-objectives optimisation without prior weighting to explore the entire conception space in the case of energy systems.

Parameters

A local electric consumption curve for 6000 typical homes over 24h along with meteorological data (wind speed and solar irradiation) are considered as inputs with a 90 minutes time-step. Decision variables of the optimisation are detailed in Table 2. Maximum value and tolerance of each variable has been defined according to the maximum power consumed during the period, which is 3.2 MW. Initial energy systems population's size has been set to 1000 inidividuals in order to represent sufficient variety and ensure that the algorithm will converge toward global optimums and not local ones. Following convergence considerations, a maximum of 300 generations was taken into account.

As described before, the two objectives considered are LCOE and energy autonomy from the main grid in order to represent the different, antagonist, categories of impact of energy systems.

Table 2: List of the decision variables with the range of their values and the tolerance associated

Decision variables Range Tolerance

Number of wind turbines (Vestas V90 -2MW) 0 -3 1 turbine Surface of PV panels 0 -64,000 m² 3000 m² Installed power of biomass power plant 0 -3.

Results

First, an optimisation with two objectives is performed. Figure 3 shows the values of the energy autonomy and the lifecycle costs of several energy sysems. The blue dots represent the Pareto front of the multi-objective optimisation. These results are compared to solutions of multi-objectives optimisations with prior weighting.

Figure 3: Pareto front of the optimisation's solutions

Single-objective optimisations present as expected a very good performance following their optimised objective but a bad one following the others, as underlined by the result of the energy autonomy optimisation. Lifecycle cost optimisation' autonomy is not too bad due to the added constraint. Weighted objectives show relatively good performances compared to the multi-objective Pareto front, with the solution of the equally weighted objectives on the front. However, it gives only one set of performances and parameters and the weights have to be chosen carefully a priori otherwise the performances of the solution varies a lot. The issue having only one solution is the risk to miss another energy system with potentially close performances but different parameters.

Some interesting energy systems were selected whose variables' values are represented on Figure 4. In the first graph, the installed power of each technology is represented, on the second the storage capacities and on the third the chosen priority order and control strategy. Two solutions of the multi-objectives optimisation, here "Multi-obj 1" and "Multi-obj 2", have been represented because they have the same performances has the weighted solution with equal priority between objectives, here 'Weighted 1'. Although they have the same power of gas installed, the capacity of the storage technologies differs along with the priority order and control strategy. Whereas the weighted solution's priority order puts first the storage devices, then the renewable energies and finally the controllable technologies, the other solutions place the controllable producers first in the priority order. Analysing technical parameters like control strategies and priority order inside a multiobjectives optimisation allows to distinguish energy systems with very close performances and thus to explore the conception space, what is impossible with a single-objective optimisation. The difference between the three weighted solutions underlines the influence of the weigthing choice. The solution giving more importance to the autonomy objective shows a greater variety of sources whereas the one with a bigger weight for the cost minimize the power installed, leading to a very low autonomy (selection has been made to avoid 0% autonomy). Nature of the optimal energy system depends on the weights chosen a priori. The multi-objectives optimisation applied here allows finding more diverse energy systems.

CONCLUSION

Power plants aging and new society's concerns are pushing for a rethink of the energy systems design. Distributed technologies with additional storage to balance their intermittency allow building local strategies in order to increase the autonomy from the main grid. However, these new options have to be compared to the existing network with economic, technical and environmental indicators. The work presented above aims at creating a decision-support tool to model, simulate and optimise energy systems. It focuses on technological parameters to model accurately power exchanges and to suggest innovative systems. In particular, management strategy is taken into account as a decision variable. To represent the complexity of energy systems design and the antagonism between some objectives, a multi-objectives optimisation algorithm without prior weighting, NSGA-II, has been performed to find different energy systems providing partial autonomy.

Non weighted multi-objectives optimisation allows finding more interesting energy systems than with a prior weighting. The exploration of the entire conception space leads to more diverse solutions. By exploring the design space without prior weighting, the developed tool offers a greater choice to the user and reveals energy systems compromises which the user would not have thought of. Especially, energy systems with very close performances but different characteristics have been found. Discriminating criteria can be applied afterward.

Robustness of the proposed mix to the variation of meteorological conditions or production capacities will be investigated in a future work to ensure a bigger reliability and more detailed models will be implemented in the tool in order to analyse the influence of local parameters.

Figure 1 :

 1 Figure 1: One day sequential simulation of an energy system -(a) Cumulated power (kW) produced/stored along with the power demand of 6000 typical homes ; (b) State of energy (kWh) of the storage devices for each time step

Figure 2 :

 2 Figure 2: Process to create a new generation with the genetic algorithm NSGA-II

  𝐿𝐶𝑂𝐸 = 𝑠𝑢𝑚 𝑜𝑓 𝑑𝑖𝑠𝑐𝑜𝑢𝑛𝑡𝑒𝑑 𝑐𝑜𝑠𝑡𝑠 𝑜𝑣𝑒𝑟 𝑙𝑖𝑓𝑒𝑡𝑖𝑚𝑒 𝑠𝑢𝑚 𝑜𝑓 𝑒𝑙𝑒𝑐𝑡𝑟𝑖𝑐𝑎𝑙 𝑒𝑛𝑒𝑟𝑔𝑦 𝑝𝑟𝑜𝑑𝑢𝑐𝑒𝑑 𝑜𝑣𝑒𝑟 𝑙𝑖𝑓𝑒𝑡𝑖𝑚𝑒 = ∑ (𝐼 𝑛 +𝑂 𝑛 +𝑀 𝑛 +𝐷 𝑛 ) (

Figure 4 :

 4 Figure 4: Content of selected energy systems solutions of various optimisations

ACKNOWLEDGMENT

The research work reported here was made possible by the Region Provence-Alpes-Côte d'Azur.


Corentin Mit

Cleo Tebby

Tristan Gueganno

Anne Bado-Nilles

Rémy Beaudouin
email: remy.beaudouin@ineris.fr

Modeling acetylcholine esterase inhibition resulting from exposure to a mixture of atrazine and chlorpyrifos using a physiologicallybased kinetic model in fish

Keywords: PBTK-TD model, Atrazine, Chlorpyrifos, Mixture, Acetylcholinesterase, Carboxylesterase

Aquatic organisms are exposed to mixtures of chemicals that may interact.

Mixtures of atrazine (ATR) and chlorpyrifos (CPF) may elicit synergic effects on the permanent inhibition of acetylcholinesterase (AChE) in certain aquatic organisms, causing severe damage. Mechanistic mathematical models of toxicokinetics and toxicodynamics (TD) may be used to better characterize and understand the interactions of these two chemicals. In this study, a previously published generic physiologically-based toxicokinetic (PBTK) model for fish was adapted to ATR and CPF. A sub-model of the kinetics of one of the main metabolites of CPF, chlorpyrifosoxon (CPF-oxon), was included, as well as a TD model. Inhibition of two esterases, AChE and carboxylesterase, by ATR, CPF and CPF-oxon, was modeled using TD modeling of quantities of total and inactive esterases. Specific attention was given to the parameterization and calibration of the model to accurately predict the concentration and effects observed in the fish using Bayesian inference and published data from fathead minnow (Pimephales promelas), zebrafish (Danio rerio) and common carp (Cyprinus carpio L.). A PBTK-TD for mixtures was used to predict doseresponse relationships for comparison with available adult fish data. Synergistic effects of a joint exposure to ATR and CPF could not be demonstrated in adult fish.

Introduction

In environmental biomonitoring programs, biomarkers are increasingly proposed as tool to detect early exposure to xenobiotics or physiological effects [1][2][3] . Among biomarkers, neuromuscular parameters have proven their interest in the study of the effects of contaminants in the environment 4 . The inhibition of those parameters by pesticides may cause severe damage to the organism. In particular, the inhibition of acetylcholinesterase (AChE), which is involved in the regulation of acetylcholine action in neuromuscular junction, is responsible for various effects, e.g. the impairment of swimming behavior 5 or the ability to feed, and may lead to the death of the organism 6 . AChE is also one of the most studied biomarkers in fish [7][8][9] . Other enzymes can also be inhibited by pesticides, such as carboxylesterases (CEs), which may not cause apparent direct toxicity to the organism 6,10 but may be involved in a detoxification process. CEs can indeed bind to pesticides and thus act as a "sink" by decreasing the amount of effective pesticide 6,11 . Atrazine (ATR) and chlorpyrifos (CPF) are among the most commonly used pesticides in the world 12 . CPF is an organophosphate (OP) insecticide that is used in both agriculture and in residential settings 13 . CPF affects the neuromuscular junction in target species by inhibiting AChE 14 , both directly, as a parent chemical, and indirectly, via one of its metabolites, CPFoxon. This metabolite is reported to be a stronger inhibitor of AChE and it has been suggested, that the effect of CPF may be mainly due to its metabolite [14][15][16] . ATR is a herbicide that belongs to the triazine family and inhibits electron transport mechanisms of photosystem II in target plants 14 . While ATR is not an AChE inhibitor since it cannot interact at the serine residue of the active site, a decrease of AChE activity is measured in fish 17,18 after exposure to ATR alone. For this reason, and also because of suspicion of endocrine disrupting effects, ATR was banned by the EU in 2003 19 . 10 years after the ban, ATR was still responsible over 70% of the undrinkable water in particular because some products of atrazine degradation are slow to disappear 20 . However, levels of atrazine dropped significantly in coastal water from more than 150 to 7 ng/L between 1991 and 2010 21 . As for CPF, the EU prohibited its use by the nonrenewal of its marketing authorization in 2020 because of its potential genotoxicity and developmental neurotoxicity 22 . Nevertheless, the mechanisms of low-dose effects of ATR and CPF as single chemicals and in mixture are still unclear.

Aquatic organisms are exposed to multiple chemicals in the aquatic environment. Certain chemicals may interact with each other and cause greater (synergism) or smaller effects (antagonism) than expected based either on the addition of concentrations of the single chemicals 23 or on independent effects 24 . A greater than additive toxicity of ATR and CPF in mixture was shown in invertebrates 14,16,17 and was also observed in zebrafish larvae (Danio rerio) and fathead minnow juveniles (Pimephales promelas) 5,25 . However, similar studies in adult fish, such as bluegill (Lepomis macrochirus), fathead minnow and common carp (Cyprinus carpio L.), were not conclusive about the synergistic effect on AChE inhibition 26,27 . Synergistic effects were suggested to be the result of an increase in synthesis of cytochrome P450 (CYP450) enzymes due to the presence of ATR 28 . This increase would lead to a higher metabolization of CPF into its two metabolites, chlorpyrifos-oxon (CPF-oxon) and trichloropyridinole (TCP) (Figure 1) in carp and zebrafish 26,27 .

Physiologically-based toxicokinetic models (PBTK) are useful tools to improve understanding of the fate of a chemical or a mixture inside an organism [29][30][31] . Recently, Grech, et al. 32 proposed a generic PBTK model for four different fish, including zebrafish, threespine stickleback (Gasterosteus aculeatus), rainbow trout (Oncorhynchus mykiss) and fathead minnow. The construction of a PBTK model with a toxicodynamic (TD) component that predicts effects at the target organs can improve the understanding of toxicity mechanisms, by testing various hypotheses, and be relevant for risk assessment 10,33,34 . Moreover, PBTK-TD models for mixtures can help understand mixture effects by quantifying the contribution of kinetics and dynamics to any possible toxicological interactions. Such tools have been applied in rats for several mixtures 33,[35][36][37] . In fish, only two models have been developed so far 10,29 . One of them successfully described interaction between melamine and cyanuric acid. It was shown to be mostly due to the chemical interaction at the target organ and demonstrated the usefulness of such models 29 . This paper aims to give new insight on the toxicokinetics and toxicodynamics of ATR and CPF by using both in silico methods and experimental data from the literature.

Materials and methods

We extended the PBTK model developed by Grech, et al. 32 in order to account for the kinetics of both ATR, CPF, and CPF-oxon, which is the most toxic metabolite of CPF. The toxicodynamic part of the model represented inhibition of AChE and CEs by ATR, CPF and CPF-oxon and was based on the equations previously proposed 10,33,34 in rainbow trout and rat. Model parameters were calibrated using experimental data relative to internal kinetics in several fish species and to enzyme inhibition in the common carp. The PBTK-TD model was then adapted to predict effects on esterase inhibition of a mixture of CPF and ATR with the aim of assessing potential induction of CPF metabolism by ATR.

Experimental data

A full description of the data used can be found in SI section 2.1. TK and TD data obtained in vivo in fish exposed to ATR, CPF, or both was collected from the literature in a total of 12 publications (see Table S1 in SI). TD dataset reported in vivo AChE and CEs inhibition in a context of exposure to ATR, CPF or mixture of both. Among these studies, Xing, et al. 17 was selected to develop our model, because they provided a complete description of the exposure scenarios as well as a study of both ATR and CPF, and a mixture of both, with the same experimental design. In addition to in vivo data, in vitro data from the literature was used for the TD sub-model (Table S2 and S3 in SI). The data had been obtained on mice, rat and trout 6, 10, 38, 39 .

Model structure

The PBTK model used is based on the generic model developed by Grech, et al. 32 . This model has been successfully applied to four different species including rainbow trout, zebrafish, fathead minnow and threespine stickleback. The PBTK model was extended to model the kinetics of both CPF and its toxic metabolite, CPF-oxon, and the kinetics of ATR.

Absorption was assumed to be branchial due to the lack of information regarding gastrointestinal absorption of ATR and CPF and because the data on kinetics was obtained using waterborne exposure to ATR of CPF. Excretion was assumed to be mediated by gills, and/or bile. Metabolization was assumed to take place only in liver and was described using a Michaelis-Menten equation (CPF and its metabolites, CPF-oxon and TCP) or first order clearance (ATR) (see Figure S1).

Equations modeling the toxicodynamics of AChE and CEs inhibition were added to the PBTK based on the work of Timchalk and Poet 33 and Abbas and Hayton 10 , on rat and trout, respectively, initially developed for organophosphates. Although, AChE inhibition by ATR was probably the result of an indirect interaction between ATR and AChE, the same equation was used to describe this process. In this case, the equation is used as an empirical description of the effects, of ATR on AChE [START_REF] Brodeur | The problem with implementing fish farms in agricultural regions: A trial in a pampean pond highlights potential risks to both human and fish health[END_REF] . Indeed, inhibition could occur before AChE protein level, maybe at the ache gene transcriptional level [START_REF] Wang | Effects of Atrazine on the Development of Neural System of Zebrafish, Danio rerio[END_REF] .

Synthesis of each enzyme was modeled in brain, muscle, and liver with a zero-order synthesis rate, Ks (nmol.d -1 ) and a first-order degradation of esterases with Kd (d -1 ), with organspecific values. Binding of CPF and CPF-oxon to AChE and CEs was assumed to be irreversible, as reactivation of inhibited AChE was assumed to be insignificant [START_REF] Wallace | Reactivation and aging of phosphorylated brain acetylcholinesterase from fish and rodents[END_REF] .The indirect inactivation of AChE by ATR was assumed to be irreversible as well. Inhibition of both esterases was modeled using a bimolecular inhibition rate constant Ki (nmol -1 .d -1 ) (Eq. 1).

𝑑𝑄 𝑒𝑠𝑡𝑒𝑟𝑎𝑠𝑒 = 𝐾 𝑠 -𝑄 𝑒𝑠𝑡𝑒𝑟𝑎𝑠𝑒 × (𝐾 𝑑 + (𝐾 𝑖 × 𝑄 𝑐ℎ𝑒𝑚𝑖𝑐𝑎𝑙 )) (1)
With Qesterase the quantity of unbound (CPF/ CPF-oxon) or inactivated (ATR) esterase in the tissue (nmol) and Qchemical the quantity of active CPF, CPF-oxon, or ATR in the tissue (nmol).

Model parameterization

Most of physiological parameters of the model for rainbow trout and fathead minnow had been collected by Grech, et al. 32 (see Table S4 in SI section 2.3). In the present paper, TD processes were modelled with additional physiological parameters, for example, the enzyme levels in tissues. Model performance was increased by calibrating a small number of TK and TD parameters, that were determinant towards whole-body concentrations, using experimental data and Bayesian methods (Monte Carlo Markov Chains, MCMC) (see SI section 2.4). TK and TD parameters were calibrated together since enzyme inactivation by the chemicals affects the kinetics of ATR, CPF, and CPF-oxon.

Toxicokinetics

Chemical-specific TK parameters from the literature were obtained from the literature or using QSAR models by Grech, et al. 32 and listed in Table 1 and2

. Parameter values for

Michaellis-Menten metabolism (Km and Vmax) of CPF into CPF-oxon and TCP were obtained in vitro [START_REF] Lavado | Microsomal biotransformation of chlorpyrifos, parathion and fenthion in rainbow trout (Oncorhynchus mykiss) and coho salmon (Oncorhynchus kisutch): Mechanistic insights into interspecific differences in toxicity[END_REF] . Partition coefficients for CPF and CPF-oxon were estimated using a QSAR method based on logKow (CPF-oxon logKow =3.5, CPF logKow =4.96, according to PubChem) and tissue composition 32,[START_REF] Pery | A Physiologically Based Toxicokinetic Model for the Zebrafish Danio rerio[END_REF] . Kinetics of CPF-oxon metabolism were assumed to be the same as the kinetics of CPF into TCP.

Calibration of the unbound fraction in blood (UF) and the blood:water partition coefficient (PCBW) was based on whole-body concentrations resulting from continuous exposure to CPF reported in fathead minnow in two different studies 15,16 , and on the quantity of non-metabolized CPF excreted by the gills which was reported to be approximately 2% of total amount excreted in guppies (Poecilia reticulata) [START_REF] Welling | Bioconcentration kinetics of the organophosphorus insecticide chlorpyrifos in guppies (Poecilia reticulata)[END_REF] (for details see SI section 2.4). Prior distribution of PCBW was based on the QSAR estimate. Tissue:blood partition coefficients were estimated using both a QSAR model and the best fit estimate of PCBW. Since no data on CPF-oxon kinetics was available, including in other species, parameter values for CPF-oxon were not calibrated.

ATR tissue:blood partition coefficients were set to in vivo measurements in whitefish (Coregonus fera) [START_REF] Gunkel | Mechanisms of bioaccumulation of a herbicide (atrazine, s-triazine) in a fresh-water mollusk (Ancylus-fluviatilis Mull) and a fish (Coregonus-fera Jurine)[END_REF] . In addition, due to the lack of data in fish, the UF of ATR in blood was set to 0.74 as reported in rat [START_REF] Lu | Salivary concentrations of atrazine reflect free atrazine plasma levels in rats[END_REF] . We considered an ATR hepatic clearance of 0.024 mL/d/g liver determined in rainbow trout based on the negligible in vitro clearance reported by Han, et al.

48

. The blood:water partition coefficient (PCBW), biliary excretion rate (Kebile) and transfer rate from bile bladder to the gastrointestinal tract lumen (KBG) were calibrated using the PBTK model for zebrafish and whole body concentrations in juvenile zebrafish during exposure and depuration phases (Table S1 in SI) [START_REF] Gorge | Kinetics and metabolism of 14c-lindane and 14c-atrazine in early life stages of zebrafish (Brachydanio-rerio)[END_REF] . Prior distributions for KBG and Kebile were noninformative; PCBW prior distribution was based on the value measured in vivo in whitefish (Coregonus fera) [START_REF] Gunkel | Mechanisms of bioaccumulation of a herbicide (atrazine, s-triazine) in a fresh-water mollusk (Ancylus-fluviatilis Mull) and a fish (Coregonus-fera Jurine)[END_REF] .

Toxicodynamics

The toxicodynamic parameters that regulate the dynamics of AChE and CEs at steady state were either collected from the literature (Table S2 and S3 in SI section 2.2) or calibrated.

The enzyme inhibition rates collected from experimental studies and used for calibration were arcsine transformed.

The degradation rate constants (Kd) of AChE in brain, muscle and liver were set to the values estimated using a dose-response modeling approach in rainbow trout [START_REF] Abbas | A physiologically based pharmacokinetic and pharmacodynamic model for paraoxon in rainbow trout[END_REF] , i.e. 1.75x10 - 3 d -1 , 1.4x10 -2 d -1 , and 1.12x10 -1 d -1 in brain, muscle, and liver respectively. The degradation rate constants of CEs were assumed to be of the same order of magnitude and were therefore set to the same values as AChE, as in the models developed in trout 10,[START_REF] Abbas | A physiologically based pharmacokinetic and pharmacodynamic model for paraoxon in rainbow trout[END_REF] and in rat 34 .

Initial quantities of AChE and CEs in brain, muscle and liver were calibrated, as well as the chemical-specific inhibition rates (Ki) of AChE and CEs for ATR, CPF, and CPF-oxon.

Calibration was based on in vivo inhibition of AChE and CEs activity in muscle and brain measured in dose-response experiments in the common carp (Cyprinus carpio L.) exposed to either ATR or CPF 17 . The PBTK model for fathead minnow was used, since both species are Cyprinids. Non-informative prior distributions were used for the Ki. Prior distributions of the initial quantities of CEs (nmol/g BW) were based on the values obtained in rats in nmol [START_REF] Maxwell | The effects of blood flow and detoxification on in vivo cholinesterase inhibition by Soman in rats[END_REF] , and transposed to fish by using the relative organ weights in rats, the relative organ weights in fish already used in our PBTK model, and the protein content in rats and fish. Prior distributions of the initial quantities of AChE were in turn estimated assuming that the ratio between AChE and CEs levels was equal to the ratio reported in rats [START_REF] Maxwell | The effects of blood flow and detoxification on in vivo cholinesterase inhibition by Soman in rats[END_REF] (for details see SI section 2.2).

In order to maintain steady state of enzymes over long periods of simulation time, the synthesis rates (Ks) were determined according to Eq. 2 for AChE and CEs in each organ.

𝐾 𝑆 = 𝑄 𝑒𝑠𝑡𝑒𝑟𝑎𝑠𝑒 𝑖𝑛𝑖𝑡𝑖𝑎𝑙 × 𝐾 𝑑 (2)

Model applications

A sensitivity analysis of AChE and CEs enzyme inhibition in brain and muscle was performed, according to the exposure scenario used with carps after 40-day exposure and after 20-day depuration 17 , and using the variance-based Sobol method [START_REF] Saltelli | Sensitivity Analysis[END_REF][START_REF] Sobol | Estimating the approximation errors when fixing unessential factors in global sensitivity analysis[END_REF] . 65 chemicalspecific parameters were studied with uniform distributions ±10% of the calibrated parameter values (see SI section 3 for details)

The calibrated PBTK-TD model was used to check how the TD (i.e. enzyme inactivation) affected internal concentrations of the xenobiotics. Furthermore, the protective role of CEs was investigated by predicting AChE inhibition dose-response curves in brain and muscle with and without modeling CEs inhibition.

Mixture effects on AChE and CEs inhibition were predicted without considering induction of CYP450 by ATR, and were compared to dose-response data obtained in the common carp 17 in order to assess the weight of this mechanism. The PBTK-TD model was used to predicted dose-response curves for equitoxic mixtures where each chemical contributes equally to the total toxicity. These dose-response curves were compared to those obtained using the PBTK-TD for single chemicals and Berenbaum's general solution for concentration addition 23 .

Software

Calculations were performed using R version 3.6.1 [START_REF] Team | A Language and Environment for Statistical Computing[END_REF] , with packages deSolve [START_REF] Soetaert | Solving Differential Equations in R: Package deSolve[END_REF] , sensitivity 56 , and GNU MCSim v6.2.0 57 .

Results

Calibration of toxicokinetic and toxicodynamic parameters

Estimates of chemical-specific kinetic and toxicodynamic parameters are reported in Table 2; estimates of initial enzyme quantities are reported in Table 3.

Calibration of ATR parameters mainly improved the elimination kinetics. The estimated PCBW was close to the mode of the prior distribution (prior values in Table S5). Predicted wholebody concentrations in the calibration dataset, obtained in juvenile zebrafish, were relatively accurate: all predictions were within a two-fold factor except for the first uptake timepoint (Figure 2).

The estimate of PCBW for CPF was within a two-fold factor of the QSAR prediction which was used as a prior, with a wide 95% confidence interval ranging from 1.8x10 3 to 7.91x10 3 .

Although the prior distribution was non-informative, the estimate of CPF UF (0.0195) was close to the value reported in Poecilia reticulata [START_REF] Welling | Bioconcentration kinetics of the organophosphorus insecticide chlorpyrifos in guppies (Poecilia reticulata)[END_REF] . As presented in the Figure 3, eight out of nine predicted whole-body concentrations were within a 3-fold factor and all were within a 10-fold factor. Predicted CPF concentrations were more accurate for Jarvinen, et al. 15 than Mehler, et al. 16 . The main difference between those two experiments was the longer exposure duration in the first one (Table S1).

The predicted internal concentrations after 40 days continuous exposure to ATR or CPF illustrate the larger bioconcentration factor of CPF: at equal exposure doses: whole body concentrations of CPF were 400-fold greater than ATR concentration levels and 30-fold greater than CPF-oxon (Table S10).

Estimated initial enzyme quantities were mostly close to the prior distributions based on data reported in rat, except for the estimated initial CEs quantity in muscle which was larger than the prior (estimates reported in Table 3 and prior distributions in Table S5)

Bimolecular rate constants (Ki) for both AChE and CEs were calibrated with noninformative prior distributions (prior in Table S5, and posterior in Table 2). The estimated bimolecular rate between ACHE and CPF-oxon (Ki AChE/CPF-oxon) was larger than Ki AChE/ATR by two orders of magnitude. The estimated CEs bimolecular rate was the largest for CPF-oxon and was of the same order of magnitude as Ki AChE/CPF-oxon. The estimated Ki CEs/CPF was particularly low. The resulting differences in responses at equivalent internal doses, which take into account the amount of each enzyme in tissues, are illustrated in Figure S13 in SI: inhibition is strongest in muscle, with comparable toxicity of ATR and CPF-oxon.

Comparison of the predicted esterase activity and observations 17 used as a calibration dataset showed that the slope of the predicted dose-responses was systematically too steep (Figure 4). Depuration was particularly badly predicted in brain: predicted activities of both AChE and CEs were almost not affected by the 20-day-depuration whereas observed data showed a clear increase in activity indicating recovery.

Single chemical model evaluation

The sensitivity analysis (SI section 3) of outputs related to AChE and CEs inhibition showed varying parameter influence depending on the chemical and the enzyme. Interestingly, some TK parameters for ATR, in particular unbound fraction, brain/muscle PC and PCBW, were

Model evaluation in mixture

Enzyme inhibition resulting from exposure to the mixtures of ATR and CPF described in Xing, et al. 17 was predicted using the PBTK-TD calibrated for single chemicals and do not account for CPF metabolism induction by ATR. The predictions were close to observations with a maximum discrepancy of 20% between observations and predictions, in particular for the highest dose (Figure 5), but as for the single chemicals, the slope of the dose-response was too steep.

The PBTK-TD model was used to predict dose-response curves for equitoxic mixtures.

Doses-response curves were also obtained using the dose-responses predicted using the PBTK-TD for single chemicals and the assumption of concentration addition 23 . The results confirm that the PBTK-TD model for mixtures prediction do not diverge significantly from the dose-responses predicted under the assumption of concentration addition (see Figure S13 in section 4.3 in SI). Minor deviations were observed at high exposure levels in muscle for CEs, and in brain for AChE.

Discussion

A PBTK-TD model in fish was successfully developed to predict the uptake and internal kinetics of two pesticides, atrazine and chlorpyrifos, and their impacts on acetylcholinesterase and carboxylesterase inhibition. This model was evaluated by using an external dataset where carps were exposed to mixtures of ATR and CPF. Although the datasets used for calibration were obtained by exposure to single chemicals, the model correctly predicted the observed inhibition caused by exposure to the mixture. Moreover, the dose-response relationship predicted by the PBTK-TD model in fathead minnow was satisfactory, although the concentration addition model used was the most parsimonious; synergistic effects of a joint exposure to ATR and CPF were not demonstrated.

Physiologically based toxicokinetic models coupled to a toxicodynamic model are still scarce in ecotoxicology. In particular, for fish, only three models have been developed so far, two in rainbow trout 10,29 and one on farmed tilapia (Oreochromis mossambicus) [START_REF] Ling | A PBTK/TD Modeling-Based Approach Can Assess Arsenic Bioaccumulation in Farmed Tilapia (Oreochromis mossambicus) and Human Health Risks[END_REF] . Two of these studies modeled effects of single chemicals, paraoxon and arsenic, whereas Tebby, et al. 29 simulated an exposure to a mixture of melamine and cyanuric acid. Because they are based on the physiology, PBTK models are interesting tools for taking into account interspecies variability. Various datasets describing toxicokinetics are available on species that were not included in the generic model in the first place, such as carp or medaka (Oryzias latipes). As shown in the present paper, certain species can be assumed to be similar, or else species-specific parameters can be modified to adapt the model to extra species.

Although they are still rare, the benefit of such PBTK-TD models in mixture is twofold: being able to propose more realistic exposure scenarios since aquatic organisms face multiple, timedependent exposures, and also being able to link internal concentrations in organs to various effects. Currently, environmental risk assessment promotes a different, more generic approach based on an empirical model, GUTS (General Unified Threshold model) [START_REF] Products | Scientific Opinion on the state of the art of Toxicokinetic/Toxicodynamic (TKTD) effect models for regulatory risk assessment of pesticides for aquatic organisms[END_REF] . Though this approach is promising by gathering under the same framework several TK-TD models, they are empirical models and fail to explain the mechanism of action of toxic: effects are only predicted on high-level functions, mainly survival. PBTK-TD models offer the possibility to explain the mechanism of action of toxicants by modeling their effects at target organ. For this reason, they represent a promising development in the future of qAOPs (quantitative adverse outcome pathway) [START_REF] Conolly | Quantitative Adverse Outcome Pathways and Their Application to Predictive Toxicology[END_REF] .

However, the requirements in experimental data while building such models are rarely met notably because experimental designs rarely consider both TK and TD. In addition, physiological data in fish are usually scarce and sometimes unreliable. For example, in our case, absorption was only considered through the gills whereas it maybe also occurs through the GIT. Likewise, complex phenomenon, such as blood-brain barrier, which would require a lot of data to be implemented in the generic model are not included although they could improve its predictive capability, notably in the brain. However, in the case of exposure to ATR and/or CPF, the availability of data in several species at different doses and at different organizational levels facilitated the construction of the first PBTK-TD model. Nonetheless, in general, experimental data does not include both variation on dose and time-course, which is a strong limitation for the domain of applicability.

A generic model developed by Grech, et al. 32 was used to predict the TK of ATR, CPF and CPF-oxon, a metabolite of CPF. Physiological parameters of the PBTK model, such as relative blood flows to organs, relative organ volumes, and tissue composition, were set to values as species-specific as possible. The chemical-specific parameterization in the model by Grech, et al. 32 was refined by fitting the model to the observations obtained in two fish species (zebrafish and fathead minnow) in three experimental studies 15,16,[START_REF] Gorge | Kinetics and metabolism of 14c-lindane and 14c-atrazine in early life stages of zebrafish (Brachydanio-rerio)[END_REF] . Although TK observations consisted in whole-body concentrations, the final predictions of esterase inhibition in specific organs were satisfactory, indicating that predicted internal concentrations in organs may be accurate, or at least compensated by the values of the organ-specific TD parameters (initial enzyme quantities). However, using different fish species in the calibration may limit the generalization of the calibrated parameters and be a source of uncertainty. For example, ATR TK parameters were calibrated on a dataset on juvenile zebrafish and may only be valid at this developmental stage. Moreover, the scarcity of the studies on the mixture of ATR and CPF lead us to use a dataset obtained in carps, which was not one of the four model species in the generic PBTK. As common carp and fathead minnow were both Cyprinids, we decided to use fathead minnow model. We assumed that though carps were larger than fathead minnows, their physiological parameters were similar, e.g. relative blood flows to organs or the relative weight of organs. Nonetheless, strong variations in metabolism can be expected, as illustrated by the difference in metabolism of CPF into CPF-oxon between two salmonids [START_REF] Lavado | Microsomal biotransformation of chlorpyrifos, parathion and fenthion in rainbow trout (Oncorhynchus mykiss) and coho salmon (Oncorhynchus kisutch): Mechanistic insights into interspecific differences in toxicity[END_REF] .

Chemical-specific parameters of CPF-oxon were not calibrated due to lack of experimental data. CPF-oxon is generally not detected in tissues [START_REF] Barron | CHLORPYRIFOS PHARMACOKINETICS AND METABOLISM FOLLOWING INTRAVASCULAR AND DIETARY ADMINISTRATION IN CHANNEL CATFISH[END_REF] because of rapid hydrolysis to TCP [START_REF] Barron | Absorption, tissue distribution and metabolism of chlorpyrifos in channel catfish following waterborne exposure[END_REF] , as suggested by in vitro microsomal biotransformation rates in mouse [START_REF] Sultatos | Kinetic analyses of the microsomal biotransformation of the phosphorothioate insecticides chlorpyrifos and parathion[END_REF] . Consequently, partition coefficients were estimated using the QSAR method 32 and metabolism for CPF-oxon was chosen to follow the same kinetics as CPF into TCP. Although CPF metabolism has been quantified in vitro in trout hepatocytes, metabolism of CPF-oxon has not been quantified. The high sensitivity of the esterase inhibition to CPF and CPF-oxon metabolism rates revealed by the sensitivity analysis underlines that metabolism is an important source of uncertainty in our TK predictions, and could be reduced by producing new experimental data. Overall, predicted kinetics of both parent chemicals using a small number of calibrated parameters were satisfactory, since most were within a 3-fold factor.

The TD part of this model is based on the first PBTK-TD model in fish 10 . This study modeled the inhibition of CEs and AChE in rainbow trout in presence of paraoxon, an organophosphate related to CPF. However, contrary to Abbas and Hayton 10 who focused on time-courses of paraoxon internal concentrations and effects, TK and TD in our model were predicted and compared to observations at several doses but at only one or two timepoints in each study.

Although the same model structure for AChE inhibition can be used for several species, the TD parameter values may vary considerably from one species to another. Bimolecular inhibition rates are also highly dependent on the species. A 3-fold variation in Ki for methylparaoxon with AChE was reported in the same study in three fish species [START_REF] Freitas | Evaluation of a Brain Acetylcholinesterase Extraction Method and Kinetic Constants after Methyl-Paraoxon Inhibition in Three Brazilian Fish Species[END_REF] where the reaction was assumed to be reversible. In our model, as in other PBTK-TD models for AChE inhibition developed in trout and rodents, the Ki were assumed to be independent of substrate concentration. This may be a limitation of the model, since in mice, the efficiency of phosphorylation appeared to decrease with increasing substrate concentration [START_REF] Kardos | Interactions of the Organophosphates Paraoxon and Methyl Paraoxon with Mouse Brain Acetylcholinesterase[END_REF] .

Enzyme synthesis, degradation rates and concentrations in organs are clearly not sufficiently documented to provide species-specific TD models for AChE in fish. Most of the data reported on AChE inhibition is reported as percent inhibition compared to control, or in the best case as nanomoles of substrate hydrolysed per unit time and weight 11,[START_REF] Sturm | Inhibition of rainbow trout acetylcholinesterase by aqueous and suspended particle-associated organophosphorous insecticides[END_REF][START_REF] Carr | Kinetic analysis of the in vitro inhibition, aging, and reactivation of brain acetylcholinesterase from rat and channel catfish by paraoxon and chlorpyrifos-oxon[END_REF][START_REF] Barron | Tissue carboxylesterase activity of rainbow trout[END_REF][START_REF] Wogram | Effects of parathion on acetylcholinesterase, butyrylcholinesterase, and carboxylesterase in three-spined stickleback (Gasterosteus aculeatus) following short-term exposure[END_REF] .

Consequently, a 15-fold variation in brain and plasma AChE activity per tissue weight has been reported across 16 species [START_REF] Chuiko | Acetylcholinesterase and butyrylcholinesterase activities in brain and plasma of freshwater teleosts: cross-species and cross-family differences[END_REF] . Moreover, variability in AChE activity has even been observed between fish of the same species maintained either in laboratory or in situ [START_REF] Sturm | Inhibition of rainbow trout acetylcholinesterase by aqueous and suspended particle-associated organophosphorous insecticides[END_REF] . The strong interspecies variability in AChE activity may add up to variability in synthesis, degradation and inhibition rates and is likely to result in large differences between species.

Though the data from in vitro experimentations helps understand differences in toxicity between chemicals or between species, it could not be directly integrated into our model to predict in vivo inhibition due to the fact that many intra-cellular mechanisms were not modeled and do not occur in in vitro cell-free experimentations. For example, modeling the dynamics of ButyrylCholinEsterase (BuChE) binding processes during detoxification 6 could increase model realism and performance. This non-target B-esterase binds stoichiometrically with organophosphates explaining their inclusion in a PBTK-TD model of brain AChE inhibition in rat 38,[START_REF] Timchalk | Pharmacokinetic and pharmacodynamic interaction for a binary mixture of chlorpyrifos and diazinon in the rat[END_REF] . Although BuChE has been detected in the plasma of some Cyprinid species, it has not been detected in many other fish [START_REF] Chuiko | Acetylcholinesterase and butyrylcholinesterase activities in brain and plasma of freshwater teleosts: cross-species and cross-family differences[END_REF][START_REF] Küster | Cholin-and carboxylesterase activities in developing zebrafish embryos (Danio rerio) and their potential use for insecticide hazard assessment[END_REF] . The lack of knowledge regarding BuChE compared to AChE in fish hinders integration of BuChE inhibition in a PBTK-TD model: BuChE tissue levels, synthesis, and degradation rates are unknown.

In our study, inhibition was modelled in the main target organs, brain and muscle, and also in liver. Indeed, the liver was thought to be a large reservoir of CEs and act as a buffer compartment through binding of CEs to metabolites 6,10 . In Abbas and Hayton 10 , the toxicodynamics were described in heart, brain and liver. In addition, other compartments were identified as the place of esterase inhibition such as plasma or muscle in general 6 . Data was too scarce to be able to include plasma esterase levels in our model even though plasmatic CEs may play an important role. For example, in rainbow trout, plasmatic CEs represents as much as 58% of total body CEs activity in adult rainbow trout [START_REF] Barron | Tissue carboxylesterase activity of rainbow trout[END_REF] . Furthermore, as plasmatic CEs is distributed round the body, it may contribute more to detoxifying than liver CEs.

One other limit of the realism and performance of this model was the lack of specific synthesis rate values for CEs, as Abbas and Hayton 10 only proposed in vivo synthesis rate values for AChE. In addition, these zero-order synthesis rates imply that synthesis is independent of the actual quantity of enzyme. However, recent studies in adult red crucian carp (Carassius auratus) observed an increase of AChE levels after 6 days of exposure to fluoxetine and a backup to normal value after 6 days of recovery [START_REF] Pan | Tissue bioconcentration and effects of fluoxetine in zebrafish (Danio rerio) and red crucian cap (Carassius auratus) after short-term and longterm exposure[END_REF] . The lack of feedback loop in the model that would modulate AChE synthesis rate could explain the fact that the predicted dose-response relationships (SI section 4.3) were systematically too steep. This may reveal a shortcoming of the model since even increasing 100-fold the degradation constant rates (and thereon increasing also 100-fold the synthesis rate) did not produce satisfactory doseresponse slopes.

MCMC calibration on in vivo data allowed us to propose bimolecular rate constants of the binding between AChE or CEs and ATR, CPF or CPF-oxon. As highlighted in Coban, et al. [START_REF] Coban | Comparison of inhibition kinetics of several organophosphates, including some nerve agent surrogates, using human erythrocyte and rat and mouse brain acetylcholinesterase[END_REF] , differences in experimental protocols resulted in a large variability in Ki determined in vitro.

Thus, Bayesian calibration could provide closer estimates of the biological value of those inhibition rates. The calibrated Ki values were largest and of the same order of magnitude for CPF-oxon and ATR: KiAChE/CPF-oxon was 2-fold higher than KiAChE/ATR and three orders of magnitude higher than KiAChE/CPF and KiCEs/CPF. Our calibrated value for KiAChE/CPF-oxon, 35.4 nmol -

1

.min -1 , was far greater than the value reported for paraoxon, 1.06 ×10 -3 nmol -1 .min -1 10 .

According to the predicted internal concentrations in muscle resulting from 40 days exposure (Table S9 in SI) internal CPF levels are greater than CPF-oxon levels and two to three orders of magnitude higher than ATR levels. CPF has been suggested to be responsible for a large part of AChE inhibition although CPF-oxon is a stronger inhibitor 16,[START_REF] Barron | Absorption, tissue distribution and metabolism of chlorpyrifos in channel catfish following waterborne exposure[END_REF][START_REF] Straus | Inhibition of acetylcholinesterase and aliesterases of fingerling channel catfish by chlorpyrifos, parathion, and S,S,S-tributyl phosphorotrithioate (DEF)[END_REF] . Given the ratio of predicted internal concentrations and the potency of each chemical regarding enzyme inhibition (Figure S14 in SI, section 4.5), our findings suggest that CPF-oxon is responsible for more AChE inhibition than CPF. In addition, CPF-oxon may be particularly responsible for the depletion of CEs. ATR has been reported to be responsible for a weak inactivation of AChE 14,16 . But according to our model predictions, at equal internal concentrations, ATR may be as potent as CPF-oxon. The difference in sensitivity of AChE activity in vivo to ATR and CPF therefore appears to be due to the difference in internal concentrations. The strong potency of ATR was unexpected because inactivation of AChE may be the result of indirect toxicity in vivo as ATR is not an AChE inhibitor.

In this paper, the PBTK-TD model was built as a tool to investigate the effect of a mixture ATR-CPF on AChE and CEs inhibition. As AChE inhibition plays an important role in neuromuscular junctions and cholinergic synapses, AChE inhibition is a relevant biomarker of exposure 8,[START_REF] Küster | Cholin-and carboxylesterase activities in developing zebrafish embryos (Danio rerio) and their potential use for insecticide hazard assessment[END_REF] . As a detoxifying enzyme, CEs would play a role of protection by binding to pesticides and act as a stoichiometric buffer 10,17,[START_REF] Küster | Cholin-and carboxylesterase activities in developing zebrafish embryos (Danio rerio) and their potential use for insecticide hazard assessment[END_REF] . Calibrating the PBTK-TD model and the sensitivity analysis showed that the quantities of enzyme alone did not have a decisive role in the TD of esterase inhibition. Rather, the rate of renewal of the enzymes may play a larger role both in the EC50 of the toxic response and in the slope of the dose-response. The data available in carp provided dose-response relationships but information about the time-course of esterase inhibition is lacking. Chemical-independent data on enzyme renewal rates under control conditions would also greatly help to develop TD models.

Although an effect of ATR on CPF metabolism had been envisioned, interactions between ATR and CPF were finally not included in the model since they did not seem necessary. The mixture that had been tested in carp was equimassic, although ATR was roughly one order of magnitude less potent than CPF regarding AChE inhibition. Interactions would be more readily identifiable in an equitoxic mixture since interactions are unlikely to appear in an equimassic mixture of dissimilarly acting chemicals [START_REF] Backhaus | Predictability of the toxicity of a multiple mixture of dissimilarly acting chemicals to Vibrio fischeri[END_REF] . Nonetheless, because the actions of inhibition by the chemicals were affecting the same pool of esterases, it would have highlighted, if any, an interaction in mixture. It is the advantage of our model to include mechanisms of effects that can interact.

Our mechanistic model predicts AChE inhibition levels equivalent, at least at low doses, to the predictions based on dose-response modeling and concentration addition of external concentrations. This mechanistic model is able to predict sigmoidal dose-response relationships without use of an empirical dose-response function, with differences in shapes that reflect differences in esterase levels in tissues, renewal rates, and inactivation of two different esterases. [START_REF] Lavado | Microsomal biotransformation of chlorpyrifos, parathion and fenthion in rainbow trout (Oncorhynchus mykiss) and coho salmon (Oncorhynchus kisutch): Mechanistic insights into interspecific differences in toxicity[END_REF] , e-Grech, et al. 30 , f-Lavado and Schlenk 39 , g-Welling and de Vries [START_REF] Wang | Effects of Atrazine on the Development of Neural System of Zebrafish, Danio rerio[END_REF] . The rest of the data was calculated using QSAR model in Grech, et al. 30 . Parameters with an X can be found in Table 2. All values presented here were not calibrated and fixed in the model. 1.29x10 -2 (1.08x10 -2 ; 2.20x10 -2 ) 1.08x10 -4 (1.69x10 -5 ; 3.34x10 -4 ) 2.46x10 -2 (2.34x10 -2 ; 6.31x10 -2 ) Bimolecular rate constant CEs (nmol -1 .d -1 ) 7.08x10 -3 (5.59x10 -3 ; 1.23x10 -2 ) 3.14x10 -6 (2.90x10 -7 ; 2.31x10 -5 ) 2.42x10 -2 (1.19x10 -2 ; 2.59x10 -2 ) Table 3. Toxicodynamic parameter estimates (mode and 95% confidence interval)

Parameter Brain Muscle Liver Initial quantity of AChE (nmol.kg BW -1 ) 4.17x10 -4 (3.47x10 -4 ; 1.67x10 -3 ) 8.31x10 -3 (1.96x10 -3 ; 9.63x10 -3 )

1.37x10 -5 (6.96x10 -6 ; 3.38x10 -5 ) Initial quantity of CEs (nmol.kg BW -1 ) 1.19x10 -2 (4.96x10 -3 ; 2.46x10 -2 ) 97.1 (74.5; 110) 6.88x10 -1 (3.31x10 -1 ; 1.50) Supporting information to "Modelling acetylcholine esterase inhibition resulting from exposure to a mixture of atrazine and chlorpyrifos using a physiologically-based kinetic model in fish"

Model parameterization 2.1. Description of the TK and TD data available from the literature

An extensive search in literature was carried out using relevant keywords, and various combinations of those keywords: "atrazine", "chlorpyrifos", "fish", "toxicokinetic", "toxicodynamic", "acetylcholinesterase inhibition", "PBTK", "in vivo", "in vitro", "metaboli*", "excretion", "accumulation", "distribution", "absorption". TK and TD data obtained in vivo in fish exposed to ATR, CPF, or both was collected from the literature in a total of 12 publications (see Table S1 in SI). The data had been obtained in 8 different fish species. Three publications reported TK data only [3][4][5] . The nine others were reported in vivo AChE inhibition in a context of to CPF, ATR, a mixture of both or pre-exposure to ATR followed by CPF [6][7][8][9][10][11][12][13] . The internal concentrations of ATR, CPF or CPF-oxon were measured at single or multiple timepoints in whole body or in organs. Activity of AChE was determined in muscle, liver or brain. In addition, CEs activity was also measured in three studies. Among these studies, Xing et al. (2010) provided a complete description of the exposure scenarios as well as a study of both ATR and CPF, and a mixture of both, with the same experimental design. Consequently, this dataset was selected to develop our model. The other datasets were not used in the present study as these authors used very small juveniles, e.g. trout weight under 5 g in 3 on 4 publications, and they do not report data on the chemical TK. Therefore, these data were considered to be outside of our application domain. Indeed, even if the model can be used with data from juveniles, those juveniles were considered to be too small, i.e. too different from adults. In addition to in vivo data, in vitro data from the literature was used for the TD sub-model (Table S2 and S3 in SI). The data had been obtained on mice, rat and trout [14][15][16][17] .

Physiological parameters in the PBTK model

Physiological parameters for toxicokinetics in rainbow trout and fathead minnow were collected by 1 and are reported for fathead minnow.

MCMC Calibration

Parameters were calibrated using MonteCarlo Markov Chains, with 3 chains of 100000 iterations each. To avoid correlations between parameters, two toxicokinetic parameters were calibrated using an intermediate ratio. Calibration was carried out on those ratio (Eq.1

for CPF and Eq.2 for ATR). Convergence was assessed in R 3.6.1 34 with the package coda 35 by checking that autocorrelations were low (i.e. that the chains were well mixed), that estimates lay well within the prior boundaries, and that the Gelman-Rubin index 36 was close to 1.

UF = PC 𝑏𝑤 × 𝑅𝑎𝑡𝑖𝑜 𝑏𝑤_𝑈𝐹 (1)
K 𝐵𝐺 = Ke 𝑏𝑖𝑙𝑒 × 𝑅𝑎𝑡𝑖𝑜 𝐾𝑏𝑔_𝐾𝑒𝑏𝑖𝑙𝑒 (2)
Table S9: Estimates of physiological TD parameters when TK and TD were calibrated separately.

Parameter Brain Muscle Liver Initial quantity of AChE (nmol.kg BW - 1 ) 5.22x10 -4 (3.45x10 -4 , 1.68x10 -3 ) 3.51x10 -3 (1.96x10 -3 ,9.46x10 -3 ) 1.20x10 -5 (6.91x10 -6 , 3.41x10 -5 ) Initial quantity of CEs (nmol.kg BW -1 )

1.21x10 -2 (5.04x10 -3 , 2.42x10 -2 ) 125 (78,115) 6.82x10 -1 (3.29x10 -1 ,1.49)

Predicted dose-responses

The dose-response relationships predicted with the PBTK-TD model and resulting from exposure to ATR, CPF, or a mixture of both were compared with observations made by Xing et al in carp (Figure S11 and Figure S12). S10.

Comparison of binding to enzymes at equal internal concentrations

Enzyme activity was predicting by considering each enzyme independently and are not representative of the inhibition predicted in the tissues when both enzymes are modelled.

Effect of binding to CEs on AChE inhibition dose-response

The calibrated PBTK model was used to investigate the effect of binding to CEs on AChE inhibition resulting from exposure to a mixture of ATR and CPF (Figure S18). The slope of the dose-responses on a log-scale of concentrations was unchanged by the presence of CEs. Binding to enzymes is modelled in such a way that the presence of CEs decreases amount of CPF by a factor (around 10 in muscle, and 2 in brain) that is independent of the exposure concentrations. The effect was particularly strong in the muscle due to large amount of CEs estimated in muscles. A B

  among the ten most influential parameters affecting esterase inhibition (FigureS2to S5 in SI).Sensitivity analysis also highlighted the importance of CPF-oxon kinetics, including CPF-oxon formation by metabolism of CPF (Michaelis-Menten parameters, Vmax and Km, FigureS8and S9) on AChE and CEs inhibition. Among TD parameters, as expected, the bimolecular rates (Ki) of a given esterase and a given chemical were often in the three most influential parameters on the esterase inhibition regardless of the exposure, time and organ. However, CEs inhibition in muscle was less influenced by Ki than by Kd of CEs in muscle for ATR and CPF (FigureS5 and S9). The effect of TD on internal concentrations was assessed by also calibrating TK and TD parameters separately (see SI section 4.1 and 4.2). It appears that parameters from TK were the ones mainly affected by the calibration of TK and TD separately (see results of separate calibration in SI section 4.1 and 4.2). In particular, the chemical rate constant from bile to GIT lumen for ATR is much lower when the calibration in done separately (TableS7in SI). Additional simulations of kinetics with and without the TD part of the model showed that binding to enzymes was responsible for a three-fold decrease in whole body concentrations of CPF-oxon at low doses (up till doses eliciting around 60% AChE activity), an almost two-fold decrease in whole body concentrations of ATR, mostly due to inactivation in muscle, and only around 4% decrease in whole body concentrations of CPF at low doses (see SI section 4.6, FigureS15in SI). The proportion of active ATR or CPF decreases when exposure doses increase.Simulations using the PBTK-TD model without modeling binding to CEs illustrates the fact that binding to CEs decreases the amount of available CPF or CPF-oxon by a 10-fold factor in muscle and by a 2-fold factor in brain, therefore contributing to a considerable part in observed toxicity (see SI section 4.7, FiguresS16 to S18). Binding to CEs causes a 300-fold decrease in ATR levels in muscle and no significant change in ATR concentrations in brain.

Figure 1 .

 1 Figure 1. Potential relation between effects of atrazine and chlorpyrifos. Potential synergistic effect is symbolized in red.

Figure 2 .

 2 Figure 2.Predicted atrazine concentrations in the whole zebrafish compared to the observation from Gorge and Nagel. The solid line represents the model predictions, dots represent the mean concentration in two juveniles (n=47) and the grey area is the 95% prediction interval, computed from the posterior distributions. These simulations were made using a sample based on every iteration of the last 3333 of each of the three MCMC chains.

Figure 3 .

 3 Figure 3. Comparison between chlorpyrifos concentrations measured in fathead minnow (µg/g fish) and model predictions. Dotted lines represent the 3-fold and 10-fold changes. Experimental datasets used are indicated in legend. Data correspond to whole body concentration. Credibility intervals were too narrow to be represented.

Figure 4 .

 4 Figure 4. Comparison between esterase inhibition (acetylcholine esterase in black and carboxylesterase in grey) reported by Xing et al. (2010) in common carp versus model predictions for 40-day-exposition (filled points) to (a, b) atrazine and to (c, d) chlorpyrifos followed by a 20-day-depuration (empty points). The inhibition in the brain, respectively in the muscle, are presented in (a) and (c), respectively (b) and (d). Solid grey line represents the identity line and 95% credibility intervals are represented, computed from a sample based on the posterior distribution.

Figure 5 .

 5 Figure 5. Comparison between esterase inhibition (acetylcholine esterase in black and carboxylesterase in grey) reported by Xing et al. (2010) in common carp versus model predictions for 40-day-exposition (filled points) to a mixture of atrazine and chlorpyrifos followed by a 20-day-depuration (empty points). The inhibition in the brain, respectively in the muscle, is presented in panel (a), respectively (b). Solid grey line represents the identity line and 95% credibility intervals are represented, computed from a sample based on the posterior distribution.

Figure 1 .

 1 Figure 1. Potential relation between effects of atrazine and chlorpyrifos. Potential synergistic effect is symbolized in red.

Figure 2 .

 2 Figure 2. Predicted atrazine concentrations in the whole zebrafish compared to the observation from Gorge and Nagel. The solid line represents the model predictions, dots represent the mean concentration in two juveniles (n=47) and the grey area is the 95% prediction interval, computed from the distributions. These simulations were made using a sample based on every iteration of the last 3333 of each of the three MCMC chains.

Figure 3 .

 3 Figure 3. Comparison between chlorpyrifos concentrations measured in fathead minnow (µg/g fish) and model predictions. Dotted lines represent the 3-fold and 10-fold changes. Experimental datasets used are indicated in legend. Data correspond to whole body concentration. Credibility intervals were too narrow to be represented.

Figure 4 .

 4 Figure 4. Comparison between esterase inhibition (acetylcholine esterase in black and carboxylesterase in grey) reported by Xing et al. (2010) in common carp versus model predictions for 40-day-exposition (filled points) to (a, b) atrazine and to (c, d) chlorpyrifos followed by a 20-day-depuration (empty points). The inhibition in the brain, respectively in the muscle, are presented in (a) and (c), respectively (b) and (d). Solid grey line represents the identity line and 95% credibility intervals are represented, computed from a sample based on the posterior distribution.

Figure 5 .

 5 Figure 5. Comparison between esterase inhibition (acetylcholine esterase in black and carboxylesterase in grey) reported by Xing et al. (2010) in common carp versus model predictions for 40-day-exposition (filled points) to a mixture of atrazine and chlorpyrifos followed by a 20-day-depuration (empty points). The inhibition in the brain, respectively in the muscle, is presented in panel (a), respectively (b). Solid grey line represents the identity line and 95% credibility intervals are represented, computed from a sample based on the posterior distribution.

Figure S1 .

 S1 Figure S1. Schematic description of the PBK model developed for rainbow trout, zebrafish, fathead minnow, and stickleback. Metabolism, uptake and excretion sites are represented in purple, blue, and green, respectively. The PBTK model used is based on the generic model developed by Grech, et al. 1 . This model has been successfully applied to four different species including rainbow trout, zebrafish, fathead minnow and threespine stickleback. It comprises twelve well-mixed blood flow limited compartments:

Figure S11 :

 S11 Figure S11: Predicted and observed AChE activity in muscle (A) and in brain (B) after 40 days continuous exposure in water. Observations are in carp, predictions are obtained with the fathead minnow PBTK-TD model

Figure S13 : 4 . 4 .

 S1344 Figure S13: Predicted dose-responses obtained using the PBTK-TD model for mixtures, using and the concentration addition model based on the dose-responses predicted with the PBTK-TD model for single compounds, in brain (A, C) and muscle (B, D) for AChE (A, B) and CEs (C, D).

Figure S14 : 4 . 6 .

 S1446 Figure S14: Predicted AChE (A) and CEs (B) activity resulting from exposure to CPF, CPF-oxon, and ATR, as a function of internal concentrations after 40 days of continuous exposure to ATR or CPF, in brain and in poorly perfused tissues (muscle) predicted with the best-fit estimates. 4.6. Effect of binding to AChE and CEs on internal levels of CPF, ATR, and CPF-oxon

Figure S15 :

 S15 Figure S15: Predicted internal concentration levels of CPF (A) CPF-oxon (B), and ATR (C) that resulted from waterborne exposure to CPF or ATR, in muscle and in whole-body, with and without the TD submodel.

Figure S16 :

 S16 Figure S16: Predicted AChE inhibition in muscle (A), and in brain (B) resulting from continuous exposure to CPF for 40d, obtained with the PBTK-TD model (red), and with the PBTK-TD model without modelling CEs (green). The shaded areas represented the 95% prediction interval.

Figure S17 :

 S17 Figure S17: Predicted AChE inhibition in muscle (A), and in brain (B) resulting from continuous exposure to ATR for 40d, obtained with the PBTK-TD model (red), and with the PBTK-TD model without modelling CEs (green). The shaded areas represented the 95% prediction interval.

Figure S18 :

 S18 Figure S18: Predicted AChE inhibition in muscle (A), and in brain (B) resulting from continuous exposure to a mixture of ATR and CPF for 40d, obtained with the PBTK-TD model (red), and with the PBTK-TD model without modelling CEs (green). The shaded areas represented the 95% prediction interval.

Table lists Table 1 .

 lists1 Compound-specific parameters used in the PBTK model

Table 2 .

 2 Compound-specific parameter estimates (mode and 95% confidence interval)

Table 3 .

 3 Toxicodynamic parameter estimates (mode and 95% confidence interval)

Table 1 .

 1 Compound-specific parameters used in the PBTK model Data were retrieved from : a-Gunkel and Streit 42 , b-Han, et al.

[START_REF] Pery | A Physiologically Based Toxicokinetic Model for the Zebrafish Danio rerio[END_REF]

, c-Nichols, et al.

[START_REF] Gunkel | Mechanisms of bioaccumulation of a herbicide (atrazine, s-triazine) in a fresh-water mollusk (Ancylus-fluviatilis Mull) and a fish (Coregonus-fera Jurine)[END_REF]

, d-Lu, et al.

Table 2 .

 2 Compound-specific parameter estimates (mode and 95% confidence interval)

	Parameter	ATR	CPF	CPF-oxon
	PC liver/blood	5.35 a	5.27	2.59
	PC brain/blood	1.94 a	3.66	1.80
	PC kidney/blood	1.47 a	8.33	4.09
	PC blood/water	X	X	122.85
	PC richly perfused tissue/blood	1.47 a	3.40	1.67
	PC poorly perfused tissue/blood	1.18 a	1.25	0.62
	PC skin/blood	1.47 a	2.31	1.14
	PC gonads/blood	1.47 a	11.03	5.42
	PC fat/blood	1.47 a	47.63	23.36
	PC GIT/blood	5.47 a	3.71	1.83
	Clhepatic (mL/d/g liver)	0.024 b	_	_
	Vmax (µg/d/g)	_	3375 f (TCP) 1305 f (CPF-oxon)	3375
	Km (µg/mL)	_	50.3	50.3
	Ke bile (d -1 )	X	0 g	_
	Ke feces (d -1 )	0.83 c	0.83 c	0.83
	KBG (d -1 )	X	0	0
	Ratio blood to plasma	1	1	1
	Unbound fraction	0.74	X	1

Table S1 . Summary on experiment data on chlorpyrifos and atrazine in fish

 S1

	Xing, et al. 13 2.2	carpio Cyprinus	20	190	200	42.8, 4.28,			static ; 20d 40 d (semi-	measurements 2-point			Brain, muscle	CES AChE
						428			depuration)
							1.16,
	Study	Species	Fis	Initial	Aquari	Dose	Dose 11.6, 116	Exposure		Timepoints	of	Tissues	where	Tissues	where	Esterase
	reference		h	weigt	um	ATR 1.13, 11.3, 113 CPF	scenario		measurements	concentration	is	esterase inhibition is	inhibition
	Jarvinen, et al. 12 al. 18 Sandahl, et	Pimephale s promelas kisutch hus Oncorhync	per dos 2 e 15	0.2 h (g) 0.7	41 volum e (L) 30	_ (µg/L )	2.68,1.21 ,0.63,0.2 ,2.50 ,1.875 5, 1.25 (µg/L) 0.00,0.62	60 (continuous) d 4 d (static)	End measurement measurement End	point point	Whole body measured		measured Brain, muscle	_ AChE
	Mehler, et al. 6 Wacksman, et al. 9 Wheelock,	ha tshawytsc hus Oncorhync	30 30 10	0.05 0.5 4	4 9.5 4	_ 0,	7, 0.12 339,201, 105,53, 25 7.3 0, 1.2,	ATR for 2 d 2 d (static) 4 d	End End measurement measurement End	point point point	_ Whole body		Brain Whole body Brain, muscle, liver	AChE AChE CES AChE
	et al. 8					1000			then CPF 2	measurement
						1000 1.03,		d					_
							8.87,		(semi-static;
							47.14	pulse-dose)
	Sandahl	Oncorhync	50	3.8	1.5		0.00,0.62	4 d (static)		End	point	_		Brain	AChE
	and Jenkins	hus					5,	1.25				measurement
	7	mykiss					,1.875
							,2.50
	Sturm, et al.		10	14.6	30		0.1, 0.3,	1		h	End	point	_		Brain, muscle	AChE
	10						1, 3.3	(continuous)	measurement
	Gorge and	Danio rerio 150 0.02	0.5	135			1	d,2	d	47-point		Whole	body,	_
	Nagel 5								(static)		measurements	repeated
	Schmidel, et		12	NS	3	0,			14		d	End	point			Brain, muscle	AChE
	al. 11					10,			(continuous)	measurement
						1000
	Barron, et	Ictalurus	21	300	40		12		4 d (static)		4-point				Whole body	CES AChE
	al. 4	punctatus										measurements
	Gunkel and	Coregonus	NS	4	300	50			1 d (static)		End	point	Blood, muscle, gall	_
	Streit 3	fera										measurement		bladder,	liver,
														stomach, intestine,
														brain, remainder

. Toxicodynamic parameter values collected from the literature Table S2. Toxicodynamic parameter values for enzyme synthesis and degradation collected from the literature. Source and value for the different parameters are listed. Only the degradation rates of AChE in rainbow trout were used in the model.

	Enzyme	Unit	brain	muscle	liver	Reference
	AChE synthesis rate in	nmol/d	11.0 a	18.9 a	7.95 b	a16 b19
	rainbow trout
	AChE synthesis rate in		0.0336			20, 21
	rat
	AChE degradation rate in	(/d)	1.75x10 -3 a 1.40x10 -2 a 0.112 b	a16 b 19
	rainbow trout
	AChE degradation rate in		0.24 c		0.24 d	c21, 22 d14, 23
	rat

Table S3 . Data used for the toxicodynamic parameterization of initial enzyme quantities

 S3

					Brain	Muscle	Liver	Reference
	Initial CEs in 250g rats (nmol)	1.60	289*	460	24, 25 . * or 65 depending on the publication
	Initial AChE in 250g rats	0.11	0.843	0.009
	(nmol)
	Relative organ weight in rat		0.0057	0.404	0.0366	26
	Relative organ weight PBTK	0.0112	0.673	0.0218
	model for fathead minnow
	Protein	content	in	rat	104 a	164 b	210 c	a27
	(mg/g tissue)			55.3 g			b 28
					60.4 h			c 29
								g 17
								h 30
	Protein	content	in	fish	91.7 d	142 e	150 f	d shortspine thornyhead (Sebastolobus
	(mg/g tissue)			39.8 i			alascanus) 31
					41.5 j			e rainbow trout 32
					38.9 k			f rainbow trout 33
					40.6 l			i trout 17
					81 m			j trout 30
								k fathead minnow 17
								l fathead minnow 30
								m Scorpaena guttata 31
	Initial CEs in fathead minnow	0.0111	1.48	0.783	calculated
	(nmol/g BW)
	Initial CEs in rainbow trout	0.00485	1.66	0.524	calculated
	(nmol/g BW)
	Initial AChE in fathead minnow	7.63x10 -4	0.00432	1.53x10 -5 calculated
	(nmol/g BW)
	initial AChE in rainbow trout	32.9	8.7	2.89	16, 19
	(nmol/min/mg prot)
	initial AChE in rainbow trout	14.8	831	6.33	16, 19 calculated using data from Abbas and
	(nmol/g BW)						protein content and relative organ weight

Table S4 : Physiological parameters used in the PBTK model for fathead minnow, as proposed by 1 . Parameters (fraction) value

 S4

	Relative weight
	Adipose tissues	0.02
	Blood	0.02
	Brain	0.01
	GIT	0.10
	Gonads	0.01
	Kidney	0.002
	Liver	0.02
	PPT	0.67
	RPT	0.03
	Skin	0.10
	Relative blood flow
	Adipose tissues	0.01
	Brain	0.04
	GIT	0.17
	Gonads	0.01
	Kidney	0.02
	Liver	0.02
	PPT	0.54
	RPT	0.13
	Skin	0.06
	PPT blood flow to kidney (αFpp) 0.60
	Skin blood flow to kidney (αFs)	0.90
	Relative water content
	Adipose tissues	0.03
	Brain	0.75
	GIT	0.77
	Gonads	0.52
	Kidney	0.49
	Liver	0.65
	PPT	0.81
	RPT	0.53
	Skin	0.76
	Relative lipid content
	Adipose tissues	1.00
	Brain	0.07
	GIT	0.07
	Gonads	0.22
	Kidney	0.17
	Liver	0.10
	PPT	0.02
	RPT	0.07
	Skin	0.05

Table S10 : Predicted internal concentrations of CPF, CPF-oxon, and ATR resulting from a continuous exposure to CPF or to ATR at 11.3µg/L in water.

 S10

			Internal concentration (µg/g)
		Brain	Muscle	Whole body
	Chlorpyrifos	1.48	0.501	1.22
	Chlorpyrifos-oxon	0.0492	0.0169	0.0426
	Atrazine	0.00605	0.000159	0.003

Supporting information

Supporting Figure S1-S18 and Supporting Table S1-S10. On Figure S2, the convergence of the three chains can be checked visually, so as the little deviation between priors and posteriors when given normal and lognormal distribution.

Each chain contains the 50000 last values of each calibrated parameters.

Sensitivity analysis

A sensitivity analysis of AChE and CEs enzyme inhibition in brain and muscle was performed, according to the exposure scenario used with carps used by 13

Results of the calibration

Toxicokinetic calibration

Toxicodynamic calibration

Table S8: Estimates of chemical-specific TD parameters when TK and TD were calibrated separately.

Parameter ATR CPF CPF-oxon Bimolecular rate constant AChE (nmol -1 .d - 1 ) 1.58x10 -2 (1.02x10 -2 , 2.06x10 -2 ) 3.93x10 -5 (9.60x10 -6 ,1.23 x10 -3) 4.14x10 -2 (2.42x10 -2 , 6.33x10 -2 ) Bimolecular rate constant CEs (nmol -1 .d -1 ) 6.88x10 -3 (5.27x10 -3, 1.17x10 -2 ) 5.67x10 -6 (1.81x10 -7 , 1.10x10 -5 )

1.78x10 -2 (1.22x10 -2 , 2.65x10 -2 )