Computational Estimation of the Acidities of Pyrimidines and Related Compounds.

Pyrimidines are key components in the genetic code of living organisms and the pyrimidine scaffold is also found in many bioactive and medicinal compounds. The acidities of these compounds, as represented by their pKas, are of special interest since they determine the species that will prevail under different pH conditions. Here, a quantum chemical quantitative structure-activity relationship (QSAR) approach was employed to estimate these acidities. Density-functional theory calculations at the B3LYP/6-31+G(d,p) level and the SM8 aqueous solvent model were employed, and the energy difference ∆EH2O between the parent compound and its dissociation product was used as a variation parameter. Excellent estimates for both the cation → neutral (pKa1, R2 = 0.965) and neutral → anion (pKa2, R2 = 0.962) dissociations were obtained. A commercial package from Advanced Chemical Design also yielded excellent results for these acidities.

1. Introduction

Pyrimidines play a central role in the terrestrial genetic code and the pyrimidine framework is also found in many bioactive compounds and medicines. As measurements of molecular properties are frequently difficult or expensive, it is of interest to develop theoretical means for estimating these properties. Among the most important and interesting properties of the pyrimidines are their acidities, as represented by their pKa values, which determine the forms of the compounds that will prevail in solution under different pH conditions. As a result, there has been a long-standing interest in estimating the pKas of chemical compounds using theoretical approaches [1,2]. This interest continues, as demonstrated by the broad range of methods employed in recent pKa studies [3,4,5,6,7,8,9,10]. In an earlier study by our group [11], we presented computational estimates of the pKas of the biologically related purines and indoles. In the present work we develop estimates for the acidities of pyrimidines and related compounds.

Three main approaches have been used in studies estimating compound acidities [1,12]. In the first approach, first-principles or absolute pKas are determined in a straight-forward manner using direct calculations, often relying on a thermodynamic cycle to separate different hypothetical stages [13,14]. The advantage of this approach is that it follows standard procedures and does not depend on prior knowledge of experimental results. Its disadvantage is that it normally requires a high level of computational effort to achieve reasonable accuracy. A second approach employs the development of an appropriate Quantitative Structure–Activity Relationship (QSAR) for the acidities of a set of compounds [2,15]. This approach requires first the collection of measured experimental pKa values for the set and then the identification of some suitable molecular parameter that is closely related to the acidity. An advantage of this approach is that it can generally achieve high accuracy while employing more modest computational efforts, and it accordingly allows the estimation of the pKas of related, unreported compounds using the same modest computational effort. Furthermore, digressions from the derived regression algorithms (i.e., outliers) can alert practitioners to compounds that may be exhibiting behaviors different from the others in the sample [16,17,18]. In a third approach, a number of commercial programs use algorithms derived from large acidity databases and employ empirical parameters such as Hammett constants to estimate pKas [19,20]. The latter two approaches were employed in the present study.

As noted, the QSAR approach relies on the discovery of some parameter or property of the compounds examined that correlates with the activity of interest, in the present case the pKa. In many cases, partial atomic charges have been employed as variation parameters for acidities and other properties. However, because the notion of an “atomic charge” in a molecule is not a proper quantum chemical observable, a variety of different approximate schemes have been developed to represent these atomic charges. Of these proposed schemes, our own group [21,22,23,24] and others [5,25,26,27] have found the natural population analysis (NPA) orbitals and charges developed by Weinhold et al. [28,29,30] to be especially useful in this role. More recently, we have also used another parameter, the energy difference ∆EH2O between the parent compound and its main dissociation product for this purpose. It is this parameter that was employed in the present study.

A particular difficulty arises in a study of the pyrimidines because these compounds do not typically appear as a single species in gas phase or solution, but rather are present as a collection of related tautomers, a situation that also prevailed in our earlier study of purines [11]. Accordingly, some accommodation must be made for this condition, as will be described in the following section.

2. Methods

Our aim in a recent series of reports has been to design and evaluate a QSAR protocol suitable for producing accurate pKa estimates for selected classes of compounds while employing relatively modest and commonly available computational tools.

The initial step in a pKa QSAR study involves the collection of reported experimental pKa values from the literature for the class of compounds of interest. In the present case, values for both the cation-to-neutral dissociation, AH2+ → AH + H+, which we will designate pKa1, and the neutral-to-anion dissociation, AH → A− + H+, which we will call pKa2, were available for a number of pyrimidines and related compounds. These values are tabulated in Table 1. Given in Table 1 also are computed pKa values obtained from the Advanced Chemical Development (ACD) [31] commercial software package.

Table 1

Reported experimental pKa1 and pKa2 values and ACD computed values for the pyrimidines and related heterocycles studied.

No.	Compound	Formula	pKa1	ACD pKa1	pKa2	ACD pKa2
1	azauracil	C3H3N3O2	-	−4.4 ± 0.2	-	7.8 ± 0.2
2	aziridine	C2H5N	7.98 [32], 8.05 [33]	8.1 ± 0.2	-	-
3	creatinine	C4H7N3O	4.8 [34]	4.8 ± 0.1	-	-
4	cytosine	C4H5N3O	4.32 [30], 4.58 [35], 4.6 [31]	4.4 ± 0.1	13 [33], 12.15 [36]	12.3 ± 0.1
5	flucytosine	C4H3FN2O2	3.26 [34]	2.6 ± 0.1	-	10.5 ± 0.1
6	imidazole	C3H4N2	7.15 [33], 6.99 [34,35]	7.2 ± 0.6	14.44 [33]	13.9 ± 0.1
7	1-methylimidazole	C4H6N2	6.95 [34]	7.0 ± 0.1	N/A	-
8	4-methylimidazole	C4H6N2	7.55 [35]	7.7 ± 0.6	-	14.3 ± 0.1
9	isocytosine	C4H5N3O	4.01 [34]	3.4 ± 0.5	9.42 [36]	9.6 ± 0.4
10	isoxazole	C3H6NO	−2.0 [34]	−2.0 ± 0.5	N/A
11	maleimide	C4H3NO2	-	−5.7 ± 0.2	9.46 [37]	8.5 ± 0.2
12	morpholine	C4H9NO	8.492 [35]	9.0 ± 0.2	N/A	-
13	piperidine	C5H11N	11.12 [35]	10.4 ± 0.1	-	-
14	piperazine	C4H10N2	9.78 [35], 9.73 [34]	9.6 ± 0.1	-	-
15	1-methylpiperazine	C5H12N2	10.19 [35]	9.3 ± 0.1	-	-
16	oxazole	C3H6NO	0.8 [34]	1.0 ± 0.1	N/A	-
17	pyrazine	C4H4N2	0.65 [34]	1.2 ± 0.1	N/A	-
18	pyrazole	C3H4N2	2.61 [35]	2.8 ± 0.1	14.21 [33]	14.0 ± 0.5
19	pyridazine	C4H43N2	2.3 [32]	3.1 ± 0.1	N/A	-
20	pyridine	C5H5N	5.23 [34]	5.2 ± 0.1	N/A	-
21	pyrimidine	C4H4N2	1.3 [32]	1.8 ± 0.1	N/A	-
22	pyrrole	C4H5N	−3.8 [34]	−0.3 ± 0.5	17.0 [33]	17.0 ± 0.5
23	pyrrolidine	C4H9N	12.10 [33], 11.31 [34,35]	10.5 ± 0.1	-	-
24	succinimide	C4H5NO2	-	−4.4 ± 0.2	9.62 [34,35], 9.68 [36]	9.6 ± 0.1
25	thymine	C5H6N2O2	-	−4.1 ± 0.4	9.9 [33], 9.79 [35], 9.44 [34]	9.2 ± 0.1
26	uracil	C4H4N2O2	-	−4.2 ± 0.1	9.43 [38], 9.45 [34,36]	8.9 ± 0.1
27	5-bromouracil	C4H3BrN2O2	-	-	7.91 [38]	6.8 ± 0.1
28	5-chlorouracil	C4H3CIN2O2	-	-	7.92 [38]	6.8 ± 0.1
29	fluorouracil	C4H4FN3O	-	-	8.04 [30], 8.00 [39], 7.93 [38]	6.7 ± 0.1
30	5-formyluracil	C5H4N2O3	-	-	6.84 [38]	7.3 ± 0.1
31	5-nitrouracil	C4H3N3O4	-	-	5.3 [38]	5.2 ± 0.1

Calculations were carried out using the Spartan’10 software package (Wavefunction, Inc., Irvine, CA). In earlier pKa studies, we found that density-functional computations at the B3LYP/6-31+G(d,p) level provided accurate accounts of molecular properties while still requiring only modest computational demands. After testing this assumption (vide infra) against available gas-phase experimental results for pyrimidines in the NIST database [40], this level of theory was adopted in the present work. For the studies in aqueous solution, the SM8 aqueous solvent model of Marenich et al. [41,42] was used. This same solvent model was also employed in our earlier study of purines and indoles [11] and other studies [16,17,43], where it was shown to perform well in helping to reproduce the experimental data.

As noted above, many of the compounds examined here exist in several tautomeric forms in solution. (For example, uracil has three cationic tautomers, six neutral tautomers, and two anionic tautomers.) Accordingly, the relative stabilities of the tautomeric forms—cationic, neutral, and anionic—of each compound were evaluated within the SM8 aqueous solvent model, and the most stable tautomer for each condition was taken as representative of that compound for computational purposes [44]. Fortunately, there frequently exists a substantial energy gap (>25 kJ/mol) between the most stable tautomer and the remaining tautomers of that species; however, ultimately the validity of this approximation must await justification from the results of the subsequent analysis. In earlier studies, we found that for neutral → anion dissociations (pKa2) the value ∆EH2O = EH2O(A−) − EH2O(AH) for the energy difference between the parent compound AH and its dissociation product A- in aqueous solution provides an excellent regression parameter for QSAR pKa studies, and after testing several other parameters this descriptor was used in the present studies. For the cation → neutral dissociation (pKa1), the analogous expression ∆EH2O = EH2O(AH) − EH2O(AH2+) was used.

We also provide a note of caution regarding directly comparing numerical results found here using Spartan’10 with those obtained using the popular Gaussian computational package (Gaussian, Inc., Wallingford, CT 06492, USA), since these programs use different basis sets for some atoms [18,45].

3. Results and Discussion

We first wished to assure that the B3LYP/6-31+G(d,p) level of computation was sufficient to provide accurate results for the dissociations in question. For this, we turned to gas-phase reaction results reported in the NIST chemical database [24]. This database contains gas-phase thermodynamic data for the Gibbs energy change ∆rG° for the anion + H+ → neutral reaction, for six of the compounds studied here. These experimental data are compared with our computed ∆rG° values in Table 2 (note that G°(H+) is −26.3 kJ/mol at 298.15 K [1]). As can be seen, there is very close agreement between the experimental and calculated ΔrG° values. With the exception of succinimide, all of the calculated values fall within the estimated experimental errors given in the NIST database. This, and the coefficient of determination of R2 = 0.998 between the computed and experimental values, suggest that the level of computation described should provide a good account of the reactions to be considered.

Table 2

Experimental and calculated gas-phase ΔrG° values (kJ/mol) for the reaction A− + H+ → AH.

Compound	Exp. ΔrG° a	Calc. ΔrG° b	Calc. ΔE b
pyridine	1601	1605	1648
pyrazine	1605	1605	1643
pyrimidine	1577	1579	1614
pyridazine	1565	1562	1601
imidazole	1433	1432	1466
succinimide	1414	1401	1436

a From the NIST database, Ref. [40]; b B3LYP/6-31+G(d,p).

As noted above, in previous studies we have found that the energy difference ΔEH2O between the parent compound and its dissociation product provides an excellent regression variable for pKa QSAR estimations. We first examined employment of the gas-phase ΔE values for this purpose. As expected, the use of ΔEgas values yielded good, but not exceptional, correlations for both the pKa1 (cation → neutral, R2 = 0.707) and pKa2 (neutral → anion, R2 = 0.874) dissociations.

We next optimized the compounds within the SM8 aqueous solvent model. Using the solvent-optimized structures and the calculated ∆EH2O values for the appropriate reactions, we obtained the following QSAR models:pK n = 17, R pK n = 12, R

In these equations, n = the number of compounds in the sample, R2 is the coefficient of determination (the fraction of the variance in the data accounted for by the model), s is the standard error of the estimate, and F is the Fisher statistic. It is evident that optimization of the structures within the solvent model significantly increases the accuracy of the model, as was also shown in earlier work [1,7,8,9].

The results for pKa1 and pKa2 are plotted in Figure 1 and Figure 2, and the calculated values are compared with the experimental values in Table 3 and Table 4. We note that several of the pKas for the cation → neutral dissociation (pKa1) fall into the difficult-to-measure negative value range and carry large uncertainties. Accordingly, these values are also not well characterized for use in this range in forming the regression Equation (1), and we prefer to recognize this uncertainty by simply indicating modestly negative (<0) or significantly negative (<<0) for the pKas of these compounds.

Figure 1

Plot of experimental pKa1s vs. calculated ΔEH2Os.

Figure 2

Plot of experimental pKa2s vs. calculated ΔEH2Os.

Table 3

Literature pKa1s and estimated pKa1s.

Compound	ΔE kJ/mol	Exp. pKa1	Calc. pKa1 a	Residual
azauracil	−1084	-	<<0	-
aziridine	−1216	8.01	7.76	0.25
creatinine	−1188	4.8	4.09	0.71
cytosine	−1183	4.5	3.43	1.07
flucytosine	−1169	3.26	1.60	1.66
imidazole	−1213	7.07	7.36	−0.29
1-methylimidazole	−1216	7.95	7.76	0.19
4-methylimidazole	−1215	7.55	7.63	−0.08
isocytosine	−1183	4.01	3.43	0.58
isoxazole	−1115	−2	<0	-
maleimide	−1035	-	<<0	-
oxazole	−872	0.8	<<0	-
piperidine	−1239	11.12	10.77	0.35
piperazine	−1226	9.76	9.07	0.69
1-methylpiperazine	−1231	10.19	9.72	0.47
pyrazine	−1166	0.65	1.21	−0.56
pyrazole	−1174	2.61	2.25	0.36
pyridazine	−1186	2.3	3.83	−1.53
pyridine	−1200	5.23	5.66	−0.43
pyrimidine	−1172	1.3	1.99	−0.69
pyrrolidine	−1239	11.71	10.77	0.94
succinimide	−1007	-	<<0	-

a pKa1s were estimated using Equation (1).

Table 4

Literature pKa2s and estimated pKa2s.

Compound	ΔE kJ/mol	Exp. pKa2	Calc. pKa2 a	Residuals
azauracil	−1179	-	6.82	-
aziridine	−1325	-	27.41	-
creatinine	−1217	-	12.18	-
cytosine	−1216	12.57	12.04	0.53
isocytosine	−1196	9.42	9.22	0.20
flucytosine	−1254	-	17.39	-
imidazole	−1216	14.4	12.04	2.36
4-methylimidazole	−1224	-	13.16	-
maleimide	−1196	9.5	9.22	0.28
piperazine	−1313	-	25.71	-
1-methylpiperazine	−1395	-	37.28	-
piperidine	−1310	-	25.29	-
pyrrole	−1243	17	15.84	1.16
pyrrolidine	−1390	-	36.57	-
thymine	−1199	9.71	9.64	0.07
uracil	−1200	9.44	9.78	−0.34
5-bromouracil	−1184	7.91	7.52	0.39
5-chlorouracil	−1182	7.92	7.24	0.68
fluorouracil	−1193	7.99	8.79	−0.80
5-formyluracil	−1180	6.84	6.96	−0.12
5-nitrouracil	−1159	5.3	4.00	1.30

a pKa2s were estimated using Equation (2).

We also tested the ability of a commercial software program, Advanced Chemical Development, Inc.’s ACD/Labs PhysChem Percepta Suite, to estimate these pKas. The results for pKa1 showed an excellent correlation:pk where n = 19, R

The results for pKa2 were also very good:pK where n = 13, R

These results encourage use of this software for studies of the pKas of these compounds.

4. Conclusions

A primary endeavor of scientific studies is to develop models of physical, chemical, and biological systems for the purpose of understanding these systems better. All models are by their very nature approximate. However, as Gauch has noted [46], in some cases—counterintuitively—a model can be more accurate than the data from which it is constructed “because it amplifies hidden patterns and discards unwanted noise” inherent in the system being examined. The QSAR equations used here take advantage of this property by “averaging through” the noise, or random errors, in the experimental pKa data. It is evident that both the QSAR Equations (1) and (2) above provide relatively simple means, via mathematical models, for estimating the pKas of the pyrimidines and related compounds. For example, in order to estimate the pKas of unmeasured compounds in this class or to check the reported pKas of measured compounds, using the QSAR equations one merely needs to determine ∆EH2O for the compound and then evaluate the pKa from the appropriate regression equation. Therefore, the equations provided should allow reasonable estimations for the pKas of other pyrimidines and compounds similar to the pyrimidines. Use of the commercial ACD/Labs program can provide a further independent and very useful check on the pKa estimates.