Re: Question on restraint behavior when part of system decouples

From: Randy J. Zauhar (r.zauhar_at_usciences.edu)
Date: Mon Jan 29 2018 - 13:24:45 CST

Dear Jerome and Brian,

Thanks for those detailed and helpful responses!

Jerome, you wrote:

I'm not sure which recommended protocol you have in mind, but translational and rotational restraints when decoupling from isotropic bulk solution would have zero effect: they do not change the excess free energy that is being calculated. Restraint in the bulk calculation will not mimic the entropic effects of binding. More typically, restraints in the "protein" leg of the alchemical transformation are there to improve numerical convergence by increasing the overlap between the end-states. And then the free energy of those restraints needs to be calculated to account for those rotational and translational degrees of freedom.

The protocol I have in mind is the one shown in the NAMD free energy tutorial, “Protein:ligand standard binding free energies…”, especially Fig. 1 and associated narrative (Paragraphs numbered 1-5, copied below). The idea is that only the conformation of the ligand is ‘pinned’ in the bulk simulation, while in the complex, ligand position, orientation AND conformation are restrained. I would think that the three types of constraint ‘capture’ three important components of small molecule binding in terms of entropy loss, and don’t understand why only one of these contributions is ‘picked out’ in the bulk simulation. Put simply, why not just annihilate the ligand in solution with no constraints?

And by the way, the tutorial implies that the purpose of the restraints is to avoid the ‘wandering ligand’ problem, where as the ligand decouples it starts moving around and gets tangled up with the protein, leading to huge energies and poor accuracy, a very real problem! - but I don’t see where it discusses entropy?

Thanks!

Randy

 To quote from the tutorial:

(1) To circumvent the wandering-ligand problem arising when the substrate is decoupled from the protein, restrain the former in the conformation, orientation and position representative of the protein:ligand complex.

(2) Determine in the bound state (“site”) the alchemical free-energy change for decoupling reversibly from the protein the ligand restrained in its native conformation, orientation and position, i.e., protein : ligand∗ .

(3) Determine in the free, unbound state (“bulk”) the alchemical free-energy change for decoupling reversibly from the bulk water the ligand restrained in its native conformation, i.e., ligand∗ .

(4) Determine in the bound state (“site”) the free-energy contribution for maintaining the ligand in its native conformation, orientation and position by means of restraining potentials (uc, uΘ, uΦ, uΨ, uθ, uφ and ur). Towards this end, the set of collective variables utilized are a root mean-square deviation with respect to the conformation of the ligand in the bound state, the three Euler angles, Θ, Φ and Ψ, the two polar angles, θ and φ, and the Euclidian distance, r, separating the substrate from the protein.

(5) Determine in the free, unbound state (“bulk”) the free-energy contribution for maintaining the ligand in its native conformation by means of a restraining potential (uc). Towards this end, the chosen collective variable is the root mean-square deviation with respect to the conformation of the substrate when bound to the protein.

On 29Jan, 2018, at 12:19 PM, Brian Radak <brian.radak_at_gmail.com<mailto:brian.radak_at_gmail.com>> wrote:

Jerome has already given a great response, but I think I can add a little bit regarding alchemical bonds and the use of extraBonds. These features were added recently (by me) and are probably not well known, nonetheless, they were not really intended for ligand binding, so you probably won't need to change any of them.

As of NAMD 2.11 bonds in the PSF (or added via extraBonds) can be affected by alchemical coupling if desired, but not be default. You can set alchBondLambdaEnd in much the same way as alchVdwLambdaEnd to control the region on which this scaling occurs and how quickly. There is also the alchBondDecouple keyword, which is somewhat similar to alchDecouple.

There are two modes:

alchBondDecouple = False (default) - Only bonds involving atoms outside the alchemical region are scaled, bonds between atoms in the alchemical region are not scaled. If you are just decoupling a ligand, then this mode does nothing different.

alchBondDecouple = True - All bonds involving alchemical atoms are decoupled. This means that the alchemical region should blow apart into non-interacting ideal gas atoms. You probably don't want this.

There is also generally a lot of confusion surrounding the alchDecouple keyword, which decouples non-bonded interactions between alchemical atoms when off (default), but keeps them when true. With the default bonded settings, alchDecouple = True means that the non-interacting ligand endpoint behaves just like the molecule in the gas phase, but with a periodic boundary if using PME. The use of this keyword can totally change what thermodynamic cycle is appropriate and what assumptions are being made regarding the unbound ligand.

HTH,
Brian

P.S. WARNING - I just discovered a subtle bug when using LJcorrection and alchVdwLambdaEnd = 0.0 - DO NOT USE THIS COMBINATION WITH CONSTANT PRESSURE. I'll hopefully have a patch up soon.

On Mon, Jan 29, 2018 at 10:53 AM, Jérôme Hénin <jerome.henin_at_ibpc.fr<mailto:jerome.henin_at_ibpc.fr>> wrote:
Hi Randy,

On 29 January 2018 at 16:19, Randy J. Zauhar <r.zauhar_at_usciences.edu<mailto:r.zauhar_at_usciences.edu>> wrote:
Hi, I have always been troubled by the following aspect of the way constraints are applied in the typical ‘alchemical’ transformation when evaluating free energy of binding of ligand to protein.

1) If ligand is in complex with protein with restraints, and some restraints involve both ligand and protein atoms (say dihedral angle involving protein and ligand) and others only ligand atoms (say ligand RMSD), and the ligand is decoupled from the complex - do all restraint energy contributions ‘disappear’ during decoupling, or only the protein-ligand part, or do NONE of the contributions disappear? In the end, is the ligand effectively unrestrained in vacuum, or are the restraints still in place?

If the restraints are implemented within the Colvars module, then they are unaffected by the alchemical lambda and will persist in the decoupled system. If they are implemented as extraBonds, I am not entirely sure, depending on the options that allow for perturbing bonded terms.

2) Similarly, when the ligand is decoupled from bulk solvent, why does the recommended protocol only restraint the conformation, why not position and and orientation as well, to mimic the loss of positional and rotational freedom that occurs upon complex formation?

I'm not sure which recommended protocol you have in mind, but translational and rotational restraints when decoupling from isotropic bulk solution would have zero effect: they do not change the excess free energy that is being calculated. Restraint in the bulk calculation will not mimic the entropic effects of binding. More typically, restraints in the "protein" leg of the alchemical transformation are there to improve numerical convergence by increasing the overlap between the end-states. And then the free energy of those restraints needs to be calculated to account for those rotational and translational degrees of freedom.

And again, if ligand in bulk decouples from solvent, is the end result a small molecule in vacuum, or a small molecule in vacuum with restraints? Those are very different situations!

They are, but the effect of those restraints in vacuum is exactly the same as in bulk solution, so their impact on the deltaG of decoupling is zero.

I had assumed that the decoupling ONLY involved the nonbonded interactions between parts of the system,

Typically, yes.

and that you needed to handle everything involving restraints via integration;

Or analytically if possible.

but if that is so, my mental picture of how the thermodynamic cycle works is really off.

There is more than one way to design a thermodynamic cycle for double decoupling, but I think the statements I've made above are general.

Jerome

Randy J. Zauhar, PhD

Prof. of Biochemistry

Dept. of Chemistry & Biochemistry
University of the Sciences in Philadelphia
600 S. 43rd Street
Philadelphia, PA 19104

Phone: (215)596-8691
FAX: (215)596-8543
E-mail: r.zauhar_at_usciences.edu<mailto:r.zauhar_at_usciences.edu>

“Yeah the night is gonna fall, and the vultures will surround you /
And when you’re lookin’ in the mirror what you see is gon’ astound you"

  — Death Cab for Cutie, “Monday Morning"

This archive was generated by hypermail 2.1.6 : Thu Sep 19 2019 - 23:19:23 CDT