2020 AMBER tutorial with PDBID 3VJK

In this tutorial, we will be modeling the dynamics of the ligand with the receptor using AMBER 16. Amber is a molecular dynamics simulation software package.

Generating Parameters for the simulation

In order to utilize Amber for molecular dynamic, parameters for the bio molecules will be needed. Luckily, there have been years of parameter development so parameters for the protein do not have to worried about. However, the small ligand does not have parameters in the standard protein force field. Consequently, we will need to generate a fcmod file specific for the ligand.

We will need to have some structures for running the simulation. This can be all stored in a directory called zzz.master. In this directory we will have the following files:

3vjkFH.pdb  3VJK_hydrogen_protein.mol2  3VJK_ligand_hydrogens.mol2

Please note that 3VJK_hydrogen_protein.mol2 had to be converted to a pdb using tleap because there were incapability issues that arose in a later step. To convert you do the following:

  tleap
  rec= loadmol2 3VJK_hydrogen_protein.mol2 
  savepdb rec 3vjkFH.pdb  
  quit

First we will make a specific directory dedicated for the generation of the parameters for the ligand:

  mkdir 000.parameters

In this directory we will run the following command:

  antechamber -i ./../zzz.master/3VJK_ligand_hydrogens.mol2 -fi mol2 -o 3vjk_ligand_antechamber.mol2 -fo mol2 -at gaff2 -c bcc -rn LIG -nc 2

notice that will be using gaff2. This stands for general amber force field 2. This will allow us to parameterize ligands for simulations. Additionally, the ligand has a charge of +2 and that was noted with the -nc flag. Once this command has run, it is beneficial to check to find if any parameters are missing:

   parmchk2 -i 3vjk_ligand_antechamber.mol2 -f mol2 -o 3vjk_ligand.am1bcc.frcmod

The output shows that the parameters that were generated were sufficient to continue.

Build system with TLeap

Before we can simulate the protein and ligand complex, we must build the whole system together. This involves adding solvent and solvent ions to the protein and ligand complex. In order to accomplish this we will be using tleap.

We will make a new directory for the system:

  mkdir 001.tleap_build

now, we will make a tleap.in file--which will contain information about building the system: vim tleap.in

    #!/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 ./../zzz.master/3vjkFH.pdb
    ##@make disulfide bonds
    bond rec.328.SG rec.339.SG
    bond rec.385.SG rec.394.SG
    bond rec.444.SG rec.447.SG
    bond rec.454.SG rec.472.SG
    bond rec.649.SG rec.762.SG
    ###load ligand frcmod/mol2
    loadamberparams ./../000.parameters/3vjk_ligand.am1bcc.frcmod  
    lig=loadmol2 ./../000.parameters/3vjk_ligand_antechamber.mol2

    ###create gase-phase complex
    gascomplex= combine {rec lig}

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

Note that the number of ions of Cl- and Na+ was set to zero. In this situation, we allowed tleap to use a grid to calculate the ions that are needed to neutralize the system. Later on, we will see that amber correctly added the about of ions to set the system equal to zero. This is an important condition to look at. The system can crash if charges are not resolved. The following files should have been created:

     3vjk.complex.parm7  3vjk.gas.complex.pdb  3vjk.gas.complex.rst7  3vjk.gas.ligand.parm7  3vjk.gas.ligand.rst7  3vjk.gas.receptor.parm7  3vjk.gas.receptor.rst7  3vjk.wet.complex.parm7  3vjk.wet.complex.pdb    3vjk.wet.complex.rst7

Now that the system is built, it is important to visualize the system to make sure that everything was generated correctly. You may do this using chimera or VMD. Please note that chimera has a special way of loading parm7 and rst7 files. You must go to trajectory and select the two files in order to visualize. Below are the images that were generated using VMD