ProteinIQ
gmx_MMPBSA icon

gmx_MMPBSA

Calculate binding free energies with MM/PBSA and MM/GBSA

What is gmx_MMPBSA?

gmx_MMPBSA calculates binding free energies from molecular dynamics trajectories using MM/PBSA (Molecular Mechanics/Poisson-Boltzmann Surface Area) and MM/GBSA (Molecular Mechanics/Generalized Born Surface Area) methods. Developed by Valdés-Tresanco and colleagues, it bridges GROMACS simulations with AMBER's MMPBSA.py engine, enabling end-state free energy calculations for protein-ligand, protein-protein, and protein-nucleic acid complexes.

These methods estimate binding affinity by calculating the free energy difference between bound and unbound states. The calculation combines molecular mechanics energies (van der Waals, electrostatic) with implicit solvation models, offering a balance between the speed of empirical scoring functions and the accuracy of rigorous alchemical methods.

How MM/PBSA and MM/GBSA work

Binding free energy (ΔGbind\Delta G_{\text{bind}}) is computed from snapshots extracted from an MD trajectory:

ΔGbind=Gcomplex−Greceptor−Gligand\Delta G_{\text{bind}} = G_{\text{complex}} - G_{\text{receptor}} - G_{\text{ligand}}

Each term comprises several energy components:

ComponentDescription
ΔEvdW\Delta E_{\text{vdW}}Van der Waals interactions (Lennard-Jones potential)
ΔEelec\Delta E_{\text{elec}}Coulombic electrostatic interactions
ΔGpolar\Delta G_{\text{polar}}Polar solvation energy (PB or GB model)
ΔGnonpolar\Delta G_{\text{nonpolar}}Nonpolar solvation, proportional to solvent-accessible surface area
−TΔS-T\Delta SConformational entropy (optional, computationally expensive)

MM/PBSA solves the Poisson-Boltzmann equation for electrostatic solvation, providing higher accuracy but at greater computational cost. MM/GBSA uses analytical Generalized Born approximations, running faster while maintaining reasonable accuracy for ranking compounds.

Generalized Born models

ModelDescription
GB-HCT (igb=1)Hawkins-Cramer-Truhlar model, fast but less accurate
GB-OBC1 (igb=2)Onufriev-Bashford-Case model, variant 1
GB-OBC2 (igb=5)OBC variant 2, recommended for most applications
GB-Neck (igb=7)Improved treatment of interstitial regions
GB-Neck2 (igb=8)Further refinements to neck region calculations

The OBC2 model (igb=5) provides a good balance of accuracy and speed for protein-ligand systems.

How to use gmx_MMPBSA online

ProteinIQ provides cloud-based access to gmx_MMPBSA, eliminating the need to configure AMBER, GROMACS, or compute environments locally.

Inputs

InputDescription
MD Trajectory fileGROMACS trajectory (.xtc, .trr) or multi-model PDB. Contains atomic coordinates across simulation frames.
Topology fileGROMACS topology (.tpr or .top). Defines force field parameters and system topology.
Index fileOptional GROMACS index file (.ndx). Specifies receptor and ligand atom groups. Auto-detected if omitted.
Input parameter fileOptional gmx_MMPBSA input file. Overrides settings with custom parameters.

Settings

Calculation settings

SettingDescription
Calculation methodMM/GBSA (faster, suitable for screening), MM/PBSA (more accurate), or Both for comparison
Energy decompositionWhen enabled, calculates per-residue energy contributions to identify key binding site residues

Trajectory settings

SettingDescription
Start frameFirst frame to analyze (1-indexed)
End frameLast frame to analyze. Use -1 to include all frames
Frame intervalAnalyze every nth frame. Set to 2 or higher to reduce computation for long trajectories

Solvation settings

SettingDescription
Salt concentrationIonic strength in molar (default 0.15 M, physiological conditions)
GB modelGeneralized Born variant for MM/GBSA. GB-OBC2 (igb=5) recommended

Output

Results include total binding free energy with component breakdown:

ColumnDescription
DELTA TOTALNet binding free energy (kcal/mol). More negative indicates stronger binding.
VDWAALSVan der Waals contribution
EELElectrostatic energy (gas phase)
EGB or EPBPolar solvation (GB or PB method)
ESURF or ENPOLARNonpolar solvation (surface area term)

When energy decomposition is enabled, per-residue contributions identify which amino acids drive binding.

Interpreting binding energies

ΔGbind\Delta G_{\text{bind}} (kcal/mol)Interpretation
< −15Very strong binding
−10 to −15Strong binding (nanomolar range)
−5 to −10Moderate binding
> −5Weak binding

These values should be interpreted as relative rankings rather than absolute affinities. MM/PBSA and MM/GBSA excel at comparing similar ligands binding to the same target, but absolute values carry significant uncertainty.

Applications

  • Lead optimization: Compare binding energies across congeneric compound series to guide medicinal chemistry
  • Binding hotspot identification: Energy decomposition reveals which residues contribute most to binding, informing mutagenesis studies
  • Selectivity analysis: Compare binding to on-target versus off-target proteins
  • Protein-protein interface analysis: Identify key residues stabilizing macromolecular complexes

Limitations

The single-trajectory approach used here assumes identical conformations for complex, receptor, and ligand. This reduces noise but misses conformational reorganization effects. Results depend heavily on the quality and length of the MD simulation—short simulations with poor sampling produce unreliable energies.

Entropy contributions are not calculated by default due to computational cost and large uncertainties. For charged ligands, electrostatic and polar solvation terms often exhibit significant cancellation, amplifying small errors. Results should guide compound prioritization rather than predict absolute binding constants.

Related tools

  • OpenMM: Run molecular dynamics simulations to generate input trajectories
  • AutoDock Vina: Fast docking for initial pose generation before MD refinement
  • DiffDock: Deep learning docking that may complement physics-based approaches
  • HADDOCK3: Information-driven docking for protein-protein complexes