Reaction Enumeration

---
name: bio-reaction-enumeration
description: Enumerates virtual chemical libraries via reaction SMARTS transformations using RDKit and Reaction templates, with explicit handling of atom mapping, template extraction (RDKit reaction mining), product validation, RECAP/BRICS fragmentation, R-group decomposition, matched molecular pair analysis (MMPA), and Free-Wilson analysis. Use when generating combinatorial libraries from building blocks, enumerating analog series, deriving structure-activity rules, or extracting transformations from reaction data.
tool_type: python
primary_tool: RDKit
---

## Version Compatibility

Reference examples tested with: RDKit 2024.09+, mmpdb 3.1+, scikit-learn 1.4+, numpy 1.26+.

Before using code patterns, verify installed versions match. If versions differ:
- Python: `pip show <package>` then `help(module.function)` 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.

# Reaction Enumeration

Generate virtual libraries by applying reaction SMARTS to building blocks, enumerate analog series via matched molecular pairs, decompose into R-groups for SAR modeling, or extract transformations from reaction data. Reaction enumeration sits at the intersection of medicinal chemistry, lead optimization, and de novo design. The two key operations: **transform** (apply known rxn to make new compounds) and **mine** (extract rules from observed analog series). RDKit's reaction SMARTS handles the former; mmpdb / Free-Wilson handle the latter.

For retrosynthetic planning (target-to-starting-material decomposition), see `chemoinformatics/retrosynthesis`. For ML-driven design, see `chemoinformatics/generative-design`. For scaffold-based design, see `chemoinformatics/scaffold-analysis`.

## Operation Taxonomy

| Operation | Goal | Tool | Fails when |
|-----------|------|------|------------|
| Forward enumeration | Apply reaction to building blocks → products | RDKit `ReactionFromSmarts` + `RunReactants` | Wrong atom mapping; missing connectivity |
| Reverse enumeration (retrosynthesis) | Product → starting materials | AiZynthFinder, Chemformer | See retrosynthesis skill |
| Template mining | Reaction database → reaction SMARTS templates | RDKit reaction mining; rxnmapper | Atom mapping ambiguous; mechanism unclear |
| RECAP fragmentation | Molecule → retro-synthetic fragments | RDKit `Chem.Recap` | Inflexible bond rules |
| BRICS fragmentation | Molecule → retro-synthetic fragments | RDKit `BRICS` module | Many false fragments |
| R-group decomposition | Set of mols + scaffold → R-group table | RDKit `Chem.rdRGroupDecomposition` | Multiple scaffolds; ambiguous attachment |
| Matched Molecular Pairs (MMPA) | Set of mols → transformation rules | mmpdb | Need ≥1k compound dataset |
| Free-Wilson | Compounds + activities → additive R-group contributions | scikit-learn linear regression | Strict additivity assumption |

## Reaction SMARTS Basics

A reaction SMARTS is `reactants >> products` with atom maps `[atom:idx]` tracking atoms through the transformation:

```python
from rdkit.Chem import AllChem, Chem

amide = AllChem.ReactionFromSmarts(
    '[C:1](=[O:2])O.[N:3]>>[C:1](=[O:2])[N:3]'
)

errors = amide.Validate()
print(errors)
```

**Atom mapping rules:**
- Atoms with the same map index `[C:1]` in both reactant and product are tracked
- Maps must be unique within each reactant/product
- Unmapped atoms are added to or removed from the product
- Bond orders may change; map index preserves identity

**Common error:** Leaving an atom unmapped causes RDKit to either lose or duplicate it.

## Common Reaction Templates

```python
REACTIONS = {
    'amide_coupling': '[C:1](=[O:2])O.[N:3]>>[C:1](=[O:2])[N:3]',
    'reductive_amination': '[C:1](=O).[NH2:2]>>[CH:1][NH:2]',
    'suzuki': '[c:1][Br].[c:2][B](O)O>>[c:1][c:2]',
    'buchwald_hartwig': '[c:1][Br].[NH:2]>>[c:1][N:2]',
    'sn2_substitution': '[CH:1][Br].[N:2]>>[CH:1][N:2]',
    'sonogashira': '[c:1][Br].[CH:2]#[C:3]>>[c:1][C:2]#[C:3]',
    'click_chemistry': '[N-:1]=[N+:2]=[N:3][CH2:4].[CH:5]#[C:6]>>[N:3]1[N:2]=[N:1][C:6]=[C:5]1[CH2:4]',
    'esterification': '[C:1](=[O:2])O.[OH:3][C:4]>>[C:1](=[O:2])[O:3][C:4]',
    'urea_formation': '[N:1]=C=O.[NH:2]>>[N:1]C(=O)[N:2]',
    'sulfonamide': '[S:1](=O)(=O)Cl.[NH:2]>>[S:1](=O)(=O)[N:2]',
}
```

These are templates; real reactions need stereo, protecting-group, and chemoselectivity considerations. For production library enumeration, use validated templates from `rxnmapper` or vendor catalogs.

## Combinatorial Library Enumeration

**Goal:** Generate every (R1, R2, ..., Rn) product combination from sets of building blocks.

**Approach:** Cartesian product of reactant lists; apply reaction SMARTS; sanitize + deduplicate.

```python
from itertools import product
from rdkit import Chem
from rdkit.Chem import AllChem

def enumerate_library(rxn_smarts, reactant_lists, mw_max=600):
    rxn = AllChem.ReactionFromSmarts(rxn_smarts)
    if rxn.Validate()[0] != 0:
        raise ValueError(f'Invalid reaction: {rxn_smarts}')

    seen = set()
    products = []
    for combo in product(*reactant_lists):
        mols = [Chem.MolFromSmiles(s) for s in combo]
        if None in mols:
            continue

        for prod_tuple in rxn.RunReactants(tuple(mols)):
            for prod in prod_tuple:
                try:
                    Chem.SanitizeMol(prod)
                    smi = Chem.MolToSmiles(prod)
                    if smi in seen:
                        continue
                    if Chem.Descriptors.MolWt(prod) > mw_max:
                        continue
                    seen.add(smi)
                    products.append(smi)
                except Exception:
                    continue
    return products
```

**Scaling:** For a 100x100x100 enumeration (1M products), parallelize with multiprocessing. For 1k x 1k x 1k (1B products), use a streaming approach + filter before materializing.

## RECAP Fragmentation

RECAP (Lewell 1998) breaks molecules at retrosynthetically reasonable bonds into reusable fragments.

```python
from rdkit.Chem import Recap

mol = Chem.MolFromSmiles('c1ccc(C(=O)Nc2ccc(F)cc2)cc1')
hier = Recap.RecapDecompose(mol)
fragments = list(hier.GetLeaves().keys())
```

RECAP bond types: amide, ester, ether, amine, urea, olefin, quaternary nitrogen, sulfonamide. Use cases: building-block library generation, scaffold-decoration enumeration.

## BRICS Fragmentation

BRICS (Degen 2008) is an extension of RECAP with more bond types. Better fragment coverage; more fragments per molecule.

```python
from rdkit.Chem import BRICS

mol = Chem.MolFromSmiles('CCN(CC)c1ccc(C(=O)NC2CCCC2)cc1')
fragments = BRICS.BRICSDecompose(mol)

builder = BRICS.BRICSBuild([Chem.MolFromSmiles(f) for f in fragments])
new_mols = [next(builder) for _ in range(10)]
```

`BRICSDecompose` produces SMILES with `[<dummy>]` attachment points; `BRICSBuild` recombines fragments at these dummies.

## R-Group Decomposition

**Goal:** Given a set of compounds sharing a scaffold, extract the R-group at each attachment point into a tabular SAR matrix.

**Approach:** Define scaffold with `[*:1]`, `[*:2]` placeholders; RDKit matches each compound and extracts R-groups.

```python
from rdkit.Chem import rdRGroupDecomposition as rgd
from rdkit import Chem

scaffold = Chem.MolFromSmiles('c1ccc(-[*:1])cc1-[*:2]')

mols = [Chem.MolFromSmiles(smi) for smi in [
    'c1ccc(C)cc1F',
    'c1ccc(CC)cc1Cl',
    'c1ccc(CCC)cc1Br',
]]

decomp, _ = rgd.RGroupDecompose([scaffold], mols, asSmiles=True)
```

`decomp` is a list of dicts `{'Core': scaffold_smi, 'R1': r1_smi, 'R2': r2_smi}`. Combined with activity column, enables Free-Wilson.

## Matched Molecular Pairs Analysis (MMPA)

MMPA (Hussain & Rea 2010) extracts SAR rules from compound pairs differing by a single transformation.

```bash
mmpdb fragment data.smi -o data.fragments
mmpdb index data.fragments -o data.mmpdb
mmpdb transform --smiles 'COc1ccccc1' data.mmpdb
```

`mmpdb` produces a database of transformations + statistics on activity changes.

| Transformation | Avg delta(pIC50) | N pairs | Confidence |
|----------------|-------------------|---------|------------|
| Me → F | +0.5 | 152 | high |
| OMe → OH | -0.3 | 89 | moderate |
| Ph → 4-pyridine | +1.2 | 23 | moderate |

**Use case:** Lead optimization. Given a hit, ask "what transformations have improved similar series?" Apply top-ranked transformations to generate analog suggestions.

**Context-based MMPA** (Awale 2024): condition rules on local chemical context (e.g., "Me→F adjacent to amide"). Outperforms classical MMPA on CYP1A2 inhibition reduction.

## Free-Wilson Analysis

**Goal:** Decompose activity into additive R-group contributions.

**Approach:** Linear regression with R-group identity as binary features.

```python
import pandas as pd
from sklearn.linear_model import Ridge

def free_wilson(decomp_results, activity_col='pIC50'):
    df = pd.DataFrame(decomp_results)
    r_groups = pd.get_dummies(df[['R1', 'R2']], prefix=['R1', 'R2'])
    X = r_groups.values
    y = df[activity_col].values
    model = Ridge(alpha=0.1).fit(X, y)
    contributions = dict(zip(r_groups.columns, model.coef_))
    return contributions, model.intercept_
```

**Trade-off:** Free-Wilson assumes additivity (R1 contribution independent of R2). Real SAR has interactions; Free-Wilson predictions for un-synthesized combinations are biased when synergy exists. Use as a *first-pass model* for analog prioritization; validate with QSAR.

## Template Extraction from Reaction Data

**Goal:** Given an atom-mapped reaction SMILES, extract a generalizable SMARTS template.

**Approach:** Use `rxnmapper` (Schwaller 2021) for atom mapping, then RDKit reaction template extraction.

```python
from rxnmapper import RXNMapper

mapper = RXNMapper()
rxns = ['CCO.OC(=O)c1ccccc1>>CCOC(=O)c1ccccc1']
results = mapper.get_attention_guided_atom_maps(rxns)
mapped_smiles = results[0]['mapped_rxn']
```

After atom-mapping, RDKit can extract a template via `ChemicalReaction.GetReactionTemplateFromMappedReaction` (custom implementation; see Coley 2019).

## Per-Tool Failure Modes

### Reaction SMARTS -- atom mapping mismatch

**Trigger:** Map index appears on reactant but not product.

**Mechanism:** RDKit treats unmapped atoms as deleted from product; an atom that was meant to be preserved disappears if its map index is missing.

**Symptom:** Products missing expected atoms; valences wrong; sanitize fails.

**Fix:** Validate with `rxn.Validate()`; manually inspect mapping; use Reaction Atom Mapping Number column 2.

### RECAP/BRICS -- over-fragmentation

**Trigger:** Highly substituted molecule with many breakable bonds.

**Mechanism:** Default bond list breaks at every retrosynthetic position; one molecule yields tens of fragments.

**Symptom:** Building-block enumeration explodes; many small irrelevant fragments.

**Fix:** Filter fragments by MW (>=80 Da), heavy atom count (>=4); use only meaningful fragments downstream.

### MMPA -- insufficient pair count

**Trigger:** mmpdb on small dataset (<500 compounds, <50 actives).

**Mechanism:** MMPA needs ≥10 pairs per transformation to yield a statistically meaningful delta(activity).

**Symptom:** Transformations report with N=1-3 pairs; effect sizes erratic.

**Fix:** Filter to transformations with N>=10; supplement with literature SAR knowledge.

### Free-Wilson -- non-additive interactions

**Trigger:** R1 and R2 interact through hydrogen bonding, steric clash, or electronic effects.

**Mechanism:** Free-Wilson is purely additive; cannot capture R1+R2 synergy.

**Symptom:** Predicted activities for un-synthesized combinations are biased low for synergistic pairs.

**Fix:** Use Free-Wilson as first-pass screen; validate predictions with QSAR (random forest, chemprop) which captures interactions.

### R-group decomposition -- multiple scaffolds

**Trigger:** Compound matches multiple scaffold templates.

**Mechanism:** `RGroupDecompose` returns the first matching scaffold; ambiguous SAR series.

**Symptom:** Same compound's R-groups differ between runs.

**Fix:** Specify scaffold unambiguously; use only compounds matching one scaffold.

### Reaction enumeration -- combinatorial explosion

**Trigger:** Large building-block sets (1k x 1k = 1M products).

**Mechanism:** Cartesian product * RunReactants is O(N^d) where d is reactant count.

**Symptom:** Memory blowup, multi-hour runtime.

**Fix:** Pre-filter building blocks; stream products to file rather than list; use mmpdb-style sparse enumeration only for valid pairings.

## Reconciliation: Free-Wilson vs MMPA

Both methods derive R-group rules but from different perspectives:
- Free-Wilson: linear regression on assembled SAR table; gives R-group contributions
- MMPA: transformation-based; gives delta(activity) for each substitution

If they agree on direction (Me→F improves activity), high confidence. If they disagree, investigate non-additive interactions or look for context dependence in MMPA.

## Common Errors

| Symptom | Cause | Fix |
|---------|-------|-----|
| `rxn.Validate()` returns errors | Bad atom mapping or invalid SMARTS | Re-check map indices; valences |
| Products contain unexpected fragments | Reactants matched in unintended way | Use more specific SMARTS; constrain with explicit ring members |
| Sanitize fails on products | Reaction breaks valence | Filter via `Chem.SanitizeMol(prod, catchErrors=True)` |
| Duplicate products | Same product from different reactant orientations | Deduplicate by canonical SMILES |
| RECAP produces single fragment | Molecule has no retrosynthetic bonds | Try BRICS for more aggressive fragmentation |
| mmpdb empty output | Insufficient dataset size or no matched pairs | Need >=1000 compounds |
| R-group decomposition wrong R | Scaffold dummy not aligned | Re-check `[*:1]` / `[*:2]` placement |

## References

- Hartenfeller et al., *J. Cheminformatics* 4:38 (2012) -- DOGS rule-based library design.
- Lewell et al., *J. Chem. Inf. Comput. Sci.* 38:511 (1998) -- RECAP.
- Degen et al., *ChemMedChem* 3:1503 (2008) -- BRICS fragmentation.
- Hussain & Rea, *J. Chem. Inf. Model.* 50:339 (2010) -- MMPA.
- Dossetter et al., *Drug Discov. Today* 18:724 (2013) -- Practical MMPA in lead optimization.
- Free & Wilson, *J. Med. Chem.* 7:395 (1964) -- Original Free-Wilson.
- Schwaller et al., *Sci. Adv.* 7:eabe4166 (2021) -- rxnmapper.
- Coley et al., *Chem. Sci.* 10:370 (2019) -- Reaction template extraction.

## Related Skills

- chemoinformatics/molecular-io - Read/write reaction SMILES
- chemoinformatics/substructure-search - SMARTS pattern matching
- chemoinformatics/scaffold-analysis - Bemis-Murcko scaffolds for R-decomp
- chemoinformatics/molecular-descriptors - Featurize products
- chemoinformatics/admet-prediction - Filter enumerated products
- chemoinformatics/retrosynthesis - Reverse direction (target → starting materials)
- chemoinformatics/generative-design - Generative alternatives to template enumeration
- chemoinformatics/qsar-modeling - Validate Free-Wilson predictions