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. Information on acpype can be accessed here.
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
If you encounter an error like Atomtype IP not found then open the topology file and put Na+ instead of IP.
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 POSRES #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 refcoord_scaling = com ; ; 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.
