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

3. Upload this file to Server and get the BNZ.gro and BNZ.top files This is topology and position files.

4. Follow the instructions from Bevan lab Lysozyme tutorial for Gromacs Step1# module load Apps/Gromacs/4.6.5 Step2# pdb2gmx -f protein_clean.pdb -o protein_processed.gro -water spce Step3# Combine BNZ.gro and protein_processed.gro to form complex.gro python combineGro_prot_lig.py protein_processed.gro BNZ.gro > complex.gro Step4# Add BNZ.itp at the top and "BNZ 1" at the bottom Step3# editconf -f complex.gro -o complex_box.gro -c -d 1.0 -bt cubic This puts protein at the center, 1nm from the edge of the box Step4# genbox -cp complex_box.gro -cs spc216.gro -o complex_box_wSPCE.gro -p topol.top This command solvates it with SPCE_water Step5# grompp -f MDP/ions.mdp -c complex_box_wSPCE.gro -p topol.top -o ions.tpr Get the net charge of system to neutralize it in next step Step6# genion -s ions.tpr -o complex_box_wSPCE_ions.gro -p topol.top -pname NA -nname CL -nn 8 Step7# grompp -f MDP/em_real.mdp -c complex_box_wSPCE_ions.gro -p topol.top -o em.tpr Step8# mdrun -v -deffnm em Step9# genrestr -f jz4.gro -o posre_jz4.itp -fc 1000 1000 1000 This generates position restrainsts to ligand Step10# Add position restraints to topol.top file right after posre.itp line Step11# make_ndx -f em.gro -o index.ndx Merge Protein and BNZ groups for NVT equilibration phase Step12# grompp -f nvt.mdp -c em.gro -p topol.top -n index.ndx -o nvt.tpr Step13# mdrun -deffnm nvt Step14# grompp -f npt.mdp -c nvt.gro -p topol.top -n index.ndx -o npt.tpr Step15# mdrun -deffnm npt Step16# grompp -f md.mdp -c npt.gro -p topol.top -n index.ndx -o md.tpr Step17# mdrun -deffnm md