Drug Mmpbsa Gbsa

drug-mmpbsa-gbsa (MM-GBSA / MM-PBSA)

Goal

To estimate relative binding free energies from MD trajectories using the single-trajectory MM-GBSA / MM-PBSA approach. For each trajectory frame, the method strips explicit solvent, evaluates potential energies of the complex, receptor, and ligand subsystems in implicit solvent, and computes:

dG = E_complex - E_receptor - E_ligand

The entropy term (-TdS) is omitted, which is standard practice when the goal is relative ranking rather than absolute binding affinity. MM-GBSA / MM-PBSA is most useful for re-ranking docked poses after MD refinement, providing an orthogonal signal to docking scores and geometric stability metrics.

Choosing a backend

Both paths give comparable GB rankings for typical drug-protein systems, but the absolute dG numbers will differ across backends because they use different GB models, radius sets, and surface-area treatments. Don't compare numbers across the two scripts.

Instructions

1. Basic usage (protein + ligand, short HTVS-style MD)

For the 1-5 ns production runs typical in the HTVS workflow:

${CLAUDE_SKILL_DIR}/../../venv/run cpu+openmm python ${CLAUDE_SKILL_DIR}/scripts/compute_mmgbsa.py \
  --topology md/system/complex_solvated.pdb \
  --trajectory md/run/production.dcd \
  --ligand_sdf md/ligand.sdf \
  --ligand_resname UNL \
  --skip_ns 0.5 \
  --stride 5 \
  --output_dir md/mmgbsa/

Key parameters:

• --topology: solvated complex PDB from drug-complex-system-builder
• --trajectory: production DCD from drug-protein-ligand-md
• --ligand_sdf: ligand SDF with correct bond orders (for OpenFF parameterization and AM1-BCC charges)
• --charge_method (compute_mmgbsa.py): ligand partial charges, am1bcc (default; needs AmberTools' sqm on PATH, not in any uv environment) or mmff94 / gasteiger (RDKit; cruder, for screening and tests). Use the same method the complex was built with.
• --ligand_resname: residue name of the ligand in the PDB topology (default: UNL)
• --skip_ns: skip the first N nanoseconds of trajectory as equilibration (default: 0.5). For longer production runs, increase proportionally: roughly 10% of production length is a reasonable rule of thumb, capped at ~5 ns for 50+ ns runs.
• --stride: evaluate every Nth frame after skipping (default: 5)
• --solute_dielectric: interior dielectric of the solute (default: 1.0). See section 3 below on when to raise this.

2. With a cofactor (e.g., NADPH) and longer MD

When the receptor has a bound cofactor that should be part of the receptor subsystem, and the MD production was long enough (>=10 ns) to justify a longer equilibration skip:

${CLAUDE_SKILL_DIR}/../../venv/run cpu+openmm python ${CLAUDE_SKILL_DIR}/scripts/compute_mmgbsa.py \
  --topology md/system/complex_solvated.pdb \
  --trajectory md/run/production.dcd \
  --ligand_sdf md/ligand.sdf \
  --cofactor_sdf md/nadph.sdf \
  --ligand_resname UNL \
  --cofactor_resname NDP \
  --skip_ns 2.0 \
  --output_dir md/mmgbsa/

3. Tuning the solute dielectric

The interior dielectric of the solute (--solute_dielectric, OpenMM's soluteDielectric parameter) is a major knob for MM-GBSA rankings and is often the difference between sensible and nonsense results. The default of 1.0 treats the solute interior as electronically non-polarizable; the explicit partial charges and atom geometries are still there, but no implicit electronic polarization is added on top. This is the most physically pure choice and works well for nonpolar binding sites with few charged residues. For polar or highly charged pockets, raise it:

Higher dielectrics damp electrostatic contributions and generally improve ranking agreement with experiment for charged systems, at the cost of losing sensitivity to directional electrostatic interactions. If you are unsure, run the rescoring at 1.0 and 2.0 on a small subset with known rank-order and pick the value that tracks better. The chosen value is recorded in mmgbsa_summary.json for provenance.

${CLAUDE_SKILL_DIR}/../../venv/run cpu+openmm python ${CLAUDE_SKILL_DIR}/scripts/compute_mmgbsa.py \
  --topology md/system/complex_solvated.pdb \
  --trajectory md/run/production.dcd \
  --ligand_sdf md/ligand.sdf \
  --ligand_resname UNL \
  --solute_dielectric 2.0 \
  --output_dir md/mmgbsa/

4. Custom selections

For non-standard residue naming or multi-chain receptors:

${CLAUDE_SKILL_DIR}/../../venv/run cpu+openmm python ${CLAUDE_SKILL_DIR}/scripts/compute_mmgbsa.py \
  --topology md/system/complex_solvated.pdb \
  --trajectory md/run/production.dcd \
  --ligand_sdf md/ligand.sdf \
  --ligand_sel "resname UNK and not name H*" \
  --output_dir md/mmgbsa/

5. Interpret results

The script outputs:

• mmgbsa_frames.csv: per-frame energies (E_complex, E_receptor, E_ligand, dG)
• mmgbsa_summary.json: mean, std, SEM of dG, and the dielectric/frame parameters used

{
    "dG_mean_kcal_mol": -42.5,
    "dG_std_kcal_mol": 8.3,
    "dG_sem_kcal_mol": 1.2,
    "n_frames_evaluated": 50,
    "solute_dielectric": 1.0,
    "solvent_dielectric": 78.5
}

More negative dG indicates stronger predicted binding. Only relative ranking within a congeneric series is interpretable. Absolute dG values from single-trajectory MM-GBSA without entropy corrections are not binding free energies in any physical sense, and specific numbers should not be reported as "predicted affinities" or compared across chemically distinct scaffolds. When ranking compounds:

• Group comparisons by formal charge state (see Constraints below).
• Compare mean dG values, but also inspect the per-frame trace in mmgbsa_frames.csv for stability. A compound with a noisy or drifting dG signal is less trustworthy than one with a tight distribution, even if the mean looks favorable.
• Use MM-GBSA as an additional signal alongside geometric stability (RMSD, contact persistence) from drug-trajectory-analysis, not as a standalone ranking.
• Compare replicate means, not within-trajectory SEMs. The honest uncertainty for a compound is the spread across independent MD replicates, not the SEM from any single one. See examples/README.md for a real CDK2 case study that demonstrates this: one compound in the example shows a 35 kcal/mol swing across three MD replicates while each individual replicate reports a tight single-run SEM. Trusting any single replicate on that compound would be deeply misleading.

6. AmberTools path: MM-PBSA / MM-GBSA via MMPBSA.py

When you need Poisson-Boltzmann (not just GB), or when you want a method that matches what reviewers expect from the MM-PBSA literature, use compute_mmpbsa.py. It builds the same dry-complex / receptor / ligand subsystems as the OpenMM path, then converts the systems to Amber prmtops via ParmEd, converts the trajectory to NetCDF via cpptraj, and hands everything to AmberTools MMPBSA.py.

${CLAUDE_SKILL_DIR}/../../venv/run cpu+openmm python ${CLAUDE_SKILL_DIR}/scripts/compute_mmpbsa.py \
  --topology md/system/complex_solvated.pdb \
  --trajectory md/run/production.dcd \
  --ligand_sdf md/ligand.sdf \
  --ligand_resname UNL \
  --skip_ns 0.5 \
  --stride 5 \
  --method both \
  --output_dir md/mmpbsa/

Key parameters:

• --method: gb (only Generalized Born), pb (only Poisson-Boltzmann), or both (default). PB is much slower than GB (~5-10x for typical systems); use gb for HTVS-style throughput and both only when you need the cross-check.

• --gb_model: igb value passed to MMPBSA.py. The script sets the prmtop's GB radius set to match the chosen model via ParmEd.

  igb	Model	Paired radius set	Notes
  1	HCT (Hawkins / Cramer / Truhlar)	mbondi
  2	OBC1 (Onufriev / Bashford / Case)	mbondi2
  5	OBC2	mbondi2	Default; most widely used for protein-ligand MM-GBSA
  7	GBn / GB-Neck	mbondi	Mongan et al. 2007 (note: the Amber manual recommends mbondi here; some downstream tools and PB-focused docs recommend bondi instead - we follow the Mongan 2007 parameterization)
  8	GBn2 / GB-Neck2	mbondi3	Nguyen et al. 2013

  Note that the OpenMM compute_mmgbsa.py path uses GBn2 (Amber-equivalent igb=8 with mbondi3-like radii), while this AmberTools path defaults to OBC2 (igb=5 with mbondi2). They are different GB models. Even on a GB-only comparison, you should expect the absolute dG to differ across the two scripts; this is the main reason for the gap, alongside the LCPO / no-SA difference noted in the comparison table.

• --pb_int_diel (default 1.0), --pb_ext_diel (default 80.0): PB dielectric constants. Same dielectric-tuning logic applies as for GB (raise --pb_int_diel to 2-4 for polar / charged pockets, see section 3). For PB specifically, raising the interior dielectric above ~2 is an empirical fitting choice rather than a physical correction: ε_int ~ 2 is the upper end of what is physically defensible (modeling electronic polarization only); larger values are commonly used in the literature but they are tuning parameters, not first-principles. Note also that the OpenMM path defaults to a solvent dielectric of 78.5 (the standard 25 °C value), while the AmberTools path defaults to AmberTools' own 80.0 - one more reason absolute numbers are not comparable across backends.

• --salt_conc (default 0.0 M): salt concentration used by both GB (saltcon) and PB (istrng). Set to 0.15 to match physiological ionic strength.

The script writes:

• mmpbsa.in: the auto-generated MMPBSA.py input file (kept for provenance).
• complex.prmtop, receptor.prmtop, ligand.prmtop, *.inpcrd: ParmEd-built Amber topologies and coords. Atom ordering is preserved against the source PDB so the trajectory aligns without any tleap-side reordering.
• trajectory_dry.nc: the stripped, NetCDF-format trajectory MMPBSA.py consumed.
• FINAL_RESULTS_MMPBSA.dat: the verbatim MMPBSA.py results file (per-method breakdown of ELE, VDW, EGB/EPB, ESURF/ECAVITY, etc.).
• mmpbsa_summary.json: parsed dG mean / std / SEM per method, plus all provenance.

{
    "method": "both",
    "results": {
        "GB": {"dG_mean_kcal_mol": -51.13, "dG_std_kcal_mol": 2.36, "dG_sem_kcal_mol": 1.36},
        "PB": {"dG_mean_kcal_mol":  -0.39, "dG_std_kcal_mol": 2.85, "dG_sem_kcal_mol": 1.64}
    },
    "gb_model": 5,
    "pb_int_diel": 1.0,
    "pb_ext_diel": 80.0,
    ...
}

PB and GB will often disagree by 30+ kcal/mol in absolute dG. This is normal: PB and GB make different approximations for the polar solvation free energy, and the nonpolar / cavity terms also differ. Within a single method, ranking is what matters; across methods, large absolute differences are diagnostic of how sensitive your system is to the implicit-solvent approximation, not a sign that one method is "wrong". When PB and GB give wildly different rank orders on the same set of compounds, neither is reliable, and you are probably out of the regime where end-point free-energy methods work, and you should consider FEP / ABFE instead.

The same caveats from sections 3 and 5 (replicates beat single-run SEM, charge state matters, electrostatic environment mismatch) apply equally to the AmberTools path.