Theoretical evaluation of pK(a) in phosphoranes: implications for phosphate ester hydrolysis

Knowledge of the pK(a) of phosphoranes is important for the interpretation of phosphate ester hydrolysis. Calculated pK(a)'s of the model phosphorane, ethylene phosphorane, are reported. The method of calculation is based on the use of dimethyl phosphate as a reference state for evaluating relative pK(a) values, and on the optimization of the oxygen and acidic hydrogen van der Waals radii to give reasonable pK(1)(a), pK(2)(a), and pK(3)(a) for phosphoric acid in solution. Density functional theory is employed to calculate the gas-phase protonation energies, and continuum dielectric methods are used to determine the solvation corrections. The calculated pK(1)(a) and p(2)(a) for the model phosphorane are 7.9 and 14.3, respectively. These values are within the range of proposed experimental values, 6.5-11.0 for pK(1)(a), and 11.3-15.0 for pK(2)(a). The mechanistic implications of the calculated pK(a)'s are discussed.     

Improved pK(a) prediction: combining empirical and semimicroscopic methods

Using three different methods we tried to compute 171 experimentally known pK(a) values of ionizable residues from 15 different proteins and compared the accuracies of computed pK(a) values in terms of the root mean square deviation (RMSD) from experiment. One method is based on a continuum electrostatic model of the protein including conformational flexibility (KBPLUS). The others are empirical approaches with PROPKA deploying physically motivated energy terms with adjustable parameters and PKAcal using an empirical function with no physical basis. PROPKA reproduced the pK(a) values with highest overall accuracy. Differentiating the data set into weakly and strongly shifted experimental pK(a) values, however, we found that PROPKA's accuracy is better if the pK(a) values are weakly shifted but on equal footing with that of KBPLUS for more strongly shifted values. On the other hand, PKAcal reproduces strongly shifted pK(a) values badly but weakly shifted values with the same accuracy as PROPKA. We tested different consensus approaches combining data from all three methods to find a general procedure for most accurate pK(a) predictions. In most of the cases we found that the consensus approach reproduced experimental data with better accuracy than any of the individual methods alone.     

Determination of the intrinsic site pK(a) value and cooperativity of the symmetrical hexadentate chelator N,N',N'-tris[2-(3-hydroxy-2-oxo-1,2-dihydropyridin-1-yl)acetamido]ethylamine

The three overlapping pK(a) values of N,N',N'-tris[2-(3-hydroxy-2-oxo-1,2-dihydropyridin-1-yl)acetamido]ethylamine, a tripodal hexadentate chelator formed from three 3-hydroxy-2(1H)-pyridinone moieties amide linked to tris-(2-aminoethyl)amine, were determined by simultaneous spectrophotometric and potentiometric titration. The data was analysed by non-linear regression with constraints to deal with (a) the highly correlated absorptivities and (b) the highly correlated pK(a) values. The three pK(a) values were optimized first from the spectrophotometric data (absorbance vs. pH) by non-linear regression to a model in which the molar absorptivity of the ith species ((i)) was constrained by the correlation equation (i) = epsilon (0) + (epsilon (3) - epsilon (0))i 3 with i = 0, 1, 2, 3, where (3) and (0) represent the molar absorptivities of the most protonated and least protonated species, respectively. The molar absorbitivity of the four species defined by three pK(a) values is, therefore, linearly related to proton stoichiometry. The pK(a) values were then optimized from the potentiometric data (pH vs. titrant volume) by non-linear regression to a model in which the three pK(a) values were constrained by the correlation equation pK(a(i)) = pK(a(int)) + b(i - 1) + (i - 2)log(3) where i = 1, 2 or 3. This expresses the three pK(a) values in terms of only two optimizable parameters, the intrinsic site pK(a) (pK(a(int))) and the interaction energy between sites (b). The fixed term (i - 2)log(3) accounts for the statistical effect on the pK(a) values of three equivalent ionizable sites. The modified analytical derivatives required for optimization of these parameters by the Gauss-Newton-Marquardt algorithm and the merits of optimizing pK(a) values with these two correlation equations are discussed. The optimized pK(a) values were 9.31 +/- 0.01, 8.75 +/- 0.01 and 8.19 +/- 0.01. The separation between pK(a) values is 0.58 comprising 0.477 for the statistical effect and 0.081 for the interaction energy while the intrinsic site pK(a) is 8.672 +/- 0.005. The tertiary amine at the centre of the tripodal backbone has a pK(a) of 5.88 +/- 0.03.     

Effects of pH on protein association: modification of the proton-linkage model and experimental verification of the modified model in the case of cytochrome c and plastocyanin.

Effects of pH on protein association are not well understood. To understand them better, we combine kinetic experiments, calculations of electrostatic properties, and a new theoretical treatment of pH effects. The familiar proton-linkage model, when used to analyze the dependence of the association constant K on pH, reveals little about the individual proteins. We modified this model to allow determination not only of the numbers of the H+ ions involved in the association but also of the pK(a) values, in both the separate and the associated proteins, of the side chains that are responsible for the dependence of K on pH. Some of these side chains have very similar pK(a) values, and we treat them as a group having a composite pK(a) value. Use of these composite pK(a) values greatly reduces the number of parameters and allows meaningful interpretation of the experimental results. We experimentally determined the variation of K in the interval 5.4 < or = pH < or = 9.0 for four diprotein complexes, those that the wild-type cytochrome c forms with the wild-type plastocyanin and its mutants Asp42Asn, Glu59Gln, and Glu60Gln. The excellent fittings of the experimental results to the modified model verified this model and revealed some unexpected and important properties of these prototypical redox metalloproteins. Protein association causes a decrease in the pK(a) values of the acidic side chains and an increase in the pK(a) values of the basic side chains. Upon association, three carboxylic side chains in wild-type plastocyanin each release a H+ ion. These side chains in free plastocyanin have an anomalously high composite pK(a) value, approximately 6.3. Upon association, five or six side chains in cytochrome c, likely those of lysine, each take up a H+ ion. Some of these side chains have anomalously low pK(a) values, less than 7.0. The unusual pK(a) values of the residues in the recognition patches of plastocyanin and cytochrome c may be significant for the biological functions of these proteins. Although each mutation in plastocyanin markedly, and differently, changed the dependence of K on pH, the model consistently gave excellent fittings. They showed decreased numbers of H+ ions released or taken up upon protein association and altered composite pK(a) values of the relevant side chains. Comparisons of the fitted composite pK(a) values with the theoretically calculated pK(a) values for plastocyanin indicated that Glu59 and Asp61 in the wild-type plastocyanin each release a H+ ion upon association with cytochrome c. Information of this kind cannot readily be obtained by spectroscopic methods. Our modification of the proton-linkage model is a general one, applicable also to ligands other than H+ ion and to processes other than association.     

Prediction of 3-hydroxypyridin-4-one (HPO) hydroxyl pKa values

As an aid in optimising the design of 3-hydroxypyridin-4-ones (HPOs) intended for use as therapeutic Fe(3+) chelating agents, various quantum mechanical (QM) and semi-empirical (QSAR) methods have been explored for predicting the pK(a) values of the hydroxyl groups in these compounds. Using a training set of 15 HPOs with known hydroxyl pK(a) values, reliable predictions are shown to be obtained with QM calculations using the B3LYP/6-31+G(d)/CPCM model chemistry (with Pauling radii, and water as solvent). With this methodology, the observed hydroxyl pK(a) values for the training set compound are closely matched by the predicted pK(a) values, with the correlation between the observed and predicted values giving r(2) = 0.98. Predictions subsequently made by this method for a test set of 48 HPOs of known hydroxyl pK(a) values (11 of which were determined experimentally in this study), gave predicted pK(a) values accurate to within ±0.2 log units. In order to further investigate the predictive power of the method, two novel HPOs were synthesised and their hydroxyl pK(a) values were determined experimentally. Comparison of these predicted pK(a) values against the measured values gave absolute deviations of 0.13 (10.18 vs. 10.31) and 0.43 (5.58 vs. 5.15).     

Prediction of acidity in acetonitrile solution with COSMO-RS

The COSMO-RS method, a combination of the quantum chemical dielectric continuum solvation model COSMO with a statistical thermodynamics treatment for realistic solvation simulations, has been used for the prediction of pK(a) values in acetonitrile. For a variety of 93 organic acids, the directly calculated values of the free energies of dissociation in acetonitrile showed a very good correlation with the pK(a) values (r(2) = 0.97) in acetonitrile, corresponding to a standard deviation of 1.38 pK(a) units. Thus, we have a prediction method for acetonitrile pK(a) with the intercept and the slope as the only adjusted parameters. Furthermore, the pK(a) values of CH acids yielding large anions with delocalized charge can be predicted with a rmse of 1.12 pK(a) units using the theoretical values of slope and intercept resulting in truly ab initio pK(a) prediction. In contrast to our previous findings on aqueous acidity predictions the slope of the experimental pK(a) versus theoretical DeltaG(diss) was found to match the theoretical value 1/RT ln (10) very well. The predictivity of the presented method is general and is not restricted to certain compound classes. However, a systematic correction of -7.5 kcal mol(-1) is required for compounds that do not allow electron-delocalization in the dissociated anion. The prediction model was tested on a diverse test set of 129 complex multifunctional compounds from various sources, reaching a root mean square deviation of 2.10 pK(a) units.     

An accurate density functional theory based estimation of pK(a) values of polar residues combined with experimental data: from amino acids to minimal proteins

We report a scheme for estimating the acid dissociation constant (pK(a)) based on quantum-chemical calculations combined with a polarizable continuum model, where a parameter is determined for small reference molecules. We calculated the pK(a) values of variously sized molecules ranging from an amino acid to a protein consisting of 300 atoms. This scheme enabled us to derive a semiquantitative pK(a) value of specific chemical groups and discuss the influence of the surroundings on the pK(a) values. As applications, we have derived the pK(a) value of the side chain of an amino acid and almost reproduced the experimental value. By using our computing schemes, we showed the influence of hydrogen bonds on the pK(a) values in the case of tripeptides, which decreases the pK(a) value by 3.0 units for serine in comparison with those of the corresponding monopeptides. Finally, with some assumptions, we derived the pK(a) values of tyrosines and serines in chignolin and a tryptophan cage. We obtained quite different pK(a) values of adjacent serines in the tryptophan cage; the pK(a) value of the OH group of Ser13 exposed to bulk water is 14.69, whereas that of Ser14 not exposed to bulk water is 20.80 because of the internal hydrogen bonds.     

A reliable and efficient first principles-based method for predicting pKa values. III. Adding explicit water molecules: Can the theoretical slope be reproduced and pKa values predicted more accurately?

A popular method for predicting pK(a) values for organic molecules in aqueous solution is to establish empirical linear least-squares fits between calculated deprotonation energies and known experimental pK(a) values. In virtually all such calculations, the empirically observed slope of the pK(a) vs. ΔE fit is significantly less than the theoretical value, 1/(2.303RT) (which is 0.73 mol/kcal at room temperature). In our own continuum solvation calculations (Zhang et al., J Phys Chem A 2010, 114, 432), the empirical slope for carboxylic acids was only 0.23 mol/kcal, despite the excellent fit to the experimental pK(a) values. There has been much speculation about the reason for this phenomenon. Although the ΔE - pK(a) relation neglects entropic effects, these are expected to largely cancel. The most likely cause for the strange behavior of the fitted slope is explicit solute-solvent (water) interactions, especially involving the ions, which cannot be described accurately by continuum solvation models. We used our previously developed pK(a) protocol (OLYP/6-311+G(d,p)//3-21G(d) with the COSMO solvation model) to investigate the effect of adding one or two explicit water molecules to the system. The slopes for organic acids (especially carboxylic acids) are much closer to the theoretical value when explicit water molecules are added to both the neutral molecule and the anion. However, explicit water molecules have almost no effect on the slopes for organic bases. Adding explicit water molecules to the ions only produces intermediate results. Unfortunately, linear fits involving explicit water molecules have much larger errors than with continuum solvation models alone and are also much more expensive. Consequently, they are not suitable for large-scale pK(a) calculations. The results compared with literature values showed that our predicted pK(a) s are more accurate.     

DFT/electrostatic calculations of pK(a) values in cytochrome c oxidase

Using classical electrostatic calculations, earlier we examined the dependence of the protonation state of bovine cytochrome c oxidase (CcO) on its redox state. Based on these calculations, we have proposed a model of CcO proton pumping that involves His291, one of the Cu(B) histidine ligands, which was found to respond to redox changes of the enzyme Fe(a)(3)-Cu(B) catalytic center. In this work, we employ combined density functional and continuum electrostatic calculations to evaluate the pK(a)() values of His291 and Glu242, two key residues of the model. The pK(a) values are calculated for different redox states of the enzyme, and the influence of different factors on the pK(a)'s is analyzed in detail. The calculated pK(a)() values of Glu242 are between 9.4 and 12.0, depending on the redox state of the protein, which is in excellent agreement with recent experimental measurements. Assuming the reduced state of heme a(3), His291 of the oxidized Cu(B) center possesses a pK(a)() between 2.1 and 4.0, while His291 of the reduced Cu(B) center has a pK(a) above 17. The obtained results support the proposal that the His291 ligand of the Cu(B) center in CcO is a proton pump element.     

Calculation of the relative acidities and oxidation potentials of para-substituted phenols. A model for alpha-tocopherol in solution

Relative acidities (Delta pK(a)) of phenols and oxidation potentials (Delta E(ox)) of the phenoxide anions have been calculated for nine para-substituted phenols using density functional theory. Solvent effects were incorporated using the conductor-like polarisable continuum method. Using the calculated Delta pK(a) and Delta E(ox) values in a thermodynamic cycle, the DeltaBDE (bond dissociation enthalpy) of the phenols were also determined with all values calculated to within 1.5 kcal mol(-1) of experiment. The Delta pK(a) and Delta E(ox) values were calculated for 6-hydroxy-2,2,5,7,8-pentamethylchroman (HPMC), a model for alpha-tocopherol for which there are no known experimental values. The acidity of this compound is raised by 2.4 pK(a) units and lowered by -0.79 V relative to phenol with a calculated Delta BDE of -14.9 kcal mol(-1). There is a negative correlation (r(2) = 0.86) between the Delta pK(a) and the Delta BDE values. A stronger and positive correlation is found between the Delta E(ox) (r(2) = 0.98) and the Delta BDE values. Using these correlations it is uncovered that hydrogen abstraction of phenols, as measured by the Delta BDE, is driven by electron transfer rather than by proton transfer.     

pKa prediction from an ab initio bond length: part 3--benzoic acids and anilines

The prediction of pK(a) from a single ab initio bond length has been extended to provide equations for benzoic acids and anilines. The HF/6-31G(d) level of theory is used for all geometry optimisations. Similarly to phenols (Part 2 of this series of publications), the meta-/para-substituted benzoic acids can be predicted from a single model constructed from one bond length. This model had an impressive RMSEP of 0.13 pK(a) units. The prediction of ortho-substituted benzoic acids required the identification of high-correlation subsets, where the compounds in the same subset have at least one of the same (e.g. halogens, hydroxy) ortho substituent. Two pK(a) equations are provided for o-halogen benzoic acids and o-hydroxybenzoic acids, where the RMSEP values are 0.19 and 0.15 pK(a) units, respectively. Interestingly, the bond length that provided the best model differed between these two high-correlation subsets. This demonstrates the importance of investigating the most predictive bond length, which is not necessarily the bond involving the acid hydrogen. Three high-correlation subsets were identified for the ortho-substituted anilines. These were o-halogen, o-nitro and o-alkyl-substituted aniline high-correlation subsets, where the RMSEP ranged from 0.23 to 0.44 pK(a) units. The RMSEP for the meta-/para-substituted aniline model was 0.54 pK(a) units. This value exceeded our threshold of 0.50 pK(a) units and was higher than both the m-/p-benzoic acids in this work and the m-/p-phenols (RMSEP = 0.43) of Part 2. Constructing two separate models for the meta- and para- substituted anilines, where RMSEP values of 0.63 and 0.33 pK(a) units were obtained respectively, revealed it was the meta-substituted anilines that caused the large RMSEP value. For unknown reasons the RMSEP value increased with the addition of a further twenty meta-substituted anilines to this model. The C-N bond always produced the best correlations with pK(a) for all the high-correlation subsets. A higher level of theory and an ammonia probe improved the statistics only marginally for the hydroxybenzoic acid high-correlation subsets.     

Acidity and alkali metal adsorption on the SiO2-aqueous solution interface

The intrinsic deprotonation constant (pK(a(2))(int)) and the intrinsic ion exchange constants (pK(Me(+))(int)) of Li(+), Na(+), and K(+) on SiO(2) were uniquely determined at 30 degrees C by using the potentiometric titration data, the Gouy-Chapman-Stern-Grahame (GSCG) model for the structure of the electrical double-layer (edl) and the double-extrapolation method. The values of these constants were pK(a(2))(int) = 6.57, pK(Li(+))(int) = pK(Na(+))(int) = pK(K(+))(int) = 5.61. The chemical meaning of intrinsic equilibrium constants and the equality in the values of pK(Li(+))(int), pK(Na(+))(int) and pK(K(+))(int) were discussed.     

Ionization behavior of amino lipids for siRNA delivery: determination of ionization constants, SAR, and the impact of lipid pKa on cationic lipid-biomembrane interactions

Ionizable amino lipids are being pursued as an important class of materials for delivering small interfering RNA (siRNA) therapeutics, and research is being conducted to elucidate the structure-activity relationships (SAR) of these lipids. The pK(a) of cationic lipid headgroups is one of the critical physiochemical properties of interest due to the strong impact of lipid ionization on the assembly and performance of these lipids. This research focused on developing approaches that permit the rapid determination of the relevant pK(a) of the ionizable amino lipids. Two distinct approaches were investigated: (1) potentiometric titration of amino lipids dissolved in neutral surfactant micelles; and (2) pH-dependent partitioning of a fluorescent dye to cationic liposomes formulated from amino lipids. Using the approaches developed here, the pK(a) values of cationic lipids with distinct headgroups were measured and found to be significantly lower than calculated values. It was also found that lipid-lipid interaction has a strong impact on the pK(a) values of lipids. Lysis of model biomembranes by cationic lipids was used to evaluate the impact of lipid pK(a) on the interaction between cationic lipids and cell membranes. It was found that cationic lipid-biomembrane interaction depends strongly on lipid pK(a) and solution pH, and this interaction is much stronger when amino lipids are highly charged. The presence of an optimal pK(a) range of ionizable amino lipids for siRNA delivery was suggested based on these results. The pK(a) methods reported here can be used to support the SAR screen of cationic lipids for siRNA delivery, and the information revealed through studying the impact of pK(a) on the interaction between cationic lipids and cell membranes will contribute significantly to the design of more efficient siRNA delivery vehicles.     

Accurate prediction of basicity in aqueous solution with COSMO-RS

The COSMO-RS method, a combination of the quantum chemical dielectric continuum solvation model COSMO with a statistical thermodynamics treatment for realistic solvation simulations, has been used for the prediction of base pK(a) constants. For a variety of 43 organic bases the directly calculated values of the free energies of dissociation in water showed a very good correlation with experimental base pK(a) values (r2 = 0.98), corresponding to a standard deviation of 0.56 pK(a) units. Thus, we have an a priori prediction method for base pK(a) with the regression constant and the slope as only adjusted parameters. In accord with recent findings for pK(a) acidity predictions, the slope of pK(a) vs. DeltaG(diss) was significantly smaller than the theoretically expected value of 1/RTln(10). The predictivity of the presented method is general and not restricted to certain compound classes, but systematic corrections of 1 and 2 pKa units for secondary and tertiary aliphatic amines are required, respectively. The pK(a) prediction method was validated on a set of 58 complex multifunctional drug-like compounds, yielding an RMS accuracy of 0.66 pK(a) units.     

Measurement of side-chain carboxyl pK(a) values of glutamate and aspartate residues in an unfolded protein by multinuclear NMR spectroscopy

A triple-resonance NMR pulse scheme is presented for measuring aspartic and glutamic acid side-chain pK(a) values in unfolded protein states where chemical shift overlap is limiting. The experiment correlates side-chain carboxyl carbon chemical shifts of these residues with the backbone amide proton chemical shift of the following residue. The methodology is applied to an (15)N, (13)C labeled sample of the N-terminal SH3 domain of the Drosophila protein drk, which exists in equilibrium between folded (F(exch)) and unfolded (U(exch)) states under nondenaturing conditions. Residue-specific pK(a) values of side-chain carboxyl groups are presented for the first time for an unfolded protein (drk U(exch) state), determined from a pH titration. Results indicate that deviations from pK(a) values measured for model compounds are likely due to local effects, while long-range electrostatic interactions appear to be of minor importance for this protein.     

Practical calculation of molecular acidity with the aid of a reference molecule

A set of linear free energy models are presented for determining the pK(a) values of amines, alcohols, and carboxylic acids. Models are determined from a series of pK(a) predictors, taken both from traditional natural atomic orbital analysis (NAO) and from a novel approach introduced here of using a reference molecule: an ammonium ion for amines and a hydrogen sulfide molecule for alcohols and carboxylic acids. Using these reference molecules, we calculate the barrier to proton transfer and show that a number of properties associated with the transition state are correlated with the pK(a). By considering 38 predictors, we obtain a four-variable model for amines and a three-variable model for oxygen-containing compounds. The model for amines is based on 145 compounds and has a root mean squared error (RMSE) of 0.45 and R(2) = 0.98. The oxygen set has 48 molecules: RMSE = 0.26, and R(2) = 0.993. Similar, linear, and multilinear models are constructed after separating the sets into chemically similar categories: alcohols, carboxylic acids, and primary, secondary, tertiary, and aromatic amines. This separation gives simpler models with relatively low RMSE values, where the most important predictor of the pK(a) is the difference in energy between transferring the proton from the reference molecular base to the conjugate acid from the data set.     

Protonation studies of modified adenine and adenine nucleotides by theoretical calculations and (15)N NMR

The acid/base character of nucleobases affects phenomena such as self-association, interaction with metal ions, molecular recognition by proteins, and nucleic acid base-pairing. Therefore, the investigation of proton-transfer equilibria of natural and synthetic nucleos(t)ides is of great importance to obtain a deeper understanding of these phenomena. For this purpose, a set of ATP prototypes was investigated using (15)N NMR spectroscopy, and the corresponding adenine bases were investigated by theoretical calculations. (15)N NMR measurements provided not only acidity constants but also information on the protonation site(s) on the adenine ring and regarding the ratio of the singly protonated species in equilibrium. Substituents of different nature and position on the adenine ring did not change the preferred protonation site, which remained N1. However, for 2-thioether-ATP derivatives a mixed population of N1 and N7 singly protonated species was observed. Reduction of basicity of 0.4-1 pK(a) units relative to ATP was also observed for all evaluated ATP derivatives, except for 2-Cl-ATP, for which K(a) was ca. 10,000-fold lower. To explain the substitution-dependent variations in the experimental pK(a) values of the ATP analogues, gas-phase proton affinities (PA), Delta Delta G(hyd), and pK(a) values of the corresponding adenine bases were calculated using quantum mechanical methods. The computed PA and Delta Delta G(hyd) values successfully explained the experimental pK(a) values. A computational procedure for the prediction of accurate pK(a) values was developed using density functional theory and polarizable continuum model calculations. In this procedure, we developed a set of parameters for the polarizable continuum model that was fitted to reproduce experimental pK(a) values of nitrogen heterocycles. This method is proposed for the prediction of pK(a) values and protonation site(s) of purine analogues that have not been synthesized or analyzed.     

Perturbation of the NH2 pKa value of adenine in platinum(II) complexes: distinct stereochemical internucleobase effects

The degree of acidification of the exocyclic N6 amino group of the model nucleobase 9-methyladenine (9MeA) in relation to the number and site(s) of Pt(II) binding has been studied in detail. It is found that twofold Pt(II) binding to N1 and N7 lowers the pK(a) value from 16.7 in the free base to 12-8. The lowest pK(a) values are observed when the resulting N6H(-) amide group is intramolecularly stabilized by an H-bond donor such as the N6H(2) group of a suitably positioned second 9MeA ligand. Deprotonation of the N6 amino group facilitates Pt migration from N1 to N6, and subsequent reprotonation of the N1 position yields a twofold N7,N6-metalated form of the rare imino tautomer of 9MeA, which has a pK(a) value of 5.03. These findings demonstrate a principle that is of potential relevance to the topic of "shifted pK(a)" values of adenine nucleobases, which is believed to be important with regard to acid-base catalysis of RNAs at physiological pH values. The principle states that a nucleobase pK(a) value can be sufficiently lowered to reach near-neutral values and that the pK(a) value of the protonated base does not necessarily have to be increased to accomplish this effect.     

pKa of acetate in water: a computational study

Several computational methods including the conductor-like polarizable continuum model, CPCM with both UAKS and UAHF cavities, Cramer and Truhlar's generalized Born solvation model, SM5.4(AM1), SM5.4(PM3), and SM5.43R(mPW1PW91/6-31+G(d)), and mixed QM/MM-Ewald simulations were used to calculate the pK(a) values of acetate and bicarbonate anions in aqueous solution. This work provided a critical and comprehensive assessment of the quality of these theoretical models in the calculation of aqueous solvation free energies for the singly charged acetate and bicarbonate ions, as well as the doubly charged acetate dianion and carbonate dianion. It was shown that QM/MM-Ewald simulations could give an accurate and consistent evaluation of the pK(a) values of acetate and bicarbonate based on both the relative and absolute pK(a) formulas, while other methods could yield satisfactory results only for certain calculations. However, this does not mean that the current QM/MM-Ewald protocol is superior to other methods. The useful information obtained in this investigation is that both the absolute and relative pK(a) formulas should better be tested in accurate calculations of pK(a) values based on any methods.