Discrepancy Between MMPBSA Binding Energy and Experimental Value

GROMACS version:2020.4
Dear Users

I have carried out a 1 μs MD simulation of a protein–peptide complex, which remained stable throughout the trajectory. Following this, I computed the binding energy using MMPBSA with the following parameters:

Sample input file for decomposition analysis
This input file is meant to show only that gmx_MMPBSA works. Althought,
we tried to used the input files as recommended in the Amber manual,
some parameters have been changed to perform more expensive calculations
in a reasonable amount of time. Feel free to change the parameters
according to what is better for your system.

&general
sys_name=“Decomposition”,
startframe=0,
endframe=2001,
forcefields=“CHARMM36”
/
&gb
igb=5, saltcon=0.150,
/
#make sure to include at least one residue from both the receptor
##and ligand in the print_res mask of the &decomp section.
##this requirement is automatically fulfilled when using the within keyword.
##Re: [AMBER] Is there any problem with this input file? from Bill Miller III on 2013-08-05 (Amber Archive Aug 2013)
&decomp
idecomp=2, dec_verbose=3
/

From this, I obtained the following binding energy result:

Δ (Complex − Receptor − Ligand): Average ΔTOTAL = −64 kcal/mol

However, the experimentally reported value is approximately −10 kcal/mol.

Am I making a mistake in my setup or analysis? Any advice or guidance would be greatly appreciated.

Thanks

Hi @Devid,

I am not an MM/PBSA expert so I cannot comment on the “correctness” of your input parameters for your system. But keep in mind that MMPBSA is a rather “rough” method to calculate binding energy. It is fast but not the most reliable one if you want to compare to experimental results. See this review for further discussion on the method. As far as I remember it becomes more reliable for ΔΔG calculations. For more accurate results you should probably switch to other methods, e.g., FEP.

Hope I could help you,

Marius

Hii @Marius

Thank you for your reply

I have gone through a few online sources and found that we need to specify the free energy term (λ, etc.) in the md.mdp file during the simulations. However, I did not include these parameters and have already simulated for 1000 ns. As per my understanding, we cannot perform FEP without this information. Is there any alternative way to calculate the free energy from the existing MD simulation, or do I need to perform the simulation again with the proper settings?

Additionally, according to the literature, the reported values for the same system are approximately ΔG_coul = −18 kcal/mol, ΔG_np = −30 kcal/mol, and ΔG_solv = +3 kcal/mol, while my calculated values are significantly higher (screenshot attached). Could this discrepancy be due to an error in my analysis, or am I misunderstanding the interpretation of these results?