Retrosynthesis and Reaction Prediction

Every compound generated by a molecular design algorithm must eventually be synthesised in a lab. Synthesis planning — working backwards from a target molecule to commercially available starting materials through a sequence of chemical transformations — is one of the most complex tasks in applied chemistry. ML-based retrosynthesis tools automate this reasoning, enabling medicinal chemists to rapidly assess the synthetic feasibility of thousands of virtual candidates generated by generative models or virtual screening.

This guide covers the two paradigms of computational retrosynthesis (template-based and template-free), how to use AiZynthFinder for end-to-end route planning, how Chemformer performs template-free retrosynthesis as a SMILES-to-SMILES translation problem, and how forward reaction prediction complements retrosynthesis for validating predicted routes.

The Retrosynthesis Problem

Given a target molecule (product), retrosynthesis identifies precursors (reactants) and the reaction type that could produce it. This is applied recursively until all precursors are commercially available or in stock. The problem is challenging because:

• A large drug-like molecule may have hundreds of possible disconnection points
• Not all chemically valid retrosynthetic steps are practically feasible (yield, selectivity, cost)
• Multi-step routes grow exponentially — a 5-step synthesis has millions of possible route trees
• Chemical literature has ~20 million reactions in USPTO, but only a fraction are drug-relevant

Retrosynthesis Approaches

AiZynthFinder: Monte Carlo Tree Search for Route Planning

AiZynthFinder (AstraZeneca, open source) uses a template-based single-step model combined with Monte Carlo Tree Search (MCTS) to plan multi-step synthesis routes. The workflow:

1. A neural network (feedforward or GNN) scores the applicability of each SMARTS reaction template to the current target molecule.
2. MCTS explores the route tree, expanding nodes (applying retrosynthetic steps) and backpropagating scores.
3. The search terminates when all leaves are commercially available compounds (in stock filter) or a maximum depth/node limit is reached.
4. Routes are ranked by a configurable scoring function (route length, stock availability, reaction feasibility).

"""
Install: pip install aizynthfinder
Pre-trained models available at: https://github.com/MolecularAI/aizynthfinder
"""
from aizynthfinder.aizynthfinder import AiZynthFinder
import json

def plan_synthesis(
    target_smiles: str,
    config_path: str = "aizynthfinder_config.yml",
    iteration_limit: int = 100,
) -> list[dict]:
    """
    Plan synthetic routes for a target molecule using AiZynthFinder.

    Args:
        target_smiles:   SMILES of the molecule to synthesise
        config_path:     Path to AiZynthFinder YAML config (model weights + stock)
        iteration_limit: MCTS search iterations

    Returns:
        List of route dicts with keys: route_id, num_steps, in_stock_fraction, score
    """
    finder = AiZynthFinder(configfile=config_path)
    finder.target_smiles = target_smiles
    finder.config.iteration_limit = iteration_limit

    finder.tree_search()
    finder.build_routes()

    routes = []
    for i, route in enumerate(finder.routes):
        routes.append({
            "route_id": i + 1,
            "num_steps": route.metadata.get("number_of_steps", "?"),
            "in_stock_fraction": round(route.metadata.get("in_stock_fraction", 0), 2),
            "score": round(route.metadata.get("score", 0), 3),
            "route_dict": route.to_dict(),
        })

    routes.sort(key=lambda r: (-r["in_stock_fraction"], r["num_steps"]))
    return routes


# Example usage
target = "CC(Cc1ccc(F)cc1)NC(=O)c1ccc(Cl)cc1"   # Arbitrary drug-like molecule
routes = plan_synthesis(target)

print(f"Found {len(routes)} routes for target molecule")
for r in routes[:5]:
    print(
        f"  Route {r['route_id']:2d}: {r['num_steps']} steps, "
        f"in-stock {r['in_stock_fraction']:.0%}, score {r['score']:.3f}"
    )

Template-Free Retrosynthesis with Chemformer

Chemformer (AstraZeneca, 2022) is a BART-based seq2seq model trained on 100 million reactions from USPTO and ZINC. For retrosynthesis, it is trained to translate a product SMILES → reactant SMILES. For forward prediction, it translates reactants → product.

from transformers import AutoTokenizer, AutoModelForSeq2SeqLM
import torch

CHEMFORMER_MODEL = "MolecularAI/chemformer-base-10M"

tokenizer = AutoTokenizer.from_pretrained(CHEMFORMER_MODEL)
model = AutoModelForSeq2SeqLM.from_pretrained(CHEMFORMER_MODEL)
model.eval()

def predict_retrosynthesis(
    product_smiles: str,
    num_beams: int = 10,
    num_results: int = 5,
) -> list[dict]:
    """
    Predict retrosynthetic precursors for a product SMILES using Chemformer.
    Uses beam search to generate top-k precursor predictions.

    Returns:
        List of dicts with keys: rank, precursors_smiles, score
    """
    # Chemformer uses special SMILES tokens; wrap product with [START] [END]
    inputs = tokenizer(
        product_smiles,
        return_tensors="pt",
        max_length=512,
        truncation=True,
    )

    with torch.no_grad():
        outputs = model.generate(
            **inputs,
            num_beams=num_beams,
            num_return_sequences=num_results,
            max_length=512,
            early_stopping=True,
            output_scores=True,
            return_dict_in_generate=True,
        )

    sequences = tokenizer.batch_decode(outputs.sequences, skip_special_tokens=True)
    # Compute token-level sequence scores
    scores = outputs.sequences_scores.exp().tolist() if hasattr(outputs, "sequences_scores") else [None] * len(sequences)

    results = []
    for i, (seq, score) in enumerate(zip(sequences, scores)):
        results.append({
            "rank": i + 1,
            "precursors_smiles": seq.strip(),
            "confidence": round(score, 4) if score is not None else None,
        })
    return results


# Example
product = "CC(=O)Oc1ccccc1C(=O)O"   # Aspirin
retros = predict_retrosynthesis(product)
print(f"Top retrosynthetic precursors for {product}:")
for r in retros:
    print(f"  {r['rank']}. {r['precursors_smiles']}  (confidence: {r['confidence']})")

Forward Reaction Prediction

Forward prediction answers "given these reactants, what product will form?" It is used to validate retrosynthetic routes: apply the predicted retro step forward and check that the product is recovered. It is also used for reaction outcome prediction in medicinal chemistry.

def predict_forward_reaction(
    reactants_smiles: str,
    reagent_smiles: str = "",
    num_results: int = 3,
) -> list[dict]:
    """
    Predict reaction products from reactant SMILES using Chemformer.

    Input format: 'reactant1.reactant2>reagent>product' (SMILES reaction format)
    For Chemformer, we encode as: 'reactant1.reactant2' (concatenated by period)

    Returns:
        List of dicts with predicted product SMILES and confidence scores.
    """
    # Standard reaction SMILES format: reactants>>product
    # For Chemformer forward prediction, input is reactant SMILES string
    input_smiles = reactants_smiles
    if reagent_smiles:
        input_smiles = f"{reactants_smiles}.{reagent_smiles}"

    inputs = tokenizer(
        input_smiles,
        return_tensors="pt",
        max_length=512,
        truncation=True,
    )

    with torch.no_grad():
        outputs = model.generate(
            **inputs,
            num_beams=num_results * 2,
            num_return_sequences=num_results,
            max_length=256,
            early_stopping=True,
            output_scores=True,
            return_dict_in_generate=True,
        )

    sequences = tokenizer.batch_decode(outputs.sequences, skip_special_tokens=True)
    scores = outputs.sequences_scores.exp().tolist() if hasattr(outputs, "sequences_scores") else [None] * len(sequences)

    return [
        {"rank": i + 1, "product_smiles": seq.strip(), "confidence": round(sc, 4) if sc else None}
        for i, (seq, sc) in enumerate(zip(sequences, scores))
    ]

Reaction Datasets

The quality of retrosynthesis models depends heavily on training data:

Evaluating Retrosynthesis Models

The standard evaluation metric on USPTO-50K is top-k exact match: the fraction of test products for which the correct reactant set appears within the top-k predictions. "Correct" means the predicted reactant SMILES, when canonicalized, matches the ground truth.

Synthetic Accessibility Scoring

Before running a full retrosynthesis, a fast synthetic accessibility score can triage compounds that are likely difficult to make. The two most common scores are:

• SA score (Ertl & Schuffenhauer): based on fragment contributions from a PubChem fragment database. Score 1 (easy) to 10 (hard). Available in RDKit.
• SCScore (Coley et al., 2018): neural network trained on reaction products vs. reactants in USPTO. Outputs 1–5. Higher = requires more synthesis steps.

from rdkit import Chem
from rdkit.Chem import RDConfig
import sys
import os

# SA score is in the RDKit Contrib folder
sys.path.append(os.path.join(RDConfig.RDContribDir, 'SA_Score'))
import sascorer

def compute_sa_score(smiles: str) -> float:
    """
    Compute Synthetic Accessibility score (SA score) for a molecule.
    Range: 1 (easy to synthesise) to 10 (very difficult).
    Drug-like compounds typically have SA score < 3.5.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None
    return round(sascorer.calculateScore(mol), 2)


# Example
test_compounds = [
    ("CC(=O)Oc1ccccc1C(=O)O", "Aspirin"),
    ("CN1C=NC2=C1C(=O)N(C(=O)N2C)C", "Caffeine"),
    ("CC12CCC3C(C1CCC2O)CCC4=CC(=O)CCC34C", "Testosterone"),
]

for smiles, name in test_compounds:
    sa = compute_sa_score(smiles)
    print(f"{name}: SA score = {sa}")

Integrating Retrosynthesis into a Design-Make-Test Cycle

Retrosynthesis tools are most valuable when integrated into a closed-loop design cycle:

1. Design: generate virtual candidates using a generative model (VAE, REINVENT) or virtual screening.
2. Triage: filter by SA score (< 3.5) and ADMET predictions to remove unsynthesisable or unsuitable compounds.
3. Route planning: run AiZynthFinder on the filtered candidates; prioritise compounds for which ≥1 route with all in-stock starting materials is found within 3 steps.
4. Route validation: run Chemformer forward prediction on each proposed route; confirm the product SMILES is recovered (round-trip validation).
5. Make: the medicinal chemist synthesises the top-ranked compound, often with the assistance of a robot synthesis platform (Chemspeed, Chemify).
6. Test: measure biological activity, feed back to the model.

Key Takeaways

• Template-based tools like AiZynthFinder find routes in seconds using MCTS and are chemically interpretable; template-free seq2seq models like Chemformer can suggest novel reaction types but at lower precision.
• Always filter virtual compound libraries by SA score (< 3.5) before retrosynthesis to avoid wasting compute on inherently unsynthesisable structures.
• Round-trip validation (retro → forward) filters hallucinated routes before they reach the chemist. A route is only reliable if the Chemformer forward model recovers the target molecule.
• The AiZynthFinder MCTS depth limit (typically 5–6 steps) and stock database coverage are the main constraints on finding viable routes. Extend stock with Enamine REAL, eMolecules, or custom in-house inventories.
• Retrosynthesis models are trained on literature reactions (mostly USPTO); their performance degrades for reactions outside standard organic chemistry (e.g., biocatalysis, exotic metal catalysis).

Try it live

Explore a related model hands-on: open AI Retrosynthesis Planner: Molecular Disconnection Search and change its parameters while it runs. Everything stays in your browser.