Stereochemistry in Drug Discovery

Stereochemistry is critical in drug development - enantiomers can have dramatically different:

• Pharmacological activity (thalidomide: one enantiomer teratogenic)

• Potency (often 10-1000x difference)

• Metabolism (different CYP450 interactions)

• Toxicity (stereoselective off-target effects)

This notebook covers:

1. Stereoisomer enumeration - generating all possible stereoisomers

2. Chiral center analysis - identifying and annotating stereocenters

3. 3D structure generation - maintaining stereochemistry

4. Enantiomer comparison - energy and geometry differences

Chemistry Background

Types of Stereoisomers

Type	Relationship	Example
Enantiomers	Mirror images, non-superimposable	R/S alanine
Diastereomers	Multiple stereocenters, not mirror images	threonine vs allo-threonine
Epimers	Differ at one stereocenter	α/β-glucose
Cis/Trans	Geometric isomers around double bonds	cis/trans-stilbene

import os
import tempfile
from pathlib import Path

import Auto3D
from Auto3D import Auto3DOptions, main
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors, rdMolDescriptors
from rdkit.Chem import Draw
from rdkit.Chem.EnumerateStereoisomers import EnumerateStereoisomers, StereoEnumerationOptions
import pandas as pd

print(f"Auto3D version: {Auto3D.__version__}")

1. Famous Stereochemistry Examples in Drugs

These examples illustrate why stereochemistry matters in drug development.

# Classic examples where stereochemistry determines activity
stereo_examples = {
    # Thalidomide - one of the most notorious examples
    # (R)-thalidomide: sedative
    # (S)-thalidomide: teratogenic
    "thalidomide_R": "O=C1CC[C@H](N2C(=O)C3=CC=CC=C3C2=O)C(=O)N1",
    "thalidomide_S": "O=C1CC[C@@H](N2C(=O)C3=CC=CC=C3C2=O)C(=O)N1",

    # Ibuprofen - only S-enantiomer is active
    "ibuprofen_S": "CC(C)CC1=CC=C(C=C1)[C@H](C)C(=O)O",
    "ibuprofen_R": "CC(C)CC1=CC=C(C=C1)[C@@H](C)C(=O)O",

    # Omeprazole vs Esomeprazole (S-omeprazole)
    "esomeprazole_S": "COC1=CC2=NC([C@H](C)S(=O)C3=NC4=CC=CC=C4N3)=NC2=CC1OC",

    # Methadone - R-enantiomer 10x more potent
    "methadone_R": "CC[C@@](C)(C(=O)CC)C(C1=CC=CC=C1)C2=CC=CC=C2",

    # Propranolol - S-enantiomer 100x more potent at β-receptors
    "propranolol_S": "CC(C)NC[C@H](O)COC1=CC=CC2=CC=CC=C12",
}

print("Stereospecific drugs:")
for name, smi in stereo_examples.items():
    mol = Chem.MolFromSmiles(smi)
    n_stereo = len(Chem.FindMolChiralCenters(mol, includeUnassigned=True))
    print(f"  {name}: {n_stereo} chiral center(s)")

2. Identifying Chiral Centers

Before generating stereoisomers, we need to identify all stereocenters.

def analyze_stereocenters(smiles):
    """
    Analyze stereochemistry of a molecule.

    Returns:
    - Number of defined stereocenters
    - Number of undefined stereocenters
    - Maximum possible stereoisomers (2^n)
    """
    mol = Chem.MolFromSmiles(smiles)

    # Find all chiral centers
    chiral_centers = Chem.FindMolChiralCenters(mol, includeUnassigned=True)

    defined = sum(1 for _, tag in chiral_centers if tag in ['R', 'S'])
    undefined = sum(1 for _, tag in chiral_centers if tag == '?')

    # Check for double bond stereochemistry
    double_bond_stereo = 0
    for bond in mol.GetBonds():
        if bond.GetBondType() == Chem.BondType.DOUBLE:
            stereo = bond.GetStereo()
            if stereo in [Chem.BondStereo.STEREOE, Chem.BondStereo.STEREOZ,
                          Chem.BondStereo.STEREOANY]:
                double_bond_stereo += 1

    total_undefined = undefined + double_bond_stereo
    max_isomers = 2 ** (defined + undefined + double_bond_stereo)

    return {
        "defined_centers": defined,
        "undefined_centers": undefined,
        "double_bond_stereo": double_bond_stereo,
        "max_stereoisomers": max_isomers,
        "chiral_centers": chiral_centers
    }


# Analyze some drug molecules
drugs_to_analyze = {
    "atorvastatin": "CC(C)C1=C(C(=C(N1CCC(CC(CC(=O)O)O)O)C2=CC=C(C=C2)F)C3=CC=CC=C3)C(=O)NC4=CC=CC=C4",
    "simvastatin": "CC[C@H](C)C(=C)[C@H]1CC[C@@H]2C1(CC[C@H](C2)OC(=O)C)C",
    "penicillin_G": "CC1([C@@H](N2[C@H](S1)[C@@H](C2=O)NC(=O)CC3=CC=CC=C3)C(=O)O)C",
}

print("Stereochemistry analysis:")
print("-" * 70)
for name, smi in drugs_to_analyze.items():
    info = analyze_stereocenters(smi)
    print(f"{name}:")
    print(f"  Defined centers: {info['defined_centers']}")
    print(f"  Undefined centers: {info['undefined_centers']}")
    print(f"  Max stereoisomers: {info['max_stereoisomers']}")
    print(f"  Chiral centers: {info['chiral_centers']}")
    print()

3. Enumerating Stereoisomers

For molecules with undefined stereochemistry, we need to enumerate all possible isomers.

def enumerate_stereoisomers(smiles, max_isomers=16):
    """
    Enumerate all stereoisomers of a molecule.

    Args:
        smiles: Input SMILES (can have undefined stereochemistry)
        max_isomers: Maximum number of isomers to generate

    Returns:
        List of (SMILES, R/S assignment) tuples
    """
    mol = Chem.MolFromSmiles(smiles)

    # Configure enumeration options
    opts = StereoEnumerationOptions(
        tryEmbedding=True,  # Verify 3D embedding is possible
        unique=True,        # Remove duplicates
        maxIsomers=max_isomers
    )

    isomers = list(EnumerateStereoisomers(mol, options=opts))

    results = []
    for iso in isomers:
        smi = Chem.MolToSmiles(iso, isomericSmiles=True)
        centers = Chem.FindMolChiralCenters(iso)
        config = "/".join([f"{idx}:{tag}" for idx, tag in centers])
        results.append((smi, config))

    return results


# Enumerate stereoisomers for a molecule with undefined centers
# Ephedrine - 2 chiral centers = 4 possible stereoisomers
ephedrine_undefined = "CC(C(C1=CC=CC=C1)O)NC"  # No stereochemistry defined

isomers = enumerate_stereoisomers(ephedrine_undefined)
print(f"Ephedrine stereoisomers ({len(isomers)} found):")
print("-" * 60)
for i, (smi, config) in enumerate(isomers):
    print(f"Isomer {i+1}: {config}")
    print(f"  SMILES: {smi}")

4. 3D Structure Generation Preserving Stereochemistry

Auto3D correctly generates 3D structures that preserve the input stereochemistry.

# Generate 3D for all stereoisomers
stereoisomer_names = {
    "1R,2S_ephedrine": "C[C@@H]([C@H](C1=CC=CC=C1)O)NC",  # (-)-ephedrine
    "1S,2R_ephedrine": "C[C@H]([C@@H](C1=CC=CC=C1)O)NC",  # (+)-ephedrine
    "1R,2R_pseudoephedrine": "C[C@@H]([C@@H](C1=CC=CC=C1)O)NC",  # (+)-pseudoephedrine
    "1S,2S_pseudoephedrine": "C[C@H]([C@H](C1=CC=CC=C1)O)NC",  # (-)-pseudoephedrine
}

# Note: ephedrine and pseudoephedrine are diastereomers!
# They have different physical properties (melting point, solubility)

with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:
    for name, smi in stereoisomer_names.items():
        f.write(f"{smi} {name}\n")
    stereo_file = f.name

print(f"Wrote {len(stereoisomer_names)} stereoisomers to {stereo_file}")

# Run Auto3D - it will preserve stereochemistry
if __name__ == "__main__":
    config = Auto3DOptions(
        path=stereo_file,
        k=1,
        optimizing_engine="AIMNET",
        use_gpu=True,
        # enumerate_isomer=False,  # Don't enumerate - we already have all isomers
    )

    stereo_output = main(config)
    print(f"Output: {stereo_output}")

5. Comparing Enantiomer Energies

In vacuum, enantiomers have identical energies. However:

• Diastereomers have different energies

• In chiral environments (proteins, enzymes), enantiomers interact differently

def compare_stereoisomer_energies(sdf_path):
    """
    Compare energies of stereoisomers.

    Enantiomers should have identical energies in vacuum.
    Diastereomers will have different energies.
    """
    mols = list(Chem.SDMolSupplier(sdf_path, removeHs=False))

    data = []
    for mol in mols:
        if mol is None:
            continue

        name = mol.GetProp("_Name")

        # Get energy
        if mol.HasProp("E_hartree"):
            e = float(mol.GetProp("E_hartree"))
        elif mol.HasProp("E_tot"):
            e = float(mol.GetProp("E_tot")) / 27.2114
        else:
            continue

        # Get stereochemistry
        centers = Chem.FindMolChiralCenters(mol)
        stereo = "/".join([tag for _, tag in centers])

        data.append({
            "Name": name,
            "Stereo": stereo,
            "E_hartree": e
        })

    df = pd.DataFrame(data)
    df["E_kcal"] = (df["E_hartree"] - df["E_hartree"].min()) * 627.509

    return df.sort_values("E_kcal")


if 'stereo_output' in dir():
    df_stereo = compare_stereoisomer_energies(stereo_output)
    print("Stereoisomer Energy Comparison:")
    print(df_stereo[["Name", "Stereo", "E_kcal"]].to_string(index=False))
    print("\nNote: Enantiomer pairs should have ~identical energies.")
    print("Diastereomers (ephedrine vs pseudoephedrine) have different energies.")

6. Double Bond Stereochemistry (E/Z)

Geometric isomers around double bonds are also important in drug design.

# Examples of E/Z isomerism in drugs
ez_examples = {
    # Tamoxifen - E and Z have opposite effects!
    "tamoxifen_Z": "CC/C(=C(\\C1=CC=CC=C1)/C2=CC=C(OCCN(C)C)C=C2)/C3=CC=CC=C3",
    "tamoxifen_E": "CC/C(=C(/C1=CC=CC=C1)\\C2=CC=C(OCCN(C)C)C=C2)/C3=CC=CC=C3",

    # Resveratrol - trans form is biologically active
    "resveratrol_trans": "OC1=CC=C(/C=C/C2=CC(O)=CC(O)=C2)C=C1",
    "resveratrol_cis": "OC1=CC=C(/C=C\\C2=CC(O)=CC(O)=C2)C=C1",

    # Retinoids - critical for vitamin A function
    "retinoic_acid_all_trans": "CC1=C(C(CCC1)(C)C)/C=C/C(C)=C/C=C/C(C)=C/C(=O)O",
}

with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:
    for name, smi in ez_examples.items():
        f.write(f"{smi} {name}\n")
    ez_file = f.name

print(f"E/Z isomers: {len(ez_examples)}")

# Generate 3D structures for E/Z isomers
if __name__ == "__main__":
    config = Auto3DOptions(
        path=ez_file,
        k=1,
        optimizing_engine="AIMNET",
        use_gpu=True,
    )

    ez_output = main(config)
    print(f"Output: {ez_output}")

# Compare E/Z isomer energies
if 'ez_output' in dir():
    df_ez = compare_stereoisomer_energies(ez_output)
    print("E/Z Isomer Energy Comparison:")
    print(df_ez[["Name", "E_kcal"]].to_string(index=False))
    print("\nNote: E and Z isomers have different energies (unlike enantiomers).")

7. Auto3D Stereoisomer Enumeration

Auto3D can automatically enumerate stereoisomers during conformer generation.

# Start with undefined stereochemistry
undefined_stereo = {
    "drug_candidate": "CC(NC(=O)C(CC1=CC=CC=C1)NC(=O)OC(C)(C)C)C(=O)O",
}

with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:
    for name, smi in undefined_stereo.items():
        f.write(f"{smi} {name}\n")
    undefined_file = f.name

# Count stereocenters
info = analyze_stereocenters(undefined_stereo["drug_candidate"])
print(f"Drug candidate has {info['defined_centers'] + info['undefined_centers']} stereocenters")
print(f"Maximum possible stereoisomers: {info['max_stereoisomers']}")

# Run with stereoisomer enumeration
if __name__ == "__main__":
    config = Auto3DOptions(
        path=undefined_file,
        k=1,                        # Top conformer per stereoisomer
        enumerate_isomer=True,      # Enumerate stereoisomers
        optimizing_engine="AIMNET",
        use_gpu=True,
    )

    enum_output = main(config)
    print(f"Output: {enum_output}")

# Analyze enumerated stereoisomers
if 'enum_output' in dir():
    mols = list(Chem.SDMolSupplier(enum_output, removeHs=False))
    print(f"Generated {len(mols)} stereoisomers")

    df_enum = compare_stereoisomer_energies(enum_output)
    print("\nStereoisomer Ranking by Energy:")
    print(df_enum[["Name", "Stereo", "E_kcal"]].to_string(index=False))

8. Best Practices for Stereochemistry in Drug Discovery

Lead Optimization

1. Early determination: Identify the active stereoisomer early

2. Synthesize pure isomers: Test individual stereoisomers, not racemates

3. Consider racemization: Some drugs racemize in vivo (ibuprofen)

Virtual Screening

1. Define stereochemistry: Use explicit @/@@/E/Z notation

2. Enumerate if unknown: Use Auto3D’s enumerate_isomer=True

3. Dock all isomers: Different stereoisomers may have different binding modes

QSAR/ML

1. Consistent representation: Use canonical SMILES with stereochemistry

2. 3D descriptors: Capture stereo-dependent properties

3. Chiral descriptors: Consider adding chirality-specific features

# Cleanup
for f in [stereo_file, ez_file, undefined_file]:
    if f in dir() and os.path.exists(f):
        os.unlink(f)

Summary

This tutorial covered:

1. Stereochemistry fundamentals - enantiomers, diastereomers, E/Z isomers

2. Chiral center analysis - identifying stereocenters in molecules

3. Stereoisomer enumeration - generating all possible isomers

4. 3D generation - Auto3D preserves input stereochemistry

5. Energy comparison - enantiomers identical, diastereomers different

Key takeaway: Always specify stereochemistry in drug discovery workflows. Undefined stereochemistry leads to ambiguous results.