Protein Contacts and Binding Energy in MMPBSA Calculations

Dear GROMACS users,

I have calculated MMPBSA for four docked protein models, keeping the receptor the same across all models. I used the following command:
gmx_MMPBSA -O -i mmpbsa.in -cs complex.tpr -ct complex.xtc -ci index.ndx -cg 18 19 -cp topol.top -o FINAL_RESULTS_MMPBSA.dat -eo FINAL_RESULTS_MMPBSA.csv*
The mmpbsa.in file contains the following script:

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=2000,
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,
print_res=“within 4”
/

I obtained the binding energy for the four models (see the attached plot).
image
However, when I calculated the number of contacts between the protein pairs, I observed that the model with the most fluctuations in the contacts vs. time plot (screenshot attached) also had the most negative binding energy (model 2).
image

Is it possible for a protein pair to exhibit such fluctuations in contacts yet have a highly negative binding energy, even when the interactions are primarily non-covalent (e.g., hydrogen bonds, hydrophobic interactions, van der Waals forces)? Or could there be an issue with how I am calculating MMPBSA?

Any suggestions or insights would be greatly appreciated.

It seems like the proteins were separated, which led to the jumping of the contact curve. This is usually due to the periodic boundary conditions (PBC), and reflected in the trajectory file. You can check what happened in the trajectory during the fluctuations.
Furthermore, the MMPBSA result will be untrustworthy if the trajectory is unstable. Therefore, only the stable part at the end of the trajectory is typically used to calculate MMPBSA.