HIgh propagated SD from gmx_mmpbsa

[Axelle Thébault]

Hello everyone,

I am trying to use gmx_mmpbsa on protein-protein complexes I have simulated using GROMACS 2025.1. 

First of all, thank you for this wonderful tool, this is exactly what I was looking for. However, I'm facing an issue regarding the output of my energy calculation. You'll find my FINAL_RESULTS_MMPBSA.dat and FINAL_RESULTS_MMPBSA.csv files attached below. 

I am working on an antibody-antigen system: my antibody (IgG) is represented by 2 chains (corresponding to one of the antibody's heavy and light 'arm') that are bound to each other, and the third chain of my system is the antigen. I conducted my simulation for 30 ns (I can also provide you the .mdp files if needed) on the system that underwent docking with HADDOCK before the MD. My MD with GROMACS performed well (I used CHARMM27 force field and TIP3P water type), my minimization (steepest descent with Fmax < 1000) and equilibration steps (NVT brought the system to 300K on average and NPT made the pressure average 1 bar) were successful and I see no issue when i visualize my protein-protein complex with PyMOL. 

But when I open my .dat file after using gmx_mmpbsa, I see that my SD(prop) is extremely high, at least twice bigger than the average dG itself (see the .dat file). There clearly is a problem here, as a dG equal to -30 more or less 80 can mean very different things, the standard deviation is bad. I used the .tpr file of the complex after the MD and the fit trajectory with PBC removed. 

It seems that the dEEL and dEGB are very high and this might be the reason why the final results are so bad and unusable, but I wonder if this issue comes from my GROMACS simulation (maybe unstable system or something else?) or from my use of your tool (did I chose bad parameters given my system?). 

I also provide you my mmpbsa.in file and my RMSD file from my MD simulation, that I used to select the frames for the energy calculation (I chose the ones where the RMSD is the more stable, so between 7ns and 15 ns). I also tried modifying the salt concentration in the input file, but it did not change the gmx_mmpbsa results much.

I am open to any hypothesis regarding this issue that I have been struggling with for quite some time now. I also performed some MD followed by gmx_mmpbsa calculation for other similar systems, and all the results share the same problems (very high SD(prop)).

[Axelle Thébault]

Sorry for the energy minimization in GROMACS, the Fmax was 100 kJ.mol-1, not 1000

[marioe911116]

That's usually the case for MMPB(GB)SA calculations, unless you use constraints in the simulation. SEM is relatively lower, though (0.22), which is usually the value that is reported... In any case, I will make sure there are no instabilities in the system, or "jump", the system has been fitted, and it should be good to go...

[Axelle Thébault]

Thanks for you quick reply. 

I thought my values were too high compared to the ones I have seen in other forum topics or even in the example output file on the gmx_mmpbsa website. How to explain such important values for the dEEL and dEGB ? I thought maybe there was a parameter I did not set well for the GB calculation.

So for you the values are good as they are and I should not consider the SD(prop) but rather the SEM(prop) ? Also, what constraints are you talking about, do you mean during the MD or the MMPBSA calculation ?

My system does not jump or break, I made sure of that by using trjconv -pbc whole, -pbc nojump, -pbc mol and -fit rot+trans on my MD trajectories.

Thanks again for your insights on the matter, I am new to these tools.

[marioe911116]

You can report SEM values... I meant restraints, sorry... to prevent the system from fluctuating too much...