Gromacs

Constant pH Simulation

    Table of contents
    1. 1. Notes

    First, check your assumptions. The ideal situation would be to have a constant-pH algorithm to perform the MD simulations. However, conventional explicit-solvent MD algorithms cannot do that, since they're constant-H+ algorithms. Furthermore, if you use a MM force field, you will never get protonation/deprotonation of sites in your solute: free H+ will just wander around, while titrable sites are stuck with their initial protonation states; if you want to overcome this you have to use a non-MM Hamiltonian which explicitly allows for proton transfer. In any case, regardless of the Hamiltonian used, the amount of H+ would never change.

    If you want to do an explicit-solvent constant-pH simulation, the simplest thing to do is to stick to conventional MD and just select the particular amount of H+ that is "typical" for your system. You will need to use orders of magnitude more water than are normally used in an MD simulation. Basically, you can neglect free H+ concentration at most pH values (even at pH=4) its concentration is orders of magnitude below that of other counterions (eg, Na+, Cl-). Furthermore, as noted above, a pure MM Hamiltonian will be unable to move protons between titrable sites (eg, unlike EVB), so you actually need to choose the protonation state of each titrable site in your molecule.

    You can get a good guess at a starting configuration by performing a standard pKa calculation with your starting structure; these are standard procedures, mostly based on continuum electrostatics for computing protonation free energies (eg, using MEAD, UHBD, DelPhi, APBS) and Monte Carlo to perform sampling of protonation states (eg, using REDTI, PETIT). Unfortunately, as the conformation of the solute changes along the MD simulation, the protonation states may become inadequate, due to the strong protonation-conformation coupling that exists in many cases. Several solutions have been proposed, but most of them are more or less heuristic attempts. The truly satisfactory solution would be to have a constant-pH MD method, and some have been proposed in the last years. This has been discussed on the GROMACS users mailing list. (Note that there is a serious theoretical problem with a method by Phil Hunenberger.) Unfortunately, constant-pH methods are recent and still under development and/or testing. Hopefully, they will become standard methods in the near future. Maybe you could one of these days just specify "pH = 7.0" in the GROMACS .mdp file and see protonation states changing during the MD run! Until then, the best solution is probably to one mentioned above: make an initial good guess with a standard pKa calculation, and eventually check it later with snapshots from the MD run.

    There are constant pH simulation models, developed by Charlie Brooks and others where you can simulate transfer of protons from one side chain to another. They only work in implicit solvent so far (and in CHARMM at that).

    There are several other constant-pH MD methods, some of them using explicit solvent. Antonio Baptista's group has developed a constant-pH MD method based on stochastic protonation changes (J. Chem. Phys. (2002) 117:4184). Although the method uses a Poisson-Boltzmann method to periodically change the protonation states, the MM/MD simulations are done with explicit solvent.

    They have actually implementated this stochastic constant-pH MD method using GROMACS (J. Phys. Chem. B (2006) 110:2927), basically following a stop-and-go approach using bash and awk scripts to interface GROMACS with MEAD (a PB solver by Don Bashford) and MCRP (a in-house program for Monte Carlo sampling of protonation states). Unfortunately, the whole thing is still too messy and hard-wired at some places, which makes it unsuitable to be submitted as a contribution to GROMACS, at least for now.

    Another explicit-solvent method, based only on MM/MD, was proposed by Hunenberger (J. Chem. Phys. (2001) 114:9706), but its theoretical basis seems to be wrong (J. Chem. Phys. (2002) 116:7766). As far as I know, all other constant-pH MD methods proposed so far used indeed implicit solvent. Stochastic approaches using implicit solvent were used by McCammon (J. Comput. Chem. (2004) 25:2038) and by Antosiewicz (Phys. Rev. E (2002) 66:051911), and Baptista's group has also proposed a fractional-charge approach (Proteins (1997) 27:523) where implicit solvent was used for computational speed; all these methods depend on some sort of simplified electrostatics-oriented method (Poisson-Boltzmann, generalized Born, etc.) to perform the protonation calculations. Implicit solvent was also used in a method proposed by Brooks (Proteins (2004) 56:738), following a different but theoretically vague approach to include protonation effects.

    For discussion of water models that can dissociate, see this paper.

     

    Notes

    The above was pasted together from various emails on the GROMACS users mailing list, but largely from some by Antonio Baptista.

    Page last modified 21:46, 9 Oct 2009 by JLemkul?