[gmx-users] Does gmx covar/gmx anaeig give <dS> or T<dS> for ligand binding?

[Billy Williams-Noonan]

Hi David,

  Thanks again for responding... Sorry if I came across the wrong way.  I'm
not trying to disprove the code, but simply understand why my values don't
make sense  I trust your knowledge on this subject too, since I suspect
you're one of the geniuses who helped to develop g_covar/g_anaeig. :)  I
know the units for entropy too.

  I should explain that I have previously performed relative FEP
calculations of ligands binding to the site of interest, and reproduced
experimental binding affinities within 1.4 kcal/mol of experiment.  Ligand
topologies came from the ATP using a GROMOS united atom force-field.  So I
know that the protocol I use for system equilibration is working.

   Using the same equilibration protocol as with the FEP protocol, and
having tried an absolute FEP calculation with restraints that failed
dismally, I have a cyclic peptide that has mM affinity for the same binding
site as the aforementioned ligands (see above paragraph).  So, using the
same protein as a model and placing the cyclic peptide in the correct
orientation as determined by the crystal structure, I am trying to use
g_mmpbsa to get an absolute binding affinity.  Of course the entropic term
from the g_mmpbsa calculation is missing, so I am using g_covar and
g_anaeig to determine the entropy.

   You're right about the size of my ligand too of course.  The cyclic
peptide is 54 atoms in size and moves quite a lot in solution.  I am used a
Parrinello-Rahman/V-rescale NPT ensemble, set to 300K and 1 bar, for the
entirety of the 100ns simulation.  And my protein is a symmetrical dimer
(two of the same protein bound to each other) so there is one ligand for
each monomer, forming 108 atoms between the two ligands.

   When I initially made this thread, the variables I was talking about
were:

<S(P.L)> = 128,886 J/mol/K
    = Entropy of one ligand bound to one side of the protein dimer, despite
another ligand being bound on the other side.  a_1-3071 was selected
(twice) in g covar to represent the P.L complex, despite there being 3125
atoms in total

<S(P)> = 153,548 J/mol/K
    = Entropy of the protein

<S(L)> = 4137 J/mol/K
    = Entropy of the cyclic peptide

   So I redefined <S(P.L)>, as <S(P.L)>', and selected a_1-3125 this
morning, to get the entropy of the dimer complexed with two cyclic
peptides, and got a value of 51,759.8 J/mol/K.  I substituted this into the
equation (1)

<dS> = <S(P.L)>' - <S(P)> - 2*<S(L)>  --------------(1)

   I multiplied the entropy of the ligand by two to account for the fact
that the beginning state now has the two ligands and the protein in
solution, while the end state has the protein dimer complexed with those
two ligands. And, the answer was -110.072 kJ/mol/K.

   So I am clearly doing something wrong and I'd like some advice on what
it is...  I doubt it's an equilibration problem, since my FEP calculations
previously worked with the same equilibration protocol.  And I doubt this
is a convergence issue too, since a <dS> this high should prohibit binding
in most cases, and my ligand definitely binds as seen by viewing the
complex simulation on VMD.

   Advice? Thoughts?  I am about to try it with just the Protein-H atoms
from the index files to see if that changes anything...

Billy

On 21 June 2016 at 21:51, David van der Spoel <spoel at xray.bmc.uu.se> wrote:

> On 21/06/16 11:26, Billy Williams-Noonan wrote:
>
>> Hi Gromacs Users,
>>
>>   I have used gmx covar and gmx anaeig to generate three ensemble average
>> entropies over 100ns: first for a ligand in solution (<S(L)>), second for
>> a
>> protein in solution (<S(P)>) and third for their respective complex in
>> solution (<S(P.L)>).
>>
>>    My understanding is that the change in entropy upon binding is given
>> by:
>>
>> <dS> = <S(P.L)> - <S(P)> - <S(L)>   -----------(1)
>>
>>    Using gmx covar/gmx anaeig I got Quasi-Harmonic entropy estimates of:
>>
>> <S(P.L)> = 128,886 J/mol/K
>>
>> <S(P)> = 153,548 J/mol/K
>>
>> <S(L)> = 4137 J/mol/K
>>
>>    As stated, these values were generated using gmx covar/anaeig by
>> selecting for the relevant biomolecule in each ensemble and ignoring the
>> effect of solvent movement.
>>
> The unit printed by the program is J/mol K, which is the normal unit for
> entropy in all handbooks. You can not prove or disprove the correctness of
> the code by an example, you will have to check the code yourself if you
> doubt it.
>
> Looking at your numbers, they are huge and the difference is huge too.
> You should probably make sure first that all you simulations are in
> equilibrium. A ligand entropy och 4137 I would expect for an organic
> molecules with close to 100 carbon atoms.
>
>
>
>>    By subbing the above-described values into (1), I got about -28
>> kJ/mol/K
>> for <dS>, which is the right answer if the units are actually kJ/mol, and
>> not kJ/mol/K.  Strangely, upon multiplying by T, I got a value of -8640
>> kJ/mol, which is quite obviously wrong.
>>
>>    So does (1) yield a value for <S> or T<dS> ?  Is anyone able to explain
>> this to me?
>>
>>    Kind regards,
>>
>> Billy
>>
>>
>>
>
> --
> David van der Spoel, Ph.D., Professor of Biology
> Dept. of Cell & Molec. Biol., Uppsala University.
> Box 596, 75124 Uppsala, Sweden. Phone:  +46184714205.
> spoel at xray.bmc.uu.se    http://folding.bmc.uu.se
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>



-- 
Billy Noonan*    |    *PhD Student    *|*    Bsci ( *Adv* ), IA Hon

*LinkedIn Profile
<http://www.linkedin.com/profile/preview?locale=en_US&trk=prof-0-sb-preview-primary-button>
**|*   +61420 382 557

Monash Institute for Pharmaceutical Sciences ( *MIPS* )
Royal Parade, Parkville, 3052