Gromacs Protein Ligand Complex Simulations

Gromacs Protein Ligand Complex Simulations

1. Preparing Protein-Ligand System

1.1 Obtaining Protein

Choose the protein of choice with ligand in it. In this example 4W52.pdb is used. It is a structure of T4 Lysozyme L99A with Benzene Bound. Since this structure has only one chain, nothing has been done. If there are more than one chain, keep only chain A/B to do MD simulations or FEP simulations further.

1.2 Cleaning & Preparing Protein for MD simulations

Code below uses dockprep functionality of Chimera to clean the protein and to add missing residues using Dunbrack Library. Hydrogens are not added to the protein as it was done automatically by pdb2gmx and to avoid any issues with differences in atom names. Hydrogens are added to Ligand as we have forcefield parameters explicitly defined for ligand atom names.

## PDB_FILE SHOULD THE COMPLETE PATH OF THE FILE
## REPLACE BNZ with LIGAND resname 
## USAGE: Chimera --nogui --script "prep_prot_lig.py 4w52.pdb BNZ" 
import chimera
from DockPrep import prep
import Midas
import sys
import os
PDB_file = sys.argv[1] 
lig_name = sys.argv[2]
os.system('grep ATOM %s > %s_clean.pdb'%(PDB_file,PDB_file[:-4]))
os.system('grep %s  %s > %s.pdb'%(lig_name,PDB_file,lig_name))
protein=chimera.openModels.open('%s_clean.pdb'%PDB_file[:-4])
ligand=chimera.openModels.open('%s.pdb'%lig_name)
prep(protein,addHFunc=None,addCharges=False)
prep(ligand)
Midas.write(protein,None,"protein_clean.pdb")
Midas.write(ligand,None,"ligand_wH.pdb")

Run this code by using Chimera. Make sure you have Chimera installed and can be called from Command Prompt.

Chimera --nogui --script "prep_prot_lig.py 4w52.pdb BNZ"

This produces protein_clean.pdb and ligand_wH.pdb that will be used further. Make sure that right number of hydrogens are added and the charge on ligand matches the expected value.

1.3 Preparing Ligand for MD simulations

Upload the ligand_wH.pdb to LigParGen server and for protein-ligand simulation, before uploading make sure that ligand residue name is changed to 1. This is because BOSS, the core of LigParGen server only works for only a certain number of residues. Download the files for Gromacs, i.e BNZ.gro (coordinate file)and BNZ.itp (topology) files.

2. Setting up Gromacs simulations

2.1 Combining Protein-Ligand coordinates

In order to simulate the protein-ligand complex, one needs to combine both coordinates and topology.

First add hydrogens to protein using this command. Please note that tutorial is made using Gromacs/4.6.5.

pdb2gmx -f protein_clean.pdb -o protein_processed.gro -water spce 

Once you have both protein_processed.gro and BNZ.gro combine them using combineGro_prot_lig.py code.

python combineGro_prot_lig.py protein_processed.gro BNZ.gro > complex.gro

2.2 Combining Protein-Ligand Topologies

First modifications is to add #include "BNZ.itp" to topol.top at the top right after the line #include "oplsaa.ff/forcefield.itp" and second you need to add the number of BNZ reisdues by adding BNZ 1 after the line saying Protein_chain_A 1

2.3 Setup MD simulations

Use the command below to center the complex and place it atleast one angstrom from the center of water box.

editconf -f complex.gro -o complex_box.gro -c -d 1.0 -bt cubic

Fill the cubic box with SPC/E water using the following command

genbox -cp complex_box.gro -cs spc216.gro -o complex_box_wSPCE.gro -p topol.top

Next step is to neutralize the system by adding ions, but one needs to get the net charge and add anions/cations accordingly. To figure out the net charge, run the command below. You can obtain ions.mdp here.

grompp -f MDP/ions.mdp -c complex_box_wSPCE.gro -p topol.top -o ions.tpr

In this case, net charge of system is +8 and to neutralize add 8 Cl- ions using the command below.

genion -s ions.tpr -o complex_box_wSPCE_ions.gro  -p topol.top -pname NA -nname CL -nn 8

Minimization

Minimize the energy using em.mdp gromacs command file.

grompp -f MDP/em_real.mdp -c complex_box_wSPCE_ions.gro -p topol.top -o em.tpr
mdrun -v -deffnm em

Before doing NVT simulations one needs to add position constraints to ligand so that it wont drift away from protein during equilibration, as the systme is heating up.

genrestr -f BNZ.gro -o posre_BNZ.itp -fc 1000 1000 1000

Once posre_BNZ.itp is generated add #include "posre_BNZ.itp" just above the line #include "oplsaa.ff/ions.itp"

To make sure that temperature is calculated by looking at both protein ligand system together you need to combine them together using following command

make_ndx -f em.gro -o index.ndx

Equilibration and Production Run

1. NVT Equilibration with Position restraints
2. NPT Equilibration with Position restraints
3. NPT Production run with no position restraints

Run the equilibration and production runs using using the LSF script as shown below.

#!/bin/bash
#BSUB -o .                       ## requests outfile be generated
#BSUB -a openmpi
#BSUB -q shared                  ## specify queue to run on, shared = default
#BSUB -J FAM           ## job-name
#BSUB -W 24:00                   ## walltime <HH:MM>= 24 hrs
#BSUB -n 4 
#BSUB -R "rusage[mem=50000]"     ## memory limit= MB count only (1024MB = 1GB)
module load Apps/Gromacs/4.6.5

grompp -f nvt.mdp -c em.gro -p topol.top -n index.ndx -o nvt.tpr
mpirun.lsf mdrun -deffnm nvt
grompp -f npt.mdp -c nvt.gro -p topol.top -n index.ndx -o npt.tpr
mpirun.lsf mdrun -deffnm npt
grompp -f md.mdp -c npt.gro -p topol.top -n index.ndx -o md.tpr
mpirun.lsf mdrun -deffnm md