Difference between revisions of "2023 AMBER tutorial 3 with PDBID 2P16"

Revision as of 17:37, 27 April 2023

Introduction

Howdy! If you've gotten this far with your MD project, you've now moved out of the realm of DOCK6 and into the world of molecular dynamics (MD). This tutorial will cover how to use the AMBER 16 software package to simulate your protein-ligand complex on a sub-microsecond timeframe and use the resulting data to calculate binding energies.

Preparing the Structure

Before you start any preparation, you need to keep yourself organized to make sure that you're running the right files and jobs for your MD simulation, otherwise, it may explode.

To do this, you can either copy the example AMBER directory to your scratch directory (this may change per year):

/gpfs/projects/AMS536/2023

Or, you can create the following directories:

mkdir 000.parameters
mkdir 001.tleap_build
mkdir 002.equilibration
mkdir 003.production
mkdir 004.analysis
mkdir zzz.master

For this tutorial, we will be copying the example AMBER setup into our folder to run our complex.

Simulation Parameters

If you're looking at this tutorial after the DOCK6 tutorial and are looking to use the original ligand, then you may want to fetch a fresh PDB file for this. Otherwise, make sure you Cartesian minimize any new ligands before sending them through the gauntlet that is MD simulation and preparation.

All of your original (stock/DOCK6 generated) files should be put into the zzz.master folder, but make sure you're in the proper directory for each step!

Receptor File Generation

To generate a fresh copy of the receptor file, go to the Protein Data Bank (PDB) website and download the PDB structure for 2P16, 2p16.pdb, which will contain the ligand-receptor complex. Open the pdb file in Chimera and go to

 Select -> Structure -> Main

Alternatively, you can instead for this protein:

 Select -> Residues -> All Standard Residues

which will select the receptor protein, and then:

 Select -> Invert (all models)

to select everything except the receptor. These will all be deleted:

 Actions -> Atoms/Bonds -> delete

Now, save your receptor as a PDB. You can save it as a mol2, but you are much more likely to bump into problems down the road:

 File -> Save PDB... 

and save as 2P16_rec.pdb. Close Chimera. You may have thought you were done, but there are still issues with your protein. Specifically, AMBER will not read disulfide bonds between cysteines properly. To fix this, you will need to manually edit your PDB file (using Notepad, vi, nano, etc.) and edit all disulfide-bonded cysteine by changing the residue from CYS to CYX. In addition, there will be some lines to add to your tleap.build.in file later on to make sure that AMBER properly recognizes the disulfide bonds. Save these changes.

Ligand File Generation

To generate the ligand file, open the 2P16.pdb file again, and go to

 Select -> Structure -> Main

Alternatively, you can instead for this protein:

 Select -> Residues -> All Standard Residues

Then you will delete the protein:

 Actions -> Atoms/Bonds -> delete

You should also delete all non-essential cofactors, ions, and waters:

 Select -> Residue -> HOH -> Actions -> Atoms/Bonds -> delete

 Select -> Residue -> CA -> Actions -> Atoms/Bonds -> delete

to delete everything except for the GG2 residue, which is the ligand. After this is completed, you will be adding hydrogens and partial charges:

 Tools > Structure editing > DockPrep

Make sure to select the same parameters as in the DOCK6 prep (ff14SB and AM1BCC). Save this file:

 File -> Save Mol2... 

and save as 2P16_lig_withH.mol2. Close Chimera. Transfer them back into your Seawulf directory (./zzz.master). The next step involved getting your ligands ready for tleap.

Force Field Parameterization

Before running TLeap, we need to generate appropriate force field parameters for the ligand. You can achieve this using antechamber- a program built specifically for parametrization of small molecules, followed by parmchk to verify if it ran correctly. Make sure that the -nc flag corresponds to the net charge of the molecule.

If your ligand happens to be a peptide, you can skip this step.

Run the following commands:

antechamber -i ../zzz.master/2p16_lig_wH.mol2 -fi pdb -o 2p16_ligand_antechamber.mol2 -fo mol2 -at gaff2 -c bcc -rn LIG -nc 0
parmchk2 -i 2p16_ligand_antechamber.mol2 -f mol2 -o 2p16_ligand.am1bcc.frcmod

TLeap

TLeap will generate files describing the topology files (prmtop) and initial parameters (rst7) for the protein, ligand, and complex in gas phase, as well as the solvated (wet) complex. Switch to the 001.tleap_build directory and open a file names tleap.build.in and type the following:

#!/bin/bash

###Load Protein force field
source leaprc.protein.ff14SB
###Load GAFF force field (for our ligand)
source leaprc.gaff
###Load TIP3P (water) force field
source leaprc.water.tip3p
####Load Ions frcmod for the tip3p model
loadamberparams frcmod.ionsjc_tip3p
###Needed so we can use igb=8 model
set default PBradii mbondi3


###Load Protein pdb file
rec=loadpdb ../zzz.master/2p16_rec.pdb

#bond rec.7.SG rec.12.SG 
#bond rec.27.SG rec.43.SG 
#bond rec.156.SG rec.170.SG  
#bond rec.181.SG rec.209.SG 
#bond rec.96.SG rec.109.SG
#bond rec.111.SG rec.124.SG
#bond rec.42.SG rec.58.SG
#bond rec.22.SG rec.27.SG

###Load Ligand frcmod/mol2
loadamberparams ../000.parameters/2p16_ligand.am1bcc.frcmod
lig=loadmol2 ../000.parameters/2p16_ligand.am1bcc.mol2

###Create gas-phase complex
gascomplex= combine {rec lig}

###Write gas-phase pdb
savepdb gascomplex 2p16.gas.complex.pdb

###Write gas-phase toplogy and coord files for MMGBSA calc
saveamberparm gascomplex 2p16.gas.complex.prmtop 2p16.gas.complex.rst7
saveamberparm rec 2p16.gas.receptor.prmtop 2p16.gas.receptor.rst7
saveamberparm lig 2p16.gas.ligand.prmtop 2p16.gas.ligand.rst7

###Create solvated complex (albeit redundant)
solvcomplex= combine {rec lig}

###Solvate the system
solvateoct solvcomplex TIP3PBOX 12.0

###Neutralize system (it will add either Na or Cl depending on net charge)
addions solvcomplex Cl- 0
addions solvcomplex Na+ 0

###Write solvated pdb file
savepdb solvcomplex 2p16.wet.complex.pdb

###Check the system
charge solvcomplex
check solvcomplex

###Write Solvated topology and coord file
saveamberparm solvcomplex 2p16.wet.complex.prmtop 2p16.wet.complex.rst7

You can then run this script and copy over the wet complex prmtop and rst7 files to your local machine:

tleap -f tleap.build.in
scp 'username@login.seawulf.stonybrook.edu:path_to_amber_directory/001.tleap_build/2p16.wet.complex.*' .

You can then visualize these files using the following Chimera commands to ensure the previous steps ran properly:

Tools → MD/Ensemble Analysis → MD Movie

Select file 2p16.wet.complex.prmtop for prmtop section, and then hit add and select the 2p16.wet.complex.rst7 file to visualize the complex.

Equilibration

Production

Analysis