Best Practices in Constant pH MD Simulations: Accuracy and Sampling

Abstract

Various approaches have been proposed to include the effect of pH in molecular dynamics (MD) simulations. Among these, the λ-dynamics approach proposed by Brooks and co-workers [Kong, X.; Brooks III, C. L. J. Chem. Phys. 1996, 105, 2414-2423] can be performed with little computational overhead and hfor each typeence be used to routinely perform MD simulations at microsecond time scales, as shown in the accompanying paper [Aho, N. et al. J. Chem. Theory Comput. 2022, DOI: 10.1021/acs.jctc.2c00516]. At such time scales, however, the accuracy of the molecular mechanics force field and the parametrization becomes critical. Here, we address these issues and provide the community with guidelines on how to set up and perform long time scale constant pH MD simulations. We found that barriers associated with the torsions of side chains in the CHARMM36m force field are too high for reaching convergence in constant pH MD simulations on microsecond time scales. To avoid the high computational cost of extending the sampling, we propose small modifications to the force field to selectively reduce the torsional barriers. We demonstrate that with such modifications we obtain converged distributions of both protonation and torsional degrees of freedom and hence consistent pK_a estimates, while the sampling of the overall configurational space accessible to proteins is unaffected as compared to normal MD simulations. We also show that the results of constant pH MD depend on the accuracy of the correction potentials. While these potentials are typically obtained by fitting a low-order polynomial to calculated free energy profiles, we find that higher order fits are essential to provide accurate and consistent results. By resolving problems in accuracy and sampling, the work described in this and the accompanying paper paves the way to the widespread application of constant pH MD beyond pK_a prediction.

Figure 1

Distributions of the λ-coordinate (A and E) and dihedrals (B–D and F–H) in constant pH MD simulation of Asp with the original (A–D) and modified CHARMM36m (E–H) force field. While the simulations were performed for an ADA tripeptide, only the central aspartic acid is shown for clarity in the inset of A. (A and D) Distributions for third-order fits for ⟨∂V/∂λ⟩_λ obtained with the original force field. Different colors correspond to independent replicas. Distributions for the λ-coordinate (A) as well as distributions for N–C_α–C_β–C_γ (B) and -C_γ-- (D) dihedrals are not identical. Right column (E–H) shows the distributions for the modified CHARMM36m force field with third-order fit for ⟨∂V/∂λ⟩_λ. Distributions are identical.

Figure 2

Distributions of dihedral angles for which the torsion potentials were corrected from standard MD simulations. (A) Dihedral distributions from publicly available trajectories^, (total simulation time ≈ 10 μs). Probability density is significant only around local minima. (B) Dihedral distributions obtained from simulations of cardiotoxin V with both the original and the modified barriers for different protonation states of the titratable residues. Local minima are preserved, and no additional configurations are observed.

Figure 3

Modification of the torsional barriers in the Asp side chain. (A) Aspartic acid and its atomic nomenclature. (B and C) Corrections added to dihedral torsions of the original force field. Corrections with two (B) and three (C) local minima were used for torsions with two and three local minima, respectively. Heights of the corrections were selected through the iterative process, aimed at achieving consistent λ-distributions without introducing additional local minima in the free energy profiles. (D–F) Original and modified torsional barriers of Asp for both protonated (H⁺) and deprotonated (H^–) states.

Figure 4

Quality of the V^MM(λ) correction potential as a function of the order of the polynomial fit to ⟨∂V/∂λ⟩_λ for Asp. (A) Fitting error as a function of fitting order (black line). Gray dashed line shows the average error in the calculated ⟨∂V^MM/∂λ⟩. (B) Distributions of λ-coordinates for the third- and seventh-order polynomial fits to ⟨∂V/∂λ⟩_λ. Whereas with the lower order fit the distribution is significantly rugged, the distribution becomes nearly flat and uniform on the [0, 1] interval of the λ-coordinate if a seventh-order fit is used.

Figure 5

Parameterization of the buffer. (A) AWH friction metric as a function of buffer charge. Since the higher the friction metric is the slower the sampling, the buffers should ideally have low charges. (B) Density distribution of buffers in a membrane system. Optimized buffers do not penetrate into the lipid bilayer, while uncharged sodium ions do. (C) Radial distribution function between buffers. When standard sodium-ion parameters are used as buffers, the tendency to form clusters is high. Optimization of buffer parameters prevents clustering. (D) Radial distribution functions of buffers with respect to the protein. With the original sodium parameters the buffers have a higher tendency to bind to the protein than with the optimized parameters.

Figure 6

Effect of the box size (A) and ionic strength (B) on the distribution of the λ-coordinate of the ADA systems (ADA, ADA₃, ADA₇, ADA_low salt, ADA_high salt).

Figure 7

Titration of cardiotoxin V. (A and B) Titration curves of the titratable residues in 1CVO₁ obtained from constant pH simulations performed with the original and modified CHARMM36 force fields, respectively. Black dots show the deprotonation ratio for the individual replica. Gray lines show the fits to the Henderson–Hasselbalch equation. Red lines are Henderson–Hasselbalch curves computed for the experimental pK_a value of the corresponding residue.^, For ASP59, the exact pK_a is not known. For each curve, the standard deviation between the calculated and the fitted deprotonation ratio, averaged over all pH values and replicas, is shown.

Figure 8

Sampling of cardiotoxin V N-terminal loop. (A) Structure of cardiotoxin V. Protein is shown in cartoon representation with HIS4, PHE10, and TYR12 shown as sticks. (B) Distributions of the distance between the centers of mass of HIS4 and TYR12 for all replicas combined. Combined distributions were calculated for constant pH MD with both the modified and the original force fields as well as for standard MD with HIS4 in all possible protonation states. (C) Distribution of HIS4 N–C_α–C_β−C_γ torsion. (D) Three observed conformations of TYR12 around HIS4. (E) Distributions of the distance between the centers of mass of HIS4 and TYR12 for individual replicas, computed for constant pH MD with the modified force field. (F) λ distributions of HIS4 protonated state for individual replicas of cardiotoxin V from the constant pH simulations with the modified force field.