Revision as of 20:43, 5 May 2024

Introduction

DOCK6 is a great tool that helps us understand ligand binding poses and relative binding energies. When DOCK6 is combined with another computational method known as molecular dynamics, we can calculate free energies of binding and probe how the ligand behaves in the binding pocket. In this tutorial we will walk through the use of AMBER16 to perform and analyze molecular dynamic simulations of PDBID 1NDV.

Getting Started

Before starting simulations it is important to create our directories so that we can keep our simulations organized. Create the following directories:

000.paramters
001.tleap_build
002.equilbration
003.production
004.analysis
zzz.master

We will use the zzz.master directory to store our initial mol2 files and other files we may use that do not belong in another directory. Once you have your directories created you can begin to prepare they system for simulation.

Before Simulation

Structure

We will need mol2 and pdb files of the ligand and receptor separately for these simulations. If you followed the 1NDV DOCK6 VS tutorial you will already have these files and can move them to the zzz.master directory. If you did not follow that tutorial the steps will be listed below, but please refer to it if you require a more in depth tutorial. Please note that the receptor file should be without charges and hydrogens, while the ligand will need the charges and hydrogens.

We will start by generating our receptor file. Open the PDB of 1NDV in Chimera and follow these steps.

1) Select an atom on the receptor and use the entire up arrow until the entire protein is selected
2) Select -> Invert (all models)
3) Actions -> Atoms/Bonds -> delete

This should leave us with just the receptor atoms. Save this as a mol2 file and a pdb file as "1ndv_protein_noH_noCharge" DO NOT add charges or hydrogens.

We will now create the ligand file in a similar fashion.

1) Select an atom on the ligand and use the entire up arrow until the entire protein is selected
2) Select -> Invert (all models)
3) Actions -> Atoms/Bonds -> delete

Before we use this file we must add hydrogens and charges.

To add Hydrogens: Tools -> Structure Editing -> Add Hydrogens
To add Charges: Tools -> Structure Editing -> Add Charges (for this ligand we will use an overall charge of 0)

Be sure to reference the PDB file to ensure that Chimera adds hydrogens in the correct places. Once charges and hydrogens are added you can save the mol2 and pdb files as "1ndv_ligand_addH_addCharge" Place both of these files in the zzz.master directory.

Force Field Parameters

Before we can use TLEap, we need to generate forcefield parameters for our ligand. We will use the antechamber program to do this. Move to the 000.paramters directory and use the following command to generate the forcefield parameter file.

antechamber -i ../zzz.master/1ndv_ligand_addH_addCharge.mol2 -fi mol2 -o 1ndv_ligand_antechamber.mol2 -fo mol2 -at gaff -c bcc -rn LIG

When running this you may encounter some warnings such as:

Warning: The number of bonds (3) for atom (ID: 1, Name: N1) does not match the connectivity (2) for atom type (N.ar) defined in CORR_NAME_TYPE.DAT.

Be sure to double check your ligand has the correct bond connectivity, and then these warnings may be safely ignored. You may also encounter some errors that prevent antechamber from running. It may require you to enter the mol2 file and edit the atom types to resolve them. Once antechamber has successfully run, run the following command to make our modified forcefield parameters:

parmchk2 -i 1ndv_ligand_antechamber.mol2 -f mol2 -o 1ndv_ligand.am1bcc.frcmod

TLEap

#!/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/4s0v_built.pdb
 
#THIS IS WHERE YOU WOULD DEFINE DISULFIDE BONDS
#NUMBERING SHOULD MATCH INPUT PDB FILE
#bond rec.Res#.SG rec.Res#.SG

###load ligand frcmod/mol2
loadamberparams ../002.parameters/1ndv_ligand.am1bcc.frcmod
lig=loadmol2 ../002.parameters/1ndv_ligand_antechamber.mol2
 
###create gase-phase complex
gascomplex= combine {rec lig}
 
###write gas-phase pdb
savepdb gascomplex 1ndv.gas.complex.pdb

###write gase-phase toplogy and coord files for MMGBSA calc
saveamberparm gascomplex 1ndv.complex.parm7 1ndv.gas.complex.rst7
saveamberparm rec 1ndv.gas.receptor.parm7 1ndv.gas.receptor.rst7
saveamberparm lig 1ndv.gas.ligand.parm7 1ndv.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 1ndv.wet.complex.pdb
 
###check the system
charge solvcomplex
check solvcomplex
 
###write solvated toplogy and coordinate file
saveamberparm solvcomplex 1ndv.wet.complex.prmtop 1ndv.wet.complex.rst7
quit

Minimization and Equilibration

Production

Analysis

RMSD

Hydrogen Bonding

MM-GBSA