Allosteric Drug Design with ML

Allostery is the phenomenon by which ligand binding at one site modulates protein activity at a distant site via conformational coupling. Unlike orthosteric inhibitors that compete with the natural substrate (often highly conserved across family members), allosteric sites show greater sequence divergence, enabling exquisite selectivity. Paradigm examples include the MEK1/2 allosteric pocket targeted by trametinib, the ABL myristoyl pocket targeted by asciminib, and the mGluR negative allosteric modulators.

The central ML challenge is that allosteric sites are often cryptic — not visible in apo crystal structures — and allosteric ligand databases are orders of magnitude smaller than orthosteric ones.

Protein Communication Network Analysis

import numpy as np
import networkx as nx
from typing import Optional


def build_residue_contact_map(
    distance_matrix: np.ndarray,
    ca_only: bool = True,
    cutoff_angstrom: float = 8.0,
) -> np.ndarray:
    """
    Convert a residue-residue distance matrix to a binary contact map.

    Args:
        distance_matrix: Square matrix of Cα-Cα distances (Å), shape (N, N)
        ca_only: Distances are already Cα-only if True
        cutoff_angstrom: Distance threshold for contact definition

    Returns:
        Binary adjacency matrix (N × N).
    """
    contact_map = (distance_matrix < cutoff_angstrom).astype(int)
    np.fill_diagonal(contact_map, 0)
    return contact_map


def compute_communication_centrality(
    contact_map: np.ndarray,
    active_site_residues: list[int],
    weight_by_coevolution: Optional[np.ndarray] = None,
) -> dict:
    """
    Compute betweenness centrality on the protein residue graph.

    Residues with high betweenness on shortest paths between the active site
    and the rest of the protein are candidate allosteric communication hubs.

    Args:
        contact_map: Binary adjacency matrix from build_residue_contact_map()
        active_site_residues: Indices of known catalytic/binding residues
        weight_by_coevolution: Optional coevolution scores for edge weights
                               (higher coevolution = stronger coupling = lower weight)

    Returns:
        Dict with betweenness centrality per residue and top allosteric candidates.
    """
    G = nx.from_numpy_array(contact_map)

    if weight_by_coevolution is not None:
        for u, v in G.edges():
            cov = weight_by_coevolution[u, v]
            G[u][v]["weight"] = 1.0 / (cov + 1e-6)

    betweenness = nx.betweenness_centrality(G, weight="weight" if weight_by_coevolution is not None else None, normalized=True)
    closeness = nx.closeness_centrality(G)

    n_residues = contact_map.shape[0]
    active_set = set(active_site_residues)
    allosteric_candidates = {
        r: betweenness[r]
        for r in range(n_residues)
        if r not in active_set and betweenness[r] > np.percentile(list(betweenness.values()), 90)
    }

    sorted_candidates = sorted(allosteric_candidates.items(), key=lambda x: x[1], reverse=True)
    return {
        "betweenness": betweenness,
        "closeness": closeness,
        "top_allosteric_candidates": sorted_candidates[:20],
    }

Cryptic Pocket Detection from MD Trajectories

import numpy as np
from dataclasses import dataclass


@dataclass
class PocketInfo:
    pocket_id: int
    centroid_residues: list[int]
    mean_volume_A3: float
    max_volume_A3: float
    druggability_score: float
    is_cryptic: bool


def detect_cryptic_pockets_from_ensemble(
    trajectory_volumes: np.ndarray,
    pocket_ids: list[int],
    apo_volume_threshold: float = 50.0,
    holo_volume_threshold: float = 300.0,
    min_open_fraction: float = 0.05,
) -> list[PocketInfo]:
    """
    Identify cryptic pockets that open transiently in MD trajectories.

    Cryptic pockets are absent or very small in the apo crystal structure
    but can be sampled in unbiased or accelerated MD simulations.

    Args:
        trajectory_volumes: Array of shape (n_pockets, n_frames) — volume per pocket per frame
        pocket_ids: Identifier for each pocket row
        apo_volume_threshold: Max volume in apo (frame 0) to classify as "cryptic"
        holo_volume_threshold: Min peak volume for the pocket to be "druggable" when open
        min_open_fraction: Minimum fraction of trajectory frames with pocket open

    Returns:
        List of PocketInfo for cryptic pockets meeting druggability criteria.
    """
    cryptic_pockets = []
    n_pockets, n_frames = trajectory_volumes.shape

    for i, pid in enumerate(pocket_ids):
        vols = trajectory_volumes[i]
        apo_vol = vols[0]
        max_vol = vols.max()
        mean_vol = vols.mean()
        n_open = np.sum(vols >= holo_volume_threshold)
        open_fraction = n_open / n_frames

        is_cryptic = apo_vol < apo_volume_threshold and max_vol > holo_volume_threshold
        meets_criteria = is_cryptic and open_fraction >= min_open_fraction

        if meets_criteria:
            druggability = min(max_vol / 500.0, 1.0) * open_fraction
            cryptic_pockets.append(PocketInfo(
                pocket_id=pid,
                centroid_residues=[],
                mean_volume_A3=round(float(mean_vol), 1),
                max_volume_A3=round(float(max_vol), 1),
                druggability_score=round(druggability, 3),
                is_cryptic=True,
            ))

    return sorted(cryptic_pockets, key=lambda p: p.druggability_score, reverse=True)


def allosteric_score_from_pocket(
    pocket_betweenness: float,
    pocket_druggability: float,
    pocket_coevolution: float,
    w_betweenness: float = 0.4,
    w_druggability: float = 0.4,
    w_coevolution: float = 0.2,
) -> float:
    """
    Composite allosteric site score combining structural, dynamic, and evolutionary evidence.

    Inputs should be normalised to [0, 1].
    """
    return (
        w_betweenness * pocket_betweenness
        + w_druggability * pocket_druggability
        + w_coevolution * pocket_coevolution
    )

Known Allosteric Drug Examples

Sequence-Based Allosteric Site Prediction

import numpy as np
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.model_selection import StratifiedKFold, cross_val_score


def encode_residue_features(
    residue_sequence: str,
    b_factors: np.ndarray,
    conservation_scores: np.ndarray,
    coevolution_max_per_residue: np.ndarray,
    window_size: int = 5,
) -> np.ndarray:
    """
    Encode per-residue features for allosteric site prediction.

    Feature set per residue (window-averaged where applicable):
    - amino acid one-hot (20)
    - normalised B-factor (crystallographic flexibility indicator) (1)
    - evolutionary conservation score (1)
    - max mutual information with other residues (coevolution) (1)
    - window-mean physicochemical properties: hydrophobicity, charge, size (3)

    Args:
        residue_sequence: Single-letter amino acid sequence
        b_factors: Per-residue Cα B-factors (or RMSF from MD)
        conservation_scores: Per-residue conservation (0=variable, 1=conserved)
        coevolution_max_per_residue: Max MI coupling for each residue vs. any other
        window_size: Half-window for local averaging (2*window_size+1 neighbours)

    Returns:
        Feature matrix of shape (N_residues, 26).
    """
    AA = list("ACDEFGHIKLMNPQRSTVWY")
    aa_to_idx = {aa: i for i, aa in enumerate(AA)}

    HYDROPHOBICITY = {"A":0.62,"C":0.29,"D":-0.90,"E":-0.74,"F":1.19,"G":0.48,"H":-0.40,
                      "I":1.38,"K":-1.50,"L":1.06,"M":0.64,"N":-0.78,"P":0.12,"Q":-0.85,
                      "R":-2.53,"S":-0.18,"T":-0.05,"V":1.08,"W":0.81,"Y":0.26}
    CHARGE = {"D":-1,"E":-1,"K":1,"R":1,"H":0.1}
    SIZE = {"G":1,"A":2,"S":2,"C":2,"T":3,"V":3,"D":3,"P":3,"N":3,"I":4,
            "L":4,"E":4,"Q":4,"M":4,"H":5,"K":5,"F":6,"R":6,"Y":6,"W":7}

    N = len(residue_sequence)
    b_nrm = (b_factors - b_factors.min()) / (b_factors.ptp() + 1e-9)
    features = []

    for i, aa in enumerate(residue_sequence):
        onehot = [1 if aa == a else 0 for a in AA]
        win = residue_sequence[max(0,i-window_size):i+window_size+1]
        hydro = np.mean([HYDROPHOBICITY.get(r, 0) for r in win])
        charge = np.mean([CHARGE.get(r, 0) for r in win])
        size = np.mean([SIZE.get(r, 4) for r in win])
        row = onehot + [b_nrm[i], conservation_scores[i], coevolution_max_per_residue[i], hydro, charge, size]
        features.append(row)

    return np.array(features, dtype=np.float32)


def train_allosteric_classifier(
    X: np.ndarray,
    y: np.ndarray,
) -> dict:
    """
    Train a residue-level allosteric site classifier on ASD/CASBench labels.

    Uses StratifiedKFold CV because allosteric residues are typically ~5%
    of all residues — heavy class imbalance requires balanced scoring.
    """
    model = GradientBoostingClassifier(
        n_estimators=300, max_depth=3, learning_rate=0.08,
        subsample=0.8, min_samples_leaf=10, random_state=42,
    )
    cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    auc_scores = cross_val_score(model, X, y, cv=cv, scoring="roc_auc", n_jobs=-1)
    model.fit(X, y)
    return {
        "model": model,
        "cv_roc_auc": round(auc_scores.mean(), 3),
        "cv_roc_auc_std": round(auc_scores.std(), 3),
        "note": "Trained on ASD 3.0 + CASBench allosteric site labels",
    }

Key Takeaways

• Allosteric selectivity arises because allosteric pockets diverge far more in sequence than the highly-conserved orthosteric sites. This makes allosteric targeting a powerful strategy for kinases, GPCRs, and nuclear receptors where achieving orthosteric selectivity is notoriously difficult.
• Protein communication network analysis (betweenness centrality on Cα contact maps) identifies allosteric "bridge" residues that mediate long-range signalling. High-betweenness residues outside the active site are primary allosteric candidate locations.
• Cryptic pocket detection requires conformational ensemble sampling — apo X-ray structures alone are insufficient. Unbiased MD, metadynamics, or cosolvent MD are standard approaches to reveal transiently open cryptic sites sufficient for small-molecule binding.
• Coevolution signals (mutual information from multiple sequence alignments) provide evolutionary evidence for allosteric coupling — residue pairs with high MI that are distant in structure are strong allosteric communication candidates.

Try it live

Explore a related model hands-on: open Molecule Screening Engine — Similarity Search Live and change its parameters while it runs. Everything stays in your browser.