Theoretical evaluation of the pka values of 5-Substitued uracil derivatives - Pham Le Nhan

Vietnam Journal of Science and Technology 55 (6A) (2017) 63-71 THEORETICAL EVALUATION OF THE pKa VALUES OF 5-SUBSTITUED URACIL DERIVATIVES Pham Le Nhan 1, * , Nguyen Tien Trung 2 1 Faculty of Chemistry, University of Dalat, 01 PhuDong Thien Vuong,Ward 8, Dalat, Viet Nam 2 Chemistry Department and Laboratory of Computational Chemistry and Modelling, Quy Nhon University, 170 An Duong Vuong, Quy Nhon, Viet Nam Email: nhanpl@dlu.edu.vn; nguyentientrung@qnu.edu.vn Revised: 15 July 2017; Accepted for publication 21 December 2017 ABSTRACT Density functional theory (DFT) calculations using numerical basis sets were employed to predict the solvation energies, Gibbs free energies and pKa values of a series of 5-substituted uracil derivatives. Obtained results show that solvation energies are not significantly different between DFT methods using the numerical (DNP) and Gaussian basis set (aug-cc-pVTZ). It is noteworthy that the independent and suitable solvation energy of proton of -258.6 kcal/mol has been proposed for the evaluation of pKa values in conjunction with the numerical basis set. In addition, the calculated pKa values suggest that the anti-conformation of 5-formyluracil is the most stable form in the aqueous solution. Keywords: numerical basis sets, Gaussian basis sets, pKa values, solvation energies, uracil. 1. INTRODUCTION There are different types of one-electron basis sets in computational chemistry including the Slater, Gaussian, numerical basis sets and others. Unlike Gaussian basis sets, numerical basis sets include basis functions which are generated numerically from the nucleus to an outer distance of 10 a.u. of each atom. Actually these functions have two parts involving spherical harmonic functions Ylm(θ,φ) as angular portions and radial functions F(r). Values of radial portions are obtained by numerically solving the atomic DFT equations [1]. Numerical basis sets have plenty of advantages in comparison to Gaussian ones and other analytical functions [1]. Although they tend to overestimate the structural information, especially bond lengths of some small molecules containing sulfur [2], numerical basis sets were proved to produce accurate geometric parameters of some hydrogen bonded systems, and it was clearly beneficial for optimization of large hydrogen bonded systems [3]. Another advantage is that numerical basis sets induce smaller basis set superposition error (BSSE) than that of Gaussian basis sets having the same size [4]. In addition, the key which is of considerable importance is the computational cost. It was demonstrated that numerical basis sets are much more efficient when compared to Gaussian basis sets [4–6]. Pham Le Nhan, Nguyen Tien Trung 64 The pKa parameters are of importance in many areas including chemistry and even in biology and pharmacy. Knowledge of pKa values could provide preliminary insights into chemical reactions and reaction mechanism, especially for the donor-acceptor proton reactions. Thus, pKa values of all types of organic compounds have been determined both experimentally and theoretically. The earlier work of Jorgensen and coworkers dates back to 1987 [7]. In this work, the authors have introduced the protocol for calculation of organic compounds pKa; all calculations used the 6-31+G(d) basis set which is one of the popular Gaussian basis sets at that time. Four years later, Lim et al. [8] presented a clearer procedure for calculation of pKa values applied to amino acids. This work was, however, conducted at the HF/6-31G(d) level of theory. To improve the theoretically calculated pKa results, some very high accurate methods for evaluating energies, for instance the use of complete basis set methods and Gaussian-n models have been applied in combination with continuum solvation methods to compute solvation energies and predict pKa values of several carboxylic acids [9]. The predicted pKa values were in excellent correlation with the experimental values. Some authors have evaluated a variety of theoretical methods using a sizable number of Gaussian basis sets [10] and found that the most accurate values of pKa were obtained at the B3LYP/6-31G(d) level. Some other researchers focused on determining pKa values of molecules which are of importance with regard to their roles in DNA [11, 12], and these studies used Gaussian basis set calculation methods as well. Until now, there are plentiful publications which have been done with the purpose of theoretically determining solvation energies and pKa values of numerous organic compounds, but these works have only been carried out using the HF, post HF and DFT methods using Gaussian basis sets. There is almost no research paying attention to numerical basis sets. For the aim of a fundamental evaluation, the present work has been conducted to calculate solvation energies and pKa values of some simple compounds with regard to the structures which are the derivatives of 5-substituted uracil. These structures are small, but they have several functional groups and additional elements beside carbon and hydrogen. The structures of these derivatives are given in Fig. 1. A couple of important facts, additionally, are that their experimental pKa values are available [13], as well as they have been studied with the DFT method using Gaussian basis sets [14] B3LYP/aug-cc-pVTZ (GBS), which is very convenient for the purpose of comparison between Gaussian and numerical basis sets. 2. MATERIALS AND METHODS 2.1. Protocol of pKa calculation It is evident that 5-substituted Uracils have two sites where the dissociation producing H + can happen [11]. These sites are at the N−H bonds marked as N1 and N3 in the structures of these derivatives depicted in Fig. 2. For the convenience, the chemical formulas of neutral compounds are written in an abbreviation way that involves the substituent only, while the chemical formulas of ions emerging from the dissociation additionally include the position number of nitrogen with the charge of ions. For example, cis-5-CHOC4H3N2 is going to be written as cis- CHO, and cis-5-CHOC4H2N2 - is denoted as cis-CHON3(-). All the shorted names of these derivatives can be seen in Table 1 and in other Tables of this text. In principle, pKa could be theoretically determined by quantum mechanical computations. Let’s consider a compound HA with ability of dissociation producing proton H+ and anion A–. The thermodynamic cycle of the dissociation in the gas phase and in the aqueous phase is illustrated in Fig. 3. To calculate pKa values of HA in the aqueous phase, the Gibbs free energy Theoretical evaluation of the pKa values of 5-substitued uracil derivatives 65 of HA dissociation in the aqueous phase ∆Gaq must be estimated. From the thermodynamic cycle in Fig. 3, ∆Gaq can be written as follows [8, 15]: ∆Gaq = ∆G(H + (aq)) +∆G(A-(aq)) − ∆G(HA(aq)) (1) where ∆G(H+(aq)), ∆G(A-(aq)) and ∆G(HA(aq)) are the Gibbs free energies of proton, anion A- and HA in solution respectively. With an arbitrary compound or ion, when transferred from gas phase into aqueous phase, the solvation happens and this process forms energy of solvation. For HA, the Gibbs free energy in the aqueous phase ∆G(HA(aq)) is given by ∆G(HA(aq)) = ∆G(HA(gas)) +∆G(HA(sol)) (2) where ∆G(HA(sol) and ∆G(HA(gas)) are the solvation energies and the Gibbs free energy in the gas phase of HA, respectively. For other terms in Eq. 1, other equations can be written as follows ∆G(A-(aq)) = ∆G(A-(gas)) +∆G(A-(sol)) (3) ∆G(H+(aq)) = ∆G(H+(gas)) +∆G(H+(sol)) (4) If the solvation energies and the Gibbs free energies in gas phase are available, ∆Gaq can be determined via Eq. 1. Thereafter, pKa of a compound is given by pKa = (1/2.303RT) ∆Gaq (5) Figure 1. Structures of Uracil’s derivatives. (a) uracil, (b) 5-fluorouracil, (c) thymine, (d) trans-5- formyluracil or syn conformation, (e) cis-5-formyluracil or anti conformation, and (f) 5-nitrouracil. Figure 2. Possible dissociation sites of HA. Pham Le Nhan, Nguyen Tien Trung 66 Figure 3. Dissociation of HA in the gas phase and aqueous phase. 2.2 Quantum mechanical calculations To calculate pKa values, we need to know the solvation energies and the Gibbs free energies of all compounds and their anions in the gas phase. These energies can be theoretically calculated. In this work, all calculations were conducted with Dmol 3 [1, 16]. Three DFT methods which are generalized gradient approximation (GGA) BLYP, B3LYP and local-density approximation (LDA) PWC have been used to predict the Gibbs free energies and the solvation energies. All the computational methods used the numerical basis set DNP (Double Numerical plus Polarization) because the DNP basis set is a typical basis set which can give accurate values of energies without computationally expensive demand [16]. For the solvation energies, the COSMO model was used to simulate the solvent (water) with dielectric constant of 78.54 at 298.15 K. In addition, in order to obtain the standard Gibbs free energy of each form in gas phase, the zero-point energy (ZPE) and the Gibbs free energies’ change from 0 to 298.15 K were also calculated by calculation of vibrational frequencies. The formula for calculating the standard Gibbs free energy presented by ∆G = E0K + ZPE +∆∆G0→298.15K (6) where the total energy of each form E0K was withdrawn from optimization of geometries. As for proton, the Gibbs free energy of proton in the gas phase ∆G(H+(gas)) was taken from literature [8, 15] given in Eq. 7. ∆G(H+(gas)) = 2.5RT − T∆S = 6.28 kcal/mol (7) Most of the solvation energies of the forms presented in the thermodynamic cycle above were extracted from quantum mechanical calculations except for that of proton. This energy was adopted from experiments. It is notable that the experimental solvation energy of proton is unclear. The value of this energy ranges from -252.6 to -262.5 kcal/mol [17–19] and this value even reaches -263.98 kcal/mol withdrawn from cluster-ion solvation data [20]. As mentioned above, this work is, basically, to analyze the capability of DFT methods using numerical basis sets in relation to those using Gaussian basis sets; therefore, the value of solvation energy of proton ∆G(H+(solv)) was chosen to be 258.32 kcal/mol which is the value used elsewhere in literature [14]. 3. RESULTS AND DISCUSSION 3.1. Solvation energies Theoretical evaluation of the pKa values of 5-substitued uracil derivatives 67 The energies of solvation are of quite importance for the calculation of pKa values. These values calculated at GGA/DNP, B3LYP/DNP and LDA/DNP levels of theory are tabulated in Table 1. It is obvious that solvation energies calculated by GGA/DNP, B3LYP/DNP and LDA/DNP are similar. These values differ from each other up to 3.5 kcal/mol, in which B3LYP/DNP gave the highest absolute values of solvation energies. These values in comparison to those calculated by GBS are systematic underestimates except for neutral trans-CHO and NO2 derivatives. Because of the highest absolute estimated values of solvation energies, B3LYP/DNP has the closest prediction to GBS, however in some cases, the differences approach 10 kcal.mol −1 for instance the predicted value of HN1(-). However, B3LYP/DNP gave -83.89 kcal/mol for the case of trans-CHON3(-), which is somehow quite different from the value calculated with GBS (80.72 kcal/mol), while two other methods gave values of around 80.5 kcal/mol. Table 1. Solvation energies (kcal/mol) of 5-substituted Uracils calculated by various DFT methods. State GGA/DNP, kcal/mol B3LYP/DNP, kcal/mol LDA/DNP, kcal/mol GBS [14], kcal/mol H -18.82 -19.85 -19.17 -20.29 HN1(-) -63.66 -65.03 -64.03 -74.31 HN3(-) -79.04 -81.30 -80.25 -87.94 F -19.83 -20.76 -19.83 -20.51 FN1(-) -60.90 -62.26 -61.65 -71.25 FN3(-) -75.52 -77.78 -76.71 -84.12 CH3 -18.03 -18.89 -18.37 -19.32 CH3N1(-) -64.26 -65.64 -64.88 -73.89 CH3N3(-) -78.71 -80.83 -79.65 -86.51 cis-CHO -20.29 -21.29 -21.21 -27.38 cis-CHON1(-) -57.20 -58.48 -57.47 -71.46 cis-CHON3(-) -70.76 -73.05 -71.14 -87.57 trans-CHO -28.22 -30.13 -28.26 -22.74 trans-CHON1(-) -62.30 -64.21 -62.37 -67.58 trans-CHON3(-) -80.40 -83.89 -80.94 -80.72 NO2 -26.21 -27.53 -26.21 -23.24 NO2N1(-) -57.60 -58.77 -57.74 -62.88 NO2N3(-) -72.82 -75.50 -73.67 -77.40 3.2. pKa values In Section 3.1, solvation energies of 18 forms emerging from the dissociation of 5 -substituted derivatives have been predicted. From these solvation energies, it is simple to estimate pKa values of these five derivatives via Eq. 5. Table 2 lists the pKa values of these compounds. Data from Table 2 reveals that both GGA/DNP and LDA/DNP gave very high overestimates in comparison with GBS and experimental values (Table 3) even double in some cases for example trans-CHON1(-) and NO2N1(-). The B3LYP/DNP method, conversely, Pham Le Nhan, Nguyen Tien Trung 68 estimated the closer values of pKa to the GBS and the experimental values. The maximum discrepancy between B3LYP/DNP and GBS is about 2. This consequence can be explained partly by the inconsistent calculation procedure between GBS and our procedure. To be more detailed, GBS calculated the solvation energies at 300 K [14] while our procedure was set up at 298.15 K, the same temperature of the experimental condition [13]. Table 2. The pKa values predicted by DFT methods using numerical basis sets. formation states of ions GGA/DNP B3LYP/DNP LDA/DNP GBS [14] HN1(-) 18.20 12.03 13.92 10.47 HN3(-) 16.48 9.67 11.88 9.34 FN1(-) 16.22 10.02 12.13 9.05 FN3(-) 14.01 7.10 9.49 7.26 CH3N1(-) 17.56 12.23 13.66 11.23 CH3N3(-) 16.25 9.96 11.61 10.04 cis-CHON1(-) 13.54 7.47 9.47 6.95 cis-CHON3(-) 14.06 7.84 9.84 7.28 trans-CHON1(-) 14.43 8.94 10.36 6.93 trans-CHON3(-) 13.11 7.32 8.93 7.96 NO2N1(-) 11.58 5.57 7.76 5.66 NO2N3(-) 12.21 5.05 7.99 6.91 3.3. New prediction of pKa and dominant conformation of 5-formyluracil in solution The B3LYP/DNP method has been proved to be a suitable alternative option in prediction of pKa values which are comparable to the results from the B3LYP/aug-cc-pVTZ. In this section, the evaluation of capability to predict pKa values independently of the proton solvation energy used by GBS methods [14] is conducted. For this purpose, a new value of the proton solvation energy must be re-determined. This value is taken from the experimental range in such a way that can accurately predict the values of pKa. Therefore, we proposed a new value of ∆G(H+(sol)) of -258.60 kcal/mol based on our preliminary estimations. Substitute this value to Eq. 4 to attain a new ∆G(H+(aq)) of -264.88 kcal/mol. By applying Eq. 5 a sets of new predicted pKa of five 5-substituted derivatives are presented in Table 3 together with the experimental values extracted from literature [13]. Only 5 experimental observations are available. It is interesting to note that most of the pKa values predicted by the B3LYP/DNP method are in good agreement with those withdrawn from the experiment except for the pKa of trans-5-formyluracil (8.73 versus 6.84). Because of large difference between two values, this value is believed not to be the pKa of the trans-5- formyluracil. Therefore, the experimental pKa of 5-formyluracil should belong to cis-5- formyluracil or in other words anti conformation is highly dominant in the real solution of 5- formyluracil. Interestingly, in 2004, Rogstad’s group also proved that the anti-conformation of 5-formyluracil is more favorable in the environment of DNA [21]. This again reinforces our idea that cis-5-formyluracil is highly populated in the aqueous solution. Theoretical evaluation of the pKa values of 5-substitued uracil derivatives 69 Table 3. The independent pKa values predicted by at B3LYP/DNP. |pKanum-pKaex.| and |pKaGBS-pKaex.|: the subtractions of experimental pKa (pKaex.) from the pKa values predicted by DFT methods with the numerical basis set (pKanum) and Gaussian basis sets pKaGBS (pKaGBS). State B3LYP/DNP Ex. [13] |pKanum-pKaex.| |pKaGBS-pKaex.| HN1(-) 11.83 × × × HN3(-) 9.46 9.42 0.04 0.08 FN1(-) 9.81 × × × FN3(-) 6.89 7.93 1.04 0.67 CH3N1(-) 12.03 × × × CH3N3(-) 9.76 9.75 0.01 0.29 cis-CHON1(-) 7.27 6.84 0.43 0.11 cis-CHON3(-) 7.63 × × × trans-CHON1(-) 8.73 6.84 1.89 * 0.09 * trans-CHON3(-) 7.12 × × × NO2N1(-) 5.36 5.3 0.06 0.36 NO2N3(-) 4.85 × × × 4. CONCLUSION A fundamental evaluation of density functional methods using the numerical basis set has been presented in this work. Our results show that the numerical basis set (DNP) is pretty suitable for the use of solvation energy calculations for 5-substituted Uracils especially in combination with B3LYP. Accordingly, B3LYP/DNP is good at prediction of pKa values for 5- substituted Uracils, which are comparable to those calculated by another high level of theory B3LYP/aug-cc-pVTZ. In addition, we have successfully introduced a novel value of -258.60 kcal/mol of proton solvation energy for a DFT method utilizing the numerical basis set to predict pKa values of 5-substituted Uracils. The predicted pKa values are in good agreement with the experimental results. And finally, we proposed that anti conformation of 5-formyluracil is dominant in the aqueous solution, which is in part consistent with the report from the literature. REFERENCES 1. Delley B. - An all-electron numerical method for solving the local density functional for polyatomic molecules, The Journal of chemical physics 92 (1) (1990) 508–517. 2. Altmann J. A., Handy N. C., and Ingamell V. E. - A study of the performance of numerical basis sets in DFT calculations on sulfur-containing molecules, International journal of quantum chemistry 57 (4) (1996) 533–542. Pham Le Nhan, Nguyen Tien Trung 70 3. Benedek N. A., Snook I. K., Latham K., and Yarovsky I. - Application of numerical basis sets to hydrogen bonded systems: a density functional theory study, The Journal of chemical physics 122 (14) (2005) 144102. 4. Inada Y. and Orita H. - Efficiency of numerical basis sets for predicting the binding energies of hydrogen bonded complexes: evidence of small basis set superposition error compared to Gaussian basis sets, Journal of computational chemistry 29 (2) (2008) 225–232. 5. Delley B. - Analytic energy derivatives in the numerical local-density-functional approach, The Journal of chemical physics 94 (11) (1991) 7245–7250. 6. Henry D. J., Varano A., and Yarovsky I. - Performance of numerical basis set DFT for aluminum clusters, The Journal of Physical Chemistry A 112 (40) (2008) 9835–9844. 7. Jorgensen W. L., Briggs J. M. and Gao J. - A priori calculations of pKa’s for organic compounds in water. The pKa of ethane, Journal of the American Chemical Society 109 (22) (1987) 6857–6858. 8. Lim C., Bashford D., and Karplus M. - Absolute pKa calculations with continuum dielectric methods, The Journal of Physical Chemistry 95(14) (1991) 5610– 5620. 9. Toth A. M., Liptak M. D., Phillips D. L., and Shields G. C. - Accurate relative pKa calculations for carboxylic acids using complete basis set and Gaussian-n models combined with continuum solvation methods, The Journal of Chemical Physics 114 (10) (2001) 4595–4606. 10. Liton M. A. K., Ali M. I. and, Hossain M. T. - Accurate pKa calculations for trimethylaminium ion with a variety of basis sets and methods combined with CPCM continuum solvation methods, Computational and Theoretical Chemistry 999 (2012) 1–6. 11. Jang Y. H., Goddard W. A., Noyes K. T., and Sowers L. C. - First Principles Calculations of the Tautomers and pKa Values of 8-Oxoguanine: Implications for Mutagenicity and Repair, Chemical research in toxicology 15(8) (2002) 1023–1035. 12. Rogstad K. N., Jang Y. H., Sowers L. C., and Goddard W. A. - First Principles Calculations of the p K a Values and Tautomers of Isoguanine and Xanthine, Chemical research in toxicology 16 (11) (2003) 1455–1462. 13. Privat E. J. and Sowers L. C. - A proposed mechanism for the mutagenicity of 5- formyluracil, Mutation Research/Fundamental and Molecular Mechanisms of Mutagenesis 354 (2) (1996) 151–156. 14. Jang Y. H., Sowers L. C., Çaǧ in T., and Goddard W. A. - First principles calculation of pKa values for 5-substituted uracils, The Journal of Physical Chemistry A 105(1) (2001) 274–280. 15. Topol I. A., Tawa G. J., Burt S. K., and Rashin A. A. - Calculation of absolute and relative acidities of substituted imidazoles in aqueous solvent, The Journal of Physical Chemistry A 101 (51) (1997) 10075–10081. 16. Delley B. - From molecules to solids with the DMol3 approach, The Journal of chemical physics 113 (18) (2000) 7756–7764. 17. Reiss H. and Heller A. - The absolute potential of the standard hydrogen electrode: a new estimate, The Journal of Physical Chemistry 89 (20) (1985) 4207–4213. Theoretical evaluation of the pKa values of 5-substitued uracil derivatives 71 18. Marcus Y. - Thermodynamics of solvation of ions. Part 5. Gibbs free energy of hydration at 298.15 K Journal of the Chemical Society, Faraday Transactions 87 (18) (1991) 2995–2999. 19. Gilson M. K. and Honig B. - Calculation of the total electrostatic energy of a macromolecular system: solvation energies, binding energies, and conformational analysis, Proteins: Structure, Function, and Bioinformatics 4 (1) (1988) 7–18. 20. Tissandier M. D., Cowen K. A., Feng W. Y., Gundlach E., Cohen M. H., Earhart A. D., Coe J. V., and Tuttle T. R. - The proton’s absolute aqueous enthalpy and Gibbs free energy of solvation from cluster-ion solvation data, The Journal of Physical Chemistry A 102 (40) (1998) 7787–7794. 21. Rogstad D. K., Heo J., Vaidehi N., Goddard W. A., Burdzy A., and Sowers, L. C. - 5- Formyluracil-induced perturbations of DNA function, Biochemistry 43 (19) (2004) 5688–5697.