- Research article
- Open Access
Predicting pK a values from EEM atomic charges
Journal of Cheminformatics volume 5, Article number: 18 (2013)
The acid dissociation constant p K a is a very important molecular property, and there is a strong interest in the development of reliable and fast methods for p K a prediction. We have evaluated the p K a prediction capabilities of QSPR models based on empirical atomic charges calculated by the Electronegativity Equalization Method (EEM). Specifically, we collected 18 EEM parameter sets created for 8 different quantum mechanical (QM) charge calculation schemes. Afterwards, we prepared a training set of 74 substituted phenols. Additionally, for each molecule we generated its dissociated form by removing the phenolic hydrogen. For all the molecules in the training set, we then calculated EEM charges using the 18 parameter sets, and the QM charges using the 8 above mentioned charge calculation schemes. For each type of QM and EEM charges, we created one QSPR model employing charges from the non-dissociated molecules (three descriptor QSPR models), and one QSPR model based on charges from both dissociated and non-dissociated molecules (QSPR models with five descriptors). Afterwards, we calculated the quality criteria and evaluated all the QSPR models obtained. We found that QSPR models employing the EEM charges proved as a good approach for the prediction of p K a (63% of these models had R2 > 0.9, while the best had R2 = 0.924). As expected, QM QSPR models provided more accurate p K a predictions than the EEM QSPR models but the differences were not significant. Furthermore, a big advantage of the EEM QSPR models is that their descriptors (i.e., EEM atomic charges) can be calculated markedly faster than the QM charge descriptors. Moreover, we found that the EEM QSPR models are not so strongly influenced by the selection of the charge calculation approach as the QM QSPR models. The robustness of the EEM QSPR models was subsequently confirmed by cross-validation. The applicability of EEM QSPR models for other chemical classes was illustrated by a case study focused on carboxylic acids. In summary, EEM QSPR models constitute a fast and accurate p K a prediction approach that can be used in virtual screening.
The acid dissociation constant p K a is an important molecular property, and its values are of interest in pharmaceutical, chemical, biological and environmental research. The p K a values have found application in many areas, such as the evaluation and optimization of candidate drug molecules [1–3], ADME profiling [4, 5], pharmacokinetics , understanding of protein-ligand interactions [7, 8], etc.. Moreover, the key physicochemical properties like lipophilicity, solubility, and permeability are all p K a dependent. For these reasons, p K a values are important for virtual screening. Therefore, both the research community and pharmaceutical companies are interested in the development of reliable and above all fast methods for p K a prediction.
Several approaches for p K a prediction have been developed [8–11], namely LFER (Linear Free Energy Relationships) methods [12, 13], database methods, decision tree methods , ab initio quantum mechanical calculations [15, 16], ANN (artificial neural networks) methods  or QSPR (quantitative structure-property relationship) modelling [18–20]. However, p K a values remain one of the most challenging physicochemical properties to predict. A promising approach for p K a prediction is to use QSPR models which employ partial atomic charges as descriptors [21–24].
The partial atomic charges cannot be determined experimentally and they are also not quantum mechanical observables. For this reason, the rules for determining partial atomic charges depend on their application (e.g. molecular mechanics energy, p K a etc.), and many different methods have been developed for their calculation. The charge calculation methods can be divided into two main groups, namely quantum mechanical (QM) approaches and empirical approaches.
The quantum mechanical approaches first calculate a molecular wave function by a combination of some theory level (e.g., HF, B3LYP, MP2) and basis set (e.g., STO-3G, 6–31G*), and then partition this wave function among the atoms (i.e., the assignment of a specific part of the molecular electron density to each atom). This partitioning can be done using an orbital-based population analysis, such as MPA (Mulliken population analysis) [25, 26], Löwdin population analysis  or NPA (natural population analysis) . Other partitioning approaches are based on a wave-function-dependent physical observable. Such approaches are, for example, AIM (atoms in molecules) , Hirshfeld population analysis  and electrostatic potential fitting methods like CHELPG  or MK (Merz-Singh-Kollman) . Another partitioning method is the mapping of QM atomic charges to reproduce charge-dependent observables (e.g., CM1, CM2, CM3 and CM4) .
Empirical approaches determine partial atomic charges without calculating a quantum mechanical wave function for the given molecule. Therefore they are markedly faster than QM approaches. One of the first empirical approaches developed, CHARGE , performs a breakdown of the charge transmission by polar atoms into one-bond, two-bond, and three-bond additive contributions. Most of the other empirical approaches have been derived on the basis of the electronegativity equalization principle. One group of these empirical approaches invoke the Laplacian matrix formalism, and result in a redistribution of electronegativity. Such methods are PEOE (partial equalization of orbital electronegativity) , GDAC (geometry-dependent atomic charge) , KCM (Kirchhoff charge model) , DENR (dynamic electronegativity relaxation)  or TSEF (topologically symmetric energy function) . The second group of approaches use full equalization of orbital electronegativity, and such approaches are, for example, EEM (electronegativity equalization method) , QEq (charge equilibration)  or SQE (split charge equilibration) . The empirical atomic charge calculation approaches can also be divided into ’topological’ and ’geometrical’. Topological charges are calculated using the 2D structure of the molecule, and they are conformationally independent (i.e., CHARGE, PEOE, KCM, DENR, and TSEF). Geometrical charges are computed from the 3D structure of the molecule and they consider the influence of conformation (i.e., GDAC, EEM, Qeq, and SQE).
The prediction of p K a using QSPR models which employ QM atomic charges was described in several studies [21–24], which have analyzed the precision of this approach and compared the quality of QSPR models based on different QM charge calculation schemes. All these studies show that QM charges are successful descriptors for p K a prediction, as the QSPR models based on QM atomic charges are able to calculate p K a with high accuracy. The weak point of QM charges is that their calculation is very slow, as the computational complexity is at least θ (E4), where E is the number of electrons in the molecule. Therefore, p K a prediction by QSPR models based on QM charges cannot be applied in virtual screening, as it is not feasible to compute QM atomic charges for hundreds of thousands of compounds in a reasonable time. This issue can be avoided if empirical charges are used instead of QM charges. A few studies were published, which give QSPR models for predicting p K a using topological empirical charges as descriptors (specifically PEOE charges) [22, 42, 43]. But these models provided relatively weak predictions.
The geometrical charges seem to be more promissing descriptors, because they are able to take into consideration the influence of the molecule’s conformation on the atomic charges. The conformation of the atoms surrounding the dissociating hydrogens strongly influences the dissociation process, and also the atomic charges.
The EEM method is a geometrical empirical charge calculation approach which can be useful for p K a prediction by QSPR. This approach calculates charges using the following equation system:
where q i is the charge of atom i; Ri,jis the distance between atoms i and j; Q is the total charge of the molecule; N is the number of atoms in the molecule; is the molecular electronegativity, and A i , B i and κ are empirical parameters. These parameters are obtained by a parameterization process, which uses QM atomic charges to calculate a set of parameters for which EEM best reproduces these QM charges. EEM is very popular, and despite the fact that it was developed more than twenty years ago, new parameterizations [39, 44–50] and modifications [47, 51, 52] of EEM are still under development. Its accuracy is comparable to the QM charge calculation approach for which it was parameterized. Additionally, EEM is very fast, as its computational complexity is θ (N3), where N is the number of atoms in the molecule.
Therefore, in the present study, we focus on p K a prediction using QSPR models which employ EEM charges. Specifically, we created and evaluated QSPR models based on EEM charges computed using 18 EEM parameter sets. We also compared these QSPR models with corresponding QSPR models which employ QM charges computed by the same charge calculation schemes used for EEM parameterization.
EEM parameter sets
In our study, we used all EEM parameters published till now. Specifically, we found 18 different EEM parameters sets, published in 8 different articles [39, 44–50]. The parameters cover two QM theory levels (HF and B3LYP), two basis sets (STO-3G and 6–31G*) and six population analyses (MPA, NPA, Hirshfeld, MK, CHELPG, AIM). Unfortunately, only some combinations of QM theory levels, basis sets and population analyses are available. On the other hand, more parameter sets were published for some combinations (i.e., 6 parameter sets for HF/STO-3G/MPA). All the parameter sets include parameters for C, O, N and H. Some sets include also parameters for S, P, halogens and metals. Most of the sets do not include parameters for C and N bonded by triple bond. Summary information about all these parameter sets is given in Table 1.
EEM charge calculation
The EEM charges were calculated by the program EEM SOLVER  using each of the 18 EEM parameter sets.
QM charge calculation
We calculated QM atomic charges for all the combinations of QM theory level, basis set and population analysis for which we have EEM parameters (see Table 1). Specifically, atomic charges were calculated for these eight QM approaches: HF/STO-3G/MPA, HF/6–31G*/MK, and B3LYP/6–31G* with MPA, NPA, Hirshfeld, MK, CHELPG and AIM). The QM charge calculations were carried out using Gaussian09 . In the case of AIM population analysis, the output from Gaussian09 was further processed by the software package AIMAll .
Data set for phenols
There are two main ways to create a QSPR model for a feature to be predicted. The first is to create as general a model as possible, with the risk that the accuracy of such a model may not be high. The second approach is to develop more models, each of them being dedicated to a certain class of compounds. Here we took the second approach, following a similar methodology as in previous studies [21–24]. Specifically, we focus on substituted phenols, because they are the most common test set molecules employed in the evaluation of novel p K a prediction approaches [21–24, 56–58]. Our data set contains the 3D structures of 74 distinct phenol molecules. This data set is of high structural diversity and it covers molecules with p K a values from 0.38 to 11.1. The molecules were obtained from the NCI Open Database Compounds  and their 3D structures were generated by CORINA 2.6 , without any further geometry optimization. Our data set is a subset of the phenol data set used in our previous work related to p K a prediction from QM atomic charges . The subset is made up of phenols which contain only C, O, N and H, and none of the molecules contain triple bonds. This limitation is necessary, because the EEM parameters of all 18 studied EEM parameter sets are available only for such molecules (see Table 1). For each phenol molecule from our data set, we also prepared the structure of the dissociated form, where the hydrogen is missing from the phenolic OH group. This dissociated molecule was created by removing the hydrogen from the original structure without subsequent geometry optimization. The list of the molecules, including their names, NCS numbers, CAS numbers and experimental p K a values, can be found in the (Additional file 1: Table S1a). The SDF files with the 3D structures of molecules and their dissociated forms are also in the (Additional file 2: Molecules).
Data set for carboxylic acids
An aspect which is very important for the applicability of the p K a prediction approach is its transferability to other chemical classes. In this work, we provide a case study showing the performance of the approach on carboxylic acids, which are also very common testing molecules for p K a prediction approaches [19–21, 43]. The data set contains 71 distinct molecules of carboxylic acids and their dissociated forms. The 3D structures of these molecules were obtained in the same way as for the phenols. The list of the molecules, including their names, NCS numbers, CAS numbers and experimental p K a values can be found in the (Additional file 3: Table S1b). The SDF files with the 3D structures of the molecules and their dissociated forms are also included in the (Additional file 2: Molecules).
p K a values
The experimental p K a values were taken from the Physprop database .
Descriptors and QSPR models for phenols
Our descriptors were atomic charges. We analyzed two types of QSPR models, namely QSPR models with three descriptors (3d QSPR models) and QSPR models with five descriptors (5d QSPR models).
The 3d QSPR models used those descriptors which proved to be the most relevant for p K a prediction in our previous study . Therefore these descriptors were the atomic charge of the hydrogen atom from the phenolic OH group (q H ), the charge on the oxygen atom from the phenolic OH group (q O ), and the charge on the carbon atom binding the phenolic OH group (qC1). These descriptors were used to establish the QSPR models by the general equation:
where p H , p O , pC1 and p are parameters of the QSPR model (i.e., constants derived by multiple linear regression). The 5d QSPR models employ the above mentioned descriptors q H , q O and qC1 and additionally also the charge on the phenoxide O- from the dissociated molecule (q O D ), and the charge on the carbon atom binding this oxygen (qC1D). Using the charges from the dissociated molecules for p K a prediction was inspired by the work of Dixon et al. . The equation of the 5d QSPR models is therefore:
where , , , , and p′ are parameters of the QSPR model.
Descriptors and QSPR models for carboxylic acids
The descriptors were again atomic charges and, similarly as for phenols, two types of QSPR models were developed and evaluated. Specifically, QSPR models with four descriptors (4d QSPR models) and QSPR models with seven descriptors (7d QSPR models). The 4d QSPR models used similar descriptors as the 3d models for phenols - the atomic charge of the hydrogen atom from the COOH group (q H ), the charge on the hydrogen bound oxygen atom from the COOH group (q O ), and the charge on the carbon atom binding the COOH group (qC1). Additionally, also the charge of the second carboxyl oxygen (qO2) is included. These 4d QSPR models are represented by the equation:
where p H , p O , pO2, pC1 and p are parameters of the QSPR model. The 7d QSPR models employ also charges from the dissociated forms, namely the charge on the carboxyl oxygens (q O D , qO2D) and the charge on the carboxylic carbon atom (qC1D). The equation of the 7d QSPR models is therefore:
where , , , , , , and p′ are parameters of the QSPR model.
QSPR model parameterization
The parameterization of the QSPR models was done by multiple linear regression (MLR) using the software tool QSPR Designer .
Results and discussion
QM and EEM QSPR models for phenols
We prepared one 3d QSPR model and one 5d QSPR model using atomic charges calculated by each of the above mentioned 18 EEM parameter sets. These models are denoted 3d or 5d EEM QSPR models. Additionally, we created one 3d and one 5d QSPR model using atomic charges calculated by each of the corresponding 8 QM charge calculation approaches (denoted as 3d or 5d QM QSPR models). The data set of 74 phenol molecules was used for the parameterization of the QSPR models, and the obtained models were validated for all molecules in the data set.
The parameterization of the 3d EEM QSPR models showed that several molecules in the data set perform as outliers. For this reason, we created also EEM QSPR models without outliers (i.e., EEM QSPR models for which parameterization was done using a data set that excluded the previously observed outliers). These models are denoted 3d EEM QSPR WO models. We classified as outliers 10% of the molecules from our data set, which had the highest Cook’s square distance. Therefore the 3d EEM QSPR WO models were parameterized using 67 molecules, and their validation was also done on the data set excluding outliers.
The quality of the QSPR models, i.e. the correlation between experimental p K a and the p K a calculated by each model, was evaluated using the squared Pearson correlation coefficient (R2), root mean square error (RMSE), and average absolute p K a error (), while the statistical criteria were the standard deviation of the estimation (s) and Fisher’s statistics of the regression (F).
Table 2 contains the quality criteria (R2, RMSE, ) and statistical criteria (s and F) for all the QSPR models analyzed. All these models are statistically significant at p = 0.01. Since our data sets contained 74 and 67 molecules, respectively, the appropriate F value to consider was that for 60 samples. Thus, the 3d QSPR models are statistically significant (at p= 0.01) when F > 4.126 and the 5d QSPR models when F > 3.339. Figure 1 summarizes the R2 of all QSPR models for ease of visual comparison, and Tables 3 and 4 provide a comparison of the models from specific points of view. The parameters of the QSPR models are summarized in the (Additional file 4: Table S2) and all charge descriptors and p K a values are contained in the (Additional file 5: Table S6). The most relevant graphs of correlation between experimental and calculated p K a are visualized in Figure 2.
Prediction of p K a using EEM charges
The key question we wanted to answer in this paper is whether EEM QSPR models are applicable for p K a prediction. For this purpose we selected a set of phenol molecules and generated QSPR models which used EEM atomic charges as descriptors. We then evaluated the accuracy of these models by comparing the predicted p K a values with the experimental ones. The results (see Tables 2 and 3, Figure 1) clearly show that QSPR models based on EEM charges are indeed able to predict the p K a of phenols with very good accuracy. Namely, 63% of the EEM QSPR models evaluated in this study were able to predict p K a with R2 > 0.9. The average R2 for all 54 EEM QSPR models considered was 0.9, while the best EEM QSPR model reached R2 = 0.924. Our findings thus suggest that EEM atomic charges may prove as efficient QSPR descriptors for p K a prediction. The only drawback of EEM is that EEM parameters are currently not available for some types of atoms. Nevertheless, EEM parameterization is still a topic of research, therefore more general parameter sets are being developed.
Improvement of EEM QSPR models by removing outliers
The quality of 3d EEM QSPR models can be markedly increased by removing the outliers. In this case, the models have average R2 = 0.911 and 83% of them have R2 > 0.9. The disadvantage of these models is that they are not able to cover the complete data set (i.e., 10% of molecules must be excluded as outliers).
On the other hand, the outliers are similar for all EEM QSPR models. For example, while 16 molecules from our data set are outliers for at least one parameter set, 10 out of these 16 molecules are outliers for five or more parameter sets. From the chemical point of view, most of the outliers contain one or more nitro groups. This may be related to reported lower accuracy of EEM for these groups . In general one limitation of the 3d EEM QSPR models is that they are very sensitive to the quality of EEM charges. Therefore, if the EEM charges are inaccurate for certain compounds or class of compounds, the 3d QSPR models based on these EEM charges will have lower performance for these compounds or class of compounds. In addition, a lower experimental accuracy of these p K a values may also be a reason for low performance in some cases. A table containing information about outlier molecules is given in the (Additional file 6: Table S3).
Improvement of EEM QSPR models by adding descriptors
Our first EEM QSPR models contained three descriptors (3d), namely atomic charges originating from the non-dissociated molecule. Nonetheless, in our study we found that using two additional charge descriptors from the dissociated molecule can markedly improve the predictive power of the EEM QSPR models. Tables 2 and 3, Figure 1 show that these new 5d EEM QSPR models provide better p K a prediction than their corresponding 3d EEM QSPR models. Specifically, adding the descriptors derived from the dissociated molecules increased the average R2 value for the EEM QSPR models from 0.876 to 0.913.
Comparison of EEM QSPR models and QM QSPR models
Another important question is how accurate the EEM QSPR models are in comparison with QM QSPR models. Table 2 and Figure 1 show that QM QSPR models provide, in most cases, more precise predictions. This is confirmed also by the average R2 values from Table 3. This is not surprising, since EEM is an empirical method which just mimics the QM approach for which it was parameterized. An interesting fact is that the differences in accuracy between QM QSPR models and EEM QSPR models are not substantial. For example, 5d EEM QSPR models have average R2 = 0.913, while 5d QM QSPR models have average R2 = 0.951. We also note that adding more descriptors to a QM QSPR model brings less improvement than adding more descriptors to an EEM QSPR model.
Influence of theory level and basis set
EEM parameters are available only for a relatively small number of theory levels (HF and B3LYP) and basis sets (STO-3G and 6–31G*). Therefore we can not perform such a deep analysis of theory level and basis set influence on p K a prediction capability from EEM atomic charges, as was done for QM QSPR models by Gross et al.  or Svobodova et al. . We can only compare the models employing HF/STO-3G and B3LYP/6–31G* charges, as these are the only combinations for which EEM parameters are available for the same population analysis (MPA). Therefore we can study only the influence of the combination of theory level / basis set, and not the isolated influence of the theory level or basis set. Our analysis revealed that B3LYP/6–31G* charges provide slightly more accurate QM QSPR models than HF/STO-3G charges (see Table 4). This is in agreement with our previous findings , and it can be explained by the fact that 6–31G* is a more robust basis set than STO-3G. However, the difference is not marked in the case of EEM QSPR models.
Influence of population analysis
Eleven EEM parameter sets were published for B3LYP/6–31G* with six different population analyses (see Table 1). Therefore it is straightforward to analyze the influence of the population analysis on the predictive power of the resulting QSPR models (see Table 4). We found that MPA and NPA provide the best QM models, while MK and CHELPG (PAs based on fitting the atomic charges to the molecular electrostatic potential) provide weak QM models. Our results are in agreement with previous studies [22, 24]. QM QSPR models based on AIM predict p K a with accuracy comparable to MPA and NPA. In the case of EEM QSPR models, we did indeed find that MPA provided the best models, but most of the other population analyses gave comparable results. This confirms previous observations that the EEM approach is not able to faithfully mimic MK charges . On the other hand, this drawback of EEM allowed the EEM QSPR models employing MK charges to predict p K a more accurately than the corresponding QM QSPR models.
Influence of the EEM parameter set
Two or more EEM parameter sets are available in literature for four combinations of theory level, basis set and population analysis (see Table 1). We found that the quality of EEM QSPR models employing the same types of charges slightly varies when using EEM parameters coming from different studies (see Table 2 and Figure 1). Even EEM parameters from the same study, but obtained by different approaches, lead to QSPR models of slightly different quality. In any case, these differences are minimal.
Comparison with previous work
QM QSPR models for p K a prediction in phenols, similar to those presented in this paper (i.e., employing similar charges) were previously published by Gross and Seybold , Kreye and Seybold  and Svobodova and Geidl . Table 5 shows a comparison between these models and the models developed in this study. Our work is the first which presents QSPR models for p K a prediction based on EEM charges. Therefore, we can not provide a comparison between EEM QSPR models, but we can compare against QSPR models based on QM charges only. It is seen therein that our 3d QM QSPR models show markedly higher R2 and F values than the models published by Gross and Seybold and Kreye and Seybold (even if some of these models employ higher basis sets) and comparable R2 and F values as models published by Svobodova and Geidl. Moreover, our 5d QM QSPR models outperform the models from Svobodova and Geidl. Our best EEM QSPR models (i.e., 5d EEM QSPR models) provide even better results than QM QSPR models from Gross and Seybold and Kreye and Seybold. These EEM QSPR models are not as accurate as the QM QSPR models published by Svobodova and Geidl or those developed in this work, but the loss of accuracy is not too high (R2 values are still > 0.91).
Our results show that 5d EEM QSPR models provide a fast and accurate approach for p K a prediction. Nonetheless, the robustness of these models should be proved. Therefore, all the 5d EEM QSPR models (i.e., 18 models) were tested by cross-validation. For comparison, also the cross-validation of all 5d QM QSPR models (i.e., 8 models) was done. The k-fold cross-validation procedure was used [64, 65], where k = 5. Specifically, the set of phenol molecules was divided into five parts (each contained 20% of the molecules). The division was done randomly, and included stratification by p K a value. Afterwards, five cross validation steps were performed. In the first step, the first part was selected as a test set, and the remaining four parts were taken together as the training set. The test and training sets for the other steps were prepared in a similar manner, by subsequently considering one part as a test set, while the remaining parts served as a training set. For each step, the QSPR model was parameterized on the training set. Afterwards, the p K a values of the respective test molecules were calculated via this model, and compared with experimental p K a values. The results are summarized in the (Additional file 7: Table S4), while the cross-validation results for the best and the worst performing 5d EEM QSPR models are shown in Table 6. The cross-validation showed that the models are stable and the values of R2 and RMSE are similar for the test set, the training set and the complete set. The robustness of EEM QSPR models and QM QSPR models is comparable.
Case study for carboxylic acids
We have shown that QSPR models based on EEM atomic charges can be used for predicting p K a in phenols. In order to evaluate the general applicability of this approach for p K a prediction, we tested the performance of such models for carboxylic acids. This case study is done for the charge schemes found to provide the best QM and EEM QSPR models in the case of phenols. Specifically, QM charges calculated by HF/STO-3G/MPA, B3LYP/6–31G*/MPA and B3LYP/6–31G*/NPA, and EEM charges calculated using the corresponding EEM parameters. Because 5d QSPR models provide the most accurate prediction for phenols, the case study is focused on their analogue for carboxylic acids, i.e., 7d QSPR models. Squared Pearson correlation coefficients of the analysed QSPR models are summarized in Figure 3, and all the quality and statistical criteria can be found in (Additional file 8: Table S5). The results show that 7d EEM QSPR models are able to predict the p K a of carboxylic acids with very good accuracy. Namely, 5 out of 12 analysed 7d EEM QSPR models were able to predict p K a with R2 > 0.9, while the best EEM QSPR model reached R2 = 0.925. Therefore, we concluded that EEM QSPR models are indeed applicable also for carboxylic acids. Again QM QSPR models perform better than EEM QSPR models, but the differences are not substantial.
We found that the QSPR models employing EEM charges can be a suitable approach for p K a prediction. From our 54 EEM QSPR models focused on phenols, 63% show a correlation of R2 > 0.9 between the experimental and predicted p K a . The most successful type of these EEM QSPR models employed 5 descriptors, namely the atomic charge of the hydrogen atom from the phenolic OH group, the charge on the oxygen atom from the phenolic OH group, the charge on the carbon atom binding the phenolic OH group, the charge on the oxygen from the phenoxide O- from the dissociated molecule, and the charge on the carbon atom binding this oxygen. Specifically, 94% of these models have R2 > 0.9, and the best one has R2 = 0.920. In general, including charge descriptors from dissociated molecules, which was introduced in our work, always increases the quality of a QSPR model. The only drawback of EEM QSPR models is that the EEM parameters are currently not available for all types of atoms. Therefore the EEM parameter sets need to be expanded to larger sets of molecules and further improved.
As expected, the QM QSPR models provided more accurate p K a predictions than the EEM QSPR models. Nevertheless, these differences are not substantial. Furthermore, a big advantage of EEM QSPR models is that one can calculate the EEM charges markedly faster than the QM charges. Moreover, the EEM QSPR models are not so strongly influenced by the charge calculation approach as the QM QSPR models are. Specifically, the QM QSPR models which use atomic charges obtained from calculations with higher basis set perform better, while the EEM QSPR models do not show such marked differences. Similarly, the quality of QM QSPR models depends a lot on population analysis, but EEM QSPR models are not influenced so much. Namely, QM QSPR models which use atomic charges calculated from MPA, NPA and Hirshfeld PA performed very well, while MK provides only weak models. In the case of EEM QSPR models, MPA performs also the best, but all other PAs (including MK) provide accurate results as well. The source of the EEM parameters also did not affect the quality of the EEM QSPR models significantly.
The robustness of EEM QSPR models was successfully confirmed by cross-validation. Specifically, the accuracy of p K a prediction for the test, training and complete set were comparable. The applicability of EEM QSPR models for other chemical classes was tested in a case study focused on carboxylic acids. This case study showed that EEM QSPR models are indeed applicable for p K a prediction also for carboxylic acids. Namely, 5 from 12 of these models were able to predict p K a with R2 > 0.9, while the best EEM QSPR model reached R2 = 0.925.
Therefore, EEM QSPR models constitute a very promising approach for the prediction of p K a . Their main advantages are that they are accurate, and can predict p K a values very quickly, since the atomic charge descriptors used in the QSPR model can be obtained much faster by EEM than by QM. Additionally, the quality of EEM QSPR models is less dependent on the type of atomic charges used (theory level, basis set, population analysis) than in the case of QM QSPR models. Accordingly, EEM QSPR models constitute a p K a prediction approach which is very suitable for virtual screening.
The concept of the study originated from JK and was reviewed and extended by RA, while the design was put together by RSV and SG and reviewed by JK and RA. SG and CMI collected and prepared the input data. SG, OS, DS and TB performed the acquisition and post-processing of data. The data were analyzed and interpreted by RSV, SG, CMI and JK. The manuscript was written by RSV and SG in cooperation with JK and CMI, and reviewed by all authors.
Radka Svobodová Vařeková and Stanislav Geidl wish it to be known that, in their opinion, the first two authors should be regarded as joint First Authors.
Atoms in Molecules
Artificial Neural Networks
Becke, three-parameter, Lee-Yang-Parr
Dynamic Electronegativity Relaxation
Electronegativity Equalization Method
Geometry-Dependent Atomic Charge
Kirchhoff Charge Model
Linear Free Energy Relationships
Multiple Linear Regression
Møller-Plesset Perturbation Theory
Mulliken Population Analysis
Natural Population Analysis
Partial Equalization of Orbital Electronegativity
Quantitative Structure-Property Relationship
Root Mean Square Error
Split Charge Equilibration
Topologically Symmetric Energy Function.
Ishihama Y, Nakamura M, Miwa T, Kajima T, Asakawa N: A rapid method for pKa determination of drugs using pressure-assisted capillary electrophoresis with photodiode array detection in drug discovery. J Pharm Sci. 2002, 91 (4): 933-942. 10.1002/jps.10087.
Babić S: Determination of pKa values of active pharmaceutical ingredients. TrAC. 2007, 26 (11): 1043-1061.
Manallack D: The pKKa distribution of drugs: application to drug discovery. Perspect Med Chem. 2007, 1: 25-38.
Wan H, Ulander J: High-throughput pKa screening and prediction amenable for ADME profiling. Expert Opin Drug Metabx Toxicol. 2006, 2: 139-155. 10.1517/17425255.2.1.139.
Cruciani G, Milletti F, Storchi L, Sforna G, Goracci L: In silico pKa prediction and ADME profiling. Chem Biodivers. 2009, 6 (11): 1812-1821. 10.1002/cbdv.200900153.
Comer J, Tam K: Pharmacokinetic Optimization in Drug Research: Biological, Physicochemical, and Computational Strategies. Switzerland: Wiley-VCH, Verlag Helvetica Chimica Acta, Postfach, CH-8042 Zürich; 2001
Klebe G: Recent developments in structure-based drug design. J Mol Med. 2000, 78: 269-281. 10.1007/s001090000084.
Lee AC, Crippen GM: Predicting pKa. J Chem Inf Model. 2009, 49: 2013-2033. 10.1021/ci900209w.
Rupp M, Körner R, Tetko IV: Predicting the pKKa of small molecules. Comb Chem High Throughput Screen. 2010, 14 (5): 307-327.
Fraczkiewicz R: In Silico Prediction of Ionization, Volume 5. 2006, Oxford: Elsevier
Ho J, Coote M: A universal approach for continuum solvent pKa calculations: Are we there yet?. Theor Chim Acta. 2010, 125 (1–2): 3-21.
Clark J, Perrin DD: Prediction of the strengths of organic bases. Q ReV Chem Soc. 1964, 18: 295-320. 10.1039/qr9641800295.
Perrin DD, Dempsey B, Serjeant EP: pKa Prediction for Organic Acids and Bases. 1981, New York: Chapman and Hall
Blower PE, Cross KP: Decision tree methods in pharmaceutical research. Curr Top Med Chem. 2006, 6: 31-39. 10.2174/156802606775193301.
Liptak MD, Gross KC, Seybold PG, Feldgus S, Shields G: Absolute pKa determinations for substituted phenols. J Am Chem Soc. 2002, 124: 6421-6427. 10.1021/ja012474j.
Toth AM, Liptak MD, Phillips DL, Shields GC: Accurate relative pKa calculations for carboxylic acids using complete basis set and Gaussian-n models combined with continuum solvation methods. J Chem Phys. 2001, 114: 4595-4606. 10.1063/1.1337862.
Hagan MT, Demuth HB, Beale M: In Neural, Network Design. 1996, Boston: PWS, MA
Jelfs S, Ertl P, Selzer P: Estimation of pKa for druglike compounds using semiempirical and information-based descriptors. J Chem Inf Model. 2007, 47: 450-459. 10.1021/ci600285n.
Dixon SL, Jurs PC: Estimation of pKa for organic oxyacids using calculated atomic charges. J Comput Chem. 1993, 14: 1460-1467. 10.1002/jcc.540141208.
Zhang J, Kleinöder T, Gasteiger J: Prediction of pKa values for aliphatic carboxylic acids and alcohols with empirical atomic charge descriptors. J Chem Inf Model. 2006, 46: 2256-2266. 10.1021/ci060129d.
Citra MJ: Estimating the pKa of phenols, carboxylic acids and alcohols from semi-empirical quantum chemical methods. Chemosphere. 1999, 1: 191-206.
Gross KC, Seybold PG, Hadad CM: Comparison of different atomic charge schemes for predicting pKa variations in substituted anilines and phenols. Int J Quantum Chem. 2002, 90: 445-458. 10.1002/qua.10108.
Kreye WC, Seybold PG: Correlations between quantum chemical indices and the pKas of a diverse set of organic phenols. Int J Quantum Chem. 2009, 109: 3679-3684. 10.1002/qua.22343.
Svobodová Vařeková R, Geidl S, Ionescu CM, Skřehota O, Kudera M, Sehnal D, Bouchal T, Abagyan R, Huber HJ, Koča J: Predicting pKa values of substituted phenols from atomic charges: Comparison of different quantum mechanical methods and charge distribution schemes. J Chem Inf Model. 2011, 51 (8): 1795-1806. 10.1021/ci200133w.
Mulliken RS: Electronic structures of molecules XI. Electroaffinity, molecular orbitals and dipole moments. J Chem Phys. 1935, 3 (9): 573-585. 10.1063/1.1749731.
Mulliken RS: Criteria for construction of good self-consistent-field molecular orbital wave functions, and significance of LCAO-MO population analysis. J Chem Phys. 1962, 36 (12): 3428-3439. 10.1063/1.1732476.
Lowdin PO: On the non-orthogonality problem connected with the use of atomic wave functions in the theory of molecules and crystals. J Chem Phys. 1950, 18 (3): 365-375. 10.1063/1.1747632.
Reed AE, Weinstock RB, Weinhold F: Natural-population analysis. J Chem Phys. 1985, 83 (2): 735-746. 10.1063/1.449486.
Bader RFW, Larouche A, Gatti C, Carroll MT, Macdougall PJ, Wiberg KB: Properties of atoms in molecules - dipole-moments and transferability of properties. J Chem Phys. 1987, 87 (2): 1142-1152. 10.1063/1.453294.
Hirshfeld FL: Bonded-atom fragments for describing molecular charge-densities. Theor Chim Acta. 1977, 44 (2): 129-138. 10.1007/BF00549096.
Breneman CM, Wiberg KB: Determining atom-centered monopoles from molecular electrostatic potentials - the need for high sampling density in formamide conformational-analysis. J Comput Chem. 1990, 11 (3): 361-373. 10.1002/jcc.540110311.
Besler BH, Merz KM, Kollman PA: Atomic charges derived from semiempirical methods. J Comput Chem. 1990, 11 (4): 431-439. 10.1002/jcc.540110404.
Kelly CP, Cramer CJ, Truhlar DG: Accurate partial atomic charges for high-energy molecules using class IV charge models with the MIDI! basis set. Theor Chem Acc. 2005, 113 (3): 133-151. 10.1007/s00214-004-0624-x.
Abraham RJ, Griffiths L, Loftus P: Approaches to charge calculations in molecular mechanics. J Comput Chem. 1982, 3 (3): 407-416. 10.1002/jcc.540030316.
Gasteiger J, Marsili M: Iterative partial equalization of orbital electronegativity - a rapid access to atomic charges. Tetrahedron. 1980, 36 (22): 3219-3228. 10.1016/0040-4020(80)80168-2.
Cho KH, Kang YK, No KT, Scheraga HA: A fast method for calculating geometry-dependent net atomic charges for polypeptides. J Phys Chem B. 2001, 105 (17): 3624-3634. 10.1021/jp0023213.
Oliferenko AA, Pisarev SA, Palyulin VA, Zefirov NS: Atomic charges via electronegativity equalization: Generalizations and perspectives. Adv Quantum Chem. 2006, 51: 139-156.
Shulga DA, Oliferenko AA, Pisarev SA, Palyulin VA, Zefirov NS: Parameterization of empirical schemes of partial atomic charge calculation for reproducing the molecular electrostatic potential. Dokl Chem. 2008, 419: 57-61. 10.1134/S001250080803004X.
Mortier WJ, Ghosh SK, Shankar S: Electronegativity equalization method for the calculation of atomic charges in molecules. J Am Chem Soc. 1986, 108: 4315-4320. 10.1021/ja00275a013.
Rappe AK, Goddard WA: Charge equilibration for molecular-dynamics simulations. J Phys Chem. 1991, 95 (8): 3358-3363. 10.1021/j100161a070.
Nistor RA, Polihronov JG, Muser MH, Mosey NJ: A generalization of the charge equilibration method for nonmetallic materials. J Chem Phys. 2006, 125 (9): 094108-094118. 10.1063/1.2346671.
Czodrowski P, Dramburg I, Sotriffer CA, Klebe G: Development, validation, and application of adapted PEOE charges to estimate pKa values of functional groups in protein–ligand complexes. Proteins Struct Funct Bioinf. 2006, 65: 424-437. 10.1002/prot.21110.
Gieleciak R, Polanski J: Modeling robust QSAR. 2. Iterative variable elimination schemes for CoMSA: Application for modeling benzoic acid pKa values. J Chem Inf Model. 2007, 47: 547-556. 10.1021/ci600295z.
Svobodová Vařeková R, Jiroušková Z, Vaněk J, Suchomel S, Koča J: Electronegativity equalization method: Parameterization and validation for large sets of organic, organohalogene and organometal molecule. Int J Mol Sci. 2007, 8: 572-582. 10.3390/i8070572.
Baekelandt BG, Mortier WJ, Lievens JL, Schoonheydt RA: Probing the reactivity of different sites within a molecule or solid by direct computation of molecular sensitivities via an extension of the electronegativity equalization method. J Am Chem Soc. 1991, 113 (18): 6730-6734. 10.1021/ja00018a003.
Jiroušková Z, Svobodová Vařeková R, Vaněk J, Koča J: Electronegativity equalization method: Parameterization and validation for organic molecules using the Merz–Kollman–Singh charge distribution scheme. J Comput Chem. 2009, 30: 1174-1178. 10.1002/jcc.21142.
Chaves J, Barroso JM, Bultinck P, Carbo-Dorca R: Toward an alternative hardness kernel matrix structure in the Electronegativity Equalization Method (EEM). J Chem Inf Model. 2006, 46 (4): 1657-1665. 10.1021/ci050505e.
Bultinck P, Langenaeker W, Lahorte P, De Proft, Geerlings P, Waroquier M, Tollenaere J: The electronegativity equalization method I: Parametrization and validation for atomic charge calculations. J Phys Chem A. 2002, 106 (34): 7887-7894. 10.1021/jp0205463.
Ouyang Y, Ye F, Liang Y: A modified electronegativity equalization method for fast and accurate calculation of atomic charges in large biological molecules. Phys Chem. 2009, 11: 6082-6089.
Bultinck P, Vanholme R, Popelier PLA, De Proft, Geerlings P: High-speed calculation of AIM charges through the electronegativity equalization method. J Phys Chem A. 2004, 108 (46): 10359-10366. 10.1021/jp046928l.
Yang ZZ, Wang CS: Atom-bond electronegativity equalization method. 1. Calculation of the charge distribution in large molecules. J Phys Chem A. 1997, 101: 6315-6321. 10.1021/jp9711048.
Menegon G, Loos M, Chaimovich H: Parameterization of the electronegativity equalization method based on the charge model 1. J Phys Chem A. 2002, 106: 9078-9084. 10.1021/jp026083i.
Svobodová Vařeková R, Koča J: Optimized and parallelized implementation of the electronegativity equalization method and the atom-bond electronegativity equalization method. J Comput Chem. 2006, 3: 396-405.
Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, Montgomery JAJr, Vreven T, Kudin KN, Burant JC, Millam JM, Iyengar SS, Tomasi J, Barone V, Mennucci B, Cossi M, Scalmani G, Rega N, Petersson GA, Nakatsuji H, Hada M, Ehara M, Toyota K, Fukuda R, Hasegawa J, Ishida M, Nakajima T, Honda Y, Kitao O, Nakai H: Gaussian 09, Revision E.01. 2004, Wallingford: Gaussian, Inc.
Keith TA: AIMAll, Version 11.12.19. 2011, USA: TK Gristmill Software, Overland Park KS, [http://aim.tkgristmill.com]
Habibi-Yangjeh A, Danandeh-Jenagharad M, Nooshyar M: Application of artificial neural networks for predicting the aqueous acidity of various phenols using QSAR. J Mol Model. 2006, 12: 338-347. 10.1007/s00894-005-0050-6.
Hanai T, Koizumi K, Kinoshita T, Arora R, Ahmed F: Prediction of pKa values of phenolic and nitrogen-containing compounds by computational chemical analysis compared to those measured by liquid chromatography. J Chromatogr A. 1997, 762: 55-61. 10.1016/S0021-9673(96)01009-6.
Tehan BG, Lloyd EJ, Wong MG, Pitt WR, Montana JG, Manallack DT, Gancia E: Estimation of pKa Using semiempirical molecular orbital methods. Part 1: Application to phenols and carboxylic acids. Quant Struct-Act Relat. 2002, 21: 457-472. 10.1002/1521-3838(200211)21:5<457::AID-QSAR457>3.0.CO;2-5.
NCI Open Database Compounds. Retrieved from [http://cactus.nci.nih.gov/] on August 10, 2010
Sadowski J, Gasteiger J: From atoms and bonds to three-dimensional atomic coordinates: Automatic model builders. Chem ReV. 1993, 93: 2567-2581. 10.1021/cr00023a012.
Howard P, Meylan W: Physical/Chemical Property Database (PHYSPROP). 1999, North Syracuse NY: Syracuse Research Corporation, Environmental Science Center
Skřehota O, Svobodová Vařeková R, Geidl S, Kudera M, Sehnal D, Ionescu CM, Koča J: QSPR designer – a program to design and evaluate QSPR models. Case study on pKa prediction. J Cheminf. 2011, 3 (Suppl 1): P16-10.1186/1758-2946-3-S1-P16.
Bultinck P, Langenaeker W, Lahorte P, De Proft, Geerlings P, Van Alsenoy, Tollenaere JP: The electronegativity equalization method II: Applicability of different atomic charge schemes. J Phys Chem A. 2002, 106 (34): 7895-7901. 10.1021/jp020547v.
Lemm S, Blankertz B, Dickhaus T, Müller KR: Introduction to machine learning for brain imaging. NeuroImage. 2011, 56 (2): 387-399. 10.1016/j.neuroimage.2010.11.004.
Organisation for Economic Co-operation and Development: Guidance Document on the Validation of (Quantitative)Structure-Activity Relationships [(Q)SAR] Models. 2007, Paris: OECD, [http://search.oecd.org/officialdocuments/displaydocumentpdf/?cote=env/jm/mono(2007)2&doclanguage=en] (accessed April 6,2013)
This work was supported by the Ministry of Education of the Czech Republic (LH13055), by the European Community’s Seventh Framework Programme (CZ.1.05/1.1.00/02.0068 to J.K. and R.S.V.) from the European Regional Development Fund and by the EU Seventh Framework Programme under the “Capacities” specific programme (Contract No. 286154 to J.K.). C.M.I. and D.S. would like to thank Brno City Municipality for the financial support provided to them through the program Brno Ph.D. Talent. This work was also supported in part by NIH grants R01 GM071872, U01 GM094612, and U54 GM094618 to R.A.. The access to computing and storage facilities owned by parties and projects contributing to the National Grid Infrastructure MetaCentrum, provided under the programme “Projects of Large Infrastructure for Research, Development, and Innovations” (LM2010005) is highly appreciated. Also the access to the CERIT-SC computing facilities provided under the programme Center CERIT Scientific Cloud, part of the Operational Program Research and Development for Innovations, reg. no. CZ. 1.05/3.2.00/08.0144 is acknowledged.
The authors declare that they have no competing interests.
Electronic supplementary material
About this article
Cite this article
Vařeková, R.S., Geidl, S., Ionescu, C. et al. Predicting pK a values from EEM atomic charges. J Cheminform 5, 18 (2013). https://doi.org/10.1186/1758-2946-5-18
- Dissociation constant
- Quantitative structure-property relationship
- Partial atomic charges
- Electronegativity equalization method
- Quantum mechanics