Binding Affinity Prediction

Binding affinity — the strength with which a drug candidate binds to its protein target — is the central quantity in drug discovery. It is experimentally measured as the inhibition constant K_i, dissociation constant K_d, or half-maximal inhibitory concentration IC₅₀, and reported as its negative log (pK_i, pK_d, pIC₅₀) so that higher values mean stronger binding. The goal of computational binding affinity prediction is to rank candidate compounds by predicted affinity before (or to replace) expensive biochemical assays.

There are two main paradigms: ligand-based models (QSAR) that learn from the activity of known compounds without using the protein structure, and structure-based models that use explicit 3D information about the protein-ligand complex. This guide focuses on structure-based models using data from protein-ligand co-crystal structures in PDBbind.

Datasets

The PDBbind core set (CASF benchmark) is the gold-standard evaluation set. Models are evaluated by Pearson R on the 290 complexes and by scoring power (ability to rank compounds for a single target).

Featurization of Protein-Ligand Complexes

Structure-based models need to encode the 3D protein-ligand interaction environment. The main approaches are:

• Interaction fingerprints (IFP): binary vector encoding which types of interactions (H-bond, hydrophobic contact, π-stacking, salt bridge) exist between specific residues and the ligand. Fast, interpretable, but loses 3D geometry.
• Atom-pair distance features (RF-Score): count atom pair contacts by type within a distance cutoff. The 36 RF-Score features (9 protein atom types × 4 ligand atom types) are the simplest structure-based features.
• 3D Voxel grids (3D-CNN): discretise the binding pocket into a 3D grid of voxels; each voxel encodes atom type channels. Captures full spatial geometry but is memory-intensive.
• Molecular graphs (GNN): represent the complex as a heterogeneous graph with protein atoms, ligand atoms, and edges for covalent bonds and proximity contacts.

Interaction Fingerprint (Protein-Ligand IFP)

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

# Interaction type definitions
INTERACTION_TYPES = [
    "hbond_donor",
    "hbond_acceptor",
    "hydrophobic",
    "aromatic_face",
    "ionic_positive",
    "ionic_negative",
    "water_bridge",
]

def compute_interaction_fingerprint(
    protein_mol: Chem.Mol,
    ligand_mol: Chem.Mol,
    distance_cutoff: float = 4.5,
    n_residue_bins: int = 50,
) -> np.ndarray:
    """
    Compute a simplified interaction fingerprint (IFP) for a protein-ligand complex.
    Each bit encodes (residue_bin, interaction_type) pairs within distance_cutoff.

    Returns a binary vector of length n_residue_bins * len(INTERACTION_TYPES).
    """
    fp_size = n_residue_bins * len(INTERACTION_TYPES)
    ifp = np.zeros(fp_size, dtype=np.uint8)

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

    for p_atom in protein_mol.GetAtoms():
        res_num = p_atom.GetMonomerInfo().GetResidueNumber() if p_atom.GetMonomerInfo() else 0
        res_bin = res_num % n_residue_bins
        p_pos = prot_conf.GetAtomPosition(p_atom.GetIdx())

        for l_atom in ligand_mol.GetAtoms():
            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 > distance_cutoff:
                continue
            # Simplified interaction type assignment
            p_sym = p_atom.GetSymbol()
            l_sym = l_atom.GetSymbol()
            if p_sym == "O" and l_sym == "N":
                itype = 0  # hbond_donor
            elif p_sym == "N" and l_sym == "O":
                itype = 1  # hbond_acceptor
            elif p_sym == "C" and l_sym == "C":
                itype = 2  # hydrophobic
            else:
                continue
            bit_idx = res_bin * len(INTERACTION_TYPES) + itype
            ifp[bit_idx] = 1

    return ifp

RF-Score: Random Forest on Atom-Pair Contacts

import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error
from scipy.stats import pearsonr
import pandas as pd

# Protein atom types for RF-Score v3
PROT_TYPES = ["C", "N", "O", "S", "P", "F", "Cl", "Br", "I"]
LIG_TYPES = ["C", "N", "O", "S"]
CUTOFFS = [2.0, 4.0, 6.0, 8.0, 10.0, 12.0]  # RF-Score v3 uses distance bins

def rfscore_v3_features(prot_mol, lig_mol) -> np.ndarray:
    """
    RF-Score v3: 36-dim feature/cutoff = 36*6 = 216-dim feature vector.
    Counts atom pair contacts for each (protein_type, ligand_type) pair
    at each distance cutoff threshold.
    """
    n_pairs = len(PROT_TYPES) * len(LIG_TYPES)
    features = np.zeros(n_pairs * len(CUTOFFS), dtype=np.float32)

    prot_conf = prot_mol.GetConformer()
    lig_conf = lig_mol.GetConformer()

    for p_atom in prot_mol.GetAtoms():
        p_sym = p_atom.GetSymbol()
        if p_sym not in PROT_TYPES:
            continue
        p_type = PROT_TYPES.index(p_sym)
        p_pos = prot_conf.GetAtomPosition(p_atom.GetIdx())

        for l_atom in lig_mol.GetAtoms():
            l_sym = l_atom.GetSymbol()
            if l_sym not in LIG_TYPES:
                continue
            l_type = LIG_TYPES.index(l_sym)
            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

            pair_idx = p_type * len(LIG_TYPES) + l_type
            for c_idx, cutoff in enumerate(CUTOFFS):
                if dist <= cutoff:
                    features[c_idx * n_pairs + pair_idx] += 1.0

    return features


# Training RF-Score on PDBbind
def train_rfscore(complexes: list[dict]) -> RandomForestRegressor:
    """
    Train RF-Score v3 on a list of protein-ligand complex dicts.
    Each dict: {'protein_mol': Chem.Mol, 'ligand_mol': Chem.Mol, 'affinity': float}
    """
    X, y = [], []
    for cpx in complexes:
        feat = rfscore_v3_features(cpx["protein_mol"], cpx["ligand_mol"])
        X.append(feat)
        y.append(cpx["affinity"])

    X, y = np.vstack(X), np.array(y)
    model = RandomForestRegressor(n_estimators=500, n_jobs=-1, random_state=42)
    model.fit(X, y)
    return model

3D-CNN Model for Affinity Prediction

3D convolutional neural networks (3D-CNNs) can learn spatial patterns of atom distributions in binding pockets directly from voxel grids. Models like GNINA's CNN, PointNet-based architectures, and Grid-based deep learning have shown Pearson R > 0.8 on PDBbind core.

import torch
import torch.nn as nn
import numpy as np
from rdkit import Chem

ATOM_TYPES = ["C", "N", "O", "S", "P", "F", "Cl", "Br", "H"]
GRID_RESOLUTION = 1.0   # Angstroms per voxel
GRID_SIZE = 24           # 24 x 24 x 24 Angstrom box

def complex_to_voxel(
    prot_mol: Chem.Mol,
    lig_mol: Chem.Mol,
    grid_size: int = GRID_SIZE,
    resolution: float = GRID_RESOLUTION,
) -> np.ndarray:
    """
    Convert a protein-ligand complex to a 3D voxel grid.
    Returns array of shape (n_channels, grid_size, grid_size, grid_size).
    n_channels = len(ATOM_TYPES) * 2  (protein channels + ligand channels)
    """
    n_channels = len(ATOM_TYPES) * 2
    grid = np.zeros((n_channels, grid_size, grid_size, grid_size), dtype=np.float32)
    center = _compute_ligand_center(lig_mol)

    for is_ligand, mol in enumerate([prot_mol, lig_mol]):
        conf = mol.GetConformer()
        channel_offset = is_ligand * len(ATOM_TYPES)
        for atom in mol.GetAtoms():
            pos = conf.GetAtomPosition(atom.GetIdx())
            sym = atom.GetSymbol()
            if sym not in ATOM_TYPES:
                continue
            ch = ATOM_TYPES.index(sym) + channel_offset
            # Convert to grid coordinates centered on ligand
            xi = int((pos.x - center[0]) / resolution + grid_size / 2)
            yi = int((pos.y - center[1]) / resolution + grid_size / 2)
            zi = int((pos.z - center[2]) / resolution + grid_size / 2)
            if 0 <= xi < grid_size and 0 <= yi < grid_size and 0 <= zi < grid_size:
                grid[ch, xi, yi, zi] += 1.0
    return grid


def _compute_ligand_center(mol: Chem.Mol) -> tuple:
    conf = mol.GetConformer()
    positions = [conf.GetAtomPosition(i) for i in range(mol.GetNumAtoms())]
    cx = sum(p.x for p in positions) / len(positions)
    cy = sum(p.y for p in positions) / len(positions)
    cz = sum(p.z for p in positions) / len(positions)
    return cx, cy, cz


class CNN3DAffinityModel(nn.Module):
    """3D-CNN for protein-ligand binding affinity prediction."""

    def __init__(self, n_channels: int = 18, grid_size: int = GRID_SIZE):
        super().__init__()
        self.encoder = nn.Sequential(
            nn.Conv3d(n_channels, 32, kernel_size=3, padding=1),
            nn.BatchNorm3d(32),
            nn.ReLU(),
            nn.MaxPool3d(2),

            nn.Conv3d(32, 64, kernel_size=3, padding=1),
            nn.BatchNorm3d(64),
            nn.ReLU(),
            nn.MaxPool3d(2),

            nn.Conv3d(64, 128, kernel_size=3, padding=1),
            nn.BatchNorm3d(128),
            nn.ReLU(),
            nn.AdaptiveAvgPool3d(2),
        )
        self.head = nn.Sequential(
            nn.Flatten(),
            nn.Linear(128 * 8, 256),
            nn.ReLU(),
            nn.Dropout(0.3),
            nn.Linear(256, 1),
        )

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        return self.head(self.encoder(x)).squeeze(-1)

CASF Benchmark Evaluation

from scipy.stats import pearsonr, spearmanr
from sklearn.metrics import mean_squared_error
import numpy as np

def evaluate_affinity_model(y_true: np.ndarray, y_pred: np.ndarray) -> dict:
    """
    Evaluate a binding affinity model on the CASF benchmark metrics.

    Args:
        y_true: Experimental pKd/pKi values
        y_pred: Predicted values

    Returns:
        Dict with Pearson R, Spearman r, RMSE, MAE, and CI (concordance index)
    """
    rmse = mean_squared_error(y_true, y_pred, squared=False)
    mae = np.mean(np.abs(y_true - y_pred))
    r_p, _ = pearsonr(y_true, y_pred)
    r_s, _ = spearmanr(y_true, y_pred)

    # Concordance Index: fraction of pairs correctly ranked
    n = len(y_true)
    concordant = 0
    total_pairs = 0
    for i in range(n):
        for j in range(i + 1, n):
            delta_true = y_true[i] - y_true[j]
            delta_pred = y_pred[i] - y_pred[j]
            if abs(delta_true) < 0.1:
                continue  # Ignore ties
            total_pairs += 1
            if delta_true * delta_pred > 0:
                concordant += 1
    ci = concordant / total_pairs if total_pairs > 0 else 0.0

    return {
        "pearson_r": round(r_p, 4),
        "spearman_r": round(r_s, 4),
        "rmse": round(rmse, 4),
        "mae": round(mae, 4),
        "concordance_index": round(ci, 4),
    }

Free Energy Methods: Context and Limits

Free Energy Perturbation (FEP) and Thermodynamic Integration (TI) are physics-based methods that estimate binding free energy by alchemically transforming one ligand into another in MD simulations. FEP+ (Schrödinger) achieves ~1 kcal/mol accuracy for congeneric series (R² ≈ 0.7–0.9) but requires:

• A high-quality protein-ligand co-crystal structure
• Molecular mechanics force field parameters for the ligand
• ~50–500 ns of simulation per perturbation (GPU days of compute)
• Congeneric series: ligands that differ by small fragment changes (~1–3 heavy atoms)

ML binding affinity models complement FEP in the workflow: ML screens millions of compounds quickly, and FEP is applied to the top handful of predicted binders to obtain reliable ΔΔG estimates before synthesis.

Key Takeaways

• RF-Score v3 (Pearson R ≈ 0.75 on CASF) is an excellent structure-based baseline that trains in minutes on PDBbind refined.
• 3D-CNN and equivariant GNN models achieve R ≈ 0.82–0.83 but require GPU training and careful voxelisation or graph featurisation of the complex.
• The practical ceiling imposed by assay noise and protein flexibility means Pearson R > 0.85 on prospective data is rarely achievable without FEP-level physics.
• ML affinity models are most powerful as ranking tools in a virtual screening funnel — not as absolute predictors of binding free energy.
• Always use the CASF-2016 core set for evaluation to enable fair comparison with published methods.

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.