Metabolite Prediction and Metabolic Stability ML

CYP3A4, CYP2D6, CYP2C9, CYP2C19, and CYP1A2 account for >90% of phase I oxidative metabolism. Site-of-metabolism (SOM) prediction identifies which atoms on a drug molecule are most likely to be oxidised, allowing medicinal chemists to block metabolic soft spots with fluorine substitution, deuterium labelling, or bioisosteric replacement without sacrificing potency.

Site of Metabolism Prediction (Atom-Level)

import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score


def encode_atom_features_for_som(
    mol: Chem.Mol,
    radius: int = 2,
    n_bits: int = 1024,
) -> tuple[np.ndarray, list[int]]:
    """
    Encode per-atom features for site-of-metabolism prediction.

    For each non-hydrogen atom, generate:
    - Atom-centred Morgan fingerprint (local environment within radius bonds)
    - Atom-level physicochemical features: atomic number, degree, total H count,
      formal charge, is_aromatic, ring membership, partial charge index

    Returns:
        (feature_matrix, atom_indices): feature matrix (N_atoms × N_features)
        and corresponding atom indices.
    """
    atom_fps = []
    atom_idxs = []

    for atom in mol.GetAtoms():
        if atom.GetAtomicNum() == 1:  # skip explicit H
            continue

        # Atom-centred Morgan environment bitvector
        env = AllChem.GetMorganFingerprintAsBitVect(
            mol, radius, n_bits, fromAtoms=[atom.GetIdx()],
        )
        fp = list(env)

        # Atom-level descriptors
        a_feats = [
            atom.GetAtomicNum() / 16.0,
            atom.GetDegree() / 4.0,
            atom.GetTotalNumHs() / 4.0,
            atom.GetFormalCharge(),
            float(atom.GetIsAromatic()),
            float(atom.IsInRing()),
            float(atom.GetHybridization().name in ("SP2", "SP")),
        ]

        atom_fps.append(fp + a_feats)
        atom_idxs.append(atom.GetIdx())

    return np.array(atom_fps, dtype=np.float32), atom_idxs


def train_som_model(
    molecules: list[Chem.Mol],
    som_labels: list[list[int]],
    cyp_isoform: str = "CYP3A4",
) -> dict:
    """
    Train an atom-level SOM classifier for a specific CYP isoform.

    Args:
        molecules: List of RDKit Mol objects
        som_labels: Per-molecule list of atom indices confirmed as SOM sites
                    (from SMARTCyp or experimental NMR/MS data)
        cyp_isoform: CYP isoform label for record-keeping

    Returns:
        Dict with trained model and CV ROC-AUC.
    """
    X_all, y_all = [], []
    for mol, som_atoms in zip(molecules, som_labels):
        feats, idxs = encode_atom_features_for_som(mol)
        som_set = set(som_atoms)
        for feat, idx in zip(feats, idxs):
            X_all.append(feat)
            y_all.append(1 if idx in som_set else 0)

    X = np.array(X_all)
    y = np.array(y_all)

    model = RandomForestClassifier(
        n_estimators=400, class_weight="balanced",
        max_depth=None, min_samples_leaf=3, random_state=42, n_jobs=-1,
    )
    cv_roc = cross_val_score(model, X, y, cv=5, scoring="roc_auc", n_jobs=-1)
    model.fit(X, y)

    return {
        "model": model, "cyp_isoform": cyp_isoform,
        "cv_roc_auc": round(cv_roc.mean(), 3),
        "positive_rate": round(float(y.mean()), 3),
    }


def predict_som(
    mol: Chem.Mol,
    som_model,
    top_n: int = 3,
) -> list[tuple[int, float]]:
    """
    Predict top-N site(s) of metabolism for a molecule.

    Returns:
        List of (atom_idx, probability) sorted by probability descending.
    """
    feats, idxs = encode_atom_features_for_som(mol)
    if len(feats) == 0:
        return []
    probas = som_model.predict_proba(feats)[:, 1]
    ranked = sorted(zip(idxs, probas.tolist()), key=lambda x: x[1], reverse=True)
    return ranked[:top_n]

Microsomal Stability Classification

import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.preprocessing import LabelEncoder
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler


STABILITY_CLASSES = {"stable": 0, "moderate": 1, "labile": 2}
# Boundaries: t½ > 60 min (HLM) = stable; 30–60 min = moderate; < 30 min = labile


def build_metabolic_stability_features(smiles: str) -> np.ndarray | None:
    """
    Compute molecular features predictive of human liver microsomal (HLM) stability.

    Covers: ECFP4-based fingerprint + 12 physicochemical properties known to
    correlate with CYP substrate recognition and metabolic clearance.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None

    fp = list(AllChem.GetMorganFingerprintAsBitVect(mol, 2, 1024))
    physchem = [
        Descriptors.MolWt(mol),
        Descriptors.MolLogP(mol),
        Descriptors.TPSA(mol),
        Descriptors.NumRotatableBonds(mol),
        Descriptors.NumHAcceptors(mol),
        Descriptors.NumHDonors(mol),
        Descriptors.NumAromaticRings(mol),
        Descriptors.FractionCSP3(mol),
        Descriptors.HeavyAtomCount(mol),
        Descriptors.NumAliphaticCarbocycles(mol),
        Descriptors.NumHeteroatoms(mol),
        Descriptors.RingCount(mol),
    ]
    return np.array(fp + physchem, dtype=np.float32)


def train_metabolic_stability_model(
    smiles_list: list[str],
    stability_labels: list[str],
) -> dict:
    """
    Train 3-class metabolic stability classifier (stable/moderate/labile).

    Training data sources: ChEMBL ADMET, Lombardo et al. datasets, AstraZeneca HLM data.
    """
    X, y = [], []
    le = LabelEncoder()
    y_encoded = le.fit_transform(stability_labels)

    for smi, label in zip(smiles_list, y_encoded):
        feats = build_metabolic_stability_features(smi)
        if feats is not None:
            X.append(feats)
            y.append(label)

    X_arr = np.array(X)
    y_arr = np.array(y)

    from sklearn.model_selection import cross_val_score
    model = Pipeline([
        ("scaler", StandardScaler()),
        ("gbt", GradientBoostingClassifier(n_estimators=300, max_depth=4, learning_rate=0.05, random_state=42)),
    ])
    cv_acc = cross_val_score(model, X_arr, y_arr, cv=5, scoring="balanced_accuracy", n_jobs=-1)
    model.fit(X_arr, y_arr)

    return {
        "model": model, "label_encoder": le,
        "cv_balanced_accuracy": round(cv_acc.mean(), 3),
        "classes": list(le.classes_),
    }

Key CYP Isoform Substrate Characteristics

Key Takeaways

• Site-of-metabolism prediction is most accurate at the atom-centred fingerprint level — Morgan environment descriptors around each candidate atom capture the local electronic and steric context that governs CYP oxidation selectivity.
• Metabolic soft-spot blocking is most effective with fluorine substitution at predicted SOM positions (C–F bonds are much stronger than C–H bonds). Deuterium labelling (kinetic isotope effect) is also used for high-value candidates.
• Reactive metabolite formation (quinones, epoxides, Michael acceptors) is the primary driver of idiosyncratic hepatotoxicity. Structural alerts (thiophenes, anilines, furans, acetylenes) are high-risk bioisosteres to avoid or protect.
• Human liver microsomal (HLM) t½ correlates poorly with in vivo clearance for highly protein-bound compounds (fu,mic correction) or transporters. ML models trained on HLM data should be supplemented with hepatocyte and organ-of-clearance data for late-stage candidates.

Try it live

Explore a related model hands-on: open Cytochrome P450 Drug Metabolism and change its parameters while it runs. Everything stays in your browser.