University of Kurdistan, Sanandaj logo

Department of Chemistry, Faculty of Science, University of Kurdistan

MM/PB(GB)SA calculations

The molecular mechanics energies combined with the Poisson–Boltzmann or generalized Born and surface area continuum solvation (MM/PBSA and MM/GBSA) methods are popular approaches to estimate the free energy of the binding of small ligands to biological macromolecules. They are typically based on molecular dynamics simulations of the receptor-ligand complex and are intermediate in accuracy and computational effort between empirical scoring and strict alchemical perturbation methods. They have been applied to many systems with varying success [1]. You can use the following procedure to perform such calculations by AMBER.

1.   Set up the receptor file as described on the AMBER-MD page

2.   Dock your ligand into the receptor

To prepare the ligand-receptor  PDB file (complex.pdb), you can follow the following web page's instructions explaining the docking procedure using AutoDock Vina.

https://prof.uok.ac.ir/m.irani/index_files/Page512.htm

 

3.   Prepare the topology and coordinate files

3.1      Solvate the complex.pdb file

Run tleap as the following and the tleap.in-MD file (see below). It solvates the complex and creates prmtop and prmcrd files for MD simulations. Include the lines that are shown in boldface in the tleap.in-MD file only if you have disulfide bridges in your system. The first and second numbers in each line represent the residue number of the corresponding cysteine group. For more information on this, see section 2.4 on the AMBER-MD page.

tleap -s -f tleap.in-MD

 

A sample tleap.in-MD file

source leaprc.ff14SB

source leaprc.gaff

source leaprc.water.tip3p

addPath /home/Mehdi/docpbsa/Parms

loadAmberPrep lig.in

loadAmberParams lig.dat

x=loadpdb complex.pdb

bond x.3.SG x.64.SG

bond x.11.SG x.57.SG

bond x.48.SG x.86.SG

bond x.75.SG x.116.SG

bond x.79.SG x.103.SG

solvateOct x TIP3PBOX 10

saveamberparm x prmtop prmcrd

savepdb x pdbout

quit

 

3.2      Create topology and coordinate files of the complex

Run tleap as the following and the tleap.compl file (see below). It creates compl.prmtop and compl.prmcrd files for MM/PB(GB)SA calculations.

 

tleap -s -f tleap.compl

A sample tleap.compl file

source leaprc.ff14SB

source leaprc.gaff

source leaprc.water.tip3p

addPath /home/Mehdi/docpbsa/Parms

loadAmberPrep lig.in

loadAmberParams lig.dat

x=loadpdb complex.pdb

bond x.3.SG x.64.SG

bond x.11.SG x.57.SG

bond x.48.SG x.86.SG

bond x.75.SG x.116.SG

bond x.79.SG x.103.SG

saveamberparm x compl.prmtop compl.prmcrd

savepdb x compl.pdbout

quit

 

3.3      Prepare topology and coordinate files for the receptor

3.3.1      Remove the ligand coordinates from the complex.pdb file and create the receptor.pdb file.

sed '/LIG/d' complex.pdb > receptor.pdb

 

3.3.1.1     Run tleap as the following and the tleap.receptor file (see below). It creates receptor.prmtop and receptor.prmcrd files for MM/PB(GB)SA calculations.

 

tleap -s -f tleap.receptor

A sample tleap.receptor file

source leaprc.ff14SB

source leaprc.gaff

source leaprc.water.tip3p

addPath /home/Mehdi/docpbsa/Parms

x=loadpdb receptor.pdb

bond x.3.SG x.64.SG

bond x.11.SG x.57.SG

bond x.48.SG x.86.SG

bond x.75.SG x.116.SG

bond x.79.SG x.103.SG

saveamberparm x receptor.prmtop receptor.prmcrd

savepdb x receptor.pdbout

quit

 

3.4      Prepare topology and coordinate files for the ligand

3.4.1      Cut the ligand coordinates from the complex.pdb file and create lig.pdb file.

grep "LIG"   complex.pdb >lig.pdb

Run tleap as the following and the tleap.lig file (see below). It creates lig.prmtop and lig.prmcrd files for MM/PB(GB)SA calculations.

 

tleap -s -f tleap.lig

A sample tleap.lig file

source leaprc.ff14SB

source leaprc.gaff

source leaprc.water.tip3p

addPath /home/Mehdi/docpbsa/Parms

loadAmberPrep lig.in

loadAmberParams lig.dat

x=loadpdb lig.pdb

saveamberparm x lig.prmtop lig.prmcrd

savepdb x lig.pdbout

quit

 

3.4.2      solvate the receptor and ligand

Optionally, you can also solvate the receptor and the ligand files, perform MD simulations independently, and use the resulting trajectory files (instead of the complex) to extract the receptor and ligand coordinates.

 

tleap -s -f tleap.lig-md

tleap -s -f tleap.receptor-md

 

A sample tleap.receptor-md file

 

source leaprc.ff14SB

source leaprc.gaff

source leaprc.water.tip3p

x=loadpdb receptor.pdb

solvateOct x TIP3PBOX 10.0

saveamberparm x receptor.prmtop-md receptor.prmcrd-md

savepdb x receptor.pdbout-md

quit

 

A sample tleap.lig-md file

source leaprc.ff14SB

source leaprc.gaff

source leaprc.water.tip3p

addpath /home/Mehdi/docpbsa/Parms/

loadAmberPrep lig.in

loadAmberParams lig.dat

x=loadpdb lig.pdb

solvateOct x TIP3PBOX 10.0

saveamberparm x lig.prmtop-md lig.prmcrd-md

savepdb x lig.pdbout-md

quit

 

4.   Run MD simulations

4.1.1      MD simulations for the solvated complex

Run MD simulations with the following commands and the sander.in files listed in section 2.6 on the AMBER-MD page.

 

pmemd.cuda  -O -i sander.in1 -o sander.out1 -r mdrest1 -p prmtop -c prmcrd  -ref prmcrd

pmemd.cuda  -O -i sander.in2 -o sander.out2 -r mdrest2 -p prmtop -c mdrest1 -ref prmcrd

pmemd.cuda  -O -i sander.in3 -o sander.out3 -r mdrest3 -p prmtop -c mdrest2 -ref prmcrd

pmemd.cuda  -O -i sander.in4 -o sander.out4 -r mdrest4 -p prmtop -c mdrest3

pmemd.cuda  -O -i sander.in5 -o sander.out5 -r mdrest5 -p prmtop -c mdrest4 -x mdcrd5

 

4.2      MD simulations for the solvated receptor and ligand

If you did point ‎3.4.2 you need to execute the following commands to run MD simulations for the solvated receptor and ligand.

Commands for MD simulations for the solvated receptor:

 

pmemd.cuda -O -i sander.in1     -o sander.out1-receptor -p receptor.prmtop-md -r mdrest1-receptor -c receptor.prmcrd-md  -ref receptor.prmcrd-md

pmemd.cuda -O -i sander.in2     -o sander.out2-receptor -p receptor.prmtop-md -r mdrest2-receptor -c mdrest1-receptor    -ref receptor.prmcrd-md

pmemd.cuda -O -i sander.in3     -o sander.out3-receptor -p receptor.prmtop-md -r mdrest3-receptor -c mdrest2-receptor    -ref receptor.prmcrd-md

pmemd.cuda -O -i sander.in4     -o sander.out4-receptor -p receptor.prmtop-md -r mdrest4-receptor -c mdrest3-receptor    -ref receptor.prmcrd-md

pmemd.cuda -O -i sander.in5     -o sander.out5-receptor -p receptor.prmtop-md -r mdrest5-receptor -c mdrest4-receptor    -x mdcrd-receptor

 

 

It is recommended to use the sander module instead of pmemd.cuda to run MD simulations of the solvated ligand.

Commands for MD simulations for the solvated ligand.

mpirun -bind-to core -np 18 sander.MPI  -O -i sander.in1-lig -o sander.out1-lig  -p lig.prmtop-md  -r mdrest1-lig  -c lig.prmcrd-md   -ref lig.prmcrd-md

mpirun -bind-to core -np 18 sander.MPI  -O -i sander.in2-lig -o sander.out2-lig  -p lig.prmtop-md  -r mdrest2-lig  -c mdrest1-lig     -ref lig.prmcrd-md

mpirun -bind-to core -np 18 sander.MPI  -O -i sander.in3-lig -o sander.out3-lig  -p lig.prmtop-md  -r mdrest3-lig  -c mdrest2-lig     -ref lig.prmcrd-md

mpirun -bind-to core -np 18 sander.MPI  -O -i sander.in4-lig -o sander.out4-lig  -p lig.prmtop-md  -r mdrest4-lig  -c mdrest3-lig     -ref lig.prmcrd-md

mpirun -bind-to core -np 18 sander.MPI  -O -i sander.in5-lig -o sander.out5-lig  -p lig.prmtop-md  -r mdrest5-lig  -c mdrest4-lig     -x mdcrd-lig

 

5.   Run MM/PB(GB)SA calculations

Run MM/PB(GB)SA calculations with the MMPBSA.py.MPI module and the mmpbsa.in file in the following. Change the number shown in the boldface (18) to the number of CPU cores of your machine.

mpirun -np 18 MMPBSA.py.MPI -O -i mmpbsa.in -o MMPBSA.dat -sp prmtop -cp compl.prmtop -rp receptor.prmtop -lp lig.prmtop -y mdcrd5

If you did point ‎3.4.2 execute the following command instead.

mpirun -np 18 MMPBSA.py.MPI -O -i mmpbsa.in -o MMPBSA.dat -sp prmtop -srp receptor.prmtop-md -slp lig.prmtop-md -cp compl.prmtop -rp receptor.prmtop -lp lig.prmtop -y mdcrd5 -yr mdcrd-receptor -yl mdcrd-lig

 

A sample mmpbsa.in file

Input file for running PB and GB

&general

   verbose=1

   interval = 1

   startframe = 1

   entropy=1,

/

&gb

  igb=2, saltcon=0.100

/

&pb

  istrng=0.100,    inp=1,   radiopt=0

/

 

6.   References

[1]      S. Genheden, U. Ryde, The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities, Expert Opin. Drug Discov. 10 (2015) 449–461. https://doi.org/10.1517/17460441.2015.1032936.