From: Aron Broom (broomsday_at_gmail.com)
Date: Mon Jun 17 2013 - 14:13:09 CDT
I'd say from what you've written that the problem is almost certainly
convergence/equilibration. I understand the windows are close together and
build off of the previous simulation, but your system (from my imagining of
what you've said) involves a lot of contacts and therefore a large number
of degrees of freedom. I really doubt that 6.75ns is enough time for those
many degrees of freedom to relax.
You mention that by using the later bits the dG decreases, but didn't get
all the way to a reasonable value. Can you extrapolate on that analysis
and see what dG you might expect for 100ns long windows (that is the
ballpark in which I would expect you to start getting reasonable results)?
Also, what force constant are you using? It seems from your window spacing
like it is extremely high, but you say the dG shouldn't be that high, so
this again suggests to me an equilibration problem. That is, you maybe
initially tried a reasonable force constant like 5 kcal/mol with 0.5
angstrom spacings, but the windows were not becoming evenly distributed
along the reaction coordinate, so you upped the force constant, when really
the problem was slow equilibration, and now that problem still exists, but
is masked by your very high force constant and low spacings.
The best quick test at this point would be to take the a single window that
is somewhere on the steepest part of the "unbinding" portion of your PMF
(that is probably ~ an angstrom from your minimum) and rerun it from a
different starting point (say, the initial equilibration, rather than the
last part of the previous window). Then, take the histogram of the
reaction coordinate from that window and the original (say 20-40 bins or
so), integrate, and see how different they are (i.e. by area). I'd bet
they are more than 5% different, in which case you would be hard pressed to
say with good confidence that they are the same.
On Mon, Jun 17, 2013 at 2:33 PM, Felipe Merino <
> Hi Aron,
> It is like this. We run ~60 ns of unbiased MD and then we started the
> Pull. The windows have a very high constant and are very close together.
> Every window is only 6.75 ns (60 windows) but the centers are spaced by
> only .05 A and every window is started with after 4.5 ns of the last window
> so i do not think the problem is really about convergence (we are really
> pulling slowly). To test it i was trying to discard different lengths as
> equilibration and while the PMF gets lower it definitely does not go to
> "reasonable" values.
> On 17.06.2013 18:51, Aron Broom wrote:
> Hi Felipe,
> A few questions:
> 1) How long are your windows?
> 2) What analysis have you done to see that all the windows have reached
> In general, assuming your initial coordinates were of the bound state,
> then you will overestimate the dG unless everything has come to
> equilibrium. If you start from an unbound state, you'll underestimate the
> dG. If you have done 400ns per window, then I would think you should have
> a lot equilibrated, but if 400ns is the total of all windows, I suspect
> equilibration may be an issue. As a reference, I have a system of a
> protein binding a lipid, and the average window length before that window
> becomes equilibrated is ~60ns (that is with explicit solvent).
> If you are in doubt about equilibration, the easiest thing to do is to
> rerun a couple of windows, from a different starting coordinate and with a
> different langevin seed. Then compare the trajectories in terms of your
> reaction coordinate, you may be surprised by how non-equilibrated many of
> the windows really are. In fact, running pairs of each window in this
> manner, while obviously costly, is a really good way to ensure you're
> results are meaningful. You can also be very extreme and run one set of
> windows starting from the unbound state, but that can be EXTREMELY time
> consuming, so I would suggest just a different bound set of coordinates
> (e.g. equilibrated for an extra 10ns) with a different langevin seed.
> On Mon, Jun 17, 2013 at 11:15 AM, Felipe Merino <
> felipe.merino_at_mpi-muenster.mpg.de> wrote:
>> Hi Giacomo,
>> Thanks for the reply. No, both domains separately have dG's around 20
>> something. The calculations are a bit expensive and that is why we did not
>> go further separating the other domain. In total it is like 400 ns per
>> domain so extending it is not particularly cheap.
>> On 06/17/2013 05:10 PM, Giacomo Fiorin wrote:
>> Hello Felipe, are you summing the free energies from two separate
>> calculations where you detach only one domain while keeping the other
>> domain still bound?
>> If you did that, this would explain why the discrepancy: those two
>> terms are not necessarily additive.
>> The correction would be to finish either simulation by also detaching
>> the bound domain, and letting the entire protein go far enough from the DNA
>> for all interactions to vanish.
>> On Mon, Jun 17, 2013 at 10:43 AM, Felipe Merino <
>> felipe.merino_at_mpi-muenster.mpg.de> wrote:
>>> Dear all,
>>> I know this is a little bit off topic but i think somebody could help us
>>> here. We have been doing some umbrella sampling simulations to calculate
>>> the binding free energy of a protein-DNA complex. We are using the minimal
>>> interatomic distance as reaction coordinate. The thing is that the protein
>>> has two domains and we are only pulling them separately, so in the end we
>>> have always an endpoint with the other domain still bound (that was
>>> planned). In the end, the binding free energies are much higher than
>>> expected (around 25 kcal/mol). The point is that it sounds to me that there
>>> should be a correction on the PMF to account for the "confinement" of being
>>> bound to the other (still attached) domain (very much line in the case
>>> where you restrain the ligand when doing FEP annihilation) but i am not
>>> sure if this is correct. Any insights for this will be highly appreciated.
> Aron Broom M.Sc
> PhD Student
> Department of Chemistry
> University of Waterloo
-- Aron Broom M.Sc PhD Student Department of Chemistry University of Waterloo
This archive was generated by hypermail 2.1.6 : Tue Dec 31 2013 - 23:23:21 CST