Protein-Protein Interactions in Drug Discovery

Protein-protein interactions (PPIs) control nearly every cellular signalling process, making them a major class of therapeutic targets. However, PPI interfaces are notoriously difficult to drug with small molecules: interfaces are typically large (1,500–3,000 Ų), flat, and lack the defined pockets that enzymes and GPCRs provide for small molecule binding.

Despite these challenges, several approved drugs target PPIs — Venetoclax (BCL2/BAX-BIM interface), Navitoclax (BCL-XL), ABT-737 — and the field has matured significantly with improved understanding of hot-spot residues and the discovery of "cryptic" pockets that open transiently in MD simulations. This guide covers the computational toolkit for PPI interface analysis and inhibitor discovery.

Interface Characterisation

The first step is to characterise the PPI interface from a co-crystal structure of the protein complex (or a homology model). Key metrics:

from Bio.PDB import PDBParser, SASA
from Bio.PDB.SASA import ShrakeRupley
import numpy as np
from typing import Optional

def compute_interface_area(
    complex_pdb: str,
    chain_a: str = "A",
    chain_b: str = "B",
    probe_radius: float = 1.4,
) -> dict:
    """
    Compute the protein-protein interface area from a PDB co-crystal structure.

    Interface area = (SASA_A + SASA_B - SASA_complex) / 2
    where SASA is solvent-accessible surface area.

    Args:
        complex_pdb: Path to PDB file with the protein-protein complex
        chain_a: Chain ID for protein A
        chain_b: Chain ID for protein B
        probe_radius: Radius of solvent probe in Angstroms (1.4 Å = water)

    Returns:
        Dict with SASA values and interface area estimate.
    """
    parser = PDBParser(QUIET=True)
    structure = parser.get_structure("complex", complex_pdb)
    model = structure[0]

    sr = ShrakeRupley(probe_radius=probe_radius)

    # Get SASA of full complex
    sr.compute(model, level="R")
    sasa_complex = sum(
        r.sasa for r in model.get_residues()
        if r.get_parent().id in [chain_a, chain_b]
    )

    # Get SASA of chain A alone
    chain_a_model = model[chain_a]
    sr.compute(chain_a_model, level="R")
    sasa_a = sum(r.sasa for r in chain_a_model.get_residues())

    # Get SASA of chain B alone
    chain_b_model = model[chain_b]
    sr.compute(chain_b_model, level="R")
    sasa_b = sum(r.sasa for r in chain_b_model.get_residues())

    interface_area = (sasa_a + sasa_b - sasa_complex) / 2.0

    return {
        "sasa_complex": round(sasa_complex, 2),
        "sasa_chain_a": round(sasa_a, 2),
        "sasa_chain_b": round(sasa_b, 2),
        "interface_area_A2": round(interface_area, 2),
        "tractability": "tractable" if interface_area < 1500 else "challenging",
    }


def identify_interface_residues(
    complex_pdb: str,
    chain_a: str = "A",
    chain_b: str = "B",
    distance_cutoff: float = 5.0,
) -> dict:
    """
    Identify residues at the protein-protein interface by proximity criterion.
    Interface residue: any residue with at least one atom within distance_cutoff
    of the partner chain.

    Returns:
        Dict with chain A and chain B interface residue lists (residue_id, residue_name).
    """
    parser = PDBParser(QUIET=True)
    structure = parser.get_structure("complex", complex_pdb)
    model = structure[0]

    def get_atoms(chain_id):
        return list(model[chain_id].get_atoms())

    atoms_a = get_atoms(chain_a)
    atoms_b = get_atoms(chain_b)

    interface_a = set()
    interface_b = set()

    for atom_a in atoms_a:
        for atom_b in atoms_b:
            dist = atom_a - atom_b  # BioPython Atom subtraction = Euclidean distance
            if dist < distance_cutoff:
                res_a = atom_a.get_parent()
                res_b = atom_b.get_parent()
                interface_a.add((res_a.get_id()[1], res_a.get_resname()))
                interface_b.add((res_b.get_id()[1], res_b.get_resname()))

    return {
        "interface_residues_chain_a": sorted(list(interface_a)),
        "interface_residues_chain_b": sorted(list(interface_b)),
        "n_interface_residues": len(interface_a) + len(interface_b),
    }

Hot-Spot Residue Prediction

Hot-spot residues contribute disproportionally to the binding free energy of the complex (their alanine-scan ΔΔG > 2 kcal/mol). Computationally, hot-spots can be predicted from structural features using ML models trained on experimental alanine scanning databases (SKEMPI, ASEdb).

import numpy as np
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.preprocessing import StandardScaler
from rdkit.Chem import rdMolDescriptors

# Hot-spot residue features (computed per interface residue)
def extract_hotspot_features(
    residue_name: str,
    position: int,
    bfactor: float,
    sasa_free: float,
    sasa_complex: float,
    n_neighbors_5A: int,
    n_neighbors_10A: int,
    is_aromatic: bool,
    is_charged: bool,
) -> np.ndarray:
    """
    Extract per-residue features for hot-spot prediction.

    Args:
        residue_name: 3-letter residue code
        position: Residue sequence position
        bfactor: B-factor / pLDDT from structure
        sasa_free: SASA in the free (unbound) state
        sasa_complex: SASA in the complex
        n_neighbors_5A: Interface residue neighbors within 5 Å
        n_neighbors_10A: Interface residue neighbors within 10 Å
        is_aromatic: Whether residue is aromatic (Phe, Tyr, Trp, His)
        is_charged: Whether residue is charged at pH 7.4

    Returns:
        1D numpy feature vector.
    """
    # Burial on binding (|ΔSASA|)
    delta_sasa = abs(sasa_free - sasa_complex)

    # One-hot encode aromatic hydrophobic residue types (known hot-spot anchors)
    # Trp, Tyr, Phe contribute most to PPI interfaces
    aromatic_types = {"TRP": [1, 0, 0], "TYR": [0, 1, 0], "PHE": [0, 0, 1]}
    residue_type_feat = aromatic_types.get(residue_name, [0, 0, 0])

    features = np.array([
        bfactor,
        sasa_free,
        sasa_complex,
        delta_sasa,
        n_neighbors_5A,
        n_neighbors_10A,
        float(is_aromatic),
        float(is_charged),
    ] + residue_type_feat, dtype=np.float32)

    return features


class HotspotPredictor:
    """
    ML model to predict if an interface residue is a hot-spot.
    Trained on SKEMPI v2 experimental alanine scanning data.
    """
    def __init__(self):
        self.model = GradientBoostingClassifier(
            n_estimators=200, max_depth=4, learning_rate=0.05,
            min_samples_leaf=5, random_state=42
        )
        self.scaler = StandardScaler()

    def fit(self, X: np.ndarray, y: np.ndarray) -> None:
        """
        Train on alanine scan data.
        y = 1 if ΔΔG > 2.0 kcal/mol (hot-spot), else 0.
        """
        X_scaled = self.scaler.fit_transform(X)
        self.model.fit(X_scaled, y)

    def predict_proba(self, X: np.ndarray) -> np.ndarray:
        """
        Return hot-spot probability for each residue.
        """
        return self.model.predict_proba(self.scaler.transform(X))[:, 1]

    def predict(self, X: np.ndarray, threshold: float = 0.5) -> np.ndarray:
        return (self.predict_proba(X) >= threshold).astype(int)

iPPI-DB Dataset and PPI Inhibitor Properties

iPPI-DB (iPPI-DB.org) is a curated database of small-molecule inhibitors and stabilisers of PPIs. It contains ~800 compounds with measured PPI modulation activity, their iPPI-score (a drug-likeness score specific to PPI inhibitors), and the target PPI class. Key observations from iPPI-DB:

from rdkit import Chem
from rdkit.Chem import Descriptors, rdMolDescriptors
import numpy as np

def ippiscore(smiles: str) -> float:
    """
    Compute the iPPI-score for a molecule (Suckling et al. adapted).
    The iPPI-score is a drug-likeness score tailored for PPI inhibitors.

    Score = sum of weighted contributions from:
      - MW (optimal 450–600 Da for PPI)
      - cLogP (optimal 3–6)
      - Aromatic rings (optimal 3–5)
      - Stereocenters (optimal 2–4)
      - Ring fraction (high preferred for PPI)

    Returns:
        Score between 0 (poor PPI inhibitor candidate) and 1 (ideal).
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return 0.0

    mw = Descriptors.ExactMolWt(mol)
    clogp = Descriptors.MolLogP(mol)
    n_arom = rdMolDescriptors.CalcNumAromaticRings(mol)
    n_stereo = len(Chem.FindMolChiralCenters(mol, includeUnassigned=True))
    n_heavy = mol.GetNumHeavyAtoms()
    n_ring_atoms = sum(atom.IsInRing() for atom in mol.GetAtoms())
    ring_frac = n_ring_atoms / n_heavy if n_heavy > 0 else 0.0

    def gaussian_score(x, mu, sigma):
        return np.exp(-0.5 * ((x - mu) / sigma) ** 2)

    # Optimal ranges centred at midpoint of iPPI-DB distributions
    score_mw = gaussian_score(mw, mu=525, sigma=100)
    score_logp = gaussian_score(clogp, mu=4.5, sigma=1.5)
    score_arom = gaussian_score(n_arom, mu=4.0, sigma=1.5)
    score_stereo = gaussian_score(n_stereo, mu=3.0, sigma=1.5)
    score_ring = min(ring_frac / 0.7, 1.0)  # reward compact rings

    # Uniform weights (could be fitted on iPPI-DB with logistic regression)
    ippiscore = np.mean([score_mw, score_logp, score_arom, score_stereo, score_ring])
    return round(float(ippiscore), 4)

Interface Fingerprint for Virtual Screening

from rdkit import Chem
from rdkit.Chem import AllChem, DataStructs
import numpy as np

def interface_mimicry_score(
    compound_smiles: str,
    hot_spot_pharmacophore_smarts: list[str],
) -> dict:
    """
    Score how well a compound mimics hot-spot residue pharmacophore patterns.

    This simplified approach counts how many hot-spot pharmacophore features
    (aromatic, H-bond donor, hydrophobic) are represented in the compound.

    Args:
        compound_smiles: SMILES of the candidate inhibitor
        hot_spot_pharmacophore_smarts: List of SMARTS patterns encoding
            the pharmacophoric features of interface hot-spot residues.
            Example: ['c1ccc[nH]1'] for Trp indole mimicry.

    Returns:
        Dict with matched features count, fraction covered, and SMARTS matches.
    """
    mol = Chem.MolFromSmiles(compound_smiles)
    if mol is None:
        return {"error": "Invalid SMILES"}

    matches = {}
    for smarts in hot_spot_pharmacophore_smarts:
        query = Chem.MolFromSmarts(smarts)
        if query:
            hit = mol.HasSubstructMatch(query)
            matches[smarts] = bool(hit)

    n_matched = sum(matches.values())
    fraction = n_matched / len(hot_spot_pharmacophore_smarts) if hot_spot_pharmacophore_smarts else 0.0

    return {
        "n_features_matched": n_matched,
        "n_features_total": len(hot_spot_pharmacophore_smarts),
        "mimicry_fraction": round(fraction, 3),
        "matched_features": {k: v for k, v in matches.items() if v},
    }


# Common hot-spot pharmacophore SMARTS for known PPI targets
COMMON_PPI_PHARMACOPHORES = {
    "Trp_indole": "c1ccc2[nH]ccc2c1",
    "Tyr_phenol":  "c1ccc(O)cc1",
    "Phe_benzene": "c1ccccc1",
    "Leu_isobutyl": "CC(C)CC",
    "Glu_carboxylate": "C(=O)[OH]",
    "His_imidazole": "c1cn[nH]c1",
    "Arg_guanidinium": "C(=N)(N)N",
}

PPI Druggability Assessment Pipeline

from dataclasses import dataclass
from typing import Optional

@dataclass
class PPIDruggabilityReport:
    ppi_name: str
    interface_area_A2: float
    n_hotspot_residues: int
    top_pocket_volume_A3: float
    top_pocket_druggability: float
    has_known_inhibitor: bool
    ippidb_max_pchembl: Optional[float]
    recommended_strategy: str
    tractability: str

def assess_ppi_druggability(
    interface_area: float,
    n_hotspots: int,
    max_pocket_volume: float,
    pockscore: float,
    has_inhibitor: bool = False,
    ippidb_pchembl: Optional[float] = None,
) -> str:
    """
    Classify PPI druggability based on structural and database evidence.

    Returns a recommended strategy string.
    """
    if has_inhibitor and ippidb_pchembl and ippidb_pchembl >= 6.0:
        return "Validated target: optimise known chemotype using iPPI-DB SAR"

    if interface_area < 1000 and pockscore > 0.5 and n_hotspots >= 3:
        return "Favourable geometry: proceed with direct docking and fragment screen"

    if max_pocket_volume > 300 and pockscore > 0.3:
        return "Cryptic pocket: consider MD-based pocket discovery (HTMD, PocketMiner)"

    if n_hotspots >= 5 and interface_area > 2000:
        return "Large flat interface: consider peptide/macrocycle mimetics or PROTACs"

    return "Challenging PPI: validate via alanine scanning + biophysical screen (SPR/NMR)"

Key Takeaways

• The single most important predictor of PPI tractability is the presence of a "hot-spot groove": an interface area <1,500 Ų with ≥3 exposed hot-spot residues (particularly Trp/Tyr) clustered around a pocket that scores >0.4 with fpocket.
• Successful PPI inhibitors in iPPI-DB are on average heavier (MW ~525), more aromatic, and more lipophilic than standard drugs — the iPPI-score captures these differences and should complement Lipinski Ro5 for PPI-specific screening.
• When no clear pocket is visible in the co-crystal structure, run 100–200 ns MD simulations and analyse trajectory frames with fpocket to detect cryptic pockets that open transiently during dynamics.
• Hot-spot prediction from the SKEMPI database (using features like ΔΔG_alanine, SASA burial, N-interacting partners) achieves MCC ~0.4–0.5, sufficient for guiding which residues to mimic or block.
• Fragment-based approaches are particularly successful for PPIs: fragments satisfy the requirement to bind in a small hot-spot sub-region, and growing/merging can yield complete interface coverage.

Try it live

Explore a related model hands-on: open Protein Folding Visualiser and change its parameters while it runs. Everything stays in your browser.