Comparative QSTR Study Using Semi-Empirical and First Principle Methods Based Descriptors for Acute Toxicity of Diverse Organic Compounds to the Fathead Minnow

Several quantum-mechanics-based descriptors were derived for a diverse set of 48 organic compounds using AM1, PM3, HF/6-31+G, and DFT-B3LYP/6-31+G (d) level of the theory. LC50 values of acute toxicity of the compounds were correlated to the fathead minnow and predicted using calculated descriptors by employing Comprehensive Descriptors for Structural and Statistical Analysis (CODESSA) program. The heuristic method, implemented in the CODESSA program for selecting the ‘best’ regression model, was applied to a pre-selection of the most-representative descriptors by sequentially eliminating descriptors that did not satisfy a certain level of statistical criterion. First model, statistically, the most significant one has been drawn up with the help of DFT calculations in which the squared correlation coefficient R2 is 0.85, and the squared cross-validation correlation coefficient RCV2 is 0.79. Second model, which has been drawn up with the help of HF calculations, has its statistical quality very close to the DFT-based one and in this model value of R2 is 0.84 and that of RCV2 is 0.78. Third and fourth models have been drawn up with the help of AM1 and PM3 calculations, respectively. The values of R2 and RCV2 in the third case are correspondingly 0.79 and 0.66, whereas in the fourth case they are 0.78 and 0.65 respectively. Results of this study clearly demonstrate that for the calculations of descriptors in modeling of acute toxicity of organic compounds to the fathead minnow, first principal methods are much more useful than semi-empirical methods.

descriptors in modeling of acute toxicity of organic compounds to the fathead minnow, first principal methods are much more useful than semi-empirical methods.

Introduction
Many of QSAR studies are based on the assumption that molecules from the same chemical domain will behave in a similar manner, so that QSAR models drawn up with the analogical molecules are hypothesized to exhibit better performance than that derived from miscellaneous data set. The traditional approach to QSARs for acute toxicity of organic compounds to the fathead minnow is the modeling of the activity of homologous or congeneric series of chemicals such as nitroaromatics [1], alkylamines [2], halogenated hydrocarbons and phenols [3], and chlorobenzenes and chloroalinines [4]. This congeneric series approach is conservative. Often, such chemicals have a single functional group or toxicophore and an alkyl moiety of variable size. Some other studies [5][6][7][8][9][10][11] by using diverse molecule sets have usually relied on dividing a molecule set into subgroups (chemical classes) by clustering the molecules based on their mode of action. Then, local QSTRs built up for each subgroup are applicable only to certain mode of action. It is worthy mention here that there has been a successful effort to draw up a global QSTR model by using a single descriptor, namely, logarithm of 1octanol/water partition coefficient LogP. This model is applicable to quite miscellaneous data set, but still counts quite a big number of molecules as outliers [12]. Hydrophobicity of a molecule is characterized by LogP which is directly related to bio-uptake of chemicals by fish or many other organisms. It has been successfully used for the modeling of acute toxicity of chemicals with different modes and mechanism of toxic action to Pimephales promelas, combined with additional parameters such as energy of lowest unoccupied molecular orbital (E LUMO ) [13] and maximum superdelocalizability (S max ), which is a molecular orbital parameter that quantifies the electro (nucleo) philicity of a molecule [14]. Developing a better QSTR for the modeling of acute toxicity of diverse chemicals is a subject of interest due to its demand by the organizations such as OECD (Organization for Economic Co-operation and Development) or EC (European Communities) to use the QSTR model for regulatory purpose.
The aim of the present study is two folds. The first one is to build QSTR multiple regression model using quantum-mechanics-based molecular descriptors that correlate and predict the LogLC 50 value of acute toxicity of 48 compounds to the fathead minnow. LC 50 (mg/l), aquatic toxicity on Pimephales promelas expressed as the chemical concentration at which 50% lethality is observed in a test batch of fish within a 96 h exposure period. Molecules used in this study are quite a diverse set and were taken from a study [12]. However, they were not strictly selected to ensure that they are sufficiently diverse. The second aim of this study is to compare the accuracy of semi-empirical and first principle methods for calculation of molecular descriptors. AM1 [15] and PM3 [16,17] are fast in computation, well suited to organic compounds, and belong to semi-empirical method family. These methods have been traditionally used to calculate the optimized 3D geometry and quantum mechanics descriptors of molecules in most of QSAR studies. Some previous comparative QSAR works [1,[18][19][20][21][22] have shown that using descriptors calculated by HF [23][24][25] or DFT [26] together with B3LYP [27] hybrid function instead of semi-empirical AM1 or PM3 methods improve the accuracy of the results that lead more reliable QSARs. On the other hand, there is an interesting comparative QSTR study relevant in this area [28]. In that study, a huge molecule set (568 molecules) has been used to establish QSTR models. These QSTR models have been built up from descriptors which were calculated using two different theory levels namely AM1 and DFT/B3LYP (6-31G**). Their study has shown that the choice of the precise but time-consuming DFT/B3LYP method does not have an advantage over AM1 method for the quality of the derived QSTRs.

Computational details
For all molecules studied here, 3-D modeling and calculations were performed using the Gaussian 03 quantum chemistry package [29]. To save in computational time, initial geometry optimizations were carried out with the molecular mechanics (MM) method using Amber force field. The lowest energy conformations of the molecules obtained by the MM method were further optimized by the DFT method by employing Becke's three-parameter hybrid functional (B3LYP) and the 6-31+G (d) basis set; their fundamental vibrations were also calculated using the same method to check if there were true minima. All the computations were carried out for the ground states of these molecules as singlet state. The lowest energy conformations of the compounds obtained using DFT were used as an input geometry for the calculations for HF/6-31+G, AM1 and PM3 methods. (CODESSA PRO) Comprehensive Descriptors for Structural and Statistical Analysis, Version 2.7.2 [30], was used for extracting descriptors of quantum mechanics and 3D geometry of the compounds from Gaussian 03 output files. CODESSA PRO enables the generation of hundreds of molecular descriptors (constitutional, topological, and quantum mechanical) from a loaded 3D geometry, and uses diverse statistical structure property/activity correlation techniques for the analysis of experimental data in combination with calculated molecular descriptors.  A QSAR/QSTR model can be developed for a given set of molecules by using a various types of descriptors. Sometimes, a model might have very good statistical parameters, but still not suffice to explore the mechanism of interaction between the ligand and receptor mechanistically. Building a model with physically interpretable descriptors is an important task for value of a QSAR/QSTR work. In this study, we aimed to draw up a QSTR model by using quantum mechanically calculated thermodynamical descriptors by virtue of which obtained models are usually mechanistically interpretable. About 50 thermodynamical descriptors depending on the number of atoms in a molecule were calculated using CODESSA PRO and Gaussian 03 packages. The heuristic method [29] implemented in CODESSA PRO was used to build up a multi-able regression model. By this method, a pre-selection of descriptors is accomplished. All descriptors are checked to ensure the following: (a) the values of each descriptor are available for each structure, and (b) there is a variation in these values. The descriptors for which values are not available for every structure in the data in question are discarded. Descriptors having a constant value for all structures in the data set are also discarded. A printout showing the values of descriptors discarded in this manner is provided. Thereafter, the oneparameter correlation equations for each descriptor are calculated. To further reduce the number in the "starting set" of descriptors, the following criteria are applied and a descriptor is eliminated if any of the following conditions are met with: (a) the F-test's value for the one-parameter correlation with the descriptor is below 1.0, (b) the squared correlation coefficient of the one-parameter equation is less than R 2 min by default 0.01, (c) the parameter's t-value is less than t 1 (where R 2 min 0.1 by default and t 1 1.5 by default are user defined values), (d) the descriptor is highly inter-correlated (above r full , where r full is a user specified value by default 0.99), with another descriptor. All the remaining descriptors are then listed in decreasing order according to the correlation coefficient of the corresponding oneparameter correlation equation. All two-parameter regression models with remaining descriptors are developed and ranked by the regression correlation coefficient R 2 . A stepwise addition of the further descriptors' scales is performed to find the best multi-parameter regression models with the optimum values of statistical criteria (highest values of R 2 , the cross-validated, R 2 CV and the F value). R 2 CV , the 'leave one out' (LOO) squared cross-validated coefficient, is a practical and reliable method for testing the predictive performance and stability of a regression model. LOO approach consists in developing a number of models with one sample omitted at a time. After developing each model, the omitted data are predicted and the differences between the experimental and predicted activity values are calculated. Then the R 2 CV is calculated according to the following formula [31]:

Theory
Among the thermodynamical descriptors, translational entropy (at 300 K), principal moment of inertia A, highest normal mode of vibrational frequency, and lowest normal mode of vibrational frequency were involved in the models that are presented in this study. Thermodynamical properties of a molecule arise from the energetics of vibrational frequencies. This connection is based upon partitioning the total energy of a macroscopic system among the constituent molecules. Translational entropy (at 300 K) is defined as [32]; where V is the volume of the system, N is the Avogadro's number, h is the Planck constant, m is the molecular mass, and kT is the Boltzman temperature. Highest normal mode of vibrational frequency and lowest normal mode of vibrational frequency are actually not frequencies; they are wavenumbers (in cm -1 unit). It is customary to call normal modes of vibration of molecule as frequency in infrared and Raman spectroscopy. Definition of normal mode of vibration arises from quantum mechanical harmonic oscillator model of a diatomic molecule. In this model, energy of the vibrational states is given as [33], where h is the Planck constant,  is the vibrational quantum number (0, 1, 2,), and v is the classical vibrational frequency given by, where k is the force constant of the chemical bond and  is the reduced mass for nuclei of two atoms.
More commonly, equation 5 is used as vibrational wavenumber () form rather than frequency form, where where  is the vibrational wavenumber, c is the velocity of light. Final descriptor involved in our model is principal moment of inertia A (I A ) that is obtained from the 3D-cooordinate of the atoms in the given molecule. Its definition is given as [34], where m i are the atomic masses and r ix denotes the distance of the i-th atomic nucleus from the main rotational axes, x. I A characterizes the mass distribution in the molecule.

Results
The assessment of toxcity of a hypothetical compound is a subject of interest. The QSAR/QSTR method saves time and cost in determining the toxicity of a series of newly synthesized compounds with the help of toxicity of previously known compounds. Forty-eight compounds have been taken in this study, and their toxicity (LogLC 50 ) and calculated logarithm of 1-octanol/water partition coefficient (LogP) values to fathead minnow have been taken from the literature [12] and are given in Table 1. Among the quantum mechanically calculated descriptors, translational entropy (at 300 K), principal moment of inertia A, and highest normal mode of vibrational frequency and lowest normal mode of vibrational frequency have been identified which are capable of modeling the toxicity and the structure of a molecule. The data matrix of these descriptors obtained from first principal (HF and DFT-B3LYP) and semi-empirical (AM1 and PM3) methods calculations are shown in Table 2 to 5. By using DFT-based descriptors, several equations were generated by using all the variables and the statistically best model that we have obtained is four-parameters equation, which is as follows: DFT-LogLC 50 = 36.37 -1.11S tr +1.65x10 -3  H -2.29x10 -3  L -0.34I A N=45, R 2 =0.85, R 2 CV =0.79, F=60.34, and s 2 =0.27 (7) where N is the number of compounds included in the model, R 2 is the squared correlation coefficient, R 2 CV is the squared cross-validation correlation coefficient, F is the Fisher test for significance of the equation and s 2 is the standard deviation of the regression. The statistical quality of the above equation is good as evident from its correlation coefficient R 2 value = 0.85 and a cross-validation coefficient R 2 CV value = 0.79. The predicted toxicity of the compounds is given in Table 2  As can be seen in equation 11, statistical quality of DFT-based model was increased dramatically by adding LogP to the model. Table 4. AM1-based descriptors and predicted toxicity of the compounds by Eq 9.   The above linear regression models obtained using the descriptors from DFT and HF calculations are much better than those obtained from semi-empirical AM1 and PM3 methods. Introducing LogP to DFT-based model, there is a rapid rising in statistical quality of the regression equation. Figure 1 gives the plots of observed LogLC 50 versus the predicted LogLC 50 by equations (7)-(11) for compounds.

Discussion
QSAR/QSTR model quality depends on the reliability of the dataset (i.e. uncertainty in toxicological and physicochemical and/or structural data). The authors usually have to rely on experimental toxicological data set that is taken from literature. The data are assumed to provide a uniform measure of toxicity for all of the compounds studied. When the uncertainties in physicochemical and/or structural data are considered, the accuracy of the descriptors is an important element for the QSAR/QSTRs. The computational level of theory is a major task for the accuracy of descriptor calculation. Above presented results clearly demonstrate the effect of used level of theory for calculations. Semi-empirical methods such as AM1 and PM3 use the empirical or experimental parameters to deal with the Schrödinger equation and omit some molecular integral calculations, so they are much faster than the first principle HF and DFT-B3LYP schemes. Therefore, they are utilized more widely in the calculations of molecular properties. But accuracy of their results is inferior to ab initio or DFT methods. The best statistical quality of the equations obtained in this study was DFTbased one with the correlation coefficient R 2 value = 0.85 and a cross-validation coefficient R 2 CV value = 0.79. The second best equation was HF-based one, and in this model R 2 is 0.84 and 2 CV R is 0.78. This result is as expected. HF does not include the effects of an instant electronic correlation, whereas DFT-B3LYP does, so that HF is inferior to DFT in theory. The statistical quality of equations obtained from the semi-empirical methods was lower than that of equations obtained from the DFT and HF methods as expected. AM1 and PM3 based models have given similar statistical fits, but PM3 based model has a lower predictive power as evident from its lower value of R 2 CV =0.66 in contrast to that of AM1 value of R 2 CV =0.72.
For all of the above mentioned models, the most representative descriptor with the highest coefficient is the translational entropy (at 300 K) S tr, which negatively correlates to LogLC 50 . It should be noted that all of the thermochemical properties of a molecule such as S tr arise from the energetics of vibrational frequencies. This connection is based upon partitioning of the total energy of a macroscopic system among the constituent molecules. Other two descriptors involved in the models are lowest normal mode of vibrational frequency,  L and highest normal mode of vibrational frequency,  H .  L correlates negatively to LogLC 50 in all the models, whereas  H correlates positively to LogLC 50 in DFT, HF, and DFT with LogP models. In several QSAR studies, fundamental vibrational frequencies of molecules have been used as descriptors [35][36][37][38][39]. They suggested that the eigen value ('EVA') descriptors are derived from fundamental IR and Raman range of molecular vibrational frequencies. Vibrational frequency of a molecule is sensitive to 3D structure. The idea behind the use of such data as descriptors was that a significant amount of information pertaining to molecular properties might be contained within the molecular vibration wave function [35]. Another descriptor involved in the models in this study is principal moment of inertia A, (I A ). It is a geometrical descriptor and originates from the rigid rotator approximation model of molecules. I A is sensitive to 3D structure and characterizes the mass distribution in the molecule. It correlates negatively to LogLC 50 in all models.
Results of present study demonstrate that quantum mechanically calculated thermo chemical descriptors, S tr ,  H , and  L jointly with geometrical descriptor, I A are capable of modeling the acute toxicity of the compounds to the fathead minnow. First principle DFT and HF methods led to statistically better models than that of semi-empirical AM1 and PM3 methods. This result is normal because DFT and HF can calculate molecular properties such as optimized geometry and spectroscopic properties more accurately than semi-empirical methods. Our results are in disagreement with the conclusions of the other comparative study [28] which has concluded that the use of DFT/B3LYP does not have an advantage over AM1 for the quality of the derived QSTRs. The disagreement between two studies may result from several reasons. It may be due to the difference of the size of molecule sets between two studies. Another possible reason is that the models were restricted to build up using only quantum chemical and thermodynamical descriptors in our study whereas Natzeva et al.
[28] has used considerably large amount of descriptors to derive their models. Finally, LogP has been inserted as an additional descriptor into the statistically best model (DFTbased one). This resulted in an increase of statistical quality of the model for the parameters (R 2 from 0.85 to 0.89, F from 60.34 to 70.82, and s 2 from 0.27 to 0.19). As mentioned in introduction section, LogP itself has been used as single descriptor for many QSTR models. Most of the compounds used in this study act in narcotic mode of action. Narcosis is a general term that describes noncovalent interaction between xenobiotics and cellular membranes. Whereas it is generally accepted that narcosis is the result of the accumulation of the compounds in cell membranes that disturbs their function, the exact mechanism is not known yet [40]. LogP characterizes the hydrophobicity of a molecule that is directly related to bio-uptake of chemicals by fish or many other organisms. Results presented in this study demonstrates that quantum mechanically calculated thermo chemical descriptors in combination with LogP are capable of modeling the acute toxicity of a quite diverse set of 48 compounds to the fathead minnow.  H , is the highest vibrational wavenumber and  L , is the lowest vibrational wavenumber.
b Residual is the differences between O-LogLC 50 and P-LogLC 50 values.

Conclusion
This QSTR study has been based on quantum mechanically calculated descriptors (such as entropy (at 300 K), principal moment of inertia A, highest normal mode of vibrational frequency, and lowest normal mode of vibrational frequency) and on the acute toxicity of the 48 organic compounds to the fathead minnow. All of these descriptors are sensitive to 3D structure of the molecules. The reliability of this study has been tested by four different methods, namely, AM1, MP3, HF, and DFT/B3LYP. A comparison of all the methods indicates that the DFT/B3LYP method is more reliable than others and has a high predictive power. Introduction of LogP as an additional descriptor into DFT-based model has resulted in an increase of statistical parameters of the model.