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@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-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
