I have docked a known inhibitor to my protein and the initial docking energy was ~ 5Kcal/mol. I ran 40ns simulation using AMBER and when I performed MMGBSA analysis on this complex, the binding energy between ligand and protein is coming ~470Kcal/mol. Throughout the simulation my inhibitor is staying pretty much inside the active site. I am unable to understand the reason for this high positive value of Binding Energy of the known inhibitor to the protein. I am giving the input file that I am using:
MMPBSA.py input file for running GB
&general
startframe=100, endframe=500, interval=5,
keep_files=0, netcdf=1,
/
&gb
igb=5, saltcon=0.1
/
&decomp
idecomp=4, dec_verbose=0, csv_format=0
Any help is gratefully appreciated. Thanks!
You need to use &pb and not &gb, an example would be
&pb
indi=4 ##Interior dielctric
exdi=80.0 ##Exterior dielectric.
1. Follow the AMBER tutorial to calculate entropy. The more snapshot you choose, the more accurate the result will be. This need huge memory and very long time. MMPBSA.MPI could make it faster.
2. Set internal dielectric constant for your system for binding energy calculation. This has nothing to do with entropy calculation. You may choose 1 to 4 according to your system. Read the following paper to be helpful to choose a proper indi parameter.
Tingjun Hou et al., “Assessing the Performance of the MM/PBSA and MM/GBSA Methods: I. The Accuracy of Binding Free Energy Calculations Based on Molecular Dynamics Simulations,” Journal of Chemical Information and Modeling 51, no. 1 (January 24, 2011): 69–82, doi:10.1021/ci100275a.
A couple of things:
1. It looks like you are only using a small part of your simulations (frames 100 to 500). You should use all of your frames and just increase your interval so you get a more diverse (and accurate) calculation (i.e. 80 frames over 40 ns is generally better than 80 frames over, say, 5 ns).
2. Check the output first to see where the source of the high positive free energy is coming from. This can give you an idea of where the error is from. It may be that there is 1 or 2 residues contributing very highly, or a bad parameter for the ligand.
3. If you calculate the entropy I would recommend using most (> 50%) of the frames. I have found this gives the most accurate results.
4. Also check out the original paper for the AMBER MMPBSA.py, 10.1021/ct300418h
this is a sample file. I have covered the entire trajectory in parts. How do i change the dielectric?
You need to use &pb and not &gb, an example would be
&pb
indi=4 ##Interior dielctric
exdi=80.0 ##Exterior dielectric.
Hi,
As you explained the inhibitor is pretty inside the binding pocket, there two possible ways to trace the error.
1. Instead of running many snapshots, try to pick 3-4 from a regions which shows stable backbone RMSD. If this gives negative BE, then do one more trial by extracting ~50 snapshots from a stable RMSD region.
2. Entropy play a major role which is significantly underestimated by MMPBSA for high computational cost. So try to use it for 5-10 snapshots as described above. If it works, then proceed further. Using many frames without tracing the error may kill ur computational as well as research time.
U can plot the BE and can see if there is any negative BE. And trace out which snapshot and frame number. Plot ur present data in grace and check the graph once.
Ref- http://www.sciencedirect.com/science/article/pii/S092777651400263X
Good luck.
- Badges
- Science topic