[gmx-users] Where to stop with TM protein PMF calculations

[Nash, Anthony]

Hi Justin,

Thanks for the swift reply. I am now convinced that the 'plateau' approach is the best. 

This brings me to a further concern. Take two systems (umbrella sampling along a reaction coordinate etc..) which are identical in composition, yet different in conformation e.g., a helical tilt/crossing angle.

I run the two systems along a reaction coordinate until I have a plateau. Both systems exhibit a plateau at around 6.5 - 7.5 nm. However, there is a significant difference in free energy (as a theoretical example 50 kj mol^-1 difference) at this region, which would bias the difference in free energy when I normalise both graphs to zero at this point. Given that both systems are atomistically identical, one would assume an identical free energy when the helices are far enough apart that they don't interact. 

By using g_wham in increments of 1 ns (0-1 ns, 0-2 ns, 0-3 ns,.....), I can see how the curve near the interfacing helices converges, yet at 4+ nm the curve is still dropping. 

Would you advice that I run additional windows near the plateau region so it is identical in value across both systems?

Many thanks
Anthony 
________________________________________
From: gmx-users-bounces at gromacs.org [gmx-users-bounces at gromacs.org] on behalf of Justin Lemkul [jalemkul at vt.edu]
Sent: 19 December 2012 13:23
To: Discussion list for GROMACS users
Subject: Re: [gmx-users] Where to stop with TM protein PMF calculations

On 12/19/12 4:12 AM, Nash, Anthony wrote:
> Good morning,
>
> A bit of a long one I am afraid.
>
> I am simulating a transmembrane dimer, and calculating the association free energy through potential of mean force calculations as a function of interhelical distance. I have got very good umbrella coverage along my reaction coordinate, however, I would like to know where I should stop calculating and normalise to zero? I am comparing two systems of identical composition but with a different conformation. Hence, the normalising step is vital.
>
> Looking in the literature the cut offs used seem to be around 2nm (two JACS papers come to mind). I've pulled as far as needed to reach a plateau in the PMF curve. This is around 6.5 nm to 8 nm. The problem is if I normalise at the plateau (which seems to be the beliefs of the professors at my institute), the comparative free energy between the two systems is very different to normalising at 2nm (note: the papers in literature do not observe a plateau, there is still an obvious upward trend at their cutoff).
>
> I need justification of where to normalise. According to literature they cut off outside of interaction range. I have decided to test umbrella windows along the reaction coordinate by decomposing the energy of the system. LJ short and Coul short between peptides at a reaction coordinate distance of 8 nm, is 0 as expected. However, I am still getting energy values on Coul. recip even when I decompose (although the gromacs literature says I can't) the system by systematically setting every charge to zero, i.e.,
>
>
> Peptide A (A), Peptide B (B), SOL, lipid and counter ions have zero charge throughout.
> E_coul_recip_A_B = E_coul_recip_(A_B,A_A,B_B) - (E_coul_recip_A_A + E_coul_recip_B_B)
>
> All charges are zero, with the exception to peptide A in E_coul_recip_A_A and peptide B in E_coul_recip_B_B. and peptide A and B in E_coul_recip_(A_B,A_A,B_B). When I add E_coul_recip_A_B to the LJ short and Coul short I do not get an energy of zero.
>
> So:
> 1) The coul-recip decomposition, is this a flawed approach? The gromacs manual says it can't be decomposed.

It's not trivial to do.  There are some posts about really detailed ways that
you might pull it off, but I don't know whether it's worth the effort.  If the
umbrella sampling simulations were conducted with PME, there is never a true
"non-interacting" state.  There is always some finite contribution to long-range
electrostatics terms.

> 2) Where along the reaction coordinate should I cut off?

I would believe the longer distance.  At 2 nm, you still may have induced order
in either lipids or water between your two proteins, thus affecting the forces
between them.  I see no reason to think that if a plateau has not been
established that the results should be correct.  It could be argued that all
you're doing is picking some random point along the reaction coordinate and
calling it "dissociated" for the sake of convenience.

> 3) Why don't authors pull until they reach a plateau?
>

I have seen several published papers that appear to have very arbitrary stopping
points, so just because it's published, doesn't mean it's perfect :)  I believe
in a plateau, a demonstration that, within the caveats inherent to the system,
the interactions between the two molecules of interest are insignificant.

-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