I have performed a 200 ns MD simulation (Gromacs - Amber forcefield) using PLUMED meta-dynamics with a distance based reaction coordinate of protein A entering protein B. I do see some noticeable minima, protein on the outside, and protein entered in the catalytic. There is a DeltaE of about 12 kcal/mol. I would like to prove that the exploration of these two minima along the reaction coordinate has converged. I understand I can:
a) normalise to zero
and,
b) output accumulative deposits of potential energy gaussian as a free energy plot to see the the free energy plot eventually settle to the final summed free energy plot.
This is plotted on a x-y (distance-energy) graph. Is there any other representation of converged data that I can use?
Many thanks
1. Plot RMS values against time and see if a stable value is reached.
2. Count number of hydrogen bonds - see if that value stabilizes.
3. Check secondary structure throughout the trajectory - see if it converges to a constant value.
4. Likewise with core packing
Since metadynamics is inherently a non equilibrium method, convergence would predominantly be measured by the free energy plot. If the free energy plot does not change significantly (within some cutoff amount) by additional simulation, then you are considered to be converged. There are several other things to look at though:
1. Errors using block averaging of the free energy surface. You can break up the simulation into blocks and determine the variance within each of these blocks. This should become relatively constant after long simulations.
2. Once converged, your system should diffuse along the reaction coordinate, meaning that it will visit all space with identical probability since you have added a biasing potential. Plot the reaction coordinate and see if you get diffusive behavior.
Hi both, thanks for your comments. After reading both comments and giving it some thought perhaps I ought to rephrase my question. Rather than "equilibrate", how can I demonstrate that my system has REVISIT conformations along the reaction coordinate?
For example, between energy wells along a distance based reaction coordinate (d(r)), visiting between sites A and B might also encompass a site C, which although has a d(r) equal to A, the overall protein conformation could be very different, thus affecting the free energy surface plot. I understand that a careful selection of multiple reaction coordinates helps prevent this e.g., a distance and an angle reaction coordinate, but adding additional reaction coordinate parameters drastically slows down the MD calculations.
As per the first commenters suggestions, I was thinking of calculating the RMSD of the hydrogen bonds over the course of the simulation at fixed d(r) along the reaction coordinate, hoping they would consistent per d(r).
Geoffrey - I particularly liked your comment regarding the equal probability for the system to diffuse along the reaction coordinate. I've just got to think of a nice way to demonstrate this.
Proteins have many degenerate states, so you will never converge in terms of all the microstates. I would really only focus on the free energy plot. If this does not change significantly with more simulation then you have sufficiently sampled orthogonal states, and are not likely to overcome additional barriers by these states.
Using accumulated (with time) FES values e.g., 0-10ns, 0-20ns, 0-30ns., etc., I have generated:
1) A distance d(r) reaction coordinate vs FES plot. The last few accumulated FES plots (see above) had 'converged'.
2) For each position along the reaction coordinate d(r), I plot the value from each accumulated FES set. All the plots i.e., each d(r) position, flat lines.
- Badges
- Science topic