5. Analysis
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> mol = loadpdb wt1mg_water.pdb
> saveAmberParm mol wt1mg_dry.parm7 wt1mg_dry.rst
> quit
The ptraj input file for cleaning up the trajectories is:
trajin wt1mg_md.crd
center :1-154
image center familiar
strip :WAT
strip :Cl-
rms first :3-152@CA
trajout wt1mg_dry.crd
go
You can also download these files from the Download page.
Here is a quick and dirty shell script I hacked to do this job:
### take snapshots at 100ps generate reference files
cat << EOF > ptraj.in
trajin wt1mg_dry.crd 1 100 10
trajout wt1mg_ref restart
go
EOF
ptraj wt1mg_dry.parm7 ptraj.in
### calculate the rmsd for each reference structure
foreach ref (`ls wt1mg_ref*`)
set ofile=$ref.rms
cat << EOF > rms.ptraj
trajin wt1mg_dry.crd 1 100
reference $ref
rms reference out $ofile :1-154
go
EOF
ptraj wt1mg_dry.parm7 rms.ptraj
end
### clean up the files and merge them into a single file
cp wt1mg_ref.2.rms rms1
set i=10
while ($i <= 100)
cut -b12-20 wt1mg_ref.$i.rms > tmp
paste rms1 tmp > rms
cp rms rms1
@ i=$i + 10
end
rm *.rms rms1 tmp
mv rms relax.rms
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:
rms first out wt1mg_ca.b :1-154@CA
atomicfluct out wt1mg.bfactor @CA byatom bfactor
go
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.
III. Visualization
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.