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

Wed Jun 22 05:40:39 CEST 2016

Sorry that was the ATB, not the ATP

On 22 June 2016 at 13:39, Billy Williams-Noonan <
billy.williams-noonan at monash.edu> wrote:

> 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
>
>


-- 
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