Beyond the contact sum#

Every score in Assessing a modelled complex and Potentials is, at bottom, a sum of a residue-pair potential over a contact list. Five things a TCR:pMHC interface does are invisible to such a sum, and tcren measures each of them with its own instrument. This page is what those instruments are for; the signatures are in tcren package.

What a contact potential can and cannot express#

A contact energy is not purely an interaction: burying a residue against any partner costs something that depends on that residue alone. decompose() separates the two exactly, and only the pair part J is beyond what a per-position model can already write down.

from tcren.potential import mj, mj1996, mj_partition_energy

d = mj1996().decompose()          # e(a,b) = mean + H(a) + H(b) + J(a,b), J double-centred
d.h("F"), d.j("F", "W")           # one-body term; the genuinely pairwise remainder
d.energy("F", "W")                # reassembles the original value

f = mj1996().hydrophobicity_fit()  # C0 + C1(q_a + q_b) + C2 q_a q_b
f.r2, f.eigenvalue_share           # 0.98 on MJ1996, 0.84 on the bundled mj
mj_partition_energy()["F"]         # 4.37 -- MJ's own one-body scale (larger = more hydrophobic)

Where a potential has that shape the interaction term is only C2 q_a q_b, so it cannot prefer one pair of side chains over another of equal hydrophobicity. Both calls refuse a directed potential: TCRen is TCR-to-peptide and must not be split this way.

Peptide conformational stability: what a contact model cannot see#

A contact potential scores whichever conformation it is handed. It cannot tell a peptide whose own side chains hold it in the TCR-facing conformation from one that merely happens to have been modelled there – both present the same contact list. tcren.mechanics.dynamics puts the backbone in motion: it samples peptide phi/psi by Metropolis Monte Carlo against DOPE and reports how far the peptide wanders, not a better pose.

from tcren import peptide_stability, stability_table

peptide_stability(structure).rmsf        # ensemble spread, A -- larger = floppier
stability_table([s1, s2])["delta_rmsf"]  # intra-peptide term ON vs OFF, paired

The hypothesis it was built to test (Sewell, 2026-08): intra-peptide interactions stabilise the productive bulge a TCR reads, so a poor binder could make many contacts and still fail to stabilise the productive peptide conformation – which would explain why an additive contact model describes some systems well and others badly.

Tested on the combinatorial-peptide-library set: about 160 best-binder and 160 worst-binder modelled complexes for each of seven clones, 2,102 structures. AUC is best-against-worst discrimination.

clone	contact energy	stability
ila1	0.348	0.862
868	0.537	0.677
sb27	0.570	0.934
mel8	0.690	0.876
4c6	0.955	0.519
1e6	0.973	0.707
mel5	0.974	0.859

Stability beats the contact energy in 4 of the 4 clones where the contact model fails, and 0 of the 3 where it works. Mean AUC over the failing clones goes 0.536 to 0.837; over the working ones the contact energy stays ahead, 0.967 against 0.695. Combining the two as a within-clone z-sum lifts the mean AUC from 0.721 to 0.826, improved in 5 of 7 – though with seven clones that paired test is underpowered (Wilcoxon p = 0.22).

The intra-peptide term is a switch, and flipping it does what the hypothesis says. Removing the peptide’s contacts with itself lets the best binders’ backbones wander further (delta rmsf = +0.021 A, s.e. 0.005, 4.4 sigma) and leaves the worst binders unchanged (+0.002 A, s.e. 0.007); best against worst p = 0.042. The same term sharpens the stability discrimination itself, by +0.024 AUC on average and in 5 of 7 clones (Wilcoxon p = 0.078).

So the mechanism is supported while the system originally guessed is not: 4c6 is one of the clones the contact model handles well here (0.955), and the ones it fails on are ila1, 868, sb27 and mel8. Caveats worth carrying: these are modelled structures, the Monte Carlo is knowledge-based rather than molecular dynamics (no solvent, no force field, no time), delta rmsf is a mechanistic signal and not a useful classifier on its own (AUC 0.526), and every clone-level test has n = 7.

Side-chain repack: what a local minimiser cannot do#

tcren.energetics.rotamers.repack() (native _relax.repack) places every side chain in the chi rotamer the DOPE potential prefers. The rigid-body refiner moves the peptide and leaves every chi where it found it, so a full-atom model whose side chains a predictor placed keeps them – which is most of why a pairwise contact energy stops discriminating on generated poses.

from tcren import repack

fixed, report = repack(structure)   # report: n_conformers, energy, p_best per residue

Like for like – same wrong-rotamer input (chi1 rotated 120 degrees), same 33 to 42 side-chain atoms, same crystal reference:

	peptide side-chain RMSD (A)	time
input (wrong chi1)	4.131	–
`repack`	2.364	6 ms
OpenMM (anchor-restrained minimisation)	4.133	3,103 ms

OpenMM leaves them where they are, and that is not a defect in OpenMM: a local minimiser cannot cross the torsional barrier between two rotamer basins, so relaxing clashes and re-sampling rotamers are different operations and only a discrete packer does the second. Over eight structures the packer recovers side-chain RMSD 3.93 to 1.66 A, 8 of 8 improved, median 6 ms.

It rotates the side chains a model has. It cannot rebuild ones tcren.refine.substitute.substitute_peptide() stripped; that is side-chain construction, and it is still open.

Footprint shape: what the contacts say before they are scored#

Every other scorer here sums over contacts. tcren.topology.footprint reads the same contact map as a shape – which of the six CDR loops touched what, and whether the resulting footprint is one connected patch. No potential, no reference structure, no fitted parameter, and no canonical orientation: every descriptor is invariant under rigid motion, so unaligned inputs are fine. Only chain typing and CDR markup are needed, which the command line does in one batched annotation pass.

Coverage is the composition over cells – the 6 CDR loops times {peptide, MHC}, optionally splitting the peptide into thirds – summarised by the normalised Shannon entropy and by the Hill numbers (Hill 1973, Jost 2006), where D2 is the effective number of engaged cells. Topology joins the contacted pMHC residues at a C-alpha threshold and builds the flag complex: fp_b0_* counts footprint patches and fp_b1_* its holes. Coverage and topology are only weakly related, which is why they belong in one channel, and why that channel is read as T, a directional score against the native crystals rather than a hand-written z-sum.

$ tcren features -s structures/ -i topology -o shape.tsv

from tcren.topology.footprint import footprint_batch, footprint_features
from tcren.reliability import t_score

row = footprint_features(structure)   # one dict, the shape descriptors at the default two radii
row["D2_pep24"], row["fp_b0_r7"], row["L_canon"]

table = footprint_batch("structures/")  # polars frame, one row per structure
T = t_score(table)                      # the shape score, fit-free, one row is enough

The cyclomatic number of the bipartite contact graph (E - V + C) is deliberately not offered: with of order thirty contacts among of order thirty residues it is dominated by E and simply tracks interface size. The patch count is scale-free instead.

Surface topology: what a TCR meets before it binds#

A contact potential scores an interface that already exists. tcren.topology.surface describes the pMHC beforehand: the peptide sits in a groove between two helices, and a TCR coming down meets one surface, so the descriptor is a height field h(x, y) over that groove with hydropathy and charge painted on. The method follows SURFMAP (surface shell, per-cell feature, 8-neighbour smoothing, Manhattan map distance, hierarchical tree) and Protein Surface Topography (centre the chart on the functional site). A flat raster rather than SURFMAP’s equal-area spherical chart, because the TCR-facing surface is an open, near-planar patch that a plane does not distort.

from tcren import surface_map, surface_stats, surface_distance, surface_tree

smap = surface_map(structure)            # channels: h, phobic, charge; source: peptide/helix/floor
surface_stats(smap)["frac_above_ridge"]  # how much peptide surface clears the MHC helix crests
ids, d = surface_distance([m1, m2, m3])  # pairwise map distance -> epitopes cluster

Two things worth knowing, because both were defects first.

The frame is refit from every structure – z from the groove-floor plane normal, y from the peptide, origin on the peptide centroid. The floor’s own principal axis is not the groove axis (its beta-strands run across the groove), which put the two helices diagonally across the map. Because the frame is intrinsic, maps compare without prealigning the inputs, which is SURFMAP’s standing caveat.

Heights come from ray casting in the groove frame, not from Shrake-Rupley surface points. Sphere sampling is fixed in global axes, so the same structure rotated gave a different map (median cell moved 1.35 A, relief by 19 %). Ray casting is exactly equivariant and needs no probe test: the highest surface in a column is by definition the one nothing is above.

“Featureless” becomes a number. Over the 374 Canonical2026 complexes (230 distinct epitopes), the epitopes the literature names as featureless and as bulged separate completely.

epitope	source	rank by `frac_above_ridge`	`frac_above_ridge`	`relief` (A)
LPEPLPQGQLTAY	EBV BZLF1 13-mer, HLA-B35 – bulged*	2 / 230	0.749	2.81
HPVGEADYFEY	HCMV pp65 11-mer, HLA-B35:08 – bulged*	5 / 230	0.562	3.59
EPLPQGQLTAY	EBV BZLF1 11-mer, HLA-B35 – bulged*	8 / 230	0.416	2.54
LLFGYPVYV	HTLV-1 Tax, HLA-A*02:01 – prominent P5-Tyr	46 / 230	0.145	2.15
GILGFVFTL	influenza M1, HLA-A02:01 – featureless*	139 / 230	0.000	1.16
TAFTIPSI	HIV RT 8-mer, HLA-B51:01 – featureless*	205 / 230	0.000	0.95

Five of the eight most-protruding epitopes are literature-named bulged HLA-B*35 epitopes; both named featureless ones have no peptide surface clearing the helix crest at all. Structure-level AUC is 1.000 on relief, peak_to_valley and frac_above_ridge (p <= 0.001, 9 featureless against 5 bulged structures) – though with two distinct epitopes per group that is a 2-against-2 comparison, so the properly powered evidence is the trend over all 279 class-I structures: frac_above_ridge rises from 0.054 (8-mers) to 0.569 (13-mers), Spearman on relief +0.414, p = 5.5e-13.

notebooks/surface_topology.py draws the elevation, charge and hydropathy maps and reproduces this comparison.

Ring stacking: the geometry an identity cannot carry#

A contact potential scores a pair by identity, so two rings face to face at 3.5 A score exactly like the same two residues brushing past edge on. tcren.stacking.ring_stacking() measures the difference and returns no energy:

from tcren import ring_stacking

ring_stacking(structure, cutoff=7.5)  # centroid_distance, interplanar_angle, vertical, lateral

interplanar_angle near 0 is face to face, near 90 edge to face; a parallel-displaced stack shows a small vertical with a few angstrom of lateral. Proline is included – its pyrrolidine ring packs face on against aromatics through CH-pi contacts.