Pharmacophore Modeling and 3D Virtual Screening

Pharmacophore models are particularly valuable for scaffold-hopping: two structurally dissimilar molecules can share the same pharmacophore and exhibit equivalent biological activity. IUPAC defines a pharmacophore as "the ensemble of steric and electronic features necessary to ensure optimal supramolecular interactions with a specific biological target and to trigger (or to block) its biological response."

The six canonical pharmacophoric feature types are: H-bond donor (HBD), H-bond acceptor (HBA), hydrophobic (HYD), positive ionisable (PI), negative ionisable (NI), and aromatic ring (AR). Exclusion volumes derived from the binding site complement the feature points to remove steric clashes.

Pharmacophoric Feature Extraction

from rdkit import Chem
from rdkit.Chem import AllChem, MolStandardize
from rdkit.Chem.Pharm3D import Pharmacophore, EmbedLib
from rdkit.Chem import ChemicalFeatures
from rdkit.Chem.rdMolChemicalFeatures import MolChemicalFeatures
import numpy as np
from dataclasses import dataclass


@dataclass
class PharmFeature:
    name: str       # HBD, HBA, HYD, PI, NI, AR
    coords: tuple[float, float, float]
    importance: float = 1.0  # weight in match scoring


def extract_pharmacophore_features(
    mol: Chem.Mol,
    conf_id: int = 0,
    feature_factory_path: str = None,
) -> list[PharmFeature]:
    """
    Extract pharmacophoric feature points from a 3D molecular conformer.

    Uses RDKit's built-in Feature Factory (BaseFeatures.fdef) which defines
    SMARTS patterns for HBD, HBA, Hydrophobe, PosIonizable, NegIonizable, Aromatic.

    Args:
        mol: RDKit Mol with at least one 3D conformer
        conf_id: Conformer index to use for 3D coordinates
        feature_factory_path: Path to custom .fdef file (None = RDKit default)

    Returns:
        List of PharmFeature objects with 3D coordinates.
    """
    from rdkit.Chem import RDConfig
    import os

    if feature_factory_path is None:
        feature_factory_path = os.path.join(RDConfig.RDDataDir, "BaseFeatures.fdef")

    factory = ChemicalFeatures.BuildFeatureFactory(feature_factory_path)
    feats = factory.GetFeaturesForMol(mol)
    conf = mol.GetConformer(conf_id)

    features = []
    for feat in feats:
        family = feat.GetFamily()
        pos = feat.GetPos()
        features.append(PharmFeature(
            name=family,
            coords=(round(pos.x, 3), round(pos.y, 3), round(pos.z, 3)),
        ))
    return features


def generate_consensus_pharmacophore(
    active_mols: list[Chem.Mol],
    min_feature_frequency: float = 0.7,
    distance_tolerance: float = 1.5,
) -> list[PharmFeature]:
    """
    Build a consensus pharmacophore from a set of active molecules.

    Algorithm:
    1. Extract pharmacophoric features from all active conformers
    2. Cluster features of the same type within distance_tolerance Å
    3. Retain only clusters present in >= min_feature_frequency fraction of actives

    Args:
        active_mols: List of 3D-embedded RDKit Mol objects (actives only)
        min_feature_frequency: Minimum fraction of actives a cluster must appear in
        distance_tolerance: Å radius for merging nearby feature centroids

    Returns:
        Consensus pharmacophore feature list.
    """
    from sklearn.cluster import DBSCAN

    all_features: dict[str, list[np.ndarray]] = {}
    source_molecules: dict[str, list[int]] = {}

    for mol_idx, mol in enumerate(active_mols):
        if mol.GetNumConformers() == 0:
            AllChem.EmbedMolecule(mol, AllChem.ETKDGv3())
        feats = extract_pharmacophore_features(mol)
        for f in feats:
            if f.name not in all_features:
                all_features[f.name] = []
                source_molecules[f.name] = []
            all_features[f.name].append(np.array(f.coords))
            source_molecules[f.name].append(mol_idx)

    n_actives = len(active_mols)
    consensus = []

    for feat_name, coords_list in all_features.items():
        coords = np.array(coords_list)
        mol_ids = source_molecules[feat_name]
        if len(coords) < 2:
            continue

        db = DBSCAN(eps=distance_tolerance, min_samples=1).fit(coords)
        for cluster_id in set(db.labels_):
            if cluster_id == -1:
                continue
            mask = db.labels_ == cluster_id
            unique_mols = len(set(m for m, keep in zip(mol_ids, mask) if keep))
            frequency = unique_mols / n_actives
            if frequency >= min_feature_frequency:
                centroid = coords[mask].mean(axis=0)
                consensus.append(PharmFeature(
                    name=feat_name,
                    coords=tuple(centroid.tolist()),
                    importance=frequency,
                ))

    return consensus

3D Pharmacophore Virtual Screening

import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem
from scipy.spatial.distance import cdist


def pharmacophore_screen(
    query_pharmacophore: list[PharmFeature],
    database_smiles: list[str],
    tolerance_A: float = 1.5,
    n_conformers: int = 10,
    required_features: list[str] = None,
    min_match_fraction: float = 0.75,
) -> list[tuple[str, float]]:
    """
    Screen a SMILES database against a 3D pharmacophore query.

    For each database molecule:
    1. Generate n_conformers with ETKDGv3
    2. Extract pharmacophoric features from each conformer
    3. Compute pharmacophore match score: fraction of required query features matched
       within tolerance_A in at least one conformer
    4. Return molecules meeting min_match_fraction threshold

    Args:
        query_pharmacophore: Consensus pharmacophore from generate_consensus_pharmacophore()
        database_smiles: List of SMILES to screen
        tolerance_A: Matching sphere radius per feature (Å)
        n_conformers: Number of conformers to generate per molecule
        required_features: Feature types that must be matched (None = all features)
        min_match_fraction: Minimum fraction of pharmacophore features to match

    Returns:
        List of (SMILES, match_score) for molecules that pass, sorted by score.
    """
    if required_features is None:
        required_features = list({f.name for f in query_pharmacophore})

    query_by_type: dict[str, list[np.ndarray]] = {}
    for f in query_pharmacophore:
        if f.name not in query_by_type:
            query_by_type[f.name] = []
        query_by_type[f.name].append(np.array(f.coords))

    hits = []
    for smi in database_smiles:
        mol = Chem.MolFromSmiles(smi)
        if mol is None:
            continue
        mol = Chem.AddHs(mol)
        params = AllChem.ETKDGv3()
        params.numThreads = 1
        AllChem.EmbedMultipleConfs(mol, numConfs=n_conformers, params=params)
        if mol.GetNumConformers() == 0:
            continue

        best_score = 0.0
        for conf_id in range(mol.GetNumConformers()):
            mol_feats = extract_pharmacophore_features(mol, conf_id=conf_id)
            mol_by_type: dict[str, list[np.ndarray]] = {}
            for f in mol_feats:
                if f.name not in mol_by_type:
                    mol_by_type[f.name] = []
                mol_by_type[f.name].append(np.array(f.coords))

            matched = 0
            total = len({f.name for f in query_pharmacophore if f.name in required_features})
            for feat_type in required_features:
                if feat_type not in query_by_type or feat_type not in mol_by_type:
                    continue
                q_coords = np.array(query_by_type[feat_type])
                m_coords = np.array(mol_by_type[feat_type])
                dist_matrix = cdist(q_coords, m_coords)
                if (dist_matrix < tolerance_A).any():
                    matched += 1

            score = matched / max(total, 1)
            best_score = max(best_score, score)

        if best_score >= min_match_fraction:
            hits.append((smi, round(best_score, 3)))

    return sorted(hits, key=lambda x: x[1], reverse=True)

Pharmacophore Feature Types

Key Takeaways

• Pharmacophore screening is complementary to docking: it is 100–1000× faster, making it effective for pre-filtering large databases before the expensive docking step. A typical workflow is pharmacophore filter → 2D similarity → docking → rescoring.
• Consensus pharmacophores from 5–10 chemically diverse actives capture the essential binding features while removing idiosyncratic group contributions. Using only 1–2 actives typically over-fits the pharmacophore to a specific scaffold and misses genuinely novel hits.
• Conformational sampling is the key bottleneck for pharmacophore screening: flexible molecules must be sampled in their bioactive conformation to match the query. Using 10–50 conformers per molecule (ETKDGv3 + clustering) reliably recovers known actives at acceptable false-negative rates.
• Exclusion volumes (derived from the protein active site shape) should always accompany pharmacophore models — without them, the model will accept molecules that match all features but clash with the protein and could never actually bind.

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.