[gmx-users] Does gmx covar/gmx anaeig give <dS> or T<dS> for ligand binding?
Billy Williams-Noonan
billy.williams-noonan at monash.edu
Wed Jun 22 06:44:43 CEST 2016
I re-did the calculation. When considering the entire biomolecule of
each ensemble:
<S(P.L)>' = 51760 J/mol K
<S(P)> = 50640 J/mol K
<S(L)> = 814 J/mol K
Resulting in <dS>= -0.51 kJ/mol K
And when just considering Protein-H, I got:
<S(P.L)>' = 31941 J/mol K
<S(P)> = 31340 J/mol K
<S(L)> = 640 J/mol K
Resulting in <dS>= -0.69 kJ/mol K
These values make more sense given my enthalpy calculation with g_mmpbsa
is likely not converged. Thank you for your time and patience. :)
Billy
On 22 June 2016 at 13:40, Billy Williams-Noonan <
billy.williams-noonan at monash.edu> wrote:
> 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
>
>
--
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
More information about the gromacs.org_gmx-users
mailing list