Mathematical Foundations

Mathematical Foundations

This page condenses the formal contract behind CYPForge's heme placement and ligand identity layers. The full proof — Theorem 11.1, "conditional consistency of the assembly contract" — lives in the Supporting Information (S3 of the manuscript). What follows is the working-mathematician's tour: definitions, the three geometric layers, the discrete addressing layer, the ligand identity layer, and the consistency theorem that ties them together.

1. The assembly contract

CYPForge produces a tuple

$$ \mathcal{K} = (P^*, \Omega_s, \Phi, \rho, \phi, \mu_L, q_L, b, \mathcal{V}), $$

where:

Symbol	Meaning
$P^*$	the standardized prepared/assembly PDB
$\Omega_s = (M_H^s, M_C^s, F_H^s)$	HEM mol2, CYP/CYM mol2, frcmod for state $s$
$\Phi$	rigid transformation from the heme template frame to the source frame
$\rho$	injective map from (chain, resid, resname) to LEaP residue index
$\phi$	heavy-atom map from SDF to complex-PDB ligand
$\mu_L$	atom-name merge map from final mol2 to complex-PDB ligand
$q_L$	ligand RESP partial charges
$b$	the explicit Fe–S bond statement
$\mathcal{V}$	validation gates (file existence, tool execution, geometric/charge checks)

The contract is conditional: the claims below hold only when every check in $\mathcal{V}$ passes. Anything else triggers WARN or FAIL.

The state set is restricted to the three predefined cases:

$$ s \in \mathcal{S}_0 = { \mathrm{IC6},\ \mathrm{DIOXY},\ \mathrm{CPDI} }. $$

Custom states are supported by the low-level scripts but do not propagate through the default ten-stage path.

---

2. Heme placement — three geometric layers

Heme placement is decomposed into three layers, each with its own degeneracy condition. If any layer's input is degenerate, the workflow halts; it does not paper over a bad geometric situation.

Layer 1: SVD plane fit

Given the macrocycle core atoms $\mathcal{P} = {\mathbf{p}_1, \dots, \mathbf{p}_N}$, center them as $\mathbf{q}_i = \mathbf{p}_i - \bar{\mathbf{p}}$, stack into $A \in \mathbb{R}^{N \times 3}$, and minimize the out-of-plane residual:

$$ \hat{\mathbf{n}} = \arg\min_{|\mathbf{n}|=1} \mathbf{n}^{\mathsf{T}} A^{\mathsf{T}} A, \mathbf{n}. $$

Lemma. Let $M = A^{\mathsf{T}} A$ with eigenvalues $\lambda_1 \geq \lambda_2 \geq \lambda_3 \geq 0$. If $\lambda_2 > \lambda_3$ (simple smallest eigenvalue), the minimizer is the eigenvector of $\lambda_3$, unique up to sign.

If $\lambda_2 = \lambda_3$ the plane is degenerate — for example, a flattened porphyrin or pathological geometry. The implementation rejects this rather than picking an arbitrary eigenvector.

Layer 2: signed normal orientation

The SVD lemma only gives an unsigned pair ${\mathbf{n}, -\mathbf{n}}$. We orient $\hat{\mathbf{n}}$ so that $+z$ points to the distal face. For the source structure, the proximal sulfur must lie on the negative side:

$$ \mathbf{n} \cdot (\mathbf{x}_{SG} - \mathbf{x}_{FE}) < 0. $$

If the initial normal violates this, flip it. For the template:

• if O1 is present (DIOXY / CPDI states): $\mathbf{n} \cdot (\mathbf{x}{O1} - \mathbf{x}{FE}) > 0$;
• if O1 is absent (IC6): use the propionate midpoint $\mathbf{m}{prop} = \frac{1}{2}(\mathbf{x}{CGA} + \mathbf{x}{CGD})$ and require $\mathbf{n} \cdot (\mathbf{m}{prop} - \mathbf{x}_{FE}) < 0$.

Proposition (sign compatibility). Under these rules, the source frame and the template frame assign $+z$ to the same chemical face. The flip ambiguity is resolved consistently.

Layer 3: in-plane phase

With the normal fixed, fix the rotational degree of freedom in the porphyrin plane. Use the propionate midpoint vector $\mathbf{v} = \mathbf{m}{prop} - \mathbf{x}{FE}$, project onto the plane:

$$ \mathbf{e}_x' = \mathbf{v} - (\mathbf{v} \cdot \mathbf{e}_z)\mathbf{e}_z. $$

Normalize, take $\mathbf{e}_y = \mathbf{e}_z \times \mathbf{e}_x / |\mathbf{e}_z \times \mathbf{e}_x|$, then re-orthogonalize. If the projection $\mathbf{e}_x'$ is numerically near zero, fall back to a single propionate-carbon anchor; if even that is degenerate, halt.

Result. The basis $B = [\mathbf{e}_x\ \mathbf{e}_y\ \mathbf{e}_z]$ satisfies $B \in SO(3)$, i.e., it is a right-handed orthonormal frame. Local and global coordinate maps,

$$ \mathrm{loc}_F(\mathbf{x}) = B^{\mathsf{T}}(\mathbf{x} - \mathbf{o}), \qquad \mathrm{glb}_F(\mathbf{u}) = \mathbf{o} + B\mathbf{u}, $$

are mutual inverses (this is just $B^{\mathsf{T}} B = I$).

Layer 4: rigid transformation

Let $F_S = (\mathbf{o}_S, B_S)$ and $F_T = (\mathbf{o}_T, B_T)$ be the source and template frames. The placement map is

$$ \Phi(\mathbf{x}) = \mathrm{glb}_{F_S}(\mathrm{loc}_{F_T}(\mathbf{x})) = \mathbf{o}_S + R(\mathbf{x} - \mathbf{o}_T), \qquad R = B_S B_T^{\mathsf{T}}. $$

Theorem. $\Phi$ is a rigid transformation: $R \in SO(3)$, and $|\Phi(\mathbf{a}) - \Phi(\mathbf{b})| = |\mathbf{a} - \mathbf{b}|$ for any two points.

Corollary. All template-internal bond lengths, bond angles, and oriented dihedrals are preserved. Bond lengths are pairwise distances; angles follow from three distances via the law of cosines; oriented dihedrals are preserved because $\det R = 1$ (a reflection, $\det R = -1$, would flip dihedral signs).

This corollary is the structural justification for not re-parameterizing the template after placement: the geometry was good before $\Phi$ and is exactly as good after.

Layer 6: propionate-side quality control

A subtle failure: even with the normal correctly oriented, an atom-matching error could mis-assign which side of the macrocycle is "proximal". Re-fit the plane using NA / NB / NC / ND nitrogens, define the signed distance $\delta(\mathbf{x}) = \mathbf{n} \cdot (\mathbf{x} - \mathbf{x}_{FE})$, and require, with tolerance $\epsilon = 0.05\ \mathrm{\AA}^2$:

$$ \delta(\mathbf{x}_{SG}) \cdot \delta(\mathbf{x}_{CGA}) &gt; \epsilon, \quad \delta(\mathbf{x}_{SG}) \cdot \delta(\mathbf{x}_{CGD}) &gt; \epsilon, \quad \delta(\mathbf{x}_{SG}) \cdot \delta(\mathbf{x}_{m}) &gt; \epsilon, $$

where $\mathbf{x}m = \frac{1}{2}(\mathbf{x}{CGA} + \mathbf{x}_{CGD})$. A positive product means the two atoms lie on the same side of the macrocycle plane. A flip changes the sign of at least one product. The third inequality is mathematically implied by the first two but is recorded in the audit manifest for redundancy.

On the bundled validation set (1R9O, 1Z10, 2FDW, 2J0D, 4EJJ, 4GQS), every structure passes with consistent sign pattern.

---

3. Addressable residue mapping

This is the discrete layer. Let $R$ be the set of (chain, resid, resname) triples in the prepared PDB. For the subset $R' \subset R$ entering assembly, define an injective map

$$ \rho: R' \hookrightarrow I \subset \mathbb{N} $$

assigning sequential LEaP residue indices in output-PDB order. The proximal cysteine is renamed CYM, the heme is added as an extra residue, and the ligand (if ligand-aware) as another. The Fe–S bond statement is then

bond mol.ρ(CYM).SG  mol.ρ(HEM).FE

Proposition. Because $\rho$ is injective, distinct (chain, resid, resname) triples land at distinct LEaP indices, and the two endpoints of the Fe–S bond are uniquely referenced. This guarantees LEaP-input uniqueness, not successful tleap execution — those are separate.

---

4. Ligand identity

Two sources, one chemistry. The SDF provides the chemical graph $G_{\mathcal{T}} = (\mathcal{T}, E_{\mathcal{T}})$ with bond orders and formal charge. The complex PDB provides coordinates and atom names; the implied graph from coordinates is $G_{\mathcal{C}} = (\mathcal{C}, E_{\mathcal{C}})$.

Strict graph isomorphism

The atom-identity map $\phi: \mathcal{T} \to \mathcal{C}$ must be edge-preserving under element and degree constraints:

$$ { (\phi(a), \phi(b)) : (a,b) \in E_{\mathcal{T}} } = E_{\mathcal{C}}. $$

If this yields a unique solution, $\phi$ is determined.

Hypergraph fallback

If strict isomorphism returns no solution or multiple solutions (common for highly symmetric ligands), CYPForge 1.1.0 invokes a heavy-atom hypergraph fallback that ranks candidates by (i) element, (ii) graph-neighborhood signature, (iii) pairwise distance terms, and (iv) Kabsch residual. The fallback admits a result only if its decision is unique (mathematically determined) or equivalent_ok (remaining ambiguities are indistinguishable under retained checks or are explicitly recorded equivalent exchanges). Anything else is rejected.

equivalent_ok is recorded as an audit boundary — it is admitted as a working identity but is not promoted to a strict uniqueness proof.

Coordinate embedding

GAFF2 atom types are assigned by Antechamber from the SDF and written to a typed mol2. Before charges are injected, parmchk2 produces the frcmod. Then the complex-PDB heavy-atom coordinates (identified by $\phi$) are written into the typed mol2; hydrogen coordinates are rebuilt from SDF/GAFF2 template geometry. The result is required to have a heavy-atom RMSD against the complex PDB of at most $0.05\ \mathrm{\AA}$.

RESP charge gate

RESP charges are fit from an HF/6-31G* wavefunction (PySCF / GPU4PySCF) via Multiwfn. The total-charge consistency check is

$$ \left| \sum_{n \in \mathcal{N}_L} q_L(n) - Q_L \right| \leq 10^{-4}, $$

where $Q_L$ is the declared integer formal charge. Any deviation rejects the parameterization. This gate verifies charge-sum consistency only; it does not certify that every per-atom charge is chemically optimal.

Atom-name merge

During ligand-aware LEaP preparation, the final mol2 and the complex-PDB ligand residue are merged by atom name, not row index. If both sides have identical, duplicate-free name sets, the merge $\mu_L(n) = n$ is unique regardless of file order. If the name sets differ, the workflow records an audit warning and falls back to coordinate/element matching.

---

5. The consistency theorem

The seven properties below jointly define what "an internally consistent CYPForge contract" means. They are not independent — see the coupling remark.

Theorem (Conditional consistency of $\mathcal{K}$). If the input satisfies the boundary conditions of §1 and every check in $\mathcal{V}$ passes, then:

1. Geometric consistency. Heme geometry is the image of the template under the rigid $\Phi$, preserving all template-internal distances, angles, and oriented dihedrals.
2. State consistency. HEM mol2, CYP/CYM mol2, and frcmod all come from the same $\Omega_s$.
3. Face-orientation consistency. Both frames assign $+z$ to the distal face, and the propionate-side QC gate passes.
4. Residue-addressing consistency. Every residue has a unique LEaP address $\rho(r)$.
5. Fe–S bond consistency. The Fe–S bond endpoints are uniquely determined by the same $\rho$.
6. Ligand identity. If strict graph isomorphism is unique, $\phi$ is determined; if the fallback is invoked, $\phi$ is admitted only under unique or equivalent_ok. When mol2 and PDB atom-name sets match and are duplicate-free, the merge $\mu_L$ is unique.
7. Ligand charge-sum consistency. $|\sum q_L - Q_L| \leq 10^{-4}$.

Coupling remark. A single state index $s$ simultaneously fixes the orientation convention (§2), the template geometry transferred by $\Phi$, and the parameter set $\Omega_s$. So properties 1–3 share one state. A single injection $\rho$ fixes both residue addresses and the two Fe–S endpoints, so properties 4–5 share one address space. It is this sharing — not the seven statements taken separately — that makes the contract internally consistent rather than merely well-formed.

---

6. Scope and limits

The contract proves conditional internal consistency. It does not prove:

• conformational optimality of any Vina-derived pose,
• chemical validity of the heme template for states outside $\mathcal{S}_0$,
• successful tleap execution or downstream MD stability,
• that per-atom RESP charges are chemically optimal (only the sum is gated),
• that an equivalent_ok fallback decision is a strict uniqueness result.

A passing contract means: given that you accept the input as chemically valid, the structural, parametric, and addressing layers are mutually compatible. It is a guarantee about CYPForge's internal logic, not about the underlying chemistry. The chemistry is on you.

---

7. Where to find this in the code

Layer	Source
SVD plane fit	src/cypforge/heme/prepare.py
Normal orientation, in-plane phase	src/cypforge/heme/prepare.py
Rigid transformation $\Phi$	src/cypforge/heme/prepare.py
Atom matching	src/cypforge/heme/mapping.py
Propionate-side QC	src/cypforge/heme/prepare.py (signed-distance gate)
Residue addressing $\rho$	src/cypforge_core/heme_mapping_leapin.py
Strict graph isomorphism + fallback	src/cypforge_core/ligand_mapping_resolver.py, ligand_heavy_hypergraph_resolver.py
Atom-name merge $\mu_L$	src/cypforge_core/ligand_mapping_leapin.py
RESP fit	src/cypforge_core/ligand_gpu4pyscf_esp.py
Charge-sum gate	src/cypforge_core/ligand_gpu4pyscf_esp.py

Next: deep dives on Heme Parameterization and Ligand Parameterization.