Toward the accurate first-principles prediction of ionization equilibria in proteins.

The calculation of pK(a) values for ionizable sites in proteins has been traditionally based on numerical solutions of the Poisson-Boltzmann equation carried out using a high-resolution protein structure. In this paper, we present a method based on continuous constant pH molecular dynamics (CPHMD) simulations, which allows the first-principles description of protein ionization equilibria. Our method utilizes an improved generalized Born implicit solvent model with an approximate Debye-Hückel screening function to account for salt effects and the replica-exchange (REX) protocol for enhanced conformational and protonation state sampling. The accuracy and robustness of the present method are demonstrated by 1 ns REX-CPHMD titration simulations of 10 proteins, which exhibit anomalously large pK(a) shifts for the carboxylate and histidine side chains. The experimental pK(a) values of these proteins are reliably reproduced with a root-mean-square error ranging from 0.6 unit for proteins containing few buried ionizable side chains to 1.0 unit or slightly higher for proteins containing ionizable side chains deeply buried in the core and experiencing strong charge-charge interactions. This unprecedented level of agreement with experimental benchmarks for the de novo calculation of pK(a) values suggests that the CPHMD method is maturing into a practical tool for the quantitative prediction of protein ionization equilibria, and this, in turn, opens a door to atomistic simulations of a wide variety of pH-coupled conformational phenomena in biological macromolecules such as protein folding or misfolding, aggregation, ligand binding, membrane interaction, and catalysis.