p<i>K</i>_a Values in Proteins Determined by Electrostatics Applied to Molecular Dynamics Trajectories

For a benchmark set of 194 measured p<i>K</i>_a values in 13 proteins, electrostatic energy computations are performed in which p<i>K</i>_a values are computed by solving the Poisson–Boltzmann equation. In contrast to the previous approach of Karlsberg⁺ (KB⁺) that essentially used protein crystal structures with variations in their side chain conformations, the present approach (KB2⁺MD) uses protein conformations from four molecular dynamics (MD) simulations of 10 ns each. These MD simulations are performed with different specific but fixed protonation patterns, selected to sample the conformational space for the different protonation patterns faithfully. The root-mean-square deviation between computed and measured p<i>K</i>_a values (p<i>K</i>_a RMSD) is shown to be reduced from 1.17 pH units using KB⁺ to 0.96 pH units using KB2⁺MD. The p<i>K</i>_a RMSD can be further reduced to 0.79 pH units, if each conformation is energy-minimized with a dielectric constant of ε_min = 4 prior to calculating the electrostatic energy. The electrostatic energy expressions upon which the computations are based have been reformulated such that they do not involve terms that mix protein and solvent environment contributions and no thermodynamic cycle is needed. As a consequence, conformations of the titratable residues can be treated independently in the protein and solvent environments. In addition, the energy terms used here avoid the so-called intrinsic p<i>K</i>_a and can therefore be interpreted without reference to arbitrary protonation states and conformations