Rapid calculation of protein pKa values using Rosetta

Abstract

We developed a Rosetta-based Monte Carlo method to calculate the pK(a) values of protein residues that commonly exhibit variable protonation states (Asp, Glu, Lys, His, and Tyr). We tested the technique by calculating pK(a) values for 264 residues from 34 proteins. The standard Rosetta score function, which is independent of any environmental conditions, failed to capture pK(a) shifts. After incorporating a Coulomb electrostatic potential and optimizing the solvation reference energies for pK(a) calculations, we employed a method that allowed side-chain flexibility and achieved a root mean-square deviation (RMSD) of 0.83 from experimental values (0.68 after discounting 11 predictions with an error over 2 pH units). Additional degrees of side-chain conformational freedom for the proximal residues facilitated the capture of charge-charge interactions in a few cases, resulting in an overall RMSD of 0.85 pH units. The addition of backbone flexibility increased the overall RMSD to 0.93 pH units but improved relative pK(a) predictions for proximal catalytic residues. The method also captures large pK(a) shifts of lysine and some glutamate point mutations in staphylococcal nuclease. Thus, a simple and fast method based on the Rosetta score function and limited conformational sampling produces pK(a) values that will be useful when rapid estimation is essential, such as in docking, design, and folding.

Figure 1

(a and b) Correlation between predicted and experimental pK_a values calculated using the standard Rosetta score function (a) without an explicit electrostatic potential and (b) with a distance-dependent Coulomb potential and calibrated solvation reference energies. In panel a, the prediction plots are flat (no shifts from the intrinsic pK_a values, denoted by + symbols) whereas in b, 78% of the pK_a predictions are within 1 pH unit from the experimental values (dashed lines) and only 4% of the predictions have errors > 2 pH units. Absolute pK_a values are plotted for clarity; see Fig. S2 for ΔpK_a correlation plots. Note that the RMSD values based on pK_a and ΔpK_a values are equivalent.

Figure 2

Conformational flexibility. (a) Interacting neighbor residues E17 and E26 from calbindin D9k (4ICB) are more accurately predicted when neighbor side-chain flexibility and protonation are allowed. (b) A structural model of RNase H (2RN2) generated using RosettaRelax (blue) compared with the native structure (cyan). D10 is closer to D70 in the model, resulting in a shift in the predicted pK_a value of the D10 from 4.0 using the native structure to 5.5 using the relaxed model. The experimental pK_a value for D10 is 6.1. (c) The structural model of H227 from phospholipase C (1GYM) has no space for a proton and thus highly favors the default (deprotonated) variant, resulting in a predicted pK_a of 2.7.

Figure 3

(a–c) Correlation between predicted and experimental pK_a values using (a) neighbor repacking and (b and c) a structural ensemble of 50 backbone conformations. (b) pK_a distribution for each generated structure in the ensemble. (c) Average pK_a value over the ensemble. The pK_a predictions are highly sensitive to local side-chain and backbone conformational changes.

Figure 4

Correlation between predicted and experimental pK_a values for SNase mutants with large pK_a shifts.