Day 1 — System Assembly

Drug discovery pipelines are notoriously manual. Docking, parameterization, and simulation setup each require domain expertise and a different toolchain. The goal here is to collapse all of that into a single Python function — one that takes a PDB ID and a SMILES string and returns a fully solvated, force-field-parameterized system ready for OpenMM simulation. No GUI. No shell scripts. No manual file wrangling.

The bigger picture

The eventual target is a batch pipeline over ten well-characterised proteins from BindingDB: carbonic anhydrase II, trypsin, HIV-1 protease, thrombin, factor Xa, CDK2, BACE1, acetylcholinesterase, c-ABL, and PTP1B. For each target, a set of ~5 high-affinity and ~5 low-affinity ligands will be pulled from BindingDB (with experimental Kᵢ or IC₅₀ values) and run through the pipeline. The end goal is a scoring function that correlates simulation observables with measured binding affinity.

Day 1 was about getting the foundation right. The stack is deliberately modern and Python-native: PDBFixer for receptor prep, RDKit for conformer generation, OpenFF Toolkit + Interchange for ligand parameterization (AM1-BCC charges, Sage 2.1 force field), and OpenMM for everything simulation-side. No AmberTools binary calls, no Antechamber wrappers.

Test case for Day 1: Carbonic Anhydrase II (1CA2) with the clinical inhibitor acetazolamide — a zinc-binding sulfonamide with a well-characterised binding mode and sub-nanomolar affinity.

Simulation

Receptor — PDB 1CA2

Carbonic anhydrase II is a zinc metalloenzyme with a deep, cone-shaped active site. Its Zn²⁺ ion is the key pharmacophoric anchor — sulfonamide inhibitors coordinate directly to it. The crystal structure (1CA2, 1.54 Å resolution) was downloaded, stripped of co-crystallized ligands and waters, and repaired with PDBFixer: missing side-chain heavy atoms were modelled and histidine protonation states were assigned at pH 7.4 using a PROPKA-style heuristic.

Ligand — Acetazolamide

Acetazolamide is the textbook CA2 inhibitor — a 1,3,4-thiadiazole sulfonamide whose –SO₂NH₂ group deprotonates and coordinates to the active-site Zn²⁺ with sub-nanomolar affinity. Its SMILES (CC(=O)Nc1nnc(s1)S(N)(=O)=O) was fed into RDKit: hydrogens were added, a 3-D conformer was generated with the ETKDGv3 distance-geometry algorithm, and the geometry was refined with MMFF94. AM1-BCC partial charges were then computed by the OpenFF Toolkit and mapped onto the OpenFF Sage 2.1 force field via Interchange, producing a full set of bonded and non-bonded parameters.

What was built today

The hard part — parameterization

Parameterization is where automated pipelines usually break. The Amber force field (amber14-all.xml) handles the protein and water (TIP3P-FB) natively inside OpenMM, but it knows nothing about arbitrary small molecules. The solution used here is the OpenFF Toolkit: it assigns AM1-BCC charges (a semi-empirical electrostatics method that balances accuracy and speed), then the Sage 2.1 general small-molecule force field fills in the bonded terms. The two systems are merged by solvating the protein first and then injecting the ligand's particles and forces into the OpenMM System programmatically — sidesteppping the issue of Amber refusing to recognise an unknown residue.

Pipeline entry point

prepare_system(
    pdb_id        = "1CA2",
    ligand_smiles = "CC(=O)Nc1nnc(s1)S(N)(=O)=O",
    ligand_name   = "acetazolamide",
    output_dir    = "outputs/",
    padding       = 1.2,          # nm
    ph            = 7.4,
    ion_concentration = 0.15      # mol/L
)
# → ~31 000 atoms, box 6.9 nm, ready for OpenMM

What's next

Day 2 will extend the pipeline to a batch mode: pull active and inactive compounds for CA2 from BindingDB, run each through prepare_system(), and extract energy and structural descriptors to build a first-pass binding score. The longer arc is to repeat this across all ten target proteins and train a model that predicts affinity from simulation observables.

Stack

Python 3.11 · OpenMM 8.5.1 · OpenFF Toolkit 0.18.0 · OpenFF Interchange 0.5.2 · RDKit 2024.03.4 · PDBFixer · Amber14-all / TIP3P-FB (protein/water) · OpenFF Sage 2.1 (ligand)