Hit-to-Lead Optimisation with ML

A hit compound typically has modest potency (IC50 1–10 µM), marginal selectivity, and poor ADMET properties. The goal of hit-to-lead is to improve potency (typically to <100 nM), establish selectivity, and confirm a tractable SAR. ML models trained on the growing primary SAR dataset can predict which structural modifications are most likely to improve a specific property — reducing the synthesis-test cycle from weeks to days.

Matched Molecular Pair Analysis

Matched Molecular Pairs (MMPs) are pairs of compounds that differ by a single well-defined structural transformation (a "cut"). Aggregating MMP data across a series reveals reliable SAR rules: "replacing H with F at position X consistently increases potency by 0.5 log units."

from rdkit import Chem
from rdkit.Chem import rdMMPFrags
import pandas as pd
from collections import defaultdict


def extract_mmp_transformations(
    smiles_list: list[str],
    activity_values: list[float],
    compound_ids: list[str],
    min_pair_count: int = 3,
    max_heavy_atom_ratio: float = 0.33,
) -> pd.DataFrame:
    """
    Extract Matched Molecular Pair transformations from a compound series.

    Args:
        smiles_list: Canonical SMILES for each compound
        activity_values: pIC50 or pKi values (log scale; higher = more potent)
        compound_ids: Compound identifiers
        min_pair_count: Minimum number of pairs to report a transformation rule
        max_heavy_atom_ratio: Max fraction of heavy atoms in the variable fragment

    Returns:
        DataFrame of transformation rules with mean delta_activity and pair count.
    """
    mols = [Chem.MolFromSmiles(s) for s in smiles_list]
    valid = [(m, a, cid) for m, a, cid in zip(mols, activity_values, compound_ids)
             if m is not None]

    # Use RDKit MMP fragmentation
    fragment_data = []
    for mol, activity, cid in valid:
        n_heavy = mol.GetNumHeavyAtoms()
        try:
            frags = rdMMPFrags.fragmentMol(mol, maxCuts=1)
            for frag in frags:
                core_smi = Chem.MolToSmiles(frag.core) if frag.core else None
                var_smi = Chem.MolToSmiles(frag.variable) if frag.variable else None
                if core_smi and var_smi:
                    var_mol = Chem.MolFromSmiles(var_smi)
                    if var_mol and var_mol.GetNumHeavyAtoms() / n_heavy <= max_heavy_atom_ratio:
                        fragment_data.append({
                            "compound_id": cid,
                            "core": core_smi,
                            "variable": var_smi,
                            "activity": activity,
                        })
        except Exception:
            continue

    df = pd.DataFrame(fragment_data)
    if df.empty:
        return pd.DataFrame()

    # Group by core, find pairs with same core but different variable
    transformation_records = []
    for core_smi, group in df.groupby("core"):
        if len(group) < 2:
            continue
        rows = group.to_dict("records")
        for i in range(len(rows)):
            for j in range(i + 1, len(rows)):
                a, b = rows[i], rows[j]
                if a["variable"] == b["variable"]:
                    continue
                delta = b["activity"] - a["activity"]
                key = (a["variable"], b["variable"]) if a["variable"] < b["variable"] else (b["variable"], a["variable"])
                transformation_records.append({
                    "from_fragment": key[0],
                    "to_fragment": key[1],
                    "delta_activity": abs(delta),
                    "direction": "+" if delta > 0 else "-",
                    "core": core_smi,
                    "compound_pair": (a["compound_id"], b["compound_id"]),
                })

    if not transformation_records:
        return pd.DataFrame()

    result_df = pd.DataFrame(transformation_records)
    rule_summary = (
        result_df.groupby(["from_fragment", "to_fragment"])
        .agg(
            mean_delta_activity=("delta_activity", "mean"),
            pair_count=("delta_activity", "count"),
            most_common_direction=("direction", lambda x: x.mode().iloc[0] if len(x) > 0 else "?"),
        )
        .reset_index()
        .query(f"pair_count >= {min_pair_count}")
        .sort_values("mean_delta_activity", ascending=False)
    )

    return rule_summary

Activity Cliff Detection

Activity cliffs are pairs of structurally similar compounds with large potency differences. They indicate highly sensitive SAR regions that are both an opportunity (small change = large gain) and a risk (fragile SAR = poor optimisability).

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


def compute_activity_cliffs(
    smiles_list: list[str],
    activities: list[float],
    compound_ids: list[str],
    similarity_threshold: float = 0.85,
    activity_cliff_threshold: float = 2.0,
    radius: int = 2,
    nbits: int = 2048,
) -> pd.DataFrame:
    """
    Identify activity cliffs: high-similarity pairs with high activity difference.

    Activity cliff criterion: Tanimoto(ECFP4) >= similarity_threshold
    AND |ΔpIC50| >= activity_cliff_threshold (default 2 orders of magnitude)

    Args:
        smiles_list: SMILES for each compound
        activities: pIC50 values
        compound_ids: Compound identifiers
        similarity_threshold: Minimum Tanimoto similarity to be an MMP candidate
        activity_cliff_threshold: Minimum |ΔpIC50| to flag as cliff

    Returns:
        DataFrame of activity cliff pairs sorted by severity score.
    """
    mols = [Chem.MolFromSmiles(s) for s in smiles_list]
    fps = []
    valid_idx = []
    for i, mol in enumerate(mols):
        if mol is not None:
            fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius, nbits)
            fps.append(fp)
            valid_idx.append(i)

    act = [activities[i] for i in valid_idx]
    cids = [compound_ids[i] for i in valid_idx]
    n = len(fps)

    cliff_pairs = []
    for i in range(n):
        for j in range(i + 1, n):
            sim = DataStructs.TanimotoSimilarity(fps[i], fps[j])
            if sim < similarity_threshold:
                continue
            delta = abs(act[i] - act[j])
            if delta < activity_cliff_threshold:
                continue
            # Severity: high similarity + high delta = most dangerous cliff
            severity = sim * delta
            cliff_pairs.append({
                "compound_a": cids[i],
                "compound_b": cids[j],
                "activity_a": act[i],
                "activity_b": act[j],
                "delta_activity": round(delta, 2),
                "tanimoto_similarity": round(sim, 3),
                "cliff_severity": round(severity, 3),
            })

    if not cliff_pairs:
        return pd.DataFrame()

    return (
        pd.DataFrame(cliff_pairs)
        .sort_values("cliff_severity", ascending=False)
        .reset_index(drop=True)
    )

Iterative SAR Model for Synthesis Prioritisation

from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import LeaveOneOut, cross_val_predict
from sklearn.metrics import mean_absolute_error
import numpy as np


def build_sar_model(
    smiles_list: list[str],
    activities: list[float],
    radius: int = 2,
    nbits: int = 2048,
) -> dict:
    """
    Train a Random Forest SAR model on ECFP4 fingerprints.

    Uses Leave-One-Out CV for small datasets (typical HtL series: 30-150 cpds).

    Args:
        smiles_list: Canonical SMILES
        activities: pIC50 values
        radius: Morgan fingerprint radius
        nbits: Fingerprint bit length

    Returns:
        Dict with trained model, LOO-CV MAE, and feature importance summary.
    """
    mols = [Chem.MolFromSmiles(s) for s in smiles_list]
    fps = []
    valid_activities = []
    for mol, act in zip(mols, activities):
        if mol is not None:
            fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius, nbits)
            fps.append(list(fp))
            valid_activities.append(act)

    X = np.array(fps)
    y = np.array(valid_activities)

    model = RandomForestRegressor(
        n_estimators=500,
        max_features="sqrt",
        min_samples_leaf=2,
        random_state=42,
        n_jobs=-1,
    )

    if len(y) <= 30:
        cv_preds = cross_val_predict(model, X, y, cv=LeaveOneOut(), n_jobs=-1)
        cv_type = "LOO-CV"
    else:
        from sklearn.model_selection import KFold
        cv_preds = cross_val_predict(model, X, y, cv=KFold(5, shuffle=True, random_state=42), n_jobs=-1)
        cv_type = "5-Fold CV"

    mae = mean_absolute_error(y, cv_preds)
    model.fit(X, y)

    return {
        "model": model,
        "cv_type": cv_type,
        "cv_mae_pIC50": round(mae, 3),
        "n_training_compounds": len(y),
        "description": "Predict new compound activities; retrain after each round.",
    }


def predict_virtual_library(
    model,
    virtual_smiles: list[str],
    top_n: int = 50,
    radius: int = 2,
    nbits: int = 2048,
) -> pd.DataFrame:
    """
    Score a virtual library and return the top_n predicted actives.

    Args:
        model: Trained sklearn model (from build_sar_model)
        virtual_smiles: SMILES list of virtual compounds to score
        top_n: Number of top compounds to return

    Returns:
        DataFrame of top-ranked virtual compounds with predicted pIC50.
    """
    fps, valid_smiles = [], []
    for smi in virtual_smiles:
        mol = Chem.MolFromSmiles(smi)
        if mol:
            fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius, nbits)
            fps.append(list(fp))
            valid_smiles.append(smi)

    if not fps:
        return pd.DataFrame()

    X = np.array(fps)
    preds = model.predict(X)

    result = pd.DataFrame({"smiles": valid_smiles, "predicted_pIC50": preds})
    return result.nlargest(top_n, "predicted_pIC50").reset_index(drop=True)

Hit-to-Lead Decision Criteria

Key Takeaways

• MMP analysis transforms raw SAR data into generalised rules. A transformation with ≥5 pairs and consistent direction is actionable; fewer pairs may be noise. Automate extraction after each data dump.
• Activity cliffs are high-value alerts. They reveal positions where single atom changes make or break activity — immediately flag these for chemical series risk assessment.
• RF models on ECFP4 with LOO-CV are the right choice for small HtL series. CV MAE below 0.5 log units is generally predictive enough to rank next-round synthesis suggestions.
• Track ligand efficiency metrics alongside raw potency. A potency gain achieved by adding MW without LE improvement is a sign of non-specific binding, not real SAR progress.

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.