Amber Workshop - tutorial 8

(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Ross Walker 2005

Case Study: All Atom Structure Prediction and Folding Simulations
of a Stable Protein (Folding Trp-Cage Peptide)

In our simulation we failed to fold up to the crystal structure, however, towards the end of our trajectory we saw a number of major structural changes occurring. Hence it is very possible that we simply haven't given the system long enough to fold to the correct structure.

So out of interest lets see what happens if we start our simulation from the NMR structure. We will start the simulation using the first structure from the NMR pdb file. You should be sufficiently fluent in AMBER now to be able to run these simulations by yourself. First things first we strip the hydrogens from the NMR structure, load it into XLeap and save a prmtop and inpcrd file. The prmtop file should be identical to our original TC5b.prmtop file. Only the inpcrd file should differ. Here it is: nmr_struct.inpcrd.

We run exactly the same simulation as we ran for our linear structure. Here is a tar file containing all of the output files: nmr_struct_output.tar.gz (124 mb)

If we consider the same temperature and energy graphs as we did for our original simulation we see that the temperature is good and that the energies are similar.

Lets try doing a backbone RMS fit of the trajectory to the first structure. This will tell us if there are any big deviations from the NMR structure. We will also calculate an RMS to our lowest energy structure.

trajin equil1.mdcrd.gz
trajin equil2.mdcrd.gz
trajin equil3.mdcrd.gz
trajin equil4.mdcrd.gz
trajin equil5.mdcrd.gz
trajin equil6.mdcrd.gz
trajin equil7.mdcrd.gz
trajin equil8.mdcrd.gz
trajin equil9.mdcrd.gz
trajin equil10.mdcrd.gz
reference nmr_struct.inpcrd
rms reference out nmr_sim_rmsd_to_nmr_struct.dat @N,CA,C

Plot showing backbone RMSD to 1st NMR structure (Black) and lowest energy structure from original folding calculation (red).

There is much more analysis that can be done in order to draw conclusions. I leave it to you the reader to try these out for yourself based on the knowledge you have gained in this tutorial. Try taking a look at the hydrogen bonding, how does this compare. Similarly try running clustering analysis on the NMR structure trajectory and then see what the structures that are closest to the various centroids give.

The fact that the system we looked at in this tutorial failed to fold to the native structure in 50 ns suggests that it may still be kinetically trapped at 325 K. This is contrary to what they state in the paper but then they only ran two simulations at 325 K so would not have had enough data to be sure of this. Thus out of interest you may wish to repeat some of the simulations in this tutorial. Some interesting things to try are:

1) Simulated Annealing
Workshop participants at the PSC Amber workshop 2005 (http://www.psc.edu/training/workshops/2005/Amber/) had great success using this approach. Here you could heat the system to around 375 K, run it for several nano seconds at this temperature and then slowly cool the system back down to 325 K. This avoids the problems with still being kinetically trapped at 325 K since the initial run at 375 K allows the system to repeatedly fold up and unfold so removing any dependence on the starting structure.

2) Langevin Dynamics
Another option you could try is to run with a Langevin thermostat (ntt=3) instead of the Berendsen thermostat. This thermostat works by simulating random collisions, as a molecule in solvent might feel, instead of simply scaling the velocities as the Berendsen thermostat does. This method equilibrates the temperature much more effectively and may allow phase space to be explored quicker. It would be interesting to see if this method yields a folded structure.

3) Replica Exchange MD
This has not been covered in this tutorial but may be added at a later date. In this approach many replicas of the system are run at the same time and every now and then the temperatures are exchanged. In this way a potential energy surface can be explored and mapped quickly.

There are many other things you could try as well such as changing the thermostat relaxation time (tautp) etc.

You could also try analysing the following trajectories. After contacting the authors of the paper I was sent the following two files: carlos.prmtop, carlos.inpcrd. These are the prmtop file used by the authors in the paper and also the structure they labelled their folded state. I have run a few simulations with these files that I provide below so that you can analyse them at your leisure.

1) 50ns from a linear starting structure using tautp=0.5ps and the paper author's prmtop file:

2) 50ns from a linear starting structure using tautp=5.0ps and the paper author's prmtop file:

3) 50ns starting from the lowest energy structure the authors obtained. Using their prmtop file and tautp=0.5ps:

You should have learnt a lot of important techniques from this tutorial. You should also have got an idea of what tricky things protein folding simulations are.

I hope that you leave this tutorial with a much better understanding of how to run more advanced MD simulations and more importantly how to carry out some advanced analysis using AMBER's tools and the MMTSB toolset.