Scaffold Analysis and Lead Optimization

Once a hit compound is confirmed in a biochemical assay, the medicinal chemistry team moves into lead optimization — systematically modifying the structure to improve potency, selectivity, solubility, metabolic stability, and safety. Computational analysis of chemical series using scaffold frameworks, structure-activity relationships (SAR), and matched molecular pair (MMP) analysis provides the data-driven backbone for these decisions.

This guide covers the key computational tools for scaffold analysis, how to build SAR tables and activity cliff maps, how to perform MMP analysis to quantify the effect of individual substituent changes, and how to use these insights for scaffold hopping.

Bemis–Murcko Scaffolds

The Bemis–Murcko scaffold decomposition, introduced in 1996, is the standard framework for classifying drug-like molecules into structural families. It reduces a molecule to its:

• Ring systems: all ring atoms and the bonds between them
• Linkers: chains connecting ring systems
• Side chains: substituents hanging off ring atoms

The scaffold is the ring systems plus linkers, stripped of all side chains. Two molecules with the same scaffold belong to the same chemical series even if their peripheral substitution patterns differ dramatically.

Computing Scaffolds with RDKit

from rdkit import Chem
from rdkit.Chem.Scaffolds import MurckoScaffold
import pandas as pd
from collections import Counter

def get_murcko_scaffold(smiles: str, generic: bool = False) -> str:
    """
    Return the Bemis-Murcko scaffold SMILES for a molecule.

    Args:
        smiles:  Input SMILES
        generic: If True, return the generic scaffold (all C, single bonds)

    Returns:
        Scaffold SMILES string, or empty string if molecule is invalid.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return ""
    if generic:
        scaffold = MurckoScaffold.MakeScaffoldGeneric(
            MurckoScaffold.GetScaffoldForMol(mol)
        )
    else:
        scaffold = MurckoScaffold.GetScaffoldForMol(mol)
    return Chem.MolToSmiles(scaffold)


# Example: analyse scaffold distribution in a ChEMBL activity set
df = pd.read_csv("chembl_activities.csv")  # columns: smiles, pIC50
df["scaffold"] = df["smiles"].apply(get_murcko_scaffold)

# Count scaffold frequency
scaffold_counts = Counter(df["scaffold"])
top_scaffolds = scaffold_counts.most_common(10)
print("Top 10 scaffolds by compound count:")
for scaffold, count in top_scaffolds:
    print(f"  {scaffold}: {count} compounds")

# Average potency per scaffold
scaffold_sar = df.groupby("scaffold")["pIC50"].agg(["mean", "std", "count"])
scaffold_sar = scaffold_sar[scaffold_sar["count"] >= 5].sort_values("mean", ascending=False)
print(scaffold_sar.head(10))

Activity Cliffs

An activity cliff is a pair of structurally similar compounds with a large difference in biological activity (typically ≥2 log units in pIC50). Activity cliffs are simultaneously the most informative and most dangerous data points in a chemical series:

• Informative: they pinpoint critical structural feature for binding
• Dangerous: ML models often fail to extrapolate across them because structural similarity does not imply activity similarity

The standard metric is the Structure-Activity Landscape Index (SALI):

SALI(i,j) = |A_i - A_j| / (1 - sim(i,j))

Where A is activity value and sim is Tanimoto similarity. High SALI scores identify the steepest activity cliffs.

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

def compute_tanimoto_matrix(smiles_list: list, radius: int = 2, n_bits: int = 2048) -> np.ndarray:
    """Compute pairwise Tanimoto similarity matrix."""
    fps = []
    for smi in smiles_list:
        mol = Chem.MolFromSmiles(smi)
        if mol is not None:
            fps.append(AllChem.GetMorganFingerprintAsBitVect(mol, radius, nBits=n_bits))
        else:
            fps.append(None)

    n = len(fps)
    matrix = np.zeros((n, n))
    for i in range(n):
        for j in range(i + 1, n):
            if fps[i] is not None and fps[j] is not None:
                sim = DataStructs.TanimotoSimilarity(fps[i], fps[j])
                matrix[i, j] = sim
                matrix[j, i] = sim
    np.fill_diagonal(matrix, 1.0)
    return matrix


def find_activity_cliffs(
    df: pd.DataFrame,
    smiles_col: str = "smiles",
    activity_col: str = "pIC50",
    sim_threshold: float = 0.6,
    activity_threshold: float = 2.0
) -> pd.DataFrame:
    """
    Find activity cliff pairs: Tanimoto >= sim_threshold but
    |pIC50 difference| >= activity_threshold.
    """
    smiles = df[smiles_col].tolist()
    activities = df[activity_col].values
    sim_matrix = compute_tanimoto_matrix(smiles)

    cliffs = []
    n = len(smiles)
    for i in range(n):
        for j in range(i + 1, n):
            sim = sim_matrix[i, j]
            act_diff = abs(activities[i] - activities[j])
            if sim >= sim_threshold and act_diff >= activity_threshold:
                sali = act_diff / (1 - sim + 1e-9)
                cliffs.append({
                    "smiles_1": smiles[i],
                    "smiles_2": smiles[j],
                    "activity_1": activities[i],
                    "activity_2": activities[j],
                    "tanimoto": round(sim, 3),
                    "activity_diff": round(act_diff, 2),
                    "SALI": round(sali, 2),
                })
    return pd.DataFrame(cliffs).sort_values("SALI", ascending=False)

Matched Molecular Pair Analysis

Matched Molecular Pair (MMP) analysis transforms the question "how does this substituent change affect activity?" into a statistical answer. An MMP is defined as two compounds that differ by exactly one structural transformation — the replacement of one fragment (context) with another (transformation). By collecting hundreds or thousands of such pairs from a dataset, you build a data-driven transformation dictionary.

Core Concepts

• Context: the shared scaffold portion (at least one ring, typically)
• Transformation: the exchanged fragment pair, e.g. H → F, CH3 → CF3, phenyl → pyridyl
• Property change Δ: the measured difference in any property (pIC50, logD, clearance)

Common insights from MMP analysis in drug discovery:

MMP Analysis with mmpdb

"""
mmpdb is the standard open-source tool for MMP analysis.
Install: pip install mmpdb

Workflow:
1. Fragment a SMILES dataset into MMP pairs
2. Load an activity CSV
3. Query property changes for specific transformations
"""
import subprocess
import pandas as pd

# Step 1: fragment the compound set
subprocess.run([
    "mmpdb", "fragment",
    "compounds.smi",         # tab-separated: SMILES  ID
    "-o", "fragments.frag"
], check=True)

# Step 2: index the fragments into an MMP database
subprocess.run([
    "mmpdb", "index",
    "fragments.frag",
    "-o", "compounds.mmpdb"
], check=True)

# Step 3: load activity data
# CSV with columns: ID, pIC50, logD, clearance
subprocess.run([
    "mmpdb", "loadprops",
    "-i", "activities.csv",
    "compounds.mmpdb"
], check=True)

# Step 4: query for all transformations affecting pIC50
result = subprocess.run([
    "mmpdb", "transform",
    "--smiles", "c1ccc(NC(=O)c2ccccc2)cc1",  # query compound SMILES
    "--property", "pIC50",
    "compounds.mmpdb"
], capture_output=True, text=True, check=True)

print(result.stdout[:2000])  # First few predicted transformations

Implementing Simple MMP Detection in Python

from rdkit import Chem
from rdkit.Chem import rdMMPA  # Matched Molecular Pair Analysis module
import pandas as pd

def extract_mmps(smiles_list: list, ids: list) -> pd.DataFrame:
    """
    Use RDKit's rdMMPA to fragment molecules and identify MMP pairs.
    Returns a dataframe of (id, smiles, core, r_group) fragments.
    """
    records = []
    for smi, mol_id in zip(smiles_list, ids):
        mol = Chem.MolFromSmiles(smi)
        if mol is None:
            continue
        # FragmentMol returns (core, r_group) pairs
        frags = rdMMPA.FragmentMol(mol, maxCuts=1, resultsAsMols=False)
        for core_smi, rgroup_smi in frags:
            records.append({
                "id": mol_id,
                "smiles": smi,
                "core": core_smi,
                "rgroup": rgroup_smi,
            })
    return pd.DataFrame(records)


# Find pairs with the same core but different R-groups
def find_mmp_pairs(df: pd.DataFrame) -> pd.DataFrame:
    """Group by core and return all pairs sharing a scaffold context."""
    pairs = []
    for core, group in df.groupby("core"):
        if len(group) < 2:
            continue
        rows = group.reset_index(drop=True)
        for i in range(len(rows)):
            for j in range(i + 1, len(rows)):
                if rows.loc[i, "rgroup"] != rows.loc[j, "rgroup"]:
                    pairs.append({
                        "core": core,
                        "id_1": rows.loc[i, "id"],
                        "id_2": rows.loc[j, "id"],
                        "smiles_1": rows.loc[i, "smiles"],
                        "smiles_2": rows.loc[j, "smiles"],
                        "rgroup_1": rows.loc[i, "rgroup"],
                        "rgroup_2": rows.loc[j, "rgroup"],
                    })
    return pd.DataFrame(pairs)

Scaffold Hopping

Scaffold hopping is the replacement of the core ring system of a compound while retaining its binding pharmacophore. Reasons to hop scaffolds include avoiding competitor IP, improving aqueous solubility, reducing CYP inhibition from the original scaffold, or escaping metabolic soft spots within the rings.

Strategies for Scaffold Hopping

• Bioisosteric replacement: replace a ring with another that has similar electronic and steric properties (benzene → thiophene, pyridine, furan)
• Ring opening/closing: macro-rings opened to acyclic linkers or vice versa
• Shape-based virtual screening: 3D overlay of candidate scaffolds to the reference compound
• ML-guided generation: constrained molecular generators (REINVENT scaffold memory) that keep pharmacophore atoms fixed
• Fragment-based approaches: replace individual ring atoms (carbon → nitrogen, bioisostere tables)

Similarity-Based Scaffold Hop Search

from rdkit import Chem, DataStructs
from rdkit.Chem import AllChem
from rdkit.Chem.Scaffolds import MurckoScaffold
import pandas as pd

def scaffold_hop_candidates(
    query_smiles: str,
    library_df: pd.DataFrame,
    smiles_col: str = "smiles",
    top_k: int = 20,
    fp_sim_min: float = 0.3,
    fp_sim_max: float = 0.6,
) -> pd.DataFrame:
    """
    Find library compounds that are moderately similar to the query
    (same pharmacophore region) but have a different Murcko scaffold.

    Low-to-moderate Tanimoto + different scaffold = scaffold hop candidate.
    """
    query_mol = Chem.MolFromSmiles(query_smiles)
    query_fp = AllChem.GetMorganFingerprintAsBitVect(query_mol, 2, 2048)
    query_scaffold = MurckoScaffold.MurckoScaffoldSmiles(mol=query_mol)

    results = []
    for _, row in library_df.iterrows():
        mol = Chem.MolFromSmiles(row[smiles_col])
        if mol is None:
            continue
        fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2, 2048)
        sim = DataStructs.TanimotoSimilarity(query_fp, fp)
        if not (fp_sim_min <= sim <= fp_sim_max):
            continue
        scaf = MurckoScaffold.MurckoScaffoldSmiles(mol=mol)
        if scaf == query_scaffold:
            continue  # Same scaffold — not a hop
        results.append({**row.to_dict(), "tanimoto": round(sim, 3), "scaffold": scaf})

    return pd.DataFrame(results).sort_values("tanimoto", ascending=False).head(top_k)

Building a SAR Table

A SAR (Structure-Activity Relationship) table organizes compounds in a chemical series to visualize the effect of each substitution position on activity. Automating SAR table generation from ChEMBL data is a common workflow at the start of a project.

import pandas as pd
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem.Scaffolds import MurckoScaffold

def build_sar_table(df: pd.DataFrame, scaffold_smi: str) -> pd.DataFrame:
    """
    For a given scaffold SMILES, extract all compounds in df that share it
    and compute their R-group substituents at each attachment point.

    Returns a dataframe with columns:
        smiles, pIC50, scaffold, rgroup_summary
    """
    target_scaffold = Chem.MolToSmiles(Chem.MolFromSmiles(scaffold_smi))
    rows = []
    for _, row in df.iterrows():
        mol = Chem.MolFromSmiles(row["smiles"])
        if mol is None:
            continue
        scaf_smi = Chem.MolToSmiles(MurckoScaffold.GetScaffoldForMol(mol))
        if scaf_smi != target_scaffold:
            continue
        rows.append(row)

    sar_df = pd.DataFrame(rows)
    sar_df = sar_df.sort_values("pIC50", ascending=False)
    return sar_df


# Example usage
df = pd.read_csv("series_activities.csv")  # smiles, pIC50, compound_id
core = "c1ccc2ncccc2c1"                    # quinoline scaffold SMILES
sar = build_sar_table(df, core)
print(sar[["compound_id", "smiles", "pIC50"]].to_string(index=False))

Activity Cliff Awareness in ML Models

Activity cliffs are the main failure mode for fingerprint-based ML models. When the training set contains a cliff pair and the model assigns them similar latent representations, predictions for cliff-related compounds will be systematically wrong. Mitigation strategies include:

• Scaffold-aware splits: never let cliff pairs end up split between train and test. Use scaffold splits (Butina clustering by scaffold membership) to ensure each scaffold appears in only one partition.
• Activity cliff flagging: precompute SALI scores and flag high-SALI pairs. Exclude them from test sets or treat them as a held-out cliff benchmark.
• Graph neural networks: GNNs operating on full molecular graphs are somewhat more cliff-aware than fingerprints because they learn continuous structural representations. But they are not immune.
• Uncertainty estimation: conformal prediction or MC Dropout can quantify model uncertainty. High uncertainty on a prediction is a reliable signal to not trust the model's output near a cliff.

from rdkit import Chem
from rdkit.Chem.Scaffolds import MurckoScaffold
from sklearn.model_selection import GroupShuffleSplit
import pandas as pd
import numpy as np

def scaffold_split(
    df: pd.DataFrame,
    test_size: float = 0.2,
    random_state: int = 42
) -> tuple:
    """
    Perform scaffold-based train/test split.
    Compounds with the same Murcko scaffold always go to the same split.
    """
    scaffolds = []
    for smi in df["smiles"]:
        mol = Chem.MolFromSmiles(smi)
        if mol is None:
            scaffolds.append("")
        else:
            scaffolds.append(
                Chem.MolToSmiles(MurckoScaffold.GetScaffoldForMol(mol))
            )

    groups = pd.Categorical(scaffolds).codes
    gss = GroupShuffleSplit(n_splits=1, test_size=test_size, random_state=random_state)
    train_idx, test_idx = next(gss.split(df, groups=groups))
    return df.iloc[train_idx], df.iloc[test_idx]


df = pd.read_csv("activities.csv")
train_df, test_df = scaffold_split(df)
print(f"Train: {len(train_df)} compounds | Test: {len(test_df)} compounds")

# Verify no scaffold overlap
train_scaffolds = set(train_df["smiles"].apply(
    lambda s: Chem.MolToSmiles(MurckoScaffold.GetScaffoldForMol(Chem.MolFromSmiles(s)))
    if Chem.MolFromSmiles(s) else ""
))
test_scaffolds = set(test_df["smiles"].apply(
    lambda s: Chem.MolToSmiles(MurckoScaffold.GetScaffoldForMol(Chem.MolFromSmiles(s)))
    if Chem.MolFromSmiles(s) else ""
))
overlap = train_scaffolds & test_scaffolds
print(f"Scaffold overlap: {len(overlap)} scaffolds")  # Should be 0

Key Takeaways

• Bemis–Murcko scaffolds are the standard way to group drug-like molecules into chemical series; compute them with RDKit's MurckoScaffold module.
• Activity cliffs (structurally similar, very different potency) are the primary failure mode for fingerprint ML. Quantify them with SALI scores before modeling.
• Matched Molecular Pair analysis turns your assay data into a transformation dictionary: "replacing H with F at this position raises pIC50 by 0.5 log units on average."
• Always use scaffold-aware train/test splits. Random splits inflate performance metrics by leaking scaffold information from train to test.
• Scaffold hopping is most productively guided by combining pharmacophore retention (shape/electrostatics match) with molecular generation or similarity-based library search.

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.