Protein pK_a calculations

Some methods are based on solutions to the Poisson–Boltzmann equation (PBE), often referred to as FDPB-based methods (FDPB stands for "finite difference Poisson–Boltzmann"). The PBE is a modification of Poisson's equation that incorporates a description of the effect of solvent ions on the electrostatic field around a molecule.

The H++ web server,[7] the pKD webserver,[8] MCCE2, Karlsberg+, PETIT and GMCT use the FDPB method to compute pK_a values of amino acid side chains.

FDPB-based methods calculate the change in the pK_a value of an amino acid side chain when that side chain is moved from a hypothetical fully solvated state to its position in the protein. To perform such a calculation, one needs theoretical methods that can calculate the effect of the protein interior on a pK_a value, and knowledge of the pKa values of amino acid side chains in their fully solvated states.[2][3][4][5]

A set of empirical rules relating the protein structure to the pK_a values of ionizable residues have been developed by Li, Robertson, and Jensen.[9] These rules form the basis for the web-accessible program called PROPKA for rapid predictions of pK_a values. A recent empirical pK_a prediction program was released by Tan KP et.al. with the online server DEPTH web server.[10]

Molecular dynamics methods of calculating pK_a values make it possible to include full flexibility of the titrated molecule.[11][12][13]

Molecular dynamics based methods are typically much more computationally expensive, and not necessarily more accurate, ways to predict pK_a values than approaches based on the Poisson–Boltzmann equation. Limited conformational flexibility can also be realized within a continuum electrostatics approach, e.g., for considering multiple amino acid sidechain rotamers. In addition, current commonly used molecular force fields do not take electronic polarizability into account, which could be an important property in determining protonation energies.

From the titration of protonatable group, one can read the so-called pK_a^1⁄2 which is equal to the pH value where the group is half-protonated (i.e. when 50% such groups would be protonated). The pK_a^1⁄2 is equal to the Henderson–Hasselbalch pK_a (pK^HH
_a) if the titration curve follows the Henderson–Hasselbalch equation.[14] Most pK_a calculation methods silently assume that all titration curves are Henderson–Hasselbalch shaped, and pK_a values in pK_a calculation programs are therefore often determined in this way. In the general case of multiple interacting protonatable sites, the pK_a^1⁄2 value is not thermodynamically meaningful. In contrast, the Henderson–Hasselbalch pK_a value can be computed from the protonation free energy via

${\displaystyle \mathrm {p} K_{\mathrm {a} }^{\mathrm {HH} }(\mathrm {pH} )=\mathrm {pH} -{\frac {\Delta G^{\mathrm {prot} }(\mathrm {pH} )}{\mathrm {RT} \ln 10}}}$

and is thus in turn related to the protonation free energy of the site via

${\displaystyle \Delta G^{\mathrm {prot} }(\mathrm {pH} )=\mathrm {RT} \ln 10\;(\mathrm {pH} -\mathrm {p} K_{\mathrm {a} }^{\mathrm {HH} })}$

The protonation free energy can in principle be computed from the protonation probability of the group ⟨x⟩(pH) which can be read from its titration curve

${\displaystyle \Delta G^{\mathrm {prot} }(\mathrm {pH} )=-\mathrm {RT} \ln \left[{\frac {\langle x\rangle }{1-\langle x\rangle }}\right]}$

Titration curves can be computed within a continuum electrostatics approach with formally exact but more elaborate analytical or Monte Carlo (MC) methods, or inexact but fast approximate methods. MC methods that have been used to compute titration curves[15] are Metropolis MC[16][17] or Wang–Landau MC.[18] Approximate methods that use a mean-field approach for computing titration curves are the Tanford–Roxby method and hybrids of this method that combine an exact statistical mechanics treatment within clusters of strongly interacting sites with a mean-field treatment of intercluster interactions.[19][20][21][22][23]

In practice, it can be difficult to obtain statistically converged and accurate protonation free energies from titration curves if ⟨x⟩ is close to a value of 1 or 0. In this case, one can use various free energy calculation methods to obtain the protonation free energy[15] such as biased Metropolis MC,[24] free-energy perturbation,[25][26] thermodynamic integration,[27][28][29] the non-equilibrium work method[30] or the Bennett acceptance ratio method.[31]

Note that the pK^HH
_a value does in general depend on the pH value.[32]

This dependence is small for weakly interacting groups like well solvated amino acid side chains on the protein surface, but can be large for strongly interacting groups like those buried in enzyme active sites or integral membrane proteins.[33][34][35]