MM/PBSA

MM/PBSA calculates the binding free energy as the difference between the free energies of the complex, receptor, and ligand in solution. The total free energy is decomposed into:

• Molecular Mechanics (E_MM): Bonded and non-bonded (van der Waals and electrostatic) terms from the force field
• Polar Solvation (G_PB): Electrostatic contribution to solvation, solved via the Poisson-Boltzmann equation
• Non-Polar Solvation (G_SA): Cavity formation and van der Waals interactions with solvent, estimated from the Solvent Accessible Surface Area (SASA)

This decomposition allows researchers to identify which energetic components drive binding, such as whether affinity is dominated by van der Waals packing or electrostatic complementarity.

MM/PBSA can be applied using two distinct protocols:

• Single-Trajectory Approach: Only the complex is simulated. The receptor and ligand snapshots are extracted from the complex trajectory. This assumes minimal conformational change upon binding and benefits from cancellation of intramolecular energy errors, yielding more stable results.
• Multi-Trajectory Approach: Separate simulations are run for the complex, receptor, and ligand. This accounts for induced fit and conformational reorganization but introduces higher statistical noise due to incomplete sampling of unbound states.

The single-trajectory protocol is the most widely used due to its computational efficiency and reduced variance, though it may underestimate the reorganization energy penalty.

MM/PBSA replaces explicit water molecules with continuum solvent models, dramatically reducing computational cost:

• Poisson-Boltzmann (PB) Model: Solves the linearized or non-linear PB equation on a grid to compute the electrostatic potential. Accounts for ionic strength via the Debye-Hückel screening parameter.
• Generalized Born (GB) Model: An analytical approximation to the PB equation that is faster but less accurate, particularly for buried charges and highly charged systems.
• Non-Polar Term: Typically calculated as G_SA = γ × SASA + b, where γ is the surface tension coefficient and b is a constant offset.

The choice of implicit solvent model and dielectric constants (commonly ε_in = 1-4 for solute, ε_out = 80 for water) critically affects the accuracy of the calculated binding free energies.

The configurational entropy change upon binding (ΔS) can be estimated using Normal Mode Analysis (NMA) or Quasi-Harmonic Analysis:

• Normal Mode Analysis: Diagonalizes the mass-weighted Hessian matrix to obtain vibrational frequencies, from which translational, rotational, and vibrational entropies are calculated using statistical mechanics formulas.
• Quasi-Harmonic Analysis: Approximates the entropy from the covariance matrix of atomic fluctuations in the MD trajectory.

Entropy calculations are computationally expensive and often the largest source of error in MM/PBSA. Many practitioners omit entropy or apply it only to a subset of snapshots. The -TΔS term can be significant for ligands with many rotatable bonds due to the loss of conformational freedom upon binding.

A powerful feature of MM/PBSA is the ability to decompose the total binding free energy into per-residue contributions:

• Identifies hot-spot residues that contribute most favorably to binding
• Reveals residues that oppose binding through electrostatic desolvation penalties
• Guides structure-based drug design by highlighting which protein-ligand interactions to optimize
• Can be extended to pairwise residue-residue decomposition to map interaction networks

This decomposition is achieved by calculating the interaction energy between the ligand and each individual receptor residue, providing a quantitative binding energy fingerprint that complements visual inspection of the binding pocket.

MM/PBSA occupies a middle ground between docking scores and rigorous alchemical methods:

• vs Docking: MM/PBSA includes solvation effects and conformational sampling from MD, providing better ranking power than static docking scores
• vs FEP/TI: Alchemical free energy perturbation (FEP) and thermodynamic integration (TI) are more theoretically rigorous but require 10-100x more computational resources due to simulating multiple intermediate lambda states
• vs LIE: The Linear Interaction Energy (LIE) method is faster but requires empirical scaling parameters calibrated to experimental data

MM/PBSA is particularly suited for lead optimization campaigns where relative binding free energies for 50-200 compounds must be estimated with moderate accuracy at reasonable computational cost.