[gmx-users] PMF Transmembrane proteins

[Nash, Anthony]

Hi Justin,

Thanks for the help. The block analysis has done wonders! I had been interested in seeing how the PMF graphs converged, i.e., 0-1ns, 0-2ns, 0-3ns. Of course this always includes the earliest time steps. Long story short, by cutting off the first 10ns so I am only sampling 10ns-40ns, the free energy curves are very different.

Regions where the peptides are close together are not that different. But between 3nm and 7.5nm the curves plateau with a lot less energy difference between the like-with-like system. I will run bootstrap analysis, get out the error, see how that goes. 

Oddly enough the curve is not as smooth now, even though there is ample histogram coverage. I may push the simulations for longer.

Thanks, I think this is now moving forwards again.

Anthony 



________________________________________
From: gmx-users-bounces at gromacs.org [gmx-users-bounces at gromacs.org] on behalf of Justin Lemkul [jalemkul at vt.edu]
Sent: 30 December 2012 20:08
To: Discussion list for GROMACS users
Subject: Re: [gmx-users] PMF Transmembrane proteins

On 12/30/12 2:57 PM, Nash, Anthony wrote:
> Dear Justin
>
> Thanks for your reply.
>
> I am studying a TM peptide, looking at how favourable the front and reverse
> face are. Hence, different orientations but the same composition of amino
> acids.
>
> I then have a duplicate of this system (call it system B), but with a couple
> of residue substitutions. I am comparing free energies of System A front
> orientation with System A back orientation (and then B with B).
>
> I have separated the peptides, and ran umbrella simulations every 1 angstrom
> apart (I experimented with windows closer together, in some regions hence the
> generalisation to half an angstrom). Papers I have looked through which
> sample up to 2 nm apart have been running umbrella windows every 0.2 of an
> angstrom. My PMF plateau is around 6.5-7.5 nm, which implies around 78
> windows per system.
>

Without seeing the actual histograms and .mdp settings, I can't comment on this
approach directly.  It seems like you should be able to use far fewer windows,
perhaps with a slightly less restrictive force constant to allow for slightly
broader distributions.  That's a hand-waving guess, though, since I've seen none
of the data.

> Taking one system as an example, if I normalise at 7.5 nm the difference at
> their global minima is around 37.63 kJ mol−1 (8.99 kCal mol−1), but then when
> I break down the PMF graph (0-1ns, 0-2ns, 0-3nm) the regions greater than 4.5
> nm have not converged. Hence, normalising anywhere in this unconverged region
> would be wildly inaccurate. If I normalise to zero at the point in which one
> of the system pairs (A&A, or B&B) converges, the free energy different is 5
> kJ mol−1.
>

I think it would be more productive to do your block analysis such that you
discard initial frames to see what effect the (presumably) unequilibrated part
of the trajectory is.  Even with prior equilibration (i.e. with position
restraints), you still need to leave some frames behind as equilibration.
Analyzing 0-40 ns, 10-40 ns, 20-40 ns, etc will probably be more informative.
If you're always considering all frames in the WHAM analysis, that can throw off
the results considerably.

> If by error estimate you mean from the inbuilt bootstrap analysis, then I am
> very confident that the error along the entire reaction coordinate is very
> small; which is confusing as I am not sure whether the slower converging
> regions of the reaction coordinate (4.5 - 7.5) should shower greater error!?
>

That's not uncommon.  At greater distance, you have two peptides floating around
rather independently, so the restraint potential may vary a bit more.  At closer
COM distances, the interactions between the peptides keep the orientations a bit
more stable.

> Yes it is heart breaking to hear that I may need to run the system for
> longer. This however, isn't a massive problem, given the time of the year the
> HPCs I am using are a bit empty. Would you suggest I need to run the 40 ns
> for longer (continuation run), or begin brand new windows from scratch and
> include those in the g_wham calculations? I guess the 40 ns one, given that
> lipid reorientation needs to be taken into account.
>

There is no need to scrap the existing data.  You may only need to extend the
simulations, but do carry out the analysis I suggested before to see if there is
undue influence from the initial frames.

> I don't have numbers in front of me, but the system box size is big (I went
> through your tutorial several times, and I have seen the PMF graph as a
> result of a box too small). The longest i.e., the reaction coordinate, is
> along the Y axes, Z is normal to the bilayer. When I calculate the PMF, I
> have pull_dim = Y Y N.
>

Good to know.  Lipid ordering effects can extend up to a few nm though, so
that's why I asked.

-Justin

--
========================================

Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================
--
gmx-users mailing list    gmx-users at gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* 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/Support/Mailing_Lists