Bio Pharmacophore Modeling

Version Compatibility

Reference examples tested with: RDKit 2024.09+, pharmer / pharmit (web service), PharmIT 1.1+, plip 2.4+ (interaction analysis).

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. Used for scaffold-hopping (find compounds with different scaffold but matching pharmacophore), virtual screening (faster than docking), and cross-target SAR transfer. Modern best practice: derive pharmacophore from co-crystal structure if available (receptor-based; apo2ph4 workflow of Heider et al 2022/2023 J Chem Inf Model 63:147-158) or align actives if no crystal (ligand-based). Diffusion-based generation (PharmacoForge, Flynn et al 2025 Front Bioinform) lets pharmacophore drive de novo design.

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.

Method Taxonomy

Decision Tree by Scenario

Ligand-Based Pharmacophore (RDKit Pharm3D)

Goal: Extract a common pharmacophore from a set of bioactive compounds.

Approach: Align actives by maximum-common-substructure or shape; extract conserved features; derive a pharmacophore signature.

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

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

active_smiles = ['CC(C)c1ccc(C(=O)NCc2ccccn2)cc1', 'CCC(C)c1ccc(C(=O)NCc2ccccn2)cc1']
active_mols = [Chem.AddHs(Chem.MolFromSmiles(s)) for s in active_smiles]
for m in active_mols:
    AllChem.EmbedMolecule(m, AllChem.ETKDGv3())
    AllChem.MMFFOptimizeMolecule(m)

# Extract pharmacophore features (family/type/atom-ids/3D-pos) per molecule
feature_lists = [factory.GetFeaturesForMol(m) for m in active_mols]

# Build a target pharmacophore from one active's features; the bounds matrix
# encodes per-feature-pair distance ranges across all conformer/active variability.
# Below uses 2 features (Aromatic + Donor) with 2.5-3.5 A and 5.0-6.0 A bands.
feature_types = ['Aromatic', 'Donor']
# bounds is symmetric 2x2 distance-range matrix; off-diagonal = (lo, hi)
import numpy as np
bounds_matrix = np.array([[0.0, 5.0], [3.5, 0.0]])  # upper / lower bounds
pharmacophore = Pharmacophore.Pharmacophore(feature_types)
pharmacophore.setLowerBound(0, 1, 3.5)
pharmacophore.setUpperBound(0, 1, 5.0)

# To search a target molecule for matches, use Pharm3D.EmbedLib.EmbedPharmacophore
# (see chemoinformatics/conformer-generation for 3D embedding fundamentals)

BaseFeatures.fdef (RDKit-shipped) defines feature SMARTS. For drug-like pharmacophores, this is the standard starting point.

Receptor-Based Pharmacophore (apo2ph4 workflow)

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

Approach: Identify hot-spots (donor / acceptor / hydrophobe regions) from protein geometry; assemble into a pharmacophore. The apo2ph4 workflow of Heider J, Kilian J, Garifulina A, Hering S, Langer T, Seidel T (2022/2023 J Chem Inf Model 63:147-158) describes the conceptual pipeline; the example below is illustrative — verify the exact CLI invocation against the published code/release before running.

# Conceptual apo2ph4-style workflow; flags shown are illustrative.
apo2ph4 -pdb receptor.pdb -site_residues 'A:100,A:101,A:104,A:108' \
        -output pharmacophore.ph4

apo2ph4 outputs a .ph4-style pharmacophore compatible with Phase, MOE, and Pharmer.

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:
        # H-bond donor / acceptor / pi-stacking / hydrophobic / salt-bridge
        feature_type = interaction.type
        feature_atom_coords = interaction.ligatom.coords

PLIP outputs interaction types per ligand atom; each interaction maps to a pharmacophore feature.

Pharmacophore Search (Pharmit / Pharmer)

For library screening, Pharmer/Pharmit are the standard tools.

# Pharmer command line (offline alternative)
pharmer search -q pharmacophore.ph4 -dbdir zinc_db -out hits.sdf

Or use Pharmit web service (https://pharmit.csb.pitt.edu) for browser-based or REST search.

import requests

url = 'https://pharmit.csb.pitt.edu/...'  # specific endpoint
response = requests.post(url, json={'pharmacophore': pharmacophore_dict})
hits = response.json()

Pharmacophore screen is orders of magnitude faster than docking; typical 100M-compound search in minutes.

Pharmacophore Quality Validation

Evaluate a pharmacophore by:

1. Retrospective enrichment: AUC-like metric on actives vs decoys (DUD-E, COCONUT)
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):
    n_active_match = sum(matches_pharmacophore(a, query_pharmacophore) for a in actives)
    n_inactive_match = sum(matches_pharmacophore(d, query_pharmacophore) for d in inactives)
    enrichment = (n_active_match / len(actives)) / (n_inactive_match / len(inactives))
    return enrichment

A good pharmacophore: enrichment >= 5x relative to random.

Diffusion-Based Pharmacophore Design (PharmacoForge)

PharmacoForge (Flynn, Shah, Dunn et al 2025 Front Bioinform) generates molecules conditioned on a pharmacophore query using a diffusion model:

# Pseudo-code; depends on PharmacoForge installation
# from pharmacoforge import PharmacophoreDiffusionGenerator
# gen = PharmacophoreDiffusionGenerator()
# generated = gen.generate(pharmacophore_query, n_molecules=100)

Trade-off: PharmacoForge produces de novo molecules that satisfy pharmacophore; lower drug-likeness than REINVENT but better novelty.

Pharmacophore vs Shape vs 2D Fingerprint

Pharmacophore is more interpretable than shape: a hit explains why it matched (donor at position X, hydrophobe at position Y).

Per-Tool Failure Modes

Ligand-based -- diverse actives confound

Trigger: Active set spans multiple scaffolds with different bound conformations.

Mechanism: No common pharmacophore exists; algorithm forces non-consensus features.

Symptom: Pharmacophore matches no actives in retrospective.

Fix: Cluster actives by scaffold first; derive per-cluster pharmacophore.

Receptor-based -- apo structure

Trigger: Protein in apo form (no bound ligand).

Mechanism: Side-chain rotamers differ between apo and holo; "binding site" geometry is wrong.

Symptom: Pharmacophore inferred from apo doesn't match holo experimental data.

Fix: Use AlphaFold3 / Boltz-1 to predict holo conformation; derive pharmacophore from predicted holo.

Pharmacophore -- single conformer bias

Trigger: Active aligned to its first generated conformer, not bioactive conformer.

Mechanism: Crystal structure not available; generated conformer may not be the bound one.

Symptom: Pharmacophore inconsistent across runs (different starting conformer chosen).

Fix: Use conformer ensemble; align all to common scaffold; choose conformer most consistent with other actives.

Tolerance too tight

Trigger: Default geometric tolerance < 0.5 Å.

Mechanism: Real bioactive conformers have flexibility; rigid pharmacophore filters most molecules out.

Symptom: Search returns zero hits.

Fix: Use tolerance 1.0-1.5 Å for drug-like; up to 2 Å for flexible peptide-like.

Pharmacophore search misses bioisostere

Trigger: Bioisostere replacement (e.g., -COOH replaced by tetrazole).

Mechanism: Tetrazole functions as acid bioisostere but RDKit features may not classify identically.

Symptom: Known bioisosteric active not found.

Fix: Use ChemAxon-style bioisosteric feature equivalence; or pharmacophore feature class expansion (acid generic vs -COOH specific).

PLIP -- water-mediated interaction missed

Trigger: Bridging water between ligand donor and protein acceptor.

Mechanism: PLIP doesn't include explicit water.

Symptom: Pharmacophore missing critical H-bond feature.

Fix: Add water-mediated interaction manually based on crystal water positions; or use PoseView for full interaction view.

Reconciliation: Ligand-Based vs Receptor-Based

For prospective scaffold-hopping, receptor-based is preferred. For ranking analogs, ligand-based is sufficient.

Common Errors

References

• Wolber & Langer, J. Chem. Inf. Model. 45:160 (2005) -- LigandScout pharmacophore.
• Heider J, Kilian J, Garifulina A, Hering S, Langer T, Seidel T 2022/2023 J Chem Inf Model 63(1):147-158 -- apo2ph4 receptor-based pharmacophore workflow (DOI 10.1021/acs.jcim.2c00814).
• Flynn JR, Shah RH, Dunn SD et al 2025 Front Bioinform -- PharmacoForge diffusion pharmacophore.
• RDKit Chem.Pharm3D framework -- documentation at rdkit.org (the unconfirmed "Stiefl 2021" citation has been removed pending verification).
• Adasme et al., Nucleic Acids Res. 49:W530 (2021) -- PLIP interaction profiler.
• Koes et al., Nucleic Acids Res. 44:W436 (2016) -- Pharmit interactive search.

Related Skills

• chemoinformatics/molecular-io - Parse molecules
• chemoinformatics/conformer-generation - Generate 3D for pharmacophore
• chemoinformatics/shape-similarity - 3D shape adjacent to pharmacophore
• chemoinformatics/virtual-screening - Pharmacophore as docking pre-filter
• chemoinformatics/scaffold-analysis - 2D scaffold-hopping context
• chemoinformatics/generative-design - PharmacoForge for de novo
• structural-biology/structure-io - PDB handling