GROMACS Protein Ligand Complex Simulations

1. Preparing Protein-Ligand System

LigParGen server provides OPLS-AAM templates with CM1A/CM1A-LBCC charges for small organic molecules. In general, molecular dynamics simulations are focused on protein/NA-ligand interactions rather than just small molecules. For this reason, in this tutorial, a robust protocol to prepare Gromacs protein/NA-ligand systems using LigParGen server will be explained in detail.

Gromacs, Chimera and python must be installed in your computer to perform this tutorial.

1.1 Obtain a complex

The tutorial uses an example structure of T4 Lysozyme L99A with a Benzene molecule bound (PDB ID: 4W52.pdb, see figure below) which contains just one chain A.

T4 Lysozyme L99A with Benzene

1.2 Cleaning & Preparing Protein for MD simulations

To prepare your system, this tutorial provides the python code below which uses Chimera dockprep functionality to clean the protein and to add missing residues using the Dunbrack rotamer Library. It is important to notice that hydrogens are not added to the protein in this step because it will be done automatically by pdb2gmx to avoid issues with atom names differences. On the other hand, for ligands, hydrogens are added using custom forcefield parameters explicitly.

## 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 the right number of hydrogens are added and the ligand charge matches the expected value.

1.3 Preparing Ligand for MD simulations

Upload the ligand_wH.pdb to LigParGen server to get the parameter files. For protein-ligand simulation, before uploading, make sure that ligand residue number is changed to 1 (BOSS, the core of LigParGen server, only works for a limited number of residues). Before to submit the structure to the server, user must specify two options: optimization and charge model (1.14*CM1A or 1.14*CM1A-LBCC). Optimization will perform a quick minimization of the submitted structure using the OPLS-AA and the charge model parameters generated by the server. If the ligand submitted is not neutral, just CM1A model can be applied (the charge must be specified by the user) and the scale factor 1.14 will not be used. Then, user can submit the structure and download the Gromacs files, 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 perform the protein-ligand complex simulation, ligand coordinates and topology files generated by the server must be combined with the protein coordinates and topology files. First, generate the coordinate file and add hydrogens to the protein using the pdb2gmx gromacs tool with the command below. Please note that this tutorial uses Gromacs 4.6.5.

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

Once you have both files (protein_processed.gro and BNZ.gro), combine them using combineGro_prot_lig.py python script.

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

2.2 Combining Protein-Ligand Topologies

Add the next line #include "BNZ.itp" to topol.top at the top right after the line #include "oplsaa.ff/forcefield.itp" to merge ligand and protein topology files. Then, add the BNZ residue number by adding BNZ 1 after the line Protein_chain_A 1 in the same topol.top file.

2.3 Setup MD simulations

Use the command below to center the complex and place it at least 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 (other water molecule model can be used such as TIP3P or TIP4P)

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 because molecular dynamics simulations require a neutral net charge of the system. User must specify the number of ions (anions and cations) to add according to the system protonation. 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

Once the system is ready, 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

Optional step: In some case, before doing NVT simulations, user needs to add position constraints to ligand so that it won't drift away from protein during equilibration, as the system is heating up.

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

Optional step: 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

Once the system set up is ready, user can perform the equilibration (NVT,NPT) and production (NPT) runs using the command lines shown below.

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