Roseline Bonnard

Thomas E Mckone

Keywords: Exposure assessment, multimedia models, probabilistic risk assessment, lead, IEUBK

Integration of the predictions of two models with dose measurements in a case study of children exposed to the emissions of a lead smelter

INTRODUCTION

Exposure models for chemicals are used around the world to assess the human health risks and to support regulatory decision making. For example, they are used to support policies for pollution control [START_REF] Bonnard | Estimation de l'impact des rejets des installations industrielles et utilisation de la démarche d'évaluation des risques à[END_REF]) such as the Hot Spots Program of California Environmental Protection Agency (OEHHA 2003), in which the risks linked to atmospheric emissions from existing facilities are estimated. They are also used in France for Registered Installations for Protection of the Environment (INERIS 2003), which requires a human health risk assessment study to be performed and included in all the environmental impact analyses for new or modified facilities. Exposure models are also employed to estimate the risks linked to manufactured substances, as required by the European regulation for new and existing substances (EC 2003(EC , 2004)). Many countries define safety limits with models for various exposure media such as food, water and air (WHO 1993(WHO , 2000 and national regulations), as well for soils, for which human health risk assessment is used to identify, prioritise and assess the need for remedial actions (US EPA, [START_REF] Ferguson | Assessing Risks from Contaminated Sites: Policy and Practice in 16 European Countries[END_REF].

Because of the critical role played by models, scientists are required to show that their predictions are reliable. Even though a "model can never be truly validated, (but only invalidated)" [START_REF] Oreskes | Evaluation (not validation) of quantitative models[END_REF], NCR 2007), it is still necessary to evaluate its fitness for use.

Part of this process is checking if models yield results matching the observations from field studies, and for the right reasons (that is by providing the correct cause-effect relationships between input and output) and capture this link with a sufficient level of precision.

In the context of industrial facilities and contaminated soils, health risk assessment studies employ multimedia exposure models to predict contaminations at a local scale and sometimes employ physiologically based pharmacokinetic (PBPK) models, when the exposure levels need to be expressed in terms of internal doses, as for lead, whose risk management threshold value is usually fixed at 100 µg/l in blood for children.

A number of different multimedia models have undergone empirical evaluation in the last decade, but these exercises were carried out at regional or continental scales [START_REF] Schwartz | Quality assurance of exposure models for environmental risk assessment of substances[END_REF] and most of the time, they were limited to environmental concentrations [START_REF] Kawamoto | Calculation of environmental concentration and comparison of output for existing chemicals using regional multimedia modelling[END_REF][START_REF] Armitage | Empirical evaluation of spatial and non-spatial European-scale multimedia fate models: results and implications for chemicals risk assessment[END_REF][START_REF] Luo | A multimedia environmental model of chemical distribution: fate, transport, and uncertainty analysis[END_REF]. In the case of the PBPK models, results have been evaluated with subjects in controlled conditions where there were either administrated doses or measured environmental concentrations [START_REF] Cohen | An arsenic exposure model: Probabilistic validation using empirical data[END_REF], Hogan et al. 1998[START_REF] Griffin | Application of a probabilistic risk assessment methodology to a lead smelter site[END_REF][START_REF] Biesiada | Blood lead levels in children: epidemiology vs. simulations[END_REF], Glorennec et al. 2007a), but none of these studies have tracked multimedia transport from a contamination source. Because of the difficulties in collecting for the same period and the same location sufficient data for the source emissions, environmental concentrations, activity patterns, biomonitoring data, and reconstructing past exposures, predictions of multimedia models have seldom been compared with biomonitoring data from the source of contamination. Such an analysis was conducted in a case of contaminated soils by polycyclic aromatic hydrocarbons [START_REF] Dor | Validation of multimedia models assessing exposure to PAH -The SOLEX study[END_REF], where urinary 1-hydroxypyren data collected from workers were compared to the levels predicted by several multimedia models. That study showed discrepancies among the models tested, in terms of magnitude of the output and of predominant exposure pathways. However, the authors carried out comparisons using a deterministic approach with a point estimate obtained from the different models. Because they did not know the precision of those point estimates, they could not determine to what extent these point estimates were truly different or even comparable.

The goal of this paper is to evaluate the performance (accuracy and precision) of two sourceto-dose modeling approaches, for a contaminated site by a secondary-lead smelter. Lead contamination was selected because lead is a contaminant often considered in health risk assessments for industrial sites. At such sites, lead can be found in various exposure media, a situation that requires consideration of multiple exposure pathways and provides the opportunity to assess exposure levels both by modeling and measuring blood lead. The comparison of lead measurements and lead modeling results can be carried out at several levels: first in environmental media (such as air or soil), second in exposure media (such as vegetables) and lastly in the human body. The intermediate points of comparison help to determine at which level discrepancies may appear between models and observations and whether the final results are correct for the right reasons. With this approach, one strives not only for good agreement between models and observations (exposure model predictions versus biomarker data) but uses the intermediate steps to confirm consistent hypotheses regarding the magnitude and variation among competing exposure pathways.

METHODS

The lead smelter is located in a French village. The local authorities have ordered many investigations in this village and have implemented a control of the environment media concentrations for several years. The data available were collected in order to carry out the model evaluations and several campaigns of measurements were performed to complete the dataset.

Two source-to-dose models, built-up from existing sub-models, were used for this study and run simultaneously by Monte-Carlo sampling from the ranges of parameter input values.

The data related to the site and the models used are presented in the following subsections.

The data not shown (measurement data and input parameters of the models) are detailed in a report available upon request [START_REF] Bonnard | Amélioration des méthodes d'estimation prospective des expositions des populations autour des ICPE[END_REF].

Study Village and Lead Emissions

The village has about 800 inhabitants. The facility producing lead began operations in 1970.

It recycles lead from engine-batteries. In relation to the facility, houses are located in a sector between the north-north-east axis to the south-west axis, with the major portion of them been located between the north-east and the south-east axes. The distance from the facility stacks ranges from 100 meters to 1 km. Currently, this smelter is the only facility in the village emitting lead in the atmosphere. From 2000 to the end of 2001, several measures were implemented to reduce the contamination linked to its activity. These measures resulted in a significant reduction of the atmospheric emissions of lead.

Child Biomonitoring Data and Exposure Survey

A biomonitoring campaigns for children's blood lead was carried out in this village in June 2002 by the local health administration. A questionnaire was filled out for each child at that sampling time to collect information on his or her exposure conditions. The questionnaire documented age, residence time in the village, places frequented, age of the home, proportion of fruit and vegetables consumed from the family garden, type of water (tap or bottled) consumed, and parents' place of work.

To protect the confidentiality of the children, we were not able to obtain the exact location of the house where each child dwelled. Only the distance to the facility was provided, with an indicator, going from 1 to 3, expressing the propensity of the wind blowing from the facility towards the child's house. In addition, a map was provided to us indicating three different sectors assigned a wind propensity value equal to 1, three others a value of 2, and two sectors a value of 3. Among all the children likely to be exposed, 67% took part into the study, namely 97 children. Among them, 19 had parents working in the facility. Because the local sanitary authorities have demonstrated a positive link between the children's blood lead and the parents working in the facility for this sample of data [START_REF] Schmitt | Nouvelle enquête sur l'imprégnation saturnine des enfants de Bourg-Fidèle[END_REF] and since no multimedia exposure model is currently able to take into account such a relationship, we only kept the lead blood data from the 78 children whose parents did not work in the facility.

Environmental Media Measurements

Air concentrations

Daily air concentration data are collected from ATMO Champagne (2002), a public association in charge of monitoring the air quality. An air sampler, measuring PM10, was placed in the village at 120 meters from the facility stacks. Air concentrations were measured by atomic absorption spectrometry (AFNOR, 1990).

Plant concentrations

A monitoring program has been conducted by a private company. Each year, it consists of the collection of about 20 samples of fruit or vegetables, from 6 or 7 volunteers' gardens, located between 80 to 600 meters from the facility. All the samples are prepared and washed before being analysed (Prost 2002[START_REF] Prost | Suivi de la qualité des sols et des végétaux autour de l'usine Métal Blanc à Bourg-Fidèle, GSC (Gestion des sites et sols contaminés), campagne 2002, campagne 2003, campagne 2004[END_REF].

Deposition to the ground surface

A campaign was conducted to measure monthly lead deposition on the ground at 15 locations in the village, using deposit gauges (ISO 1989). These measurements were made in October and November 2005. One duplicate device was used to check the reproducibility and 7 of the samples collected were distributed on the North-East axis relating to the facility, which corresponds to the preferential wind direction. The results of this measurement campaign, as reported in [START_REF] Bonnard | Amélioration des méthodes d'estimation prospective des expositions des populations autour des ICPE[END_REF], showed a clear decreasing level of deposition with the distance.

Soil concentrations

Several campaigns of soil lead concentration measurements have been performed in the village since 1998. However, because of the location and the depth of the samples, few data could be considered as representative of the population exposure. Thus, there were two additional measurement campaigns in the summer of 2006. The first consisted of 48 surface soil samples and 12 samples taken from the first 20 centimeters. All were composite samples assembled by mixing five sub-samples from the same area. The samples were not distributed according to a regular grid but taken in the village from locations close to dwellings and selected according to accessibility. The surface samples were collected in open areas and in the inhabitants' gardens, where the soil had not been disturbed. All of the the first-twentycentimeter samples came from the inhabitants' kitchen gardens. The soils were digested by aqua regia and analysed by inductively coupled plasma optical emission spectrometry (ICP/OES) according to the European standard analytical method (ISO 1996, CEN 2000).

The second campaign consisted of in-situ measurements made with a portable X-ray fluorescence device. Measurements were performed at 133 locations. In order to evaluate correlation of results between the two methods, some of the fluorescence measurements were made at the location where soil surface analytical samples were collected and 24 additional measurements by X-ray fluorescence were performed on soil samples prepared at the laboratory for extraction before ICP/OES analysis. During the program that monitored lead in vegetables cultivated by inhabitants, measurements of lead in the first twenty centimeters of soil were performed. They were added to our database. A geostatistical study, with a set of 156 data corresponding to the lead concentration in the surface soil and with a set of 39 data corresponding to the first-twenty-centimeter soil samples, was achieved.

Lead is characterized by a low mobility in soils and tends to accumulate in the surface layers [START_REF] Adriano | Trace elements in the terrestrial environment[END_REF], Kabatia-Pendias 1992, Alloway 1995). However, even if the quantities of lead collected in each gauge during our two-month-measurement campaign (while the stack emissions were higher than the 2002 to 2006 average ones) had been deposited on the ground for four years (between 2002 and 2006) and had not moved, that deposition of lead would correspond to an increase of the concentration in the first two centimeter layer of less than 5%. Given this low accumulation rate, we assumed that the values measured in the soils in 2006 are representative of the concentrations in 2002. Subject to the condition that no other source was added to the soils and the surface soil was not mixed with other soil layers, this assumption seems reasonable.

Source to Dose Modeling

Source to dose modeling was carried out with two alternative integrated models. Both of them link a multimedia fate model with a pharmacokinetic model. The first integrated model is process based and its purpose is to provide a regional mass balance to track lead from a source to receptor using generic parameters. The second model uses an empirical relationship to estimate the lead deposition on the ground and is intended to be both site and receptor specific.

Description of integrated model A

Model A consists of the CalTOX multimedia fate and multi-pathway exposure model [START_REF] Mckone | A multimedia total exposure model for hazardous-waste sites[END_REF][START_REF] Mckone | A multimedia total exposure[END_REF] integrated with the exposure-uptake-biokinetic (IEUBK) model (US EPA 1994, White et al 1998). We set up this model system in EXCEL© files, and ran probabilistic assessments with the software Crystal Ball [START_REF] Decisoneering | Crystal Ball 7.2[END_REF].

CalTOX is a fugacity-based mass balance model in which each environmental medium (air, soil, water, etc) is represented by a homogeneous compartment. IEUBK is the model recommended by the United States Environmental Protection Agency (US EPA) to estimate children's blood lead exposure. It permits calculation of the lead intake and uptake from the exposure media concentrations and then determines the blood lead levels using a biokinetic sub-model. To address the inability of IEUBK to propagate uncertainty and variability from lead intake, Syracuse Research Corporation developed the Integrated Stochastic Exposure (ISE) model (SRC, 2003) using the same conceptual model as IEUBK but allowing for stochastic simulations. We did the same, by developing our own version of IEUBK using EXCEL© and Visual Basic, so that we could link it directly to the stochastic outputs of CalTOX. We checked that our implementation of IEUBK gives exactly the same results as the original code. With this approach our blood lead prediction model could consider temporal variation of lead intake and probabilistic distributions for some inputs defined as point estimates in the ISE model.

We accounted for intake of lead through inhalation, soil ingestion, ingestion of local fruits and vegetables, and background exposure from non-local food consumption (Glorennec 2007b).

Based on atmospheric emissions and the distribution of lead concentrations measured in the first twenty centimeters of soil, model A was used with Crystal Ball to characterize the distributions of lead in the village air, in the surface soil, in the above-ground fruits and vegetables and the below-ground vegetables grown in the village, and finally in the blood of the children of the village.

Description of integrated model B

Model B has three components: (1) a regression model that calculates lead deposition according to the location of the household in the village, (2) a multimedia exposure model that calculates local lead concentrations in environmental and exposure media based on the deposition rate, the surface soil concentration and the soil concentration in the first twenty centimeters and (3) our version of IEUBK recoded in EXCEL©.

Model B was built to estimate lead exposures for precise locations, an approach often required for contaminated sites to distinguish the areas requiring remediation from those which do not. Unlike model A, which is aimed to capture the distribution of the concentrations over the whole village, model B enables to represent the exposure of individuals. Because the exact location of each child's house was not known, we had to identify it basing on the concordance among three sources of information: (1) the child's distance from the facility, (2) the sector where the house was located and (3) a bird's-eye view zonal map of the village. In some cases, only one location of household corresponded to information derived from these three sources, in other there were several possibilities.

Because of this difficulty, blood lead level was modeled for a subset of nine children for which we had identified only one or two possibilities of housing location. If there were two possible locations for one child's house, two predictions were made for that child's blood lead based on a different soil concentration and a different lead deposition on the ground. The children of this subset had contrasting exposure conditions (high, low or no consumption of home-grown vegetables, close or far from the facility, low or high blood lead level) (Table 1).

We use the regression model in place of a classical gaussian air dispersion model because we found the accuracy of these models, which predict deposition and air concentration, at each point of a grid, to be inappropriate in regards of the estimation of fugitive emissions (see section called Source term definition). Instead we fitted deposition data collected during the sampling campaign (see section called Deposition to the ground surface) using a statistical relationship, taking into account the distance (D i in meters) of a specific location i from the facility, the frequency (W ij unitless) of the wind blowing from the facility towards a specific location during month j and the magnitude of stack emissions (E j in g/month). We used Statistica (Statsoft 1999) for this analysis. We obtained the higher correlation coefficient and the lower residual standard deviation with the following equation by ordinary least squares:

Ŷ ij = log dep ij = 0.35 -1.7 10 -3 D i + 1.31W ij + 1.4 10 -5 E j ( 2 = 0.93)
where dep ij = total deposition of lead on the ground at location i during month j (in mg/m 2 /month).

We used a Student's T-test to test if a coefficient of this relationship is significantly different from zero. Although the equation above is based on data collected for only two months, all coefficients were found to be statistically significant. In a good regression model, residuals (difference between the observed values and the predicted values) are normal and homoscedastic (residuals have the same variance). The residuals of the above equation appeared normal and homoscedastic.

We characterized the uncertainty of the predicted deposition (Ŷ) using a Student's distribution, with 26 degrees of freedom, multiplied by the estimated deviation of Ŷ [START_REF] Foucart | L'analyse des données, Mode d'emploi[END_REF]. We have 26 degrees of freedom because we have 30 data to define a relationship that has 4 coefficients. The estimated deviation of Ŷ is a T distribution based on the residual estimation error from the regression model. The output of this relationship along with the error were used as inputs in the multimedia exposure model.

The media concentrations were calculated from the equations described in the Human Health Risk Assessment Protocol (HHRAP) document (US EPA 2005). However, a supplemental soil layer (between 0 to 10 cm from the surface) for predicting the forage contamination from

Values used for model parameters and inputs

Here we describe how we obtained model inputs. First we consider the period for which we had to model the doses adsorbed by children to estimate their blood lead level at the time of the biomonitoring campaign. Next we consider the case of source terms, lead chemical properties, environmental parameters, intake and uptake parameters. We developed probability distributions for each input parameter of both models.

Period for reconstructing children's intake

The half-life of lead in blood is estimated to be one month. Although lead is also stocked in bones and slowly released in blood from there, it has been shown that a reduction of child exposure brings out a rather fast reduction of the blood lead (Declercq andLadrière 2004, Khoury et al 2003). The conceptual model implemented in the IEUBK software also yields a rapid reduction of the blood lead, once intake is stopped. We ran the IEUBK model for the case of a child with a cumulative ingestion dose 3 times higher than the estimated background intake in France for 3 years and then a dose equal to the background intake (Glorennec 2007b). We found that one year after reducing the intake, the blood lead would be only 2 percent higher than the one calculated with the background intake for 4 years.

Therefore, the efforts to reconstruct the intake doses were focused on a period starting from spring 2001 to the date of the campaign of blood lead measurements. Beforehand, the exposure conditions were assumed to be the same as those between January and June 2001.

Source term definition

Model A

A lognormal distribution was fitted to the observed lead concentrations in the first twenty centimeters of soil using a Z-plot chart. The resulting distribution with an arithmetic mean of 261 mg/kg and a standard deviation of 249 mg/kg was assigned to the root zone soil concentration parameter of CalTOX.

We estimate smelter stack emissions from monthly sampling of particule emissions, the measurement of the lead concentration in the emitted particles and the duration of operation. The smelter operators try to limit the fugitive emissions by keeping the pressure inside the buildings lower than the outside pressure, washing outdoor surfaces daily and stocking raw materials products and by-products indoors. Nevertheless, preliminary calculations performed with an air dispersion model (CERC 2000) showed that the air concentrations and particle fallout recorded downwind could not be explained purely by the recorded stack emissions [START_REF] Bonnard | Amélioration des méthodes d'estimation prospective des expositions des populations autour des ICPE[END_REF]. Therefore, we had to include fugitive emissions. We estimated them using The atmospheric emissions were estimated from the statistical model described above (Equation 1) for each location and period used to represent the atmospheric deposition on the different kinds of home-grown plants.

Chemical properties for lead

Model A With the exception of the particle-water partition coefficients, we used the default distributions provided by CalTOX for chemical properties-water solubility, partition coefficients, and bioconcentration factors. We found that the particle-water partition coefficient (Kd) in CalTOX was high relative to ranges found in other papers. Since Kd for lead is known to depend on the soil pH, we used the recommendations from the Office of Air and Radiation (US EPA 1999) and the pH values collected in the garden soils of the village to define a triangular distribution. The measured pH going from 5.3 to 8.3, the minimum and maximum values for the closest ranges of soil pH, given by the Office of Air and Radiation, were used as minimum and maximum of the distribution (150 L/kg and 23,279 L/kg respectively). The mode was assigned a value of 1731 l/kg, which was calculated from the modal value of measured pH and the Office of Air and Radiation's relationship between lead Kd and pH.

Model B

The parameters required for model B are the soil-plant transfer coefficients. To obtain values, we made a review of the primary literature and fitted the resulting dataset of selected values to a lognormal distribution (Table 2,Bonnard 2008).

Environmental parameters

Model A

For model A we used the default distributions provided by CalTOX for all environmental parameters except area, rainfall, wind speed, and temperature. We replaced distributions for these parameters in CalTOX by local data. In addition we replace the CalTOX root zone soil depth with a mean value equal to 0.2 m to be consistent with the measured data for the village.

Model B

Most of environmental parameters required in model B are different from those of CalTOX.

They mainly concern the description of the various categories of plants. Distributions were defined from the literature (Table 2, [START_REF] Bonnard | Amélioration des méthodes d'estimation prospective des expositions des populations autour des ICPE[END_REF].

Intake parameters

From birth to the age of 7, we defined plant consumption and lead exposure from local food consumption for seven age classes based on French surveys [START_REF] Boggio | Consommation alimentaire des nourrissons et des enfants en bas âge en France en 1997[END_REF][START_REF] Volatier | Enquête I individuelle et nationale sur les consommations alimentaires, Editions Technique & Documentation[END_REF].

For soil ingestion, we used the cumulative distribution provided by ISE with a median equal to 135 mg/d (US EPA 1994EPA , 1999)). We note that other distributions with lower median values were available from empirical data published by [START_REF] Thompson | Parametric distributions for soil ingestion by children[END_REF] and [START_REF] Stanek | Daily soil ingestion for children at Superfund site[END_REF][START_REF] Stanek | Soil ingestion distributions for Monte Carlo risk assessment in chidren[END_REF]. However, the IEUBK model was calibrated with a point estimate equal to 135 mg/day. On the other hand, the lognormal distribution given in the ISE's model for soil ingestion, has a standard deviation that appears too low relative to what is available in the literature. For example, the ratio of the 95 th percentile to the 5 th percentile in the ISE distribution is less than 2, indicating almost no inter-individual variability. We did use the ISE age-specific weighting factors for soil ingestion in our recoded version.

For model B, we use observed individual data for each child from the study population to obtain the age, the seasonal variation in the percentage of home-grown fruits and vegetables consumed, and the date of moving into the village.

Uptake parameters

The bioavailability of lead in food and the relative bioavailability of lead in soil were revised to better reflect data documented in IEUBK for these factors.

The range reported by the USEPA in the IEUBK guidance (US EPA 1994) for absorption of lead from diet is 42 to 53%. These values correspond to the 40 th and 63 th percentiles of the lognormal distribution defined in the ISE model for this parameter (mean = 0.5, standard deviation = 0.2). A lower value was used for the standard deviation (0.05 instead of 0.2) to get a distribution more in accordance with the range of data given above (in these conditions, the low and high values given by the US EPA correspond to the 5 th and 75 th percentiles of the new distribution defined for absorption of lead from diet).

For the absolute bioavailability of soil lead, ISE gives a point estimate of 0.3, but no distribution is defined. We used the values collected by [START_REF] Ruby | Advances in evaluating the oral bioavailability of inorganics in soil for use in human health risk assessment[END_REF] for lead uptake from soils (data from wastes and mines were eliminated), to define a normal distribution with a mean equal to 0.6 and a standard deviation equal to 0.2 for the relative biovailability of soil lead. The mean of the soil lead absolute bioavailability obtained is then close to the point estimate used in the IEUBK and ISE models.

Probabilistic computation

All of our calculations were based on Monte Carlo simulations using 5000 runs with latinhypercube sampling. We used the Monte Carlo results to develop cumulative probability distributions that can be compared to cumulative distributions of measured concentrations in the media.

RESULTS AND DISCUSSION

Model A

Because no significant statistical relationship between the blood level and the age could be identified [START_REF] Schmitt | Nouvelle enquête sur l'imprégnation saturnine des enfants de Bourg-Fidèle[END_REF], we analysed the model's output as the mean of the blood lead from 6 months to 7 years old.

Figure 1 shows the probability plot charts for the (base 10) logarithms of the observed and predicted concentrations in the different exposure media. On these charts, concentrations in the y-axis are represented against the cumulative probability expressed as the number of standard deviation from the geometric mean. For blood lead, figure 1e shows the concentrations calculated by the integrated model and those predicted with the new version of IEUBK from the measurements performed in the samples of surface soil, fruit and vegetables collected in the village.

The observed values in air, ground soil, exposed produce and blood lead tend to follow lognormal distributions. However, concentrations above 500 mg/kg in the ground soil deviate from a lognormal distribution (these values correspond to a small area close to the facility and without dwelling). The distributions of the calculated concentrations in air and ground soil also follow lognormal distributions, whereas the distributions of the predicted concentrations for the protected produce, the exposed produce and in particular for the blood lead are characterized by two different slopes, indicating bimodal distributions.

Except for the exposed produce, the median values predicted by the integrated model are in good agreement with the observed medians. The predicted median-to-observed median ratios are between 0.5 and 1.3 for these media. But for the exposed produce, the ratio is equal to 18.

Regarding distributions, the observed data for ground soil are well captured until 500 mg/kg. The ranges of the observed and predicted values for the air concentrations are close and the two distributions are quite similar. But for protected produce, the distribution of the observed data is underestimated and above all the distribution of the observed concentrations for exposed produce is overestimated. For blood lead, the variability of the calculated concentrations is higher than that of the observed concentrations. In the case of concentrations predicted from the measurements performed in the exposure media, blood lead concentrations are underestimated up to the 90 th percentile, then they are overestimated. The values predicted from the integrated model are higher and the difference between the two distributions of calculated concentrations increased with percentiles.

According to figure 2, even though emissions are primarily to air, inhalation of smelter emissions is not a major exposure pathway, because most of the lead emitted to air deposits to soil and from there is transferred to the children by various ingestion pathways.

When blood lead level is predicted from the exposure media concentrations, soil ingestion appears as the predominant pathway in the upper percentiles of the distribution of the predicted blood lead concentrations, as in former studies conducted near smelters [START_REF] Roels | Exposure to lead by the oral and the pulmonary routes of children living in the vicinity of a primary lead smelter[END_REF][START_REF] Gulson | Source and pathways of lead in humans from the Broken Hill mining community; an alternative use of exploration methods[END_REF][START_REF] Landrigan | Lead levels, home dust, and proximity to lead smelters[END_REF]). However, when blood lead level is calculated by the integrated model, the ingestion of exposed produce appears as the main pathway in the last quarter of the blood lead distribution. Therefore, the overestimation of the exposure produce concentrations by the integrated model contributes to the overestimation of the upper part of the distribution of blood lead concentrations calculated by this model.

Analysis of the contributions to exposed produce contamination showed that in the first part of the distribution, transfer from air and from surface soil to plants are higher than transfer from root-soil, and then the situation is reversed. In order to understand the reasons for the overestimation of the exposed produce concentrations, we carried out a sensitivity analysis with Crystal Ball. The parameters with the highest contributions to the variance are listed in Table 3 in descending order.

The distributions of the first two parameters with the highest sensitivity (particle-water partition coefficient and root-soil concentration) were defined by taking into account specific data from the studied site. The third parameter is the leaf wet density (rho_leaf). The default value, which was used in the computations, is equal to 820 kg/m 3 . With a volume fraction of water in leaf (beta_leaf) equal to 0.5, as provided in CalTOX, we get a value of 0.61 for the mass water fraction (0.5 x 1000 / 820), whereas the data for the exposed produce (leafy vegetables) collected in the gardens ranged from 0.77 to 0.98. So, we determined a new value for the leaf wet density on the basis of the mean of the leaf mass water fraction measured in the collected samples of exposed produce (0.91). We found a value equal to 550 kg/m 3 . It remains in accordance with the typical values used in models for leaf wet density, which are comprised between 500 and 900 kg/m 3 [START_REF] Riederer | Estimating partitioning and transport of organic chemicals in the foliage atmosphere system[END_REF][START_REF] Paterson | A fugacity model of chemical uptake by. plants from soil and air[END_REF][START_REF] Bacci | Chlorinated dioxins -Volatilization from soils and bioconcentration in plant leaves[END_REF][START_REF] Trapp | Plant contamination: modelling and simulation of organic chemical processes[END_REF][START_REF] Maccrady | Uptake and photodegration of 2,3,7,8 tetrachlorodibenzodioxins sorbed to grass foliage[END_REF].

As a result, the distribution obtained for the exposed produce concentrations (figure 3a,curve 3) is closer to that of the observed concentrations. The predicted concentrations were divided nearly by a factor 2, except at the lower tail of the distribution. The air, soil and protected concentration distributions were not altered by this change, but the variance in the predicted blood-lead concentrations (calculated by the integrated model from the predicted concentrations in produce) was reduced (figure 3b, curve 3).

For the other parameters pointed out in the sensitivity analysis, we found no new data justifying any modifications of the default values and significantly improving the prediction of the exposed produce concentrations.

But, we needed to consider that fruits and vegetables had been washed before the analyses and recognized that the impact of washing is not considered in the model. Unfortunately, the extent to which such a treatment reduces the contribution of air deposition and rain-splash to fruit and vegetables contamination is not well known. In our model, if the contributions of both these pathways are divided arbitrarily by ten, it appears clearly that the lower part of the predicted distribution (figure 3a, curve 4) is closer to that of the observed concentrations data.

That means that the variability of the measurements is then captured better by the model. The model over-predicts the concentration in exposed fruits and vegetables by a factor 4 at the median, under these assumptions. The distribution of the blood lead concentrations predicted from the integrated model is then very close to that predicted by the new version of IEUBK from the exposure media measurements (figure 3b, curves 4 and 5).

Nevertheless, these two distributions still deviate from the distribution of the observed blood lead at the upper percentiles and according to the integrated model, the ingestion of homegown plants remains the main contributor to the blood lead level above the 90 th percentile. This is not in accordance with the result of the univariate analysis of observed data showing a higher mean of blood lead concentration for no consumers of homegrown plants (m=65 µg/l, n=46) than for consumers (m=45 µg/l, n=32), even if for the group of homegrown plant consumers, the blood lead level increases with the fraction of homegrown plants in the diet. The discrepancy between observed and predicted data may be due to the fact that in the sample of children tested, those with the highest consumption of kitchen garden foods were those living the farthest away from the facility, thereby with the lowest lead concentrations in plants. In other words, the actual consumption of home-grown produce in the village and its contribution to blood lead level at the upper percentiles is overestimated by the model, which assumes a homogeneous distribution of the home-grown produce consumptions within the village.

Model B

The results of model B are provided in Figure 4 for nine children with one or two possible housing locations. We can see that the measured blood lead levels are between the 5 th and the 95 th percentiles of the model results, with the exception of child 43. The measured blood-lead levels correspond to various percentiles of the distribution given by the model, going from the 25 th to the 88 th percentile according to the different children. Measured values of blood lead above the predicted median ratios are between 0.7 and 1.7 (excluding the results corresponding to child 43).

Model B was able to capture the blood-lead-level ranges for the various conditions tested, except for child 43. The discrepancy between the model output and the blood lead level measured for this child may be explained by the fact that this child lived in an old house built before the ban of lead for interior paint. The univariate analysis of the data of the child biomonitoring and of the exposure survey showed a higher mean of blood lead concentration for children living in houses built before 1948 than for children living in more recent houses (DDASS 2002). But, this exposure pathway is not taken into account by the model. It has to be noted that using either the [START_REF] Thompson | Parametric distributions for soil ingestion by children[END_REF] or [START_REF] Stanek | Daily soil ingestion for children at Superfund site[END_REF][START_REF] Stanek | Soil ingestion distributions for Monte Carlo risk assessment in chidren[END_REF] distributions for soil ingestion, instead of a distribution centered on the value used to calibrate IEUBK for soil ingestion, results in lower estimates of blood lead concentrations, especially when the soil ingestion is a significant contributor to overall exposure. For example soil ingestion dominates when the children consume few or no homegrown foods, and when the soil lead concentration is high. In the sub-sample tested here, the geometric mean of the measured median above the calculated value ratio would increase from 1.1 to 1. Many studies showed a correlation between blood lead level and soil lead [START_REF] Xintaras | Analysis Paper : Impact of lead-contaminated soil on public health[END_REF][START_REF] Landrigan | Lead levels, home dust, and proximity to lead smelters[END_REF][START_REF] Mielke | Soil is an important pathway of human lead exposure[END_REF]. However, for instance, for child 71, who did not consume kitchen garden plants and who lived close to the facility, model B indicates that exposure due to ingestion of no local food is slightly higher than exposure due to soil ingestion. This is certainly due to the fact that soil lead concentration at the location of his house is not high (200-230 mg/kg) compared to values measured in these studies. For child 28, who ate a large amount of home-grown plants and lived relatively far from the facility, the ingestion of homegrown produce appears as the largest contributor to the blood lead level. Wilhem et al. (2005), who studied the dietary intake of lead by two groups of children living in an industrial region, found no difference between those eating a large part of produce from the families own vegetable garden and those consuming exclusively food from the supermarket. But in this case, comparison was made at a region scale without an identified local source of lead. Here, measures of concentrations in plants from child 28's garden would be necessary to confirm the role played by the ingestion of homegrown plants.

The uncertainties in the modeled blood lead predictions for model B are reflected by the ratios of the 95 th percentile to the 5 th percentile, which are between 2.6 and 4.9. The widest 90% interval is given for child 28. In this case, the root soil concentration is the highest contributor to the variance. This high contribution is due to a lack of precision in the estimation of the root soil concentration in the child's living area. Increasing the sampling efforts in soils of the farthest areas of the village could permit the reduction of the overall uncertainty. Nevertheless, model B is able to predict blood lead distributions with lower variance than model A (the ratio of the 95 th percentile to the 5 th percentile is equal to 10.2 with model A). Indeed, model A has to integrate inter-individual variability in the input distributions, while model B took it into account by a better description of the influence of the distribution of lead deposition and individual characteristics.

Limitations of the comparison between predicted data and observed data

Finally, we acknowledge that the exercise of comparison made with both models may be limited by the quality of data used to perform it.

Soil lead concentration appeared as the most sensitive parameter for the blood lead predictions performed. However, it was measured in bulk soil and not in the finest particle fraction, more likely representative of lead concentration in soil ingested by children.

Campaigns for soil lead and blood lead concentrations were not conducted at the same time.

(so, there is a possibility that soil was moved in some areas of the village between the two campaigns) and soil lead was not measured in children's garden soil, but estimated by a geostatistical study. Air emissions used in model A were also estimated from indirect methods.

On the other hand, the lack of measurement data in the village for tapwater and lead bioavailability, which is usually recognized as a determining factor for estimating blood lead level [START_REF] Mushak | Uses and limits of empirical data in measuring and modeling human lead exposure[END_REF], should not represent a major fault. For drinking water, the analysis of observed data showed that there is no significant difference between blood lead level of children drinking tapwater and those drinking bottled water. Lead bioavailability is not either a very sensitive parameter in the previous calculations. The absolute lead bioavailability and the relative bioavailability of soil lead together contribute to 3% of variance in model A and to 7% in model B at the maximum.

CONCLUSIONS

Both models, built for estimating human health risks linked to future emissions of facilities, may provide useful insight regarding source-to-dose relationships for the exposure population. Model A may be used to address general trends and model B is more focused on capturing individual exposures. In considering the results presented above, there is an issue that cuts across our evaluation of the individual models and data. In this process, we recognized that not only did we gain insight into the two models, but by making a systematic comparison of two models with a specific case study, we gained important insights about lead exposure in these communities that could not be obtained from either model or the environmental/biomonitoring data used in isolation. Thus we conclude that using the two models together with population and site-specific measurements of blood lead and environmental conditions provides for this population's lead exposure patterns a level of understanding that could not be achieved with either model alone or with just the blood surveillance measurements.

Model A provides an ability to capture the distributions of lead in the media and the blood lead concentrations at a local scale. It offers the opportunity to identify the potential pathways of concern and to characterize also the distribution of the exposure level. From the initial results of the model, mainly used with default values, the information from the sensitivity analysis and the pathway contributions evaluation help to identify the parameter improvements needed to calibrate the model to this case. However, predicted concentrations in exposed vegetables may be a potential source of overestimation of blood lead level.

Model B provides outputs that can be targeted to specific areas and exposure conditions. It was set up to answer some questions, such as "in which zones or from which distance may the consumption of home-grown produce raise a health concern ?" Because of this feature, model B can predict blood lead distributions with lower variance than model A. Precision of assessment is important for discriminating those situations that require management measures from those which do not, especially if the risk management threshold is low. Although, the results provided by model B, for contrasting exposure conditions, are in accordance with the measurements (apart from a case that may involve an exposure pathway that is not considered in the model), the low number of cases tested here does not enable to conclude on the performances of model B at the moment.

(Plantago lanceolata L.). Contaminated soils:

International conference on the biogeochemistry of trace elements, 15/05/95, Paris, France. Proceedings 85: 325-335

  soil and additional classes of plants (fruit, leafy vegetables and fruit vegetables) were introduced, as well as equations for plant contamination by irrigation and deposition of resuspended soil particles. The equations used for these pathways are those given in the model ERMYN (US DOe 2003a). The concentrations in air, surface soil and plants yielded by the multimedia exposure model were used as input of the new version of IEUBK to obtain the child blood lead distributions.

  three methods (inverse modeling, predictions based on the emission factors and the production rate, and a calculation based on the indoor air concentration in the buildings of the factory and the ventilation rate) and obtained consistent estimates. With inverse modeling, we calculated the volumetric flux required to get the best correlation between the air dispersion model and the lead deposition measured during the campaign conducted in October and December 2005 and obtained a fugitive emissions estimate between 20 and 60 kg/month. Emissions factors for this type of facility and information on its production rate give fugitive emissions between 25 and 73 kg/month. Data on lead concentration measured in indoor air of the facility building with an assumption of a ventilation rate of 10 times per hour give a fugitive emissions estimate of 34 kg/month.To cover the whole growing period of the local fruit and vegetables consumed till the date of the campaign of blood lead measurements, stack lead emissions were estimated from the data recorded from April 2001 to June 2002. As information on fugitive emissions is poor, a triangular distribution was assigned to the total atmospheric emissions from the estimates obtained for the stack and the fugitive emissions during this period. The mode is equal to 1,200 g/d, the minimum to 750 g/d and the maximum to 2,800 g/d.For comparing the predicted vegetation concentrations with the observed ones, the emissions were estimated over the vegetation exposure period to fallout from the facility (from April toJuly, in 2002July, in to 2005)). The total atmospheric emissions were defined by a triangular distribution with a mode of 780 g/d, a minimum of 450 g/d and a maximum of 1,450 g/d.Model B Spatial soil mean concentration in the land associated to each child's house was estimated by kriging. Cumulative distributions were defined for the surface soil and the root-zone soil owing to the turning band method(Chilès 1999[START_REF] Bonnard | Amélioration des méthodes d'estimation prospective des expositions des populations autour des ICPE[END_REF]).

  3, if the distribution used for the quantity of soil ingested was replaced by Thompson et al.'s distribution and it would reach a higher value with Stanek et al.'s distribution.

Figure 1 :

 1 Figure 1 : Distributions of the measured concentrations and of the concentrations calculated

Table 1 :

 1 Cases investigated with the model B (because of uncertainty in the location of some 796 children's houses, there are two entries for some children based on two different assumptions 797 for soil concentrations and lead deposition around their houses) 798

	child number	age	blood lead	distance to	% home-	% of home-	surface soil lead
	identification	(month)	(µg/l)	the stacks	grown	grown	concentration in the
				(m)	fruit	vegetables	location of the child's
							home (mg/kg)
		33	67.1	250	0	50	173
		47	83.1	900	75	75	96
		47	83.1	900	75	75	84
		74	41.5	900	62.5	75	96
		74	41.5	900	62.5	75	84
		52	31.9	450	10	20	112
		52	31.9	450	10	20	103
		65	44.7	150	0	0	203
		65	44.7	150	0	0	227
		37	44.7	300	0	0	261
		37	44.7	300	0	0	174
		81	22.3	1000	0	0	66
		21	108.7	450	0	0	112
		21	108.7	450	0	0	103
		71	99.1	200	0	0	493
	799
	800

Table 3 :

 3 Input parameters with the highest sensitivity on exposure produce concentrations in

	model A
	Parameters	Symbols	Rank correlation coefficients between the
			input parameter and exposure produce
			concentration
	Particle-water partition	Kd	-0.44
	coefficient
	Concentration in root-zone soil	Cs	0.41
	Leaf wet density	rho_leaf	0.37
	Primary production dry	veg_prod	0.31
	vegetation
	Rainsplash	rainsplash	0.21
	Leaf surface erosion half-life	Thalf_le	0.14
	Stem wet density	rho_stem	-0.12
	Annual average precipitation	rain	0.11
	Ambient environmental	Temp	-0.09
	temperature
	Land surface runoff	runoff	0.09

ACKNOWLEDGEMENTS

The authors thank the local administration (Direction Régionale de l'Industrie, de la Recherche et de l'Equipement de Champagne-Ardennes and Direction Départementale des Affaires Sanitaires et Sociales des Ardennes), the director's facility, ATMO Champagne and the company Gestion des sites et des sols contaminés (GSC) for communicating their data, as well as the studied village inhabitants who permitted us to take samples from their garden and the mayor for his help. We also thank F. Bois of Institute of Industrial Environment and Risks for his review of the paper and M. Sohn of Lawrence Berkeley National Laboratory for his guidance and input on this work. R. Bonnard's work was supported by the French Ministry of Ecology and Sustainable Development (Budget Civil de la Recherche et du Développement, grant BRCD-03-DRC-06). T. McKone was supported in part by the US Environmental Protection Agency through Interagency Agreement DW-988-38190-01-0 carried out through the US Department of Energy contract Grant No. DE-AC02-05CH11231. McKone was also supported by Cooperative Agreement Number U19/EH000097-03 from the US Centers for Disease Control and Prevention (CDC).

 -3, 3.1 E-3 , 6.1 E-4, 1.9 E-2) Root-vegetables 4.5 E-2) a U (x1, x2) specifies a uniform distribution with minimum x1 and maximum x2. e LN (x1,x2,x3,x4) specifies a lognormal distribution with mean x1, standard deviation x2, lower truncation limit x3 and upper truncation limit x4.


Guillaume Fayet

Patricia Rotureau
email: patricia.rotureau@ineris.fr

Laurent Joubert

Carlo Adamo

On the prediction of thermal stability of nitroaromatic compounds using quantum chemical calculations

Keywords: Thermal stability, explosibility, property prediction, Quantitative Structure-Property Relationship (QSPR), quantum chemistry

This work presents a new approach to predict thermal stability of nitroaromatic compounds based on quantum chemical calculations and on Quantitative Structure-Property Relationship (QSPR) methods. The data set consists of 22 nitroaromatic compounds of known decomposition enthalpy (taken as a macroscopic property related to explosibility) obtained from differential scanning calorimetry. Geometric, electronic and energetic descriptors have been selected and computed using Density Functional Theory (DFT) calculation to describe the 22 molecules. First approach consisted in looking at their linear correlations with the experimental decomposition enthalpy.

Molecular weight, electrophilicity index, electron affinity and oxygen balance appeared as the most correlated descriptors (respectively R²= 0.76, 0.75, 0.71 and 0.64). Then multilinear regression was computed with these descriptors. The obtained model is a six-parameter equation containing descriptors all issued from quantum chemical calculations. The prediction is satisfactory with a correlation coefficient R 2 of 0.91 and a predictivity coefficient R² cv of 0.84 using a cross validation method.

Introduction

An important initial step for the management of industrial risks consists in identifying and determining as soon and as accurately as possible explosive abilities of molecules, whatever their final use. Moreover, the explosive intrinsic property of a substance ranks at top of physicochemical hazards that may be feared from the use of a given chemical [START_REF] Etchells | Why reactions run away[END_REF]. This is a reason why, in addition to the complexity with which this hazardous property is triggered in real case, the experimental approach has remained the "golden" way to assess hazards. Keeping on this unique approach would however be a real burden for the industry with the requirements of the new regulatory framework REACH 1 and its tool the GHS 2 given the tremendous number of substances (up to 30000) that might require a new assessment of hazardous properties. Another interest for predictive methods is to address needs of screening processes of substances that are applied to search for a targeted chemical activity (e.g. phytotoxicity, medical efficiency for a given illness, reactivity…) at the R&D level.

A number of early works are worth being mentioned in the field of hazard prediction. At first, some methods of prediction based on thermodynamic concepts have been developed. Significant weaknesses in those methods lie in the chemical thermodynamic and energy release evaluation (CHETAH) [START_REF] Shanley | A review of ASTM CHETAH 7.0 hazard evaluation criteria[END_REF] and the calculated adiabatic reaction temperature (CART) [START_REF] Melhem | On the estimation of hazard potential for chemical substances[END_REF]. Grewer [START_REF] Grewer | The influence of chemical structure on exothermic decomposition[END_REF][START_REF] Grewer | Prediction of Thermal Hazards of Chemical Reactions[END_REF] proposed another way to predict the thermal stability by considering the influence of the chemical structure on the decomposition of nitro compounds. Keshavarz's works have also to be noted as the elemental 1 REACH for "Registration, Evaluation and Authorization of Chemicals" is a new European regulation published by the European Commission in December 2006 and entered into force on 1st June 2007. All manufacturers, importers and suppliers of chemicals must identify and manage risks linked to the substances they manufacture and market. http://echa.europa.eu/ 2 GHS for "Globally Harmonized System of classification and labeling of chemicals" is a United Nations system which aims to identify hazardous chemicals and to inform users about these hazards through composition of substances is used to predict various properties of energetic materials such as the impact sensitivity [START_REF] Keshavarz | Prediction of impact sensitivity of nitroaliphatic, nitroaliphatic containing other functional groups and nitrate explosives[END_REF]. Saraf et al. [START_REF] Saraf | Prediction of reactive hazards based on molecular structure[END_REF] outlined the pertinence (in terms of way of investigation) of screening tools based on the identification of relationships between chemical structure and thermal stability for nitroaromatic compounds and we have made up our mind to explore this route.

An alternative tool for the prediction of chemical hazards is the Quantitative Structure Activity/Properties Relationship (QSAR/QSPR) methods. Nowadays, they are generally limited to the scope of toxic property screening (i.e. the nitrobenzene molecule [START_REF] Agrawal | QSAR prediction of toxicity of nitrobenzenes[END_REF]).

Indeed, the first applications of these methods based on statistical analyses have mainly concerned biology [START_REF] Gao | Comparative QSAR analysis of estrogen receptor ligands[END_REF][START_REF] Winkler | The role of quantitative structure-activity relationships (QSAR) in biomolecular discovery[END_REF], toxicology [START_REF] Selassie | Comparative QSAR and the radical toxicity of various functional groups[END_REF][START_REF] Bradbury | Quantitative structure-activity relationships and ecological risk assessment: an overview of predictive aquatic toxicology research[END_REF] or drug design [START_REF] Garg | Comparative Quantitative Structure-Activity Relationship Studies on Anti-HIV Drugs[END_REF][START_REF] Grover | Quantitative structure-property relationships in pharmaceutical research -Part 1[END_REF][START_REF] Grover | Quantitative structure-property relationships in pharmaceutical research -Part 2[END_REF]. However, their interest has been growing up in recent years for other physico-chemical activities [START_REF] Katritzky | QSPR: the correlation and quantitative prediction of chemical and physical properties from structure[END_REF][START_REF] Katritzky | QSPR as a means of predicting and understanding chemical and physical properties in terms of structure[END_REF][START_REF] Taskinen | Prediction of physico-chemical properties based on neural network modelling[END_REF]. In such methodology, computational chemistry may help to describe the molecular electronic structure and the decomposition reaction at quantum chemical level.

In this paper, we decided to explore the abilities of molecular modeling to predict the explosibility of nitroaromatic compounds (known as potentially explosive chemical substances and presenting complex decomposition channels [START_REF] Brill | Kinetics and mechanisms of thermal decomposition of nitroaromatic explosives[END_REF]). A series of geometric, electronic and energetic descriptors were computed using quantum chemical calculations. Correlations between these data and macroscopic properties related to explosibility were determined. In particular, explosibility is related to detonation and deflagration performances and to sensitivity to mechanical (impact, shock or friction) and thermal stresses or to electric discharges. Here, we focused on the thermal stability property and more precisely on the experimental decomposition enthalpy (or heat of decomposition). Experimental characterization is well defined, particularly by calorimetric measurements [START_REF] Grewer | Thermal Hazards of Chemical Reactions[END_REF] but, until now, only few approaches concern the standard symbols and phrases on the packaging labels and through safety data sheets (SDS). http://ecb.jrc.ec.europa.eu/classification-labelling/ prediction of the thermal stability properties based on QSPR models [START_REF] Yu | Prediction of the thermal decomposition property of polymers using quantum chemical descriptors[END_REF][START_REF] Keshavarz | Predicting activation energy of thermolysis of polynitro arenes through molecular structure[END_REF][START_REF] Keshavarz | Simple method for prediction of activation energies of the thermal decomposition of nitramines[END_REF]. In most of the cases, they are dedicated to specific (and small) classes of compounds (e.g. chromophores [START_REF] Figueiredo | Chemometric analysis of nonlinear optical chromophores structure and thermal stability[END_REF] or ionic liquids [START_REF] Kroon | Quantum chemical aided prediction of the thermal decomposition mechanisms and temperatures of ionic liquids[END_REF]). In the framework of chemical hazards management, Saraf and coworkers [START_REF] Saraf | Prediction of reactive hazards based on molecular structure[END_REF] proposed, to our knowledge, the only existing model concerning the heat of decomposition with an average absolute error of about 6%:

2 75 ) / ( NO n mol kcal H × = ∆ - (1)
where n NO2 is the number of nitro groups in the molecule.

Our work presents the first significantly correlated model for the prediction of this experimental property using an original approach combining QSPR methodology with quantum chemical calculations.

QSPR methodology

The quantitative structure-property relationship methodology (QSPR) consists in correlating quantitatively an experimental property with the molecular structures of the considered compounds. Thus, the relationship has the following general form between macroscopic and microscopic properties:

Property = f (Descriptors) (2)
An experimental data set provides the property values. Then different descriptors can be calculated to characterize the molecular structure of the compounds of the data set. The method can be based on neural networks [START_REF] Gasteiger | Neural networks in chemistry[END_REF], genetic algorithms [START_REF] Leardi | Genetic algorithms in chemometrics and chemistry: a review[END_REF] or statistical analyses such as multilinear regressions. We chose the last approach in this study.

The model obtained with such a regression has the following form:

n n X A X A X A A Y + + + + = ... 2 2 1 1 0 (3)
where Y is the property to predict, X i are the molecular descriptors and A i the corresponding regression constants. The reliability of this model is estimated with the coefficient R², which characterize the fitting of the calculated values with the experimental ones. To ensure the validation of the models, the use of external data is recommended. Nevertheless, in this study, the set of available experimental data contained only 22 molecules making its division into training and validation sets impossible. For this reason, the coefficient R² cv , using a cross-validated method, has been considered to characterize the predictivity of the model. [START_REF] Karelson | Molecular Descriptors in QSAR/QSPR[END_REF] is based on a leave-one-out procedure. For each molecule in the data set, a new multilinear regression is recalculated without changing the descriptors for all compounds in the data set, except this one. Then, the corresponding property value is predicted from this model equation. Finally, R² cv is calculated by correlating the obtained values for each compound with the experimental ones.

The cross-validation technique

Once a reliable equation is validated, this model can be used to predict the property, for other similar compounds to the selected ones, not yet measured and maybe not yet produced. The results can also help to understand the mechanisms leading to the studied property.

Data set selection

The choice of the training set of experimental data is a critical point in a QSPR analysis.

Experimental conditions may have a strong influence on the studied properties.

Therefore, all experimental values used in the fitting procedure have to be obtained in the same conditions. Differential Scanning Calorimetry (DSC) is often used to characterize the thermal stability of explosive compounds [START_REF] Frurip | Effective use of differential scanning calorimetry in reactive chemicals hazard evaluation[END_REF]. The experimental property studied in this work is the decomposition enthalpy taken from literature [START_REF] Duh | Chemical incompatibility of nitrocompounds[END_REF].

Our data set is composed of 22 nitroaromatic compounds, known as potentially explosive chemicals [START_REF] Gustin | Runaway reaction hazards in processing organic nitro compounds[END_REF]. The 22 molecules and their corresponding experimental decomposition enthalpy values are presented in Figure 1 and Table 1.

Descriptors

Different types of descriptors (geometric, electronic and energetic) have been selected to describe the 22 molecules. In Table 2 are presented the 14 descriptors calculated among the large number and diversity of those actually used in QSPR methodology [START_REF] Karelson | Quantum-chemical descriptors in QSAR/QSPR studies[END_REF].

Local descriptors

In nitro compounds, the carbon-nitrogen bond dissociation is currently considered as a rate-determining step of decomposition [START_REF] Saraf | Application of transition state theory for thermal stability prediction[END_REF]. For this reason, we decided to describe this carbon-nitrogen bond and its attached nitro group using geometric, electronic and energetic descriptors.

First, the length of the carbon-nitrogen bond d CN and the charge on the NO 2 functional group were calculated.

Secondly, the mid-point potential V mid , also used as a descriptor [START_REF] Murray | A relationship between impact sensitivity and the electrostatic potentials at the midpoints of C-NO 2 bonds in nitroaromatics[END_REF][START_REF] Rice | A quantum mechanical investigation of the relation between impact sensitivity and the charge distribution in energetic molecules[END_REF] in the field of energetic materials, was evaluated. This descriptor is an approximation of the electrostatic potential at the midpoint of the carbon-nitrogen bond [START_REF] Owens | Computational analysis of some properties associated with the nitro groups in polynitroaromatic molecules[END_REF].

CN N C mid d Q Q V 5 . 0 + = (4)
where Q C and Q N are respectively the atomic charges on carbon and nitrogen atoms.

The C-NO 2 dissociation bond energy (E diss ) has already been correlated to the impact sensitivity [START_REF] Politzer | Relationships between dissociation energies and electrostatic potentials of C-NO 2 bonds: applications to impact sensitivity[END_REF][START_REF] Rice | Density functional calculations of bond dissociation energies for NO 2 scission in some nitroaromatic molecules[END_REF] and considered for the estimation of decomposition temperatures [START_REF] Theerlynck | Towards improved models to rationalize and estimate the decomposition temperatures of nitroalkanes, nitramines and nitric esters[END_REF]. To calculate this energy, an homolytic dissociation of the R-NO 2 molecule is considered:

• • → - 2 2 NO + R NO R (5)
Finally, the corresponding dissociation energy E diss is calculated as the energy difference between products and reactants:

) ( ) ( ) ( 2 RNO E NO E + R E E 2 diss - = • • (6)
It is worth being noted that, for polynitroaromatic compounds, these descriptors were extracted for the weakest carbon-nitrogen bond in the molecules considering their dissociation energy.

Global descriptors

If the previous local descriptors as above defined are specific for nitro compounds, global descriptors present the advantage to be more suitable for extended data sets with compounds without any nitro group. Among such descriptors, the molecular weight M w was considered. It has also been used in the calculation of the oxygen balance as defined by Shanley et al. [START_REF] Shanley | The oxygen balance criteria for thermal hazards assessment[END_REF].

( ) [ ] w M Z Y X OB - + - = 2 2 1600 (7)
X, Y and Z are respectively the number of carbon, hydrogen and oxygen atoms in the molecule. This is a traditional empirical descriptor used in hazard prediction related to energetic materials [START_REF] Médard | Les explosifs occasionnels -Tomes 1 et 2[END_REF].

The atomization energy is the energy needed to break all bonds in the molecule. It can be calculated from the following equation:

) ( ) ( molecule E i atom E E natom i i n atomizatio - = ∑ (8)
The electronic structure can be described by the dipole moment (DM) and the mean polarizability (α). Dipole moment and polarizability characterize the influence of an external field on the electronic density. The dipole moment is a scalar entity whereas the polarizability is a tensor. The mean polarizability is calculated from the polarizability matrix:

( )

zz yy xx 3 1 α α α α + + = (9)
where α ii are the diagonal components of the polarizability matrix.

Conceptual density functional theory descriptors

Global electronic descriptors can also be estimated from conceptual density functional theory [START_REF] Geerlings | Conceptual density functional theory[END_REF][START_REF] Chermette | Chemical reactivity indexes in density functional theory[END_REF]. This methodology allows to redefine classical chemical reactivity concepts, e.g. electronegativity, in the framework of the density functional theory [START_REF] Pearson | Chemical Hardness[END_REF].

The ionisation potential (IP) and the electron affinity (EA) are calculated from the energies of the highest occupied and the lowest unoccupied molecular orbital ε HOMO and ε LUMO , according to the Koopmans theorem [START_REF] Koopmans | Über die Zuordnung von Wellenfunktionen und Eigenwerten zu den einzelnen Elektronen eines Atoms[END_REF].

HOMO IP ε - = ( 10
) LUMO EA ε - = (11)
The electronegativity (χ) characterizing the electron donating property of the system was identified to be the negative of the chemical potential (µ) [START_REF] Mulliken | A new electroaffinity scale; together with data on valence states and on valence ionization potentials and electron affinities[END_REF] and therefore:

µ ε ε χ - = + - = + = 2 ) ( 2 ) ( LUMO HOMO EA IP ( 12
)
The hardness (η) [START_REF] Parr | Absolute hardness: companion parameter to absolute electronegativity[END_REF] has been defined similarly to express the resistance of the system to the change in the number of electrons.

HOMO LUMO EA IP ε ε η - = - = (13)
These quantities (eq. 12 and 13) are already used as descriptors for different physicochemical properties such as heats of formation or boiling points [START_REF] Thanikaivelan | Application of quantum chemical descriptor in quantitative structure activity and structure property relationship[END_REF].

More recently, Parr et al. [START_REF] Parr | Electrophilicity index[END_REF] constructed the electrophilicity index (ω) to measure the loss in energy for a maximal electron flow from donor to acceptor. This index is proportional to the square of the chemical potential divided by the hardness.

η µ ω 2 ² = (14)
This descriptor has already been applied to the prediction of biological activity [START_REF] Parthasarathi | Electrophilicity Index as A Possible Descriptor of Biological Activity[END_REF].

Computational details

The molecular structures of 22 nitroaromatic compounds have been calculated. For each one, the structures of their phenyl radicals were also computed to access to the carbonnitrogen bond dissociation energy, the NO 2 radical being of course also calculated. All calculations presented in this work were performed using the Gaussian03 package [START_REF] Frisch | Gaussian 03[END_REF],

employing the density functional theory (DFT) with the parameter-free PBE0 [START_REF] Adamo | Toward reliable density functional methods without adjustable parameters: the PBE0 model[END_REF] hybrid functional. Geometry optimizations of stable species were performed with a 6-31+G(d,p) basis set to describe correctly the geometric structure. Vibrational frequencies were obtained at the same level of theory. The nature of the stationary points was checked by showing no imaginary frequency for our stable structures. The electronic structure of these molecules has been investigated using the natural population analysis (NPA) [START_REF] Reed | Natural population analysis[END_REF].

The Codessa software [54] has been used to obtain correlation coefficients and the QSPR model (using the integrated Best Multi Linear Regression analysis) for the prediction of the experimental decomposition enthalpy from the computed descriptors.

The significance of each descriptor in the equations was validated by performing a Student t-test validation at a 95% confidence level.

Results and Discussion

Linear correlations

The molecular descriptors, presented above, have been calculated for each molecule of the data set. These data, reported in Table 3, were analyzed to obtain a relationship between the molecular descriptors and the experimental explosibility.

Simple linear fittings were computed between each descriptor and the experimental enthalpy change. Correlations were appreciated from the coefficient R² (in Table 3) as shown for the descriptors that provide the best and lowest correlations (see Figure 2), with R² of 0.76 for molecular weight and less than 0.01 for mean polarizability respectively.

The local descriptors related to the C-NO 2 bond (presented in paragraph 4.1) do not exhibit any significant correlation with the decomposition enthalpy (R²<0.5). Therefore, a simple and direct breaking of the carbon-nitrogen bond appears to be not sufficient for the description of the decomposition enthalpy. More complex processes might occur [START_REF] Brill | Kinetics and mechanisms of thermal decomposition of nitroaromatic explosives[END_REF] and in particular other decomposition paths could exist as experimentally observed [START_REF] Brill | Influence of the substituent on major decomposition channels of the NO 2 group in para-substituted nitrobenzenes: a tandem mass spectrometric study[END_REF] and evidenced computationally [START_REF] Fayet | Theoretical study of the decomposition reactions in substituted nitrobenzenes[END_REF]. Furthermore, major decomposition paths involving interactions between the nitro group and ortho-substituents have been calculated [START_REF] Chen | A computational study on the kinetics and mechanisms for the unimolecular decomposition of o-nitrotoluene[END_REF][START_REF] Cohen | Mechanism of thermal unimolecular decomposition of TNT (2,4,6-trinitrotoluene): a DFT study[END_REF].

More substantial correlations are exhibited for the molecular weight (R²=0.76), the electrophilicity index (0.75), the electron affinity (0.71) and the oxygen balance (0.64).

It can be noted that these descriptors are auto-correlated. Indeed, the molecular weight and the electron affinity are used in the calculation of the oxygen balance and the electrophilicity index respectively.

Even if significant correlations appear, a single-descriptor approach is not sufficient for the prediction of the experimental property. For instance, considering the molecular weight or oxygen balance, the values of these descriptors do not vary with the position of the substituent. Indeed, ortho, meta and para nitrotoluenes present the same molecular weight (137 g/mol) whereas their decomposition enthalpies are different, i.e.

329, 284 and 318 kJ/mol, respectively. The use of at least one more descriptor is needed to characterize the influence of the substituent position in this case. Therefore, in a next step, multivariable regressions are investigated.

Multilinear model

All the descriptors previously studied have been integrated in a multivariable analysis using a multilinear regression. The best QSPR model is estimated as the most predictive in term of R² cv .

The obtained model is a six-parameter equation composed with the hardness, the electrophilicity index, the mean polarizability, the ionization potential, the dipole moment and the dissociation energy.

= = + - + - - + = ∆ - cv diss R R E DM IP H α ω η (15)
First we note that the six descriptors selected by this multilinear fitting procedure are different from the four global descriptors presenting the best linear correlations with the experimental decomposition enthalpy (see paragraph 6.1). The only exception is the electrophilicity index. The other parameters (Eq.15) are less linearly correlated with the experimental property.

Obviously, 22 molecules are not sufficient to obtain a robust predictive model. However, these results (with a R² of 0.91) are promising. For instance, when applied to the evaluation of the decomposition enthalpy of a molecule not included in the training set, as 2,4,6-trinitrotoluene (TNT), the predicted value is 920 kJ/mol which is close to its experimental value (998 kJ/mol [START_REF] Grewer | Thermal Hazards of Chemical Reactions[END_REF]). Work is in progress in this direction, but here we stress, once again the difficulty to obtain homogenous experimental data. Moreover, the model gives interesting indications for the future exploration of larger data sets. Indeed, the presence of descriptors arising from the so-called conceptual density functional theory (i.e. η, ω, IP) has to be noticed, all being related to the molecular reactivity. Hence, the characterization of the C-NO 2 bond, through E diss in Eq 15, and the molecular reactivity properties are important parameters for the prediction of decomposition properties like thermal stability. Furthermore, contrary to classical constitutional descriptors (for instance, the oxygen balance OB), the selected descriptors in the model are able to distinguish between isomers. Hence, a model based only on the number of nitro groups (like Eq. 1) is not sufficient to completely characterize thermal stability of nitroaromatic compounds.

Nevertheless this parameter is obviously fundamental. Indeed, the presence of chemical groups indicating explosive properties (e.g. nitro) in chemicals is a pre-evaluation element in chemical safety regulations for substances which may have explosive properties [START_REF]Recommendations on the Transport of Dangerous Goods: Manual of Tests and Criteria[END_REF]. Besides, mono-, di-and tri-nitroaromatic compounds can be clearly distinguished on figure 3, which represents the calculated values versus the experimental ones. So our model is consistent with this empirical consideration.

Conclusion

Nitroaromatic compounds have been modeled using the density functional theory. The aim was to observe correlation between the molecular structure of such compounds and an experimental property of explosibility, the decomposition enthalpy. The electronic and geometric structures were characterized with 14 molecular descriptors. The Sketchs of the considered nitroaromatic molecules (see table 1 for the exact nomenclature).

  Fig. 1.Sketchs of the considered nitroaromatic molecules (see table1for the

Fig. 2 .

 2 Fig. 2.Plot of (a) molecular weight and (b) mean polarizability with the experimental decomposition enthalpy.

Fig. 3 .

 3 Fig. 3.Calculated vs experimental decomposition enthalpies (kJ/mol).

Table 1 :

 1 Experimental decomposition enthalpies (-∆H exp /kJ/mol) of the 22 substituted

	nitrobenzene molecules from [30].
			-∆H exp
	1	nitrobenzene	339
	2	1,2-dinitrobenzene	518
	3	1,3-dinitrobenzene	586
	4	1,4-dinitrobenzene	622
	5	2-nitrotoluene	329
	6	3-nitrotoluene	284
	7	4-nitrotoluene	318
	8	2,6-dinitrotoluene	576
	9	3,4-dinitrotoluene	666
	10	2,4-dinitrotoluene	596
	11	2-nitroaniline	307
	12	3-nitroaniline	314
	13	4-nitroaniline	279
	14	2-nitrobenzoic acid	297
	15	3-nitrobenzoic acid	298
	16	4-nitrobenzoic acid	304
	17	2-nitrophenol	345
	18	3-nitrophenol	316
	19	4-nitrophenol	300
	20 1-chloro-4-nitrobenzene	360
	21	2,4-dinitrophenol	662
	22	2,4,6-trinitrophenol	1173

Table 2 :

 2 Descriptors.

	d CN	carbon nitrogen distance in Ǻ
	Q NO2	charge on the nitro group
	V mid	mid-point potential in Ǻ -1
	E diss	carbon nitrogen dissociation energy in a.u.
	M w	molecular weight in g/mol
	OB	oxygen balance in percents
	E atom	atomization energy in a.u.
	DM	dipole moment in D
	α	mean polarizability in Ǻ 3
	IP	ionization potential in a.u.
	EA	electron affinity in a.u.
	χ	electronegativity in a.u.
	η	hardness in a.u.

ω

electrophilicity index in a.u.

Table 3 :

 3 Descriptors calculated for 22 nitro compounds under study. See table 1 for acronyms and figure1for molecules nomenclature.Distances in Ǻ, dipole moments in D, oxygen balance in %, molecular weight in g/mol, polarizability in Ǻ 3 , Energies( ‡) in kcal/mol or(#) in u.a.

		d CN a	Q NO2	b	V mid a	E diss a ( ‡)	M w	OB	E atom a ( ‡)	DM a	α a	IP a(#)	EA a(#)	χ a(#)	η a(#)	ω a(#)
	1	1.467	-0.26	1.64	71.2	123	-163	1497.3	4.832	0.40	0.299	0.097	0.198	0.202	0.097
	2	1.467	-0.20	1.74	60.4	168	-95	1666.5	7.000	0.45	0.314	0.117	0.216	0.197	0.118
	3	1.470	-0.23	1.65	68.5	168	-95	1676.5	4.458	0.39	0.329	0.122	0.226	0.207	0.123
	4	1.472	-0.23	1.71	68.2	168	-95	1676.4	0.000	0.42	0.327	0.136	0.231	0.191	0.140
	5	1.446	-0.26	1.65	68.7	137	-181	1778.5	4.542	0.39	0.286	0.091	0.188	0.195	0.091
	6	1.467	-0.26	1.66	71.4	137	-181	1778.0	5.196	0.40	0.285	0.094	0.190	0.191	0.094
	7	1.462	-0.26	1.61	72.4	137	-181	1778.5	5.567	0.38	0.289	0.092	0.190	0.197	0.092
	8	1.468	-0.23	1.68	63.8	182	-114	1950.6	3.071	0.41	0.311	0.110	0.211	0.201	0.111
	9	1.468	-0.20	1.78	60.5	182	-114	1948.1	7.722	0.46	0.304	0.112	0.208	0.192	0.113
	10	1.469	-0.24	1.68	65.9	182	-114	1955.5	5.177	0.39	0.317	0.115	0.216	0.201	0.116
	11	1.441	-0.33	1.45	75.6	138	-151	1663.2	5.026	0.30	0.243	0.087	0.165	0.155	0.088
	12	1.468	-0.26	1.71	71.6	138	-151	1660.1	5.904	0.42	0.244	0.090	0.167	0.154	0.091
	13	1.448	-0.30	1.49	75.7	138	-151	1662.6	7.472	0.34	0.248	0.080	0.164	0.169	0.080
	14	1.466	-0.22	1.77	62.7	167	-120	1862.3	4.259	0.47	0.302	0.101	0.201	0.201	0.101
	15	1.468	-0.24	1.63	70.1	167	-120	1870.3	2.815	0.39	0.312	0.106	0.209	0.206	0.106
	16	1.470	-0.24	1.57	69.5	167	-120	1869.9	3.762	0.36	0.310	0.115	0.212	0.195	0.116
	17	1.460	-0.25	1.49	67.4	139	-132	1589.4	6.091	0.34	0.271	0.086	0.179	0.185	0.086
	18	1.469	-0.25	1.70	70.5	139	-132	1594.9	6.066	0.42	0.270	0.096	0.183	0.174	0.096
	19	1.455	-0.28	1.54	73.8	139	-132	1596.8	5.504	0.36	0.275	0.090	0.182	0.185	0.090

c correlation coefficient for the linear regression with the experimental decomposition enthalpies in table

1