As I have mentioned earlier, analysis of a molecular dynamics trajectory is a highly custom job. It is not practical to anticipate all possible methods of analysis and implement them in a shrink-wrap package (although I have to admit that it does sound like an attractive idea). In the previous section, we have already seen how to use ptraj to process the trajectory files and calculate RMSD.
I. Determination of the initial relaxation time by trace-back RMSD analysis
In this section, we'll revisit the issue of equilibration. Initially, we took an educated guess and equilibrated the system for 100ps. We then examined the potential energy and rmsd and decided that the trajectory was stable and it was okay to begin collecting data. In the paper J. Chem. Phys. Vol 109(23), 10115, Stella et al. pointed out that the common practice of using root mean square deviation from the initial structure to determine the initial relaxation time often over-estimates the true value and resulting in the waste of valuable simulation data. They suggested a better approach to determine how much of the trajectory really belong to the initial relaxation phase and, hence, offers a way to recover some of the data we might have thrown away initially. Their approach consists of using a number of later stage structures as references to calculate the rmsd values and then compare them to the rmsd values calculated using the initial structure as a reference. By tracing the curves backward in time, we can better identify the point at which the trajectory reached a stationary state. We will apply this analysis to the 1-ns trajectory we collected including the first 100-ps of constant presure-constant temperature equilibration trajectory. I have already concatenated the trajectory files into one single crd file and re-imaged it for you here. To make the file size smaller, I have also stripped off water and counter ions from the trajectory. If you are going to do this by yourself, don't forget to use tleap to generate a new topology file for the stripped down trajectory. You can do this in 3 easy steps:1. take the initial topology and restart file and generate a pdb file from them:
2. edit the wt1mg_water.pdb file and delete all the water molecules and Cl- ions.3. use tleap and load up this watered down pdb file and save new topology file
The ptraj input file for cleaning up the trajectories is:
You can also download these files from the Download page.
Here is a quick and dirty shell script I hacked to do this job:
After running the above script, you should get an output file called "relax.rms". Look at the plot below and make your own conclusion ;-)
II. B-factor calculation
Another common analysis is to calculate the B-factors from the MD trajectory and compare it to the crystallographic B-factors. This is really easy to do with ptraj
The ptraj input file is:
This will give you the B-factor of each C-alpha atom. I plotted the result against the experimental B-factors taken from the pdb file 1QS4, chain B (this is the structure we used as our initial model). You probably noticed that the B-factors from the MD simulation is much lower than the crystallographic B-factors, except for a few key residues. This is to be expected given that 1-ns isn't really does not give enough sampling compare to crystal structures.
Finally, what's an MD project without a movie?
There are several free programs that can play back AMBER format trajectories (PyMol, VMD, etc.). I have found that VMD does a very good job of playing back AMBER trajectory as movies and offers quite a number of useful tools for interactively working with the data. I have used VMD to generate a short movie of the 1-ns trajectory here for your enjoyment.
What insights can you get from the MD simulation? Like digging for gold, it's all hidden in there and waiting for you to find out.