(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 2006

MM-PBSA - SECTION 1

MM-PBSA

By Ross Walker & Thomas Steinbrecher

1) Build the starting structure and run a simulation to obtain an equilibrated system.

The system we will model in this simulation is the complex between the human H-Ras protein and the Ras-binding domain of C-Raf1 (Ras-Raf), which is central to the signal transduction cascade. Here is a partially equilibrated, pre-prepared pdb file of the RAS-RAF complex.

ras-raf.pdb

This structure contains the ras and raf proteins and also a physiologically necessary GTP nucleotide as illustrated in the figure below:

For the purposes of this tutorial and for the sake of simplicity we will avoid treating the GTP molecule in the calculation since this would require the setup of new parameters for this compound and is beyond the scope of this tutorial. Thus we will simply remove it from the calculation by erasing it from the pdb file. While not strictly correct this approximation is somewhat reasonable since, as can be seen from the figure above, GTP is not directly involved in the binding interface.

There is also a magnesium ion in this protein that is essentially bound to the GTP molecule so we will remove this as well. Hence you should remove residues 243 and 244 from the pdb file.

The next step is to split this pdb file into the two separate structures such that you have a ras-raf.pdb, a ras.pdb and a raf.pdb. We will use these three structures to create three gas phase prmtop and inpcrd file pairs for the MM-PBSA calculation as well as one for the solvated complex which will be used to run the MD simulations:

> $AMBERHOME/bin/tleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99

Caution: For AMBER 14 please use -f $AMBERHOME/dat/leap/cmd/oldff/leaprc.ff99 in the tleap call for loading the ff99 force field.

com = loadpdb ras-raf.pdb
ras = loadpdb ras.pdb
raf = loadpdb raf.pdb

Make sure you select the correct radii for the calculation method you intend to use. For details see the set PBRadii paragraph in the LEaP section of the manual and the recommendations provided here.

set default PBRadii mbondi2

saveamberparm com ras-raf.prmtop ras-raf.inpcrd
saveamberparm ras ras.prmtop ras.inpcrd
saveamberparm raf raf.prmtop raf.inpcrd

Now before you quit tleap you should create the solvated complex for running the MD simulation:

charge com
> Total unperturbed charge: -0.000000
> Total perturbed charge: -0.000000    (Hence there is no need to add counter ions)
solvatebox com TIP3PBOX 12.0
saveamberparm com ras-raf_solvated.prmtop ras-raf_solvated.inpcrd
quit

Here are the files: ras-raf.prmtop, ras-raf.inpcrd, ras.prmtop, ras.inpcrd, raf.prmtop, raf.inpcrd, ras-raf_solvated.prmtop, ras-raf_solvated.inpcrd

1.1) Equilibrate the solvated complex

We will equilibrate the solvated complex by carrying out a short minimisation, 50ps of heating and 50 ps of density equilibration with weak restraints on the complex followed by 500ps of constant pressure equilibration at 300K. All simulations will be run with shake on hydrogen atoms, a 2 fs time step and langevin dynamics for temperature control. The input files are as follows:

 

Caution: In the examples in this tutorial we do not change the value of the random seed used for the random number generator. This is controlled by the namelist variable ig. This is largely for issues of reproducibility of the results within a tutorial setting. However, when running production simulations, especially when using ntt=2 or 3 (Anderson or Langevin thermostats) it is essential that you change the random number seed from the default on EVERY MD restart. If you are using AMBER 10 (bugfix.26 or later) or AMBER 11 or later you can do this automatically by setting ig=-1 in the cntrl namelist. Otherwise you can specify a positive random number of your choosing for ig each time you restart a calculation. For more details on the pitfalls of not doing this you should refer to the following publication: Cerutti DS, Duke, B., et al., "A Vulnerability in Popular Molecular Dynamics Packages Concerning Langevin and Andersen Dynamics", JCTC, 2008, 4, 1669-1680

You should run all 4 of these simulations using commands along the lines of:

$AMBERHOME/bin/sander -O -i min.in -o min.out -p ras-raf_solvated.prmtop -c ras-raf_solvated.inpcrd \
-r min.rst -ref ras-raf_solvated.inpcrd
$AMBERHOME/bin/sander -O -i heat.in -o heat.out -p ras-raf_solvated.prmtop -c min.rst \
-r heat.rst -x heat.mdcrd -ref min.rst
gzip -9 heat.mdcrd
$AMBERHOME/bin/sander -O -i density.in -o density.out -p ras-raf_solvated.prmtop -c heat.rst \
-r density.rst -x density.mdcrd -ref heat.rst
gzip -9 density.mdcrd
$AMBERHOME/bin/sander -O -i equil.in -o equil.out -p ras-raf_solvated.prmtop -c density.rst \
-r equil.rst -x equil.mdcrd
gzip -9 equil.mdcrd

This takes approximately 5 hours on 16 processors of a 1.7GHz IBM P690.

Here are the output files: equil.tar.gz

Before we proceed with the MM-PBSA production MD run we need to verify that the system has equilibrated. For this we will look at temperature, density, total energy and RMSD. We can start by using the following perl script (process_mdout.pl) that will extract useful information from the output files.

./process_mdout.pl heat.out density.out equil.out

Since the first heating run was performed under constant volume conditions the density data is not recorded. Hence you will need to edit the summary.DENSITY file and remove the first 50 lines. (Because xmgrace is stupid and just gets confused otherwise).

xmgrace summary.DENSITY
xmgrace summary.TEMP
xmgrace summary.ETOT

Additionally we will examine the protein backbone RMSD with respect to the minimized structure to see if conformational stability has been achieved during the equilibration. This can be done using ptraj or cpptraj with the following script:

xmgrace equil.rmsd

The density, temperature and total energy plots have all clearly converged by the end of our equilibration period. The RMSD while seeming to begin to level does not appear to have completely converged but for the purposes of this tutorial is acceptable. In a real calculation you would want to run, depending on your system, significantly more equilibration time. We are now ready to perform the production runs.

---

CLICK HERE TO GO TO SECTION 2

---

(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 2006