Constant pH molecular dynamics with proton tautomerism

Abstract

The current article describes a new two-dimensional lambda-dynamics method to include proton tautomerism in continuous constant pH molecular dynamics (CPHMD) simulations. The two-dimensional lambda-dynamics framework is used to devise a tautomeric state titration model for the CPHMD simulations involving carboxyl and histidine residues. Combined with the GBSW implicit solvent model, the new method is tested on titration simulations of blocked histidine and aspartic acid as well as two benchmark proteins, turkey ovomucoid third domain (OMTKY3) and ribonuclease A (RNase A). A detailed analysis of the errors inherent to the CPHMD methodology as well as those due to the underlying solvation model is given. The average absolute error for the computed pKa values in OMTKY3 is 1.0 pK unit. In RNase A the average absolute errors for the carboxyl and histidine residues are 1.6 and 0.6 pK units, respectively. In contrast to the previous work, the new model predicts the correct sign for all the pKa shifts, but one, in the benchmark proteins. The predictions of the tautomeric states of His12 and His48 and the conformational states of His48 and His119 are in agreement with experiment. Based on the simulations of OMTKY3 and RNase A, the current work has demonstrated the capability of the CPHMD technique in revealing pH-coupled conformational dynamics of protein side chains.

FIGURE 1

Acid-base and proton tautomeric equilibria of histidine (A) and aspartic acid (B).

FIGURE 2

Two-dimensional potential of mean force (PMF) map along the titration coordinate λ and tautomeric interconversion coordinate x for the blocked His (left) and Asp (right) residues. The color bars shown on the right indicate that the PMF increases from dark blue to dark red.

FIGURE 3

2-ns standard GBSW simulation of the blocked doubly protonated histidine. The top plot shows the time series of the distance from the backbone carbonyl oxygen to ND1 (solid) and that to NE2 (shaded). The middle and bottom plots show the time series of the χ₁ and χ₂ angles, respectively.

FIGURE 4

Coupling between Asp⁷ and Glu¹⁰ in a 2-ns simulation of OMTKY3 at pH 3. The top plot shows the time series of the running average distances between Asp⁷ and Glu¹⁰ (solid), Asp⁷ and Lys³⁴ (shaded dashed), Glu¹⁰ and Lys³⁴ (dashed), and Glu¹⁰ and Arg²¹ (dotted) computed from 100-ps time windows. The definitions of the residue-residue distances are given under Table 4. The bottom plot shows the running unprotonated fraction for Asp⁷ (solid) and Glu¹⁰ (dashed) computed from 100-ps time windows.

FIGURE 5

Coupling between His¹² and His¹¹⁹ in a 2-ns simulation of RNase A at pH 5. The top plot shows the time series of the distance between NE2 of His¹² and ND1 of His¹¹⁹. The bottom plots shows the running fraction of neutral states for His¹² (solid) and His¹¹⁹ (dashed) computed from 100-ps time windows.

FIGURE 6

Snapshots from the 2-ns simulation of RNase A showing the salt-bridge and hydrogen-bond interactions involving Glu², Asp⁸³, Asp¹²¹, and CT-Val under pH conditions of 2, 0, 0, and 0, respectively.