Bio Pharmacophore Modeling

Version Compatibility

Reference examples tested with: RDKit 2024.09+, Pharmit web service, and PLIP 2.4+ (interaction analysis). Verify the deployed Pharmit/Pharmer interface and query format before automation.

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show rdkit then help(rdkit.Chem.Pharm3D) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Pharmacophore Modeling

Build 3D pharmacophore queries that capture the essential interaction features of a ligand-target binding event. A pharmacophore is the spatial arrangement of pharmacophore features (donor, acceptor, hydrophobe, aromatic, charged) sufficient for activity, abstracted from any specific chemotype. Use pharmacophores for scaffold hopping, virtual-screening prefilters, and cross-target SAR transfer. Derive interaction features directly from a co-crystal when available, use apo2ph4 to derive models from an apo pocket (Heider et al. 2023), or align known actives for a ligand-based model. PharmacoForge generates candidate 3D pharmacophores conditioned on a protein pocket; those pharmacophores can then retrieve matching molecules from a library (Flynn et al. 2025).

For 2D scaffold-based searches, see chemoinformatics/scaffold-analysis. For 3D shape similarity, see chemoinformatics/shape-similarity. For protein-ligand interaction analysis, see chemoinformatics/virtual-screening.

Pharmacophore Feature Types

Tolerances are pharmacophore-feature distance windows in the search. Tighter tolerances = fewer hits but more specific.

The ranges in this table are repository starting heuristics, not universal feature tolerances. Set final bounds from aligned-feature variability, coordinate uncertainty, and retrospective validation for the selected search engine.

The one-letter labels above are human-readable shorthand, not RDKit API codes. RDKit's shipped BaseFeatures.fdef uses family names such as Donor, Acceptor, Hydrophobe, Aromatic, PosIonizable, and NegIonizable. Its default feature definitions do not provide every halogen-bond or metal-coordination model; add and validate project-specific feature definitions when those interactions matter.

Method Taxonomy

Decision Tree by Scenario

Ligand-Based Pharmacophore (RDKit Pharm3D)

Goal: Derive a common pharmacophore from aligned bioactive conformers, then apply that established model to candidate molecules.

Approach: Consensus derivation is a separate modeling step: select or generate plausible bioactive conformers, align them using a documented method, identify conserved feature correspondences, and estimate distance bounds or tolerances. RDKit does not provide a single EmbedPharmacophore call that performs those steps. EmbedPharmacophore instead generates conformations of a molecule that satisfy an already defined pharmacophore.

from rdkit import Chem, Geometry
from rdkit.Chem import ChemicalFeatures
from rdkit.Chem.Pharm3D import EmbedLib, Pharmacophore
from rdkit.RDPaths import RDDataDir
import os

fdef_file = os.path.join(RDDataDir, 'BaseFeatures.fdef')
factory = ChemicalFeatures.BuildFeatureFactory(fdef_file)

# This is an already defined model. Coordinates and bounds must come from a
# validated consensus-derivation workflow or another justified source. RDKit
# requires FreeChemicalFeature objects, not feature-family strings.
query_features = [
    ChemicalFeatures.FreeChemicalFeature(
        'Aromatic', Geometry.Point3D(0.0, 0.0, 0.0)),
    ChemicalFeatures.FreeChemicalFeature(
        'Donor', Geometry.Point3D(4.0, 0.0, 0.0)),
]
pharmacophore = Pharmacophore.Pharmacophore(query_features)
pharmacophore.setLowerBound(0, 1, 3.5)
pharmacophore.setUpperBound(0, 1, 5.0)

target = Chem.AddHs(Chem.MolFromSmiles('c1ccc(cc1)CCN'))
can_match, feature_matches = EmbedLib.MatchPharmacophoreToMol(
    target, factory, pharmacophore)
if can_match:
    atom_match = tuple(tuple(matches[0].GetAtomIds())
                       for matches in feature_matches)
    _, embeddings, n_failed = EmbedLib.EmbedPharmacophore(
        target, atom_match, pharmacophore, randomSeed=23, silent=True)

BaseFeatures.fdef (RDKit-shipped) defines feature SMARTS and is a useful starting feature taxonomy. The code above demonstrates applying an existing two-feature model; it does not infer a consensus model from active compounds.

Receptor-Based Pharmacophore (apo2ph4 workflow)

Goal: Derive a pharmacophore from a protein binding-pocket structure without requiring a bound ligand.

Approach: Identify donor, acceptor, and hydrophobic hot spots from apo-pocket geometry, cluster them, and assemble candidate pharmacophores. Heider et al. describe apo2ph4 in J. Chem. Inf. Model. 63:101-110 (2023). Use the source release's documented scripts and environment rather than assuming a packaged apo2ph4 command: the published workflow writes LigandScout PML output, not a generic .ph4 file. Treat conversion to Pharmit, Pharmer, MOE, or Phase as a separate, explicitly validated step because pharmacophore formats are not interchangeable.

When a co-crystal ligand is available, derive pharmacophore directly from the ligand binding pose: each ligand feature in contact with a complementary protein residue is part of the pharmacophore.

from plip.basic import config
from plip.structure.preparation import PDBComplex

mol_complex = PDBComplex()
mol_complex.load_pdb('complex.pdb')
mol_complex.analyze()

for site in mol_complex.interaction_sets.values():
    for interaction in site.all_itypes:
        # Objects are interaction-class-specific. Inspect the documented fields
        # for HydrophobicContact, HydrogenBond, PiStacking, SaltBridge, etc.;
        # there is no universal `.type` or `.ligatom.coords` interface.
        interaction_class = type(interaction).__name__
        print(interaction_class, interaction)

PLIP exposes typed interaction records with class-specific ligand/protein atoms and coordinates. Map those records to pharmacophore features explicitly and retain the interaction class and source atom identifiers.

Pharmacophore Search (Pharmit / Pharmer)

For library screening, configure feature types, centers, radii, and optional shape constraints in Pharmit, or use a Pharmer database and query produced in the format required by the installed release. Do not pass LigandScout PML or a vendor .ph4 file directly unless the selected interface documents that import path. Pharmit reported searching millions of conformers in seconds to minutes; actual runtime depends on query selectivity, database size, and deployment (Sunseri & Koes 2016).

Pharmacophore Quality Validation

Evaluate a pharmacophore by:

1. Retrospective enrichment: a stated metric on target-relevant actives and inactives/decoys. DUD-E can provide a benchmark with known decoy-construction biases; COCONUT is a natural-products collection, not a target-specific active/decoy benchmark.
2. Geometric tightness: feature distance variance across actives
3. Selectivity: false positives in inactive set should be low
4. Specific consistency: pharmacophore matches each active's bioactive conformer

def pharmacophore_enrichment(query_pharmacophore, actives, inactives,
                             matches_pharmacophore):
    """Return active/inactive match-rate enrichment for a supplied matcher."""
    if not actives or not inactives:
        raise ValueError('actives and inactives must both be non-empty')
    n_active_match = sum(
        bool(matches_pharmacophore(mol, query_pharmacophore))
        for mol in actives)
    n_inactive_match = sum(
        bool(matches_pharmacophore(mol, query_pharmacophore))
        for mol in inactives)
    active_rate = n_active_match / len(actives)
    inactive_rate = n_inactive_match / len(inactives)
    return float('inf') if inactive_rate == 0 else active_rate / inactive_rate

For this repository, enrichment >=5x may be used as a starting triage heuristic only after the active/decoy construction and matching policy are documented. Report the full metric and uncertainty, and calibrate the acceptance threshold on the project dataset.

Pocket-Conditioned Pharmacophore Generation (PharmacoForge)

PharmacoForge (Flynn et al. 2025) applies a diffusion model to a protein pocket and generates candidate 3D pharmacophores. It does not directly generate molecular structures from an input pharmacophore. The validated workflow is:

1. Prepare the protein pocket in the representation required by the published PharmacoForge release.
2. Sample and rank pocket-conditioned pharmacophores.
3. Convert a selected pharmacophore into the query representation used by the search engine.
4. Retrieve matching, purchasable compounds and evaluate them with docking, strain, and physical-validity checks.

The paper compares pharmacophore and downstream retrieval performance with other pocket-based approaches; it does not support a drug-likeness or novelty comparison with REINVENT.