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

[Justin Lemkul]

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

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