Stochasticity in physiologically-based kinetics models: Implications for cancer risk assessment 2

Keywords: PBPK, butadiene, systems biology, stochasticity, cancer

In case of low dose exposure to a substance, its concentration in cells is likely to be stochastic.

Assessing the consequences of this stochasticity in toxicological risk assessment requires the coupling of macroscopic dynamics models describing whole body kinetics with microscopic tools designed to simulate stochasticity. In this paper, we propose an approach to approximate stochastic cell concentration of butadiene in the cells of diverse organs. We adapted the dynamics equations of a physiologically-based pharmacokinetic (PBPK) model and used a stochastic simulator for the system of equations we derived. We then coupled kinetics simulations with a deterministic hockey stick model of carcinogenicity. Stochasticity induced substantial modifications relative to dose-response curve, compared to the deterministic situation. In particular, there was non-linearity in the response and stochastic apparent threshold was lower than the deterministic one. The approach we developed could easily be extended to other biological studies to assess the influence at macroscopic scale of stochasticity for compounds dynamics at cell level.

INTRODUCTION

The dynamics of chemical species at the level of the cell are rather discrete and stochastic than continuous and deterministic, in so far as they are determined by the action of only a few molecules [START_REF] Rao | Stochastic chemical kinetics and the quasi-steady-state assumption: application to the Gillespie algorithm[END_REF]. Stochastic resonance in biology is a well known phenomenon, able to enhance detection and improve the transmission efficiency of weak information in nonlinear systems [START_REF] Hänggi | Stochastic resonance in biology. How noise can enhance detection of weak signals and help improve biological information processing[END_REF]. In particular, it explains why deterministic threshold can be exceeded, with a rate in relation to the intensity of the signal, even if the mean value of this signal throughout time is below the threshold. The concepts apply to toxicology: For low exposures to a chemical, stochasticity is likely to play a role in the occurrence of toxic effects at the cell level, even if the mean cellular concentration is below some deterministic threshold, inferred for instance, through mechanistic considerations [START_REF] Lovell | Dose response and threshold mediated mechanisms in mutagenesis: statistical models and study design[END_REF]. Studying stochasticity in toxicant cellular concentration is particularly relevant for carcinogenesis which might be induced at very low doses of exposure. In this paper, we present a modelling framework to simulate stochastic concentrations at cell level and to derive consequences for cancer risk assessment. Simulations are performed based on actual exposure levels and a physiologically-based toxicokinetic model for butadiene.

1,3-Butadiene is a highly volatile four-carbon chemical mostly made from the processing of petroleum. It can be detected in urban air pollution, cigarette smoke and gasoline vapors.

Butadiene is an established carcinogen, as a DNA-reactive chemical leading to production of DNA adducts in rodents in liver, lung and tissue [START_REF] Preston | Cancer risk assessment for 1,3-butadiene: data integration opportunities[END_REF]. It is hypothesized that butadiene carcinogenicity is a consequence of the genotoxicity of its metabolites [START_REF] Albertini | The use of non-tumor data in cancer risk assessment: reflections on butadiene, vinyl chloride, and benzene[END_REF]. It is not clear whether or not there is a threshold for effects in humans [START_REF] Preston | Cancer risk assessment for 1,3-butadiene: data integration opportunities[END_REF].

Physiologically based toxicokinetic (PBPK) models propose a realistic even if simplified description of the mechanisms of absorption, distribution, metabolism and elimination of chemicals in the body. In these models, the body is subdivided into various compartments representing specific organs or homogeneous groups of tissues linked and irrigated by blood vessels. Compartments are characterized by a set of parameters of physiological relevance (e.g., volume or blood perfusion rate) which play a crucial role in explaining the behavior of chemical substances in the body, and represent invariants across substances. A threecompartment physiologically based pharmacokinetics (PBPK) model has been proposed to describe the distribution of butadiene and the production of its first metabolite following oxidation, 1,2-epoxy-3-butene [START_REF] Brochot | Use of a chemical probe to increase safety for human volunteers in toxicokinetics studies[END_REF]. We use here an extension of this model with 23 compartments. The model parameters (organ volumes etc.) correspond to those of an adult man.

These macroscopic models can describe the distribution of chemicals in the different organs but are unable to capture the stochasticity at cell level. In contrast, a recent discipline, systems biology, aims at studying the dynamics of the components of a cell, and tools to study the influence of stochasticity at cell level have being developed. The main objective of this paper is to develop a methodology for coupling microscopic and macroscopic dynamics models to assess the consequences at organ level of stochastic chemical concentration and effects at cell level. Coupling PBPK models and systems biology is a step forward to develop integrated approaches able to relate information obtained at cell level, like for instance "omics" data, and effects on health.

MATERIALS AND METHODS

PBPK model for butadiene

Our PBPK model contains 23 compartments (See Figure 1). Compartments are mainly connected by blood circulation, by air exchange at the lung level, excretion to urine and feces, and metabolism. Concentrations C, in µg/L, are obtained at any time by dividing the quantity of butadiene, Q (in µg) by the compartment volume (supposed constant in time). Volumes are in L, time in min, flows and rates in L/min. Q and C depend on time but we omit the time argument when possible for simpler notation. A list of parameters and their values can be found in Tables I-IV. We assume that butadiene is only eliminated through metabolisation or exhalation, and that intake only occurs through inhalation.

For adipose tissue, adrenals, bone marrow, brain, breast, heart, kidneys, muscles, other organs and tissues, pancreas, skin, spleen, testes and thyroid, the differential equation giving the rate of change for the quantity of BD is:

        - × = ∂ ∂ i i art i i P C C F t Q (1)
where F i are blood flows (values given in Table 2), and P i are tissue over blood partition coefficients. These partition coefficients were calculated by multiplying the partition coefficients for fat (Table III) by fat content of each organ (Table IV). Indeed, for highly lipohilic organic chemicals, P i values can be approximated by the ratio of lipids in adipose tissues and blood [START_REF] Haddad | Relative lipid content as sole mechanistic determinant of the adipose tissue:blood partition coefficients for highly lipophilic organ chemicals[END_REF].

( )

        - - - = ∂ ∂ air over lung lung urt alv urt air pul urt P C C F C C F t Q _ _ (2)
where F pul designates the pulmonary ventilation rate (9 L/min for an average human) and F alv the alveolar ventilation rate (6 L/min) (ICRP, 2002).

For the lung:

        - +         - = ∂ ∂ lung lung ven total air over lung lung urt alv lung P C C F P C C F t Q _ _ (3)
where F total designates the sum of the F values for the following organs or tissues: adipose tissue, adrenals, bone marrow, brain, breast, heart, kidneys, liver, muscles, other organs and tissues, skin, testes and thyroid.

For the liver, blood comes from the arterial pool, spleen, pancreas, stomach and gut, and there is metabolism occuring in the liver (see further):

F P C F C F t Q × - × - × + × + × + × + × = ∂ ∂ (4)
where F liver , is:

stomach gut pancreas spleen eport liver F F F F F F + + + + = (5)
For arterial blood:

        - × = ∂ ∂ art lung lung total art C P C F t Q (6)
For venous blood:

ven total i i i i ven C F P C F t Q × -       × = ∂ ∂ ∑ (7)
where i designates the adipose tissue, adrenals, bone marrow, brain, breast, heart, kidneys, liver, muscles, other organs and tissues, skin, testes and thyroid.

For epoxy metabolite (EB) in the liver:

EB liver met liver liver met met P C F Q K t Q _ - × = ∂ ∂ (8)
For simplicity, the possibility of EB molecules to come back to the liver through blood circulation has been neglected, assuming susbtantial elimination through secondary metabolism and exhalation. To obtain more relevance, incorporating EB kinetics may be proposed later, based on the PBPK model for butadiene and major metabolites proposed by [START_REF] Brochot | Development of a physiologically based toxicokinetic model for butadiene and four major metabolites in humans: global sensitivity analysis for experimental design issues[END_REF].

Exposure and epidemiological data

We used data from [START_REF] Higashino | Exposure and risk assessment of 1,3-butadiene in Japan[END_REF] on human exposures to 1,3-butadiene in Japan.

The average concentration in the general environment is 0.25 µg/m 3 , with a background concentration 0.06 µg/m 3 in unpolluted areas. Exposure concentrations above 0.8 µg/m 3 are only found in vicinity of industrial activities. Lifetime excess cancer risk level is estimated at 10 -5 for an exposure concentration of 1.7 µg/m 3 [START_REF] Higashino | Exposure and risk assessment of 1,3-butadiene in Japan[END_REF]. In Japan, 0.03% of the total population (that is 36 000 persons) are exposed to concentrations exceeding this value. With a molar mass of 54.09 g/mol for butadiene, 1.7 µg/m 3 corresponds to 18.7 10 12 molecules/L.

Simulation Software used

We used the Systems Biology Workbench (SBW) [START_REF] Sauro | Next Generation Simulation Tools: The Systems Biology Workbench and BioSPICE Integration[END_REF], version 2.0.39. We implemented our PBPK model in Jdesigner. The model was then converted to an SBML file and input in the Dizzy 1.11.4 software [START_REF] Ramsey | Dizzy:stochastic simulation of large-scale genetic regulatory networks[END_REF], which is able to perform stochastic simulations of chemical kinetics. We used the Gillespie stochastic algorithm [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF] for stochastic simulations. It is an algorithm for modeling the kinetics of a set of coupled chemical reactions, taking into account stochastic effects from low copy numbers of the chemical species. With the same differential equation, the deterministic approach regards the time evolution as continuous, whereas the stochastic approach regards the time evolution as a kind of random-walk process. In Gillespie's approach, chemical reaction kinetics are modelled as a Markov process in which reactions occur at specific times separated by Poisson-distributed intervals. The mean interval is recomputed each reaction time. At each reaction time, a specific chemical reaction occurs, randomly selected from the set of all possible reactions with a probability given by the individual reaction rates.

Stochastic simulations at steady state

Exact stochastic simulators track the actual number of molecules involved in a set of reaction in a given portion of space. It is impossible for them to handle as many molecules as can be found in an entire organ. We first computed the steady state concentration values in each organ using the deterministic (ordinary differential equation) simulator of the SBW, with continuous inhalation exposure concentrations of 0.25 and 1.7 mg/m 3 of BD. Initial condition for the number of molecules in a cell was set to the deterministic value from the PBPK model.

We consider that a human has typically about 10 14 cells. In our PBPK model, the total volume of a man is 75 L, which results in a mean density of 1.37 10 12 cells/L.

Except for liver, the uptake of butadiene is performed from the arterial blood, in which butadiene concentration was set at steady state deterministic value from PBPK model. The kinetic equation for one cell is similar to that of the organ to which it belongs except for the flux value:

        - = organ cell art cell cell P C C F dt Q d ) ( (9) with organ organ cell V density Cell F F * _ =
This equation cannot be computed directly in the softwares we use because of rounding errors due to the extremely small value of F cell (less than 10 -13 for most organs). We reformulate equation ( 9):

        - × × ℵ = organ cell butadiene art organ organ cell P N M density Cell C V F dt N d _ ) ( ( 10
)
where ℵ is the Avogadro's number, M butadiene is the molar weight of butadiene and N cell the number of molecules in the cell. This equation can be studied with Jdesigner and Dizzy with a fixed C art value. Equation (4) was reformulated in a similar way.

We performed simulations with five organs: spleen, which has the lowest partition coefficient over blood (0.77), fat with the highest one ( 22), marrow, which has a low scaled flux value (0.0786 min -1 ), kidney with a high scaled flux value (3.9838 min -1 ) and liver, for which we studied the concentrations of both butadiene and its first metabolite.

That choice of organs is also relevant relative to cancer risk of butadiene. Mutagenicity has been shown to occur for spleen and marrow of rodents [START_REF] Preston | Cancer risk assessment for 1,3-butadiene: data integration opportunities[END_REF]. The same author reported butadiene-induced lymphomas so as lung and liver tumours in mice.

We ran 100 simulations at steady state. For spleen, kidney and liver, we reported the number of molecules at time 100 min, which is large compared to the time needed to reach equilibrium in the deterministic model for these organs. For marrow and fat, we reported this number at 500 min and 10000 min respectively due to a longer time to reach equilibrium.

Coupling with cancer PD model

To assess how dose-response curves are affected by stochasticity at the cellular level, we coupled the stochastic PBPK model described above to a linear model with a threshold for effect, which belongs to the family of the "hockey stick models" introduced into carcinogenic risk assessment by [START_REF] Cornfield | Carcinogenic risk assessment[END_REF]. For a given cell at a given time, the probability R of carcinogenesis is supposed to be proportional to max(N(t)-N 0 , 0) where N(t) is the number of molecules of BD metabolites per cell at time t, N 0 the threshold number of molecules to get an effect. We chose 1 as the threshold number of molecules per cell able to initiate liver carcinogenesis with a non-zero probability. This is the minimum reasonable number and it corresponds to approximately twice the general environment exposure in Japan [START_REF] Higashino | Exposure and risk assessment of 1,3-butadiene in Japan[END_REF]. This is just a choice for simulations. It is not based on particular knowledge about butadiene mechanisms of carcinogenicity. It is also worth noting that a purely linear doseresponse would not be affected by stochasticity. We adapted the proportionality factor between R and excess number of molecules so that the lifetime excess cancer risk level estimated for a continuous exposure to 1.7 µg/m 3 BD is 10 -5 , as in the study by [START_REF] Higashino | Exposure and risk assessment of 1,3-butadiene in Japan[END_REF].

We performed simulations for exposure concentrations from 0 to 1.7 µg/m 3 , with a step of 0.1 µg/m 3 . For each concentration, we simulated, at steady state, the number of BD metabolites molecules in 100 cells exposed during 100 minutes.

Stochastic simulations in time-varying conditions

To study the system in a non steady-state situation (e. g. for a time-varying change in the exposure level) the exposure scenario simulated 9 hours a day to a BD concentration of 1.7 µg/m 3 , followed by 15 hours a day to a BD concentration of 0.25 µg/m 3 . This is typically what a factory worker would be exposed to in Japan [START_REF] Higashino | Exposure and risk assessment of 1,3-butadiene in Japan[END_REF]. For simplicity, holidays and week-ends were not accounted for.

We added a component (accounting for a liver cell) in parallel to the liver in the PBPK model.

It was not possible to implement time-varying deterministic concentrations in Jdesigner and

Dizzy and use the same methods as in section 2.4. Therefore, we had to simulate with Dizzy simultaneously the dynamics in all organs and in a cell. However, we could not track all the molecules in the body even for an environmental exposure level of BD. Therefore, we had to adapt the PBPK model equations implemented in Jdesigner and Dizzy, so that simulated organ concentrations were approximately the deterministic ones and so that Gillespie algorithm was unchanged for the cell.

Equation ( 10) is formally equivalent, relative to dynamics at the cell level, to the following equation: , fix at 0.01 the volume of the liver cell (both operations largely minimise stochasticity at organ level so as the contribution of the liver cell to the whole system) and divide F organ by 100, we have approximately the same kinetics at the organs level, without affecting the kinetics for Q_cell. The Gillespie algorithm is unchanged, because the reaction probability density function is unchanged. We indeed checked with Jdesigner and Dizzy that, once the model was adapted, the stochastic values of all organ concentrations differed by less than 5 % from their deterministic values during the simulations we performed. Using 10 and 0.1 instead of 100 and 0.01 in equation ( 13) would have led a large difference between deterministic values and stochastic ones. In contrast, using 1000 and 0.001 would have led to huge calculation times for Dizzy.

        × - × × ℵ × × =
We focused here only on EB in the liver. We performed simulations for 20 cells during the 15 hours in the general environment and reported cell concentration. The concentrations at time 0 (time at which the subject is just leaving industrial vicinity) are extracted from simulations with Jdesigner over 50 days. For the last 10 days, the values at time 0 only differed by 0.1 %, i.e. dynamics steady state was reached. Cancer risk was computed at each time step (every minute).

RESULTS

Stochastic simulations at steady state

The mean number of molecules in a cell, the standard deviation and coefficient of variation per organ for exposure concentration 1.7 µg/m 3 are presented in Table V. Stochastic mean values equal deterministic ones.

The stochastic variability is mainly influenced by the partition coefficient (which determines the mean number of molecules at steady state). The fact that marrow and kidney have comparable standard deviations but very different mean number of molecules suggests that low blood flow to the cell also tends to decrease variability relative to the mean value. The results obtained for metabolites in the liver are presented in Figure 2. There is a large variability of the number of metabolites in a liver cell (mean 2.91, standard deviation 1.71 resulting in coefficient of variation of 0.588). Simulations for exposure concentration 0.25 µg/m 3 led to the same conclusions relative to organ characteristics and coefficient of variation (data not shown).

Coupling with cancer model

Figure 3 shows that accounting for stochasticity when coupling the PBPK model and cancer model increased the excess of risk even below the deterministic threshold of 0.6 µg/m 3 . The dose-response curve appears to have two distinct regimes: a "quadratic" looking region and a

"linear" looking one, closed to the deterministic response.

Stochastic simulations in time-varying conditions

Figure 4 shows the results at dynamics steady state for metabolites in one liver cell for 15 hours at exposure concentration 0.25 µg/m 3 after 9 hours at exposure concentration 1.7 µg/m 3 .

The metabolite mean concentration quickly reaches low levels (half of the initial value is reached after 13 minutes), at which stochasticity is high and kinetics is slow.

The cancer model was used to assess the instantaneous excess risk of liver cancer during this period. There was no increased risk for times over 65 minutes, according to the deterministic approach. The sum over time of instantaneous excess risk values during this period outside the industrial vicinity corresponded to 12 minutes at steady state exposure at the high concentration, which is low compared to the 9 hours exposure. In contrast, when accounting for stochasticity, the mean sum of excess risk corresponded to 85 minutes of high exposure.

DISCUSSION

The simple algorithms outlined in the works by Gillespie [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF] permit the modeling of microscopic stochastic phenomena [START_REF] Haseltine | Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics[END_REF]). Gillespie's algorithms can be so expensive computationally that alternatives have been proposed to approximate exact simulations [START_REF] Haseltine | Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics[END_REF][START_REF] Rathinam | Stiffness in stochastic chemically reacting systems: the implicit tau-leaping method[END_REF]. For the present work, we chose to use the exact algorithm, for simulations could be performed in a reasonable amount of time with modern calculation stations, thanks to the progress that computer technology has still made in the past few years.

Coupling our PBPK model to a model of excess risk of cancer affected dose response curve. It looks curvilinear, with a breakpoint dose lower than the theoretical one. Mutagenesis data from [START_REF] Elhajouji | Indications for thresholds of chromosome nondisjunction versus chromosome lagging induced by spindle inhibitors in vivo in human lymphocytes[END_REF] showed experimentally comparable profile, with a slight increase in mutagenicity, then a highly significant steep increase for compounds likely to have a threshold for effects. We do not pretend to fully assess cancer risk associated to exposure to butadiene. Our study is a simulation study aiming at assessing whether or not stochasticity should be considered when assessing risk. Therefore, at least in the case of exposure to low doses of compounds having a threshold for carcinogenicity, stochasticity in kinetics should not be neglected.

In a paper on dose-response and threshold-mediated mechanisms in mutagenesis, [START_REF] Lovell | Dose response and threshold mediated mechanisms in mutagenesis: statistical models and study design[END_REF] points that "absolute" threshold are difficult to estimate from toxicity data due to background noise as a consequence of stochasticity. The concept of absolute threshold may therefore make no physical sense at all. Lovell is in favour of "pragmatic" thresholds, constructed through dose-response models (in particular biology-based dose-response (BBDR) models) coupled with knowledge about what level of response is biologically important, and what level is not. We believe that dose-response models construction should seriously consider stochasticity to get the best accurate representation of reality at low doses.

BBDR models would permit to derive a deterministic approximate threshold, then stochasticity studies would permit to derive a realistic one based on stochastic simulations with BBDR models coupled with PBPK models.

To our knowledge, this study is the first time that the dynamics of compounds at the cell level is studied in parallel to the dynamics at the body level. Coupling microscopic and macroscopic dynamics is a real challenge, because of the very large scale difference between numbers of molecules in cell and in organ. Here, we achieved the coupling of PBPK models with tools developed in the framework of systems biology through two different approaches.

In case of steady state for organs, concentrations in PBPK compartments were fixed and exact dynamic equations for butadiene in the cell was derived. In case of time varying concentrations in the organs, we adapted the equations to get an approximate simulation of the dynamics at the cell level. We were then able to assess the influence of stochasticity relative to dose-response in toxicant risk assessment, with constant or time-varying exposure concentration. This coupling is not limited to toxicology. It can be immediately generalized to the study of stochastic concentration of many compounds in the cell, like for instance pharmaceuticals and hormones which may be effective at extremely low concentrations [START_REF] Gurevich | Stochastic description of the ligand-receptor interaction of biologically active substances at extremely low doses[END_REF].

ACKNOWLEGMENTS

The authors acknowledge Sandrine Micallef for her useful comments on our work. We would also like to thank two anonymous reviewers who greatly helped to improve the quality of the manuscript. This publication was supported by the European Commission 6th Framework Bergerova (1983) and Van der Molen (1996). Default value is 0.049, which corresponds to "remaining organs" in Van der Molen (1996).

  , in the PBPK model equations, we multiply the concentrations in all compartments but the liver cell by 

Figure 1 .

 1 Figure 1. PBPK model for butadiene, as implemented in Jdesigner.

Figure 2 .Figure 3 .

 23 Figure 2. Distribution of BD metabolites number per cell in a liver cell at steady state at time 100 for 100 simulations performed with Dizzy. Exposure concentration is 1.7 µg/m 3 .

Figure 4 .

 4 Figure 4. Kinetics of butadiene metabolites in one liver cell for a man leaving industrial

  Figure 1.

Table III .

 III computed using cardiac output, percent blood flows per tissue mass and organ weights given in ICRP 2002 Pub 89 [14] (Table2.8 p18-19, Table2.39 p28, Table2.40 p29) or provided BD and EB-specific parameters. These have been taken from[START_REF] Brochot | Use of a chemical probe to increase safety for human volunteers in toxicokinetics studies[END_REF] for BD parameters and[START_REF] Brochot | Development of a physiologically based toxicokinetic model for butadiene and four major metabolites in humans: global sensitivity analysis for experimental design issues[END_REF] for EB parameter.

	by William & Leggett [15]. Parameter		Symbol	Value	Unit
	BD fat over blood partition coefficient			22	-
	Tissue or organ Symbol BD lung over air partition coefficient	Value P lung_over_air	0.653	-
	Adipose EB liver over blood partition coefficient F adip	0.564	P liver_EB	0.59	-
	Adrenals Metabolisation rate for BD into EB F adrenal	0.02	K met	0.3	L/min
	Brain	F brain	0.78
	Breast	F breast	0.00
	Gut	F gut	0.98
	Heart	F heart	0.35
	Kidney	F kidney	1.23
	Liver	F eport	0.45
	Lung	F total	6.72
	Bone marrow	F marrow	0.29
	Muscles	F muscle	1.11
	Others	F other	0.19
	Pancreas	F pancreas	0.065
	Skin	F skin	0.33
	Spleen	F spleen	0.19
	Stomach	F stomach	0.065
	Testes	F testes	0.004
	Thyroid	F thyroid	0.094

Table IV .

 IV Fat content for the various organs or tissues. These have been taken from Fiserova-

Table V .

 V Mean number of molecules in a cell, standard deviation and coefficient of variation per organ at steady state at time 100 (5000 for fat, 500 for marrow) for 100 simulations performed with Dizzy. Exposure concentration is 1.7 µg/m 3 .

	Organ Mean number of molecules in a cell Standard deviation Coefficient of variation
	Fat	379.3	10.0	0.026
	Marrow 82.3	4.49	0.0545
	Kidney 22.9	4.52	0.197
	Liver	15.2	3.2	0.211
	Spleen 13.6	3.86	0.284
	Tissue or organ	Value
	Adipose		0.859
	Adrenals	default
	Brain		0.11
	Breast		default
	Gut		0.065
	Heart		0.083
	Kidney		0.052
	Liver		0.049
	Lung		0.017
	Bone marrow	0.186
	Muscles		0.064
	Others		default
	Pancreas	0.105
	Skin		0.15
	Spleen		0.03
	Stomach		default
	Gonads		default
	Thyroid		default
	22

standard man of 1.76 m and 73 kg, using the organ weights given by the ICRP 2002 Pub 89 (p.18(T 2.8) p.19(T 2.9)). Density for the organs is supposed equal to 1 excepted for adipose tissues (density 0.9) and bones (density 2)