Binding Affinity Prediction

Overview

This skill predicts protein-ligand binding affinity from docked poses — converting structural information into estimated ΔG (kcal/mol), pKd, and Kd (nM). It complements the molecular-docking skill's interaction analysis (score.py) which counts contacts but does NOT predict binding strength in energy units.

Key capabilities:

• Empirical scoring: descriptor + contact-based affinity prediction using RDKit and BioPython
• Ligand energy inspection: explicitly selected RDKit heuristic; it does not calculate receptor-dependent binding energy
• Full MM/GBSA: use a separately validated external method; the bundled script does not implement it
• Consensus scoring: combine multiple scoring methods with rank-based normalization
• Batch virtual screening: efficiently score large compound libraries

Validation Warning

All predictions from this skill are computational estimates, NOT experimentally validated measurements.

• Empirical scoring (predict.py): typical error is 1-2 log units pKd (~10-100x in Kd)
• The bundled ligand heuristic cannot establish binding affinity or receptor-dependent ranking
• Use for prioritizing compounds for experimental testing, not for making clinical claims

The scripts include uncertainty ranges and confidence flags to help calibrate expectations.

Output Integrity

Script outputs are RAW computational estimates. The agent MUST NOT:

• Apply "calibration" or scaling to raw pKd/Kd values
• Adjust values to match known experimental data
• Add fields like "calibrated_pKd" not produced by the script
• Present approximate methods (simplified MM/GBSA) as full implementations

The raw output IS the prediction. Report it exactly as produced.

Successful script outputs are logged to _script_manifest.jsonl. A failed rescoring run may retain a JSON artifact with explicit failed records and exits nonzero; it does not write a success manifest entry. Verify the exit status and record status before using a score.

When to Use This Skill

Use the binding-affinity skill when you need to:

• Predict how tightly a ligand binds to a protein (after docking)
• Rank docked poses by estimated binding affinity
• Inspect ligand conformer energies with the explicit ligand heuristic; obtain MM/GBSA from a validated external workflow when required
• Screen compound libraries for binding potential
• Combine multiple scoring methods into a consensus ranking

Trigger phrases: "predict binding affinity", "estimate Kd", "score binding strength", "rescore with MM/GBSA", "rank compounds by affinity", "virtual screening"

Do NOT use this skill for:

• Finding binding pockets (use pocket-detection)
• Generating docked poses (use molecular-docking)
• Predicting ADMET properties (use admet-prediction)
• Free energy perturbation (requires specialized MD simulations)

Related Skills

• molecular-docking: Generate docked poses first. This skill scores them.
• pocket-detection: Find binding pockets before docking.
• admet-prediction: Filter affinity hits by drug-likeness and safety.

Installation

Required Dependencies

# Core (required for all modes)
pip install rdkit-pypi biopython numpy scipy

MM/GBSA capability boundary

rescore.py does not implement full MM/GBSA. Its former OpenMM path omitted the protein–ligand complex and is unavailable. Installing OpenMM does not enable it. For full MM/GBSA, choose a validated external workflow that parameterizes receptor, ligand and complex, records the executed solvent/force-field method, and has independent numerical reference checks. Do not relabel the ligand heuristic or old openmm_mmgbsa files as evidence from such a workflow.

Quick Verification

python -c "from rdkit import Chem; print('RDKit OK')"
python -c "from Bio.PDB import PDBParser; print('BioPython OK')"
python -c "import numpy; print('NumPy OK')"

Core Workflows

Workflow 1: Predict Binding Affinity

Score docked poses using the empirical descriptor-based model.

python scripts/predict.py \
    --protein prepared_protein.pdb \
    --poses docking_results/poses.sdf \
    --output affinity.json

Workflow 2: Ligand Energy Inspection

The retained RDKit fallback is a ligand-only heuristic. It uses no receptor or complex energy, and its output is not a binding free energy or a target-specific score. Select it explicitly; old invocations without a method fail with guidance.

python scripts/rescore.py \
    --method ligand-heuristic \
    --poses poses.sdf \
    --output ligand_scores.json

The output identifies method: ligand_energy_heuristic, receptor_used: false and ligand_score in heuristic units. It records invalid molecules as failed and exits nonzero if any pose failed. Minimization and GB-model flags are unsupported and rejected instead of being silently ignored.

Workflow 3: Consensus Scoring

Combine multiple scoring methods for robust ranking.

python scripts/consensus.py \
    --scores affinity.json ligand_scores.json \
    --docking-scores docking_results/scores.csv \
    --interactions interactions.json \
    --output consensus.json \
    --top-n 10

Workflow 4: Batch Virtual Screening

Screen a compound library against a target.

python scripts/batch.py \
    --protein prepared_protein.pdb \
    --library compounds.sdf \
    --output screening_hits.csv \
    --top-n 50 \
    --threshold 6.0

Workflow 5: Full Pipeline

# 1. Detect pockets
python ../pocket-detection/scripts/detect.py \
    --input protein.pdb --output pockets.json

# 2. Dock ligand
python ../molecular-docking/scripts/dock.py \
    --protein protein.pdb --ligand ligand.sdf \
    --output-dir dock_results/ --method vina \
    --center_x 10 --center_y 20 --center_z 15

# 3. Interaction analysis
python ../molecular-docking/scripts/score.py \
    --protein protein.pdb --poses dock_results/poses.sdf \
    --output interactions.json

# 4. Predict affinity
python scripts/predict.py \
    --protein protein.pdb --poses dock_results/poses.sdf \
    --output affinity.json

# 5. Optional ligand-only energy inspection (not binding affinity)
python scripts/rescore.py \
    --method ligand-heuristic --poses dock_results/poses.sdf \
    --output ligand_scores.json

# 6. Consensus
python scripts/consensus.py \
    --scores affinity.json ligand_scores.json \
    --docking-scores dock_results/scores.csv \
    --interactions interactions.json \
    --output final_ranking.json --top-n 5

Script Reference

Output Format

Affinity JSON (from predict.py)

{
  "protein": "protein.pdb",
  "method": "descriptor",
  "n_poses": 5,
  "note": "Empirical estimate. Typical error: 1-2 log units pKd (~10-100x in Kd). Use for relative ranking only.",
  "predictions": [
    {
      "pose_id": 1,
      "pose_name": "ligand_pose_1",
      "predicted_pKd": 7.2,
      "pKd_uncertainty": 1.5,
      "pKd_range": [5.7, 8.7],
      "predicted_dG_kcal": -9.8,
      "predicted_Kd_nM": 60,
      "confidence": "moderate",
      "features": {
        "mw": 342.4,
        "logp": 2.1,
        "n_hbonds": 4,
        "n_hydrophobic": 12,
        "burial_fraction": 0.65
      }
    }
  ]
}

Consensus JSON (from consensus.py)

{
  "n_poses": 5,
  "sources": ["predict.py", "rescore.py", "dock.py", "score.py"],
  "agreement_tau": 0.72,
  "agreement_class": "high",
  "rankings": [
    {
      "pose_id": 1,
      "pose_name": "ligand_pose_1",
      "consensus_score": 0.85,
      "consensus_rank": 1,
      "individual_ranks": { "predict": 1, "rescore": 2, "docking": 1, "interactions": 3 }
    }
  ]
}

Output Interpretation

pKd Values

Critical: These are computational estimates with ~1-2 log unit uncertainty. A predicted pKd of 7.2 means the true value is likely somewhere between 5.7 and 8.7 (Kd between ~2 nM and 2 uM).

Confidence Levels

Ligand Heuristic Scores

ligand_score is a ligand-only inspection heuristic. Lower scores are ordered first for inspection; this does not imply stronger binding. Consensus retains the ligand_heuristic source name and rejects legacy rescore files with unverified MM/GBSA attribution. Agreement between heuristic rankings is not validation.

Troubleshooting

• All poses get similar scores: The ligands may be too similar, or the scoring function may not discriminate well for this target class.
• Negative confidence: Check if molecules are drug-like (MW 200-600, LogP -1 to 5). Non-drug-like molecules get unreliable scores.
• Full MM/GBSA requested: the bundled script fails explicitly. Use a validated external workflow; do not substitute a ligand-only score.
• Very large library (>10K molecules): Use batch.py with --threshold to filter early.

References

• Wang, R. et al. "The PDBbind database." J. Med. Chem. 47, 2977-2980 (2004).
• Ballester, P.J. & Mitchell, J.B.O. "A machine learning approach to predicting protein-ligand binding affinity." Bioinformatics 26, 1169-1175 (2010).
• Hou, T. et al. "Assessing the performance of the MM/PBSA and MM/GBSA methods." J. Chem. Inf. Model. 51, 69-82 (2011).
• Li, H. et al. "Improving AutoDock Vina Using Random Forest." J. Chem. Inf. Model. 55, 1291-1299 (2015).