[gmx-users] Does gmx covar/gmx anaeig give <dS> or T<dS> for ligand binding?
David van der Spoel
spoel at xray.bmc.uu.se
Thu Jun 23 14:16:17 CEST 2016
On 22/06/16 08:22, Billy Williams-Noonan wrote:
> Yes, still with two ligands... So I assume that <dS> should be halved to
> about -0.25 kJ/mol K, which would give T<dS>= -75.
>
> I found out recently that we have access to a node with 1TB of RAM.
>
> So the solvent is still relevant? If the amount of memory I need for the
> entropy calculation is given by *( 3*N )^2 * 4 *where N is the number of
> atoms, I should be able to calculate the entropy of systems (including
> solvent) on that node. If I have to multiply by the number of frames then
> I can't.
Number of frames should not be in the formula.
However including the solvent in the analysis will not help I think.
you need a PMF to get anything close to reliable.
Forget about MMPBSA.
>
> Is the formula above correct? Because according to it, with about 50,000
> atoms in my system, I would need 90 GB of RAM, which is odd, because GMX
> covar crashed when I used a node with 128GB RAM and one core.
>
> Billy
>
>
>
> On 22 June 2016 at 16:05, David van der Spoel <spoel at xray.bmc.uu.se> wrote:
>
>> On 22/06/16 06:44, Billy Williams-Noonan wrote:
>>
>>> 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
>>>
>> Still with two ligands? This still corresponds to a binding entropy change
>> of -150 kJ/mol. Of course you are ignoring the entropy change of the water
>> which is probably almost the same magnitude and with opposite sign. If you
>> want more quantitative results you could consider doing a PMF but your
>> ligand is very large so that will be difficult to converge as well.
>>
>>
>>
>>>
>>> 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
>>>>
>>>>
>>>>
>>>
>>>
>>
>> --
>> 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.
>>
>
>
>
--
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
More information about the gromacs.org_gmx-users
mailing list