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.
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.
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 lig_name = sys.argv 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"
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.
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
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.
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
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
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
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
define = -DPOSRES ; position restrain the protein and ligandfrom the NVT and NPT equilibration scripts.