Theoretical and Computational Chemistry

Best practices in constant pH MD simulations: accuracy and sampling



Various approaches have been proposed to include the effect of pH in Molecular Dynamics (MD) simulations. Among these, the lambda-dynamics approach proposed by Brooks and co-workers can be performed with little computational overhead and hence be used to routinely perform MD simulations at microsecond timescales, as shown in the accompanying paper. At such timescales, however, the accuracy of the molecular mechanics force field and the parameterization becomes critical. Here, we address these issues and provide the community with guidelines on how to set up and perform long timescale 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 timescales. 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 pKa 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.


Thumbnail image of 220517_constant_ph_best_practices.pdf

Supplementary material

Thumbnail image of 220517_constant_ph_best_practices_si.pdf
Best practices in constant pH MD SI
lambda- and angle distributions for Glu, His, Lys, C- and N-termini, torsion corrections for Glu, His, and C-terminus, effects of force field modifications on dihedral dynamics and energy profiles, effects of fitting order on lambda-distributions, and conformations of HIS4 of cardiotoxin V are presented in SI.
Thumbnail image of
Zip file with all the parameters used in the simulation
In this file, all the information needed to reproduce our simulations is provided: force field parameters, input files, etc.

Supplementary weblinks

Gromacs constant pH gitlab page
The GROMACS constant pH MD code, manual and tutorials