GPU-Accelerated All-Atom Particle-Mesh Ewald Continuous Constant pH Molecular Dynamics in Amber

Abstract

Constant pH molecular dynamics (MD) simulations sample protonation states on the fly according to the conformational environment and user specified pH conditions; however, the current accuracy is limited due to the use of implicit-solvent models or a hybrid solvent scheme. Here, we report the first GPU-accelerated implementation, parametrization, and validation of the all-atom continuous constant pH MD (CpHMD) method with particle-mesh Ewald (PME) electrostatics in the Amber22 pmemd.cuda engine. The titration parameters for Asp, Glu, His, Cys, and Lys were derived for the CHARMM c22 and Amber ff14sb and ff19sb force fields. We then evaluated the PME-CpHMD method using the asynchronous pH replica-exchange titration simulations with the c22 force field for six benchmark proteins, including BBL, hen egg white lysozyme (HEWL), staphylococcal nuclease (SNase), thioredoxin, ribonuclease A (RNaseA), and human muscle creatine kinase (HMCK). The root-mean-square deviation from the experimental pK_a's of Asp, Glu, His, and Cys is 0.76 pH units, and the Pearson's correlation coefficient for the pK_a shifts with respect to model values is 0.80. We demonstrated that a finite-size correction or much enlarged simulation box size can remove a systematic error of the calculated pK_a's and improve agreement with experiment. Importantly, the simulations captured the relevant biology in several challenging cases, e.g., the titration order of the catalytic dyad Glu35/Asp52 in HEWL and the coupled residues Asp19/Asp21 in SNase, the large pK_a upshift of the deeply buried catalytic Asp26 in thioredoxin, and the large pK_a downshift of the deeply buried catalytic Cys283 in HMCK. We anticipate that PME-CpHMD will offer proper pH control to improve the accuracies of MD simulations and enable mechanistic studies of proton-coupled dynamical processes that are ubiquitous in biology but remain poorly understood due to the lack of experimental tools and limitation of current MD simulations.

Figure 1:. Nonlinear fitting of the mean forces to obtain the model PMF parameters for His titration.

a) and b) Fitting $\langle \partial U/\partial \theta \rangle$ at ${\theta }_x=0$ (a) or at ${\theta }_x=\pi /2$ (b) to $2A_0({\sin }^2\theta -B_0)\sin 2\theta$ gives $A_0$ and $B_0$ (a) or $A_1$ and $B_1$ (b), respectively. c) Fitting $\langle \partial U/\partial {\theta }_x\rangle$ at $\theta =\pi /2$ to $2A_{10}\left({\sin }^2{\theta }_x-B_{10}\right)\sin 2{\theta }_x$ gives $A_{10}$ and $B_{10}$. The fitting equations are derivatives of Eqs. 20, 21, and 22. The red curves are the best fits.

Figure 2:. Simulated titration plots of model peptides ACEAAXAANH₂ (X=Asp, Glu, His, Cys, and Lys) at independent pH conditions.

Top panel: unprotonated fractions of Asp, Glu, and His at different pH. Bottom panel: unprotonated fractions of Cys and Lys at different pH. At each pH, three simulation runs were performed starting from different initial velocity seeds. The $\text{p}{K}_{\text{a}}$, Hill coefficient (n), and fitting error are given. The boot strap errors are given in Table 2. The fitting was performed on all data points using the generalized Henderson-Hasselbalch equation. Performing the fits against the Henderson-Hasselbalch equation yields identical $\text{p}{K}_{\text{a}}$ values and error estimates.

Figure 3:. Comparison between the calculated and experimental pKa′s and pKa shifts of the benchmark proteins.

a) Calculated $\text{p}{K_{\text{a}}}^{\prime }\text{s}$ vs. experimental $\text{p}{K_{\text{a}}}^{\prime }\text{s}$. b) Calculated vs. experimental $\text{p}{K}_{\text{a}}$ shifts with respect to the experimental model peptide $\text{p}{K_{\text{a}}}^{\prime }\text{s}$ (Table 1). The data for Asp, Glu, His, and Cys are shown in magenta, cyan, blue, and orange, respectively. Pearson’s correlation coefficient $(r)$ and RMSE are given. The solid black lines represent the linear regression. The shaded region indicates the calculated $\text{p}{K_{\text{a}}}^{\prime }\text{s}$ within the overall RMSE (0.76 units) of the experimental values. To guide the eye, the dashed diagonal line $(x=y)$ is shown. c) Histograms of the deviations between the calculated and experimental $\text{p}{K_{\text{a}}}^{\prime }\text{s}$ for Asp (left), Glu (middle), and His (right) residues.

Figure 4:. Protonation of His166 in BBL is correlated with the pH-dependent increase in solvent exposure.

a) Structure of BBL protein with His166 explicitly shown (PDB ID: 1W4H). b) Deprotonated fraction of His166 at different pH. c) Buried ratio of His166 at different pH, defined as 1-fSASA, where fSASA (fraction of solvent accessible surface area) was calculated as SASA of the sidechain atoms relative to that in the model pentapeptide.

Figure 5:. Factors influencing the pKa′s of the catalytic dyad in HEWL.

(a) A representative snapshot at pH 7.5 showing the h-bonding and salt bridge environment of Glu35 and Asp52. (b) The unprotonation fractions of Glu35 and Asp52 at different pH. (c) Fractions of the sidechain solvent exposure (SASA value relative to the model pentapeptide) of Glu35 and Asp52 at different pH. (d) The h-bond and salt bridge occupancies of Glu35 and Asp52 at different pH.

Figure 6:. Linked titration of Asp19 and Asp21 in SNase.

(a) A snapshot from the pH 4 simulation showing the h-bonding environment of Asp19 and Asp21. (b) Total number of protons of Asp19/Asp21 at different pH. The stepwise $\text{p}{K_{\text{a}}}^{\prime }\text{s}$ are obtained from the best fit (black curve) to the two-proton coupled equation (Eq. 27). (c) The pH-dependent probabilities of four microscopic states: two protons (HH, red); proton on Asp19 (H−, magenta) or Asp21 (−H, cyan); zero proton (—, blue).

Figure 7:. Protonation of His12 in RNaseA is correlated with the decreased solvent exclusion and hydrogen bonding.

a) A zoomed-in view of the hydrogen bonding between His12 and Asn11 in RNase A. The snapshot was taken from the simulation at pH 7. b) Unprotonated fraction of His12 at different pH. c) Buried fraction of His12 at different pH. Definition of the buried fraction is given in the caption of Fig. 4. d) Occupancy of the h-bond between His12 and Asn11 at different pH.

Figure 8:. Factors influencing the pKa of Cys283 in HMCK.

(a) A zoomed-in view of the h-bond environment of Cys283 (from simulation at pH 7.5). Cys283 thiolate can form h-bonds with the backbones and sidechains of Ser285 and Asn286. (b) Unprotonated fraction of Cys283 at different pH. (c) Exposure fraction (SASA relative to that of the model pentapeptide) at different pH. d) Occupancy of the total h-bond formation of Cys283 thiolate with Ser285 and Asn286 at different pH.

Figure 9:. Effect of box size on the calculated pKa′s.

The errors of the raw (a) and finite-size corrected (b) $\text{p}{K}_{\text{a}}$ of SNase with different solvent cushion spaces, 10 (magenta), 12 (orange), 14 (green), and 18 $\text{Å}$ (cyan). The corresponding RMSE values are shown next to the legends.