pKa calculations in solution and proteins with QM/MM free energy perturbation simulations: a quantitative test of...

The accuracy of biological simulations depends, in large part, on the treatment of electrostatics. Due to the availability of accurate experimental values, calculation of pKa provides stringent evaluation of computational methods. The generalized solvent boundary potential (GSBP) and Ewald summation electrostatic treatments were recently implemented for combined quantum mechanical and molecular mechanics (QM/MM) simulations by our group. These approaches were tested by calculating pKa shifts due to differences in electronic structure and electrostatic environment; the shifts were determined for a series of small molecules in solution, using various electrostatic treatments, and two residues (His 31, Lys 102) in the M102K T4-lysozyme mutant with large pKa shifts, using the GSBP approach. The calculations utilized a free energy perturbation scheme with the QM/MM potential function involving the self-consistent charge density functional tight binding (SCC-DFTB) and CHARMM as the QM and MM methods, respectively. The study of small molecules demonstrated that inconsistent electrostatic models produced results that were difficult to correct in a robust manner; by contrast, extended electrostatics, GSBP, and Ewald simulations produced consistent results once a bulk solvation contribution was carefully chosen. In addition to the electrostatic treatment, the pKa shifts were also sensitive to the level of the QM method and the scheme of treating QM/MM Coulombic interactions; however, simple perturbative corrections based on SCC-DFTB/CHARMM trajectories and higher level single point energy calculations were found to give satisfactory results. Combining all factors gave a root-mean-square difference of 0.7 pKa units for the relative pKa values of the small molecules compared to experiment. For the residues in the lysozyme, an accurate pKa shift was obtained for His 31 with multiple nanosecond simulations. For Lys 102, however, the pKa shift was estimated to be too large, even after more than 10 nanosecond simulations for each lambda window; the difficulty was due to the significant, but slow, reorganization of the protein and water structure when Lys 102 was protonated. The simulations support that Lys 102 is deprotonated in the X-ray structure and the protein is highly destabilized when this residue is protonated.