Day 2 — What the Pipeline Gets Wrong
Day 1 ended with a working pipeline. A PDB ID and a SMILES string go in; a fully solvated, force-field-parameterized OpenMM system comes out. Every test passes. The simulation runs without crashing.
Day 2 is about asking a harder question: how correct is it? A pipeline can be internally consistent — no crashes, no missing atoms, valid force-field parameters — and still be modeling the wrong physics. This post is a detailed account of where version one falls short, using carbonic anhydrase II as the test case. CA2 is chosen deliberately: it is well-studied enough that the literature spells out exactly what a rigorous preparation looks like, which makes it a good mirror to hold up against what the code actually does.
The standard is not what you might expect
PDB files are not standardized input. The same protein — the same amino acid sequence, the same active site — can appear in hundreds of entries in the RCSB, each reflecting a different crystallographer's interpretation of electron density, a different ligand, a different pH in the crystallization buffer, a different resolution, a different organism. For CA2 alone, there are over 600 deposited structures.
This means that the "default preparation" — download, strip waters, add hydrogens, solvate — is actually a set of implicit choices masquerading as neutral operations. Which waters to keep? Which protonation states to assign? Which crystallographic artifact is a real feature? None of that has a universal answer. It depends on the protein, the ligand, and the scientific question being asked.
The current pipeline makes these choices silently. The goal of this post is to surface them.
The active site — a quick orientation
Before listing the limitations, it helps to know what CA2's active site looks like. The catalytic center is a Zn²⁺ ion held in place by three histidine residues (His94, His96, His119) in a distorted tetrahedral geometry. The fourth coordination site is occupied by either a water molecule or a hydroxide ion — the distinction is pH-dependent (pKa ≈ 6.8). His64 shuttles protons out of the active site during catalysis and adopts two distinct rotamers depending on pH.
|
His96 — Zn²⁺ — His119
|
H₂O / OH⁻ (pKa ≈ 6.8)
↕
His64 (proton shuttle)
When a sulfonamide inhibitor like acetazolamide binds, it deprotonates, displaces the Zn-bound water, and coordinates directly to Zn²⁺ through its nitrogen. The binding event is inseparable from the protonation chemistry. Any pipeline that ignores that chemistry is modeling a different system from the one that was assayed.
Seven things the pipeline gets wrong
Amber14 treats Zn²⁺ as a bare divalent ion. It has no bonded terms, no angle terms, no awareness of the His94/His96/His119 tetrahedral coordination shell. During minimization and equilibration, Zn will drift — there is nothing to hold it in place. Published CA2 workflows derive custom Zn partial charges from ab initio electrostatic potential fitting and add explicit harmonic bond and angle terms for each Zn–His interaction. The current pipeline has none of that.
At pH 7.4, the fourth Zn ligand is genuinely ambiguous: it could be H₂O or OH⁻.
receptor.py calls addHydrogens(pH=7.4) via PDBFixer,
which assigns protonation globally using a local hydrogen-bonding heuristic.
It cannot distinguish Zn–OH⁻ from Zn–H₂O. Getting this wrong changes the
electrostatics of the entire binding site and determines whether the sulfonamide
displaces a neutral or anionic ligand — a physically meaningful difference.
receptor.py passes keepWater=False, deleting every
crystallographic water. For most proteins this is a defensible default.
For CA2 it is documented to degrade results: Mukherjee et al. (2010) benchmark
CA2 crossdocking with and without active-site waters and show a clear performance
drop when they are absent. The waters within ~10 Å of Zn — including the
Thr199-mediated network — are not bulk solvent. They are part of the binding
pharmacophore.
The SMILES CC(=O)Nc1nnc(s1)S(N)(=O)=O represents the neutral
sulfonamide. The bioactive form is the deprotonated sulfonamidate anion
(–SO₂NH⁻), which coordinates directly to Zn²⁺ with sub-nanomolar affinity.
AM1-BCC charges computed on the wrong protonation state produce incorrect
electrostatics throughout the active site. Any scoring function trained on
these simulations will inherit this systematic error.
His64 is CA2's proton shuttle. Its rotamer preference — "in" (toward Zn) below pH ≈ 7, "out" above — is pH-dependent and affects access to the active site. PDBFixer will assign whatever protonation state is suggested by the local crystal environment, without any knowledge of the intended simulation pH or the rotameric context. For binding-affinity calculations near physiological pH, this is a meaningful source of variance.
After solvate.py injects the ligand particles and builds the OpenMM
System, there is no check that the active-site geometry is reasonable.
Zn–His distances, the coordination angle, and the presence of the fourth ligand
are not inspected. Errors here would propagate silently through minimization
and equilibration, producing plausible-looking trajectories that are modeling
a broken active site.
The pipeline hard-codes one PDB ID per protein. CA2 has over 600 deposited structures — different resolutions, crystal forms, pH conditions, alternate conformers, different organisms. The best structure for a given ligand is not always the one with the canonical PDB ID. Choosing the right starting structure requires reasoning about the specific assay conditions under which the affinity data was collected — the kind of context-sensitive judgment that is difficult to encode as a fixed rule.
Why this matters for the project
The goal of this pipeline is not just to run stable simulations — it is to extract observables that correlate with measured binding affinities. That requires modeling the correct binding state. For CA2 and sulfonamide inhibitors, the correct binding state involves deprotonated ligand, defined Zn coordination, a specific water network, and a resolved His64 rotamer. Version one gets none of those right.
This is not a CA2-specific problem. Every metalloenzyme in the target set (there are others beyond CA2) will have analogous metal-coordination issues. Every active site with titratable residues will have analogous protonation ambiguities. The "default preparation" is a floor, not a solution.
What gets built next
Closing thought
Version one is a foundation, not a finished pipeline. The machinery works: PDB in, solvated system out, simulation runs. But there is a gap between "the code runs" and "the code models the right physics," and Day 2 is about being honest about how wide that gap is.
The encouraging thing is that most of these limitations are not obstacles to building an automated system — they are arguments for it. A human expert reviewing a CA2 preparation would check the Zn coordination, inspect His64, and verify the ligand protonation state before running a single nanosecond of simulation. That expert reasoning is what the agent layer needs to encode. The hard part is not the MD — it is knowing what questions to ask before you run it.
A detour worth noting
Something became clear while working through these seven problems: the more interesting challenge is not establishing a general workflow from protein + ligand → simulate, but learning how to prepare specific systems correctly. Each target has its own physics — its own metal coordination, its own protonation ambiguities, its own crystallographic waters that matter. A generic pipeline will always be wrong about something important.
That realization points toward a different next step: a paper hunt paired with a batch of simulations drawn from BindingDB, using published protocols as ground truth to audit where the pipeline diverges and why. It is a more grounded way to make progress than continuing to generalize prematurely.
There is something satisfying about this kind of drift — starting with "automate the boring parts" and ending up deep in Zn²⁺ coordination chemistry and protonation state thermodynamics. Trying to make one piece of the pipeline more robust has a way of revealing exactly how much the other pieces were glossing over.
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)
Simulation trajectory
The 1CA2 + acetazolamide complex was equilibrated with an NVT ensemble (1 checkpoint frame) followed by NPT production (10 frames). The solvated box contains ~30 k atoms; only the protein and ligand are shown here. Drag to rotate, scroll to zoom. The animation cycles through NVT → NPT frames.
The ligand here is the anionic sulfonamidate form of acetazolamide
(CC(=O)Nc1nnc(s1)S([NH-])(=O)=O) — the species that actually
coordinates to the catalytic Zn²⁺. Using the neutral SMILES produces the
wrong protonation state and a poor Zn–N distance.