48.1 Drug Discovery Pipeline

1. The Drug Discovery Pipeline

Every year, pharmaceutical companies abandon late-stage drug candidates after investing hundreds of millions of dollars, often because a flaw that was detectable at the screening stage slipped through unnoticed. Understanding exactly where candidates fail, and why, is the prerequisite for building AI systems that catch those flaws early.

Of the roughly 10⁶⁰ drug-like molecules that are theoretically possible, fewer than 10⁴ have ever reached a pharmacy shelf, and each of those survivors navigated a gauntlet that discards 99.9% of candidates at every stage. Understanding this pipeline is essential before applying AI to any part of it, because the value of a computational prediction depends entirely on where in the pipeline it sits.

The pipeline has five major stages:

1. Target identification: finding the biological molecule (usually a protein) whose modulation will treat the disease.
2. Hit finding: discovering chemical matter that interacts with the target, typically through high-throughput screening (HTS) or virtual screening.
3. Lead optimization: iteratively modifying hits to improve potency, selectivity, and drug-like properties.
4. Preclinical development: evaluating toxicity, pharmacokinetics, and formulation in animal models.
5. Clinical trials: testing safety (Phase I), efficacy (Phase II), and comparative effectiveness (Phase III) in humans.

Figure 48.1 maps these five stages as a narrowing funnel, highlighting where AI methods currently deliver the greatest acceleration. Figure 48.1.1 illustrates Drug discovery pipeline funnel with AI intervention points.

AI has its greatest current impact in stages 1 through 3. Preclinical and clinical stages involve regulatory constraints, biological variability, and ethical considerations that limit purely computational approaches (though clinical AI, covered in Section 48.4, addresses trial design and patient stratification). In short: AI's highest leverage is at the cheapest, most uncertain end of the funnel, where a better computational filter today prevents a billion-dollar clinical failure years from now.

2. Target Identification with Genomics and Knowledge Graphs

Target identification asks: which protein, when inhibited or activated, will alter the disease phenotype? Traditional approaches relied on phenotypic screens (assays that measure a compound's effect on cells or organisms without knowing the molecular target) and literature curation. Modern AI-driven target ID combines three data sources: genome-wide association studies (GWAS), transcriptomic differential expression, and biomedical knowledge graphs.

Knowledge graphs for target identification connect genes, diseases, pathways, drugs, and side effects in a heterogeneous graph. Link prediction on this graph identifies gene-disease associations not yet reported in the literature. The approach from Chapter 38 applies directly. Train a graph embedding model (TransE, RotatE, or a GNN (graph neural network)-based encoder) on known associations, then rank candidate gene-disease links by predicted score.

import numpy as np
from pykeen.pipeline import pipeline
from pykeen.triples import TriplesFactory

def train_target_identification_model(
    triples_path: str,
    embedding_dim: int = 256,
    num_epochs: int = 200,
):
    """Train a knowledge graph embedding model for target identification.

    The triples file contains (gene, relation, disease) tuples from
    curated databases like DisGeNET, OMIM, and DrugBank.
    """
    tf = TriplesFactory.from_path(triples_path)
    training, testing, validation = tf.split([0.8, 0.1, 0.1])

    result = pipeline(
        training=training,
        testing=testing,
        validation=validation,
        model="RotatE",
        model_kwargs={"embedding_dim": embedding_dim},
        optimizer="Adam",
        optimizer_kwargs={"lr": 1e-3},
        training_kwargs={
            "num_epochs": num_epochs,
            "batch_size": 1024,
        },
        evaluation_kwargs={"batch_size": 256},
    )

    return result.model, result.metric_results


def predict_novel_targets(model, disease_entity: str, top_k: int = 50):
    """Predict novel gene targets for a disease.

    Returns genes ranked by predicted association score that are
    NOT in the training set (truly novel predictions).
    """
    predictions = model.predict_target(
        head=disease_entity,
        relation="associated_with",
    )
    return predictions.topk(top_k)

Mendelian randomization provides causal evidence for target-disease associations by using genetic variants as instrumental variables (factors that influence the outcome only through the gene of interest, creating a natural experiment). If a variant that reduces expression of gene \(G\) also reduces risk of disease \(D\), this constitutes causal evidence that inhibiting \(G\)'s protein product could treat \(D\). The causal inference methods from Chapter 31 formalize this reasoning.

3. Hit Finding: Virtual Screening

Once a target protein is identified, the next step is finding molecules that bind to it. Physical high-throughput screening tests hundreds of thousands of compounds in robotic assays. Virtual screening replaces physical testing with computational prediction, enabling screening of billions of molecules. (In reported prospective studies, AI-guided virtual screening has typically achieved hit rates of 10 to 12%, compared to roughly 0.1% from random physical screening: a 100x enrichment that turns a needle-in-a-haystack search into a targeted retrieval.)

Virtual screening evaluates large chemical libraries computationally to identify molecules likely to bind a biological target. No compound is synthesized or physically tested at this stage. Physical screening of even a million compounds costs millions of dollars and weeks of robotic assay time. Virtual screening evaluates billions of candidates on a GPU cluster in hours for a fraction of that cost. The process has two stages. First, each candidate molecule is represented numerically (as a fingerprint, a graph, or a 3D coordinate set). Then a scoring function ranks molecules by predicted binding strength or structural similarity to known actives. Use virtual screening when a validated target structure or known active ligands exist; when neither is available, phenotypic screening or fragment-based approaches remain necessary starting points.

Two complementary approaches dominate virtual screening:

3.1 Ligand-Based Virtual Screening

When known active molecules exist, ligand-based methods search for similar compounds. The key representation is the molecular fingerprint: a fixed-length binary or count vector encoding the presence of chemical substructures. Each molecule is specified as a SMILES (Simplified Molecular-Input Line-Entry System) string, a compact text notation where atoms are letters and bonds are symbols. RDKit provides several fingerprint types; the Morgan (circular) fingerprint with radius 2 is the most widely used:

from rdkit import Chem
from rdkit.Chem import AllChem, DataStructs
import numpy as np

def compute_morgan_fingerprint(smiles: str, radius: int = 2, n_bits: int = 2048):
    """Compute Morgan circular fingerprint from SMILES string.

    The Morgan fingerprint encodes circular substructures up to
    the specified radius. Each atom's neighborhood is hashed to
    a bit position, creating a fixed-length representation of
    the molecule's local chemical environments.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        raise ValueError(f"Invalid SMILES: {smiles}")
    fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius, nBits=n_bits)
    return np.array(fp)


def tanimoto_similarity(fp1: np.ndarray, fp2: np.ndarray) -> float:
    """Tanimoto coefficient between two binary fingerprints.

    T(A, B) = |A ∩ B| / |A ∪ B|

    Values above 0.4 indicate meaningful structural similarity.
    Values above 0.7 suggest similar biological activity.
    """
    intersection = np.sum(fp1 & fp2)
    union = np.sum(fp1 | fp2)
    if union == 0:
        return 0.0
    return intersection / union


def virtual_screen_ligand_based(
    query_smiles: str,
    library_smiles: list[str],
    threshold: float = 0.4,
) -> list[tuple[str, float]]:
    """Screen a compound library by similarity to a query molecule.

    Returns compounds with Tanimoto similarity above threshold,
    sorted by decreasing similarity.
    """
    query_fp = compute_morgan_fingerprint(query_smiles)
    hits = []
    for smiles in library_smiles:
        try:
            lib_fp = compute_morgan_fingerprint(smiles)
            sim = tanimoto_similarity(query_fp, lib_fp)
            if sim >= threshold:
                hits.append((smiles, sim))
        except ValueError:
            continue
    return sorted(hits, key=lambda x: x[1], reverse=True)

The Tanimoto coefficient between two fingerprints \(A\) and \(B\) is defined as:

For binary fingerprints, this simplifies to the ratio of shared "on" bits to the total number of "on" bits in either fingerprint.

3.2 Structure-Based Virtual Screening (Docking)

When the target protein's 3D structure is available (from X-ray crystallography, cryogenic electron microscopy (cryo-EM), or AlphaFold prediction), molecular docking predicts how a small molecule binds in the protein's active site. Docking scores approximate binding free energy:

where the terms capture van der Waals interactions, electrostatics, hydrogen bonds, desolvation penalties (the energy cost of stripping water molecules from both the ligand and the binding pocket before they can interact), and torsional strain. More negative scores indicate stronger predicted binding.

import subprocess
import os
from dataclasses import dataclass

@dataclass
class DockingResult:
    ligand_smiles: str
    score: float  # kcal/mol (more negative = better)
    pose_file: str

def dock_molecule(
    ligand_smiles: str,
    receptor_pdbqt: str,
    center: tuple[float, float, float],
    box_size: tuple[float, float, float] = (25.0, 25.0, 25.0),
    exhaustiveness: int = 32,
    work_dir: str = "/tmp/docking",
) -> DockingResult:
    """Dock a single molecule against a receptor using AutoDock Vina.

    The binding site is defined by a center coordinate and box size
    in Angstroms. Exhaustiveness controls search thoroughness.
    """
    os.makedirs(work_dir, exist_ok=True)

    # Convert SMILES to 3D structure and prepare for docking
    from rdkit import Chem
    from rdkit.Chem import AllChem
    mol = Chem.MolFromSmiles(ligand_smiles)
    mol = Chem.AddHs(mol)
    AllChem.EmbedMolecule(mol, randomSeed=42)
    AllChem.MMFFOptimizeMolecule(mol)

    ligand_sdf = os.path.join(work_dir, "ligand.sdf")
    Chem.MolToMolFile(mol, ligand_sdf)

    # Convert to PDBQT format (required by Vina)
    ligand_pdbqt = os.path.join(work_dir, "ligand.pdbqt")
    subprocess.run([
        "mk_prepare_ligand.py", "-i", ligand_sdf, "-o", ligand_pdbqt
    ], check=True)

    # Run Vina docking
    output_pdbqt = os.path.join(work_dir, "output.pdbqt")
    result = subprocess.run([
        "vina",
        "--receptor", receptor_pdbqt,
        "--ligand", ligand_pdbqt,
        "--center_x", str(center[0]),
        "--center_y", str(center[1]),
        "--center_z", str(center[2]),
        "--size_x", str(box_size[0]),
        "--size_y", str(box_size[1]),
        "--size_z", str(box_size[2]),
        "--exhaustiveness", str(exhaustiveness),
        "--out", output_pdbqt,
    ], capture_output=True, text=True, check=True)

    # Parse best score from Vina output
    best_score = _parse_vina_score(result.stdout)

    return DockingResult(
        ligand_smiles=ligand_smiles,
        score=best_score,
        pose_file=output_pdbqt,
    )


def _parse_vina_score(vina_output: str) -> float:
    """Extract the best binding affinity from Vina output."""
    for line in vina_output.split("\n"):
        if line.strip().startswith("1"):
            parts = line.split()
            return float(parts[1])
    raise ValueError("Could not parse Vina output")

4. Lead Optimization

Screening billions of molecules to find thousands of hits solves only the first half of the problem; the harder challenge is transforming a weak, imperfect binder into a compound safe and potent enough to enter a patient.

Lead optimization iterates on a hit's chemical structure to improve potency (binding affinity), selectivity (avoiding off-target binding), and drug-likeness (oral bioavailability, metabolic stability, low toxicity). This multi-objective optimization problem maps directly to the methods from Chapter 45.

4.1 ADMET Prediction

ADMET stands for Absorption, Distribution, Metabolism, Excretion, and Toxicity. These pharmacokinetic properties determine whether a molecule that binds a target in a test tube will work as a drug in a human body. AI models predict ADMET properties from molecular structure, enabling computational triage before synthesis.

from rdkit import Chem
from rdkit.Chem import Descriptors, Lipinski

@dataclass
class ADMETProfile:
    molecular_weight: float
    logp: float
    hbd: int   # hydrogen bond donors
    hba: int   # hydrogen bond acceptors
    tpsa: float  # topological polar surface area
    rotatable_bonds: int
    lipinski_violations: int
    drug_likeness_score: float

def compute_admet_profile(smiles: str) -> ADMETProfile:
    """Compute ADMET-relevant molecular descriptors.

    Combines Lipinski's Rule of Five with additional descriptors
    (TPSA, rotatable bonds) that predict oral absorption.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        raise ValueError(f"Invalid SMILES: {smiles}")

    mw = Descriptors.MolWt(mol)
    logp = Descriptors.MolLogP(mol)
    hbd = Lipinski.NumHDonors(mol)
    hba = Lipinski.NumHAcceptors(mol)
    tpsa = Descriptors.TPSA(mol)
    rot_bonds = Lipinski.NumRotatableBonds(mol)

    # Count Lipinski violations
    violations = sum([
        mw > 500,
        logp > 5,
        hbd > 5,
        hba > 10,
    ])

    # Quantitative Estimate of Drug-likeness (QED)
    from rdkit.Chem.QED import qed
    qed_score = qed(mol)

    return ADMETProfile(
        molecular_weight=mw,
        logp=logp,
        hbd=hbd,
        hba=hba,
        tpsa=tpsa,
        rotatable_bonds=rot_bonds,
        lipinski_violations=violations,
        drug_likeness_score=qed_score,
    )

4.2 Multi-Objective Lead Optimization

Lead optimization requires balancing competing objectives. A molecule with nanomolar binding affinity is useless if it is metabolically unstable or toxic. We formalize this as a multi-objective optimization problem:

subject to constraints on molecular weight, log P, and synthesizability score (a computed estimate of how difficult a molecule would be to manufacture in a chemistry lab, where lower scores indicate easier synthesis). The Pareto optimization framework from Section 45.2 identifies the set of non-dominated solutions: molecules where no single property can be improved without degrading another.

from rdkit.Chem import RDConfig
import os
import sys

# Multi-property optimization with molecular generation
def optimize_lead(
    seed_smiles: str,
    target_protein_pdbqt: str,
    binding_site_center: tuple[float, float, float],
    n_iterations: int = 50,
    population_size: int = 100,
) -> list[dict]:
    """Multi-objective lead optimization using genetic algorithm.

    Objectives:
    1. Maximize predicted binding affinity (docking score)
    2. Maximize drug-likeness (QED)
    3. Minimize synthetic accessibility (SA score)

    Uses matched molecular pair transformations to generate
    neighboring molecules from the seed compound.
    """
    from rdkit.Chem import AllChem, RDConfig
    from rdkit.Chem.SA_Score import sascorer

    seed_mol = Chem.MolFromSmiles(seed_smiles)
    population = [seed_smiles]

    # Generate initial population via random modifications
    for _ in range(population_size - 1):
        mutant = _random_molecular_mutation(seed_mol)
        if mutant is not None:
            population.append(Chem.MolToSmiles(mutant))

    pareto_front = []
    for iteration in range(n_iterations):
        # Evaluate all objectives
        scored = []
        for smi in population:
            try:
                profile = compute_admet_profile(smi)
                dock_result = dock_molecule(
                    smi, target_protein_pdbqt, binding_site_center
                )
                mol = Chem.MolFromSmiles(smi)
                sa_score = sascorer.calculateScore(mol)

                scored.append({
                    "smiles": smi,
                    "docking_score": dock_result.score,
                    "qed": profile.drug_likeness_score,
                    "sa_score": sa_score,
                })
            except Exception:
                continue

        # Non-dominated sorting (NSGA-II style, where NSGA-II is a
        # multi-objective genetic algorithm that ranks solutions by
        # Pareto dominance and crowding distance)
        pareto_front = _non_dominated_sort(scored)

        # Generate next population from Pareto-optimal parents
        population = _breed_next_generation(pareto_front, population_size)

    return pareto_front

Library Shortcut: TDC (Therapeutics Data Commons)

The Therapeutics Data Commons provides standardized benchmarks and datasets for every stage of the drug discovery pipeline. Instead of curating your own ADMET datasets and splitting them for evaluation, TDC provides ready-to-use data loaders with standard train/validation/test splits:

from tdc.single_pred import ADME

# Load the Caco-2 permeability dataset (intestinal absorption)
data = ADME(name="Caco2_Wang")
split = data.get_split(method="scaffold")  # Scaffold-based split

train = split["train"]  # SMILES + permeability values
valid = split["valid"]
test = split["test"]

print(f"Train: {len(train)}, Valid: {len(valid)}, Test: {len(test)}")
# Train: 728, Valid: 91, Test: 91

Listing 48.4b: Loading a standardized ADMET benchmark from the Therapeutics Data Commons. Scaffold-based splitting ensures that training and test molecules have distinct chemical scaffolds, preventing inflated accuracy from structural memorization.

TDC covers 66 datasets across 22 tasks (circa 2022; the collection has since grown) spanning target discovery, activity modeling, ADMET prediction, and drug-drug interaction. Using standard benchmarks rather than custom datasets takes approximately 5 lines instead of 200 lines of data curation, and ensures your results are comparable to published baselines.

5. Molecular Representations for Deep Learning

The choice of molecular representation determines what patterns a model can learn. Three representations dominate:

• Simplified Molecular-Input Line-Entry System (SMILES) strings: text representations processed by sequence models (LSTMs, transformers). Simple but fragile; small edits to a SMILES string can produce invalid molecules.
• Molecular graphs: atoms as nodes, bonds as edges, processed by graph neural networks (GNNs). Naturally capture molecular topology.
• 3D conformations: atom coordinates in 3D space, processed by equivariant neural networks. Required for modeling steric effects and binding poses.

import torch
from torch_geometric.data import Data

def smiles_to_graph(smiles: str) -> Data:
    """Convert a SMILES string to a PyTorch Geometric graph.

    Node features: atomic number, degree, formal charge,
    number of hydrogens, aromaticity.
    Edge features: bond type (single/double/triple/aromatic),
    conjugation, ring membership.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        raise ValueError(f"Invalid SMILES: {smiles}")

    # Node features: one per atom
    atom_features = []
    for atom in mol.GetAtoms():
        features = [
            atom.GetAtomicNum(),
            atom.GetDegree(),
            atom.GetFormalCharge(),
            atom.GetTotalNumHs(),
            int(atom.GetIsAromatic()),
            int(atom.IsInRing()),
        ]
        atom_features.append(features)

    x = torch.tensor(atom_features, dtype=torch.float)

    # Edge index and features: one per bond (bidirectional)
    edge_indices = []
    edge_features = []
    bond_type_map = {
        Chem.rdchem.BondType.SINGLE: [1, 0, 0, 0],
        Chem.rdchem.BondType.DOUBLE: [0, 1, 0, 0],
        Chem.rdchem.BondType.TRIPLE: [0, 0, 1, 0],
        Chem.rdchem.BondType.AROMATIC: [0, 0, 0, 1],
    }

    for bond in mol.GetBonds():
        i, j = bond.GetBeginAtomIdx(), bond.GetEndAtomIdx()
        bt = bond_type_map.get(bond.GetBondType(), [0, 0, 0, 0])
        feat = bt + [int(bond.GetIsConjugated()), int(bond.IsInRing())]
        # Add both directions
        edge_indices.extend([[i, j], [j, i]])
        edge_features.extend([feat, feat])

    edge_index = torch.tensor(edge_indices, dtype=torch.long).t().contiguous()
    edge_attr = torch.tensor(edge_features, dtype=torch.float)

    return Data(x=x, edge_index=edge_index, edge_attr=edge_attr)

6. Where AI Delivers and Where It Falls Short

With molecular representations and scoring functions in hand, the natural question becomes: how much of the drug discovery pipeline can these tools actually accelerate today, and where do fundamental gaps remain?

Assessing AI in drug discovery requires separating demonstrated impact from aspirational claims:

Current Scorecard

Demonstrated impact:

• Virtual screening enrichment: AI-guided screening typically achieves 10 to 100x enrichment over random selection in reported prospective studies.
• ADMET prediction: models trained on public datasets (ChEMBL, a curated database of bioactive molecules with drug-like properties maintained by the European Bioinformatics Institute, and TDC) predict solubility, permeability, and cytochrome P450 (CYP) inhibition with accuracy that, in benchmark evaluations, is generally sufficient to eliminate clearly unsuitable candidates before synthesis.
• Structure prediction: AlphaFold2/3 provides experimental-quality structures for targets without crystallographic data, enabling structure-based drug design on previously intractable targets.
• Synthesis planning: retrosynthesis models (covered in Chapter 49) identify viable synthetic routes, reducing the time from hit to synthesized compound.

Persistent challenges:

• Binding affinity prediction: absolute binding free energy calculations remain unreliable. Docking scores rank compounds roughly but cannot predict nanomolar affinity differences between closely related molecules.
• Off-target effects: predicting which of the body's thousands of proteins a molecule will inadvertently bind remains an open problem.
• Clinical translation: no AI-discovered drug has yet completed Phase III trials from a fully AI-driven discovery pipeline (as of mid-2026), though several are in Phase II.
• Data quality: public bioactivity databases contain measurement errors, assay variability, and reporting biases that limit model accuracy.

7. Discovery Workbench Integration

The Discovery Workbench gains a DrugDiscoveryPipeline component that orchestrates target identification, virtual screening, ADMET filtering, and lead optimization as a single workflow. Each stage logs its inputs, parameters, and outputs to the experiment registry (Chapter 47), enabling full provenance tracking from initial target hypothesis to optimized lead compound.

from discovery_workbench import Pipeline, ExperimentRegistry

class DrugDiscoveryPipeline(Pipeline):
    """End-to-end drug discovery pipeline for the Discovery Workbench."""

    def __init__(self, registry: ExperimentRegistry):
        self.registry = registry
        self.stages = [
            "target_identification",
            "virtual_screening",
            "admet_filtering",
            "lead_optimization",
        ]

    def run(
        self,
        disease: str,
        compound_library: list[str],
        target_structure: str,
        binding_site: tuple[float, float, float],
    ) -> dict:
        """Execute the full pipeline with provenance tracking."""
        run_id = self.registry.start_run(
            pipeline="drug_discovery",
            params={"disease": disease, "library_size": len(compound_library)},
        )

        # Stage 1: Target identification (uses KG model)
        targets = predict_novel_targets(self.kg_model, disease, top_k=10)
        self.registry.log_stage(run_id, "target_id", targets)

        # Stage 2: Virtual screening
        hits = virtual_screen_ligand_based(
            query_smiles=targets[0].known_ligand,
            library_smiles=compound_library,
            threshold=0.35,
        )
        self.registry.log_stage(run_id, "screening", {"n_hits": len(hits)})

        # Stage 3: ADMET filtering
        filtered = []
        for smiles, sim in hits:
            profile = compute_admet_profile(smiles)
            if profile.lipinski_violations <= 1 and profile.drug_likeness_score > 0.3:
                filtered.append((smiles, sim, profile))
        self.registry.log_stage(run_id, "admet", {"n_passed": len(filtered)})

        # Stage 4: Docking and ranking
        docked = []
        for smiles, sim, profile in filtered[:500]:  # Top 500 by similarity
            result = dock_molecule(smiles, target_structure, binding_site)
            docked.append({
                "smiles": smiles,
                "similarity": sim,
                "docking_score": result.score,
                "qed": profile.drug_likeness_score,
            })

        ranked = sorted(docked, key=lambda x: x["docking_score"])
        self.registry.log_stage(run_id, "ranked_leads", ranked[:50])
        self.registry.end_run(run_id)

        return {"run_id": run_id, "top_leads": ranked[:50]}