From Docking to Dynamics: How to set up a docking generated protein-ligand complex for Molecular Dynamics simulation

Create a separate folder where the tasks will be performed.

You should have following softwares installed
UCSF Chimera
AmberTools15
Acpype
Gromacs 4.6.2 or later versions

Use the docked complex to generate separate rec.pdb and big.mol2 files.

Open the complex.pdb with UCSF Chimera and Select > Residue > LIG and then remove the ligand by Actions > Atoms/Bonds > Delete . Now save the receptor in PDB format and give it a name rec.pdb

Now again open the Complex and Save the ligand in the .mol2 format. Make sure your ligand has hydrogens.

2. H++ server http://biophysics.cs.vt.edu/

H++ server assigns correct protonation state for the receptor. Choose the protonation state at your chosen pH.

Please uncheck “correct orientation of ASN, GLN and HIS groups, add H atoms, and assign HIS H atoms to the δ or ε O, based on van der Waals contacts and H-bonding.” this option if you don’t wanna flip the orientation of some residues.

Save the .top and .crd file from the output.

3. Create a PDB with correct protonation state:

~/amber14/bin/ambpdb -p 0.15_80_10_pH6.5_rec.top <0.15_80_10_pH6.5_rec.crd> rec_final.pdb

the rec_final.pdb will be the receptor file which we will use.

4. Create the Ligand topology and .frcmod files:

-nc is the charge of the ligand.

antechamber -i lig.mol2 -fi mol2 -o LIG.mol2 -fo mol2 -j 4 -at gaff -c bcc -nc 0

parmchk -i LIG.mol2 -f mol2 -o LIG.frcmod

5. Use tleap of AmberTools15

source leaprc.ff14SB
loadAmberParams frcmod.tip4pew
loadAmberParams frcmod.ionsjc_tip4pew
source leaprc.gaff
LIG = loadmol2 LIG.mol2
loadamberparams LIG.frcmod
check LIG
receptor = loadPDB rec_final.pdb
complex = combine {receptor LIG}
set default PBRadii mbondi2
saveAmberParm LIG LIG.top LIG.crd
receptor = loadPDB rec_final.pdb
saveAmberParm receptor rec.top rec.crd
saveAmberParm complex com.top com.crd
savepdb complex complex_gas.pdb
charge complex
addIons2 complex Na+ 0
solvateBox complex TIP4PBOX 10.0
saveAmberParm complex com_solvated.top com_solvated.crd
savepdb complex com_solvated.pdb
quit

Use either Cl- or Na+ based on the counter ions need to neutralise the system.

Use the following command:

tleap -s -f tleap.all

You will get the output com_solvated.top and com_solvated.crd from tleap.

6. Transform the Amber .top .crd to Gromacs compatible .to and .gro files:

Use the acpype

acpype -p com_solvated.top -x com_solvated.crd -b complex

it will create two files complex_GMX.top and complex_GMX.gro

7. Manipulate the topology file

Acpype can’t make a Gromacs compatible .top file with TIP4P water model hence we need to add some parts in the complex_GMX.top file to make it Gromacs compatible.

Add the following line in the atom_types section

[ atomtypes ]


HW_tip4pew HW_tip4pew 0.000 0.0000 A 0.00000e+00 0.00000e+00
OW_tip4pew OW_tip4pew 0.00 0.0000 A 3.16435e-01 6.80946e-01
 MW MW 0.00000 0.00000 A 0.00000e+00 0.00000e+00 ; 0.00 0.0000

and make sure you have the following part at the end of the topology file

[ moleculetype ]
 ; molname nrexcl
 NA+ 1

[ atoms ]
 ; id_ at type res nr residu name at name cg nr charge mass
 1 Na+ 1 NA+ NA+ 1 1 22.9898

[ moleculetype ]
; molname nrexcl ; TIP3P model
 WAT 2

[ atoms ]
; id at type res nr res name at name cg nr charge mass
 1 OW_tip4pew 1 SOL O 1 0 16.00000
 2 HW_tip4pew 1 SOL H1 1 0.52422 1.00800
 3 HW_tip4pew 1 SOL H2 1 0.52422 1.00800
 4 MW 1 SOL EPW 1 -1.04844 0.00000


#ifndef FLEXIBLE

[ settles ]
; i funct doh dhh
1 1 0.09572 0.15139


#else
[ bonds ]
; i j funct length force.c.
1 2 1 0.09572 502416.0 0.09572 502416.0
1 3 1 0.09572 502416.0 0.09572 502416.0

[ angles ]
; i j k funct angle force.c.
2 1 3 1 104.52 628.02 104.52 628.02

#endif


[ virtual_sites3 ]
; Vsite from funct a b
4 1 2 3 1 0.106676721 0.106676721


[ exclusions ]
1 2 3 4
2 1 3 4
3 1 2 4
4 1 2 3


; The position of the virtual site is computed as follows:
;
; O
;
; V
;
; H H
;
; Ewald tip4p:
; const = distance (OV) / [ cos (angle(VOH)) * distance (OH) ]
; 0.0125 nm / [ cos (52.26 deg) * 0.09572 nm ]
; then a = b = 0.5 * const = 0.106676721
;
; Vsite pos x4 = x1 + a*(x2-x1) + b*(x3-x1)


[ system ]
 complex

[ molecules ]
; Compound nmols
 complex 1 
 NA+ 8 
 WAT 9297

8. Now the first step is of energy minimisation:

Use the em.mdp file to grompp and generate a .tpr file

#em.mdp
; to test
; grompp -f em.mdp -c test_GMX.gro -p test_GMX.top -o em.tpr -v
; mdrun -v -deffnm em
 #Grompp minimization
integrator = steep
nsteps = 500
nstlist = 10
energy-grps = protein LIG
cutoff-scheme = verlet
vdw-type = cut-off
energy-grps = protein LIG
rvdw = 1.0
coulombtype = pme
rcoulomb = 1.0

grompp -f em.mdp -c complex_GMX.gro -p complex_GMX.top -o em

In order to minimise use the following command

mdrun -v -deffnm em

You should get em.gro , em.log , em.trr and em.edr files as output.

9. Put the restraint in the protein C-alpha and ligand heavy atoms during first step of equilibration:

#Generate index file
echo '13 & !a H*
17 | 3
q' | make_ndx -f em.gro

Please check which option is suitable for your protein

#Generate restraints

echo 18 | genrestr -f em.gro -n index.ndx -fc 500 500 500

After that open the complex_GMX.top file and add the following lines

#ifdef POSRE
#include "posre.itp"
#endif

Now it is time for grompp using the following eq.mdp file

~/gromacsbin462/bin/grompp -f eq.mdp -c em.gro -p complex_GMX.top -o eq

You should get eq.tpr as output. Now you can submit an equilibration run.

eq.mdp

; Input file
;
define = -DPOSRES
; integrator
integrator = md 
nsteps = 100000000000
dt = 0.002
;
; removing CM translation and rotation
comm_mode = Linear 
nstcomm = 1000
;
; output control 
nstlog = 500
nstenergy = 500
nstxout = 50000 
nstvout = 50000 
nstfout = 0
; group definition
nstxtcout = 500
xtc-grps = protein LIG
energy-grps = protein LIG
;
; neighbour searching
nstlist = 20
ns_type = grid
pbc = xyz
rlist = 1.2
;
; electrostatic
rcoulomb = 1.2
coulombtype = pme
fourierspacing = 0.16
;
; vdw
vdw-type = Cut-off 
rvdw = 1.2 
;
; constraints
constraints = all-bonds 
constraint-algorithm = lincs
lincs_iter = 4
;
; temperature
Tcoupl = v-rescale 
tc_grps = system 
tau_t = 0.1
ref_t = 298.0
;
; pression
;Pcoupl = no
Pcoupl = berendsen
Pcoupltype = isotropic
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0 
;
; initial velocities
gen_vel = yes 
gen_temp = 298.0 
gen_seed = -1

10. Running a non-restraint second NPT equilibration:

Take the last snapshot of the NPT restraint simulation and grompp again

~/gromacsbin462/bin/grompp -f eq.mdp -c confout.gro -p complex_GMX.top -o eq

Use the following eq.mdp
Note: there is no -DPOSRES in the define part

; Input file
;
define = 
; integrator
integrator = md 
nsteps = 100000000000
dt = 0.002
;
; removing CM translation and rotation
comm_mode = Linear 
nstcomm = 1000
;
; output control 
nstlog = 500
nstenergy = 500
nstxout = 50000 
nstvout = 50000 
nstfout = 0
; group definition
nstxtcout = 500
xtc-grps = protein LIG
energy-grps = protein LIG
;
; neighbour searching
nstlist = 20
ns_type = grid
pbc = xyz
rlist = 1.2
;
; electrostatic
rcoulomb = 1.2
coulombtype = pme
fourierspacing = 0.16
;
; vdw
vdw-type = Cut-off 
rvdw = 1.2 
;
; constraints
constraints = all-bonds 
constraint-algorithm = lincs
lincs_iter = 4
;
; temperature
Tcoupl = v-rescale 
tc_grps = system 
tau_t = 0.1
ref_t = 298.0
;
; pression
;Pcoupl = no
Pcoupl = berendsen
Pcoupltype = isotropic
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0 
;
; initial velocities
gen_vel = yes 
gen_temp = 298.0 
gen_seed = -1

11. Running a NVT production:

Combine the ener.edr file produced from NPT non-restraint equilibration run and plot the volume component

g_energy -f ener.part0001.edr -o vol.xvg

Take a point in the vol.xvg and dump a .gro file

trjconv -s plm.tpr -f traj.part0001.trr -dump 1800 -o conf1800.gro

Now use conf1800.gro file and grompp using the md.mdp file (NVT ensemble).

~/gromacsbin462/bin/grompp -f md.mdp -c conf1800.gro -p complex_GMX.top -o md

md.mdp

; Input file
;
define = 
; integrator
integrator = md 
nsteps = 100000000000
dt = 0.002
;
; removing CM translation and rotation
comm_mode = Linear 
nstcomm = 1000
;
; output control 
nstlog = 500
nstenergy = 500
nstxout = 50000 
nstvout = 50000 
nstfout = 0
; group definition
nstxtcout = 500
xtc-grps = protein Ion
energy-grps = protein Ion
;
; neighbour searching
nstlist = 20
ns_type = grid
pbc = xyz
rlist = 1.2
;
; electrostatic
rcoulomb = 1.2
coulombtype = pme
fourierspacing = 0.16
;
; vdw
vdw-type = Cut-off 
rvdw = 1.2 
;
; constraints
constraints = all-bonds 
constraint-algorithm = lincs
lincs_iter = 4
;
; temperature
Tcoupl = v-rescale 
tc_grps = system 
tau_t = 0.1
ref_t = 298.0
;
; pression
;Pcoupl = no
;Pcoupl = berendsen
;Pcoupltype = isotropic
;tau_p = 0.5
;compressibility = 4.5e-5
;ref_p = 1.0 
;
; initial velocities
gen_vel = yes 
gen_temp = 298.0 
gen_seed = -1

How to fix the PBC issue in Gromacs

gmx trjcat -f traj.part*.xtc -o trajall.xtc
echo 16 | gmx trjconv -s md.tpr -f trajall.xtc -pbc mol -o trajwhole.xtc -n index.ndx
echo 16 16 | gmx trjconv -s md.tpr -f trajwhole.xtc -pbc cluster -o trajclust.xtc -n index.ndx
echo 16 | gmx rmsf -f trajclust.xtc -ox av.pdb -s md.tpr -n index.ndx
echo 3 0 | gmx trjconv -s av.pdb -f trajclust.xtc -fit rot+trans -o trajfit_sparse20.xtc -dt 20

Choose your index number carefully.