Introduction¶

The MM/PB(GB)SA method can be used to calculate the binding free energies of noncovalently bound complexes.

The binding free energy of a complex can be estimated as follows:

∆𝐺_{𝑏𝑖𝑛𝑑} = 〈𝐺_𝐶𝑂𝑀〉−〈𝐺_𝑅𝐸𝐶〉−〈𝐺_𝐿𝐼𝐺〉

(1)

where each term on the right-hand side is given by:

〈𝐺_𝑥〉 = 〈𝐸_𝑀𝑀〉 + 〈𝐺_𝑠𝑜𝑙〉 − 〈𝑇𝑆〉

(2)

In turn, ∆𝐺_{𝑏𝑖𝑛𝑑} can also be represented as:

∆𝐺_{𝑏𝑖𝑛𝑑} = ∆𝐻 − 𝑇∆𝑆

(3)

where ∆𝐻 corresponds to the enthalpy of binding and −𝑇∆𝑆 to the conformational entropy after ligand binding. When the entropic term is omitted, the computed value is an enthalpy-like effective binding estimate. It can be useful for a defined relative-comparison protocol, but omitting entropy is an approximation whose adequacy depends on the systems, sampling, and scientific question; it is not generally sufficient by itself for relative affinity claims.

The enthalpy, ∆𝐻, can be decomposed into different terms:

∆𝐻 = ∆𝐸_𝑀𝑀 + ∆𝐺_𝑠𝑜𝑙

(4)

where:

∆𝐸_𝑀𝑀 = ∆𝐸_{𝑏𝑜𝑛𝑑} + ∆𝐸_{𝑎𝑛𝑔𝑙𝑒} + ∆𝐸_{𝑑𝑖ℎ𝑒𝑑} + ∆𝐸_𝑣𝑑𝑊 + ∆𝐸_𝑒𝑙𝑒 + ∆𝐸_{1-4 VDW} + ∆𝐸_{1-4 EEL}

(5)

The gas-phase contributions are calculated by sander within AmberTools according to the force field and method. The ordinary terms above are supplemented by UB, IMP, and CMAP for applicable CHARMM calculations and by ESCF for QM/MMGBSA. These terms are included when the output parser reports them. In ST, component differences for matching topologies can cancel; that cancellation does not redefine the component totals or apply automatically to MT.

The ∆𝐺_𝑠𝑜𝑙 is given by:

∆𝐺_𝑠𝑜𝑙 = ∆𝐺_𝑝𝑜𝑙 + ∆𝐺_{𝑛𝑜𝑛−𝑝𝑜𝑙} = ∆𝐺_{𝑃𝐵/𝐺𝐵} + ∆𝐺_{𝑛𝑜𝑛−𝑝𝑜𝑙}

(6)

where:

∆𝐺_{𝑛𝑜𝑛−𝑝𝑜𝑙𝑎𝑟} = 𝑁𝑃_{𝑇𝐸𝑁𝑆𝐼𝑂𝑁} ∗ ∆𝑆𝐴𝑆𝐴 + 𝑁𝑃_{𝑂𝐹𝐹𝑆𝐸𝑇}

(7)

or,

∆𝐺_{𝑛𝑜𝑛−𝑝𝑜𝑙} = ∆𝐺_{𝑑𝑖𝑠𝑝} + ∆𝐺_{𝑐𝑎𝑣𝑖𝑡𝑦} = ∆𝐺_{𝑑𝑖𝑠𝑝} + (𝐶𝐴𝑉𝐼𝑇𝑌_{𝑇𝐸𝑁𝑆𝐼𝑂𝑁} ∗ ∆𝑆𝐴𝑆𝐴 + 𝐶𝐴𝑉𝐼𝑇𝑌_{𝑂𝐹𝐹𝑆𝐸𝑇})

(8)

In the above equations, ∆𝐸_𝑀𝑀 corresponds to the molecular mechanical energy changes in the gas phase. ∆𝐸_𝑀𝑀 includes ∆𝐸_{𝑏𝑜𝑛𝑑𝑒𝑑}, also known as internal energy, and ∆𝐸_{𝑛𝑜𝑛𝑏𝑜𝑛𝑑𝑒𝑑}, corresponding to the van der Waals and electrostatic contributions. The solvation energy is determined differently depending on the method employed. In the 3D-RISM model, both the polar and nonpolar components of the solvation energy are calculated. However, the PB and GB models estimate only the polar component of the solvation energy. The nonpolar component is usually assumed to be proportional to the molecule's total solvent-accessible surface area (SASA), with a proportionality constant derived from experimental solvation energies of small nonpolar molecules (Eq. 7). Alternatively, a modern approach that separates nonpolar solvation free energies into cavity and dispersion terms can be used. In this approach, SASA is used to correlate the cavity term only, while a surface-integration method is employed to compute the dispersion term (Eq. 8).

Furthermore, the entropic component can be estimated with normal-mode analysis (NMODE). NMODE is Hessian-based: the energy is minimized and a mass-weighted Hessian is diagonalized around the minimized structure to obtain vibrational modes and frequencies. It is therefore distinct from quasi-harmonic (QH) analysis, which estimates fluctuations from a coordinate covariance matrix over a sampled trajectory. NMODE can be computationally expensive, although truncated systems can reduce the cost. New QH calculations are not supported in 1.7.0; historical QH results remain readable for compatibility only.

Interaction Entropy (IE) estimates an entropic contribution from the fluctuation of the interaction energy along an MD trajectory and has low additional post-processing cost. Its numerical behavior depends strongly on the distribution, fluctuations, and convergence of the sampled interaction energies; it is not universally superior to NMODE and should be checked with block or cumulative convergence diagnostics before interpretation. Multiple-trajectory IE/C2 use is experimental in this release because the bound and unbound trajectories are independently sampled. See the GROMACS normal-mode reference and Ekberg and Ryde (2021) for methodological context.

Typically, MM/PB(GB)SA calculations use one of two approaches: the single-trajectory protocol (STP) or the multiple-trajectory protocol (MTP). In STP, both the receptor and ligand trajectories are extracted from the complex trajectory. This approach is valid when the bound and unbound states of the receptor and ligand are similar. It is computationally less expensive than the MTP approach since only a simulation of the complex is required. Additionally, the potential internal terms (e.g., bonds, angles, and dihedrals) cancel exactly in STP since these terms are the same in both bound and unbound states. On the other hand, the MTP is a more realistic approach because it considers separate trajectories for the complex, receptor, and ligand. However, subtracting energies from independently sampled conformations can introduce substantial uncertainty. In practice, the system must be studied carefully to select the appropriate approach.

Literature¶

Further information can be found in Amber manual:

• MMPBSA.py
• The Generalized Born/Surface Area Model
• PBSA
• Reference Interaction Site Model
• Generalized Born (GB) for QM/MM calculations

and the foundational papers:

• Srinivasan J. et al., 1998
• Kollman P. A. et al., 2000
• Gohlke H., Case D. A. 2004

as well as some reviews and expert opinions:

• Genheden S., Ryde U. 2015
• Wang et. al., 2018
• Wang et. al., 2019
• Tuccinardi, 2021