Cleo Tebby
email: cleo.tebby@ineris.fr

Enrico Mombelli

Pascal Pandard

Alexandre R R Péry

Exploring an ecotoxicity database with the OECD (Q)SAR Toolbox and DRAGON descriptors in order to prioritise testing on algae, daphnids, and fish

Keywords: Cross-species comparison, Acute aquatic toxicity, Sensitivity, Quantitative Structure-Activity Relationships (QSAR), OECD (Q)SAR Toolbox

The European regulation on chemicals (REACh) places emphasis on reduction of systematic toxicity testing, thus fostering development of alternative methods. Consequently, we analysed acute toxicity data gathered by the Japanese Ministry of Environment for three species belonging to three different trophic levels (i.e., Pseudokirchneriella subcapitata 72-hour EC50, Daphnia magna 48-hour EC50 and Oryzias latipes 96-hour LC50). This paper investigates the relationships between the chemical structure and both the toxicity of the chemicals and the cross-species differences in sensitivity. The physicochemical properties of the chemicals were represented by the categories they belonged to in several widely-used categorisation schemes implemented by the freely available OECD (Q)SAR Toolbox and by quantitative molecular descriptors using DRAGON software. The outputs of these software products were analysed and compared in terms of quality of prediction and biological interpretation. Amongst the categorisations implemented by the OECD Toolbox, those focussing on bioaccumulation or biotransformation appeared to be the most interesting in terms of environmental prediction on a whole set of chemicals, in particular as the predicted biotransformation half-life is strongly dependent on hydrophobicity. In predicting toxicity towards each species, simple linear regression on logP performed better than PLS regression of toxicity on a very large set of molecular descriptors. However, the predictions based on the interspecies correlations performed better than the QSAR predictions. The results in terms of cross-species comparisons encourage the use of test strategies focussing on reducing the number of tests on fish.

Introduction

According to the European Community Regulation No 1272/2008 (CLP) on classification, labelling and packaging of substances and mixtures, substances are classified according to their acute aquatic toxicity and environmental fate. Acute aquatic toxicity data are usually determined using an algal species 72-hour EC50, a crustacean species 48-hour EC50, and a fish 96-hour LC50. Most of the time, these tests are carried out according to the Organisation for Economic Cooperation and Development (OECD) test guidelines (i.e. TG201 (OECD, 2006), TG202 (OECD, 2004), and TG203 (OECD, 1992)). Countries outside Europe have also adopted similar regulations: in Japan for instance, the Ministry of the Environment is responsible for evaluating toxicity of chemicals on ecosystems, and recommends the use of tests involving aquatic organisms such as Oryzias latipes (fishes) or Daphnia magna (daphnia), in addition to algae data [START_REF] Furuhama | Development of an ecotoxicity QSAR model for the KAshinhou Tool for Ecotoxicity (KATE) system, March 2009 version[END_REF]. The choice of these three trophic levels (primary producers, primary and secondary consumers) is considered to be relevant in order to protect aquatic ecosystems [START_REF] Wei | A new method for evaluating biological safety of environmental water with algae, daphnia and fish toxicity ranks[END_REF]. Furthermore, algae, daphnids, and fish are often used for cross-species comparisons of sensitivity [START_REF] De Roode | QSARs in ecotoxicological risk assessment[END_REF][START_REF] Dom | Aquatic multi-species acute toxicity of (chlorinated) anilines: Experimental versus predicted data[END_REF][START_REF] Henegar | What can be learnt from an ecotoxicity database in the framework of the REACh regulation?[END_REF][START_REF] Lessigiarska | QSAR investigation of a large data set for fish, algae and Daphnia toxicity[END_REF][START_REF] Tremolada | Quantitative inter-specific chemical activity relationships of pesticides in the aquatic environment[END_REF][START_REF] Weyers | Acute toxicity data: a comprehensive comparison of results of fish, daphnia, and algae tests with new substances notified in the European Union[END_REF].

The European regulation on chemicals (REACh) places emphasis on the reduction of systematic toxicity testing, thus fostering the development of alternative methods, such as testing strategies or statistical methods based on existing data. In this context, Quantitative Structure-Activity Relationships (QSAR) methods relate the physicochemical properties of chemicals with their toxicity on the basis that similar compounds have similar biological activities or properties [START_REF] Tropsha | Best Practices for QSAR Model Development, Validation, and Exploitation[END_REF].

Many QSAR models predict toxicity with specific descriptors such as the n-octanol/water partition coefficient (logP) alone [START_REF] Hsieh | Quantitative structure-activity relationships for toxicity of nonpolar narcotic chemicals to Pseudokirchneriella subcapitata[END_REF], associated with a chemical class as in EcoSAR software, or, depending on the mechanism of toxic action [START_REF] Duchowicz | Quantitative Structure-Toxicity Models for Heterogeneous Aliphatic Compounds[END_REF], associated with other parameters such as the energy of the lowest unoccupied molecular orbital [START_REF] Chen | Toxicity of substituted anilines to Pseudokirchneriella subcapitata and quantitative structure-activity relationship analysis for polar narcotics[END_REF][START_REF] Huang | Toxicity and quantitative structure-activity relationships of nitriles based on Pseudokirchneriella subcapitata[END_REF], or the dissociation constant, Hammett Co constant, index of valence molecular connectivity, perimeter of the efficient cross-section of molecule, and melting point [START_REF] Shigeoka | Acute toxicity of chlorophenols to green algae, Selenastrum capricornutum and Chlorella vulgaris, and quantitative structure-activity relationships[END_REF]. Other QSAR models rely on generating large sets of descriptors and using statistical methods to reduce dimensionality and identify underlying structural factors influencing toxicity. QSAR models have been made available in various software packages [START_REF] De Roode | QSARs in ecotoxicological risk assessment[END_REF].

Structure-Activity Relationships (SAR) can take the physicochemical properties into account by considering the functional or chemical class of molecules. Besides considering the overall level of toxicity of a chemical towards each species, studying species' relative sensitivity is particularly relevant in the regulatory context of environmental risk assessment and when considering the avoidance of tests on vertebrates [START_REF] Hoekzema | Alternative approaches can greatly reduce the number of fish used for acute toxicity testing[END_REF][START_REF] Hutchinson | A strategy to reduce the numbers of fish used in acute ecotoxicity testing of pharmaceuticals[END_REF][START_REF] Jeram | A strategy to reduce the use of fish in acute ecotoxicity testing of new chemical substances notiWed in the European Union[END_REF]. [START_REF] Tremolada | Quantitative inter-specific chemical activity relationships of pesticides in the aquatic environment[END_REF] developed Quantitative Inter-specific Chemical Activity Relationships (QICAR) for pesticides on algae, daphnids, and fish. The toxicity was related to the functional and chemical class of the compounds and the physicochemical properties also contributed to predicting the toxicity towards one species based on the toxicity towards a different species.

In this study, we analysed an ecotoxicity database produced by the Japanese Ministry of Environment and available online (Japanese Ministry of Environment, 2010).

Additional data on physicochemical properties of these chemicals were collected in two different ways. The chemicals were classified according to several widely-used classifications implemented by the OECD (Q)SAR Toolbox. Quantitative physicochemical properties were obtained by calculating a set of molecular descriptors using DRAGON software and the results were compared with those obtained using the OECD Toolbox. Since the latter is freely available for download and a free version of the DRAGON descriptors is available online [START_REF] Tetko | Virtual computational chemistry laboratorydesign and description[END_REF]VCCLAB. Virtual Computational Chemistry Laboratory, 2005) our methodology can be easily adopted and reproduced by all the stakeholders involved in regulations on chemicals.

The aims of this paper are two-fold. Firstly, it studies structure-activity relationships between the acute toxicity of the chemicals for three aquatic species and the physicochemical properties. Secondly, it relates the cross-species differences in sensitivity with the physicochemical properties.

Material and methods

Data Set

The data in the ecotoxicity database were gathered by the Japanese Ministry of Environment and available online (Japanese Ministry of Environment, 2010). The tests were conducted according to the OECD test guidelines performed under Good Laboratory Practices (GLP). We selected the toxicity data of organic chemicals that had been collected in tests in or after 2002. Data from tests using vehicles to dissolve the test substance were discarded.

We studied the acute toxicity expressed as the EC50 or LC50 for three aquatic species from three different trophic levels, one alga, one daphnid and one fish. The species selected were Pseudokirchneriella subcapitata, previously known as Selenastrum capricornutum (72-hour EC50, calculated from average specific growth rate), Daphnia magna (48-hour EC50), and Ozyrias latipes (96-hour LC50) also known as the Japanese Medaka. P. subcapitata and D. magna are two widely used and relatively sensitive species; the sensitivity of O. latipes is comparable to that of other fish species [START_REF] Wei | A new method for evaluating biological safety of environmental water with algae, daphnia and fish toxicity ranks[END_REF]. The toxicity information for all three species was available for 85 chemicals in the database. The Chemical Abstract Service (CAS) registration numbers and names of these chemicals are reported in Table 1 of the supplementary material. The final dataset is chemically heterogeneous, containing aliphatic and aromatic organic chemicals, and excluding mixtures. At least eight of the chemicals are pesticides . There were no highly hydrophobic chemicals, as the MlogP (Moriguchi's octanol-water partition coefficient obtained by multiple regression analysis on 13 structural parameters [START_REF] Moriguchi | Simple method of calculating Octanol/Water Partition Coefficient[END_REF]) ranged from -1.4 to 4.5 and the AlogP (Ghose-Crippen octanol-water partition coefficient based on hydrophobic atomic constants measuring the lipophilic contribution of atoms in the molecule [START_REF] Viswanadhan | Atomic physicochemical parameters for 3 dimensional structure directed quantitative structure -activity relationships .4. Additional parameters for hydrophobic and dispersive interactions and their application for an automated superposition of certain naturally-occurring nucleoside antibiotics[END_REF]) ranged from -0.9 to 5.8.

The two-dimensional structure of the chemicals was coded in SMILES notation [START_REF] Weiniger | SMILES, a chemical language and information system. 1. Introduction to methodology and encoding rules[END_REF][START_REF] Weiniger | SMILES. 2. Algorithm for Generation of Unique SMILES Notation[END_REF] in the database. This was used as an input in the freely available OECD Q(SAR) Toolbox version 2.0 (OECD, 2010) and used for profiling the chemicals according to several predefined classification schemes implemented by the Toolbox. The chemicals were thus profiled according to their organic functional groups, the EcoSAR Classification, the mutagenicity/Carcinogenicity alerts by Benigni/Bossa, the US EPA Categorization, the Lipinski Rule, the OASIS Acute Toxicity Mode Of Action, the OASIS protein binding, the superfragment profiling, the OASIS DNA Binding, the Estrogen Receptor binding, the bioaccumulation metabolism alerts (BioWIN MITI fragments, developed to predict biodegradation by microorganisms [START_REF] Tunkel | Predicting ready biodegradability in the Japanese ministry of international trade and industry test[END_REF], with some additional classes, US EPA), and the primary biotransformation half-lives in fish (US EPA), developed partly to improve design of bioaccumulation and toxicity testing [START_REF] Arnot | A quantitative structureactivity relationship for predicting metabolic biotransformation rates for organic chemicals in fish[END_REF]. These categorisations were used in order to analyse the relationships between the presence of particular structural fragments and the toxicity of chemicals. Attention is drawn to the fact that, even inside a given categorisation scheme, a chemical can belong to several classes. In view of further analysis, the data extracted from the Toolbox was thus recoded into a table where each column was a binary variable indicating whether a given chemical belonged or not to a given class (one column for each chemical class).

The three-dimensional models of the molecules were built with hydrogens and their stereochemistry coincides with the chirality information displayed by the website ChemIDplus lite (http://chem.sis.nlm.nih.gov/chemidplus/chemidlite.jsp). Molecules with less than 6 rotatable bonds were submitted to systematic conformational analysis otherwise they were submitted to a conformational search by Boltzman jump (50000 steps, T=300 K, RMSD=60 degrees). At the end of this search, a conjugate gradients minimization was performed with the MM2 force field and the MOPAC/AM1 semiempirical quantum mechanical calculations were used to generate atom partial charges and further optimize molecular conformations. All the described procedures were carried out in vacuum thanks to the software VegaZZ [START_REF] Pedretti | VEGA -An open platform to develop chemo-bio-informatics applications, using plug-in architecture and script programming[END_REF]. A set of 797 two-and three-dimensional molecular descriptors were calculated using the software DRAGON version 5.3 licensed by Talete srl, Milan (Italy).

Statistical analysis

The analyses aimed to answer three main questions:

-In the database, are there overall cross-species differences in sensitivity? -How are classes defined according to chemical structure related to the toxicity of chemicals and to the cross-species differences? -When considering a large set of DRAGON descriptors, how are these related to the categorisations implemented by the OECD (Q)SAR Toolbox, to toxicity, and to cross-species differences in sensitivity?

Data set verification

Basic properties of the EC50 and LC50 distributions were investigated to ensure the relevance of our statistical analysis. Distribution of log(EC50) or log(LC50) of each species was checked for large gaps (exceeding 15% of the total range) between two consecutive values [START_REF] Tropsha | Best Practices for QSAR Model Development, Validation, and Exploitation[END_REF]. 15% is a relatively high threshold but accounts for the small size of the dataset.

The EC50 or LC50 values were compared with the estimated solubility (logS) calculated with software ALOGPS version 2.1 [START_REF] Tetko | Estimation of aqueous solubility of chemical compounds using E-state indices[END_REF]VCCLAB, 2005).

Cross-species comparison

The interspecies toxicity correlations were measured by the square of Pearson's correlation coefficient between the log(EC50) or log(LC50) of all pairs of species.

The most and the least sensitive species towards each chemical were identified. χ² tests on the counts of most or least sensitive species were used to determine whether overall sensitivity was significantly different between species. The proportion of chemicals which were more toxic towards P. subcapitata or D. magna than towards O. latipes was also computed.

Relationships between chemical class and toxicity

Either Student's t-test or the Mann-Whitney test on the log(EC50) or log(LC50) was used to determine whether the toxicity of a chemical class was significantly different from that of the rest of the chemicals. Student's t-test was used when the log(EC50) or log(LC50) was distributed normally in each population; the Mann-Whitney test was used if this were not the case. Normality in subpopulations was assessed by the Shapiro-Wilk test.

The chemical class was cross-tabulated against the most sensitive species and against the least sensitive species, thus producing one contingency table for each chemical class. Fisher's exact test was used to determine whether sensitivity to a given chemical class was significantly different between trophic levels, taking into account the possible overall cross-species differences in sensitivity. This test was repeated for each chemical class. The expected counts under the null hypothesis of no association between the chemical classes and the species' relative sensitivity were compared with the observed counts in order to identify which species (either one or two) were the most sensitive.

Relationships between molecular descriptors and toxicity

The molecular descriptors were first of all explored by Principal Component Analysis (PCA) in order to identify the main structural differences between the chemicals in the dataset that are highlighted by the set of descriptors, and to relate these with the categorisations implemented by the OECD (Q)SAR Toolbox (with Mann-Whitney's test), with the toxicity (test of nullity of Pearson's correlation coefficient), and the species' relative sensitivity (Fisher's exact test).

Either χ² tests or Fisher's exact tests were used to identify the chemical classes that were characteristic of specific subgroups revealed by the PCA. Fisher's exact test was used when the χ² approximations could have been incorrect.

The database was randomly split into a training set and a test set (for external validation) with a three to one size ratio (63 chemicals in the training set and 22 chemicals in the test set, Table 1 of the supplementary material), using random sampling without replacement. The distribution of log(EC50) or log(LC50) in both sets was then compared and a Mann-Whitney test was performed for each endpoint to check that the values in the test set were not significantly different from the ones in the training set. Partial Least Squares (PLS) regression was performed on the training set on the whole set of molecular descriptors. Prediction quality was evaluated with the cross-validated R² with 7 segments and with the external validation Q² ext calculated using the average log(EC50) or log(LC50) on the test set [START_REF] Schüürmann | External Validation and Prediction Employing the Predictive Squared Correlation Coefficient -Test Set Activity Mean vs Training Set Activity Mean[END_REF]. Linear regression of toxicity on two different calculations of the logP (MlogP and AlogP) and their squares was performed. The quality was evaluated by calculating the coefficient of determination (R²), both on the training and on the test set (Q² ext ).

All calculations and graphs were made with R software version 2.12.0 (R Development Core Team, 2010), using contributed packages FactoMineR and pls.

Results

Data set verification

72-hour EC50 for P. subcapitata ranged from 0.0064 mg/L to 82 mg/L. 48-hour EC50 for D. magna ranged from 0.016 mg/L to 87 mg/L. 96-hour LC50 for O. latipes ranged from 0.031 mg/L to 90 mg/L. All the chemicals could thus be ranked in one of the classes "very toxic", "toxic", or "harmful" defined for the classification and labelling of chemicals according to the European Community Regulation No 1272/2008 (CLP). Subsequent calculations were performed on the logarithm of EC50 or LC50.

For all three species, there were no gaps in log(EC50) or log(LC50) greater than 15% of the total log(EC50) and log(LC50) ranged.

The log(EC50) and log(LC50) expressed in mol/L were compared with the estimated logS for 81 of the chemicals (the logS could not be estimated for the four remaining chemicals). The calculated EC50 or LC50 values were greater than the estimated solubility for two chemicals for P. subcapitata (0.43 log units in the worst case) and for two chemicals for O.latipes (0.12 log units in the worst case).

According to these comparisons, we found no reason to discard any of the substances due to solubility aspects.

General cross-species comparisons

The interspecies toxicity correlation was highest between D. magna and O. latipes (R²=68.3%). The correlation was smaller between the toxicity towards P. subcapitata and D. magna (R²=53.8%) and between the toxicity towards P. subcapitata and O. latipes (R²=41.3%) (Figure 1).

Insert Figure 1

For each chemical, the EC50 or LC50 were compared between the three species in order to identify the most and least sensitive species. There was one tie in the most sensitive species and two ties in the least sensitive species; in those cases the calculations were done for all possible combinations (two combinations for the test on the most sensitive species and four combinations for the test on the least sensitive species).

All three species were about equally often the most sensitive species (χ² =2.85 or 2.14, p-value=0.24 or 0.34, depending on the assignment in case of ties). D. magna was rarely the least sensitive species (to 14 or 15 chemicals) (χ² =10.54 to 13.01, p-value=0.0015 to 0.0051, depending on the assignment in case of ties). Toxicity towards P. subcapitata or D. magna was higher than towards O. latipes for 75.3% of the chemicals.

Chemical class and toxicity

The OECD (Q)SAR Toolbox identified a total of 243 classes. Only classes that comprised at least 10% of the whole dataset were selected: as there were 85 chemicals in the dataset, classes with 9 to 76 chemicals were selected for further analysis. Some categorisations did not yield classes with sufficient chemicals: the Lipinski rule and the superfragment profiling were not utilised in this study. In the dataset, 80 out the 85 chemicals were bioavailable according to the Lipinski rule and 77 had no superfragment. Some of the classes were identical in several categorisations, such as phenols or anilines (with some differences in the US EPA classification).

Degree of toxicity according to the chemical class

Few chemical classes show a difference in toxicity towards any of the species considered independently. The ratios of EC50 or LC50 in each class on EC50 or LC50 out of the class are reported for each species in the columns entitled "Toxicity of chemical class" in Table 1. Ratios greater than 1 indicate that the class is less toxic towards a given species than the other chemicals in the data set. The ratios also allow cross-species comparisons to a certain extent, as one class of chemicals can be more toxic than the other classes towards one particular species.

Insert Table 1

Only the categorisations regarding bioaccumulation and biotransformation yielded classes of chemicals that were both represented by a sufficient number of chemicals to be tested and were significantly more or less toxic than the other chemicals in the dataset (i.e. presented a significant difference in log(EC50) or log(LC50)). The chemicals with methyl substituents on an aromatic ring (identified in the bioaccumulation metabolism alerts) were less toxic towards all species than other chemicals (EC50 or LC50 2.7 to 6.0 times greater). Chemicals with an unsubstituted phenyl group were less toxic towards D. magna and O. latipes. The chemicals with several fused 6C aromatic rings (bioaccumulation metabolism alert) were more toxic towards P subcapitata (EC50 11 times smaller) and D. magna (EC50 7.1 times smaller) than other chemicals. Seven out of nine of these chemicals were naphthalene derivatives. The biotransformation half-lives were undetermined for two chemicals, and were moderate (1 to 10 days), fast (0.1 to 1 day) or very fast (under 0.1 day) for all other chemicals. Chemicals with moderate biotransformation half-lives (i.e. longer half-lives than the other chemicals) were more toxic towards all species (EC50 or LC50 3.4 to 4.7 times smaller) and chemicals with very short half-lives were less toxic towards D. magna and O. latipes (EC50 or LC50 2.8 times greater). Closer inspection of the 19 chemicals with very short biotransformation half-lives revealed that 8 of them belonged to the aniline class (in the EcoSAR classification), and none of them were neutral organics (in the EcoSAR classification) or basesurface narcotics (OASIS Acute Toxicity Mode Of Action).

A more specific toxicity of chemicals with methyl groups towards O. latipes was also shown. Also, chemicals with an alkyl substituent on an aromatic ring were less toxic than other chemicals towards P. subcapitata.

Cross-species differences in sensitivity

Each chemical class was cross-tabulated separately against the most sensitive species and against the least sensitive species, thus producing one contingency table for each chemical class. Fisher's exact test was applied to each contingency table to determine whether sensitivity to a given chemical class was significantly different between the species. To account for ties, the test was carried out on all combinations for each chemical class. The number of tests with a p-value smaller than 0.05 are reported in Table1 in two separate columns for the test on the most sensitive species and on the least sensitive species. When p<0.05, the observed counts were compared with the counts expected under the null hypothesis of no associations between the structural categories and ecotoxicological sensitivity in order to determine which species was (or were) significantly the most sensitive. The significantly most sensitive species to each chemical class are reported in Table 1 in the third to last and in the last column.

When considering the most sensitive species, there was a difference in sensitivity between the three species for 7 out of the 40 chemical classes investigated. O. latipes was the most sensitive species towards chemicals with an impaired OH or NH 2 group (Estrogen Receptor binding) or a methyl metabolism alert (BioWIN MITI fragment) (most sensitive species to 41% or 45% of the chemicals with methyl alerts and most sensitive to only 15% of the chemicals without methyl alerts). D. magna was the most sensitive species towards arenes (most sensitive species to 45% of the arenes and most sensitive to only 12.5% or 17% of the non-arenes). In the organic functional group classification, arenes are chemicals with one to three six-membered aromatic rings bearing diverse substituents. D. magna was also the least sensitive towards chemicals with a methylene functional group (most sensitive species to only 14% of the chemicals with methylene and most sensitive to 44% or 46% of the chemicals without methylene). P. subcapitata was less sensitive than the other species towards chemicals with an alkyl substituent on an aromatic ring (BioWIN MITI fragment) (most sensitive species to only 11% of the chemicals with an alkyl substituent on an aromatic ring and most sensitive to 45% of the other chemicals).

When considering the least sensitive species, there was a difference in sensitivity between the three species for phenols, which were identified as a class by four of the categorisations investigated. P. subcapitata was shown to be less sensitive towards phenols. Because of ties in EC50 or LC50, it was unclear whether P. subcapitata was less sensitive towards chemicals with a methyl functional group and whether O. latipes was less sensitive to anilines (in the US EPA classification) and more sensitive towards chemicals with a carbon with four single bonds and no hydrogen (BioWIN MITI fragment).

The chemical classes for which the toxicity towards P. subcapitata or D. magna was notably higher than towards O. latipes in over 85% of cases were identified. This was observed only for a small number of chemical classes: chemicals with aryl halide organic functional groups (aromatic chloride BioWIN MITI fragments for example), and chemicals with several fused 6-carbon aromatic rings (BioWIN MITI fragments).

Molecular descriptors and toxicity

Exploratory analysis

Amongst the set of 797 molecular descriptors, the ten that correlated most with the log(EC50) or log(LC50) of each species are shown in Table 2. The signification of the descriptors is reported in Table 2 of the supplementary material. The descriptor that correlated best with the toxicity was the squared logP (AlogP² with P. subcapitata: r=-0.43; MlogP² with D. magna: r=-0.57; MlogP² with O. latipes: r=-0.53). AlogP was amongst the top ten for P. subcapitata and both MlogP and AlogP were amongst the top ten for D. magna and O. latipes. As expected, toxicity of chemicals increased with hydrophobicity.

Insert Table 2

Principal Component Analysis (PCA) was applied to the 797 standardised molecular descriptors in order to identify the main structural differences between the chemicals in the dataset that are highlighted by the set of descriptors, without focussing on toxicity. Pearson's correlation coefficient was calculated between the four first principal components and the two different calculations of the logP and their squared values that were included in the molecular descriptors. The correlations between MlogP, MlogP², AlogP, and AlogP² and the first and fourth principal components were all significantly different from zero (p<0.05) and relatively high with the fourth component (Pearson's r ranged from 0.34 to 0.52 with the first component and ranged from 0.50 to 0.56 for the fourth component). This shows that hydrophobicity is one of the main characteristics reflected by the DRAGON molecular descriptors and that differentiates the chemicals in this database.

The coordinates of the chemicals of each class on the four first principal components were compared with the coordinates of chemicals that were not in the class (Table 2 of the supplementary material).

Most of the classes of chemicals were represented. The two first principal components discriminated chemicals with methyl functional groups (Figure 2). The chemicals in this class were more hydrophobic (Student's t test, p<0.05 for MlogP and AlogP). This could be related to the fact that methyl groups are small hydrophobic units. The second principal component strongly discriminated chemicals with arene functional groups (aromatic hydrocarbons) from non-cyclic structures (Figure 2).

Insert Figure 2

The differences in biotransformation half-life were represented on the fourth component, therefore the relationship between logP and biotransformation half-life was investigated. There was a significant relationship between logP and biotransformation half-life (one-factor analysis of variance without the data with an undefined half-life; p-value= 4.705e-11): the chemicals with a very short half-life were less hydrophobic (Figure 3).

Insert Figure 3

Some chemical classes that were not represented on the first five principal components were actually pointed out earlier in this paper as being important regarding overall toxicity (chemicals with a methyl substituent on an aromatic ring) or relative sensitivity of species (chemicals with an impaired OH or NH 2 group). A number were not represented on the first three. Also there were no significant correlations between log(EC50) or log(LC50) and the two first PCA components (test of nullity of Pearson's correlation coefficient, p<0.05). The 2D and 3D DRAGON molecular descriptors used were thus representative of several structural and chemical properties that are related to toxicity, but not of all the structural fragments pointed out earlier in this paper.

The most and the least sensitive species to each chemical are shown on the scatter plot in Figure 4. In case of ties, the most/least sensitive species to a chemical was randomly assigned. The patterns tended to show that relative toxicity towards the various species were at least related to the two first principal components. Approximately one quarter of the scatter plot contained 30 chemicals that were mostly more toxic towards D. magna or O. latipes and less toxic towards P. subcapitata or O. latipes. This area (coordinate on 1 st component smaller than 1 and coordinate on 2 nd component smaller than -1.85) was examined by a series of Fisher's exact tests on the chemical class (results in Table 4 of the supplementary material). The chemicals in this area were all arenes (organic functional group classification) but displaying a wide variety of substituents; compared to the rest of the dataset, a large proportion had an aryl halide functional group with chlorine or were classified as anilines in the US EPA classification. None of them had a methylene functional group. The chemicals in this area were on average as hydrophobic as the other chemicals (Mann-Whitney test on MlogP, p=0.051), however the range and variance of logP were smaller (MlogP ranged from 1.2 to 3.8; Levene's Test for Homogeneity of Variance on MlogP: p=0.0040). The chemicals in this area had noticeably lower values for some DRAGON descriptors belonging to the groups of Edge Adjacency indices Eigenvalues, of lowest eigenvalues of the Burden matrix, 3rd component size or accessibility directional WHIM index as well as a lower spherosity (SPH descriptor values were characteristic of aromatic chemicals).

Insert Figure 4

The first principal components are thus related to cross-species differences in sensitivity as well as to some chemical classes which in turn are related to the degree of toxicity of the chemicals. Therefore the molecular descriptors might be predictors of the acute toxicity of the chemicals. The relationships between DRAGON descriptors on one side and toxicity and relative sensitivity on the other side were further investigated by PLS regression during QSAR modelling. Regarding external validation quality of prediction, PLS modelling did not perform as well as a simple linear regression of log(EC50) or log(LC50) on either calculations of the logP or their squares (R², R² adjusted for the number of variables, and Q² are reported in Table 3), where the Q² ext ranged from 0.12 (MlogP, P. subcapitata) to 0.61 (AlogP², D. magna) although the R² values obtained on the training set were smaller. These predictions were also of poor predictive validity. Multiple linear regression on both the logP and logP² [START_REF] Patrick | Quantitative Structure-Activity Relationships[END_REF] did not perform better regarding adjusted R² or Q² ext (results not shown).

QSAR models of toxicity

Insert Table 3

Compared to linear regression on AlogP², the predictions obtained considering pairwise interspecies correlations (Table 3) were better for all three species. For P. subcapitata, Q² ext with linear regression on AlogP² was equal to 0.47, whereas on log(EC50) D.magna , Q² ext was equal 0.55. For D. magna, the increase in Q² ext was greater: on AlogP², Q² ext was equal to 0.61 whereas on log(LC50) O. latipes , Q² ext was equal to 0.80. For O. latipes, on AlogP², Q² ext was equal to 0.57; on log(EC50) D. magna , Q² ext reached 0.77. Predicting toxicity towards one species with toxicity from both other species (instead of only one) only improved the predictions for D. magna when considering the adjusted R² (R² adj =0.64 on log(LC50) O. latipes ; R² adj =0.72 on both log(EC50) P subcapitata and log(LC50) O. latipes ). Adding the AlogP² measure of hydrophobicity to the later model did not further improve the predictions when considering the adjusted R² (Table 3).

PLS discriminant analysis of the most sensitive species according to the physico-chemical descriptors performed poorly (results not shown).

Discussion

The dataset analysed was chemically heterogeneous, containing aliphatic and aromatic organic chemicals, but there were no highly hydrophobic chemicals (with AlogP>6) and no persistent chemicals. Its representativeness is quite limited, considering the physicochemical properties of chemicals notified in Europe.

The interspecies correlations between fish and Daphnia toxicity values were higher than the correlation with algae, thus confirming previous research on other ecotoxicity databases [START_REF] Henegar | What can be learnt from an ecotoxicity database in the framework of the REACh regulation?[END_REF][START_REF] Lessigiarska | QSAR investigation of a large data set for fish, algae and Daphnia toxicity[END_REF][START_REF] Tremolada | Quantitative inter-specific chemical activity relationships of pesticides in the aquatic environment[END_REF][START_REF] Weyers | Acute toxicity data: a comprehensive comparison of results of fish, daphnia, and algae tests with new substances notified in the European Union[END_REF][START_REF] Zhang | Interspecies correlations of toxicity to eight aquatic organisms: Theoretical considerations[END_REF]. [START_REF] Zhang | Interspecies correlations of toxicity to eight aquatic organisms: Theoretical considerations[END_REF] suggest the difference is due to bio-uptake and mode of action. The magnitude of the correlations are similar to those calculated by [START_REF] Tremolada | Quantitative inter-specific chemical activity relationships of pesticides in the aquatic environment[END_REF] for pesticides between Daphnia and trout (R²=59%), and for unspecific pesticides between algae and Daphnia (R²=52%) and algae and trout (R²=49%).

In an analysis of 164 chemicals from a database provided by the French Ministry of Ecology and Sustainable Development, [START_REF] Henegar | What can be learnt from an ecotoxicity database in the framework of the REACh regulation?[END_REF] showed that fish were less sensitive than algae and Daphnia. [START_REF] Weyers | Acute toxicity data: a comprehensive comparison of results of fish, daphnia, and algae tests with new substances notified in the European Union[END_REF] showed that alga was the most sensitive trophic level in a database of 694 chemicals, extracted from the New Chemicals Database which comprises information about new chemical substances manufactured or imported in the European Union since 1981. In a study on 1,439 chemicals, [START_REF] Jeram | A strategy to reduce the use of fish in acute ecotoxicity testing of new chemical substances notiWed in the European Union[END_REF] showed that alga was the most sensitive and fish was the least sensitive.

However, their database included many herbicides, which are likely to target algae. The database analysed in this study was smaller and contained only eight pesticides, of which one herbicide, one insecticide, and several non-specific pesticides. Part of our results suggests that daphnids might be the most sensitive in general, but the small size of the dataset leads to limit the significance of this conclusion. Moreover, the other published studies do not support this conclusion.

OECD, based on the several studies [START_REF] Hoekzema | Alternative approaches can greatly reduce the number of fish used for acute toxicity testing[END_REF][START_REF] Hutchinson | A strategy to reduce the numbers of fish used in acute ecotoxicity testing of pharmaceuticals[END_REF][START_REF] Jeram | A strategy to reduce the use of fish in acute ecotoxicity testing of new chemical substances notiWed in the European Union[END_REF], proposed a threshold approach to reduce the number of fish for the determination of acute aquatic toxicity in a regulatory context. The authors report that algae or daphnia were more sensitive than fish for 80% to 88% of the chemicals in their databases, compared to 75% in our dataset. The testing strategy based on that developed by [START_REF] Hutchinson | A strategy to reduce the numbers of fish used in acute ecotoxicity testing of pharmaceuticals[END_REF] would however still have enabled reducing the number of fish used, by designing the fish tests on the basis of the results obtained with algae or daphnia and proceeding sequentially (limit test in a first step and full test if lethal effects are observed at the lowest EC50 determined in the daphnia and algae tests).

When using the OECD (Q)SAR Toolbox to identify structural fragments related to toxicity, only the categorisations focussing on bioaccumulation or biotransformation yielded classes of chemicals that were both represented by a sufficient number of chemicals to be tested and were significantly more or less toxic than the other chemicals in the dataset. In particular, chemicals with very short biotransformation half-lives in fish were less toxic to all species. In this data set, the chemicals with a very short biotransformation half-life were less hydrophobic, as the QSAR model which was used to predict the biotransformation half-life [START_REF] Arnot | A quantitative structureactivity relationship for predicting metabolic biotransformation rates for organic chemicals in fish[END_REF] was partly based on the logP. Linear regression of log(EC50) or log(LC50) on logP confirmed previous research showing that toxicity of chemicals increased with hydrophobicity [START_REF] Duchowicz | Quantitative Structure-Toxicity Models for Heterogeneous Aliphatic Compounds[END_REF][START_REF] Hsieh | Quantitative structure-activity relationships for toxicity of nonpolar narcotic chemicals to Pseudokirchneriella subcapitata[END_REF][START_REF] Shigeoka | Acute toxicity of chlorophenols to green algae, Selenastrum capricornutum and Chlorella vulgaris, and quantitative structure-activity relationships[END_REF]. It is unclear whether the difference in toxicity observed in the dataset is due to biotransformation half-life or to hydrophobicity.

Linear regression on a measure of logP performed better than PLS regression on a large set of DRAGON molecular descriptors. According to [START_REF] Golbraikh | Beware of q²![END_REF], two of the necessary conditions for a model to be considered as predictive are that R²>0.6 and Q²>0.5. None of the QSAR models developed in this paper by PLS regression or by linear regression on a measure of logP satisfy these criteria. The results from the PLS regression do not encourage the selection of any particular molecular descriptors apart from the logP. Interspecies correlations between toxicity tests were better than between toxicity and logP, as shown by [START_REF] Weyers | Acute toxicity data: a comprehensive comparison of results of fish, daphnia, and algae tests with new substances notified in the European Union[END_REF] in a study on part of the New Chemicals Database. The external validation Q² values were often greater than the R² values and therefore the distribution of the chemical classes in the original data set in the training set and test set was checked. In the test set, there were a slightly higher proportion of aromatic chlorides and chemicals with moderate biotransformation half-lives, but this did not result in any significant difference in toxicity.

The conclusions in terms of relative sensitivity of species are not straightforward, as the results obtained with the test on the most sensitive species are not confirmed by the test on the least sensitive species. On the other hand, tests generally focus on the most sensitive species rather than on the least sensitive, meaning that the data concerning the most sensitive species are more reliable. Basing their research on the Verhaar classification, Vaal et al. (1997b) concluded that chemicals with an unspecific mode of action were the least toxic and demonstrated smaller differences in interspecies sensitivity.

They also concluded (1997a) that, except for anilines, compounds with the highest overall toxicity also have the largest variation in toxicity for different species. In our dataset however, the most toxic chemical classes (i.e. chemicals with the longest biotransformation half-lives) did not display any cross-species differences in toxicity.

Our analysis with the OECD Toolbox outputs was repeated on the same dataset with data collected on 141 additional chemicals before 2002, which might be less reliable as in some cases a dispersant was used, thus increasing almost three-fold the number of chemicals. The results obtained on the larger database confirmed the differences in sensitivity towards arenes (D. magna being more sensitive than the other species), methylene functional groups (D. magna being less sensitive), and phenols (P. subcapitata being less sensitive). Some of the classes identified by the Toolbox, such as anilines (aromatic amines), are subject to differences between categorisation schemes (US EPA versus organic functional groups or EcoSAR). On the other hand, arenes appear to be easily identifiable and not subject to differences between categorisation schemes. It has to be noted that arenes are a rather heterogeneous class of chemicals and further research focussing on this class would be required to better define the structural characteristics that could make daphnids more sensitive.

Anilines (aromatic amines) are pinpointed by a number of the classifications implemented by the OECD (Q)SAR Toolbox as an important class of chemicals regarding toxicity and have been shown to be particularly toxic towards daphnids due to a specific mode of action [START_REF] Chen | Toxicity of substituted anilines to Pseudokirchneriella subcapitata and quantitative structure-activity relationship analysis for polar narcotics[END_REF][START_REF] Dom | Aquatic multi-species acute toxicity of (chlorinated) anilines: Experimental versus predicted data[END_REF][START_REF] Neuwoehner | QSAR-analysis and mixture toxicity as diagnostic tools: Influence of degradation on the toxicity and mode of action of diuron in algae and daphnids[END_REF]Vaal et al., 1997a). However, in this dataset, chemicals belonging to the aniline class were not significantly more toxic towards D. magna.

Overall, the results regarding classes of chemicals representing structural alerts were obtained by performing a large number of tests with a high significance threshold (5%) in an exploratory analysis.

The number of significant relationships is relatively small, and due to the fact that several classifications yield some identical or very similar classes, it is difficult comparing the proportion of significant results to the chance level. However, the results of the analysis of the larger database confirmed the results in terms of relative sensitivity of species and also generalised to all three species the effect of biotransformation half-lives observed with the 85-chemical dataset.

Conclusions

Although this dataset is relatively small, it is nevertheless valuable in particular as 85 chemicals were tested with the same methodology, in a relatively short span of time, and on three aquatic species belonging to three different trophic levels. Also the main results obtained on this small dataset were confirmed by the analysis of a larger dataset. The effort put into sharing the data by the Japanese Ministry of Environment is thus very valuable and should encourage the publication of this type data in general as it could increase the relevance of comparative studies. Datasets relative to more homogeneous chemicals would also be welcome.

Exploratory analyses both on presence of structural fragments and on quantitative molecular descriptors illustrated relationships between toxicity and cross-species differences in sensitivity on the one hand and physico-chemical characteristics, especially hydrophobicity, on the other hand. In predicting toxicity towards each species, simple linear regression on logP performed better than PLS regression of toxicity on a very large set of molecular descriptors, thus confirming the importance of hydrophobicity in predicting ecotoxicity. The predictions of toxicity based on the interspecies correlations performed better than the QSAR predictions. In the scope of reducing the number of in vivo tests, the cross-species differences in sensitivity obtained with the outputs of the OECD (Q)SAR Toolbox suggest that aromatic chemicals should be tested first of all on D. magna, as part of a testing strategy. The concentrations tested on fish could then be adjusted according to the toxicity towards D. magna, in order to avoid unnecessary tests and thus reduce the number of vertebrates used in ecotoxicity testing.

Though a number of statistically significant associations were revealed in this paper, the interpretation in terms of biological mode of action is not straightforward due to a lack of knowledge. Though the OECD Toolbox does implement a categorisation based on the mode of action, the classes obtained did not differ in toxicity. Amongst the categorisations implemented by the OECD Toolbox, those focussing on bioaccumulation or biotransformation appeared to be the most interesting in terms of environmental prediction on a whole set of chemicals, partly as they yielded classes of chemicals that were both represented by a sufficient number of chemicals to be tested. However, some categorisations may pinpoint particularly relevant structural fragments but may not be sufficiently represented in this dataset to allow drawing any conclusions. The use of the categorisation schemes implemented in the OECD (Q)SAR Toolbox and of DRAGON molecular descriptors showed that the molecular descriptors are related to some well-defined chemicals classes, such as arenes or chemicals with a methyl functional group, which are relevant considering toxicity or relative sensitivity of species.

Further research is required to investigate the relationship between hydrophobicity, biotransformation half-life and toxicity towards aquatic species. Further developments in QSAR methods could include nonlinear models to improve predictions or the integration of uncertainty on the experimental values (EC50or LC50) in order to estimate the uncertainty on the predictions. The integration of molecular reactivity represented by quantum descriptors in the QSAR models might improve the predictions.

  set. The log(EC50) or log(LC50) in the training set and test set were not significantly different for any of the species (Mann-Whitney test, p>0.05), although the lowest value for P. subcapitata and for D. magna was in the test set. A visual inspection of the position of the test set chemicals on the first two axes of the PCA with all chemicals did not reveal any differences in distribution of the test set and training set.PLS1 regression was performed for each of the three endpoints on the training set with 7-segment cross-validation on the 782 descriptors that were non-constant both on the training and test set. The PLS models had negative cross-validated coefficients of determination (Q² CV ) and very low predictive validity (results not shown).

Figure 1 :

 1 Figure 1: Pairwise interspecies correlations between experimental EC50 or LC50 values for O. latipes,

Figure 2 :

 2 Figure 2: Scatter plots of the two first principal components, highlighting chemicals with a methyl

Figure 4 :

 4 Figure 4: Scatter plots of the two first principal components, showing most sensitive species to each

  Table 1 of the supplementary material shows which chemicals were in the training set and in the test

Table 1 :

 1 Ratios of mean EC50 or LC50 of chemical class over mean EC50 or LC50 of other chemicals; number of cases where the Fisher's exact test p-value was smaller than 0.05; most sensitive species.Highlighted cells relate to p-values under 0.05 meaning that the chemical class is either more or less toxic than the rest of the chemicals (Student tests or Mann-Whitney tests depending on normality of EC50 or LC50 distributions). Tests on most sensitive species were carried out twice, allowing for the two possible cases due to one tie. Tests on least sensitive species were carried out four times, allowing for the four possible cases due to two ties.

	Toxicity of chemical class: ratio EC50 or LC50	Species' relative sensitivity Fisher's exact test on most Fisher's exact test on least
		sensitive species	sensitive species

Table 2 :

 2 Top ten molecular descriptors with highest Pearson correlation coefficient (r) with log(EC50) Signification of the abbreviations are given in Table 2 of the supplementary material.

	or log(LC50) values
		log(EC50)		log(EC50)		log(LC50)
		P. subcapitata	D. magna		O. latipes
		Molecular		Molecular		Molecular
	Rank		r		r		r
		descriptor		descriptor		descriptor
	1	ALOGP² -0.43	MLOGP² -0.57	MLOGP² -0.53
	2	MLOGP² -0.41	MLOGP -0.55	ALOGP² -0.52
	3	Mor27u	-0.4	BLTF96	0.55	ALOGP -0.49
	4	D/Dr10 -0.39	ALOGP² -0.54	MLOGP -0.47
	5	L/Bw -0.38	ALOGP	-0.5	BLTF96	0.47
	6	nR10 -0.35	TPSA(NO)	0.46	EEig14d	0.39
	7	piPC09 -0.35	nHAcc	0.45	Hy	0.37
	8	IC1	0.35	nO	0.43	MAXDN	0.36
	9	ALOGP -0.35	SEigv -0.41	TPSA(NO)	0.36
	10	EEig14d	0.34	TPSA(Tot)	0.41	Ms	0.35

Table 3 :

 3 Coefficient of determination (R², first row), adjusted coefficient of determination (2 nd row) on training set, and externally validated coefficient of determination (Q² ext , 3 rd row) for linear regression models of log(EC50) or log(LC50) for each species on logP (two different calculations by Dragon and their squares), on each other species, on both other species, and on both other species and AlogP².

					Log(EC50) or	Log(EC50)or
		Log(EC50)	Log(EC50)	Log(LC50)
		P. subcapitata	D. magna	O. latipes	log(LC50)	log(LC50)	MlogP MlogP² AlogP AlogP²
					2 species	2 species + AlogP²
			0.525	0.382	0.529	0.529	0.084	0.166	0.037	0.068
	Log(EC50) P. subcapitata		0.518	0.371	0.513	0.505	0.069	0.153	0.022	0.052
			0.555	0.477	0.560	0.560	0.117	0.117	0.258	0.474
		0.525		0.644	0.729	0.731	0.269	0.299	0.129	0.119
	Log(EC50) D. magna	0.518		0.638	0.720	0.717	0.257	0.288	0.115	0.104
		0.540		0.761	0.799	0.816	0.318	0.363	0.465	0.615
		0.382	0.644		0.647	0.655	0.138	0.185	0.111	0.132
	Log(LC50) O. latipes	0.371	0.638		0.635	0.638	0.124	0.172	0.097	0.118
		0.474	0.775		0.777	0.784	0.352	0.459	0.464	0.575

Acknowledgements

This study was supported by the French Ministry in charge of Ecology and Sustainable Development, within the framework of Programme 181 and Programme 190 and by the National Research Agency (ANR) within the project AMORE (Contract number : no. 2009 CESA 15 01).


Guillaume Fayet
email: guillaume.fayet@ineris.fr.

Patricia Rotureau

Laurent Joubert

Carlo Adamo

Development of QSPR models on thermal stability of nitroaromatic compounds considering their decomposition mechanisms

Keywords: Quantitative Structure-Property Relationship (QSPR), nitroaromatic compounds, heat of decomposition, REACH regulation, Density Functional Theory (DFT)

est destinée au dépôt

Development of a QSPR model for predicting thermal stabilities of nitroaromatic compounds taking into account their decomposition mechanisms

Introduction

After initiation from various external stimuli (impact, electric discharge, heat), energetic materials can undergo to decomposition reactions leading to the release of large amount of energy. So, the comprehension of their thermal stability properties is very important in order to evaluate not only their explosive power, but also the related hazards during storage, transportation and handling [START_REF] Grewer | Thermal Hazards of Chemical Reactions[END_REF]. Indeed, many reported (industrial) accidents have been caused by the lack of knowledge about decomposition process [START_REF] Etchells | Why Reactions Run Away[END_REF], particularly those concerning nitro compounds [START_REF] Gustin | Runaway Reaction Hazards in Processing Organic Nitro Compounds[END_REF].

To evaluate the amount of the energy released during the decomposition, calorimetric analyses are performed on chemicals. A typical experimental screening test is the differential scanning calorimetry (DSC) [START_REF] Grewer | Thermal Hazards of Chemical Reactions[END_REF][START_REF] Jones | Evaluation of systems for use in DSC measurements on energetic materials[END_REF][START_REF] Chervin | Method for estimating decomposition characteristics of energetic chemicals[END_REF][START_REF] Yoshida | Safety of Reactive Chemicals[END_REF], which provides heats of decomposition with an error of about 5-10% [START_REF] Ando | Analysis of differential scanning calorimetric data for reactive chemicals[END_REF].

Such measurement represents a pre-selection step for substances which may have explosive properties according to the Recommendations on the Transport of Dangerous Goods [START_REF]Recommendations on the Transport of Dangerous Goods: Manual of Tests and Criteria[END_REF] which are used in REACH (for "Registration, Evaluation, Authorisation and Restriction of Chemicals") [START_REF]the European Parliament and of the Council of 18 December 2006 concerning the Registration, Evaluation, Authorisation and Restriction of Chemicals (REACH)[END_REF] and CLP (for "Classification, labelling and packaging of substances and mixtures") [START_REF]the European Parliament and of the Council of 16 December 2008 on classification, labelling and packaging of substances and mixtures[END_REF] regulations for classification purpose. Within this new European regulatory framework dedicated to chemicals, the evaluation of a tremendous number of substances may be required in a restricted calendar. Unlikely, the complete experimental characterization of chemicals causes not only time, cost and ethical (e.g. animals testing for toxicology and ecotoxicology) problems but also hazards for the case of potentially explosive compounds, (e.g. nitroaromatic compounds) [START_REF] Medard | Les explosifs occasionnels[END_REF]. In this context, the development of alternative tools, notably those issuing from computational chemistry [START_REF] Lewis | Integrating process safety with molecular modeling-based risk assessment of chemicals within the REACH regulatory framework: Benefits and future challenges[END_REF], may be useful and is even recommended in these regulations for the screening and the prioritization of chemicals for experiments.

Among the available methods, Quantitative Structure-Property Relationships (QSPR) represent powerful tools to predict different properties of chemicals. Widely used in 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] Bradbury | Quantitative structure-activity relationships and ecological risk assessment: an overview of predictive aquatic toxicology research[END_REF][START_REF] Selassie | Comparative QSAR and the Radical Toxicity of Various Functional Groups[END_REF] and drug design [START_REF] Grover | Quantitative structure-property relationships in pharmaceutical research -Part 2[END_REF][START_REF] Grover | Quantitative structure-property relationships in pharmaceutical research -Part 1[END_REF], their applications for physico-chemical properties increased since many years [START_REF] Cocchi | Development of Quantitative Structure-Property Relationships Using Calculated Descriptors for the Prediction of the Physicochemical Properties (n D , ρ, bp, ε, η) of a Series of Organic Solvents[END_REF][START_REF] Katritzky | QSPR: the Correlation and Quantitative Prediction of Chemical and Physical-Properties from Structure[END_REF], particularly for the properties of energetic materials [START_REF] Katritzky | QSPR Analysis of Flash Points[END_REF][START_REF] Katritzky | QSPR modeling of flash points: An update[END_REF][START_REF] Keshavarz | Simple method for prediction of activation energies of the thermal decomposition of nitramines[END_REF][START_REF] Keshavarz | Predicting activation energy of thermolysis of polynitro arenes through molecular structure[END_REF][START_REF] Rice | A Quantum Mechanical Investigation of the Relation between Impact Sensitivity and the Charge Distribution in Energetic Molecules[END_REF][START_REF] Saraf | Prediction of reactive hazards based on molecular structure[END_REF][START_REF] Theerlynck | Towards improved models to rationalize and estimate the decomposition temperatures of nitroalkanes, nitramines and nitric esters[END_REF][START_REF] Afanas'ev | Comparative characteristics of some experimental and computational methods for estimating Impact Sensitivity Parameters of Explosives[END_REF][START_REF] Toghiani | Prediction of physicochemical properties of energetic materials[END_REF][START_REF] Politzer | Computational Prediction of Standard Gas, Liquid, and Solid-Phase Heats of Formation and Heats of Vaporization and Sublimation[END_REF][START_REF] Politzer | An electrostatic interaction correction for improved crystal density prediction[END_REF][START_REF] Pospisil | A possible crystal volume factor in the impact sensitivities of some energetic compounds[END_REF]. Their principle consists in developing a mathematic relationship connecting a macroscopic property of a compounds series to microscopic descriptors derived from their molecular structures, using a reliable experimental data set. Once a model is developed and validated, it can be used to predict the macroscopic property for other molecules, with close structures, not yet characterized and even maybe not yet synthesized. Moreover, it might help to understand the investigated phenomena at a molecular scale. To encourage the development and use of QSAR/QSPR models, the Organization for Economic Co-operation and Development (OECD) introduced recently 5 principles for the validation, for regulatory purposes, of these models [START_REF]Principles for the Validation, for Regulatory Purposes, of (Quantitative) Structure-Activity Relationship Models[END_REF]. They require a defined endpoint, an unambiguous algorithm, a defined domain of applicability, appropriate measures of goodness-of-fit, robustness and predictivity and a mechanistic interpretation, if possible. Until now, only few QSPR models were investigated under these principles [START_REF] Papa | Development, Validation and Inspection of the Applicability Domain of QSPR Models for Physicochemical Properties of Polybrominated Diphenyl Ethers[END_REF]. Some works were devoted to predict the thermal stability of chromophores [START_REF] Figueiredo | Chemometric analysis of nonlinear optical chromophores structure and thermal stability[END_REF], polymers [START_REF] Yu | Prediction of the thermal decomposition property of polymers using quantum chemical descriptors[END_REF] or ionic liquids [START_REF] Kroon | Quantum chemical aided prediction of the thermal decomposition mechanisms and temperatures of ionic liquids[END_REF]. Within the framework of energetic materials, some trends connecting the presence of particular functional groups to the decomposition temperatures (T onset ) of potentially explosive compounds were highlighted by Grewer in the early nineties [START_REF] Grewer | The influence of chemical structure on exothermic decomposition[END_REF]. The popular CHETAH software, dedicated to the prediction of reactivity hazards, estimates a maximum heat of decomposition of chemicals using a Benson's group contribution method [START_REF] Benson | Thermochemical kinetics. 2nd edn[END_REF]. Nevertheless, this quantity cannot be directly compared to experimental heats of decomposition since it considers a maximal decomposition of the molecule. More recently, T onset was correlated to the dissociation energies of the weak bonds in nitro molecules [START_REF] Theerlynck | Towards improved models to rationalize and estimate the decomposition temperatures of nitroalkanes, nitramines and nitric esters[END_REF][START_REF] Saraf | Application of Transition State Theory for Thermal Stability Prediction[END_REF].

To the best of our knowledge, the first QSPR type analysis related to the thermal stability of nitroaromatic compounds was realized by Saraf [START_REF] Saraf | Prediction of reactive hazards based on molecular structure[END_REF], based on 19 differential scanning calorimetric data. This study proposed the number of nitro groups (n NO2 ) as the only variable to estimate the decomposition enthalpy (-∆H) with an average error about 8%, an incertitude close to the experimental error. More recently, Keshavarz proposed two QSPR models based on constitutional descriptors to predict the activation energy of thermolysis of nitroaromatics [START_REF] Keshavarz | Predicting activation energy of thermolysis of polynitro arenes through molecular structure[END_REF] and nitramines [START_REF] Keshavarz | Simple method for prediction of activation energies of the thermal decomposition of nitramines[END_REF] with both R²=0.87 in correlation.

In previous works, we proposed preliminary QSPR models to predict the heat of decomposition of 22 nitroaromatic compounds [START_REF] Fayet | On the use of descriptors arising from the conceptual density functional theory for the prediction of chemicals explosibility[END_REF][START_REF] Fayet | On the prediction of thermal stability of nitroaromatic compounds using quantum chemical calculations[END_REF][START_REF] Fayet | QSPR Modeling of Thermal Stability of Nitroaromatic Compounds: DFT vs. AM1 Calculated Descriptors[END_REF], the most robust one presenting good correlation with experimental data (R²=0.98). Nevertheless, the lack of experimental data (only 22 molecules) did not allow the estimation of predictive power and applicability domains of the models. This paper deals with the application of an original approach combining the QSPR methodology with the Density Functional Theory (DFT) for the prediction of the heat of decomposition of a larger series of nitroaromatic compounds (77 molecules), following the OECD recommendations for the validation of such predictive models. The chosen quantum chemical level of theory (DFT) allows for the optimization of molecular structures and, overall, the accurate calculations of a series of chemically-comprehensive descriptors of molecular reactivity (e.g. conceptual DFT descriptors) [START_REF] Fayet | On the use of descriptors arising from the conceptual density functional theory for the prediction of chemicals explosibility[END_REF][START_REF] Fayet | QSPR Modeling of Thermal Stability of Nitroaromatic Compounds: DFT vs. AM1 Calculated Descriptors[END_REF][START_REF] Fayet | On the prediction of thermal stability of nitroaromatic compounds using quantum chemical calculations[END_REF]. Furthermore, the importance of the training set selection is demonstrated and its definition will be in relation with the subjacent microscopic decomposition mechanisms.

Materials and methods

Experimental data set

The choice of the experimental data set is, of course, a critical point of any QSPR analysis.

Since experimental conditions may have a strong influence, all experimental values used in the fitting procedure should be obtained in the same conditions to ensure, at best, their reliability and compatibility.

In this study, a data set of 77 nitroaromatic compounds was considered. It consists in nitrobenzenes derivatives including mono, di and trinitrobenzenes and a large variety of substituents (e.g. nitro, carboxylic acids or halogens). All the heats of decomposition (in tables 1 and 2) were extracted from a single reference [START_REF] Ando | Analysis of differential scanning calorimetric data for reactive chemicals[END_REF] to ensure they were obtained using a single protocol. In the present case, a pressure DSC apparatus was used on 1-2 mg samples in aluminium cells with a pin-hole (i.e. open sample cells) with a heat rate of 10 K/min.

Molecular structures

All molecular structures were calculated using the density functional theory (DFT) with the Gaussian03 package [START_REF] Frisch | Gaussian Inc[END_REF]. Geometry optimizations were performed using the parameter-free PBE0 hybrid functional [START_REF] Adamo | Toward reliable density functional methods without adjustable parameters: The PBE0 model[END_REF] and the 6-31+G(d,p) basis set. Vibrational frequencies were computed at the same theoretical level to ensure that all stable species corresponds to energy minima. These molecular structures were then loaded into CodessaPro software [START_REF] Codessapro | University of Florida[END_REF], including information about geometry, atomic charges, molecular orbitals energies and vibrational frequencies. More than 300 descriptors were calculated and, in addition, some particular external descriptors have been included since their interest has been shown in previous work [START_REF] Fayet | On the use of descriptors arising from the conceptual density functional theory for the prediction of chemicals explosibility[END_REF][START_REF] Fayet | On the prediction of thermal stability of nitroaromatic compounds using quantum chemical calculations[END_REF]. This is notably the case for the number of nitro groups and some conceptual DFT descriptors that have been already correlated successfully with thermal stability of similar nitroaromatic compounds [START_REF] Fayet | On the use of descriptors arising from the conceptual density functional theory for the prediction of chemicals explosibility[END_REF][START_REF] Fayet | On the prediction of thermal stability of nitroaromatic compounds using quantum chemical calculations[END_REF].

The considered molecular descriptors can be divided into several classes: constitutional, topological, geometrical and quantum chemical descriptors. Constitutional descriptors characterize the presence and number of specific atoms, groups or bonds in the molecule (e.g. number of O atoms, single bonds). Topological indices, like Wiener index, are based on the atomic connectivity and give information about the size and branching degree of molecules.

Geometric descriptors, contrary to the previous classes calculated from 2D-structures, stem from the 3D-structure (e.g. molecular volume). Quantum chemical descriptors gather binding, formation, molecular orbital energies, thermodynamic and electronic information, like dipole moment or polarizability, and descriptors characterizing the charge distribution within the molecule (e.g. partial charges). Detailed definitions and information about descriptors can be found in [START_REF] Karelson | Molecular Descriptors in QSAR/QSPR[END_REF].

Statistical analysis

Various data mining tools can be used to develop QSPR models, e.g. artificial neural networks [START_REF] Gasteiger | Neural Networks in Chemistry[END_REF] or genetic algorithms [START_REF] Leardi | Genetic algorithms in chemometrics and chemistry: a review[END_REF]. In this study, multilinear regressions have been computed and the final model has the following general mathematic form:

∑ = + = n i i i X a a Y 1 0 (1)
where Y is the property to predict, X i are the molecular descriptors and a i are the corresponding regression constants.

Equations ( 1) have been obtained using the "best multilinear regression" (BMLR) technique, described in [START_REF] Karelson | Molecular Descriptors in QSAR/QSPR[END_REF] and implemented into CodessaPro software.

The first step of the BMLR analysis consists in reducing the initial set of descriptors by rejecting all descriptors with insignificant variance, to ensure that descriptors are not included by chance in the model whereas they are not related to the investigated property. At the same time, if two descriptors are highly correlated together, only the one presenting the best correlation with the property is kept. This step not only prevents against the introduction of inappropriate descriptors, but also makes the analysis faster, since fewer variables have to be treated.

After that, starting from pairs of orthogonal (i.e. not intercorrelated) descriptors, higher rank models are computed by including orthogonal descriptors successively as soon as an increase in correlation is observed. To the end, the BMLR analysis selects the best models at each rank and the final model must be chosen between them. It has to be sufficiently correlated and, at the same time, ensure against any over-parameterization, which leads to a loss in predictive power for molecules outside the training set.

In this paper, the "breaking point" rule has been used to manage this problem. This method, already successfully used in previous works [START_REF] Katritzky | QSPR modeling of flash points: An update[END_REF][START_REF] Katritzky | QSAR study of mosquito repellents using Codessa Pro[END_REF][START_REF] Fayet | QSPR Modeling of Thermal Stability of Nitroaromatic Compounds: DFT vs. AM1 Calculated Descriptors[END_REF], consists in analyzing the correlation improvement with the number of variables in the model. By plotting the R² values as functions of the number of descriptors, an asymptotic behavior was observed and the improvement of correlation became less significant after a certain rank (∆R² < 0.02-0.03). At this point ("breaking point"), the model is considered as optimal, representing the best compromise between correlation and parameterization.

The robustness and the stability of the models have been evaluated through the square correlation coefficient (R²) and the mean absolute error (MAE). The choice of the descriptors was confirmed by performing a student's t-test at a 95% confidence level. Moreover, they have been validated internally using the cross-validation technique (R² cv ), and externally, using a validation set of data (R² valid ).

Finally, the applicability domain of the models, i.e. the domain in which predictions are reliable, was investigated. In the present paper, the applicability domain estimation was performed based on Euclidean distances in the descriptors space (after a principal component analysis pre-processing) using Ambit Discovery software [51]. The applicability domain was defined to include 95% of the training set molecules and all predictions within this domain are expected to be reliable. This last statement was checked by calculating the correlation coefficients (R² in ) for the molecules of the validation set inside this applicability domain.

Results and discussion

Complete data set

In a first step, the whole data set of 77 molecules has been investigated. This data set is sufficiently large to be divided into a training and a validation set. In order to keep similar distributions for the two ensembles, the entire set was sorted from the smallest to the highest values of heat of decomposition and the validation set was composed by the 3 rd , 7 th , etc. molecules. Thus, distributions of both training and validation sets (58 and 19 molecules, respectively) are quite similar to the entire one, as shown on figure 1. This ratio makes both sets sufficiently important for robust development and validation of the model. Then, the BMLR method was applied on the training set (see table 1). Models including up to 19 descriptors were built. The improvement in R² with the number of descriptors included in the model is represented in figure 2. Based on the "breaking point" rule, the following fourparameter model was considered as the best compromise between correlation and number of descriptors.

-∆H = -282.3 + 333.5 n NO2 -1214.5 E C,avg + 7.4 α -275.6 0 IC avg [START_REF] Etchells | Why Reactions Run Away[END_REF] where n NO2 is the number of nitro groups, E C,avg the average 1-electron reactivity index for a C atom, α the mean polarizability and 0 IC avg the average information content (order 0), a topological index. This model presents the advantage to include descriptors directly related to thermal stability. In fact, the amount of energy released during decomposition is expected to be linked to the loss of nitro groups [START_REF] Brill | Kinetics and mechanisms of thermal decomposition of nitroaromatic explosives[END_REF]. It is then pertinent to find the number of these nitro groups in the model. Moreover, E C,avg characterizes the reactivity of carbon atoms in the molecule and the one connected to leaving nitro groups is expected to be the most reactive in the molecule. Besides, nitro groups also influence the electronic properties of nitro groups, notably the polarizability.

This model is significantly correlated with R²=0.84 and quite stable (R² cv =0.81). If it is less correlated than our previous work based on 22 molecules (R²=0.98) [START_REF] Fayet | QSPR Modeling of Thermal Stability of Nitroaromatic Compounds: DFT vs. AM1 Calculated Descriptors[END_REF], the distribution of experimental data was not as homogeneous as that considered in the present work and the single trinitro molecule of the previous set strongly influenced the regression.

To evaluate the predictive power of this new model, heats of decomposition of the validation set were calculated. As shown on figure 3 At this point, standard experimental uncertainty is not sufficient to explain the poor predictive power of the model. Another factor affecting the accuracy of the model is related to the chemical reactivity. Indeed, the decomposition of nitroaromatic compounds is complex.

Whereas the direct breaking of the carbon nitro bond was sometimes considered as the initiation step of decomposition for nitro compounds [START_REF] Saraf | Application of Transition State Theory for Thermal Stability Prediction[END_REF], numerous studies demonstrated that more complex reaction paths could be specifically involved in nitroaromatic molecules [START_REF] Brill | Kinetics and mechanisms of thermal decomposition of nitroaromatic explosives[END_REF], such as ortho-substituted nitrobenzene derivatives (e.g. o-nitrotoluenes [START_REF] Chen | A Computational Study on the Kinetics and Mechanism for the Unimolecular Decomposition of o-Nitrotoluene[END_REF] and its derivatives [START_REF] Fayet | A theoretical study of the decomposition mechanisms on substituted ortho-nitrotoluenes[END_REF] including 2,4,6-trinitrotoluene [START_REF] Cohen | Mechanism of Thermal Unimolecular Decomposition of TNT (2,4,6-Trinitrotoluene): A DFT Study[END_REF]).

For this reason, models have been developed considering non-ortho and ortho nitroaromatic compounds separately, i.e. a first model is applicable for compounds presenting no substituent in ortho position to the nitro group and a second one for compounds presenting one. So, two models have been developed following the protocol previously used for the entire set with same division into training and validation sets (in tables 1 and 2). Besides, Storm [START_REF] Storm | Sensitivity Relationships in Energetic Materials[END_REF] and Kamlet [START_REF] Kamlet | The relationship of Impact Sensitivity with Structure of Organic High Explosives. II. Polynitroaromatic explosives[END_REF] already considered the importance of the substituents in this position by distinguishing nitroaromatic explosives with an alpha C-H linkage when reporting experimental impact sensitivities and their correlation with the oxygen balance.

Non-ortho compounds

A four-parameter model was developed from the 31 non-ortho molecules of the training set:

-∆H = 0.8 G -3.8 WPSA1 -4255.1 Q max + 26.8 RPCS -251.2

where G is the gravitation index, WPSA1 the weighted positive surface area (from Mulliken calculated charges [START_REF] Mulliken | Electronic Population Analysis on LCAO-MO Molecular Wave Functions[END_REF]), Q max the maximal partial charge in the molecule (calculated according to Gasteiger's method [START_REF] Gasteiger | Iterative partial equalization of orbital electronegativity--a rapid access to atomic charges[END_REF]) and RPCS the relative positive charged surface area (from Zefirov calculated charges [START_REF] Karelson | Molecular Descriptors in QSAR/QSPR[END_REF]). If Q max is related to the nitrogen atoms in nitro groups, the other descriptors are more difficult to link to the decomposition process. G characterizes the molecular shape and the mass distribution within the molecule. WPSA1 and RPCS are charged partial surface area descriptors, developed to encode the features responsible for the polar interactions between molecules. Nevertheless, these last descriptors are related, in a certain way, to the distribution of charge within the molecule which is very influenced by nitro groups. So, they are indirectly related to the properties of this group, which is central in the decomposition of nitroaromatic compounds.

This model is strongly correlated with experimental data (R²=0.90) with an average deviation of 12%, close to the experimental uncertainty (see figure 4 and table 1). The internal validation, by cross-validation, is satisfying with R² cv =0.86. Moreover, it presents good predictivity since calculated values are close to experimental ones for the 11 non-ortho molecules of the validation set (R² valid =0.84, see table 2). Regarding to experimental uncertainty, the performance of this model is satisfying since the predictive power of the model reaches even to 0.86 when excluding the molecule determined to be out of the applicability domain of the model (see table 3).

Ortho compounds

For ortho nitroaromatic compounds, the final model, developed on the 27 molecules of the training set, is another four-parameter equation:

-∆H = 4.1 PNSA1 -3298.8 RPCG -56228 N C,min -1245.5 S ZX/ZX + 1117.8

where PNSA1 and RPCG are the partial negative charged surface area and the relative positive charge (from Zefirov calculated charges), N C,min is the minimum nucleophilic reactivity index for a C atom and S ZX/ZX (ZX shadow / ZX rectangle) characterizes the molecular shape within a ZX plane. The charge related descriptors are the main descriptors in this equation. The only descriptor being directly related to the thermal stability is N C,min , which is influenced by the electronic properties of nitro groups, so to their reactivity with the aromatic ring.

Concerning the performance of the model, the correlation with experimental data is high (R²=0.94) with an average deviation of 17% and the cross-validation procedure exhibits a good robustness (R² cv =0.91) but the predictivity is low since calculated values for the 8 ortho molecules of the validation set deviate about 44% from experiment (R² valid =0.42).

In fact, a reliable model was more expected for the non-ortho than for the ortho molecules since many different reaction paths exist upon the nature of the substituent in ortho position from the nitro group (nitro, alcohol, amino…) [START_REF] Brill | Kinetics and mechanisms of thermal decomposition of nitroaromatic explosives[END_REF] whereas all non-ortho nitroaromatic compounds follow the same decomposition process, C-NO 2 homolysis being evidenced in previous theoretical studies [START_REF] Brill | Influence of the substituent on the major decomposition channels of the NO2 group in para-substituted nitrobenzenes: a tandem mass spectrometric study[END_REF][START_REF] Fayet | Theoretical Study of the Decomposition Reactions in Substituted Nitrobenzenes[END_REF] as the main reaction path (without influence of the substituent nature in meta and para position to the nitro group). Besides the applicability domain analysis confirms this observation since the predictive power (R 2 in ) of models 2 (entire set) and 4 (ortho molecules) are low in their applicability domain, compared to the one of model 3 (non-ortho molecules).

To compare the models developed in this study to previous ones, the heats of decomposition of the validation molecules have been calculated using Saraf's model [START_REF] Saraf | Prediction of reactive hazards based on molecular structure[END_REF] (based on the only number of nitro groups of 19 nitrobenzene derivatives) and to our previous model [START_REF] Fayet | QSPR Modeling of Thermal Stability of Nitroaromatic Compounds: DFT vs. AM1 Calculated Descriptors[END_REF] (based on 22 molecules), as shown in table 2. No satisfactory correlation was exhibited since the average deviations are 36% and 46% (R²= 0.53 and 0.37) for Saraf's model and our previous work, respectively. This may be due to the lack of homogeneity in the data sets, which moreover consider no or only one trinitro compounds. Besides, these models did not distinguished ortho and non-ortho compounds. Finally, the predictive powers of these models are very low, even in their respective applicability domains (as shown in table 3). These previous models present the same limit as Eq. 2 and 4 by considering different decomposition mechanisms in the same data set. So, the knowledge of molecular decomposition paths is a critical point for the prediction of nitroaromatics thermal stabilities.

Accordance with OECD principles

Our best QSPR model, developed for non-ortho compounds (Eq. 3) follows the five OECD principles of validation of QSAR/QSPR models for regulatory use [START_REF]Principles for the Validation, for Regulatory Purposes, of (Quantitative) Structure-Activity Relationship Models[END_REF]: Principle 1: The endpoint is well defined as the heat of decomposition, performed using a pressure DSC apparatus with a 10 K/min heat rate on 1-2mg samples in aluminum cells with pin-hole [START_REF] Ando | Analysis of differential scanning calorimetric data for reactive chemicals[END_REF].

Principle 2: The model is transparent since it consists in a simple multi-linear equation, including four parameters, calculated from well defined DFT-optimized structures at the PBE0/6-31+G(d,p) level.

Principle 3: The model is applicable to all nitrobenzene derivatives, without substituent in ortho position to the nitro group, that are included into an applicability domain defined by interpolation of the training set Principle 4: The performance of the model was estimated by its correlation in the training set (R²=0.90), its robustness by cross validation (R² cv =0.86) and its predictive power into its domain of applicability on an external validation set of molecules (R² in =0.86).

Principle 5: If no mechanistic interpretation can be easily provided from descriptors used in the model, the molecular mechanisms were considered during its development and form part of the defined applicability domain, since the model is not applicable to ortho-substituted compounds (see Principle 3). Indeed, the molecular mechanism involved in the decomposition of the target non-ortho substituted nitroaromatic compounds has been characterized from density functional theory calculations in previous work [START_REF] Kamlet | The relationship of Impact Sensitivity with Structure of Organic High Explosives. II. Polynitroaromatic explosives[END_REF].

As this model satisfies all requirements of OECD principles, it could be used as an efficient alternative to experimental characterization as a first screening test to evaluate if a target nitroaromatic compound may have explosive properties.

Conclusion

A set of 77 nitroaromatic compounds was considered aiming to develop a reliable QSPR model for the prediction of their heats of decomposition. Molecular structures, calculated at a DFT level of theory, were described using more than 300 descriptors.

Considering the entire set of data led to significant correlation in the training step (R²=0.84) but it failed in predicting the heats of decomposition of external molecules (about 32% in deviation). If experimental uncertainty (e.g. use of open sample cells) contributed to the lack of predictive power, this paper demonstrated that the data set selection had to take into account chemical reactivity during the decomposition process. Indeed, considering a set of nitro compounds without any ortho substituent to the nitro group led to a performant model with significant correlation between calculated and experimental heats of decomposition of molecules out of the training set (R² valid =0.84). Besides, within its own domain of applicability, the predictivity of this model is very high (R² in =0.86). The improvement of accuracy, compared with the model for the entire set, is related to the fact that various specific decomposition mechanisms involve when substituents are in ortho position from the nitro groups whereas all non-ortho compounds decomposed following the same reaction path (homolysis of the C-NO 2 bond).

Finally, this model is, to our knowledge, the most reliable QSPR model predicting the heat of decomposition within the class of nitroaromatic compounds and the first dedicated to energetic compounds following all OECD requirements for regulatory use.

Figure captions Figure 1 .Figure 2 .

 captions12 Figure captions

Figure 3 .

 3 Figure 3. Experimental versus calculated heats of decomposition (in kJ/mol) of nitroaromatic compounds according to Eq. 2 (the molecules of the training and validation sets are represented in wide and plain circles respectively).

Figure 4 .

 4 Figure 4. Experimental versus calculated heats of decomposition (in kJ/mol) of non-ortho nitroaromatic compounds according to Eq. 3 (the molecules of the training and validation sets are represented in wide triangles and plain circles respectively).

  the applicability domain of the model (see table3). So, the predictive power of the model is low within its own applicability domain (R² in =0.43).

	and table 2, the model gives only low correlation
	with an average deviation of 32% with experiments (R² valid =0.43). In particular, the calculated

values for 4-fluoro-2-nitrotoluene, 2,4-and 2,6-dinitrobenzoic acids present large errors: 64, 83 and 206 %, respectively (see table 2). Nevertheless, no molecule of this validation set was excluded from

Table 1 .

 1 Experimental and calculated heats of decomposition (-∆H) in kJ/mol for the

	molecules of the training set
	molecules	Exp [7]	Eq. 2	Eq. 3	Eq. 4
		non-ortho compounds
	2-amino-4-nitrophenol	130	238	173	-
	3-nitrotoluene	149	238	212	-
	2-amino-5-nitrophenol	153	201	239	-
	4-nitrotoluene	213	253	192	-
	4-nitrophenol	232	183	235	-
	3-nitroanisole	243	223	288	-
	3-nitrobenzoic acid methyl ester	256	334	277	-
	2,6-dichloro-4-nitroaniline	264	394	284	-
	4-nitrophenetole	270	347	249	-
	4-nitrophenylhydrazine	277	389	279	-
	3-nitrophenol	283	165	227	-
	3-nitrobenzoic acid	289	267	372	-
	4-nitroacetophenone	291	348	343	-
	4-nitrobenzyl alcohol	292	223	272	-
	4-nitrobenzoic acid methyl ester	302	329	264	-
	4-nitro-2-toluidine	306	287	315	-
	4-nitrobenzamide	319	275	321	-
	4-nitrobenzyl chloride	337	675	333	-
	4-nitroaniline	347	298	308	-
	3-nitroaniline	350	302	317	-
	3-nitrophenylacetic acid	358	397	347	-
	4-nitrobenzhydrazide	362	415	335	-
	3-nitroacetoanilide	369	289	394	-
	2-amino-4-nitroanisole	375	339	325	-
	4-nitroacetoanilide	387	343	372	-
	4-nitrobenzoyl chloride	408	463	303	-
	3-nitrocinnamic acid	414	314	417	-
	4-nitrobenzaldehyde	421	495	394	-
	3,5-dinitrobenzonitrile	654	699	698	-
	3,5-dinitrobenzoic acid	674	658	679	-
	3,5-dinitrobenzylchloride	711	682	673	-
		ortho compounds
	5-chloro-2-nitrobenzotrifluoride	40	96	-	7
	2-nitrophenol	123	176	-	172
	2-nitrophenylacetic acid	175	288	-	266
	2-nitrotoluene	182	223	-	211
	2-nitroanisole	230	203	-	276
	2-nitrobenzamide	256	317	-	341
	2-nitrobenzoic acid	271	212	-	247
	2-nitrobenzoic acid methyl ester	274	284	-	285