Frédéric Y Bois
email: frederic.bois@ineris.fr

Yves Frédéric

Bois

Physiologically-Based Modelling and Prediction of Drug Interactions

Physiologically based modeling and prediction of drug interactions

Introduction

The mathematical modelling of interactions between drugs or toxicants is essentially motivated by a search for efficiency. While assessing or predicting the toxicity of single substances is difficult enough, considering them by pairs, triplets, etc. is daunting. Even though much advances have been made on in vitro systems [START_REF] Weaver | Assessment of drug-drug interactions concepts and approaches[END_REF][START_REF] Pelkonen | Prediction of drug metabolism and interactions on the basis of in vitro investigations[END_REF][START_REF] Pelkonen | Inhibition and induction of human cytochrome P450 enzymes: current status[END_REF], the experimental route suffers the "curse of dimensionality":

Beyond pairs of interacting substances, factorial experimental designs of the order of size 2 n are required to assess simply the sign and order of interactions between n drugs or environmental stressors [START_REF] Bois | Optimal pollution control strategies in the presence of interacting toxicants[END_REF]. Modelling, with the aid of computers, has the potential of being much faster, provided that adequate models are available. In that case, "adequate" does not only mean "correctly fitting the data", but also "correctly predicting future observations", a much stricter requirement which early on oriented research toward mechanistic models. The domain of validity of a mechanistic model is conditioned by its structure and parametric assumptions, while the same domain for empirical (data-fitting) models is restricted to the range of the available data.

Therefore, this short review will leave aside empirical statistical models, such as population models [START_REF] Barrett | Application and impact of population pharmacokinetics in the assessment of antiretroviral pharmacotherapy[END_REF], or semi-empirical dose-response models [START_REF] Jonker | Towards a mechanism-based analysis of pharmacodynamic drug-drug interactions in vivo[END_REF], however interesting they are. We will not review either qualitative or semi-qualitative datamining approaches [START_REF] Bonnabry | Quantitative drug interactions prediction system (Q-DIPS) -A dynamic computer-based method to assist in the choice of clinically relevant in vivo studies[END_REF][START_REF] Isoherranen | Qualitative analysis of the role of metabolites in inhibitory drug-drug interactions: literature evaluation based on the metabolism and transport drug interaction database[END_REF], even though they can be a useful first screen. If we simply examine mechanistic models, two classes emerge: A first group of mechanistic models aims at predicting interactions at the enzyme or receptor level, through molecular modelling or QSAR approaches, or a combination of both [START_REF] Refsgaard | In silico prediction of cytochrome P450 inhibitors[END_REF][START_REF] Schuster | Predicting drug metabolism induction in silico[END_REF][START_REF] Vedani | The challenge of predicting drug toxicity in silico[END_REF]. A second group of models, which will be the subject of this review, integrates pharmacokinetic and enzyme or receptor pathways modelling. As this "physiological" or "systems biology" approach is the most comprehensive, it tends to couple with advances and tools or concepts from the above empirical and QSAR models [START_REF] Ekins | A combined approach to drug metabolism and toxicity assessment[END_REF]. To further reduce the scope of this paper without loosing generality, while still capturing some well-known sources of interactions among chemical substances, we will focus here on metabolic interactions. Those occur because a relatively limited number of enzyme species process a majority of xenobiotics and endogenous chemicals. For any of those reactions, the number of enzyme molecules is itself limited, leading to predictable bottlenecks in the detoxification or, conversely, in the activation of chemicals [START_REF] Pelkonen | Inhibition and induction of human cytochrome P450 enzymes: current status[END_REF][START_REF] Isoherranen | Qualitative analysis of the role of metabolites in inhibitory drug-drug interactions: literature evaluation based on the metabolism and transport drug interaction database[END_REF].

We mentioned above that experimental methods suffer from the curse of dimensionality, but to be fair modelling also faces the same problem. However, new solutions are emerging from pathways' modelling, which is proceeding at a fast pace (see http://systems-biology.org) [START_REF] Demin | Kinetic Modelling in Systems Biology[END_REF]. Although quite complex, pathways' modelling can actually reduce the dimensionality of the problem, as we will see in the following.

Physiologically based pharmacokinetic (PBPK) models

PBPK models are increasingly used in drug development and regulatory toxicology to predict the kinetics and metabolism of substances in the body [START_REF] Leahy | Progress in simulation modelling for pharmacokinetics[END_REF][START_REF] Parrott | Application of full physiological models for pharmaceutical drug candidate selection and extrapolation of pharmacokinetics to man[END_REF][START_REF] Barton | Characterizing uncertainty and variability in physiologically-based pharmacokinetic (PBPK) models: state of the science and needs for research and implementation[END_REF][START_REF] Bouvier D'yvoire | Physiologically-based kinetic modelling (PBK modelling): meeting the 3Rs agenda -The report and recommendations of ECVAM Workshop 63a[END_REF][START_REF] Loizou | Development of good modelling practice for physiologically based pharmacokinetic models for use in risk assessment: the first steps[END_REF][START_REF] Jamei | The Simcyp population-based ADME simulator[END_REF]. In these models, the body is represented by a set compartments corresponding specific organs or tissues, and the transfers or transports of drugs are dictated by various physiological flows (blood, bile, pulmonary ventilation, etc) [START_REF] Leahy | Progress in simulation modelling for pharmacokinetics[END_REF]. A system of differential equations can be written, with parameters representing blood flow rates, organ volumes etc., for which information is available in the published scientific literature or may be obtainable in vitro. Numerical integration of that differential system computes the quantity and concentration of the drug considered in each compartment, as a function of time and exposure dose. Indeed, such a description of the body is approximate, if not rough, but a balance has to be found between precision (and therefore complexity) and simplicity (ease of use).

The first report of a PBPK model able to describe simultaneously administration, kinetics and metabolic interactions was for warfarin and bromosulfophtalein [START_REF] Luecke | Drug elimination interactions: analysis using a mathematical model[END_REF]. The application to drug-drug interactions did not really take on at that time, in part because PBPK model parameterisation was difficult and did not adhere to the standards of statistical model fitting that was imposed to pharmacokinetic data analysis.

Interest in PBPK modelling then shifted to the regulatory toxicology arena, where data were scarce and predictive modelling (for low dose or inter-species extrapolations) much needed [START_REF] Menzel | Physiological pharmacokinetic modeling[END_REF]. That shift led to sporadic applications of PBPK models to interactions between pairs of substances [START_REF] Andersen | Quantitative evaluation of the metabolic interactions between trichloroethylene and 1,1dichloroethylene in vivo using gas uptake methods[END_REF][START_REF] Purcell | In vivo metabolic interactions of benzene and toluene[END_REF][START_REF] Tardif | Physiologically based modeling of the toxicokinetic interaction between toluene and m-xylene in the rat[END_REF][START_REF] Krishnan | Physiologically based pharmacokinetic analyses of simple mixtures[END_REF][START_REF] Tardif | Physiologically-based pharmacokinetic modeling of a mixture of toluene and xylene in humans[END_REF][START_REF] El-Masri | Physiologically based pharmacokinetic/pharmacodynamic modeling of chemical mixtures and possible applications in risk assessment[END_REF][START_REF] El-Masri | Exploration of an interaction threshold for the joint toxicity of trichloroethylene and 1,1-dichloroethylene: utilization of a PBPK model[END_REF][START_REF] Yu | Physiologically based pharmacokinetic modeling of metabolic interactions between n-hexane and toluene in humans[END_REF]. For reviews of the state of the art at that time, see [START_REF] Simmons | Application of physiologically based pharmacokinetic modelling to combination toxicology[END_REF][START_REF] Haddad | Physiological modeling of toxicokinetic interactions: implications for mixture risk assessment[END_REF]. The study of more complex mixtures, up to the five-component mixture of dichloromethane, benzene, toluene, ethylbenzene, and xylene was an interesting advance [START_REF] Haddad | Physiological modeling of the toxicokinetic interactions in a quaternary mixture of aromatic hydrocarbons[END_REF][START_REF] Haddad | A PBPK modeling-based approach to account for interactions in the health risk assessment of chemical mixtures[END_REF]. But, the exercise is quite painstaking and was not extended significantly afterward.

The beginning of this century witnessed a renewed interest in PBPK modelling in the context of drug development, including for drug-drug interactions [START_REF] Chien | Physiological approaches to the prediction of drug-drug interactions in study populations[END_REF][START_REF] Kato | Evaluation of methods for predicting drug-drug interactions by Monte Carlo simulation[END_REF][START_REF] Rostami-Hodjegan | In silico' simulations to assess the 'in vivo' consequences of 'in vitro' metabolic drug-drug interactions[END_REF], albeit without much awareness of the work being pursued in environmental toxicology [START_REF] Dobrev | Assessing interaction thresholds for trichloroethylene in combination with tetrachloroethylene and 1,1,1trichloroethane using gas uptake studies and PBPK modeling[END_REF][START_REF] Lee | Comparison of pharmacokinetic interactions and physiologically based pharmacokinetic modeling of PCB 153 and PCB 126 in nonpregnant mice, lactating mice, and suckling pups[END_REF][START_REF] Isaacs | Visualization-based analysis for a mixedinhibition binary PBPK model: determination of inhibition mechanism[END_REF][START_REF] Campbell | A PBPK modeling assessment of the competitive metabolic interactions of JP-8 vapor with two constituents, m-xylene and ethylbenzene[END_REF]. However, this is changing, but still, even recent publications focus primarily on binary drug-drug interactions [START_REF] Einolf | Comparison of different approaches to predict metabolic drug-drug interactions[END_REF][START_REF] Chenel | Drug-drug interaction predictions with PBPK models and optimal multiresponse sampling time designs: application to midazolam and a phase I compound. Part 1: comparison of uniresponse and multiresponse designs using PopDes[END_REF][START_REF] Buur | A physiologically based pharmacokinetic model linking plasma protein binding interactions with drug disposition[END_REF].

Systems biology pathways modelling

A problem faced by all the above investigators was to efficiently describe the cascades of reactions for the set of substances considered, in addition to the eventual transport and elimination of the metabolites. Assume that we describe each metabolic reaction by a set of four differential equations (for the formation of the enzyme-substrate complex and its dissociation into enzyme plus product), for 50 substances leading to 50 metabolites, all of them distributed primarily in 10 organs or tissues of the body. The corresponding model would consist in about 1200 equations and at least a thousand parameters. The difficulty or even the impossibility to manually develop and check such a model is obvious.

Dealing with a large number of biological reactions is a problem faced daily by systems biologists. To solve it, a set of tools has been developed, and in particular the Systems Biology Markup Language (SBML), a conventional syntax for representing biochemical pathway models [START_REF] Hucka | The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models[END_REF]. SBML can code for models of metabolism, signalling, transcription etc. An international community of software developers and users maintains SBML since the mid-2000 (cf. http://sbml.org). SBML defines a uniform symbolic coding system that enumerates reactants, products, parameters, reactions, their formulae etc., needed to define and simulate unambiguously the specified reactions. The model representation is symbolic and does not directly refer to a set of differential equations or any specific mathematical representation. That facilitates the interpretation of the model, and its automatic translation and integration into diverse software. SBML is de facto becoming a standard for model specification and exchange in systems biology [45].

Automated coupling of PBPK and systems biology modelling

To be able to deal with more complex mixtures, streamlining the process and taking advantage of the emerging tools of systems biology is a research avenue worth exploring. However, that exercise is still difficult when starting from a Henry-Michaelis-Menten description of metabolic reactions, because the number of parameters needed increases with the square of the number substances considered. In those conditions, modelling is not much more advantageous than experiments.

The challenge is to efficiently describe together cascades of reactions and the eventual transport and elimination of intermediate and final metabolites. The associated model can be very complex, even if limited to pharmacokinetic and metabolic interactions, because the transport of each metabolite requires a set of differential equation of its own. The main steps of a newly proposed "divide and conquer" approach [START_REF] Bois | Modeling physiological of metabolic interactions[END_REF] are:

-Develop for each substance and each of its metabolic reaction an SBML-coded model. It is necessary at this step to avoid simplifications of the Henry-Michaelis-Menten type, and to account for the balance of the various quantities of enzymes and enzyme-substrate complex involved, hence using the individual reaction micro-constants, rather than aggregated Henry-Michaelis-Menten constants.

-Define a general PBPK model template, applicable to all the substances considered (for example, the model shown on Figure 1). Such a model is needed to describe the relationship between personal exposure and internal concentrations at the site metabolism, which condition the intensity of interactions between substances; -Automatically build a global model of transport and metabolism of all the substances in the body, on the basis of the previously defined elements.

-Use that model to perform predictive simulations of substances' interactions.

The global model can be used to study the behaviour of complex mixtures, representative of the reality of human population exposures.

The above steps can be performed with the help of version 5.3.1 of the GNU MCSim software. (http://www.gnu.org/software/mcsim) [START_REF] Bois | MCSim: a simulation program[END_REF][START_REF] Bois | GNU MCSim: Bayesian statistical inference for SBML-coded systems biology models[END_REF]. GNU MCSim can use a PBPK model as a template for the description of the pharmacokinetics of each drug (and metabolites) involved in a reaction network. The template equations are copied for each chemical species and aggregated to form a large system of differential equations, able to describe the simultaneous pharmacokinetics and interactions of all the substances considered. Physiological parameters (for example the liver volume) keep the same value for all chemical species (see Table 1), but other parameter values may differ between chemical species. GNU MCSim recognises compartments defined in SBML models, maps them to the PBPK template compartments, and automatically locates the specified biochemical reactions in the correct organs. The substances placed (in SBML) outside of any physiological compartment are automatically distributed in all the compartments of the template PBPK model.

Analysis of a PBPK-coupled metabolic network

To illustrate the automated approach just described, and to explore the interaction properties of a sufficiently complex network, a set of primary metabolic reactions were simulated for a list of 50 hypothetical parent chemicals. For each of those, one or two metabolic reactions involving one of 10 possible enzymes were randomly generated, each with equal probability. The differential equations corresponding to a given substrate (S) and enzyme (E) pair are:

ES K t dP ES K K E S K t dES ES K K E S K t dE ES K E S K t dS 3 3 2 1 3 2 1 2 1 (1)
where ES represents the transient enzyme-substrate complex, P the reaction product, and the K i the reactions' micro-constants. The product species were randomly drawn among the 50 parents (thereby generating cross-linked reactions) plus another 150 other species (purely terminal products). In addition 5 parent species (S18 to S20) were supposed to induce an increased synthesis of enzymes 7 and 9, according to the following equations:

S K E K t dE 5 4 1 (2)
The random draws mentioned above were all from discrete uniform distributions. The enzymes and enzyme-substrate complexes were supposed to stay local to the liver. Parent species and metabolites were defined as circulating, and therefore distributed according to the generic PBPK model of Figure 1. The metabolic network generated by that procedure is presented graphically in Figure 2. The network is relatively complex and highly interconnected, like natural metabolic networks can be. That network was used as a basis for all the following simulations and results.

For simplicity, and without loss of generality, all substances were supposed to be non-volatile and excreted by the kidney, in addition to being eventually metabolised. Substances S1 to S20 (group A) were supposed to be natural A first step is to evaluate the "extent" of the possible interactions between substances of the network presented on Figure 2. In that example case, every substance affects the metabolism of any other one either directly or indirectly (there is a path connecting every two nodes of Figure 2). Yet, the previous statement is qualitative and solely based on the structure of the metabolic network. In reality, to paraphrase Paracelsus, the dose makes the interaction. If we turn to quantitative analyses, Figure 3 shows the simulated time-course of venous concentrations of all substances during the co-administration of groups A and B substances. These results were obtained by numerical integration of the global model for joint pharmacokinetics and metabolism of the substances. Parameters were set the values given in Tables 1 and2. The initial quantities of substances in the body were set to zero, except for the enzymes, set at an arbitrary initial liver concentration of 1 mM each. Ingestion rates per minute for group A and B substances were set at a high value (0,1 mM) for 4 hours and to zero afterward.

Since the parameter values for all substances were the same in that simulation, the noticeable differences observed on Figure 3 between substances are entirely due to the cascading network structure and to interaction effects. The global model is able to predict the kinetics of parent and daughter species in any compartment of the associated PBPK model for any exposure schedule. Such a simulation only takes a fraction of second on an i686processor.

To have an overview of the magnitude of the possible interactions, it is possible to compute an "interaction index" for each chemical species. For a given species, we define that index as the ratio of the concentration resulting from a joint exposure to groups A and B substances over the sum of the concentrations resulting from separate exposures to the same substances, at the same dose. In the absence of interactions the index value is 1, it is superior to 1 if the joint exposure leads to a disproportionate increase in the concentration of the substance considered, and lower than 1 otherwise. Figure 4 gives that index for all species, as a function of the ingestion rate of substances of groups B (the same rate was used for all substances of that group), 24 hours after the start of 4-hour exposures.

Substances of group A were administered at the rate of 0.1 mM/min. During exposure to the mixture, the internal concentrations of groups A and B substances were in general increased (we will call that "synergy"). On the contrary, metabolites (group C) had reduced concentrations on average (an effect we will term "antagonism"). Synergy (which would correspond to an increased toxicity if the substances were toxic) was observed in 60% of the cases. The concentration increase at 24 h, due to interactions, reached a factor 7 (substance S6).

Antagonism (by inhibition of metabolite formation) lowered concentrations by a factor 2 at most (for S39, metabolite of S20). At lower doses of group B substances, interaction effects dampened. Maximal synergies went down to a factor 2 at most (still for S6). Antagonisms reached at most a factor 1,5 (for S59 and S38), but interactions with group A substances (which are supposed to represent natural substrates, ingested at relatively high doses) were still observed, so that interactions indices did not return to 1. Note that for a few substances, the sign of the interaction changed with dose.

The above results were obtained with a model whose metabolic parameter values were all set to the arbitrary values given in Table 2. They are simply results of a stylised case study, aimed at exploring the possible extent of interactions and their behaviour in a realistically large metabolic network. In the real world, the parameters characterising each drug or substance might be very different and exposure levels and routes will be highly variable. To simulate such a variability, it is possible to perform Monte-Carlo simulations [START_REF] Bois | Precision and sensitivity analysis of pharmacokinetic models for cancer risk assessment: tetrachloroethylene in mice, rats and humans[END_REF], by randomly sampling parameter values. Besides, the metabolic network studied here is only a plausible example and our results depend on it. A specific real world network would probably respond differently, and this should be the subject of specific investigations.

Advantages of a micro-constants approach

The potential complexity of interaction effects between drugs, food constituents, and environmental chemicals, is an open problem since many decades. From a sheer mathematical point of view it is clear that combinatorial approaches (e.g., exploring all possible pairs, triplets etc.), be they experimental, statistical or mechanism-based, are impractical as soon as the complexity of the mixtures considered exceeds a small handful of substances. Modelling approaches do also suffer that "curse of dimension" unless they decompose the problem and take advantage of their power for synthesis. Our approach tends toward that goal and takes advantage of recent advances in bioinformatics and biomathematics at three levels: Automatic reconstruction of a metabolic network on the basis of individual reactions; Automatic coupling to an ADME prediction model; Fewer parameters through the use of reaction micro-constants (instead of the Henry-

Michaelis-Menten approximations).

To illustrate the latest point, crucial to the feasibility of the proposed approach, Figure 5 shows the number of model parameters required by various approaches as a function of the number of interacting substances considered. When using reaction micro-constants, the number of parameters required increases only linearly. The use Michaelis-Menten inhibition constants imposes at least a quadratic increase (when dealing only with pairs of chemicals, that is with order two interactions) and quickly becomes prohibitive.

The above advantages translate in feasibility for large networks of substances, ease of development and rapidity. The development of the template PBPK model of Figure 1 took about an hour and that of SBML models about three hours (using JDesigner version 2.1c). The compilation of all SBML models with the PBPK template took about a second. Using the global model and the GNU MCSim software, 10000 stochastic simulations can be run in 6 minutes on a i686 processor under GNU/Linux. The analysis and interpretation of the results take much longer.

However, to be fully operational, the systems biology approach presented here still requires gathering qualitative and quantitative data on the individual metabolism of the target substances. In that context, going back to experimentation becomes useful, with the advantage of being better targeted and therefore more efficient. In particular, the reaction micro-constant required (K 1 , K 2 and K 3 , in Table 2) are seldom available in the published literature, as everyone got used to Michaelian V max , K m , and K i constants (which suffice when considering simple binary interactions). However, at least for some enzymes, simple modifications of standard activity measurements allow an easy determination of the micro-constants in vitro (work in progress in our laboratory.) We are also exploring quantum chemistry predictions of those microconstants.

Perspectives

The ab initio modelling of metabolic interactions in cocktails of drugs or substances is feasible and even promising. Computing can be fast and the development of the model is limited to that of individual drug sub-models, reusable each time that that substance is to be included in a mixture. The latest developments of the systems biology approach presented here requires slightly different parameters than the Henry-Michaelis-Menten constants usually available, and its application to concrete cases will require a specific experimental effort, probably worth the investment on a medium to long term.

This system biology approach can obviously be extended to the description of pharmacodynamic interactions arising from signalling and regulation cascades.

However, while that approach gives access to the study of the behaviour of complex networks, the number of such networks and of realistic mixtures of drugs, food-borne chemicals and environmental contaminants is extremely large.

The future in this area, with the new tools available, may well lie in formulating general and robust statements about complex mixtures, with specific investigations of a few very relevant ones. Urinary excretion takes place in the kidney, part of viscera. Metabolism is not figured, and is specified in separate Systems Biology Markup Language models.

For symbols, see Table 1.

Muscle and skin

Viscera

Fat

Liver Urinary excretion takes place in the kidney, part of viscera. Metabolism is not figured, and is specified in separate Systems Biology Markup Language models.

f m F t f v F t f f F t f l F t V m V v V f V l V t
For symbols, see Table 1.

CellDesigner 4 .

 4 0.1 was used to generate Figure 2. The other graphs were plotted with Kaleidagraph v4.02 and R v2.8.1.

Figure 1 :

 1 Figure 1: Graphical representation of a physiologically based

Figure 2 :

 2 Figure 2: Computer-generated hepatic metabolic network for 50 hypothetical

Figure 3 :

 3 Figure 3: Simulated blood kinetics of the group A, B and C substances. The

Figure 4 :

 4 Figure 4: Interaction indices for the substances of groups A, B and C, as a

Figure 5 :

 5 Figure 5: Number of model parameters required by various approaches as a

Figure 1 :

 1 Figure 1: Graphical representation of a physiologically based

Figure 2 :Figure 4 :

 24 Figure 2: Computer-generated hepatic metabolic network for 50 hypothetical

Figure 5 :

 5 Figure 5: Number of model parameters required by various approaches

  chemicals, such as those present in food (phenols, etc.), absorbed daily at relatively high dose. Substances S21 to S50 (group B) were drugs absorbed at lower rate. The group C substances, with numbers between 51 and 200, correspond to terminal metabolites. There are only 46 of those metabolites because the random drawing of substances for possible reactions did not select all of the 150 possible choices.

	GNU MCSim version 5.3.1 (http://www.gnu.org/software/mcsim) was used to
	generate and simulate the global interaction model. That model had 466 state
	variables and differential equations for the quantities of parent species and
	metabolites in the various compartments of the PBPK model, and quantities of

The PBPK model template shown in Figure

1

was used to describe the kinetics of group A, B and C substances in the human body. It corresponds to four differential equations (the lung exchange kinetics are modelled with a simple algebraic equation under an instantaneous equilibrium hypothesis) and ten physiological parameters (see Table

1

). The template model does not include the metabolic reactions affecting the substances to which the body is exposed, those were automatically generated, using equations 1 and 2 above, for all reactions of the network. The substance-specific parameters of the PBPK template (partition coefficients, urinary elimination clearance) and the metabolic rate constants were set to the arbitrary values given in Table

2

.

enzymes and enzyme-substrate complexes in the liver, and 818 parameters.

Acknowledgements

I thank Karen Watanabe and Stefanie Stöckl for helpful discussions. This work was supported by the European Commission, 6 th FP Priority 6, project 2-FUN [contract #036976] and the French MEEDDAD (PRG189_08_DRC03).

Tables

 Table 1:Physiological parameter values for the PBPK model shown in Figure 1.

Parameter (a) Symbol Value Total blood flow F t 6.4 Fractional blood flows Fat f f 0.09 Liver f l 0.24 Muscle and skin

Total body volume

(a) Units: volumes are in litres, total blood flow in litres/min. (b) The fraction of blood flow to the viscera is computed by difference between 1 and the sum of the fractional flows to the other compartments.

(c)

Those volumes correspond to those of a 30 years old woman [START_REF] Van Der Molen | A generic toxicokinetic model for persistent lipophilic compounds in humans: an application to TCDD[END_REF]. (d) The volume of the muscle and skin compartment is computed by difference between the total body volume, 10% of it (corresponding to bones) and the other volumes.


Alexandre R R Pery
email: alexandre.pery@ineris.fr

Adrien Troise

Sylvie Tissot

Jean-Marc M Vincent

Alexandre R R Péry

Comparison of models to analyze mortality data and derive

Keywords: accidental risk, modeling, lethality data, parameters estimation, probit

Introduction

To prevent risks related to accidental releases of dangerous chemical substances in the atmosphere, risk managers need acute toxicity thresholds in association with accident scenarios to produce safety reports and design emergency plans. In France, they usually determine the zones of lethal, irreversible and reversible effects relative to the location of plants storing, producing or using toxic substances, especially for land use planning. The lethality threshold is related to a certain percentage of death occurring during the experimental test or in the following fourteen days post exposure (including animal sacrificed for ethical reasons). The "irreversible effects" correspond to the persistence over time of a lesion or a functional damage induced by an exposure. Three types of irreversible effects are pointed out, lesion without functional repercussions, lesion with functional repercussions (like bronchopneumopathy, pulmonary fibrosis, necrosis of olfactory epithelium with anosmia) and the functional irreversible impairment (like asthma). The "reversible effects" correspond to a return to the level of health prior to exposure (immediately or in a reasonable time after).

It is therefore crucial to evaluate properly the thresholds, because they determine the distances of effects. Indeed, if the thresholds are overestimated, distances are overprotective with economic impact. In contrast, if the thresholds are underestimated, the health of the exposed population is threatened.

A French methodology was developed to evaluate the quality of the data available and to deduce from these data acute toxicity thresholds. It comprises many steps. The first step is the selection of an experimental key study (mainly, an animal study) for each effect, based on the method developed by [START_REF] Klimisch | A systematic approach for evaluating the quality of experimental toxicological and ecotoxicological data[END_REF] for quality assessment and on expert judgment for relevance of the observed effects related to the type of threshold. The second step consists in finding the relevant critical toxic effect for the two types of effects which are considered in addition to lethal effects, i.e. irreversible and reversible effects. The third step is the data analysis based on a statistical model. The fourth step considers the extrapolation from animals to humans (with or without uncertainty factors). Here we focus on lethality data only.

Usually, the data to analyze are rodent lethality measurements for different exposure concentrations and different exposure durations, observed after at least a fourteen-day period.

It can however happen that the information is available only for single exposure duration.

Standard dose-response models are generally based on the Haber's law, or its generalizations (ten [START_REF] Ten Berge | Concentration--time mortality response relationship of irritant and systemically acting vapours and gases[END_REF]. These generalizations state that the effect for exposure to concentration C during a period of time t is a function of (C n t), named fixed effect level [START_REF] Jarabek | Considerations of temporal toxicity challenges current default assumptions[END_REF], where n is called the Haber constant. Thus, in the standard probit model, the probability of death equals:

) ) ) log( ) log( ( ) (  µ t C n F death P   
where F is the cumulative distribution function of a normal distribution with mean µ and standard deviation σ [START_REF] Diack | Pharmacokinetic-pharmacodynamic models for categorical toxicity data[END_REF]. In this formulation, each individual has a given threshold, log-normally distributed. If he is exposed to an effect level exceeding this threshold, death will occur.

There are many limits to this model. First, it is not possible to analyze data for which only one exposure time is available, because only two combinations of parameters (n/σ and (log(t)µ)/σ) could be identified in this case. This model has also been criticized for its lack of biological parameters by [START_REF] Diack | Pharmacokinetic-pharmacodynamic models for categorical toxicity data[END_REF] who have proposed an alternative model (the so-called PKPD model). However, the model they propose has even more limitations in a risk assessment perspective.

First, the gain in terms of realism compared to the standard probit approach is far from obvious. Indeed, the proposed equation for kinetics is the following one:

keQ kaC dt dQ n  
with Q the internal quantity of substance in the tissue, C the exposure concentration, k a the pulmonary ventilation rate, k e the elimination rate, and n an unknown parameter analogous to the Haber constant. This is a very unusual kinetics equation, because there is no reason for the parameter n to be different from 1, i.e. for the intake rate of the substance not to be proportional to its concentration in the air. Authors try to justify their equation with some biological consideration in their discussion, but their reasoning would be acceptable only for very rapid kinetics.

Second, there are five parameters to estimate (ka, ke, n, but also µ and σ as in the probit model), compared to the three ones in the probit approach. In practice, when analyzing data for accidental risk assessment, with both the PKPD model and the standard probit model, we found that the common parameters have similar estimates but that the remaining parameters have large confidence intervals. This observation indicates an excessive number of parameters to estimate. In particular, the estimation of the kinetics parameters is not feasible when only data for one time of exposure are available.

Here, we propose to evaluate three alternative models. These models have the same numbers of parameters (three) as the standard probit model. The first two models incorporate chemical kinetics through a one-compartment model to add realism relative to compound uptake and elimination. They differ in the toxicodynamics part. The first one, which is called DEBtox (Kooijman and [START_REF] Bedaux | Statistical analysis of bioassays based on hazard modelling[END_REF], is based on a threshold approach which has already been used in ecotoxicology. The second one is based on the Haber's law and will be here referred to as Haber TKTD. The third approach we propose is an approximation of the probit standard by using a loguniform statistical distribution for F instead of a normal one. The comparison of the models we propose is based on the analysis of datasets generated from simulations with each of the models and on actual data which have already been used to derive toxicity threshold for accidental risk assessment.

Materials and Methods

Mathematicals models to analyse survival data

Standard probit

As already presented in the introduction, in the standard probit model, the probability of death equals:

) ) log( ) log( ( ) (      t C n F death P
where F is the cumulative distribution function of a normal distribution with mean µ and standard deviation σ. It consequently assumes that the individual threshold for response follows a lognormal distribution. When the product C n t exceeds its threshold, the individual dies.

Loguniform model

In this model, the probability of death equals:

) ) log( ) log( ( ) ( a b a t C n F death P    
where F is the cumulative distribution function of an uniform distribution with lower and upper bounds 0 and 1, and where a and b are bounds of the response. As for the standard probit, when the product C n t exceeds the threshold of an individual, the individual dies. There is no individual having a threshold below a and above b, which means that no individual are expected to die for C n t values below a and no individual are expected to survive for C n t values over b. The Loguniform model can be seen as a linearization of the standard probit model. The comparison will permit to assess the relevance of the choice of the statistical distribution for F.

DEBtox model

DEBtox model is a mathematical model that has first been developed to analyse aquatic ecotoxicological survival data. It has been proposed by Kooijman and Bedaux (1996). It is particularly adapted for the analysis of toxicological data obtained under time varying exposures (Pery et al., 2001;Pery et al., 2002).

In DEBtox model, a kinetics module, accounting for the dynamics of compound body concentration is coupled with an effects module. To keep the same number of parameters as for the standard probit, we assumed, as in the standard DEBtox model, that all individuals have a common threshold for effects. Once this threshold is exceeded by the dose at target organ level, the probability to die is not 100% but increases linearly as a function of this dose. Dose at target level is described using the following linear one-compartment kinetics model:

)) ( ) ( ( ) ( t c t C k t dt dc i e i  
where k e is the elimination rate, C is the concentration measured in exposure air, and c i the scaled body concentration. This parameter corresponds to the ratio of the amount of compound in the whole body to the body volume. It is scaled by the bioconcentration factor (a constant corresponding to the ratio of the concentration in the target organ to the concentration in air at toxicokinetics steady state) in order to ensure the feasibility of parameters estimation. Toxic effects occur only when c i (t) exceeds a threshold, the NEC (No Effect Concentration), which corresponds to the maximal toxicant concentration at target organ level that can be handled by regulation systems without generating detectable effects on mortality. Survival probability in exposed organisms is described based on c i (t), which drives toxicodynamics. Death being assumed to be a stochastic process, the probability q(t) to survive until time t is defined as:

   t d h t q 0 ) ) ( exp( ) (  
where h(τ) is the hazard rate at time τ. For a small time interval dτ, h(τ) dτ represents the probability of dying between τ and τ+dτ, for an organism which has survived until time τ.

Toxicant hazard linearly increases with the difference between c i (t) and the threshold value when this threshold is exceeded:

        0 ) ( ), ) ( ( ) ) ( ( ) ( ), ) ( ( t h NEC t c if NEC t c b t h NEC t c if i i i
where b is the killing rate, a descriptor of the intensity of effects. The effects data permit both to estimate the parameters b and NEC, but also the kinetics parameter. For instance, under constant exposure, it is related to the time at which the level of effects reaches a constant value.

Haber-TKTD model

This model has the same kinetics module as the DEBtox model. It differs from the latter model in the effects module, with the following definition for h:

n t bC t h ) ( ) ( 
where b is a proportionality constant, and n is the Haber constant.

It is noteworthy to realize that, in case of n=1, we have, for high values of t, corresponding for instance to observation times (typically 14 days, whereas exposure is limited to 8 hours maximum): (-bCt) whatever the value for k e. This means that, in this case, effects are related to the Haber's law, regardless of the actual compound kinetics.

q(t)=exp

Data analyzed in this study

Simulated data

To test the ability of the models for data analysis, we simulated datasets according to each model and analyzed these data with all the models. We aimed at producing realistic data, with small effects for low exposure times and low exposure concentrations and high effects for high exposure times and high exposure concentrations. We used four concentrations (50, 100, 200 and 300 ppm) and four times of exposure (30, 60, 120 and 240 minutes). We generated twenty different data sets per model. For each condition in each dataset, we simulated the results for 10 individuals. First, we calculated the mean probability for an individual to be dead at the observation time. Using a random generator (from the commercial software Excel), we simulated mortality data by comparing, for each individual, the number generated (between 0 and 1) and the probability to be dead. If the random number was below this probably, the outcome for this individual was death.

For better comparison of performance of the different models, we chose parameters so that the theoretical value for the 5% letal concentration (5%LC 30 min , which is the concentration for which the mean expectation is 5% mortality after 30 minutes of exposure) is 100 ppm for all the simulations. The 5%LC 30 min is of great interest in accidental risk assessment. We then evaluated if the models were able to fit correctly the data, to estimate the parameters, and to estimate properly the 5%LC 30 min .

Real data

We analyzed data obtained for rats exposed to chemicals, with many times of exposure available to be able to use all our models, including those with kinetics. We had data for arsine (3 times of exposure, 5 to 7 concentrations, 10 female rats per condition) (IRDC, 1985a;IRDC, 1985b;IRDC, 1985c), cyanhydric acid (1 to 5 times of exposure, 4 to 39 concentrations, 10 rats per condition) [START_REF] Ballantyne | Acute inhalation toxicity of hydrogen cyanide vapor to the rat and rabbit[END_REF][START_REF] Blagden | Hydrogen cyanide: multiple exposure time acute inhalation toxicity study in the rat[END_REF][START_REF] Higgins | Acute toxicity of brief exposures to HF, HCl, NO2 and HCN with and without CO[END_REF], sulfuric acid (4 times of exposure, 3 or 4 concentrations, 10 to 14 mice per condition) [START_REF] Runckle | The toxicity of H2SO4 aerosols to CD-1 mice and Fisher-344 rats[END_REF], hydrogen sulfide (1 to 4 times of exposure, 3 to 7 concentrations, 10 mice per condition) [START_REF] Vernot | Toxic hazards research unit annual technical report[END_REF][START_REF] Zwart | Determination of concentration time mortality relationships to replace LC50 values[END_REF], nitrogen dioxide (10 times of exposure, 7 concentrations, 5 to 14 mice per condition) [START_REF] Hine | Pulmonary changes in animals exposed to nitrogen dioxide, effects of acute exposures[END_REF], vinyl chloride (2 to 4 times of exposure, 3 to 14 concentrations, 5 to 90 mice per condition) [START_REF] Mastromatteo | Acute inhalation toxicity of vinyl chloride to laboratory animals[END_REF][START_REF] Prodan | Experimental acute toxicity of vinyl chloride (monochloroethene)[END_REF] and sulphur dioxide (3 to 4 times of exposure, 3 concentrations, 14 to 70 mice per condition) [START_REF] Bitron | Delayed mortality of mice following inhalation of acute doses of CH2O, SO2Cl2, and Br2[END_REF]. These data have already been used, with the standard probit approach, to derive thresholds used in accidental risk management.

Estimation of parameters

Parameters were estimated by maximizing the logarithm of the likelihood of the data,

)) 1 ln( ) ( ) ln( ( 1 j j j k j j j p n m p n L      
where p j is the mean probability to be dead at the time of observation (14 days) for concentration c j , calculated with the model selected for data analysis, k the number of tested concentrations, n j the number of organisms dead at the time of observation for concentration c j , m j the number of exposed organisms for concentration c j , .

To perform that maximization, we coded a program in C++, which calculates the likelihood for combinations of large ranges of parameters, and selects the parameters providing the highest value. It is not time-saving but very robust, so that we avoid the selection of local maximums.

Results

Simulated data

Table 1 presents the results of the analysis of the simulated data. Except in two cases (data simulated with Loguniform model n=1 and DEBtox with rapid kinetics), regardless of the model we used to simulate data, the standard probit model provides the best fit of the data.

Moreover, the estimate of n is never significantly different from the theoretical value (Wilcoxon test, p<0.05). Regarding the estimate of 5%LC 30 min , it is significantly different from the theoretical value only for one simulation (DEBtox and slow kinetics).

DEBtox and Haber-TKTD models provide correct fits and parameters estimation only when they have generated the data they analyze. It is important to note that, for data generated by DEBtox in case of rapid kinetics, it outperforms all other models. For the other simulations, fits are of very bad quality (with more than twice the sum of the squares of the deviations (SSD) of standard probit model).

Despite a relative good fit of the data (low SSD values), the Loguniform model provides biased estimates (significant differences between theoretical value and estimated values, for the 5%LC 30 min in 6 of the 8 simulations), even when data have been generated with this model (case n=3).

Other simulations with different σ and a values (data not presented here) did not modify our conclusions.

Real data

Table 2 and 3 present the results of the analysis of the real data. To permit to compare the performance of fits between experiments, the SSD has been normalized by the number of animals used in the study. Except for cyanhydric acid, the fit is better using the standard probit model than the other ones. In contrast, Haber-TKTD model always provide the worst fits. The second best model is DEBtox in 5 cases over 7 ones and the Loguniform model in 2 cases over 7 ones. These results are in accordance with the analysis of simulated data. The value of n estimated with the Loguniform model is close to the one estimated with the standard probit model except for sulfuric acid. From this point of view, the Loguniform model performs better than our expectations based on the study with simulated data.

Possible use of the models in case of one exposure time only available

As we showed in the introduction, it is not possible to estimate the parameters of the standard probit model for data characterized by only one time of exposure. This is also the case for the Loguniform model. Even setting the lower bound at a fixed value (corresponding, for instance, to an available NOAEL) would not permit to estimate the two remaining parameters.

In contrast, for models with kinetics module, like DEBtox or Haber-TKTD, the parameters could be estimated in case of only one time of exposure, provided there are several observation times. We performed two additional datasets simulations, for each model, with 6 exposure concentrations (50, 100, 150, 200, 250 and 300 ppm), only one time of exposure (60 min) but two observation times (60 min and 15 days). We used the same parameters as for the previous simulations in case of slow kinetics and produced 20 datasets per model. The estimation of the parameters of the Haber-TKTD model was not satisfactory with a huge variability. In contrast, DEBtox parameters were adequately estimated. The NEC estimate was 25 ± 7 ppm, and the 5%LC 30 min estimate was 99 ± 8 ppm (the theoretical value is 100 ppm), which is a better estimate than the one obtained for only one time of observation but several times of exposure (Table 1).

Discussion

When the data for many times of exposure are available, the standard probit model performed at best relative to goodness of fit. The description was the closest for real and simulated data, even, in the latter case, when the data have been simulated with other models. In contrast, the Haber-TKTD model was the worst of all the models we tested. DEBtox model was also not able to provide satisfactory fits in most simulated cases but was often the second best model when analyzing real data. The Loguniform model provided fits with a quality close to the one of the standard probit model, but it suffered from biases in the estimation.

For the analysis of the data usually used in an accidental risk assessment context, i.e. data with only one observation time but several times of exposure, we can conclude that better fits are obtained when accounting for variability of response between individuals than for kinetics during exposure and post-exposure period. This is coherent with our remark that the mathematical description of effects with the Haber-TKTD model does not depend on k e when n equals 1 (See 2.1.4), as a consequence of measuring the data only when there is no more compound in the body. Models like DEBtox perform well when there are many measurements for different times of exposure, even if exposure concentrations are time-varying (Pery et al., 2001;Pery et al., 2002). For rodent lethality, measuring effects just after exposure and much later may be enough to improve the estimation of the kinetics parameter. Unfortunately, such data are seldom available. In the data we have, only those relative to Sulfur dioxide comprise several times of observation. We analyzed these data with DEBtox model. Sulfur dioxide was the chemical for which DEBtox performed at worst with data characterized by only one time of observation. However, with several observations throughout time, the sum of the squares of the deviations between mortality predictions and actual data normalized by dividing by the total number of exposed organisms in the study was divided by 1.8 compared to Table 2. Moreover, the estimate for the 5%LC 30 min was 986 ppm, which is very similar to the estimate value using the standard probit model (1025 ppm).

The standard probit model is not able to analyze data for which only one time of exposure is available. In this case, the methodology used in France recommends the use of the probit standard to derive dose response for the only time of exposure available and the extrapolation of the results using Haber's law with n=1 for longer duration than the one available, and n=3 for shorter duration [START_REF] Ten Berge | Concentration--time mortality response relationship of irritant and systemically acting vapours and gases[END_REF].

Among the models we studied, only DEBtox model would be able to analyze data characterized by only one time of exposure, provided mortality is reported at least two times, typically just after and long after exposure. However, it will be necessary to perform comparisons between DEBtox and standard probit models with real data and at least two observation times and several exposure times. This will permit to assess if the thresholds needed for accidental risk assessment are estimated with the same accuracy when applying DEBtox to data with several times of observation or applying the standard probit model to data with several times of exposure. Unfortunately, we could not find such data.

To conclude, it is recommended to use the standard probit model to analyze lethality data for accidental risk assessment, when data are available for many times of exposure. In contrast, there is no model that could replace the extrapolation with n=1 and n=3 to analyze data for only one time of exposure and extrapolate to other times if only one time of observation is available.

Table 1 Mean and standard deviations for the estimates of 5%LC 30 min (in ppm) and sum over all the simulated tests of the squares of the deviation (SSD) between mortality predictions and actual data. The different columns account for the different models used for the simulations, and the different lines for the different models used to analyze the simulated data. For the models in which the Haber constant is involved, the estimates for n are presented.

Standard Probit n=3

Standard

For

  each model, two sets of parameters values were used. The parameters values we selected for the simulations were n = 3, µ = 7.43, σ = 1, and n = 1, µ = 3.43, σ = 1, for the simulations based on standard probit model. This permits to cover the two options generally applied when fixing arbitrarily a value for n to extrapolate between doses (ten Berge et al., 1986). With the same rationale, we chose n = 3, a = 5.38, b = 7.38, and n = 1, a = 1.38, b = 3.38 for the Loguniform model. For DEBtox model, we chose k e = 0.01 min -1 , b = 0.0003 min -1 ppm -1 , NEC = 25 ppm, and k e = 10 min -1 , b = 0.0000228 min -1 ppm -1 , NEC = 25 ppm, to have a simulation with slow kinetics and one with rapid kinetics. With the same rationale, we chose k e = 0.01 min -1 , b = 0.0000022 min -1 ppm -1 , n = 2, and k e = 10 min -1 , b = 0.000000304 min - 1 ppm -1 , n = 2, for the Haber-TKTD model.

Table 3

 3 Estimates of n using standard probit or Loguniform models, with for the latter either all the data available or only those for 60 minutes of exposure.

	Probit n=1	Loguniform n=3	Loguniform n=1

Acknowledgements

Authors would like to thank Enrico Mombelli, Céline Brochot and two anonymous reviewers for their comments on the study. This work was supported by the French Ministry in charge of Ecology and Sustainable Development for this study, within the framework of Programme 181.

Conflicts of interest

The authors declare that there are no conflicts of interest.