Constant-pH simulation

[jchodera]

I wanted to restart the discussion about how we should best implement constant-pH methods into OpenMM, and most importantly, which layer this should be implemented into.

The best strategy is to combine the Mongan framework [1] for use in implicit solvent with the Stern framework [2] for explicit solvent within an NCMC move [3]. This effectively means that the proton will be switched on or off (with a set of atomic partial charges changing as this happens, but no LJ changes) over the course of a user-selected N steps, followed by a Metropolis-like acceptance/rejection step.

For the question of which layer(s) to implement this, there are several possibilities:

In the OpenMM API, where the user could, say, create a MonteCarloTitration force, with an interface like this:

force = simtk.openmm.MonteCarloTitration(temperature, pH, nsteps)
titration_group_index = force.addTitratableGroup(atom_indices)
force.addTitrationStates(titration_group_index, pKref, relative_energy, charges, proton_count)
force.addTitrationStates(titration_group_index, pKref2, relative_energy2, charges2, proton_count2)
...
system.addForce(force)

It is probably straightforward to implement these methods, along with some others to query the protonation states (setSwitchingSteps, getSwitchingSteps, getNumTitratableGroups, getNumTitrationStates, getTitrationState, getTitrationStateTotalCharge, setTitrationState, update [to perform an MC update of titration state], getAcceptanceProbability [to query the acceptance probability], getNumAttemptsPerUpdate, setNumAttemptsPerUpdate). Nearly all of these are implemented in the Python example code I emailed earlier.

This could be implemented at the C++ level in a Platform-independent manner, thus avoiding the need to code any GPU kernels.

At the OpenMM app layer, we could add information to the ffxml files about which residues are titratable forms of the same residue, providing reference pKas (pKref), proton counts, and potentially estimated relative energies. We would also want to implement an automated "calibration" routine that, given a set of parameters the system of interest (forcefield, nonbonded treatment), runs short alchemical free energy calculations on individual titratable residues present in the protein to compute the free energy difference between titratable states to use for the relative_energy of each state. The results from this could be cached to avoid the overhead of having to do this each time.

We could also provide routines to read an AMBER-generated .cpin file for users who wanted to depend on AmberTools for setup. My script gives an example of that.

References:

[1] Mongan J, Case DA, and McCammon JA. Constant pH molecular dynamics in generalized Born implicit solvent. J Comput Chem 25:2038, 2004.
http://dx.doi.org/10.1002/jcc.20139

[2] Stern HA. Molecular simulation with variable protonation states at constant pH. JCP 126:164112, 2007.
http://link.aip.org/link/doi/10.1063/1.2731781

[3] Nonequilibrium candidate Monte Carlo is an efficient tool for equilibrium simulation. PNAS 108:E1009, 2011.
http://dx.doi.org/10.1073/pnas.1106094108

[jchodera]

Note that this scheme can simply be extended to handle small-molecule titratable groups (with the multiple protonation states built externally and rolled into the ffxml file), as well as tautomers (same deal).

[jchodera]

Looping in @VijayPande @peastman @jwagoner

[jchodera]

The actual update scheme deserves more comment, since this may actually make it less suitable as a Force than something else.

During a protonation state update:

• One or more titratable groups is selected to be changed to some other titration state. There should likely be default probabilities that can be set by the user.
• The initial total energy (kinetic + potential) is recorded.
• For each step of nsteps, the charges are interpolated between old and new states, and a step of velocity Verlet (or another symplectic integrator) is taken.
• The final total energy is recorded.
• A Metropolis-like acceptance criterion is used to determine whether to accept the new protonation state, charges, configuration, and velocity; or to reject them all and return to the initial state. Note that there is some subtlety here in that the velocity must be inverted (multiplied by -1) on either success or rejection. The rule can be either (e.g. always invert velocity on rejection), but it must be uniformly applied.

This means that the update() call to MonteCarloTitration would actually need to run some simulation and change the coordinates, velocities, and charges, which may break how Force objects are currently conceptualized. To me, this still sounds like the most straightforward interface, however.

Instead, it's possible that an Integrator subclass (or CustomIntegrator) could be used for NCMC updates for explicit solvent, though since we cannot currently have more than one Integrator associated with a Context, this would require reinitializing the Context every time we wanted to try an MC protonation state change. That would not be ideal.

Other ideas?

[peastman]

I'll take a look at the papers. I agree that having the Force run its own integration inside updateContextState() sounds like a reasonable approach. That way it isn't dependent on what integrator you're using. And it can ensure that no other Force (like a thermostat or barostat) will modify the state in the middle, which would of course ruin the results.

How frequently would MC steps be attempted?

Is it ok for this to only work with NonbondedForce? Or would you want it to be able to handle custom nonbonded forces, AmoebaMultipoleForce, etc.?

[jchodera]

@peastman Thanks!

The Mongan and Stern papers are probably the best to start with. Our paper just provides a general framework to knit the two together, and the presentation is complex for generality, but the application to this special case is simple. We can help with that.

The frequency of MC steps should be either set by the user (e.g. every 500 steps for implicit; every 5000 steps for explicit). If this scheme is used, we need to also be able to specify a probability that a single titratable group is changed, or two groups simultaneously, or perhaps three or more (e.g. 50%, 30%, 20%).

A simpler alternative could be to set a probability per timestep that a protonation state move will be attempted for each titratable group, potentially something that the user could set. This way, there will naturally be some 2- or 3-group moves attempted at the same time.

I think we eventually want it to work with CustomNonbondedForce and AmoebaMultipoleForce as well so that we can use this feature for binding free energy calculations. That makes things considerably more complex, I realize, so I'll think a bit harder as to whether we can at least avoid needing this feature for CustomNonbondedForce (since there isn't even a concept of parameters for charges there). The one way this could work is for it to apply to parameters with a specific name in Custom*Force terms, such as "charge".

[jchodera]

One further complication: scaled 1-4 interactions involving charges would need to change when the titratable group charges are modified.

[jchodera]

Another complication: For explicit solvent constant-pH, charge creation/deletion will turn on the PME background charge, which could cause complications that necessitate corrections of the sort David Mobley is working out.

We may want to instead neutralize/create a counterion when a titratable group has its charge changed. Jason Wagoner suggests this is most easily done by turning a water molecule into an ion (or vice versa), which can potentially be done by simply linearly interpolating the water/ion Lennard-Jones parameters in addition to the charges. That suggests we would want to allow the specification of both charges and Lennard-Jones parameters (sigma, epsilon) for the different titration states.

We are experimenting with this in Python at the moment and will make sure it actually works well before we commit to the API, but it would still be useful to discuss what the API would look like over the next week or two.

[VijayPande]

I like that idea. I wonder what ion makes sense. using OH- and OH3+ in a water molecule near the proton which we change in the protein would be cute but maybe opens us up for unexpected problems.

Sent from my Phone. Sorry for the brevity or unusual tone.

On Oct 23, 2013, at 3:01 PM, John Chodera [email protected] wrote:

Another complication: For explicit solvent constant-pH, charge creation/deletion will turn on the PME background charge, which could cause complications that necessitate corrections of the sort David Mobley is working out.

We may want to instead neutralize/create a counterion when a titratable group has its charge changed. Jason Wagoner suggests this is most easily done by turning a water molecule into an ion (or vice versa), which can potentially be done by simply linearly interpolating the water/ion Lennard-Jones parameters in addition to the charges. That suggests we would want to allow the specification of both charges and Lennard-Jones parameters (sigma, epsilon) for the different titration states.

We are experimenting with this in Python at the moment and will make sure it actually works well before we commit to the API, but it would still be useful to discuss what the API would look like over the next week or two.

—
Reply to this email directly or view it on GitHub.

[VijayPande]

PS One obvious problem with what I just suggested is that once the water is changed, it could diffuse away and we'll start having a mix of water + OH- + OH3+ which would be weird (and would require us to handle protonation of the water dynamically too). I guess while we're at it, we could do that, but that just seems like opening up a huge can of worms.

Sent from my Phone. Sorry for the brevity or unusual tone.

On Oct 23, 2013, at 3:59 PM, Vijay Pande [email protected] wrote:

I like that idea. I wonder what ion makes sense. using OH- and OH3+ in a water molecule near the proton which we change in the protein would be cute but maybe opens us up for unexpected problems.

Sent from my Phone. Sorry for the brevity or unusual tone.

On Oct 23, 2013, at 3:01 PM, John Chodera [email protected] wrote:

Another complication: For explicit solvent constant-pH, charge creation/deletion will turn on the PME background charge, which could cause complications that necessitate corrections of the sort David Mobley is working out.

We may want to instead neutralize/create a counterion when a titratable group has its charge changed. Jason Wagoner suggests this is most easily done by turning a water molecule into an ion (or vice versa), which can potentially be done by simply linearly interpolating the water/ion Lennard-Jones parameters in addition to the charges. That suggests we would want to allow the specification of both charges and Lennard-Jones parameters (sigma, epsilon) for the different titration states.

We are experimenting with this in Python at the moment and will make sure it actually works well before we commit to the API, but it would still be useful to discuss what the API would look like over the next week or two.

—
Reply to this email directly or view it on GitHub.

[jchodera]

@VijayPande: There's obviously lots of room for experimentation and research here, which is why it's good to discuss an API that will meet all our research and production needs for constant pH and potentially also tautomers, GCMC ion concentrations, and even GCMC ligand design.

The simplest approach is to use monovalent counterions like Na+ and Cl- (or whatever is around). The reference calculations (terminally-blocked amino acids) would also use the same ions, so the exact ion identity would mostly cancel out (though some additional corrections from David Mobley's paper may still be required).

Jason Wagoner and I can look into the model systems for this using my Python model implementations over the next couple of weeks as we work to agree on an API for OpenMM.

[VijayPande]

I really like that approach. With a sufficiently general API, people can start to play around and it will be interesting to see what comes up.

Thanks,

Vijay

Sent from my Phone. Sorry for the brevity or unusual tone.

On Oct 23, 2013, at 4:05 PM, John Chodera [email protected] wrote:

@VijayPande: There's obviously lots of room for experimentation and research here, which is why it's good to discuss an API that will meet all our research and production needs for constant pH and potentially also tautomers, GCMC ion concentrations, and even GCMC ligand design.

The simplest approach is to use monovalent counterions like Na+ and Cl- (or whatever is around). The reference calculations (terminally-blocked amino acids) would also use the same ions, so the exact ion identity would mostly cancel out (though some additional corrections from David Mobley's paper may still be required).

Jason Wagoner and I can look into the model systems for this using my Python model implementations over the next couple of weeks as we work to agree on an API for OpenMM.

—
Reply to this email directly or view it on GitHub.

[swails]

We may want to instead neutralize/create a counterion when a titratable group has its charge changed. Jason Wagoner suggests this is most easily done by turning a water molecule into an ion (or vice versa), which can potentially be done by simply linearly interpolating the water/ion Lennard-Jones parameters in addition to the charges.

Perhaps this is jumping the gun a little bit? What leads you to believe that maintaining a net neutral cell is required? I've seen plenty of evidence that for systems of typical size, having a net, non-zero charge does not cause actual harm, and the one paper I've seen that tried to assert the importance of charge neutrality did not even come close to convincing me.

I agree it would be nice to have the flexibility to scale the Lennard-Jones parameters, but I don't think that should be an initial concern, IMO...

[jchodera]

What leads you to believe that maintaining a net neutral cell is required?

David Mobley has a lengthy paper coming out (with Philippe Hünenberger) that gives a thorough accounting of all of the details for charged ligand binding. We anticipate a number of these corrections / issues will also apply to constant pH in explicit solvent with PME.

[swails]

David Mobley has a lengthy paper coming out (with Philippe Hünenberger) that gives a thorough accounting of all of the details for charged ligand binding.

I'd be very interested to read this account when it comes out (I assume it's not out yet).

We anticipate a number of these corrections / issues will also apply to constant pH in explicit solvent with PME.

In the method I developed and implemented, maintaining cell neutrality was unimportant. Admittedly it was a hybrid method and the protonation state moves were evaluated in implicit solvent, but in the tests I ran I saw no discernible effect on the dynamics compared to what I was seeing with relevant constant protonation runs in a neutral unit cell.

Even if the protonation state move is evaluated in pure explicit solvent a la Stern et al., I expect that the use of the model compound will mitigate the effect of the fluctuating charge on the computed pKa shift. While I have data suggesting that using counterions is important even in hybrid GB/PME method (where they are ignored for the protonation state move), they didn't have to morph to maintain charge neutrality.

The one paper that suggested charge neutralization is important (http://dx.doi.org/10.1063/1.4766352) is misleading. They were using small boxes with tiny compounds with GRF (not PME), and only 1 other counterion. The titration process led to massive fluctuations in the local ionic strength around the solute as the protonation state(s) changed. Had the authors used a larger box with more (fixed) ions, I'm quite confident they would have seen zero effect from their charge-leveling procedure. Of course they were also not using a long-range method formally dependent on a net neutral cell so your proposed simulations may suffer in a way that the existing continuous protonation state methods do not.

So far I know of 3 groups doing continuous CpHMD in explicit solvent based on lambda-dynamics. The Shen group claims charge leveling is important, but they use GRF (where charge neutrality is not nearly as important as Ewald-based methods) and seem to be observing an obvious artifact unrelated to the issue at hand. The Brooks group claims charge leveling is not important and they seem to get good results with proteins and nucleic acids, but they neglect long-range electrostatics altogether and use a force-shifting approach instead. The Grubmuller group, on the other hand, actually uses PME and no charge neutralization, although they haven't done anything larger than a small test system that I'm aware of.

I've heard hearsay of pathological cases (particularly involving bilayers and micelles) where a net-charge on a unit cell can lead to unphysical phenomena, but in my experience with constant pH MD methods (my entire grad school career) I would suggest that charge neutralization is not an urgent concern.

[jchodera]

I'd be very interested to read this account when it comes out (I assume it's not out yet).

Feel free to email him! It's been accepted, but he's only distributing preprints to people that ask for a copy.

[jchodera]

Even if the protonation state move is evaluated in pure explicit solvent a la Stern et al., I expect that the use of the model compound will mitigate the effect of the fluctuating charge on the computed pKa shift.

My expectations is that most of the effects will indeed be accounted for in the model compound simulations as well, though in discussions with David, it appears that not all of the effects will be accounted for. We will need to establish which effects remain and whether their magnitude is large enough to warrant inclusion.

[jchodera]

So far I know of 3 groups doing continuous CpHMD in explicit solvent based on lambda-dynamics.

I (and others) are suspicious of continuous-lambda protonation states, since the expectation of some property of a continuous protonation state is not guaranteed to be equivalent to the expectation of that property evaluated over a mixture of physical protonation states.

I've heard hearsay of pathological cases (particularly involving bilayers and micelles) where a net-charge on a unit cell can lead to unphysical phenomena, but in my experience with constant pH MD methods (my entire grad school career) I would suggest that charge neutralization is not an urgent concern.

I'm hoping that we can integrate all input (especially yours!) to ensure that the API changes support several styles of constant-pH dynamics. My thought was that we should be able to support keeping changes charge-neutral by creating/destroying counterions, but we would by no means require it.

[jwagoner]

I share John's concerns of lambda-dynamics for constant pH, and will also lay my vote in favor of versatility. Our interest in incorporating insertion and deletion movesets for ions extends well beyond this discussion, so I don't consider it to be more or extraneous work. It will be great if the background charge causes no issue--the resulting movesets will be simpler and have higher acceptance rates. But, it's worth discussing and exploring, starting with David Mobley's results on the matter.

[VijayPande]

+1

Thanks,

Vijay

Sent from my phone. Sorry for any brevity or unusual tone.

On Oct 23, 2013, at 11:41 PM, jwagoner [email protected] wrote:

I share John's concerns of lambda-dynamics for constant pH, and will also lay my vote in favor of versatility. Our interest in incorporating insertion and deletion movesets for ions extends well beyond this discussion, so I don't consider it to be more or extraneous work. It will be great if the background charge causes no issue--the resulting movesets will be simpler and have higher acceptance rates. But, it's worth discussing and exploring, starting with David Mobley's results on the matter.

—
Reply to this email directly or view it on GitHub.

[peastman]

When evaluating the PME reciprocal space terms, we ignore the zero frequency component. This is formally equivalent to filling all of space with a uniform charge density to neutralize the system. So strictly speaking, your system is always neutral when using PME.

[jchodera]

Excellent point, @peastman. I think @swails and I are discussing whether or not it is important to be sure to add/delete a counterion to ensure that the "uniform neutralizing charge density" remains zero, since the creation of this uniform neutralizing charge may require a number of corrections to be added to the Monte Carlo acceptance criteria that could make the scheme much more complex. These corrections are required when computing alchemical free energies involving the creation/deletion of charged ligands.

[peastman]

How is that case any different from adding an ion?

[swails]

How is that case any different from adding an ion?

An ion does not spread uniformly over the PME grid; its contribution to the charge grid is localized around its location. This may be an important distinction.

EDIT: Apologies, forgot the first not in my original post.

[jchodera]

As @swails points out, a counterion is a physical localized species, while the uniform neutralizing charge is spread over the entire box. These could lead to different effects on the configurations sampled at equilibrium.

[peastman]

Given all the potential complications here (many different Force classes to work with, both Coulomb and LJ parameters, 1-4 interactions, etc.), I wonder if it would be best to drive this from the Python level? That would allow a completely generic interface that can set arbitrary parameters on arbitrary objects. The only part that is really performance critical is the integration, so that's the only part that needs to be in lower level code.

[jchodera]

@VijayPande : Would you be comfortable with this feature being only implemented at the Python level?

Alternatively, we can do some profiling and tuning with the Python version and see if it still makes sense to implement the whole feature (or more of it) at the C++ level.

@peastman : I think we would still need some way to support steps of velocity Verlet integration at the C++ level so that we don't need to create a new Integrator and Context each time we want to try a protonation state change, because this will result in very poor efficiency. We can still prototype and debug the method in parallel without C++ support, but I'd be curious to hear your thoughts on how to support this integration at the C++ level.

To review, each attempted protonation state change includes these steps:

• Retrieve the coordinates
• Compute the initial potential energy
• Take N steps of velocity Verlet, where the steps of dynamics are interleaved with updates to some charges (and potentially LJ radii)
• Compute the final potential energy
• Accept or reject the new configuration and parameters, or return to the old ones

If you don't think Context creation would be enormously expensive in addition to these operations, we may not need C++ support...

[VijayPande]

I'd strongly prefer this at the C++ level for several reasons, most of which is that many of our key new apps (FAH Core17 and the new MSMBuilder) are C++ only apps.

I would be fine with developing first in python and they porting something to C++. I could see this as a powerful scheme to play around with several possibilities in Python and then port then one we like best to C++.

Sound good?

Sent from my Phone. Sorry for the brevity or unusual tone.

On Oct 31, 2013, at 12:52 PM, John Chodera [email protected] wrote:

@VijayPande : Would you be comfortable with this feature being only implemented at the Python level?

Alternatively, we can do some profiling and tuning with the Python version and see if it still makes sense to implement the whole feature (or more of it) at the C++ level.

@peastman : I think we would still need some way to support steps of velocity Verlet integration at the C++ level so that we don't need to create a new Integrator and Context each time we want to try a protonation state change, because this will result in very poor efficiency. We can still prototype and debug the method in parallel without C++ support, but I'd be curious to hear your thoughts on how to support this integration at the C++ level.

—
Reply to this email directly or view it on GitHub.

[jchodera]

What if we

1. finish our prototype in Python, tune it, and figure out what features are essential to move to C++
2. Implement essential components at the C++ layer and app layer. I could imagine this might include support for only the NonbondedForce and protein protonation states, with ligand promets and tautomers provided by a plugin that makes of Schroedinger or OpenEye tools for predicting these, generating ffxml files. You couldn't do alchemical calculations that involved protonation or tautomeric states because CustomNonbondedForce isn't supported (unless we use a trick, which I can elaborate on later) so that will be in Python in Yank. We can make sure whatever C++ support is implemented also enables more elaborate Python versions.

Would that work?

[VijayPande]

that works for me

Sent from my Phone. Sorry for the brevity or unusual tone.

On Oct 31, 2013, at 1:44 PM, John Chodera [email protected] wrote:

What if we

1. finish our prototype in Python, tune it, and figure out what features are essential to move to C++
2. Implement essential components at the C++ layer and app layer. I could imagine this might include support for only the NonbondedForce and protein protonation states, with ligand promets and tautomers provided by a plugin that makes of Schroedinger or OpenEye tools for predicting these, generating ffxml files. You couldn't do alchemical calculations that involved protonation or tautomeric states because CustomNonbondedForce isn't supported (unless we use a trick, which I can elaborate on later) so that will be in Python in Yank. We can make sure whatever C++ support is implemented also enables more elaborate Python versions.

Would that work?
—
Reply to this email directly or view it on GitHub.

[jchodera]

I've created a repository for experimenting with constant-pH methods in a pure Python implementation:

https://github.com/choderalab/openmm-constph

Currently, it just contains the result of my previous experiments, but I'll update this over the next few days to be a more complete testbed for experimenting with functionality we should pull into OpenMM.

[jchodera]

@VijayPande : I earlier alluded to a "trick" where we could avoid the need to do constant-pH simulations for alchemically modified states. This should be workable if we impose some sort of additional field that has the effect of confining the system to a single fixed protonation state. One way to do this is to simply fix an extreme protonation state of +$\infty$ or -$\infty$ and compute the free energy for changing this pH for both protein+ligand and protein // ligand physical endpoints. The free energy calculation would then proceed with a fixed extreme protonation state. In this way, there is no requirement for the ability to combine constant-pH and alchemically-modified Custom*Force simulations.

This may have some issues: The alchemical path may be inefficient for something with high charge, and the additional neutralizing counterions or high charge could cause convergence issues. Instead, we could create a nonphysical ensemble where there is not a single global pH, but each titratable group has its own pH that we individually extremize to keep the system in a fixed (but reasonable) protonation state.

[VijayPande]

That's an interesting idea. I am also interested in const pH more generally, especially as we move into systems where dynamic protonation is starting to be an issue.

Sent from my Phone. Sorry for the brevity or unusual tone.

On Nov 7, 2013, at 1:39 PM, John Chodera [email protected] wrote:

@VijayPande : I earlier alluded to a "trick" where we could avoid the need to do constant-pH simulations for alchemically modified states. This should be workable if we impose some sort of additional field that has the effect of confining the system to a single fixed protonation state. One way to do this is to simply fix an extreme protonation state of +$\infty$ or -$\infty$ and compute the free energy for changing this pH for both protein+ligand and protein // ligand physical endpoints. The free energy calculation would then proceed with a fixed extreme protonation state. In this way, there is no requirement for the ability to combine constant-pH and alchemically-modified Custom*Force simulations.

This may have some issues: The alchemical path may be inefficient for something with high charge, and the additional neutralizing counterions or high charge could cause convergence issues. Instead, we could create a nonphysical ensemble where there is not a single global pH, but each titratable group has its own pH that we individually extremize to keep the system in a fixed (but reasonable) protonation state.

—
Reply to this email directly or view it on GitHub.

[jchodera]

This is under active development here:
http://github.com/choderalab/openmm-constph

[peastman]

Since this is being implemented elsewhere, should we close this issue?

[jchodera]

Still working on an update on this with concrete proposals for OpenMM extensions to help support this; stay tuned. Please keep this open for now?

[peastman]

Just wondering what's the status of this? I notice the other repo you linked to has had no commits in well over a year.

[jchodera]

@pgrinaway and I have some modifications we haven't merged yet---I'm still hoping to have some time to do that immediately after sorting out the core 21 stuff needed to get any of our projects running.

My plan was to first get us all running in Python, but maybe it's worth restarting the discussion for the C++ API in case this was a desired feature for OpenMM 7 (which would be awesome).

[jchodera]

@bas-rustenburg is now picking up the constant-pH work.

[jchodera]

Maybe mark this as "enhancement" since we are working on it?

[peastman]

Done.