Setschenow Constant Prediction Based on the IEF-PCM Calculations

Article pubs.acs.org/IECR

Setschenow Constant Prediction Based on the IEF-PCM Calculations Xinliang Yu†,‡ and Ruqin Yu*,† †

State Key Laboratory of Chemo/Biosensing and Chemometrics, College of Chemistry and Chemical Engineering, Hunan University, Changsha, Hunan 410082, People’s Republic of China ‡ College of Chemistry and Chemical Engineering, Hunan Institute of Engineering, Xiangtan, Hunan 411104, People’s Republic of China S Supporting Information *

ABSTRACT: The Setschenow constant Ksalt of a compound in NaCl solution is an important parameter. The solvent environment influences the geometries, energies, charge distributions, and other properties of solutes. The integral equation formalism polarizable continuum model (IEF-PCM) for solvent effects was used to optimize molecular geometries with the density functional theory method combined with Becke’s three-parameter hybrid functional and Lee−Yang−Parr’s gradientcorrected correlation functional at 6-31G(d) level. Single-point energy calculations and natural bond orbital analyses were carried out with the same method. After 1672 molecular descriptors generation, four descriptors were selected to develop models for Ksalt of 101 organic compounds, by using the genetic algorithm method together with multiple linear regression (MLR) technique. The optimal MLR and support vector machine models of Ksalt have the mean root-mean-square errors of 0.0287 and 0.0227, respectively. Compared with previous models, the two models in this work have better predictive performance. Results of the study suggest that calculating molecular descriptors from IEF-PCM to predict the Setschenow constants Ksalt of organic compounds in NaCl solution is feasible.

■

the previous model,3 by applying several chemometric tools including multiple linear regression (MLR),6,7 artificial neural networks (ANNs),6,7 and support vector machines (SVMs).6 But the SVM model of Ksalt had the squared correlation coefficient of 0.7511 and the root-mean-square (rms) error of 0.0385 for the test set, which were poorer than the MLR and ANN models.6 The descriptors used to develop quantitative structure− activity relationship (QSAR) models for Ksalt were calculated from a single molecule in vacuum (i.e., a gas-phase molecule).3,6,7 Solute molecules in the gas-phase and in solution might be different in the geometries, energies, and charge distributions because of the solvent effects. These effects play a key role in physical process and chemical reactions in solution. In this work, the integral equation formalism polarizable continuum model (IEF-PCM) for solvent effects was applied to molecular geometry optimizations, single-point energy calculations and natural bond orbital analyses.8−10 Molecular descriptors derived from optimized molecular structures were then used to develop MLR and SVM models predicting Setschenow constants Ksalt.

INTRODUCTION The aqueous solubility of organic compounds is of importance in the fields of chemical processes and reactions, environmental hazard assessment, and pharmaceutical and pesticide design. For example, the aqueous solubility of pharmaceuticals strongly affects their transport, release, absorption, and tissue distribution.1,2 With the advent of combinatorial chemistry and highthroughput synthesis, the demand for fast methods for the accurate assessment of molecular properties prior to organic synthesis increases. Therefore, it is very significant to develop a reliable predictive model for aqueous solubility of compounds.1−3 When there are strong electrolyte salts in an aqueous solution, the solubility of organic compounds in water can increase or decrease depending on the polarity of both the solute and the salts. NaCl salt solution increases the polarity of water, which increases the “squeezing out” effect of water on nonpolar solutes according to the “like dissolves like” principle of solubility. This “salting-out” effect on the solubility of organic compounds in an aqueous salt solution can be described by the Setschenow equation:4 log(Ssalt /Swater) = −K saltCsalt

■

MATERIALS AND METHODS Experimental Data. Supporting Information, SI, Table 1S shows 101 organic compounds and their respective experimental Setschenow constants Ksalt in NaCl solution, which are taken from the literature.2 The data set have been studied by

(1)

where Ssalt and Swater are the solubilities of the organic compound in an aqueous salt solution and in water, respectively, Csalt is the molar concentration of an electrolyte, and Ksalt is the Setschenow constant. Theoretical treatments to predict Ksalt have been carried out by using the molar volume,5 intrinsic solubility,1,2 octanol− water partition coefficient,1,2 connectivity index,3 molecular weight,3 and other molecular descriptors.6,7 The models developed by Xu et al. were tested to be more accurate than © 2013 American Chemical Society

Received: Revised: Accepted: Published: 11182

January 1, 2013 June 19, 2013 July 23, 2013 July 23, 2013 dx.doi.org/10.1021/ie400001u | Ind. Eng. Chem. Res. 2013, 52, 11182−11188

Industrial & Engineering Chemistry Research

Article

atom (QMH), and the most negative Mulliken charge (qM−). The total dipole moment μ reflects the global polarity of a molecule. EHOMO characterizing the susceptibility of the compound toward attack by electrophiles is directly related to the ionization potential. Correspondingly, ELUMO characterizing the susceptibility of the compound toward attack by nucleophilesis is directly related to the electron affinity. The total energy ET can reflect the molecular size and represent a bulk measure of the solute’s ability to participate in nonspecific interactions with the stationary phase. Charge-based descriptors can denote molecular chemical reactivity or weak intermolecular interactions. Commonly, the minimum (most negative) and maximum (most positive) atomic partial charge in the molecule or the minimum or maximum partial charges for particular types of atoms are used. Thus, natural atomic charges (QH and q−) and Mulliken charges (QMH and M−), respectively, based on the natural population analysis and the Mulliken population analysis were calculated in this work. In addition, we calculated 1664 molecular descriptors from each optimized molecular structure using Dragon software (Version 5.4).16 In total, 1672 descriptors were calculated for each molecule in SI Table 1S. In order to do a comparison with previous models based on respective 4 and 5 descriptors,6,7 4-descriptor models were developed in this work. Stepwise MLR in SPSS 11.5 and the genetic algorithm (GA) in BuildQSAR17 were used to develop MLR models. BuildQSAR17 can deal with a relatively small number of descriptors. First, stepwise MLR in SPSS 11.5 was used to investigate the correlation between 101 Ksalt in SI Table 1S and the 1664 descriptors calculated with Dragon software (Version 5.4)16 and obtain a descriptor subset (Subset I) comprising 15 descriptors. Then GA method, together with MLR technique in the BuildQSAR17 was applied to select Subset II from Subset I and the eight quantum chemical descriptors. During the process of GA analysis, following parameters were controlled: total variables in models being 4; number of models per generation being 30; mutation probability being 35%; and number of generations being 200. In addition, we used stepwise MLR in SPSS 11.5 to directly analyze the 101 Ksalt in SI Table 1S and 1672 descriptors, and obtained another Subset III. Support Vector Machine. Support vector machine (SVM) is a powerful state-of-the-art data mining tool and possesses many attractive features and promising empirical performance.18−23 The main advantage of SVM originates from the algorithm based on guaranteed risk bounds of statistical learning theory, i.e., the so-called structural risk minimization principle. In SVM method, systems map training data in the input space into a high dimensional feature space and subsequently carry out the linear regression in the feature space. SVMs approximate the regression function in the following form:6,19,22,23

other researchers.2,3,6,7 These compounds include alkanes, alkenes, acetates, ketones, acids, aniline, amides, aromatics, and so on. The reported Ksalt values ranged from −0.068 (cystine, No. 15) to 0.354 (1,2-benzanthracene, No. 4). To ensure the test set for Ksalt (SI Table 1S) covering the full range as the data used for training, we randomly split samples into a training set (51 compounds) and a test set (50 compounds). Molecular Descriptor Derivation. Polarizable continuum models (PCMs) are widely used to model quantum effects of molecules in liquid solutions.8,9 In PCMs, the solute molecule is located inside a cavity, which is surrounded by a homogeneous dielectric (the solvent). The solute−solvent interactions between the charge distributions which compose the solute and the dielectric are limited to those of electrostatic origin. IEF-PCM for solvent effects is a convenient approach to resolve the electrostatic problem between the solute and polarization charges. This method creates the solute cavity via a set of overlapping spheres and can be used in all cases where the Green’s function for the considered environment is known.8−10 Although it is not always possible to have analytical Green functions, the Green function may be effectively built numerically, and thus the IEF approach can be generalized to many other environments.11 Actually, IEF is a method more general than other members of the PCM family, such as DPCM (D stands for dielectric), or COSMO (conductor-like screening model).8 Generally, it can describe accurately the charge distribution of solute. In addition, density functional theory (DFT) is becoming popular in the prediction of a broad range of properties of molecules in the ground state because of two important reasons. One is its lower computational cost when compared with conventional ab initio methods with inclusion of electron correlation. The other is the fact that DFT includes the effects of electron correlation at some reasonable level by using a functional of electron density to model exchange and correlation.12 We used IEF-PCM and DFT because of the combination of lower computational cost with reasonable accuracy. The molecular structures were sketched using ChemBioDraw Ultra 11.0, then preliminarily optimized using ChemBio3D Ultra 11.0 based on molecular mechanics (MM2) until the rms of gradient value became smaller than 0.1 kcal/molÅ. The energy minimized molecules were subjected to reoptimization using DFT with IEF-PCM in the self-consistent reaction field (SCRF) framework implemented in Gaussian 09 program,13 with the default solvent (water) and solvent parameters. Becke’s three-parameter hybrid exchange functional combined with the Lee−Yang−Parr correlation functional (B3LYP),14 with 6-31G(d) basis set,15 was used to reoptimize molecular geometries. Besides the natural bond orbital analyses, frequency calculations were carried out to ensure the nature of the stationary points found (no imaginary frequencies) and to calculate the zero point or thermal energy by employing IEFPCM/B3LYP/6-31G(d), after full geometry optimization. Eight quantum chemical descriptors were derived for each optimized structure of 101 molecules in SI Table 1S. They were the total dipole moment (μ), the energies of the highest occupied molecular orbital (HOMO), and the lowest unoccupied molecular orbital (LUMO) (EHOMO and E LUMO), the total energy (ET), the most positive natural atomic charge (nuclear charge minus the summed natural populations of the natural atomic orbitals on the atom) on hydrogen atom in a molecule (QH), the natural atomic charge of the most negative atom (q−), the most positive Mulliken charge on hydrogen

n

f (x ) =

∑ ϕ(xi)w + b

(2)

i

where n is the total number of input−output pairs, φ(x) is a mapping function, x is the input vector, f(x) is the output. The coefficients w and b are estimated by the following minimization: min J(w , ξ , ξ*, b) =

w , b , ξ , ξ*

11183

1 || w ||2 + C ∑ (ξi + ξi*) 2 i

(3)

dx.doi.org/10.1021/ie400001u | Ind. Eng. Chem. Res. 2013, 52, 11182−11188


Article

characteristics are listed in Table 1. The optimal MLR model was obtained from the training set in SI Table 1S:

subject to the following: yi − ϕT (xi)w − b ≤ ε + ξi

(4)

ϕT (xi)w + b − yi ≤ ε + ξi*

(5)

Table 1. Characteristics of Descriptors Selected in MLR Model Based on 101 Compounds

where C is a penalty parameter to determine the trade-off between the training error and the model flatness. ε is a prescribed parameter of the ε-insensitive loss function. If the predicted value is within the tube, then the loss is zero and the loss function ignores errors. While if predicted point is outside the tube, the loss is the magnitude of the difference between the predicted value and the radius ε of the tube, and the loss function penalizes errors.6,19 ξ and ξ * are positive slack variables for the data points. Usually, support vector regression (SVR) uses the εinsensitive loss function to measure the empirical risk (training error): ⎧|f (x) − y| − ε ,(|f (x) − y| ≥ ε) |f (x) − y|ε = ⎨ ⎩ (|f (x) − y| < ε) 0

i

i

(6)

(7)

(8)

QH

0.113 0.012 0.000 9.435 2.893

−0.107 0.031 0.001 −3.502 1.497

(10)

q2 > 0.5

(11)

R2 > 0.6

(12)

(R2 − R 0 2)/R2 < 0.1

or

0.85 ≤ k ≤ 1.15

0.85 ≤ k′ ≤ 1.15

or

(R2 − R′0 2 )/R2 < 0.1

(13) (14)

where q2 is the internal correlation coefficient q2int for the training set (or external q2ext for the test set); R2 is the squared the correlation coefficient between the predicted and the observed activities; R02 and R′02 are determination coefficients of the predicted vs the observed values and of the observed vs the predicted values, respectively; k and k′ are slopes of regression lines of the predicted vs the observed values and those of the observed vs the predicted values, respectively. The formulas of these criteria can be found in refs 25 and 26 Statistical results evaluated externally with the test set in SI Table 1S were obtained: q2ext = 0.859, R2 = 0.868, R02 = 0.804, R′02 = 0.859, k = 0.976, and k′= 1.007, which suggest that the proposed MLR model satisfies above accepted conditions (eqs 11−14). The calculated values of Ksalt from eq 10 are listed in SI Table 1S and plotted versus experimental values in Figure 1. The rms errors for the training and test sets are 0.0288 and 0.0286, respectively. Descriptor Discussion. As seen from Table 1, all of the significance-values (i.e., Sig.-values or P-values) are less than 0.05, which suggests that the four descriptors all are significant descriptors.6 Moreover, the four descriptors are weakly correlated with each other because all of the values of variance inflation factor (VIF) in this work are less than 10. According to the t-test (in Table 1), the most significant descriptor is the descriptor ALOGP, which is positively correlated with Ksalt. The molecular descriptor ALOGP is Ghose-Crippen octanol−water partition coefficient (log P).27,28 The descriptor ALOGP is related to Setschenow constants Ksalt in accord with previous reports.1,2

where s is the number of input data having nonzero values (of αi and αi*). Thus, it is smaller than the total number of input− output pairs (n). In this study, we used the Gaussian radial basis function (RBF) as the kernel function of SVM models, whose expression is shown as follows: K (xi , xj) = exp( −γ || xi − xj ||2 )

GATS1m

1.00 × 10−4 9.77 × 10−6 0.000 10.247 1.763

where n is the number of samples, R2 is the square of the correlation coefficient, se is the standard error of estimation, F is the Fischer ratio. Then the leave-one-out (LOO) cross-validation procedure was adopted to test the internal predictive ability of the MLR model. The internal correlation coefficient q2int for training set is 0.875. According to refs 25 and 26, a QSPR model is successful if it satisfies the following criteria:

s

∑ (ai − ai*)K (x , y) + b

ET

0.079 0.004 0.000 20.446 3.293

n = 51, R2 = 0.904, se = 0.0303, F = 107.741

where αi and αi* are the introduced Lagrange multipliers given by the solution of a corresponding quadratic optimization problem. Through selecting the appropriate kernel function k(x, y), the entire problem can be solved in the input space itself: f (x ) =

ALOGP

0.082 0.015 0.000 5.355 /

− 0.066Q H

n

∑ (ai − ai*)φ(xi)·φ(x) + b

constant

K salt = 0.068 + 0.079ALOGP + 0.0001E T + 0.115GATS1m

Thus, eq 2 becomes the following form: f (x ) =

descriptor coefficient standard error significance t-test VIF

(9)

The kernel width γ plays a major role in the performance of the kernel and should be carefully tuned to the problem at hand. A small γ value may cause under-fitting, and large one overfitting. The optimal kernel width γ is merely selected based on the trade-off between under-fitting and overfitting of regression models.24 Similarly, too small or too large parameter C value will lead to under-fitting or overfitting, respectively. In ε-SVR models, the three parameters C, ε, and γ should be tuned. All SVM models from the present work were obtained with winSVM (http://www.cs.ucl.ac.uk/staff/M.Sewell/ winsvm/).23

■

RESULTS AND DISCUSSION MLR Models. There is a slightly stronger correlation between Ksalt and Subset II than Subset III for 101 compounds in SI Table 1S, thus the four descriptors in Subset II were selected. Subset II comprises the molecular descriptor ALOGP (Ghose-Crippen octanol−water partition coefficient),16 Geary autocorrelation GATS1m,16 the total energy (ET), and the most positive natural atomic charge on hydrogen atom(QH). These descriptor values are listed in SI Table 1S and their 11184



Article

SVM parameters, the regularization constant C, the Gaussian function parameter γ, and ε in the ε -insensitive loss function, were optimized by performing the grid search.23 The overall performance of SVM was evaluated in terms of the rms error of the test set in SI Table 1S. The curve of rms errors versus each parameter with the other two fixed is shown in Figure 2a−c.

Figure 1. Plot of the experimental versus predicated Ksalt values with the MLR model.

In ref 2, the octanol−water partition coefficient log Kow values were calculated by using ClogP 4.0 software (Bio Byte Corp., Claremont, CA), which adopted another computational method for estimating log P values of organic compounds.28 The correlation coefficient between log Kow and Ksalt for 101 compounds (in SI Table 1S) is 0.7717,2 which is smaller than that in our paper (R = 0.843). This means that there is a stronger correlation between the Setschenow constants and ALOGP than that calculated with respect to log Kow.2 By t-test, the second significant descriptor in Table 1 is the molecular descriptor GATS1m (Geary autocorrelation − lag 1/ weighted by atomic masses). Geary autocorrelation GATSkW is calculated by applying the Geary coefficient to the molecular graph:16 GATSkW =

1 2Δ

nSK

nSK

· ∑i = 1 ∑ j = 1 δij·(Wi − Wj)2 1 (nSK − 1)

nSK

· ∑i = 1 (Wi − W̅ )2

(15)

where Wi is any atomic property such as atomic mass (m), k is the lag varied from 1 to 8, W̅ is its average value on the molecule, nSK is the number of non-hydrogen atoms, and Δ is defined as follows: Δ = Σ dij, where dij is the Kronecker delta (dij = 1 if dij = k, zero otherwise, dij being the topological distance between two considered atoms). It should be noted that GATSkW is the generalized form of GATS1m, and it becomes GATS1m when k = 1 and W = m. The molecular descriptor GATS1m describes how atomic mass is distributed along a topological molecular structure and denotes the polarity of a molecule. Therefore, GATS1m is related to the Setschenow constants Ksalt. The next significant descriptor is the total energy ET. Previous studies have tested that the molecular weight correlates well with Ksalt.3 A large molecule has a more negative descriptor ET and a large molecular weight, which results in a small Ksalt. Thus, ET is positively correlated with Ksalt. The last significant descriptor in Table 1 is QH (the most positive net natural charge on hydrogen atom). According to the classical chemical theory, all chemical reactions are resulted from electrostatic (or orbital) interactions. Atomic charges are used widely in describing molecular polarity.29,30 Therefore, it is easy to understand that the descriptor QH is correlated with Setschenow constants Ksalt. SVM Model. The four descriptors stated above were used to develop SVM models from the 51 compounds in the training set (see SI Table 1S), by applying the winSVM program. Three

Figure 2. (a) The rms error versus C parameter (with γ = 0.5 and ε = 0.01); (b) rms error versus γ parameter (with C = 12 and ε = 0.01); and (c) rms error versus ε parameter (with C = 12 and γ = 0.5).

Figure 2a suggests that the SVR model reaches the best performance with the minimal rms error (0.0279) (with γ = 0.5 and ε = 0.01) and the optimal C is 12. Figure 2b shows the optimal γ corresponding to the minimal error (rms = 0.0279) is set to 0.5 (with C = 12 and ε = 0.01). Figure 2c suggests that the best SVR model (rms = 0.0263) is obtained with the optimal ε being 0.001 (with C = 12 and γ = 0.5). Thus, the optimal SVM parameters, C, γ, and ε, are 12, 0.5 and 0.001, respectively. 11185



Article

09 program, and rebuilt a SVR model using the similar four descriptors: ALOGP, ET, GATSkw, and QH. The rebuilt model (C = 12, γ = 0.5, and ε = 0.001) based on gas-phase descriptors produces the mean rms error of 0.0269 for 101 compounds in SI Table 1S, which is larger than that (rms = 0.0227) in the optimal SVM model based on IEF-PCM. Similar phenomenon can be found in the comparison between MLR models. For example, the test set error of the recalculated MLR model from gas-phase descriptors is 0.0293, which is larger than the corresponding error of 0.0285 of the MLR model based on IEF-PCM. Therefore, it is feasible to calculate descriptors from IEF-PCM and to develop SVR models for Ksalt. Applicability Domain. By showing standardized residuals versus leverages, the Williams plots, Figure 4A,B, were used to

Similarly, the optimal SVM model was validated both internally and externally. Statistical results validated were obtained: q2int = 0.626, q2ext = 0.880, R2 = 0.887, R02 = 0.869, R′02 = 0.827, k = 0.957 and k′ = 0.942, which suggest that the proposed SVM model (C = 12, γ = 0.5, and ε = 0.001) also satisfies above accepted conditions, i.e., eqs 11−14. The calculated values of Ksalt with the optimal SVM model (C = 12, γ = 0.5, and ε = 0.001) are listed in SI Table 1S and plotted versus experimental values in Figure 3.

Figure 3. Plot of the experimental versus predicated Ksalt values with the SVM model.

Model Discussion. Table 2 shows the statistical performance comparison of models in this work with previous models. In this study, the rms errors for the training and test sets of the SVM model are 0.0186 (R = 0.980) and 0.0263 (R = 0.942), respectively. The mean rms error and correlation coefficient for 101 compounds are 0.0227 and 0.965, respectively, which are more accurate than previous models.3,6,7 Moreover, the test set in present SVM model possesses more samples (50 compounds) than that in reference models (30 compounds).3,6,7 These mean our SVM model has better prediction ability than all of the previous models. In addition, our MLR has less mean rms error than previous MLR models.3,6,7 Thus, the optimal descriptors subset (ALOGP, ET, GATSkw, and QH) in this work is superior to the descriptors subsets in previous models.3,6,7 To make a head-to-head comparison between models constructed with descriptors from the IEF-PCM calculations versus those using similar gas-phase descriptors, we also optimized the 101 compounds in a gas-phase with the DFT method at B3LYP level with 6-31G(d) basis set using Gaussian

Figure 4. (A) Williams plot for MLR model with a warning leverage of 0.294. (B) Williams plot for SVM model with a warning leverage of 0.294.

Table 2. Statistical Performance Comparison training set

a

test set

model

algorithm

Pa

rms′b

nc

rms

R

n

rms

R

Zhong model3 Xu model (1)7 Xu model (2)7 Xu model (3)6 Xu model (4)6 Xu model (5)6 Current model (1) Current model (2)

MLR MLR ANN MLR ANN SVM MLR SVM

3 5 5 4 4 4 4 4

0.041 0.0293 0.0265 0.0306 0.0289 0.0288 0.0287 0.0227

71 71 71 71 71 71 51 51

0.042 0.0320 0.0284 0.0319 0.0292 0.0236 0.0288 0.0186

0.887 0.936 0.950 0.932 0.943 0.964 0.951 0.980

30 30 30 30 30 30 50 50

0.040 0.0216 0.0215 0.0273 0.0282 0.0385 0.0286 0.0263

0.872 0.948 0.950 0.946 0.943 0.867 0.932 0.942

P is the number of variables. brms′ is the mean rms for the total set. cn is the number of samples. 11186



Article

21190041) and the Hunan Provincial Natural Science Foundation of China (Contract No. 12JJ6011).

visualize the applicability domain for both the MLR and SVM models. Predictions for only those compounds that fall into this domain may be considered reliable.26,31,32 The leverage of a compound in the original variable space is defined as follows:31,32 hi =

xiT (XT X )−1xi

■

(1) Ni, N.; El-Sayed, M. M.; Sanghvi, T.; Yalkowsky, S. H. Estimation of the Effect of NaCl on the Solubility of Organic Compounds in Aqueous Solutions. J. Pharm. Sci. 2000, 89, 1620−1625. (2) Ni, N.; Yalkowsky, S. H. Prediction of Setschenow Constants. Int. J. Pharm. 2003, 254, 167−172. (3) Li, Y. J.; Hu, Q. H.; Zhong, C. L. Topological Modeling of the Setschenow Constant. Ind. Eng. Chem. Res. 2004, 43, 4465−4468. (4) Setschenow, J. Z. Ü ber Konstitution der Salzlösungen auf Grund ihres Verhaltens zu Kohlensäure. Z. Physik. Chem. 1889, 4, 117−125. (5) Xie, W. H.; Shiu, W. Y.; Mackay, D. A Review of the Effect of Salts on the Solubility of Organic Compounds in Seawater. Marine Environ. Res. 1997, 44, 429−444. (6) Xu, J.; Wang, L.; Wang, L.; Shen, X.; Xu, W. QSPR study of Setschenow Constants of Organic Compounds Using MLR, ANN, and SVM Analyses. J. Comput. Chem. 2011, 32, 3241−3252. (7) Xu, J.; Wang, L.; Wang, L.; Liang, G.; Shen, X.; Xu, W. Prediction of Setschenow Constants of Organic Compounds Based on a 3D Structure Representation. Chemom. Intel. Lab. Sys. 2011, 107, 178− 184. (8) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999−3094. (9) Miertuš, S.; Scrocco, E.; Tomasi, J. Electrostatic Interaction of a Solute with a Continuum. A Direct Utilization of ab Initio Molecular Potentials for the Prevision of Solvent Effects. Chem. Phys. 1981, 55, 117−129. (10) Cancès, M. T.; Mennucci, B.; Tomasi, J. A New Integral Equation Formalism for the Polarizable Continuum Model: Theoretical Background and Applications to Isotropic and Anisotropic Dielectrics. J. Chem. Phys. 1997, 107, 3032−3041. (11) Frediani, L.; Cammi, R.; Corni, S.; Tomasi, J. A Polarizable Continuum Model for Molecules at Diffuse Interfaces. J. Chem. Phys. 2004, 120, 3893−3907. (12) Zhan, C. G.; Dixon, D. A. A Density Functional Theory Approach to the Development of Q-e Parameters for the Prediction of Reactivity in Free-Radical Copolymerizations. J. Phys. Chem. A 2002, 106, 10311−10325. (13) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, Revision A.02; Gaussian Inc.: Wallingford CT, 2009. (14) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785−789. (15) Petersson, G. A.; Al-Laham, M. A. A Complete Basis Set Model Chemistry. II. Open-Shell Systems and the Total Energies of the FirstRow Atoms. J. Chem. Phys. 1991, 94, 6081−6090. (16) Todeschini, R. Consonni, V. Mauri, A. Pavan, M. DRAGON software for the Calculation of Molecular Descriptors, Revision 5.4 for Windows; TALETE srl: Milan, 2006. (17) de Oliveira, D. B.; Gaudio, A. C. BuildQSAR: A New Computer Program for QSAR Analysis. Quant. Struct.-Act. Relat. 2000, 19, 599− 601. (18) Kovatcheva, A.; Golbraikh, A.; Oloff, S.; Xiao, Y. D.; Zheng, W.; Wolschann, P.; Buchbauer, G.; Tropsha, A. Combinatorial QSAR of Ambergris Fragrance Compounds. J. Chem. Inf. Mod. 2004, 44, 582− 595. (19) Camps-Valls, G.; Chalk, A. M.; Serrano-López, A. J.; MartínGuerrero, J. D.; Sonnhammer, E. L. Profiled Support Vector Machines for Antisense Oligonucleotide Efficacy Prediction. BMC Bioinformatics 2004, 5, 135. (20) Tang, H.; Wang, X. S.; Huang, X.-P.; Roth, B. L.; Butler, K. V.; Kozikowski, A. P.; Jung, M.; Tropsha, A. Novel Inhibitors of Human Histone Deacetylase (HDAC) Identified by QSAR Modeling of Known Inhibitors, Virtual Screening, And Experimental Validation. J. Chem. Inf. Model. 2009, 49, 461−476.

(16)

where xi is the descriptor vector of the considered compound, and X is the model matrix derived from the training set descriptor values. The warning leverage h* is defined as follows:31,32 h* = 3 × (p + 1)/n

(17)

where p is the number of variables, n is the number of samples in training set. From Figure 4A,B, one can see that there are two chemicals (1,2,4-trichlorobenzene, No. 1 and cystine, No. 15 in SI Table 1S in the training set having leverage values greater than the warning leverage h* (= 0.294). This suggests that the two compounds show some peculiarity in their structure and derived descriptor values and are poorly represented in the training set. So their presence strongly influences the variable selection and modeling. But their standardized residuals (SR) from both the MLR and SVM models are less than 3, which means that they can stabilize both models. Similarly, the one chemical (lindane, No. 95) with h (= 0.593) > h* and SR (= −0.270 for MLR model and 0.059 for SVM model) < 3 in the test set suggests that both the MLR and SVM models have good generalizability.6,31,32

■

CONCLUSIONS Linear and nonlinear QSAR models based on four molecular descriptors (ALOGP, ET, GATSkw, and QH) were successfully developed for Setschenow constants Ksalt in NaCl solution. The descriptors subset obtained from IEF-PCM calculation, at the DFT/B3LYP level of theory with 6-31G(d) basis set, is superior to the descriptors subsets of molecules in gas-phase. Our SVM model having the mean rms error of 0.0227 for Ksalt of 101 compounds shows better statistical results than all the previous models. In addition, our MLR has less mean rms error than previous MLR models. The feasibility of calculating molecular descriptors from IEF-PCM to predict the Setschenow constants Ksalt of organic compounds in NaCl solution has been demonstrated.

■

ASSOCIATED CONTENT

* Supporting Information S

Table S1 showing the molecular descriptors selected and Ksalt values for 101 compounds. This material is available free of charge via the Internet at http://pubs.acs.org.

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*Phone: +86 731 88822577. Fax: +86 731 88822577. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work has been supported by the Major Program of National Natural Science Foundation of China (Contract No. 11187



Article

(21) Polanski, J.; Bak, A.; Gieleciak, R.; Magdziarz, T. Modeling Robust QSAR. J. Chem. Inf. Model. 2006, 46, 2310−2318. (22) Kumar, R.; Kulkarni, A.; Jayaraman, V. K.; Kulkarni, B. D. Structure−activity Relationships Using Locally Linear Embedding Assisted by Support Vector and Lazy Learning Regressors. Int. Electron. J. Mol. Des. 2004, 3, 118−133. (23) Yu, X. L.; Yi, B.; Wang, X. Y.; Chen, J. F. Predicting Reaction Rate Constants of Ozone with Organic Compounds from Radical Structures. Atmos. Environ. 2012, 51, 124−130. (24) Keerthi, S. S.; Lin, C. J. Asymptotic Behaviors of Support Vector Machines with Gaussian Kernel. Neural Comput. 2003, 15, 1667− 1689. (25) Golbraikh, A.; Tropsha, A. Beware of q(2)! J. Mol. Graphics Modell. 2002, 20, 269−276. (26) Tropsha, A.; Gramatica, P.; Gombar, V. K. The Importance of Being Earnest: Validation Is the Absolute Essential for Successful Application and Interpretation of QSPR Models. QSAR Comb. Sci. 2003, 22, 69−77. (27) Ghose, A. K.; Crippen, G. M. Atomic Physicochemical Parameters for Three-Dimensional Structure-Directed Quantitative Structure-Activity Relationships. 1. Partition Coefficients As a Measure of Hydrophobicity. J. Comput. Chem. 1986, 7, 565−577. (28) Ghose, A. K.; Viswanadhan, V. N.; Wendoloski, J. J. Prediction of Hydrophobic (Lipophilic) Properties of Small Organic Molecules Using Fragmental Methods: An Analysis of ALOGP and CLOGP Methods. J. Phys. Chem. A 1998, 102, 3762−3772. (29) Karelson, M.; Lobanov, V. S.; Katritzky, A. R. Quantum− Chemical Descriptors in QSAR/QSPR Studies. Chem. Rev. 1996, 96, 1027−1043. (30) Yu, X. L.; Yi, B.; Wang, X. Y. Quantitative Structure−Property Relationships for the Reactivity Parameters of Acrylate Monomers. Eur. Polym. J. 2008, 44, 3997−4001. (31) Eriksson, L.; Jaworska, J.; Worth, A. P.; Cronin, M. T. D.; McDowell, R. M.; Gramatica, P. Methods for Reliability and Uncertainty Assessment and for Applicability Evaluations of Classification and Regression-Based QSARs. Environ. Health Perspect. 2003, 111, 1361−1375. (32) Wang, Y. N.; Chen, J. W.; Li, X. H.; Wang, B.; Cai, X. Y.; Huang, L. P. Predicting Rate Constants of Hydroxyl Radical Reactions with Organic Pollutants: Algorithm, Validation, Applicability Domain, And Mechanistic Interpretation. Atmos. Environ. 2009, 43, 1131−1135.

11188


Setschenow Constant Prediction Based on the IEF-PCM Calculations

Recommend Documents