Constant pH Dynamics

During my 4-week visit as an OpenMM Visiting Scholar, I implemented a constant pH algorithm for implicit solvent.  A brief overview of the algorithm is given in this video:

Below, I describe the algorithm in much greater detail and also highlight key features of OpenMM that enabled me to quickly implement the algorithm. 


Constant pH dynamics is a simulation procedure whereby the protonation states of residues and molecules are allowed to fluctuate.  The role that this process plays in structural biology can be surprisingly difficult to model with classical molecular dynamics.  Jacobson et al [1] provide an excellent review of some of these cases.  Protonation states often are inferred with programs like MCCE [2], but in a static structural context.  With this approach, a set of stable initial protonation states can then be studied and compared, either by docking or simulation.

As noted by Jacobson, however, pH plays a much more dynamic role.  Notable examples include DNA binding, allostery, fibril formation, and other mechanisms.  These systems can exhibit cascades of proton transport, requiring an in situ approach to the simulation.  To account for these effects properly remains one of the great challenges in biomolecular simulation.  Many useful approaches have been proposed [3-6], and we propose to build on this impressive body of work.

Implementation of Constant pH Dynamics

In my implementation of constant pH dynamics, I built from an existing alchemical formulation that uses a fractional parameter  to describe the protonation state of the system.  Using this formulation, we can gradually morph molecules from one protonation state to another.  The gradual parameter adjustment allows the surrounding environment to adjust to the new protonation state.  The new state is then accepted or rejected based on a Monte Carlo criterion.  Below, I describe the new energy formulation and the Monte Carlo criterion used for accepting the trial protonation states.  

Energetics:  Consider a system with coordinates $r$ and protonation state(s) ${\lambda}$ , and an energy $U_{SYS}$, which includes bonded and nonbonded terms, as well as solvation.  The full energy term $U_T$ includes the Henderson-Hasselbach energy, giving

$$U_{T}(r,\lambda)=U_{SYS}\left(r,p(\lambda) \right) + G_{HH}(\lambda)$$

where $G_{HH}(\lambda)=\sum_{i} G_{HH}(\lambda_i) $ is the sum over all titratable residues and molecules, with the term for the $i$th residue follow the form given by Brooks [4] as

$$\beta G_{HH}(\lambda_i))=2.3\cdot \beta \lambda_i \cdot \left( pK_a(i)-pH \right )$$

The energy of the system is currently computed using a linear mixing rule for the parameter states, computed as 

$$p(\lambda)=(1-\lambda)p_0 + \lambda p_1 $$

Where $p_0$ and $p_1$ are the initial and final parameter states, respectively, which can include partial charges, Van der Waals radii, and other forcefield parameters.  The current implementation is in implicit solvent, and explicit solvent implementations will be available in the near term.

Sampling:  The perturbation procedure is a Monte Carlo technique that collaborators and I have developed known as Nonequilibrium Candidate Monte Carlo (NCMC)[7].  The idea behind this procedure is to generate a trial move that gradually adjusts the  parameter to allow the local environment to relax to the new parameter state.  We then accept this trial with Monte Carlo criterion that ensures correct statistics.  If we use a Verlet integrator in the trial move, the Metropolis acceptance criterion is

$$\textrm{acc}(r_0 \rightarrow r_1)=\min(1,e^{-
\beta \Delta H_T}) ,$$

where $\Delta H_T$ is the change in the total Hamiltonian moving from state 0 to state 1 through the nonequilibrium trial.

OpenMM Experience and Tips

Implementing this algorithm in OpenMM was very straightforward, given that the data structures can readily be accessed through the Python API.  This type of prototype development can often be very difficult in large molecular dynamics codes, but OpenMM is clearly designed with rapid prototyping in mind. It took me about a week to refresh my Python skills and learn OpenMM, and another couple of weeks to implement the algorithm using OpenMM’s Python API and begin testing the algorithm.  We still have much to do, but the basic machinery seems to be in place and working as expected, including the GPU acceleration.

For the scientist whose interest is primarily developing prototypes and testing hypotheses, the Python API is an ideal development platform.  Python is very easy to use and debug.  It is also a very legible language, so that others can easily see your logic.  If you leverage the object-oriented structure, it is easy to see how your prototype code might be implemented in C++ once the prototype is stable.

Since the Python layer is an API to OpenMM’s GPU-accelerated code, it also runs extremely fast, so that generating validation data is much easier.  We found that using special functions that allow for the rapid updating of forcefield parameters without needing to fully update the context on the GPU was extremely useful for our particular algorithm, as we needed to make frequent communications to the GPU.

The development team has provided extremely high quality documentation with great examples in the user manual.  The code is also very well written and organized, and I found the online class diagrams to be particularly helpful when looking for routines to incorporate in the code.

Future Work

I hope to continue to develop this code as part of my research program moving forward.   Having a stable constant pH code for implicit and explicit solvent would be our first effort.  Another goal is to port this procedure onto the C++ layer so that it can be readily accessible to Folding@Home, which generates many simulation trajectories that can be pieced together using Markov State models via MSMBuilder.  This would allow for microsecond to millisecond timescale studies of pH dependent folding events (for example).  OpenMM has made it easy to test out my idea, and I am excited to see what new research results come of it.


1.            Schönichen, A., et al., Considering protonation as a posttranslational modification regulating protein structure and function. Annual review of biophysics, 2013. 42: p. 289-314.

2.            Gunner, M., X. Zhu, and M.C. Klein, MCCE analysis of the pKas of introduced buried acids and bases in staphylococcal nuclease. Proteins: Structure, Function, and Bioinformatics, 2011. 79(12): p. 3306-3319.

3.            Bürgi, R., P.A. Kollman, and W.F. van Gunsteren, Simulating proteins at constant pH: an approach combining molecular dynamics and Monte Carlo simulation. Proteins: Structure, Function, and Bioinformatics, 2002. 47(4): p. 469-480.

4.            Khandogin, J. and C.L. Brooks, Constant pH molecular dynamics with proton tautomerism. Biophysical journal, 2005. 89(1): p. 141-157.

5.            Mertz, J.E. and B.M. Pettitt, Molecular dynamics at a constant pH. International Journal of High Performance Computing Applications, 1994. 8(1): p. 47-53.

6.            Mongan, J. and D.A. Case, Biomolecular simulations at constant pH. Current opinion in structural biology, 2005. 15(2): p. 157-163.

7.            Nilmeier, J.P., et al., Nonequilibrium candidate Monte Carlo is an efficient tool for equilibrium simulation. Proceedings of the National Academy of Sciences, 2011. 108(45): p. E1009-E1018.

8.            Kong, X. and C.L. Brooks, λ‐dynamics: A new approach to free energy calculations. The Journal of chemical physics, 1996. 105(6): p. 2414-2423.

9.            Stern, H.A., Molecular simulation with variable protonation states at constant pH. The Journal of chemical physics, 2007. 126(16): p. 164112.

10.          Wagoner, J.A. and V.S. Pande, A smoothly decoupled particle interface: New methods for coupling explicit and implicit solvent. The Journal of chemical physics, 2011. 134(21): p. 214103.

11.          Wagoner, J.A. and V.S. Pande, Reducing the effect of Metropolization on mixing times in molecular dynamics simulations. The Journal of chemical physics, 2012. 137(21): p. 214105.

12.          Börjesson, U. and P.H. Hünenberger, Explicit-solvent molecular dynamics simulation at constant pH: methodology and application to small amines. The Journal of chemical physics, 2001. 114(22): p. 9706-9719.