Announcement

BIRLA INSTITUTE OF SCIENTIFIC RESEARCH
Bioinformatics Workshop
on
"Bioinformatics Application in modelling and Drug Designing"
15 - 17 February 2018


Day3

Go Back

Section - 1

Molecular Simulation Using Gromacs

Check if GROMACS is installed.

How? Type: gmx on terminal.

Downloading Protein structure file

We must download the protein structure file we will be working with. For this tutorial, we will utilize T4 lysozyme L99A/M102Q (PDB code 3HTB). Go to the RCSB website and download the PDB text for the crystal structure.

Clean the PDB file

Command to extract protein
grep '^ATOM' 3htb.pdb > protein.pdb

Command to extract ligand
grep '^HETATM' 3htb.pdb | grep 'JZ4' > ligand.pdb

If you want to make a complex pdb cat protein.pdb ligand.pdb > complex.pdb


Prepare the protein topology with pdb2gmx

Now that the crystal waters are gone and we have verified that all the necessary atoms are present, the PDB file should contain only protein atoms, and is ready to be input into the first GROMACS module, pdb2gmx. The purpose of pdb2gmx is to generate three files:
1. The topology for the molecule.
2. A position restraint file.
3. A post-processed structure file.

Execute pdb2gmx by issuing the following command:

gmx pdb2gmx -f protein.pdb -o protein_proc.gro -water spc

    Prepare the ligand topology using external tools

    Go to the PRODRG website:
    (a) http://davapc1.bioch.dundee.ac.uk/cgi-bin/prodrg
    (b) Upload ligand.pdb file.

    Chirality (yes/no): maintain chiral centers in the input molecule (yes), or reconstruct them (no). In our case, we have no chirality, so leave this option set to its default.

    Charges (full/reduced): full charges are for condensed-phase systems (43A1 force field); reduced are for "vacuum" simulations (43B1 force field). Use full charges for nearly all applications. The accuracy of 43B1 is debatable, anyway.

    EM (yes/no): perform energy minimization (yes), or leave coordinates alone (no). In our case, we want to construct our protein-ligand complex from existing coordinates, so we do not want PRODRG to minimize our molecule in vacuo. Choose "no" for this option.

    Edit the topology and itp of ligand generated.

    Including the parameters for the JZ4 ligand in the system topology is very easy. Just insert a line that says #include "drg.itp" into topol.top after the position restraint file is included. The inclusion of position restraints indicates the end of the "Protein" moleculetype section.

    Concatenate protein and ligand topologies into complex, after add 15 to complex.gro header (2nd line).
    cat protein_proc.gro jz4.gro > complex.gro Insert #include "drg.itp" into topol.top after the position restraint file is included. After line number 10647.

    TOP

    Define the box

    Let's define the box using editconf:

    At this point, the workflow is just like any other MD simulation. We will define the unit cell and fill it with water.

    gmx editconf -f 1AKI_processed.gro -o 1AKI_newbox.gro -c -d 1.0 -bt cubic

    Now that we have defined a box, we can fill it with solvent (water). Solvation is accomplished using solvate:
    gmx solvate -cp 1AKI_newbox.gro -cs spc216.gro -o 1AKI_solv.gro -p topol.top

    Adding charges to the system

    Download em.mdp

    Use grompp to assemble a .tpr file, using any .mdp file.
    gmx grompp -f em.mdp -c solv.gro -p topol.top -o ions.tpr

    Now we have an atomic-level description of our system in the binary file ions.tpr. We will pass this file to genion:
    gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -nn 6

    Energy Minimization Step


    Now that the system is assembled, create the binary input using grompp using this input parameter file:

    First download em_real.mdp

    Run command
    gmx grompp -f em_real.mdp -c solv_ions.gro -p topol.top -o em.tpr


    Run the minimization

    We are now ready to invoke mdrun to carry out the EM:


    gmx mdrun -v -deffnm em
    Observe the convergence of the system

    Equilibration


    Equilibrating our protein-ligand complex will be much like equilibrating any other system containing a protein in water. There are a few special considerations, in this case:
    Applying restraints to the ligand
    Treatment of temperature coupling groups

    Restraining the Ligand

    PRODRG did not generate a file analogous to posre.itp for our ligand, but GROMACS provides the means to do so with the genrestr module. Simply run genrestr on the jz4.gro file that we obtained from PRODRG:
    gmx genrestr -f jz4.gro -o posre_jz4.itp -fc 1000 1000 1000

    Now, we need to include this information in our topology. We can do this in several ways, depending upon the conditions we wish to use. If we simply want to restrain the ligand whenever the protein is also restrained, add the following lines to your topology in the location indicated: ; Include Position restraint file #ifdef POSRES #include "posre.itp" #endif ; Include ligand topology #include "drg.itp" ; Ligand position restraints #ifdef POSRES #include "posre_jz4.itp" #endif ; Include water topology #include "gromos43a1.ff/spc.itp"

    Thermostats

    Proper control of temperature coupling is a sensitive issue. Coupling every moleculetype to its own thermostatting group is a bad idea.

    Make index using commands

    gmx make_ndx -f em.gro -o index.ndx

    Merge the "Protein" and "JZ4" groups with the following, where ">" indicates the make_ndx prompt: > 1 | 13 > q We can now set tc_grps = Protein_JZ4 Water_and_ions to achieve our desired "Protein Non-Protein" effect.

    Equilibration, Part 2

    Once the NVT simulation is complete, proceed to NPT with this .mdp file:
    gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p topol.top -n index.ndx -o npt.tpr
    gmx mdrun -deffnm npt

    Copy nvt.mdp file
    Run command gmx grompp -f nvt.mdp -c em.gro -p topol.top -n index.ndx -o nvt.tpr
    gmx mdrun -v -deffnm nvt

    Once the NVT simulation is complete, proceed to NPT with this .mdp file:


Go Back

TOP