ORIGINAL PAPER Benchmarking pK a prediction methods for Lys115 in acetoacetate decarboxylase Yuli Liu 1 & Anand H. G. Patel 1 & Steven K. Burger 1 & Paul W. Ayers 1 Received: 1 December 2016 /Accepted: 17 March 2017 # Springer-Verlag Berlin Heidelberg 2017 Abstract Three different pK a prediction methods were used to calculate the pK a of Lys115 in acetoacetate decarboxylase (AADase): the empirical method PROPKA, the multiconformation continuum electrostatics (MCCE) method, and the molecular dynamics/thermodynamic integration (MD/ TI) method with implicit solvent. As expected, accurate pK a prediction of Lys115 depends on the protonation patterns of other ionizable groups, especially the nearby Glu76. However, since the prediction methods do not explicitly sam- ple the protonation patterns of nearby residues, this must be done manually. When Glu76 is deprotonated, all three methods give an incorrect pK a value for Lys115. If protonated, Glu76 is used in an MD/TI calculation, the pK a of Lys115 is predicted to be 5.3, which agrees well with the experimental value of 5.9. This result agrees with previous site-directed mutagenesis studies, where the mutation of Glu76 (negative charge when deprotonated) to Gln (neutral) causes no change in K m , suggesting that Glu76 has no effect on the pK a shift of Lys115. Thus, we postulate that the pK a of Glu76 is also shifted so that Glu76 is protonated (neutral) in AADase. Keywords Protein pKa . Acetoacetate decarboxylase . Net-event kinetic Monte Carlo Introduction Of the many factors that can affect enzyme substrate binding and catalysis, the protonation state of its ionizable residues is quite significant. Thus, when studying enzymatic mecha- nisms, it is generally useful to have accurate pK a values for the active site residues. Although the pK a values of the free amino acids in solution are well known, the influence of the protein environment, including the protonation state of other surrounding residues, can change proton affinities dramatical- ly. To predict residue pK a values, there are several classes of available computational methods, including empirically based methods [1–3], Poisson-Boltzmann (PB) equation based methods [4–6], and molecular dynamics (MD) or quantum mechanics/molecular mechanics (QM/MM) free energy sim- ulation methods [7]. Empirically based methods, such as PROPKA, rely on pa- rameterization of factors that are believed to influence protein pK a values using experimentally available pK a s[1, 2]. PROPKA 2.0 determines the pK a of ionizable groups empir- ically by parameterizing the following determinants: global (GlobalDes) and local (LocalDes) desolvation effects, side- chain (SC-HB) and backbone (BB-HB) hydrogen bonds, and Coulomb interactions between charged groups (chg- chg), as shown in Eq. (1): ΔpK a ¼ ΔpK GlobalDes þ ΔpK LocalDes þ ΔpK SC-HB þ ΔpK BB-HB þ ΔpK chg-chg ð1Þ Although they are very fast, the accuracy of these methods depends strongly on the fitting procedure and the choice of reference pK a values used in parameterization. Thus, these methods can be unreliable for residues in environments sig- nificantly different from those in the fit set. This paper belongs to Topical Collection Festschrift in Honor of Henry Chermette Electronic supplementary material The online version of this article (doi:10.1007/s00894-017-3324-x) contains supplementary material, which is available to authorized users. * Paul W. Ayers ayers@mcmaster.ca; ayers@chemistry.mcmaster.ca 1 Department of Chemistry & Chemical Biology, McMaster University, Hamilton, Ontario, Canada J Mol Model (2017) 23:155 DOI 10.1007/s00894-017-3324-x