Scaffold Morphing and Bioisosteric Replacement

A series built around a single scaffold may hit IP wall (competitor patent), face synthetic intractability for key analogues, or show class-wide liability (e.g., photoreactivity of a benzimidazole). Scaffold morphing generates structurally distinct but pharmacophore-equivalent compounds — opening new IP space while maintaining binding interactions.

Modern approaches combine rule-based bioisostere databases with generative ML models. The ML models learn to make pharmacophore-conservative transformations from large chemical databases, producing suggestions beyond the finite bioisostere lookup tables.

Murcko Scaffold Analysis and Scaffold Tree

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


def get_murcko_scaffold(smiles: str, generic: bool = False) -> str | None:
    """
    Get the Murcko scaffold (ring systems + linkers) for a compound.

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

    Returns:
        Scaffold SMILES or None if parsing fails.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None
    scaffold = MurckoScaffold.GetScaffoldForMol(mol)
    if generic:
        scaffold = MurckoScaffold.MakeScaffoldGeneric(scaffold)
    return Chem.MolToSmiles(scaffold)


def scaffold_series_analysis(
    smiles_list: list[str],
    compound_ids: list[str],
    activities: list[float] | None = None,
) -> pd.DataFrame:
    """
    Analyse scaffold diversity and potency distribution in a compound series.

    Args:
        smiles_list: Compound SMILES
        compound_ids: Compound identifiers
        activities: Optional pIC50 values for per-scaffold activity summary

    Returns:
        DataFrame with scaffold SMILES, member count, and mean activity.
    """
    records = []
    for smi, cid in zip(smiles_list, compound_ids):
        scaffold = get_murcko_scaffold(smi)
        records.append({"compound_id": cid, "smiles": smi, "scaffold": scaffold})

    df = pd.DataFrame(records)
    if activities:
        df["activity"] = activities

    if activities:
        summary = (
            df.groupby("scaffold")
            .agg(
                n_compounds=("compound_id", "count"),
                mean_activity=("activity", "mean"),
                max_activity=("activity", "max"),
                compound_ids=("compound_id", list),
            )
            .sort_values("n_compounds", ascending=False)
            .reset_index()
        )
    else:
        summary = (
            df.groupby("scaffold")
            .agg(n_compounds=("compound_id", "count"), compound_ids=("compound_id", list))
            .sort_values("n_compounds", ascending=False)
            .reset_index()
        )

    return summary

Classical Bioisostere Lookups

from typing import Optional

# Curated classical bioisostere replacements
# Based on medicinal chemistry knowledge and BIOSTER database patterns
BIOISOSTERE_TABLE: dict[str, list[str]] = {
    # Carboxylic acid bioisosteres
    "C(=O)O":        ["C(=O)[NH]S(=O)(=O)", "c1cnco1", "C(=O)NS(=O)(=O)C",
                       "C(=O)NHC(=O)", "[C@@H](O)C(=O)O"],
    # Amide bioisosteres
    "C(=O)N":        ["C(=S)N", "C(=O)O", "C(=N)N", "c1ocnc1", "C=NN"],
    # Ester bioisosteres
    "C(=O)OC":       ["C(=O)NC", "C(=O)SC", "c1ccoc1"],
    # Phenyl ring bioisosteres
    "c1ccccc1":      ["c1ccncc1", "c1cccn1", "c1cncs1", "C1CCCCC1", "c1ccco1"],
    # Chloro substitution
    "Cl":            ["F", "CF3", "C#N", "Br"],
    # Methyl group
    "C":             ["CF3", "C#N", "cyclopropyl"],
    # Hydroxyl bioisosteres
    "O":             ["S", "NH", "F"],
    # Ketone bioisosteres
    "C(=O)C":        ["C(=S)C", "C(=NH)C", "c1ccoc1"],
    # Thioether bioisosteres
    "CSC":           ["COC", "CNC", "C-C"],
}


def find_bioisosteres(
    smiles: str,
    max_replacements: int = 5,
) -> list[dict]:
    """
    Find bioisosteric replacements for substructures in a compound.

    Args:
        smiles: Input compound SMILES
        max_replacements: Maximum number of replacements to return per substructure

    Returns:
        List of {original, replacement, product_smiles} dicts.
    """
    from rdkit.Chem import AllChem, RWMol
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return []

    results = []
    for query_smi, replacements in BIOISOSTERE_TABLE.items():
        query = Chem.MolFromSmarts(query_smi) if ":" not in query_smi else Chem.MolFromSmarts(query_smi)
        if query is None:
            continue
        if not mol.HasSubstructMatch(query):
            continue

        for replacement_smi in replacements[:max_replacements]:
            try:
                rxn_smarts = f"[{query_smi}:1]>>[{replacement_smi}:1]"
                # Simple SMILES replacement via AllChem reaction
                product_smi = smiles.replace(query_smi, replacement_smi, 1)
                product_mol = Chem.MolFromSmiles(product_smi)
                if product_mol and product_smi != smiles:
                    results.append({
                        "original_fragment": query_smi,
                        "replacement_fragment": replacement_smi,
                        "product_smiles": Chem.MolToSmiles(product_mol),
                    })
            except Exception:
                continue

    return results


def property_filter_bioisosteres(
    candidates: list[dict],
    original_smiles: str,
    max_mw_change: float = 50.0,
    max_logp_change: float = 1.5,
) -> list[dict]:
    """
    Filter bioisostere candidates to those with similar physicochemical properties.

    Args:
        candidates: Output of find_bioisosteres()
        original_smiles: Reference compound SMILES
        max_mw_change: Maximum MW change allowed (Da)
        max_logp_change: Maximum cLogP change allowed

    Returns:
        Filtered and property-annotated candidate list.
    """
    from rdkit.Chem import Descriptors
    orig_mol = Chem.MolFromSmiles(original_smiles)
    if orig_mol is None:
        return candidates

    orig_mw = Descriptors.MolWt(orig_mol)
    orig_logp = Descriptors.MolLogP(orig_mol)

    filtered = []
    for cand in candidates:
        prod_mol = Chem.MolFromSmiles(cand["product_smiles"])
        if prod_mol is None:
            continue
        delta_mw = abs(Descriptors.MolWt(prod_mol) - orig_mw)
        delta_logp = abs(Descriptors.MolLogP(prod_mol) - orig_logp)
        if delta_mw <= max_mw_change and delta_logp <= max_logp_change:
            cand["delta_mw"] = round(delta_mw, 1)
            cand["delta_logp"] = round(delta_logp, 2)
            filtered.append(cand)

    return filtered

Common Bioisostere Pairs in Drug Discovery

Key Takeaways

• Bioisosteric replacement is most effective when the original group has a known liability: metabolic hotspot, hERG pharmacophore element, IP-blocked fragment, or flat SAR indicative of non-specific binding.
• Classical bioisostere lookup tables (BIOSTER, SwissBioisostere) should be the first step — they encode expert knowledge in well-validated transformation pairs with predictable property changes.
• ML generative approaches (REINVENT Mol2Mol, fragment VAE) are best for exploring beyond classical replacements, especially in scaffold-hopping where the solution space is too large for manual enumeration.
• Always apply physicochemical filters after generating bioisostere candidates: MW change <50 Da and cLogP change <1.5 are practical guardrails that prevent inadvertently generating ADMET-problematic structures.

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.