View source for 2015 AMBER tutorial with PARP


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===

[http://ambermd.org/ Amber] - '''A'''ssisted '''M'''odel '''B'''uilding with '''E'''nergy '''R'''efinement - is a multi-program suite for macromolecular simulations. Amber14 is the most recent version of the software and it includes new force fields such as ff14SB. In addition, in this release, more features from sander have been added to pmemd for both CPU and GPU platforms, including performance improvements, and support for extra points, multi-dimension replica exchange, a Monte Carlo barostat, ScaledMD, Jarzynski sampling, explicit solvent constant pH, GBSA, and hydrogen mass repartitioning. Support is also included for the latest Kepler, Titan and GTX7xx GPUs expanded options for Poisson-Boltzmann solvation calculations, accelerated molecular dynamics, additional features in ''sander'' pmemd code, and expanded methods for free energy calculations. Our lab is set up with Amber14 and the latest update of AmberTools14 which contains the programs such as antechamber and tleap to set up your simulation. 

The [http://ambermd.org/doc12/Amber14.pdf Amber 14 Manual] is available to get started with using Amber14. You can search the document for keywords such as "tleap" if you use Adobe Acrobat to view the file. Additionally, [http://ambermd.org/doc12/AmberTools13.pdf AmberTools Reference Manual] is another reference for the programs available under Amber tools. 

Here below are some of the programs available in both Amber and AmberTools:

#'''LEaP''': LEaP is an X-windows-based program that provides for basic model building and Amber coordinate and parameter/topology input file creation. It includes a molecular editor which allows for building residues and manipulating molecules.
#'''ANTECHAMBER''': This program suite automates the process of developing force field descriptors for most organic molecules. It starts with structures (usually in PDB format), and generates files that can be read into LEaP for use in molecular modeling. The force field description that is generated is designed to be compatible with the usual Amber force fields for proteins and nucleic acids. 
#'''SANDER''': Sander is short for Simulated annealing with NMR-derived energy restraints. This allows for NMR refinement based on NOE-derived distance restraints, torsion angle restraints, and penalty functions based on chemical shifts and NOESY volumes. Sander is also the "main" program used for molecular dynamics simulations, and is also used for replica-exchange, thermodynamic integration, and potential of mean force (PMF) calculations. Sander also includes QM/MM capability. 
#'''PMEMD''': This is an extensively-modified version (originally by Bob Duke) of the sander program, optimized for periodic, PME simulations, and for GB simulations. It is faster than sander and scales better on parallel machines.
#'''PTRAJ and CPPTRAJ''': These are used to analyze MD trajectories, computing a variety of things, like RMS deviation from a reference structure, hydrogen bonding analysis, time-correlation functions, diffusional behavior, and so on. 
#'''MM_PBSA and MM_PBSA.py''': These are scripts that automate post-processing of MD trajectories, to analyze energetics using continuum solvent ideas. It can be used to break energies energies into "pieces" arising from different residues, and to estimate free energy differences between conformational basins.

There is also a [http://lists.ambermd.org/mailman/listinfo/amber mailing list] available as an additional resource. What you can do with it is: you document your questions and sent to this mail address, some specialists of Amber will be assigned to reply your email and help you.

===PARP===
Poly (ADP-ribose) polymerase (PARP) is a family of proteins involved in a number of cellular processes involving mainly DNA repair and programmed cell death. PARP is composed of four domains of interest: a DNA-binding domain, a caspase-cleaved domain (see below), an auto-modification domain, and a catalytic domain. The DNA-binding domain is composed of two zinc finger motifs. The PARP family comprises 17 members (10 putative). They have all very different structures and functions in the cell. We're going to use PARP5b in this tutorial.

===Organizing Directories===

It makes things easier to organize your files in a clean and logical way. The following directory structure and naming scheme is a convenient way to organize your files. We could make these directories first before doing anything further. 

 ~username/AMS536/Amber-Tutorial/001.system.prep/
                                 002.tleap/ 
                                 003.pmemd/       
                                 004.ptraj/
                                 005.mmgbsa/

==II. Structural Preparation==


====Antechamber, Parmchk, tLeap====

Before beginning the Molecular Dynamics protocol using AMBER, you must first set up your files.
In your 001.chimera folder, you will add the following files:

  4TKG.lig.mol2
  4TKG.rec.noH.pdb
  4TKG.rec.noH.ZIN.pdb

To prepare the first two files, please see the 2025 DOCK tutorial at the following link: [http://ringo.ams.sunysb.edu/index.php/2015_DOCK_tutorial_with_Poly(ADP-ribose)_polymerase_(PARP)]

Note: delete any headers before the atoms/helix information.

In order to prepare 4TKG.rec.noH.ZIN.pdb, first open 4TKG.rec.noH.pdb and use the command "SHIFT + g" to reach the bottom of the pdb.  The last lines of your files should look like this:

  ATOM   1648  OE1 GLU B 205      23.839 -23.190  57.747  1.00  0.00           O
  TER    1649      GLU B 205
  HETATM 1650 ZN   ZN  B 206      28.130  -3.467  55.482  1.00  0.00          Zn
  END

To make 4TKG.rec.noH.ZIN.pdb, you will need to change the "ZN" atom ID to "ZIN" so that AMBER can read the atom type. 

  ATOM   1648  OE1 GLU B 205      23.839 -23.190  57.747  1.00  0.00           O
  TER    1649      GLU B 205
  ATOM   1650  ZIN ZIN B 206      28.130  -3.467  55.482  1.00  0.00          Zn
  END

==III. Simulation using pmemd==

==IV. Simulation Analysis==

===Ptraj===


===MM-GBSA Energy Calculation===
Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) is a great method to calculate or estimate relative binding affinity of a ligand(s) to a receptor. The binding energy calculated from this method are also known as free energies of binding, where the more negative values indicate stronger binding. For this section, the topology files for the ligand, receptor and complex are needed. 

Create a new directory:
  mkdir 005.MMGBSA
Create an input file name
  vim gb.rescore.in
Enter the following into the input file:
   Single point GB energy calc
 &cntrl
  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,
  imin   = 5,        maxcyc = 1,        ncyc   = 0,
  /
Create a tcsh/bash/csh script (run.sander.rescore.csh) with the following information:
  #! /bin/tcsh
  #PBS -l nodes=1:ppn=1
  #PBS -l walltime=48:00:00
  #PBS -o zzz.qsub.out
  #PBS -e zzz.qsub.err
  #PBS -V
  #PBS -N mmgbsa
  
  set workdir = /nfs/user03/kbelfon/amber_tutorial/005.mmgbsa
  
  cd $workdir
  
  sander -O -i gb.rescore.in \
  -o gb.rescore.out.com \
  -p ../002.tleap/4TKG.com.gas.leap.prm7 \
  -c ../002.tleap/4TKG.com.gas.leap.rst7 \
  -y ../004.ptraj/4TKG.com.trj.stripfit \
  -r restrt.com \
  -ref ../002.tleap/4TKG.com.gas.leap.rst7 \
  -x mdcrd.com \
  -inf mdinfo.com
 
  sander -O -i gb.rescore.in \
  -o gb.rescore.out.lig \
  -p ../002.tleap/4TKG.lig.gas.leap.prm7 \
  -c ../002.tleap/4TKG.lig.gas.leap.rst7 \
  -y ../004.ptraj/4TKG.lig.trj.stripfit \
  -r restrt.lig \
  -ref ../002.tleap/4TKG.lig.gas.leap.rst7 \
  -x mdcrd.lig \
  -inf mdinfo.lig 
  
  sander -O -i gb.rescore.in \
  -o gb.rescore.out.test.rec \
  -p ../002.tleap/4TKG.rec.gas.leap.prm7 \
  -c ../002.tleap/4TKG.rec.gas.leap.rst7 \
  -y ../004.ptraj/4TKG.rec.trj.stripfit \
  -r restrt.rec \
  -ref ../002.tleap/4TKG.rec.gas.leap.rst7 \
  -x mdcrd.rec \
  -inf mdinfo.rec 
  
  exit

Execute this script on the seawulf cluster or machine(s) of your preference
  qsub run.sander.rescore.csh
Three output files will be generated once the job is completed:              
'''gb.rescore.out.com''', '''gb.rescore.out.lig''', and '''gb.rescore.out.rec'''
These files represent the single point energy calculation results for the complex (.com), the ligand (.lig) and the receptor (.rec). The energy will be output by the program Sander for each frame specified in the input file. The final results for one frame in one of the three files should look as the following:

                                    FINAL RESULTS
   
   
    
    NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER
       1       5.9132E+03     2.0005E+01     1.2640E+02     C         159
   
  BOND    =      661.8980  ANGLE   =     1751.7992  DIHED      =     2581.7692
  VDWAALS =    -1696.6585  EEL     =   -13958.9335  EGB        =    -3125.9524
  1-4 VDW =      747.0185  1-4 EEL =     7750.8118  RESTRAINT  =        0.0000
  ESURF   =    11201.4791
 minimization completed, ENE= 0.59132314E+04 RMS= 0.200047E+02

==V. Frequently Encountered Problems==