Difference between revisions of "2021 AMBER tutorial 1 with PDBID 1HW9"

Revision as of 23:49, 6 May 2021

Introduction

In this tutorial, we'll be exploring the wonders of molecular dynamics through the use of AMBER 16. AMBER 16 we'll help study the binding characteristics of the ligand to the target protein!

Directory Setup

It's important to stay organized during this process as we'll be generating a number of files during the process. Overall, we'll be working within five directories. They should be created as follows:

mkdir 001_structure
mkdir 002_parameters
mkdir 003_leap
mkdir 004_equil
mkdir 005_production

Before proceeding, double-check that all directories have correctly been made.

1HW9 Structures

Let's move into our directory titled 001_structure.

cd 001_structure

Receptor

We need to retrieve a fresh copy of the receptor before proceeding forward. Even if you've done the DOCK tutorial, you should still retrieve a fresh .pdb just in case Chimera did something wrong during the preparation process. In order to do so, head over to the PDB website and type in 1HW9 into the header. It should take you straight to the webpage. Download the system as a pdb and open it in Chimera. Our focus will be on Chain A of the complex. In order to visualize the chain:

 Select -> Chain -> A

Notice, however, that this does not remove Chains B, C, and D from our view. We need to perform this action manually.

 Select -> Invert (Selected Molecules) -> Actions -> Atoms/Bonds -> Delete

After deletion, the only thing remaining should be Chain A of the original file. It is within your best interest to save this session as you'll need to come back to this file when preparing the ligand. Doing so will prevent you from having to do the above deletion step again.

 File -> Save Session As ->

Now that the session has been saved, we can go ahead and isolate the receptor. While we've isolated Chain A, you'll notice that Simvastatin is still located on the molecule. In order to remove everything besides HMGR, we need to do the following:

 Select -> Residue -> ADP -> Actions -> Atoms/Bonds -> Delete
 Select -> Residue -> HOH -> Actions -> Atoms/Bonds -> Delete
 Select -> Residue -> SIM -> Actions -> Atoms/Bonds -> Delete

With the receptor isolated, go ahead and save the file as

1HW9_fresh.pdb

Ligand

Open the session in Chimera we previously saved that consisted of Chain A with both the receptor and Simvastatin ligand. This time around, we'll be using the session file to isolate the ligand instead of the receptor.

 Select -> Residue -> SIM -> Select -> Invert (Selected Molecules) -> Actions -> Atoms/Bonds -> Delete

We'll need to properly protonate the ligand in order for the AMBER calculations to function:

Tools -> Structure Editing -> Add H
Tools -> Structure Editing -> Add Charge -> (have Amber ff14SB and AM1-BCC selected) -> Ok

Once done, make sure to save your file as:

1HW9_lig_wH.mol2

NOTE: Before proceeding, make sure both the ligand and receptor file you just created have been copied from your local computer to the 001_structure directory in your Seawulf account.

Simulation Parameters

Let's move into the second directory:

cd 002_parameters

Here, we'll be generating force field parameters for our ligand so that AMBER can use them in later calculations. We'll be doing this using antechamber and parmchk2:

antechamber -i ../001_structure/1HW9_lig_wH.mol2 -fi mol2 -o 1HW9_ligand_antechamber.mol2 -fo mol2 -at gaff2 -c bcc -rn LIG -nc -1

Notice the -nc flag at the end of the command line. Make sure that this value corresponds to the protonation state of your ligand! Once the output file is generated, we need to run parmch2 to generate a frcmod file:

parmchk2 -i 1HW9_ligand_antechamber.mol2 -f mol2 -o 1HW9_ligand.am1bcc.frcmod

TLeap System

Let's now move into the third directory:

cd 003_leap

In this stage, we'll be creating files to simulate the ligand and receptor as one system. This way, calculations like binding affinities and energies can be performed on the system as a whole. Fortunately, tleap will do this for us by generating parameter (parm7) and restart (rst7) files. In order to get started, we'll need to create the input file.

vi leap.in

Copy the following code into your new input file:

#!/usr/bin/sh 
###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 ../001_structure/1HW9_fresh.pdb
###load ligand frcmod/mol2
loadamberparams ../002_parameters/1HW9_ligand.am1bcc.frcmod
lig=loadmol2 ../002_parameters/1HW9_ligand_antechamber.mol2
###create gase-phase complex
gascomplex= combine {rec lig}
###write gas-phase pdb
savepdb gascomplex 1HW9.gas.complex.pdb
###write gase-phase toplogy and coord files for MMGBSA calc
saveamberparm gascomplex 1HW9.complex.parm7 1HW9.gas.complex.rst7
saveamberparm rec 1HW9.gas.receptor.parm7 1HW9.gas.receptor.rst7
saveamberparm lig 1HW9.gas.ligand.parm7 1HW9.gas.ligand.rst7
###create solvated complex (albeit redundant)
solvcomplex= combine {rec lig}
###solvate the system
solvateoct solvcomplex TIP3PBOX 12.0
###Neutralize system
addions solvcomplex Cl- 0
addions solvcomplex Na+ 0
#write solvated pdb file
savepdb solvcomplex 1HW9.wet.complex.pdb
###check the system
charge solvcomplex
check solvcomplex
###write solvated toplogy and coordinate file
saveamberparm solvcomplex 1HW9.wet.complex.parm7 1HW9.wet.complex.rst7
quit