Prediction of pKa values using the PM6 semiempirical method

Comparison of pKa values predicted using PM6 and ab initio methods

Sastre et al. (2012) have computed the pKa values using isodesmic reactions and a several ab initio method for a variety of molecules containing six types of ionizable groups. Table 1 lists the molecules from Sastre et al. (2012) used in this study. The molecules in the first row are the reference molecules (ref) with the corresponding pKa^ref value in parenthesis. Molecules containing chlorine have been eliminated because PM6 calculations for this elements involves d-integrals, which have not yet been implemented in GAMESS.

Columns 2–4 of Table 2 lists mean absolute differences (MADs) and maximum absolute differences (Max AD) relative to experiment for pKa values calculated by Sastre et al. (2012) using B3LYP and M05-2X/6-311++G(d,p) as well as the CBS-4B3* composite method (Casasnovas et al., 2010) and the SMD solvation method. The data shows that all three ab initio methods perform roughly equally well, with all three methods giving a MAD below 1 pH unit, with the exception of alcohols where the MAD ranges from 1.0–1.3 pH units. The Max ADs are lowest for amines (0.6–0.8 pH units) and highest for alcohols (2.3–2.9 pH units).

The fifth column lists the corresponding values computed using PM6-D3H+ with the SMD solvation method. For pyridines, alcohols, phenols, and benzoic acids the overall accuracy of PM6-D3H+ is comparable to the ab initio methods: the MADs are within 0.5 pH units of the ab initio values and while the Max ADs range from 0.4 (pyridines) to 2.4 (phenols). For carboxylic acids the results are dominated by a 3.5 pH unit error for trimethylacetic acid, without which the MAD is 1.0 pH units. Thus, different reference molecules should be used to predict pKa values for carboxylic acid groups bonded to secondary and tertiary carbons, using PM6 based methods. For amines the MAD and Max AD is 1.2 and 3.9 pH units, respectively. If only primary amines, which are most similar to the reference compound, are considered the MAD and Max AD drops to 0.5 and 1.2 pH units, respectively. We investigate this point further in the next subsection.

The sixth column of Table 2 lists PM6-D3H+/SMD* pKa values computed with the G_RRHO^°(X) term in Eq. (6) removed (denoted by the “*”). In all cases the change in MAD and Max AD is ≤ 0.2 and 0.3 pH units, respectively. This small change is not surprising the use of isodesmic reactions and approach has been used in pKa prediction before (Li, Robertson & Jensen, 2004). Neglecting the dispersion correction (PM6/SMD*) has an even smaller effect on the pKa values, changing the MAD and Max AD by at most 0.1 pH units. It is important to note that the molecules used in this part of the study are relatively small and contain only one functional group. The effect of neglecting vibrational free energies and dispersion corrections may have a bigger effect on the pKa values computed for larger molecules with, for example, intramolecular interactions where both dispersion and vibrational effects can play an important role.

The final column of Table 2 lists PM6/COSMO* pKa values. The pKa values for alcohols, phenols, and benzoic acids are very similar to PM6/SMD with MAD and Max ADs changing by at most 0.1 pH units. In the case of pyridines and carboxylic acids Max AD changes by 0.5 and −1.0 pH units, respectively although this only changes the MAD by at most 0.2 pH units. In the case of pyridines the PM6/SMD* and PM6/COSMO* Max AD is observed for 2,3-dimethylpyridine and 2,4-dimethylpyridine, respectively, while in the case of carboxylic acids the Max AD is observed for trimethylacetic acid. In the case of amines, the accuracy of PM6/SMD* and PM6/COSMO* is very similar for primary amines, but the error for di- and trimethylamine is reduced by 1.9 and 2.2 pH units, respectively, by using the COSMO solvation method implemented in MOPAC. To understand these differences, we look more closely at dimethylamine and compare the results to corresponding M05-2X/6-311++G(d,p)/SMD calculations, which is one of the methods used by Sastre et al. (2012), but used here without the G_RRHO^°(X) contribution to make the results directly comparable to PM6/SMD* and PM6/COSMO*. Both M05-2X/6-311++G(d,p)/SMD* and PM6/COSMO* yield pKa values for dimethylamine that are virtually identical in accuracy: 10.1 and 11.2 compared to the experimental value of 10.6 pH units. In the case of M05-2X/6-311++G(d,p)/SMD ΔE_ele (which replaces ΔH_f in Eq. (6)) and ΔΔG_solv^° the values are 11.4 and −10.7 kcal/mol, while the corresponding values for PM6/COSMO* are 3.5 and −4.2 kcal/mol. Taking M05-2X/6-311++G(d,p)/SMD* as a reference, the good performance of PM6/COSMO* is thus a result of significant error cancellation. The corresponding ΔΔG_solv^° computed using PM6/SMD* is −6.8 kcal/mol. While this value is closer to the M05-2X/6-311++G(d,p)/SMD* value it leads to worse error cancellation with the electronic energy contribution and therefore a less accurate pKa prediction (8.2 pH units).

In summary, the pKa values of the pyridines, alcohols, phenols, and benzoic acids considered in this study can generally be predicted with PM6 and ab initio methods to within the same overall accuracy, with average MADs for these four functional groups are 0.7–0.8 and 0.6–0.7 pH units, for the ab initio and PM6-based predictions. Similarly, the corresponding Max ADs ranges are 1.6–1.7 and 1.3–1.5 pH units, respectively. For carboxylic acids the PM6-based results are dominated by 2.3–3.5 pH unit errors for trimethylacetic acid, without which the MAD is 0.7–1.0 pH units and comparable to the corresponding ab initio results (0.6–0.7 pH units). Similarly, for amines the PM6-based results are dominated by a 1.9–4.1 pH unit errors for di- and trimethylamine, without which the MAD is 0.5–0.6 pH units and comparable to the corresponding ab initio results (0.2–0.3 pH units). For these simple molecules, dispersion corrections and vibrational free energy make a negligible contribution to the predicted pKa values.

Following Seybold & Shields (2015), Table 3 summarizes the overall statistics for the predictions presented in Table 2 (labeled “Sastre”) where outliers have been removed using the Modified Thompson τ method. As expected from our discussion above, the ab initio predictions are slightly better than the semiempirical results with r² and standard errors of 0.963–0.970 and 0.7–0.8 pH units compared to 0.933–0.951 and 1.0–1.2 pH units. PM6/COSMO is seen to perform slightly better than the other semiempirical approaches. The outliers are identified in Figs. 1A and 1B for the ab initio and semiempirical predictions. Trimethylamine and trimethyl acetic acid are outliers for all three SMD-based semiempirical predictions and the PM6/SMD* energy contributions are compared to the corresponding M05-2X/SMD* values to gain further insight. In both cases, the differences between PM6 and DFT is largest for the change in the gas phase deprotonation energy: 5.6 vs 1.5 and 3.7 vs 0.5 kcal/mol for trimethylamine and trimethylacetic acid, respectively.

Secondary and tertiary amines

Here we investigate whether the accuracy of PM6-based predictions of amines can be improved by using different reference molecules for primary, secondary, and tertiary amines. Table 4 lists experimental and predicted pKa values for six secondary and tertiary amines shown in Fig. 2 using di- and tri-ethylamine as respective reference. The accuracy of the predicted pKa values for secondary amines is slightly worse compared to primary amines (Table 2): the MADs and Max ADs are 0.5–1.0 and 1.0–1.6 pH units, respectively, compared to 0.5–0.6 and 1.2–1.4 pH units. The lowest MAD and Max AD is observed for PM6/COSMO*. The contributions of vibrational and dispersion effects are larger than for primary amines, with respective changes of upto 0.8 and 0.9 pH units—both observed for diallylamine. This is presumably due to the fact that the secondary amines are structurally more different from the reference (diethylamine) than for the primary amines. For example, if piperidine is taken as a reference for the prediction of the pKa of morpholine and piperazine then the effects of vibrations and dispersion contribute at most 0.1 pH units. For the SMD-based predictions, the lowest MAD is observed for PM6-D3H+ without vibrational contributions.

The accuracy of the predicted pKa values for tertiary amines is significantly worse than for primary and secondary amines with MADs and Max ADs of 1.0–2.8 and 2.1–4.4 pH units, respectively. As observed for secondary amines, the lowest and next-lowest MAD is observed for PM6/COSMO and PM6-D3H+/SMD*. For these two methods, the largest error is observed for DABCO: 3.2 and 2.1 pH units for PM6-D3H plus;/SMD* and PM6/COSMO, respectively. With the exception of diisopropylmethylamine, both methods underestimate the pKa values, and using the 2 pH unit correction proposed by Eckert & Klamt (2005) reduces the MAD and Max AD to 0.7 and 1.2 for PM6-D3H+/SMD* for these molecules, although the Max AD increases to 3.8 pH units if diisopropylmethylamine is included. Alternatively, the accuracy can be improved by changing the reference molecule. For example, using quinuclidine as a reference, the pKa of DABCO is predicted to within 0.9 and 0.5 pH units using PM6-D3H+/SMD* and PM6/COSMO, respectively.

In summary, the large errors observed for secondary and tertiary amines in Table 2 (i.e. di- and tri-ethylamine) can be decreased by using di- and tri-ethylamine as a reference. The MAD and Max AD for secondary amines (0.5–1.0 and 1.0–1.6 pH units) are only a little larger than those observed for primary amines (0.5–0.6 and 1.2–1.4). The MAD and Max AD for tertiary amines (1.0–2.5 and 2.1–4.5 pH units) are significantly larger than those observed for primary amines and secondary amines. As observed by Eckert & Klamt (2005) the pKa values tend to be underestimated and the error can be reduced somewhat by adding a 2 pH unit correction factor. Alternatively, the error can be reduced for individual molecules by choosing reference molecules with similar structures. PM6/COSMO results in the lowest errors, followed by PM6-D3H+/SMD* for both secondary and tertiary amines.

Table 3 summarizes the overall statistics for the primary amines in Table 2 and the amines in Table 4 (labeled “Amines”) where outliers have been removed using the Modified Thompson τ method. As expected from our discussion above, the predictions for amines are significantly worse than for the molecules in Table 1. In particular, the slopes deviate significantly from 1.0 and the intercept is in the range 4.0–6.2 pH units, due in part to the underestimated pKa values of tertiary amines. Interestingly, the standard error is 0.5–0.6 pH units, suggesting that reasonably accurate pKa predictions may be possible with the chosen reference molecules if the slope and intercept determined here are transferable to other systems. The outliers are identified in Fig. 1C. DABCO is an outlier for all four semiempirical predictions and the PM6/SMD* energy contributions are compared to the corresponding M05-2X/SMD* values to gain further insight. Again, the difference between PM6 and DFT is largest for the change in the gas phase deprotonation energy: 5.9 vs 4.2 kcal/mol.

Application to a drug-like molecule

We explore the effect of using different reference molecules further for compound 1 shown in Fig. 3. Settimo, Bellman & Knegtel (2013) have shown that the Chemaxon pKa predictor predicts a pKa value of 9.1 for compound 1, which is significantly higher than the experimental value of 4.2, i.e. Chemaxon predicts that 1 is charged as physiological pH when, in fact, it is neutral. Table 5 list the pKa values for 1 predicted using PM6-based methodologies using three different reference molecules (cf. Table 3). The absolute errors range from 1.7–8.5 with the error being smallest for PM6/SMD using tri-ethylamine as a reference. Given the size of compound 1 we expect that dispersion effects will make important contributions to intramolecular interactions and the difference in pKa values predicted with PM6-D3H+/SMD* and PM6/SMD* is indeed substantial (9.2–10.5 kcal/mol). The low error observed for PM6/SMD* is therefore likely fortuitous and, indeed, the error increases for reference molecules more closely related to 1, while the opposite is seen for PM6-D3H+/SMD(*). Furthermore, the PM6-D3H+/SMD(*) results are consistent with the near systematic pKa-underestimation observed for the tertiary amines in Table 4 and if the 2 pH unit correction suggested by Eckert & Klamt (2005) is used the error decreases to 3.7–4.1 pH units when benzylpyrrolidine or heliotridane are used as references. While these errors are substantial they lead to the correct qualitative prediction that 1 is neutral at physiological pH. However, whether PM6-based pKa predictions are sufficiently accurate to be useful in drug-design will require a great deal of additional study (see Summary and Outlook section for further information).

The computational cost of computing the free energy of a single conformation of 1 is ca 5 min on a single Intel Xeon 2.67 GHz X5550 core processor with the time roughly equally split between geometry optimization and vibrational frequency calculations. Thus, if the vibrational contributions to the standard free energy can be neglected the time requirement drops to 2–3 min per conformer per core processor. For 1 we computed the free energy of roughly 200 conformers. Thus, PM6-based pKa prediction is computationally efficient enough to be used for high throughput screening using on the order of 100 core processors.