Skip to main page content
U.S. flag

An official website of the United States government

Dot gov

The .gov means it’s official.
Federal government websites often end in .gov or .mil. Before sharing sensitive information, make sure you’re on a federal government site.

Https

The site is secure.
The https:// ensures that you are connecting to the official website and that any information you provide is encrypted and transmitted securely.

Access keys NCBI Homepage MyNCBI Homepage Main Content Main Navigation
. 2022 Oct 11;18(10):6148-6160.
doi: 10.1021/acs.jctc.2c00516. Epub 2022 Sep 21.

Scalable Constant pH Molecular Dynamics in GROMACS

Affiliations

Scalable Constant pH Molecular Dynamics in GROMACS

Noora Aho et al. J Chem Theory Comput. .

Abstract

Molecular dynamics (MD) computer simulations are used routinely to compute atomistic trajectories of complex systems. Systems are simulated in various ensembles, depending on the experimental conditions one aims to mimic. While constant energy, temperature, volume, and pressure are rather straightforward to model, pH, which is an equally important parameter in experiments, is more difficult to account for in simulations. Although a constant pH algorithm based on the λ-dynamics approach by Brooks and co-workers [Kong, X.; Brooks III, C. L. J. Chem. Phys.1996, 105, 2414-2423] was implemented in a fork of the GROMACS molecular dynamics program, uptake has been rather limited, presumably due to the poor scaling of that code with respect to the number of titratable sites. To overcome this limitation, we implemented an alternative scheme for interpolating the Hamiltonians of the protonation states that makes the constant pH molecular dynamics simulations almost as fast as a normal MD simulation with GROMACS. In addition, we implemented a simpler scheme, called multisite representation, for modeling side chains with multiple titratable sites, such as imidazole rings. This scheme, which is based on constraining the sum of the λ-coordinates, not only reduces the complexity associated with parametrizing the intramolecular interactions between the sites but also is easily extendable to other molecules with multiple titratable sites. With the combination of a more efficient interpolation scheme and multisite representation of titratable groups, we anticipate a rapid uptake of constant pH molecular dynamics simulations within the GROMACS user community.

PubMed Disclaimer

Conflict of interest statement

The authors declare no competing financial interest.

Figures

Figure 1
Figure 1
Illustration of the additional λ-dependent potential energy terms (B–D). Panel A shows the protonation free energy of a titratable residue in its reference state obtained at the force field level, ΔGiMMi) . To compensate for the shortcomings of the force field and obtain a zero free energy difference between the two protonation states (A + B), we add the correction potential, ViMM(λ), shown in panel B. A biasing potential, Vibias(λ), is introduced to avoid sampling of nonphysical states (panel C). To model the proton chemical potential (pH), we add a pH-dependent term, VpHi) (panel D). For a titratable residue at pH ≠ pKa, the total λ-dependent potential, including the interpolated force field functions and the three additional terms, is illustrated in panel (A + B + C + D).
Figure 2
Figure 2
Multisite representation illustrated for a histidine side chain. Each possible protonation state is represented by its own λ-coordinate. HSP refers to doubly protonated histidine, HSD, and HSE to histidine singly protonated at the Nδ or Nϵ, respectively. HS0 is a common, nonphysical state of the residue. To restrict the sampling to the plane connecting the physical states, a constraint λ1 + λ2 + λ3 = 1 is applied (gray plane). A biasing potential is also applied to enhance sampling at the end states, where one of the λ-coordinates is close to one, while the other coordinates are close to zero.
Figure 3
Figure 3
Titrations of tripeptide amino acids (Glu, Asp, and His) in water. The top and bottom rows show titrations performed with the modified AA CHARMM36 and CG Martini force fields, respectively. In all simulations, neutrality was maintained by including ten buffer particles in combination with the charge constraint. Dots show the fraction of frames in which the residue is deprotonated, and the dashed lines represent the fits to the Henderson–Hasselbalch equation. For His, the blue color represents the macroscopic pKa, while yellow and red represent the microscopic pKa values for HSD (proton on Nδ) and HSE (proton on Nϵ), respectively. In the Martini 2.0 model, HSD and HSE are indistinguishable and hence only the macroscopic titration curve is shown. Errors of Sdeprot were estimated from the standard error of the mean for the different replicas. From the fits, the pKa values were estimated and listed in Table 1.
Figure 4
Figure 4
Titration curves of the cardiotoxin V protein obtained from constant pH MD simulations with the CHARMM36 (top) and Martini 2.0 force fields (bottom). For each of the four titratable residues in this protein, the dots show the fraction of frames in which the residue is deprotonated. Errors of Sdeprot were estimated from the standard error of the mean for the different replicas. The lines show the best fits to the Henderson–Hasselbalch equation. The pKa values for each titratable residue were estimated from these fits and listed in Table 1. The right panel shows the protein structure with the four titratable residues highlighted in stick representation.
Figure 5
Figure 5
Titration curves of the HEWL protein obtained from constant pH MD simulations with the CHARMM36 (top) and Martini 2.0 force fields (bottom). For each of the ten titratable residues, the dots show the fraction of frames in which that residue is deprotonated. Errors of Sdeprot were estimated from the standard error of the mean for the different replicas. The lines show the best fit to the Henderson-Hasselbach equation. The pKa values for each titratable residue were estimated from these fits and listed in Table 1. The right panel shows the protein structure with the ten titratable residues highlighted in stick representation.
Figure 6
Figure 6
(A) Relative performance of interpolating potentials in a previous implementation of CpHMD into a fork of GROMACS 3.3 release (blue) and of charge interpolation in our new implementation (red) as a function of the number of titratable sites. The simulations were performed for the turkey ovomucoid inhibitor protein (PDB ID: 2GKR(53)), shown in the inset, where the titratable sites are highlighted in stick representation. (B) Comparison of the performance between CPU-only and CPU+GPU implementations for the ligand-gated ion channel GLIC (PDB ID: 4HFI(52)) with 185 titratable sites. In total, the GLIC system contained 292135 atoms.

Similar articles

Cited by

References

    1. Hollingsworth S. A.; Dror R. O. Molecular Dynamics Simulation for All. Neuron 2018, 99, 1129–1143. 10.1016/j.neuron.2018.08.011. - DOI - PMC - PubMed
    1. Groenhof G.; Modi V.; Morozov D. Observe while it happens: catching photoactive proteins in the act with non-adiabatic molecular dynamics simulations. Curr. Opin. Struct. Biol. 2020, 61, 106–112. 10.1016/j.sbi.2019.12.013. - DOI - PubMed
    1. Warshel A.; Sussman F.; King G. Free energy of charges in solvated proteins: microscopic calculations using a reversible charging process. Biochemistry 1986, 25, 8368–8372. 10.1021/bi00374a006. - DOI - PubMed
    1. Alexov E.; Mehler E. L.; Baker N.; Baptista A. M.; Huang Y.; Milletti F.; Erik Nielsen J.; Farrell D.; Carstensen T.; Olsson M. H. M.; Shen J. K.; Warwicker J.; Williams S.; Word J. M. Progress in the prediction of pKa values in proteins. Proteins: Struct., Funct., Bioinf. 2011, 79, 3260–3275. 10.1002/prot.23189. - DOI - PMC - PubMed
    1. Chen W.; Morrow B. H.; Shi C.; Shen J. K. Recent development and application of constant pH molecular dynamics. Mol. Simul. 2014, 40, 830–838. 10.1080/08927022.2014.907492. - DOI - PMC - PubMed

LinkOut - more resources