2012 AMBER Tutorial with Biotin and Streptavidin

From Rizzo_Lab
Revision as of 10:34, 14 March 2012 by Stonybrook (talk | contribs) (MMGBSA Energy Calculation)
Jump to: navigation, search

For additional Rizzo Lab tutorials see AMBER Tutorials.

In this tutorial, we will learn how to run a molecular dynamics simulation of a protein-ligand complex. We will then post-process that simulation by calculating structural fluctuations (with RMSD) and free energies of binding (MM-GBSA).

I. Introduction


Amber - Assisted Model Building with Energy Refinement - is a suite of about multiple programs for perform macromolecular simulations. Amber11, the current version of Amber, includes newly released functionality such as PMEMD, particle mesh Ewald MD and soft-core Thermodynamics Integration MD. For the tutorial, we are using an older version which is AMBER10.

The Amber 10 Manual is the primary resource to get started with Amber10. (Tip: Using Adobe Acrobat to view the file, you can simply search the document for keywords such as the name of a simulation parameter, which saves much time.) In addition, Amber Tools User's Manual serves as another reference while using Amber tools. You can also read the manual for Amber11 on Amber11 and AmberTools Users' Manuals

Here are some programs in Amber

  1. LEaP: an preparing program for constructing new or modified systems in Amber. It consists of the functions of prep, link, edit, and parm for earlier version of Amber.
  2. ANTECHAMBER: in additional to LEap, this main Antechamber suite program is for preparing input files other than standard nucleic acids and proteins.
  3. SANDER: according to the Amber 10 manual, it is 'a basic energy minimizer and molecular dynamics program' that can be used to minimize, equilibrate and sample molecular conformations. And this is the program we mainly use in this tutorial to generate trajectory files of the molecular system.
  4. PMEMD: version of SANDER that has improved parallel scaling property and optimized speed.
  5. PTRAJ: an analysis program for processing trajectory files. One can use ptraj to rotate, translate the structures, evaluate geometrical features and so on.

There is a mailing list you could sign-up for, as an additional resource.

Biotin and Streptavidin

For information of the Biotin and Streptavidin system, see 2012 DOCK tutorial with Streptavidin.

Organizing Directories

While performing MD simulations, it is convenient to adopt a standard directory structure / naming scheme, so that files are easy to find / identify. For this tutorial, we will use something similar to the following:


II. Structural Preparation

Preparation in Chimera

In this AMBER tutorial, we will use the same system with previous DOCK part. Chimera can directly get the structure by its PDB ID 1DF8. To begin with, we need three files under directory 001.CHIMERA.MOL.PREP.

1. To lower computational cost and make the system clear, we remove chain B of the dimer. Select - chain - B, Action - Atoms/Bonds - delete
2. Remove the water molecules. Select - residue - HOH, Action - Atoms/Bonds - delete

Then we separate the receiver and its ligand.

Select the ligand and delete it. Select - Residue - BTN, Action - Atoms/Bonds - delete

1. Select the protein and delete it. Select - Residue - BTN, Select - Invert, Action - Atoms/Bonds - delete
2. Use Dock Prep to add hydrogens and charges(AM1-BCC) to the ligand. Tools - Structure Editing - Dock Prep


A antechamber input file requires all the atom names to be unique and it only uses the first 3 characters as the name. So if we use 1DF8.lig.chimera.mol2 as the input file, it will cause errors("H102" and "H103" will have the same name "H10").

    14 O3         26.9770   10.6020   12.2050 O.2       1 BTN201     -0.6531
    15 N2         28.6480   12.1210   11.8210 N.pl3     1 BTN201     -0.4789
    16 C4         28.9670   13.2060   10.9010 C.3       1 BTN201      0.0816
    17 H102       32.0358   17.4077   15.9597 H         1 BTN201      0.0272
    18 H103       30.6885   18.0473   15.0102 H         1 BTN201      0.0219
    19 H92        30.5603   15.6243   15.3130 H         1 BTN201     -0.0094
    20 H93        32.1288   15.4384   14.4982 H         1 BTN201      0.0384
    21 H82        29.6624   16.7125   13.2433 H         1 BTN201      0.0275

We need to manually rename the atoms. One way is to use the first column numbers to be the atom names. If you are using Vim, the visual block mode can help by selecting a rectangular section of text. We rename the atoms and save the file as 1DF8.lig.mol2 under directory 001.CHIMERA.MOL.PREP.

    14 O14       26.9770   10.6020   12.2050 O.2       1 BTN201     -0.6531
    15 N15       28.6480   12.1210   11.8210 N.pl3     1 BTN201     -0.4789
    16 C16       28.9670   13.2060   10.9010 C.3       1 BTN201      0.0816
    17 H17       32.0358   17.4077   15.9597 H         1 BTN201      0.0272
    18 H18       30.6885   18.0473   15.0102 H         1 BTN201      0.0219
    19 H19       30.5603   15.6243   15.3130 H         1 BTN201     -0.0094
    20 H20       32.1288   15.4384   14.4982 H         1 BTN201      0.0384
    21 H21       29.6624   16.7125   13.2433 H         1 BTN201      0.0275


To begin with, go to 002.ANTE.TLEAP directory.

To make sure we have access to the three programs that we want to run (antechamber, parmchk and tleap) and we are using the correct version of amber, we can use the which command, type:

which antechamber
which parmchk
which tleap

You should have results similar to this:


Please note that you need to be in the tcsh shell for running the which command shown above. Also, if you are using a different version of amber, you might want to change to the correct version by editing the .cshrc file and source it.

Copy parameters of ions to your working directory from the following resource:

scp -r ~lingling/AMS536/AMBER_Tutorial/002.ANTE.TLEAP/rizzo_amber7.ionparms .

Antechamber is a set of tools to generate files for organic molecules, which can then be read into LEaP. The antechamber program itself is the main program of Antechamber package. It can perform many file conversions, and can also assign atomic charges and atom types. In this tutorial, we use antechamber to convert our input mol2 file into files ready for LEaP. In the command line, type:

antechamber -i ../001.CHIMERA.MOL.PREP/1DF8.lig.mol2 -fi mol2 -o 1DF8.lig.ante.pdb -fo pdb

-i input file name; -fi input file format; -o output file name; -fo output file format. You will have an output file:


Similarly, we can use antechamber to change the fomat of 1DF8.lig.mol2 file to prep file:

antechamber -i ../001.CHIMERA.MOL.PREP/1DF8.lig.mol2 -fi mol2 -o 1DF8.lig.ante.prep -fo prepi

You will get a set of output files:

1DF8.lig.ante.prep         ANTECHAMBER_BOND_TYPE.AC0  NEWPDB.PDB

These files give detailed information of the ligand such as atom coordinates, bond angle, bond length, dihedral angles, etc. For example:


You can see this file contains information related to the ring system.

Parmchk is another program included in the Antechamber package. After running antechamber, we run parmchk to check the parameters. If there are missing parameters after antechamber is finished, the frcmod template generated by parmchk will help us in generating the needed parameters:

parmchk -i 1DF8.lig.ante.prep -f prepi -o 1DF8.lig.ante.frcmod

Then open the output file 1DF8.lig.ante.frcmod, you will see something like this:

remark goes here
c3-o -c -o          1.1          180.0         2.0          General improper torsional angle (1 general atom type)
n -n -c -o         10.5          180.0         2.0          General improper torsional angle (2 general atom types)

This means all the needed force field parameters are avalaible except there are two missing dihedral parameters, which were assigned a default value.

Next, we need create 3 new input files – “tleap.lig.in”, “tleap.rec.in” and “tleap.com.in”, for the ligand, the receptor, and the complex, respectively. These input files will be used by LEaP to create parameter/topology files and initial coordinate files.

vim tleap.lig.in 

In this new ligand input file, type:

set default PBradii mbondi2                                      #set default PBradii                                         
source leaprc.ff99SB                                             #load a force field
loadoff rizzo_amber7.ionparms/ions.lib
loadamberparams rizzo_amber7.ionparms/ions.frcmod                #load an AMBER format parameter file
source leaprc.gaff
loadamberparams rizzo_amber7.ionparms/gaff.dat.rizzo
loadamberparams 1DF8.lig.ante.frcmod
loadamberprep 1DF8.lig.ante.prep                                 #load an AMBER PREP file
lig = loadpdb 1DF8.lig.ante.pdb                                  #load the ligand pdb file
saveamberparm lig 1DF8.lig.gas.leap.prm7 1DF8.lig.gas.leap.rst7  #save the ligand gas phase AMBER topology and coordinate files
solvateBox lig TIP3PBOX 10.0                                     #solvate the receptor using TIP3P, solvent box radii 10 angstroms
saveamberparm lig 1DF8.lig.wat.leap.prm7 1DF8.lig.wat.leap.rst7  #save the ligand water phase AMBER topology and coordinate files
charge lig                                                       #charge the ligand

Similarly, create the tleap.rec.in and tleap.com.in files:

set default PBradii mbondi2
source leaprc.ff99SB
loadoff rizzo_amber7.ionparms/ions.lib
loadamberparams rizzo_amber7.ionparms/ions.frcmod
REC = loadpdb ../001.CHIMERA.MOL.PREP/1DF8.rec.noh.pdb
saveamberparm REC 1DF8.rec.gas.leap.prm7 1DF8.rec.gas.leap.rst7
solvateBox REC TIP3PBOX 10.0                                     
saveamberparm REC 1DF8.rec.wat.leap.prm7 1DF8.rec.wat.leap.rst7
charge REC
set default PBradii mbondi2
source leaprc.ff99SB
loadoff rizzo_amber7.ionparms/ions.lib
loadamberparams rizzo_amber7.ionparms/ions.frcmod
source leaprc.gaff
loadamberparams rizzo_amber7.ionparms/gaff.dat.rizzo
REC = loadpdb ../001.CHIMERA.MOL.PREP/1DF8.rec.noh.pdb
loadamberparams 1DF8.lig.ante.frcmod
loadamberprep 1DF8.lig.ante.prep
LIG = loadpdb 1DF8.lig.ante.pdb
COM = combine {REC LIG}
saveamberparm COM 1DF8.com.gas.leap.prm7 1DF8.com.gas.leap.rst7
solvateBox COM TIP3PBOX 10.0
saveamberparm COM 1DF8.com.wat.leap.prm7 1DF8.com.wat.leap.rst7
charge COM

Use tleap command:

tleap -s -f tleap.lig.in > tleap.lig.out
tleap -s -f tleap.rec.in > tleap.rec.out
tleap -s -f tleap.com.in > tleap.com.out

you will see the following output files:


An alternative way to use antechamber, parmchk and tleap is to make a csh script, which contains the commands you want to do, and run the script.

You can visualize the output files in VMD, and use the output files to run SANDER.


Amber Tools User's Manual

2011 AMBER Tutorial with Biotin and Streptavidin

Visualization in VMD

First, copy “002.ANTE.TLEAP”directory back to Herbie. On Herbie,type the following commands:

cd ~/AMS536/AMBER_Tutorial/
scp –r sw:~/AMS536/AMBER_Tutorial/002.ANTE.TLEAP/ .

Open VMD software:


Then, we can load the gas phase ligand file generated by tleap: VMD Main->File->New Molecule (Please note that in order to visualize in VMD, we need first load the .prm7 file and then load the .rst7 file)

Load file for: New Molecule
File name: 1DF8.lig.gas.leap.prm7
Determine file type: AMBER7 Parm
Load file for: 1DF8.lig.gas.leap.prm7
File name: 1DF8.lig.gas.leap.rst7
Determine file type: AMBER7 Restart

You can see the ligand in the gas phase looks like this:


You can compare this with ligand in the water phase by opening 1DF8.lig.water.leap.prm7, 1DF8.lig.water.leap.rst7 at the same time.


In order to see the ligands clearly, we can visualize without hydrogen atoms. VMD Main->Graphics->Representations->

In the Selected Atoms command line, type: all and noh->Click Apply


You can actually see the ligand in water is very different from that in gas phase.

III. Simulation using sander

Minimization and equilibration


Although we run minimization and production in the same script, they represent distinct stages

The final stage now we have equilibrated our system is to run a production simulation. We should run the production phase of the simulation using the same conditions as the final phase of equilibration to prevent an abrupt jump in the potential energy due to a change in simulation conditions.

First create the file 10md.in:

 10md.in: production (500000 = 1ns)
    imin = 0, ntx = 5, irest = 1, nstlim = 500000,
    temp0 = 298.15, tempi = 298.15, ig = 71287,
    ntc = 2, ntf = 1, ntt = 1, dt = 0.002,
    ntb = 2, ntp = 1, tautp = 1.0, taup = 1.0,
    ntwx = 500, ntwe = 0, ntwr = 500, ntpr = 500,
    scee = 1.2, cut = 8.0, iwrap = 1,
    ntr = 1, nscm = 100,
    restraintmask = ':1-118@CA,C,N', restraint_wt = 0.1,

Then create similar files with the modified parameters:

 06md.in: ntx = 1, irest = 0, nstlim = 50000, restraintmask = ':1-119 & !@H=',
          restraint_wt = 1.0, dt = 0.001,
          ntwx = 1000, ntwr = 1000, ntpr = 1000,
 07md.in: nstlim = 50000, restraintmask = ':1-119 & !@H=', restraint_wt = 0.5, dt = 0.001,
          ntwx = 1000, ntwr = 1000, ntpr = 1000,
 08md.in: nstlim = 50000, restraintmask = ':1-118@CA,C,N', restraint_wt = 0.1, dt = 0.001,
          ntwx = 1000, ntwr = 1000, ntpr = 1000,
 09md.in: nstlim = 50000, restraintmask = ':1-118@CA,C,N', restraint_wt = 0.1, dt = 0.001,
          ntwx = 1000, ntwr = 1000, ntpr = 1000,
 10md.in: nstlim = 500000, restraintmask = ':1-118@CA,C,N', restraint_wt = 0.1, dt = 0.002,
 11md.in: nstlim = 500000, restraintmask = ':1-118@CA,C,N', restraint_wt = 0.1, dt = 0.002,

In The above code we see that we are running for 50000 steps (nstlim = 50000), with a time step of dt where dt of 0.001 is 1fs. We are starting the run with no previous velocity information ntx=1, irest=0. We will print energy output every 500 steps (ntwr=500) and save coordinates every 500 (ntwx=500). We want the initial temperature to be 298.15K (tempi=298.15) and the reference temperature to be 300K (temp0=298.15). Initially, we restraint everything except hydrogens with a restraint weight of 1.0, and with later runs we decrease the restraint weight and only the backbone is restrained. With each future run, we are using velocity information from the previous restart file (ntx=5, irest=1), and switched on constant pressure (ntb=2) with isotropic position scaling (ntp=1).

IV. Simulation Analysis



To run the Ptraj analysis, you should first create a separate folder for it under the root directory if you don't have one.

mkdir  004.PTRAJ

Woo Suk

Finally, you want to perform MMGBSA analysis to obtain energetic information about the system. To do this, you need a trajectory file of the gas phase complex (which we already created above and saved as 1DF8.com.trj.stripfit), and want to create two more trajectory files containing just the receptor and just the ligand. To do so, you have to run two more ptraj using the following ptraj input files:


trajin 1DF8.com.trj.stripfit 1 2000 1
trajout 1DF8.rec.trj.stripfit
strip :119


trajin 1DF8.com.trj.stripfit 1 2000 1
trajout 1DF8.lig.trj.stripfit
strip :1-118

To execute, type in the following commands:

ptraj  ../002.ANTE.TLEAP/1DF8.com.gas.leap.parm  ptraj.4.in  > ptraj.4.log
ptraj  ../002.ANTE.TLEAP/1DF8.com.gas.leap.parm  ptraj.5.in  > ptraj.5.log

To visualize trajectories, we have to move the trajectory files to Herbie with scp or rsync commands and run VMD.

rsync -arv 004.PTRAJ/ ID@compute.mathlab.sunysb.edu:~AMS536/AMBER_tutorial/004.PTRAJ

In VMD, we open one of the .prm7 files in 002.ANTE.TLEAP. If you want to visualize the whole complex in gas state, you can open 1DF8.com.gas.leap.prm7 with AMBER7 Parm from 002.ANTE.TLEAP and then 1DF8.com.trj.stripfit with AMBER coordinates from 004.PTRAJ. With these files, you can look at the real-time movement of the complex in the gas state like following snapshot.


You can repeat this procedure to observe the real-time movement of the complex in the water state. Only different thing is opening 1DF8.com.wat.leap.prm7 instead of 1DF8.com.gas.leap.prm7. The following snapshot is that of movement in the water phase.


In the snapshot above, the water molecules are shown as points.

MMGBSA Energy Calculation

MM/GBSA is the acronym for Molecular Mechanics/Poisson-Boltzmann Surface Area. This part of AMBER combines molecular mechanics (MM) with both the electrostatic (PB) and nonpolar (SA) contribution to solvation . Topology files are needed for the receptor (R), ligand (L), and receptor-ligand complex (RL). The trajectory files (trj) generate the coordinates. Therefore, molecular dynamics is used to generate a set of snapshots taken at fixed intervals from the trajectories. These snapshots are processed to remove solvent and generate the free energy of binding.

In the AMBER_Tutorial directory, create a new directory:

mkdir 005.MMGBSA

In this new directory, create the file gb.score.in containing:

Single point GB energy calc
 ntf    = 1,       ntb    = 0,       ntc  = 2,
 idecomp= 0,       
 igb    = 5,       saltcon= 0.00,
 gbsa   = 2,       surften= 1.0,
 offset = 0.09,    extdiel= 78.5,
 cut    = 99999.0, nsnb   = 99999,
 scnb   = 2.0,     scee   = 1.2,     
 imin   = 5,       maxcyc = 1,       ncyc   = 0,

Then create a csh script, run.sander.rescore.csh, that contains the following lines of command:

#! /bin/tcsh
#PBS -l nodes=1:ppn=2
#PBS -l walltime=24:00:00
#PBS -o zzz.qsub.out
#PBS -e zzz.qsub.err
set workdir = “${HOME}/AMBER_Tutorial/005.MMGBSA”
cd ${workdir}

sander -O -i gb.rescore.in \
-o gb.rescore.out.com \
-p ../002.TLEAP/1df8.com.gas.leap.prm7 \
-c ../002.TLEAP/1df8.com.gas.leap.rst7 \
-y ../004.PTRAJ/1df8.com.trj.stripfit \
-r restrt.com \
-ref ../002.TLEAP/1df8.com.gas.leap.rst7 \
-x mdcrd.com \
-inf mdinfo.com

sander -O -i gb.rescore.in \
-o gb.rescore.out.lig \
-p ../002.TLEAP/1df8.lig.gas.leap.prm7 \
-c ../002.TLEAP/1df8.lig.gas.leap.rst7 \
-y ../004.PTRAJ/1df8.lig.trj.stripfit \
-r restrt.lig \
-ref ../002.TLEAP/1df8.lig.gas.leap.rst7 \
-x mdcrd.lig \
-inf mdinfo.lig

sander -O -i gb.rescore.in \
-o gb.rescore.out.rec \
-p ../002.TLEAP/1df8.rec.gas.leap.prm7 \
-c ../002.TLEAP/1df8.rec.gas.leap.rst7 \
-y ../004.PTRAJ/1df8.rec.trj.stripfit \
-r restrt.rec \
-ref ../002.TLEAP/1df8.rec.gas.leap.rst7 \
-x mdcrd.rec \
-inf mdinfo.rec


Then this script should be sent to the queue, i.e., qsub the script using the commands:

chmod +x run.sander.rescore.csh
qsub run.sander.rescore.csh

When the job is complete, you will obtain the following output files:

gb.rescore.out.com gb.rescore.out.lig gb.rescore.out.rec

In these files, the single point energy calculation results will be written for each individual frame. It will be found in the results section and the output file will have an infrastrucutre that is similar to the following:

                    FINAL RESULTS

   NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER
      1       7.9662E+03     1.9504E+01     1.3272E+02      C        	3403

BOND    =      836.7674  	ANGLE   =     2343.1208  	DIHED      =     2917.9356
VDWAALS = -2299.0570  	EEL     =   -19265.6215  	EGB        =    -3539.4050
1-4 VDW =     1071.1234  	1-4 EEL =    11631.2887  	RESTRAINT  =        0.0000
ESURF   =    14270.0746
minimization completed, ENE= 0.79662270E+04 RMS= 0.195036E+02
minimizing coord set #     2

V. Frequently Encountered Problems


1.If you have a problem opening the structure of your molecule using VMD, be certain that the prm7 file corresponds to the rst7 file. For example, 1df8.lig.wat.leap.prm7 must be matched to 1df8.lig.wat.leap.rst7.

2. The VMD Console window (often hidden behind the visualization screen of your molecule) can be used to determine where VMD has found an error in your typed commands if you are unable to visualize your molecule.

3. If your "rescore" files do not appear when you are running the script: run.sander.rescore.csh, it may be the result of there not being a free node. You job may be sitting on the queue. If you type qstat -u username, you will be able to see if the job is sitting there because a Q will appear.

Woo Suk