Heme Parameterization

Heme Parameterization

CYP450 catalysis cycles through several distinct electronic states of the heme. Each one has a different Fe oxidation state, axial coordination, formal charge, and partial-charge distribution. CYPForge ships parameters for three of them and the geometric machinery to place the right one in your pocket.

The three supported states

State	What it is	Fe oxidation	Axial ligand (distal)	Net charge	Use case
IC6	Resting Fe(III), low-spin or high-spin	+3	water / open / vacant	overall heme +1 with porphyrin / propionate	Substrate-bound but uncatalyzed; most X-ray structures show this state
DIOXY	Oxy-ferrous Fe(II)–O₂ adduct (Compound 0 precursor)	+2 with bound O₂	end-on dioxygen	similar carrier framework, modified Fe / O block	Modeling the O₂-bound resting state before electron transfer
CPDI	Compound I: oxo-ferryl Fe(IV)=O with porphyrin π-cation radical	formally +4 (or +3 with porphyrin radical)	Fe=O (oxo)	shifted charges on Fe, oxo-O, porphyrin N	Modeling the active oxidant; rate-limiting H-abstraction step

The state is declared at init time with --heme-state IC6|DIOXY|CPDI. Picking the wrong state is the single highest-impact undetected error in CYP450 MD setup — the Amber topology will build, minimization will run, and only careful inspection of Fe–N distances or trajectory dynamics will reveal that the partial charges and axial coordination are wrong.

Why these three?

The Shahrokh et al. 2012 paper derives consistent AMBER-compatible parameters for the full P450 catalytic cycle (resting Fe(III), reduced Fe(II), oxy-ferrous, Compound 0, Compound I, Compound II, hydroxide/hydroxo, etc.). CYPForge bundles the three states that cover the majority of CYP MD use cases:

• IC6 — what you have when you read a substrate-bound crystal structure off the PDB.
• DIOXY — what you have when you want to model the oxygen-binding step.
• CPDI — what you have when you want to model the active oxidant that abstracts hydrogen from the substrate.

The remaining Shahrokh states (Compound 0, Compound II, Fe(II) reduced, ferric hydroxide) are not in the default library; if you need them you can fit your own and load via the CUSTOM state hook in cypforge module heme prepare, but they do not propagate through the default ten-stage workflow.

What ships in each state directory

src/cypforge/data/heme_params/
├── IC6/
│   ├── HEM.mol2          # ferric heme, atom names + GAFF-style types + Shahrokh charges
│   ├── CYP.mol2          # axial cysteinate (CYM-equivalent), state-specific charges
│   └── IC6.frcmod        # bonded + nonbonded parameters for this state's Fe block
├── DIOXY/                # same triad for Fe-OO
│   ├── HEM.mol2
│   ├── CYP.mol2
│   └── DIOXY.frcmod
├── CPDI/                 # same triad for Fe=O Compound I
│   ├── HEM.mol2
│   ├── CYP.mol2
│   └── CPDI.frcmod
└── PROVENANCE.json       # which file came from which paper / figure / table

The .mol2 files carry the atom names, types, and partial charges; the .frcmod files carry the bonded force constants and any van der Waals overrides. Together they fully describe the heme block for one electronic state. CYPForge never mixes files across states — that's the state consistency property in Theorem 11.1.

Heme placement in the experimental pocket

Once the state is chosen, the template heme (in its own reference frame) has to be placed into the source structure's pocket. This is where most hand-built CYP450 systems silently fail. CYPForge does this with a six-layer construction; see Mathematical Foundations for the proofs. In one paragraph:

1. Plane fit — least-squares macrocycle plane from the 4 pyrrole nitrogens + 20 macrocycle carbons (SVD).
2. Normal orientation — the smallest eigenvector is unsigned; sign it so $+z$ is the distal face (SG below, O1 / propionate convention above).
3. In-plane phase — fix the rotational degree of freedom using the A/D propionate carbons.
4. Rigid transformation $\Phi$ — composition of template-frame and source-frame coordinate maps; provably $\Phi \in SE(3)$.
5. Atom matching — cylindrical-coordinate sort to pair source atoms with template atoms (diagnostic; final coordinates come from $\Phi$, not from per-atom matching).
6. Propionate-side QC — signed-distance product gate that detects residual proximal/distal flips before parameterization.

The output of Core 1 is a standardized PDB where the heme has the right coordinates, the right atom names, and the right state-specific block; the axial cysteine has been renamed CYM; and the Fe–S bond is recorded with explicit residue references for tleap.

Axial cysteine treatment

The proximal cysteine in CYP450s is deprotonated (thiolate, not thiol) — the SG ligates Fe directly. In Amber this requires:

• the residue to be named CYM (not CYS), so tleap selects the CYM library entry with no HG hydrogen,
• the explicit Fe–S bond statement,
• the state-specific charges from the corresponding CYP.mol2.

CYPForge handles all three. The axial_cys_resid (or the validated axial_cys_resname=CYM/CYP) is declared at init; src/cypforge/cys/axial_identification.py verifies the Fe–S distance is in the chemically reasonable 2.0–3.0 Å range (the bundled benchmark structures land at 2.3–2.6 Å); src/cypforge/cys/proximal_rewrite.py performs the rename. The Fe–S bond is then emitted to the tleap input as

bond mol.<CYM-index>.SG  mol.<HEM-index>.FE

where the indices are the LEaP residue numbers from the injective map $\rho$ described in Mathematical Foundations §3.

Why hand-stitching this is risky

Every individual step above is doable by hand:

• you can identify the right state from the literature,
• you can find the Shahrokh parameters and put them on disk,
• you can use PyMOL to align the heme, rename the Cys, and write the bond statement.

The problem is that doing all of them correctly and in a way you can audit six months later is much harder. Common silent failure modes:

What you forgot	What happens
Renaming CYS to CYM	tleap builds it as protonated cysteine, no Fe–S bond, Fe drifts off
Adding the bond statement	Same as above
Using Shahrokh IC6 charges with a Compound I geometry	Partial charges don't match the actual electronic state
Hand-aligning the porphyrin and rotating it 180°	Proximal/distal flip; substrate ends up where O₂ should be
Using two different heme states for two subunits	Asymmetry undetected until trajectory analysis

CYPForge does not eliminate the chemistry — you still have to decide which state matches your experimental context. But once that decision is made, the rest is mechanical and audited.

Validation set

The placement and QC pipeline has been exercised against six CYP450 structures spanning all three states:

PDB	Family / variant	Heme state	Notes
1R9O	CYP2C9	IC6	warfarin-bound
1Z10	CYP2A6	IC6	coumarin-bound (bundled in benchmark/)
2FDW	CYP24A1	IC6
2J0D	CYP3A4	IC6
4EJJ	CYP46A1	CPDI	bundled; Compound I model
4GQS	CYP4	IC6

Every structure passes the propionate-side QC gate with consistent sign pattern. This is not a calibration set in the machine-learning sense — it is a sanity-check ledger that confirms the geometry construction does what the math says it does on real X-ray structures.

Citation

When you publish results that use the bundled heme parameters, cite:

Shahrokh K, Orendt A, Yost GS, Cheatham TE III. Quantum mechanically derived AMBER-compatible heme parameters for various states of the cytochrome P450 catalytic cycle. J. Comput. Chem. 2012, 33(2): 119–133. doi:10.1002/jcc.21922 · PMID 21997754 · PMCID PMC3242737.

Per-file provenance: src/cypforge/data/heme_params/PROVENANCE.json.

Next: Ligand Parameterization — the SDF → GAFF2 → RESP pipeline.