[gmx-users] Re: Re: (dH/dl) calculation

[David Mobley]

Mauricio,

> Exactly: I`d want to calculate "dG/dl for some particular component of
> the energy". Why? I did 2 FEPs, so I have two Delta_Gs for a mutations in
> 2 states: folded and unfolded. Both Delta_Gs are positive values. The
> difference between them (Delta_Delta_G) is the value I was looking for.
> This value agrees with experimental value of Delta_G of unfolding between
> 2 mutant proteins. Now, I want to know what is stabilizing or
> destabilizing each state with respect the other. Analysis of each
> component (e.g. Coul-SR:Prot-Prot) is impossible: it happens that I'm
> searching a delta_E of < 10 Kcal between a magnitude of 10E6. So I
> thought that if I could obtain a <dE/dl> for a component and then
> integrate them along lambda (TI method), I could get some idea about its
> magnitude.

OK, I have two comments about this, and then I'll give you an
alternative suggestion, which is somewhat better, but probably still
not great. First, the comments:
(1) To obtain component-wise dG/dl values, you would need to modify
the source code (probably fairly extensively) to get GROMACS to save
components of the dG/dl, not just the total. You are of course free to
do this, but I can't help.
(2) This would probably not be particularly informative, as you'll
probably find that most of the difference is in Lennard-Jones and
Coulomb interactions -- in the unfolded state, between the protein and
the solvent, and in the folded state, between the protein and itself
(and possibly also solvent). Maybe you learn something from that, but
I doubt you'll learn anything terribly surprising, given the amount of
work it will take (from (1)).

Something you could do that would be somewhat less work than the
above, and probably about as informative, is to break the
transformation into several parts, where you (for example) first
change the charges involved in your mutation, and then change the
Lennard-Jones interactions. Then you'll get several component free
energies -- a free energy of changing the electrostatics, and a free
energy of changing the LJ interactions. These are, of course,
path-dependent, but still, comparing the charging component, for
example, between the unfolded and folded states may give you some idea
of whether the difference is largely electrostatic or not.

Otherwise, you may want to look at structural order parameters to
figure out what is going on -- i.e. the radial distribution function
of water around the residue of interest, number of hydrogen bonds,
things like that.

One other comment on the approach of looking at components of the
potential energy: Keep in mind, the force field was parameterized to
get the total potential energy right (as a function of the
coordinates), *not* to give correct individual components. Personally
(although others may disagree) if I were reviewing a paper that tried
to make a significant case by looking at the breakdown of a free
energy by components in the way you are describing, I would want to
see some validation on a simple system where the right answer is known
in order to show that this approach really works. In other words, I
wouldn't trust that the results are right just because you can do it.

Best wishes,
David


> Mauricio
>
> > Mauricio,
> >
> > > That was exactly what I wanted to know. From one hand, because I just
> > > wanted to know it and, from the other hand, because I think that if I
> > have
> > > some idea about how this calculus is carried out, I was able to
> > estimate a
> > > Delta_E for some energy components, for example,
> > Coul-SR:Protein-Protein
> > > by means of calculating the corresponding <dE/dl> (even when I know
> > that
> > > the value is not a state function as stated in the Mark & v Gunsteren
> > > paper). Maybe I should have started from this point. Is it posible?
> >
> > I'm sorry, I have no idea what you're asking here. Are you talking
> > about trying to calculate a dG/dl for some particular component of the
> > energy? That doesn't sound useful, to me, since the components
> > generally are all interdependent. And if you are trying to separate
> > out the contribution of different components to the total dG/dl, this
> > would require source code modifications, in general. Again, if you can
> > be more specific about what exactly you want to do (what problem are
> > you trying to solve?) people may be able to be more helpful.
> >
> > David
> >
> >
> > > > Message: 6
> > > > Date: Mon, 6 Nov 2006 10:19:05 -0800
> > >
> > > > Mauricio,
> > > >
> > > > I'm somewhat confused by your question and notation. However, I
> > think
> > > > the basic answer is something like this: In molecular dynamics, you
> > > > know the Hamiltonian from which you are sampling; call it H(x,p,
> > l),
> > > > where x denotes all of the positions, p the momentums, and l
> > lambda.
> > > > This, of course, is closely linked to the potential energy. Anyway,
> > at
> > > > any snapshot, you can simply take the derivative dH(x,p,l)/dl, and
> > you
> > > > have dH/dl at that snapshot. This is usually straightforward since
> > you
> > > > know the dependence of all of the terms in your Hamiltonian on
> > lambda,
> > > > so you actually have the functional form for dH/dl as well -- so it
> > > > just involves taking the appropriate combination of positions,
> > > > momentums, etc. This is of course all handled internally by the
> > code.
> > > > <dH/dl>, then, is just the time-average of dH/dl, which can be
> > > > evaluated every step by the code.
> > > >
> > > > I am not sure if that's helpful at all, as I'm not entirely sure
> > what
> > > > problem you're having. After all, whenever you do TI calculations
> > in
> > > > GROMACS, the code gives you back dG/dl (or dH/dl, or dA/dl) for
> > every
> > > > snapshot in an xvg output file. Are you just confused about how the
> > > > code gets this (I think I just answered that above), or are you
> > trying
> > > > to figure out how to use it? If you're confused about how to use
> > it,
> > > > try to ask a question that relates to the specific issue you're
> > > > confused about.
> > > >
> > > > Best wishes,
> > > > David Mobley
> > > > UCSF
> > > >
> > > >
> > > > On 11/5/06, Mauricio Sica <msica at unq.edu.ar> wrote:
> > > > > Dear experts
> > > > >
> > > > > I am doing FEP (thermodynamic integration method) simulations.
> > > > > I have a questions about <dH/dl> calculation in GROMACS.
> > > > > Take in mind equation 3.77 from the GROMACS 3.3 manual.
> > > > > There, dA/dl is calculated as
> > > > >
> > > > > dA/dl = SS{ (dH/dl) exp()dp dq } / SS{ exp()dp dq = <dA/dl>NVT;l
> > }
> > > > >
> > > > > where SS are doble integrals (sorry for the notation).
> > > > >
> > > > > My question is: how is (dH/dl) (in the middel-term of the
> > equation)
> > > > > calculated?
> > > > > My idea is that the difference V(L=1)-V(L=0) is calculated for
> > every
> > > > time
> > > > > step (irrespective of the lambda value of the simulation) and
> > <dG/dl>
> > > > is
> > > > > the time average of that difference.
> > > > >
> > > > > <dG/dl> = < V(L=1)(i)-V(L=0)(i)/1 >
> > > > >
> > > > > Is this correct?
> > > > >
> > > > >
> > > > > Thanks
> > > > >
> > > > >
> > >
> > >
> > >
> >
>
>
>
> _______________________________________________
> gmx-users mailing list    gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>