Why CYPForge

Why CYPForge

Setting up a CYP450 protein–heme–ligand system for Amber MD by hand is one of those tasks where every individual step looks tractable but the coupling between steps quietly produces incorrect topologies. This page enumerates the pain points and shows where CYPForge intervenes.

The manual workflow

A typical hand-stitched setup of a CYP450 + heme + drug substrate complex involves, at minimum:

1. Determine the heme catalytic state. Resting Fe(III) (IC6), oxy-ferrous / Fe-OO (DIOXY), or Compound I oxo-ferryl (CPDI) — each has a different formal charge, spin multiplicity, and partial-charge distribution on Fe and the porphyrin nitrogens. Picking the wrong state is silent: tleap won't complain.
2. Locate compatible heme parameters. The Amber distribution does not ship a complete heme set covering all catalytic states. The community standard is Shahrokh et al. 2012, but you have to find them, verify them, and place them on disk.
3. Place the template into the experimental pocket. The bundled heme template lives in its own reference frame; the PDB pocket lives in another. Aligning them by eye in PyMOL is error-prone and produces an "upside-down" porphyrin in ~30% of attempts (proximal vs. distal face swap).
4. Rename the axial cysteine to CYM and declare the Fe–S bond. Forgetting either step gives tleap an unbonded Fe and a deprotonated cysteine floating in space — the topology builds, energy minimization runs, and only the trajectory tells you something is off.
5. Parameterize the ligand. SDF → GAFF2 typing → RESP charges. Now you have three files describing the same molecule (SDF, mol2, complex-PDB), each with its own atom ordering, atom names, and possibly differing hydrogen treatment. If atom mapping is done by row index, the charges land on the wrong atoms.
6. Set protonation states. Histidines (HID/HIE/HIP), glutamates/aspartates near the heme (GLH/ASH), and any non-standard residues each need a decision. Default tleap behavior is not always what you want for a CYP active site.
7. Solvate, ionize, equilibrate. Standard, but the cumulative cost of every prior decision shows up here as crashes during minimization or weird Fe–S distances after equilibration.

Every one of these steps is recoverable in isolation. The problem is that a failure introduced at step 3 will not surface until step 7, after hours of CPU and GPU time.

What goes wrong, concretely

Failure mode	Where it originates	Where it surfaces	How long to detect
Heme upside-down (proximal/distal flip)	Hand placement (step 3)	Fe–S distance > 4 Å after minimization	Hours
Wrong heme state used	Visual judgment (step 1)	Spin/charge mismatch — not detected	Never (silent)
Axial Cys still named CYS	Forgot to rename (step 4)	tleap builds an unbonded Fe	Visible in pdb4amber, ignored by tleap
Ligand RESP on wrong atoms	Row-index atom mapping (step 5)	Charges sum is correct, per-atom is wrong	Often never
tleap exit 0 ≠ correct topology	Implicit trust in zero exit code	Long-time-scale artifact	Days

What CYPForge changes

CYPForge is built around one observation: the failure modes above are all silent under a "command succeeded" criterion, but they all leave detectable evidence in the structure. So the shell is organized to generate that evidence at every step and gate on it.

The orchestration layer

Ten stages, fixed order, each one producing a JSON manifest with a hard PASS / WARN / FAIL:

environment_check
  → core1_prepare_heme_cym
  → core2_prepare_ligand_resp_gaff2
  → core3_finalize_protonation
  → core3_solvate_ionize
  → core3_render_pre_md
  → core3_run_pre_md
  → global_audit
  → equilibration_decision
  → production_readiness_check

FAIL stops the workflow. WARN pauses unless you set --auto-accept-warn at init. There is no implicit "carry on with the previous result" — every transition is an audited decision.

The math layer

• Heme placement is a rigid-body transformation $\Phi$ between the template frame and the source frame, constructed from an SVD plane fit + signed normal orientation + propionate in-plane phase. The proof in S3 shows $\Phi \in SE(3)$, so all template bond lengths, angles, and oriented dihedrals are preserved exactly. A proximate-side quality gate then catches any residual proximal/distal flip before the structure enters parameterization.
• Heme state is a discrete index $s \in \{IC6, DIOXY, CPDI\}$ that simultaneously fixes the template geometry, the HEM mol2, the proximal CYP/CYM mol2, and the frcmod. A "mix-and-match" inconsistency is structurally impossible if the state is declared.
• Ligand atom identity is recovered by strict graph isomorphism on the SDF / complex-PDB graph; if strict isomorphism returns no unique solution, a heavy-atom hypergraph fallback ranks candidates by element, neighborhood, pairwise distance, and Kabsch residual. Identity is never inferred from row order.
• RESP charge sum is gated at $|\sum q_i - Q_{\text{formal}}| \leq 10^{-4}$. A non-passing fit is rejected, not warned about.
• Residue addressing uses an injective map $\rho$ from (chain, resid, resname) to LEaP indices. The Fe–S bond mol.ρ(CYM).SG mol.ρ(HEM).FE statement is therefore globally unambiguous.

Details and proofs: Mathematical Foundations.

What this buys you, in time

Activity	Manual (typical)	CYPForge
Heme placement	10–60 min by hand in PyMOL	seconds (SVD + Φ)
Heme/state parameter lookup	Email collaborators, hunt papers	declared at init
Ligand parameterization (GAFF2 + RESP)	~1 hour of antechamber + Multiwfn juggling	one script, gated
Verifying topology correctness	run MD, look at trajectory	per-stage JSON manifests
Reproducing a prior setup 6 months later	ad-hoc notes	run_config.json + manifests

The wall-clock saving is real but not the main point. The main point is that what CYPForge produces is auditable — if a referee asks "how did you choose the heme state, the axial Cys ID, the ligand protonation?", the answer is in a JSON file you can attach as Supplementary Material.

What CYPForge does not replace

• It does not replace Amber/AmberTools. tleap, cpptraj, antechamber, parmchk2, and pmemd.cuda are all called as subprocesses; CYPForge just orchestrates them.
• It does not replace human judgment on protonation. Core 3 requires an explicit protonation_decision.json for non-default residue states — the shell will not guess HID vs HIE in the active site.
• It does not run production MD. The deepest stage is a 20 ns free NPT equilibration. Production sampling is yours.

Continue with Architecture for the workflow detail, or Mathematical Foundations for the proof structure.