PROTAC Design

PROTACs (Proteolysis-Targeting Chimeras) are bifunctional molecules that simultaneously bind a target protein of interest (POI) and an E3 ubiquitin ligase, recruiting the ligase to the POI and triggering its ubiquitination and subsequent proteasomal degradation. Unlike classical inhibitors that occupy an active site, PROTACs degrade the entire protein — including scaffolding and allosteric functions — and are catalytic (the PROTAC is released after ubiquitin transfer and can degrade another molecule of POI).

The key computational challenge in PROTAC design is modelling the ternary complex (E3 ligase : PROTAC : POI) and optimising the linker that connects the two warheads so that the protein surfaces are oriented productively for ubiquitin transfer.

PROTAC Architecture

PROTAC Property Filters

PROTACs violate the Lipinski Rule of Five (MW typically 700–1100 Da). They pass through cells via active transport or endocytosis, and special property rules apply:

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

def protac_property_filter(smiles: str) -> dict:
    """
    Apply PROTAC-specific property filters (Bemis et al. / Ciulli lab guidelines).

    Recommended PROTAC property space:
    - MW: 700–1100 Da
    - cLogP: 2–6
    - HBD: ≤ 5
    - HBA: ≤ 15
    - RotBonds: 15–30 (flexibility needed for ternary complex)
    - PSA: 150–300 Ų (accepts larger values vs. drugs)
    - TPSA: ≤ 300 Ų for cell permeability

    Args:
        smiles: PROTAC SMILES string

    Returns:
        Property dict with pass/fail annotations.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return {"error": "Invalid SMILES"}

    mw = Descriptors.ExactMolWt(mol)
    clogp = Descriptors.MolLogP(mol)
    hbd = Lipinski.NumHDonors(mol)
    hba = Lipinski.NumHAcceptors(mol)
    rot = rdMolDescriptors.CalcNumRotatableBonds(mol)
    tpsa = rdMolDescriptors.CalcTPSA(mol)
    n_heavy = mol.GetNumHeavyAtoms()

    violations = []
    if not (700 <= mw <= 1100):
        violations.append(f"MW={mw:.0f} outside 700–1100 Da")
    if not (2 <= clogp <= 6):
        violations.append(f"cLogP={clogp:.2f} outside 2–6")
    if hbd > 5:
        violations.append(f"HBD={hbd} > 5")
    if hba > 15:
        violations.append(f"HBA={hba} > 15")
    if not (10 <= rot <= 35):
        violations.append(f"RotBonds={rot} outside 10–35")
    if tpsa > 300:
        violations.append(f"TPSA={tpsa:.0f} > 300 Ų")

    return {
        "smiles": smiles,
        "mw": round(mw, 1),
        "clogp": round(clogp, 2),
        "hbd": hbd,
        "hba": hba,
        "rotatable_bonds": rot,
        "tpsa": round(tpsa, 1),
        "n_heavy_atoms": n_heavy,
        "violations": violations,
        "passes_protac_rules": len(violations) == 0,
    }

E3 Ligase Warhead Library

from rdkit import Chem
from rdkit.Chem import AllChem
import pandas as pd

# Key E3 ligase warheads for PROTAC design
E3_WARHEADS = {
    # CRBN (Cereblon) binders — thalidomide analogs
    "thalidomide": {
        "smiles": "O=C1CCC(=O)N1C1CCC(=O)NC1=O",
        "e3_ligase": "CRBN",
        "pchembl_range": (6.0, 8.0),
        "notes": "Classic CRBN binder; base scaffold for IMiDs",
    },
    "pomalidomide": {
        "smiles": "Nc1cccc2C(=O)N(C3CCC(=O)NC3=O)C(=O)c12",
        "e3_ligase": "CRBN",
        "pchembl_range": (6.5, 8.5),
        "notes": "More potent CRBN binder; amine for linker attachment",
    },
    # VHL (von Hippel-Lindau) binders
    "vhl_ligand_vh032": {
        "smiles": "CC(C)(C)OC(=O)NC1CC(O)(C(=O)NC2CC(c3ccc(Cl)cc3)N(C(=O)c3ccccc3)C2)C1",
        "e3_ligase": "VHL",
        "pchembl_range": (7.0, 9.0),
        "notes": "VH032 scaffold; hydroxyproline mimetic essential for binding",
    },
    # IAP (Inhibitor of Apoptosis Proteins) binders
    "lcl161_smac_mimetic": {
        "smiles": "CC1CN(C(=O)c2cc(Cl)ccc2NC(=O)C2CCCN2C)CC1C",
        "e3_ligase": "cIAP1",
        "pchembl_range": (6.0, 7.5),
        "notes": "Smac mimetic; AVPI tetrapeptide analog",
    },
}

def get_e3_warhead(name: str) -> dict:
    """Retrieve E3 ligase warhead information by name."""
    if name not in E3_WARHEADS:
        raise KeyError(f"Unknown warhead '{name}'. Available: {list(E3_WARHEADS.keys())}")
    entry = E3_WARHEADS[name].copy()
    mol = Chem.MolFromSmiles(entry["smiles"])
    if mol:
        entry["mw"] = round(Chem.Descriptors.ExactMolWt(mol), 1)
        entry["n_heavy_atoms"] = mol.GetNumHeavyAtoms()
    return entry

Linker Enumeration and Geometry

The linker must span the distance between the E3 ligase surface and the POI surface in the productive ternary complex. Typical linker lengths for CRBN-based PROTACs are 8–20 atoms; for VHL-based PROTACs, shorter linkers (4–12 atoms) often suffice due to the smaller VHL–POI interface distance.

from rdkit import Chem
from rdkit.Chem import AllChem
import numpy as np
from itertools import product

# Linker building blocks available in PROTAC-DB and commercial libraries
LINKER_BLOCKS = {
    "peg1": ("OCC", "PEG unit 1"),
    "peg2": ("OCCO", "PEG unit 2"),
    "peg3": ("OCCCCO", "PEG unit 3"),
    "alkyl2": ("CC", "ethyl"),
    "alkyl3": ("CCC", "propyl"),
    "alkyl4": ("CCCC", "butyl"),
    "piperazine": ("N1CCNCC1", "piperazine ring"),
    "pip_peg": ("N1CCN(CC)CC1", "N-ethylpiperazine"),
    "rigid_phen": ("c1ccncc1", "pyridine spacer"),
    "amide": ("C(=O)N", "amide bond"),
}

def assemble_protac(
    poi_warhead_smiles: str,
    e3_warhead_smiles: str,
    linker_smiles: str,
    poi_attachment_smarts: str = "[NH2]",
    e3_attachment_smarts: str = "[NH2]",
) -> str | None:
    """
    Assemble a PROTAC from POI warhead, linker, and E3 warhead using amide coupling.

    This simplified version assumes both warheads have a free amine for linker coupling
    and the linker has carboxylic acids at both ends.

    Returns PROTAC SMILES or None if assembly failed.
    """
    # Amide coupling: ArNH2 + HOOC-linker-COOH + ArNH2
    # For a real pipeline use RDKit reaction SMARTS or DeLinker/PROTAC-DB tools
    rxn_smarts = "[C:1](=O)[OH].[NH2:2]>>[C:1](=O)[N:2]"
    rxn = AllChem.ReactionFromSmarts(rxn_smarts)

    try:
        poi_mol = Chem.MolFromSmiles(poi_warhead_smiles)
        e3_mol = Chem.MolFromSmiles(e3_warhead_smiles)
        linker_mol = Chem.MolFromSmiles(linker_smiles)
        if None in (poi_mol, e3_mol, linker_mol):
            return None

        # Step 1: couple linker COOH (end 1) to POI amine
        products_1 = rxn.RunReactants((linker_mol, poi_mol))
        if not products_1:
            return None
        intermediate = products_1[0][0]
        Chem.SanitizeMol(intermediate)

        # Step 2: couple remaining linker COOH (end 2) to E3 amine
        products_2 = rxn.RunReactants((intermediate, e3_mol))
        if not products_2:
            return None
        protac = products_2[0][0]
        Chem.SanitizeMol(protac)
        return Chem.MolToSmiles(protac)
    except Exception:
        return None


def estimate_linker_end_distance(
    e3_exit_vector: np.ndarray,
    poi_exit_vector: np.ndarray,
) -> dict:
    """
    Estimate required linker length from two exit vector coordinates.

    Args:
        e3_exit_vector: 3D coordinates of E3 warhead exit atom
        poi_exit_vector: 3D coordinates of POI warhead exit atom

    Returns:
        Dict with distance, suggested linker atom count, and linker type hints.
    """
    distance = float(np.linalg.norm(e3_exit_vector - poi_exit_vector))
    # Extended chain: ~1.25 Å per heavy atom in extended conformation
    n_atoms_extended = int(np.ceil(distance / 1.25))
    # Recommend slightly longer to allow geometric flexibility
    n_atoms_recommended = max(4, n_atoms_extended + 2)

    if n_atoms_recommended <= 6:
        linker_type = "short_alkyl_or_peg1"
    elif n_atoms_recommended <= 12:
        linker_type = "medium_peg_or_mixed_alkyl_peg"
    else:
        linker_type = "long_peg_or_piperazine_chain"

    return {
        "exit_vector_distance_A": round(distance, 2),
        "min_linker_atoms": n_atoms_extended,
        "recommended_linker_atoms": n_atoms_recommended,
        "linker_type_suggestion": linker_type,
    }

Ternary Complex Modelling and Cooperative Binding

The PROTAC must bridge two protein surfaces in a geometry that: (1) positions the POI's ubiquitination-acceptor lysines towards the E3 ligase's ubiquitin-transfer surface, (2) avoids steric clashes, and (3) has a sufficiently large burial area to give a cooperative ΔΔG.

The standard computational approach is:

1. Dock the POI warhead into the POI crystal structure (or predicted structure)
2. Dock the E3 warhead into the E3 crystal structure
3. Use protein-protein docking (ClusPro, HADDOCK, ZDOCK) to enumerate possible relative orientations of POI and E3 ligase
4. Score each orientation by: linker length compatibility + buried interface area + ubiquitin transfer site accessibility

import numpy as np
from scipy.spatial.distance import cdist

def score_ternary_geometry(
    poi_ub_acceptor_positions: np.ndarray,   # Lys residues on POI available for ubiquitination
    e3_ub_donor_position: np.ndarray,        # Active site Cys of E2/E3 machinery (3D coords)
    linker_end_distance: float,              # Estimated from docked PROTAC
    max_ub_transfer_distance: float = 12.0, # Max acceptable Lys-E2Cys distance for transfer
    buried_interface_area: float = 1000.0,  # Estimated from protein-protein docking
    min_interface_for_cooperativity: float = 800.0,
) -> dict:
    """
    Score the viability of a ternary complex geometry for productive ubiquitination.

    Returns a dictionary of geometric quality metrics.
    """
    # Distances from each POI lysine to E3 ubiquitin donor site
    poi_positions = poi_ub_acceptor_positions.reshape(-1, 3)
    e3_position = e3_ub_donor_position.reshape(1, 3)
    lys_distances = cdist(poi_positions, e3_position, metric="euclidean").flatten()

    n_reachable_lys = int((lys_distances <= max_ub_transfer_distance).sum())
    min_lys_dist = float(lys_distances.min())

    # Cooperativity: is the buried interface large enough?
    cooperative = buried_interface_area >= min_interface_for_cooperativity

    # Overall ternary complex score (0–1)
    reachability_score = min(n_reachable_lys / 3.0, 1.0)  # saturates at 3+ lysines
    interface_score = min(buried_interface_area / 1500.0, 1.0)
    linker_score = 1.0 if linker_end_distance < 25.0 else max(0.0, 1.0 - (linker_end_distance - 25.0) / 15.0)

    total_score = np.mean([reachability_score, interface_score, linker_score])

    return {
        "n_reachable_ubiquitination_lysines": n_reachable_lys,
        "min_lys_to_e3_dist_A": round(min_lys_dist, 2),
        "buried_interface_area_A2": round(buried_interface_area, 1),
        "is_cooperative": cooperative,
        "reachability_score": round(reachability_score, 3),
        "interface_score": round(interface_score, 3),
        "linker_geometry_score": round(linker_score, 3),
        "ternary_complex_score": round(float(total_score), 3),
    }

PROTAC-DB: Known PROTACs Dataset

Key Takeaways

• PROTACs operate catalytically — one molecule can degrade multiple copies of the POI — which enables sub-stoichiometric dosing and potentially overcoming drug resistance mutations in the binding site.
• Property filters for PROTACs differ substantially from Lipinski: MW 700–1100 Da, HBA ≤ 15, RotBonds 10–35. Standard Ro5 filtering will incorrectly reject all viable PROTACs.
• Ternary complex modelling (protein-protein pose sampling + linker geometry) is the key bottleneck. Tools like HTMD, PyDock, and HADDOCK can enumerate poses; score by: accessible Lys count, buried interface area, and linker end-to-end distance compatibility.
• In practice, a small focused PROTAC library (20–50 compounds varying linker length by 1–2 PEG/alkyl units) is synthesised and evaluated by SPR ternary binding + cellular degradation (Western blot / nanobret). Computational tools narrow this set, not replace it.
• Molecular glues (e.g., lenalidomide, indisulam) differ from PROTACs: a single small molecule simultaneously touches both E3 and POI without a linker, stabilising a neo-interface. ML models predicting molecular glue activity are an active research area (Molecular Glue DB, Magic Methylation).

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.