Molecular Docking and ML Scoring Functions

Molecular docking is a computational method that predicts the preferred binding pose of a small molecule (ligand) inside a protein binding site, and estimates the free energy of that interaction. It is the workhorse of structure-based virtual screening: given a protein structure and thousands of candidate compounds, docking ranks them by predicted binding affinity, focusing experimental resources on the most promising binders.

Classical docking scoring functions have significant limitations — they are parameterised on limited training sets and make simplifying assumptions about solvation and entropy. Machine learning scoring functions, trained on large collections of protein-ligand co-crystal structures and binding data, can substantially improve pose ranking and binding affinity prediction.

Docking Fundamentals

A docking calculation has two components:

• Search algorithm: explores the conformational space of the ligand inside the binding pocket. Techniques include genetic algorithms (AutoDock), Monte Carlo methods, and fragment-growing.
• Scoring function: evaluates each pose by estimating the free energy of binding.

The output is a ranked list of poses, each with a docking score. Better scores (more negative in most programs) indicate more favorable predicted binding.

Types of Scoring Functions

AutoDock Vina: The Standard Open-Source Docking Engine

AutoDock Vina uses a quasi-Newton gradient optimization search combined with an empirical scoring function. It is free, fast (seconds per ligand), and has been extensively validated. Vina's scoring function sums repulsion, hydrophobic, and hydrogen bonding terms.

Vina Docking Pipeline

"""
Prerequisites:
    pip install vina rdkit meeko
    AutoDock Vina binary in PATH
"""
import subprocess
from pathlib import Path
from rdkit import Chem
from rdkit.Chem import AllChem
import meeko  # Ligand preparation for AutoDock4/Vina

def prepare_ligand_pdbqt(smiles: str, output_path: str) -> bool:
    """
    Convert a SMILES to a 3D-conformer PDBQT file for Vina docking.
    Uses RDKit for 3D generation and Meeko for PDBQT conversion.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return False
    mol = Chem.AddHs(mol)
    result = AllChem.EmbedMolecule(mol, AllChem.ETKDGv3())
    if result != 0:
        return False
    AllChem.MMFFOptimizeMolecule(mol)

    preparator = meeko.MoleculePreparation()
    preparator.prepare(mol)
    preparator.write_pdbqt_file(output_path)
    return True


def run_vina(
    receptor_pdbqt: str,
    ligand_pdbqt: str,
    output_pdbqt: str,
    center_x: float,
    center_y: float,
    center_z: float,
    size_x: float = 20.0,
    size_y: float = 20.0,
    size_z: float = 20.0,
    exhaustiveness: int = 8,
    num_modes: int = 9,
) -> list[dict]:
    """
    Run AutoDock Vina and return a list of pose scores.

    Returns:
        List of dicts with keys: rank, affinity_kcal_per_mol, rmsd_lb, rmsd_ub
    """
    cmd = [
        "vina",
        "--receptor", receptor_pdbqt,
        "--ligand", ligand_pdbqt,
        "--out", output_pdbqt,
        "--center_x", str(center_x),
        "--center_y", str(center_y),
        "--center_z", str(center_z),
        "--size_x", str(size_x),
        "--size_y", str(size_y),
        "--size_z", str(size_z),
        "--exhaustiveness", str(exhaustiveness),
        "--num_modes", str(num_modes),
    ]
    result = subprocess.run(cmd, capture_output=True, text=True)

    poses = []
    for line in result.stdout.split("\n"):
        # Vina output lines look like: "   1        -8.5      0.000      0.000"
        parts = line.split()
        if len(parts) == 4 and parts[0].isdigit():
            poses.append({
                "rank": int(parts[0]),
                "affinity_kcal_per_mol": float(parts[1]),
                "rmsd_lb": float(parts[2]),
                "rmsd_ub": float(parts[3]),
            })
    return poses

High-Throughput Docking Pipeline

import pandas as pd
from pathlib import Path
import tempfile
import os

def screen_library(
    receptor_pdbqt: str,
    compound_df: pd.DataFrame,
    smiles_col: str,
    id_col: str,
    box_center: tuple,  # (x, y, z)
    box_size: tuple = (20, 20, 20),
    output_dir: str = "docking_results",
    exhaustiveness: int = 4,  # Lower for initial triage
) -> pd.DataFrame:
    """
    Dock a compound library against a single receptor.
    Returns a DataFrame with compound IDs and best docking scores.
    """
    os.makedirs(output_dir, exist_ok=True)
    results = []

    for _, row in compound_df.iterrows():
        mol_id = row[id_col]
        smiles = row[smiles_col]

        with tempfile.TemporaryDirectory() as tmpdir:
            lig_pdbqt = os.path.join(tmpdir, "ligand.pdbqt")
            out_pdbqt = os.path.join(output_dir, f"{mol_id}_out.pdbqt")

            success = prepare_ligand_pdbqt(smiles, lig_pdbqt)
            if not success:
                results.append({"id": mol_id, "vina_score": None, "status": "prep_failed"})
                continue

            poses = run_vina(
                receptor_pdbqt=receptor_pdbqt,
                ligand_pdbqt=lig_pdbqt,
                output_pdbqt=out_pdbqt,
                center_x=box_center[0],
                center_y=box_center[1],
                center_z=box_center[2],
                size_x=box_size[0],
                size_y=box_size[1],
                size_z=box_size[2],
                exhaustiveness=exhaustiveness,
                num_modes=3,
            )

            if poses:
                best_score = poses[0]["affinity_kcal_per_mol"]
                results.append({"id": mol_id, "vina_score": best_score, "status": "ok"})
            else:
                results.append({"id": mol_id, "vina_score": None, "status": "no_poses"})

    return pd.DataFrame(results).sort_values("vina_score")

ML-Based Rescoring with GNINA

GNINA (Graph Neural Network In AutoDock) is an open-source docking engine that replaces Vina's empirical scoring with a 3D convolutional neural network trained on the PDBbind dataset. It can be used as:

• A drop-in replacement for Vina (GNINA does both docking and scoring)
• A rescorer applied to poses generated by any docking program

GNINA's CNN affinity score is significantly better than Vina at separating active from inactive compounds in prospective screening campaigns.

"""
Install GNINA: download binary from https://github.com/gnina/gnina/releases
                pip install gnina (Python wrapper)
"""
import subprocess
import pandas as pd

def rescore_with_gnina(
    receptor_pdb: str,
    ligand_sdf: str,
    cnn_scoring: str = "rescore",  # Options: 'rescore', 'refinement', 'docking'
) -> pd.DataFrame:
    """
    Rescore docking poses with GNINA's CNN.

    Args:
        receptor_pdb:  Receptor in PDB format (no ligand, no water)
        ligand_sdf:    Multi-conformer SDF file with docking poses
        cnn_scoring:   How to use CNN ('rescore' = evaluate poses without re-docking)

    Returns:
        DataFrame with CNN affinity and pose CNNscore for each pose.
    """
    result = subprocess.run(
        [
            "gnina",
            "-r", receptor_pdb,
            "-l", ligand_sdf,
            "--score_only",   # Evaluate poses without re-docking
            "--cnn_scoring", cnn_scoring,
            "--out", "/dev/null",
        ],
        capture_output=True,
        text=True,
    )

    rows = []
    for line in result.stdout.split("\n"):
        if line.startswith("CNNscore"):
            parts = line.split()
            rows.append({
                "CNNscore": float(parts[1]),
                "CNNaffinity": float(parts[3]),
                "vinascore": float(parts[5]),
            })
    return pd.DataFrame(rows)

RF-Score: Random Forest Rescoring

RF-Score is a classic ML scoring function that replaces VINA's empirical terms with a Random Forest trained on protein-ligand interaction counts from PDBbind. The feature vector counts interactions between 9 protein atom types and 4 ligand atom types within a 12 Å cutoff, giving a 36-dimensional fingerprint of the binding pose.

import numpy as np
from rdkit import Chem
from rdkit.Chem import rdMolTransforms

# Protein atom type indices (simplified RF-Score v3)
PROTEIN_TYPES = ["C", "N", "O", "S", "P", "F", "Cl", "Br", "I"]
LIGAND_TYPES = ["C", "N", "O", "S"]
DISTANCE_CUTOFF = 12.0  # Angstroms

def compute_rfscore_features(
    protein_mol: Chem.Mol,
    ligand_mol: Chem.Mol,
    interaction_cutoff: float = DISTANCE_CUTOFF
) -> np.ndarray:
    """
    Compute 36-dimensional RF-Score feature vector from a protein-ligand complex.
    Counts protein-ligand atom pair interactions by type within cutoff.
    """
    n_p = len(PROTEIN_TYPES)
    n_l = len(LIGAND_TYPES)
    features = np.zeros(n_p * n_l, dtype=np.float32)

    prot_conf = protein_mol.GetConformer()
    lig_conf = ligand_mol.GetConformer()

    for p_atom in protein_mol.GetAtoms():
        p_symbol = p_atom.GetSymbol()
        if p_symbol not in PROTEIN_TYPES:
            continue
        p_idx_type = PROTEIN_TYPES.index(p_symbol)
        p_pos = prot_conf.GetAtomPosition(p_atom.GetIdx())

        for l_atom in ligand_mol.GetAtoms():
            l_symbol = l_atom.GetSymbol()
            if l_symbol not in LIGAND_TYPES:
                continue
            l_idx_type = LIGAND_TYPES.index(l_symbol)
            l_pos = lig_conf.GetAtomPosition(l_atom.GetIdx())

            dist = ((p_pos.x - l_pos.x)**2 +
                    (p_pos.y - l_pos.y)**2 +
                    (p_pos.z - l_pos.z)**2) ** 0.5

            if dist <= interaction_cutoff:
                feature_idx = p_idx_type * n_l + l_idx_type
                features[feature_idx] += 1.0

    return features

Evaluating Docking Performance

The two key evaluation metrics for docking in virtual screening are:

import numpy as np
from sklearn.metrics import roc_auc_score

def enrichment_factor(y_true: np.ndarray, scores: np.ndarray, fraction: float = 0.01) -> float:
    """
    Compute enrichment factor at a given fraction of the ranked library.

    Args:
        y_true:    Binary array (1=active, 0=inactive)
        scores:    Docking scores (lower = better binding in Vina convention)
        fraction:  Top fraction to evaluate (0.01 = top 1%)

    Returns:
        Enrichment factor (EF_x%)
    """
    n_total = len(y_true)
    n_actives = y_true.sum()
    n_top = max(1, int(n_total * fraction))

    # Sort by score — more negative = better binding for Vina
    sorted_idx = np.argsort(scores)
    top_idx = sorted_idx[:n_top]
    n_actives_top = y_true[top_idx].sum()

    ef = (n_actives_top / n_top) / (n_actives / n_total)
    return round(ef, 2)


# Example
docking_results = pd.read_csv("dud_e_docking_results.csv")  # id, score, active
y = docking_results["active"].values
# Vina scores are negative; negate so higher = better for sklearn
scores_neg = -docking_results["score"].values

print(f"EF 1%:   {enrichment_factor(y, -scores_neg, 0.01)}")
print(f"EF 5%:   {enrichment_factor(y, -scores_neg, 0.05)}")
print(f"ROC-AUC: {roc_auc_score(y, scores_neg):.3f}")

Practical Docking Workflow for Virtual Screening

1. Receptor preparation: remove water molecules (keep structural ones), add hydrogens at physiological pH, assign partial charges, export to PDBQT (Vina) or MOL2/PDB (GNINA).
2. Binding site definition: define the docking box around the known active site. If the structure has no co-crystallised ligand, use SiteMap or fpocket to predict binding pockets.
3. Ligand preparation: standardise SMILES (remove salts), enumerate ionisation states at target pH with Epik or RDKit MolStandardize, generate 3D conformers, assign partial charges.
4. ADMET pre-filter: apply Lipinski/Veber rules, PAINS filter, and catalytic group filters before docking to reduce library size 4–10× and save compute.
5. Docking and scoring: dock with Vina (fast, open-source) or Glide (commercial, higher accuracy). Optionally rescore top poses with GNINA CNN.
6. Post-processing: cluster poses by binding mode, visualise top 50–100 in PyMOL/UCSF ChimeraX, check interactions (H-bonds, hydrophobic contacts, key water-mediated interactions).

Key Takeaways

• AutoDock Vina is the standard open-source docking engine; its empirical scoring function is fast and widely validated for initial triage.
• GNINA's CNN-based scoring substantially outperforms Vina's empirical score for separating actives from decoys, especially for scaffold-diverse libraries.
• Always apply ADMET pre-filters before docking — reducing the library from 1 million to 100,000 compounds saves 10× compute with no loss of chemically interesting candidates.
• Report enrichment factor at 1% (EF1%) and BEDROC, not just ROC-AUC. A docking campaign succeeds if the top-ranked 1% is highly enriched for actives, not if the overall ranking is perfect.
• Pose accuracy (RMSD ≤ 2 Å from crystal structure) must be validated separately from enrichment. A scoring function can enrich actives without predicting the correct binding pose.

Try it live

Explore a related model hands-on: open Molecular Docking: Ligand Binding Search and change its parameters while it runs. Everything stays in your browser.