[gmx-users] PMF Transmembrane proteins

[Justin Lemkul]

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

========================================