Historical contact-map Potts model

These fitted models and their empirical outputs are retained for historical and methodological reproduction. They are excluded from final fit-free scoring. The local benchmark contact GLM is distinct from this coupled model.

Every scoring path elsewhere in tcren reads a contact map that a structure has. This one models the map itself: which residue pairs could have contacted, which did, and how the answers at neighbouring cells depend on each other.

A site \(a = (i, j)\) is an available pair — a receptor residue \(i\) and a partner residue \(j\) whose C_α atoms lie within radius (15 Å by default) — and \(\sigma_a = 1\) iff a heavy-atom contact formed within cutoff (5 Å, the TCRen contact definition, unchanged). A whole contact map is the configuration \(\sigma\), and the model is a Boltzmann distribution over it:

\[E(\sigma) = -\sum_a \eta_a \sigma_a - \tfrac{1}{2}\sum_{a,b} A_{ab}\,\sigma_a \sigma_b, \qquad P(\sigma) = \frac{e^{-E(\sigma)}}{Z}, \qquad Z = \sum_{\sigma \in \{0,1\}^{n}} e^{-E(\sigma)}\]

Here \(n\) is the number of available sites, \(a,b\) index sites, \(\eta_a\) is the contact log-odds field, and \(A_{ab}\) is a symmetric zero-diagonal coupling. \(E\) is a dimensionless statistical energy; \(\log Z\) and entropy use natural-log units (nats). These are not measured binding free energies.

Site tables retain physical chain.id.rec and chain.id.par plus heavy-atom measurement probe (Å). Partner offsets and lengths belong to a chain and region together, so equal class-II groove-region names cannot merge chains. Orientation lookups index the same compact Cα rows used to compute their matrix. Unavailable orientations remain NaN; a model fitted with orientation rejects missing geometry explicitly instead of assigning it a bin. Legacy peptide-only tables remain readable. Legacy MHC tables without partner chain identity must be rebuilt from their structures. Frozen coefficient arrays are unchanged; class-II predictions require fresh assessment with corrected sites.

The one-body term is additive over categorical blocks,

\[\eta_a = \alpha + h^{\mathrm{rec}}(a_i) + h^{\mathrm{par}}(b_j) + J(a_i, b_j) + g_{\mathrm{dist}} + g_{\mathrm{region}} + g_{\mathrm{role}} + g_{\mathrm{class}}\]

so the fields carry single-residue propensity — “the backbone put you in reach; did your side chain engage?” — J carries the pair chemistry, and g the geometry and annotation the fields must be adjusted for. Every block is in the zero-sum (Ising) gauge, which makes J directly comparable with a double-centred Potential.

Why a reference state of available pairs

A TCRen potential is a Boltzmann inversion conditioned on a contact existing, so a residue that could have reached the peptide and declined contributes nothing at all. Here that non-event is the observable, which is what lets the fields separate reach from chemistry.

Why the couplings

Contacts on a real interface are not independent. On the 362 αβ Native2026 crystals, at fixed C_α distance, the chance of a contact runs 0.155 → 0.401 when the same receptor residue’s neighbouring peptide residue contacts (odds ratio 5.37 over 48,921 pairs). Three coupling families capture it, all defined on sequence offsets so no extra coordinates are needed:

K(di, dj)

within-loop, \(|di| \le 2\), \(|dj| \le 2\), both sites in one receptor loop and one partner chain — 12 classes after the \(K(d) = K(-d)\) symmetry;

L(|dj|, same chain?)

cross-loop, the two sites in different hypervariable loops — 6 classes;

M

cross-class, the same receptor residue against both the peptide and the MHC groove — 1 class, joint models only;

X(same receptor residue, |dj| > 2)

long-offset exclusion, the complement of K’s reach on one receptor residue — 1 class, far=True only. New in 3.8.0.

Everything else is asserted uncoupled: the model’s statement that linkage falls to zero across loops and beyond two residues.

Note

X exists because that statement is not neutral at long offset. K(+0,+1) and K(+0,+2) are strongly positive — a receptor residue gripping peptide position j also grips j+1 and j+2 — and nothing in the model then says the run stops, though a side chain gripping one end of the groove is plainly not also at the other. Fitted on 367 crystals with J pinned to tcren2, the class comes out at K = −0.457, s.e. 0.049, z = −9.34, worth +81.8 nats of pseudo-log-likelihood for one parameter.

It is off by default and every shipped model carries far=False, so nothing already written changes. Where it bites is a predicted Cα geometry: on the crystal’s own coordinates the field is already sharp and the couplings have little left to move (Jensen–Shannon similarity 0.7133 → 0.7162 over 367 crystals), while under a geometry averaged over reference crystals it improves both placement and level on three of the four Native2026 receptors with held-out priors — js 0.4846 → 0.4957, 0.5226 → 0.5350, 0.4834 → 0.4924, with the expected contact count moving 36.80 → 34.11 against 32 observed, 33.81 → 31.40 against 33, and 29.97 → 27.82 against 25. On the fourth (3qdg) it is flat, and that receptor is already under-counted.

On the crystals every axial class comes out positive and every off-axis class negative — a made contact recruits its own sequence neighbours onto the same partner residue and suppresses the diagonal one. That sign pattern is not visible in the raw data; it appears only once the axial terms and the distance profile are held fixed.

Fitting needs no partition function

The conditional \(P(\sigma_a = 1 \mid \sigma_{-a})\) is logistic in \(\eta_a + \sum_k K_k n_k(a)\) with \(n_k(a)\) the count of contacting neighbours in class k, so the coupled fit is an ordinary weighted-binomial GLM with a handful of extra integer covariates, and it stays concave. That is Besag’s pseudolikelihood; consistency for this model class is Ravikumar, Wainwright & Lafferty (arXiv:1010.0311), and plmDCA (arXiv:1211.1281) is the same recipe on Potts sequence variables. The design is over-parametrised and identified by an \(\ell_2\) ridge, then projected to the zero-sum gauge — penalise then project, because an \(\ell_2\) penalty otherwise picks its own gauge.

PottsModel.n_parameters() counts stored fitted coefficients, including orientation and redundant or unobserved one-body levels. It is not the number of identifiable degrees of freedom. pseudo_loglik stores the unpenalized weighted value at the penalized optimum. IRLS validates its numerical domain and backs off Newton steps that decrease the penalized objective; failure to converge raises.

Scoring does, and it is exact at the reference

Z is estimated by annealed importance sampling (Neal, Stat. Comput. 11:125–139, 2001), annealing only the coupling term. At \(\beta = 0\) the model is the uncoupled one, whose partition function is exact and closed form, \(\log Z_0 = \sum_a \log(1 + e^{\eta_a})\) — so the reference is a verified model rather than an approximation, and the estimator is unbiased in Z. Transitions are block Gibbs on a greedy colouring of the actual coupling graph; same-colour sites are conditionally independent, which is asserted against the real edge lists rather than argued.

Always read ais_ess in the output. It is the effective sample size of the AIS weights out of --particles; close to --particles means the annealing schedule was long enough, and a small value means it was not.

Bound versus unbound, for the whole interface

A single site has two states and \(\eta_a\) is the free energy between them. The same contrast for the whole interface needs a macrostate, and the contact count \(N(\sigma) = \sum_a \sigma_a\) defines one. Because \(E(\varnothing) = 0\) exactly, the observed map’s log-odds against the empty configuration is just its negated energy, and the scores already emitted by score_sites() decompose:

\[\underbrace{-E(\sigma^{\mathrm{obs}})}_{\text{binding log-odds}} \;=\; \underbrace{\log Z}_{\text{capacity}} \;+\; \underbrace{\mathcal{L}(\sigma^{\mathrm{obs}})}_{\text{typicality}}\]

bound_unbound() gives three readings of that contrast, all from one Gibbs pass, because every tilt in \(N\) is an exponential family:

df_empty

\(\log[P(N \ge 1)/P(N = 0)] = \log(Z - 1)\), exact, from the AIS log Z.

df_threshold

\(\log[P(N \ge x)/P(N < x)]\) from the sampled histogram — Z cancels, so no AIS is needed, but it is only finite while x sits inside the sampled range.

mu_star

the chemical potential at which \(\langle N\rangle_\mu\) equals the observed count. Positive means the structure made more contacts than its fields and couplings warrant. nan outside the sampled support, where reweighting would be extrapolating.

The two are not competing estimates of one number. The unbound basin of a docked pose is astronomically improbable — the model is conditioned on an available set that already holds the receptor against the peptide — so no sampler reaches \(N = 0\), and only the log Z route gives it. count_profile() returns the pooled \(F(N) = -\log p(N)\) profile so the landscape can be inspected before a threshold is chosen.

Constraining a statistic of the whole configuration

Jaynes’ construction lifts a constraint from a single site to any statistic of the configuration: fix \(\langle O_k\rangle\) and its multiplier enters the Hamiltonian as \(-\lambda_k O_k(\sigma)\). Tkačik et al. (PLoS Comput. Biol. 10(1):e1003408, 2014) is this with O the total activity — the “K-pairwise” model — and it is why gibbs() takes an observer callback: a statistic of whole configurations can be accumulated during sampling without materialising every draw.

Because such a statistic depends on \(\sigma\) only through a low-dimensional summary, the tilted expectation is an importance-weighted average over draws taken at \(\lambda = 0\), so one sampling pass serves every step of a moment-matching fit. tilt_mean() and mu_star() implement that reweighting for the contact count. A linear tilt is exactly a constant added to every field, \(E - \mu N = -(\eta + \mu)\cdot\sigma\), which is what makes the reweighting identity checkable against direct simulation rather than assumed.

Shipped models

Key	Sites	Contacts	Provenance
potts_tcr_peptide	64,622	7,865	The default. TCR:peptide over the 362 αβ Native2026 crystals, redundancy-balanced on both axes. Reproduce with tcren potts fit -s data/Native2026 -o potts_tcr_peptide.json --balance both.
potts_tcr_peptide_sc	64,064	7,799	The same interface with the pair term pinned on tcren.potential.tcren2_em() rather than on TCRen2, over the 359 crystals that survive the withdrawn-entry rule; beta_matrix +0.7629. Reproduce with tcren potts selfcon (see Potentials) — it is the model that loop leaves behind, not a separate fit.
potts_tcr_mhc	239,093	15,451	TCR:MHC groove over the same crystals — twice as many contacts as the peptide interface. Reproduce with tcren potts fit -s data/Native2026 --partner mhc --balance both -o … (needs the allele reference: tcren build-mhc-ref).
potts_pep_mhc	195,674	23,492	The presentation arm: the groove’s grip on the peptide, over the same crystals. Here the groove is the receptor, so region.rec is a groove region and the model carries tcren.potts.MHC_RECEPTOR_REGIONS rather than the TCR loop set. Reproduce with tcren potts fit -s data/Native2026 --receptor mhc --partner peptide --balance both -o ….

Load any of them with tcren.potts.PottsModel.bundled(), or by name on the command line — -m/--model on potts score, contacts, map, ensemble and scan takes a bundled name as well as a JSON path:

tcren potts contacts -s complex.pdb.gz -m potts_tcr_peptide_sc -o contacts.tsv

The three arms are three separate Hamiltonians, not one field read three ways, and the fit says so: the peptide:MHC coupling block correlates with the shipped TCRen2 at r = +0.017 and with Miyazawa-Jernigan at r = +0.042. The groove’s chemistry is its own. Its per-position contact profile is the textbook one — on 1ao7, peaks of 10, 9 and 10 contacts at P1, P2 and P9 against 1 at the TCR-facing bulge.

Which potential belongs on which interface

Fixing J to one scale on a bundled potential (--coupling-matrix) gives every candidate an identical parameter count and an identical design, so their pseudo-log-likelihoods compare directly. Because the matrix is double-centred it contributes nothing to the one-body marginals — the fields carry those, refitted freely — so the comparison is about pair structure alone.

On the Native2026 crystals the ranking inverts between the two interfaces: TCRen2 beats MJ by 103.3 nats on TCR:peptide, and MJ beats TCRen2 by 35.5 nats on the TCR:MHC groove. TCRen2’s fitted scale falls 5.4-fold across that move (\(\beta\) = +1.131 → +0.209) while MJ’s barely changes (+0.803 → +0.974). This is the measurement behind tcren’s long-standing default of scoring Phi_tcr_mhc with Miyazawa–Jernigan and reserving TCRen for TCR:peptide.

Three limits of one free energy

The interface free energy is closed form, \(\Phi = \sum_a \operatorname{softplus}(\eta_a)\), and its derivative is the contact marginal, \(\partial\Phi/\partial\eta_a = p_a\). So \(\Phi\) is an interaction sum weighted by how often each pair is in contact, and the fixed contact map every other scoring path in tcren reads is the special case \(p_a \in \{0, 1\}\). Three limits follow, and they are not competing models — each is the one the interface’s own physics selects.

limit	weight on each pair	the interface it fits
hard contact	\(\mathbf{1}[d \le 5\,\text{Å}]\)	reproduces tcren.ddg.reference_delta() exactly, when J is pinned with pin_centred=False
smoothed	\(p_a \in (0, 1)\)	a plastic interface, where an ensemble average is a different number from one snapshot
saturated	\(p_a \to \{0, 1\}\)	a frozen interface, where the two above coincide and a fixed map is already exact

Which limit an interface is in is measurable, not a matter of taste. Over four deposited 100 ns trajectories, collapsed onto the peptide and restricted to the pairs inside the availability radius, 36 of the 38 engaged peptide positions (94.7%) reach a maximum groove-pair contact frequency above 0.98, against 14 of 38 (36.8%) on the receptor side — at the groove a peptide residue is in contact in essentially every frame or not at all, while at the receptor loops the same maximum runs the whole range from 0.00 to 1.00. Among pairs that touch at all the groove mean is 0.468 against 0.322, and 28.6% of groove pairs sit above 0.98 against 17.4% of receptor pairs.

That is why the groove takes a fixed map and Miyazawa–Jernigan: it is already in the saturated limit, so smoothing it can only add noise. The receptor loops are not, which is what the smoothed limit is for.

What the marginals settle, and what they do not

The interaction energy of a contact map is linear in the map, \(\Phi(\sigma) = \sum_a \sigma_a \varphi_a\), so by linearity of expectation

\[\langle\Phi\rangle \;=\; \sum_a \langle\sigma_a\rangle\, \varphi_a\]

for any joint distribution over \(\sigma\), coupled or not. Dependence between sites does not break that identity. What it does is fix its two conditions, and both are easy to violate:

1. \(\langle\sigma_a\rangle\) must be the marginal of the coupled model — p_model, by block Gibbs — and not p_independent, the factorised \(\sigma(\eta_a)\).

2. It must be read at the sequence being scored. ddg(weights=...) and cpl.predict_response(tcr_weights=...) take one weight vector computed on the template and apply it to every threaded candidate, giving \(\sum_a p_a(\text{template})\,\varphi_a(x)\). A candidate residue enters \(\eta_a\) twice, through \(h^{\mathrm{par}}\) and through \(J\), so the contact ensemble itself moves under substitution and the two are different numbers.

tcren.potts.threaded_ensemble() recomputes the marginals per candidate. It accepts peptide-partner sites only and validates each whole candidate against par.len, including unrepresented tail residues. Later candidates start from the previous retained matrix of chain configurations and use candidate_burn; candidate ordering remains part of the sampling protocol. tcren.potts.energy_ensemble() carries the read-outs no set of marginals can determine:

read-out	units	what it is
phi_var	(potential units)²	\(\operatorname{Var}\Phi\), which needs the two-point function \(\langle\sigma_a\sigma_b\rangle\). How much the energy depends on which of the maps the interface admits is the one realised.
phi_free	potential units	\(-\beta^{-1}\log\langle e^{-\beta\Phi}\rangle\), a soft minimum over configurations. \(\beta \to 0\) recovers \(\langle\Phi\rangle\); \(\beta \to \infty\) the best map the geometry admits.
map_entropy	nats	\(\langle E\rangle + \log Z\), how many contact maps the interface is choosing between.
hamming_obs	count	mean distance from the map the structure itself made.

On the 14 αβ complexes in tests/assets/pdb the spread is not a rounding effect: \(\operatorname{Var}\Phi\) runs 0.341 to 1.286 squared TCRen units, a 3.8-fold range, so how plastic an interface is varies as much between complexes as its mean energy does.

A configuration needs coarse-graining before it has a probability

tcren.potts.map_ensemble() returns the contact configurations themselves, ranked, each with its probability. The granularity is not a free choice, and it is measured rather than assumed. On 1ao7’s TCR:peptide interface — 174 available sites, 6,400 draws — the raw site-level map recurs never: 6,400 distinct configurations out of 6,400 draws, every one at \(\hat p = 1.6\times10^{-4}\), because the state space is \(2^{174}\). Coarse-grained to which peptide positions are engaged at all, the same draws collapse to 41 distinct patterns with the most likely at \(\hat p = 0.3219\) — and that modal pattern is the one the crystal itself made.

by	groups on 1ao7	distinct patterns / 6,400	reading
"pair"	174	6,400, top \(\hat p\) 0.0002	unresolvable; use tcren.potts.exact_log_z() on an interface small enough to enumerate
"loop"	39	3,308, top \(\hat p\) 0.0216	the granularity the molecular-dynamics validation works at; at the edge of resolvable
"position"	9	41, top \(\hat p\) 0.3219	which peptide positions the receptor engages — the default

Across those same 14 complexes the modal engagement pattern is the observed one in 6 of 14, and n_distinct_total runs 37 to 556.

When the contact indicator is not data

tcren.potts.fit_potts() treats \(\sigma\) as observed. In a crystal it mostly is — but a contact is called by thresholding the closest heavy-atom distance at 5 Å, and a pair sitting within half an ångström of that line is a coin toss recorded as a fact. tcren.potts.em treats those sites as latent and runs a Gibbs-imputation/pseudolikelihood hybrid around the existing fit. The full-posterior imputation and conditional-pseudolikelihood update do not share the ordinary likelihood-EM monotonicity theorem.

E-step. Sample \(\sigma_{\mathcal L}\) from \(P(\sigma_{\mathcal L} \mid \sigma_{-\mathcal L})\) by block Gibbs with every confident site clamped to what the structure did. The couplings are what make this worth doing: a borderline pair beside three made contacts is imputed differently from one standing alone.

M-step. The existing penalised pseudolikelihood on the completed data, draws completions stacked at weight \(1/\text{draws}\) so the effective sample size — and therefore what the ridge means — is unchanged. Each completion carries a #m suffix on its pdb.id, because every coupling edge in tcren.potts.kernel joins on that column and unsuffixed copies would link sites across draws that never coexisted.

Warning

The latent band has to reach both ways. available_pairs measures d_heavy out to a requested probe=cutoff, so every pair that missed carries a null distance. The current default probe is the availability radius. A band built on that table is one-sided: on the 17 bundled test structures it selects 44 sites, every one of them a contact, and an EM run on it could remove contacts and never add them. Build the table with available_pairs(..., probe=cutoff + margin) — which populates d_heavy past the cutoff without moving sigma — and the same band selects 115 sites, 44 contacts and 71 near misses. tcren.potts.latent_mask() checks the recorded probe against cutoff + margin; an all-contact band alone does not establish truncation. Legacy tables without probe provenance need regeneration for publication. margin=0 explicitly makes every site observed.

This adds a fit; it does not replace one. The shipped pseudolikelihood models are unchanged and are the point of comparison. tcren.potts.consistency() reports how far a refit moved per parameter block, scaled by the reference fit’s own cluster-robust standard errors where it has them (the coupling kernel) and by the block’s own standard deviation where it does not. A block with no spread — g_class on a single-partner model — is reported as degenerate rather than divided by something near zero.

Note

Monte-Carlo EM with a fixed draws does not converge to a point. It random-walks in a ball around the optimum whose radius is set by the Monte-Carlo noise, so max_shift has a floor that tol may sit below. Read max_shift_vs_init beside it: a run that is merely jittering has a flat one. Raise draws to lower the floor. Every sampler is seeded, so a run is exactly reproducible either way.

tcren potts fit -s <structures> -o potts_em.json --em --em-draws 5 --em-steps 10 \
    --em-trace-out em_trace.tsv

Retired geometry experiments

The geometry prior, coordinate-chain model and cross-product experiment are archived in the manuscript research archive (2026-10-audit). They are not shipped APIs. The contact model on supplied coordinates and its fitting recipes remain supported.

Scoring the agreement: placement and level

tcren.potts.map_agreement() reads any scored site table — from contact_probabilities() or tcren potts contacts — and reports two numbers per map, because no single one answers both questions:

ratio

The level, \(\sum_a p_a / \sum_a \sigma_a\). One is the contact count reproduced.

js

The placement, Jensen–Shannon similarity of the two maps read as distributions over sites, \(q_a = \sigma_a/\sum\sigma\) against \(r_a = p_a/\sum p\), in bits:

\[\mathrm{js} = 1 - \tfrac12 D_{\mathrm{KL}}(q\|m) - \tfrac12 D_{\mathrm{KL}}(r\|m), \qquad m = \tfrac12(q+r)\]

Bounded in \([0, 1]\), symmetric, and scale-free in the level — normalising divides the total out, which is why ratio sits beside it rather than being folded in.

tcren potts agreement -i predicted.tsv
tcren potts agreement -i contacts.tsv --column p_independent -o agreement.tsv

Read ``js`` against the same model on the structure’s own geometry, not against 1. \(\mathrm{js} = 1\) requires the model to put all its mass on exactly the contacting cells in equal proportion. A probabilistic map \(p=c\sigma\), \(0<c<1\), also reaches 1 after normalization; ratio records its different level. On the four receptors of the manuscript’s Figure 2c the crystal’s own C\(\alpha\) geometry reads 0.75–0.81 and the predicted geometry 0.42–0.58, so the first is the ceiling the second is measured against.

R² is deliberately not offered. On a binary outcome 1 - SS_res/SS_tot is the Brier skill score — the two agree to three decimals on every map checked here — so reporting both is reporting one number twice. It is also dominated by the many easy zeros and scales with prevalence, which makes two maps of different site counts incomparable under it; that is exactly the comparison the predicted and measured geometries require, since the predicted site set is roughly twice the size.

Reproducing a referenced score

centred_potential double-centres by default, which is right for ranking potentials against each other — with the one-body terms removed the comparison is about pair structure alone. It is wrong for reproducing a referenced contact-map score. The peptide-referenced energy is a difference of one-body sums, so double-centring re-injects a burial-scaled composition term \(n_i\,c(a)\): the position’s contact count times the potential’s partner-residue column mean. On TCRen2 that mean has s.d. 0.0668 and runs −0.212 to +0.027 over the 19 residues the potential observes, and \(n_i\) runs 1 to 54.

Pass pin_centred=False (or centre=False to centred_potential) and the coupling is the potential, so any linear read-out of the field reduces to the potential’s own score up to the fitted scale \(\beta\).

Command line

# fit a model (the alpha-beta TCR:pMHC HARD RULE applies, as in derive-potential)
tcren potts fit -s structures/ -o potts.json --balance both

# energy, log Z and the likelihood of each structure's observed contact map
tcren potts score -s structures/ -o scores.tsv          # bundled model by default

# per-residue-pair contact probability
tcren potts contacts -s complex.pdb -o contacts.tsv

# close those onto the grids an experiment measures
tcren potts map -s complex.pdb --by loop     -o map.tsv    # CDR loop x peptide position
tcren potts map -s complex.pdb --by position -o import.tsv # peptide residue importance

# historical log-Z0 substitution effects (nats; not calibrated binding energy)
tcren potts scan -s complex.pdb -o scan.tsv

score emits one row per structure: n_sites and n_contacts (the available pairs and how many of them engaged), energy and neg_energy (the Hamiltonian of the observed map and its negation, lower energy being more favourable), log_z and log_z0, log_lik and psi (the log-likelihood, and the same per available pair so it compares across interfaces of different size), pseudo_log_lik/psi_pseudo as the MCMC-free cross-check, and ais_ess.

contacts emits one row per site with three probabilities, whose differences are the couplings: p_independent from the one-body model alone, p_model the marginal of the full coupled model by block Gibbs (the one to use), and p_conditional = \(P(\sigma_a = 1 \mid \text{the observed rest})\).

map groups pairs onto a loop/position grid. Its p_any uses the independence approximation

\[P_{\mathrm{ind}}(N \ge 1) = 1 - \prod_j (1-p_j),\]

where \(p_j\) is the coupled marginal of pair \(j\). Distinct pairs need not be independent: coupled marginals alone do not determine the probability of their union. This approximation can be compared empirically with a trajectory contact-frequency map, but is not its exact coupled model expectation. No numerical definition changed in the audit.

p_expected is \(\sum_j p_j\), the expected number of contacts; linearity makes this valid without independence, and it need not lie in [0, 1]. Other columns are n_pairs, n_observed and observed. p_any lies in [0, 1]; the product is accumulated in log space. Neither this probability approximation nor the expected count is a binding free energy.

score, contacts and map all take --workers (default: every core). The per-structure numbers are functions of (seed, pdb.id) alone, so splitting the work changes nothing.

Historical substitution log-partition and linear response

map reads how engaged a position is expected to be before any residue identity is scored. peptide_free_energy() reads what happens when the identity changes. The partner residue enters \(\eta\) in exactly two places, so for a site \(s\) and a candidate residue \(a\)

\[\eta_s(a) = r_s + h^{\mathrm{par}}(a) + J(\mathrm{rec}_s, a)\]

with \(r_s\) — the intercept, the receptor field, the distance profile, region, role and class — independent of \(a\). Threading \(a\) through partner position \(i\) moves every site carrying that position, and the uncoupled log-partition changes with it:

\[\Phi^{\mathrm{Potts}}(x) = \log Z_0\big(\eta(x)\big) = \sum_s \log\!\big(1 + e^{\eta_s(x)}\big), \qquad \Delta\Phi_i(a) = \Phi^{\mathrm{Potts}}(x_{i \to a}) - \tfrac{1}{20}\sum_b \Phi^{\mathrm{Potts}}(x_{i \to b})\]

Higher is more favourable, and the reference is the equimolar one — the mean over the twenty residues at that position rather than the residue the structure carries — which is the null a positional-scanning library holds its other positions at. coupled=True takes the linear response about the observed sequence, \(\Delta \log Z \approx \sum_s p_s\,\Delta\eta_s\), since \(\partial \log Z/\partial \eta_s = \langle\sigma_s\rangle\): one Gibbs pass, then a dot product per cell. Because \(\log Z_0\) is a sum over independent sites the result is additive over positions, so one \(L \times 20\) table scores a single substitution and any whole partner sequence alike. Only aa.par changes — the backbone, the Cα distances, the receptor residues and the partner roles stay the structure’s own, the same fixed-backbone approximation every threading score in the package makes. Custom radii/cutoffs must also agree with the model’s fitted categorical blocks; changing the radius alone does not create missing field levels. A partner position carrying two different residues has no sequence to substitute into and is rejected rather than averaged.

The positive \(\log Z_0\) is dimensionless and measured in nats. The uncoupled statistical free energy is \(-\log Z_0\); neither is a calibrated binding free energy. With coupled=True, the preserved log_z0 column contains the frozen- marginal linear surrogate, not coupled \(\log Z\). Finite substitutions do not make the first-order approximation exact. Existing names and values are unchanged.

Only the full-field centred F_full quantity overlaps the uncoupled scan. The centred F_full_p quantity overlaps its fixed-native-marginal path when identical sites, parameters and marginal values are supplied. Benchmark composition-removed, pair-only, observed-support, normalized-contact and within-clone standardized combinations have no complete equivalent CLI route. Do not replace them with potts scan. Their equations, original results and losses remain historical evidence; no new fitted replacement is authorized. See the issue-32 disposition record in the benchmark research archive (16 October 2026).

The two terms it separates are different quantities, and which one a task needs is an empirical question. \(h^{\mathrm{par}}\) is composition — how much a residue engages an available partner at all, wherever it is put — while \(J\) is complementarity, the pair chemistry, and is the block \(\beta_\Phi\) ties to TCRen2 above.

\[\text{eta} \;=\; \underbrace{h^{\mathrm{rec}} + h^{\mathrm{par}}}_{\text{composition}} \;+\; \underbrace{J}_{\text{complementarity}} \;+\; \underbrace{g}_{\text{geometry}}\]

Closing the loop: the potential re-derived from the contacts the model predicts

Everything above reads a finished potential. tcren.derive_tcren() tallies the amino-acid pairs of the contacts a crystal shows, and tcren.potts.fit_potts() reads the resulting matrix back as the pair term in its field. Nothing carries the model’s own view of which contacts formed into the tally. tcren.potts.selfcon closes that loop.

\[\begin{split}\eta_a^{\mathrm{post}} &= \eta_a(\theta^{(t)}) + \operatorname{clip}\!\left(\frac{c - d_a}{m},\ \pm K\right) \\ p_a^{(t)} &= \langle \sigma_a \rangle_{P(\sigma \mid \mathrm{obs},\, \theta^{(t)})} \\ N^{(t)}(a,b) &= \sum_s w_s \sum_{a\,:\,\mathrm{pair}(a)=(a,b)} p_a^{(t)} \\ e^{(t+1)}(a,b) &= -\ln \frac{N^{(t)}(a,b)\,N_{\mathrm{tot}}} {N_{\mathrm{row}}(a)\,N_{\mathrm{col}}(b)}, \qquad \theta^{(t+1)} = \arg\max\ \mathrm{PL}\big(\theta \mid e^{(t+1)}\big)\end{split}\]

with \(c\) the 5 Å cutoff, \(d_a\) the site’s closest-heavy-atom distance, \(m\) the margin, \(K\) the logit clip and \(w_s\) the structure’s redundancy weight.

Weighting each conformation by its probability costs nothing. The tally over whole contact configurations, \(\sum_\delta P(\delta \mid \mathrm{obs})\, n_{ab}(\delta)\), equals \(\sum_a p_a\) by linearity of expectation — the same identity that governs tcren.potts.energy_ensemble(). So the configurations never have to be enumerated to get the counts; tcren.potts.map_ensemble() lists them as a read-out, not as a dependency.

The estimator does not change; only the counts it is fed. tcren.derive_tcren() already accepts a per-contact fractional weight, documented there as the way to let a pair that “will only sample an interacting geometry some of the time” cast a fractional vote. The log-odds, the contact-fraction reference state and the pseudocount are the shipped ones, untouched.

The observation enters as one additive field. The distance sigmoid \(\operatorname{sigmoid}((c-d_a)/m)\) defines a multiplicative observation field. A Bayesian posterior interpretation requires its logit to be the observation likelihood ratio, or the sigmoid to use an equal-odds reference prior; an arbitrary posterior probability given distance does not suffice. The field adds a constant to \(\eta\), so the E-step is the existing block-Gibbs sampler called with a modified vector — no clamping machinery, no second code path. One knob spans the family: margin → 0 pins the crystal map, margin = 0.5 Å releases the threshold band alone, and a large margin lets the crystal enter only through the fit.

The pair term is pinned, not free. J carries one coefficient — beta_matrix, the scale on the current matrix — rather than 400 cells. That is what closes the loop: a free J would absorb the pair structure and leave the re-derived matrix no route back into \(P(\sigma)\). It also makes beta_matrix reportable as a statistical scale relating TCRen natural-log units to contact log-odds. It is not a calibrated physical temperature. tcren.potts.fit_potts() accepts a Potential object for coupling_matrix, not only one of the five bundled names, and a model pinned that way stores the matrix itself so it round-trips.

Warning

The heavy-atom layer has to reach past the cutoff, exactly as for tcren.potts.em. Build the table with available_pairs(..., probe=cutoff + 3*margin), or d_heavy is null for every non-contact and the observation field can only ever remove contacts. tcren potts pairs --probe and tcren potts selfcon both do this.

Note

This is EM-shaped but it is not EM. The E-step is a genuine conditional expectation, but the potential update is the shipped log-odds re-derivation applied to expected counts, not the maximiser of any one objective — so the pseudo-log-likelihood need not rise, and on the 17 bundled crystals it does not (\(-337.94\) then \(-335.05\), drifting to \(-335.74\) over six steps at damping=0.5). Convergence is declared on max_abs_de, the largest cell movement in natural-log units, which contracts at the damping rate. pseudo_loglik rides along as a diagnostic.

The reduction gate. At margin = 0 the crystal’s distance read pins every site, the posterior is the observed map exactly, and the probability-weighted tally reproduces the plain unweighted one bit for bit — on the bundled crystals, max |soft − hard| = 0 over all 380 cells and a posterior mass of 309.000000 against 309 observed contacts. contact_mass and mass_from_uncertain sum unweighted site probabilities; structure weights enter the fit and potential derivation separately. That is the check that the two halves of the loop are wired to the same contacts, and it is a test rather than an approximation.

Note

That gate is exact on one site table, and the Potts availability set is not quite the derivation’s contact set. Over Native2026 the 15 Å Cα radius holds 7,799 of the deposited derivation’s 7,816 TCR:peptide contacts (99.78 %); the 17 it misses are spread over 17 structures at one or two each, side chains touching while their Cα sit further apart than the radius. So margin = 0 reproduces the derivation from the availability set, not the shipped TCRen2_potential.csv — which is derived over the full contact map.

tcren potts selfcon -s <structures> --balance both --margin 0.5 \
    -o TCRen_sc.csv --model-out potts_sc.json --trace-out sc_trace.tsv

API

A coupled Potts model over the TCR:pMHC contact map.

The configuration is the contact map itself. Sites are the residue pairs that could have contacted — Cα within radius — and sigma_a = 1 iff a heavy-atom contact formed:

E(sigma) = - sum_a eta_a sigma_a - 1/2 sum_{a,b} A_ab sigma_a sigma_b
P(sigma) = exp(-E(sigma)) / Z

eta is an additive one-body model (fields for each side, a pair coupling J, and the geometry the fields must be adjusted for); A couples neighbouring cells of the map. Fitting is penalised pseudolikelihood and needs no Z; scoring gets Z by annealed importance sampling from the uncoupled model, whose partition function is exact.

Typical use:

from tcren.potts import PottsModel, available_pairs, score_sites, contact_probabilities

pairs = available_pairs(structure)                     # the sites
model = PottsModel.bundled()                           # or fit_potts(pairs)
score_sites(pairs, model)                              # energy, log Z, likelihood
contact_probabilities(pairs, model)                    # per-site contact probability
contact_map(pairs, model, by="loop")                   # loop x peptide-position frequency map
energy_ensemble(pairs, model)                          # <Phi>, Var Phi, free energy, entropy
map_ensemble(pairs, model)                             # the configurations, with probabilities
threaded_ensemble(pairs, model, candidates)            # the same per candidate sequence

The retired sequence-only geometry experiment is preserved in the repository attic.

The marginals settle the mean interaction energy and nothing else: Phi is linear in sigma, so <Phi> = sum_a p_a phi_a exactly — but only with p_model, the coupled marginal, and only read at the sequence being scored. tcren.potts.ensemble carries the read-outs a set of marginals cannot determine, and re-reads the marginals per candidate.

The CLI mirrors it: tcren potts fit, tcren potts score, tcren potts contacts, tcren potts map, tcren potts ensemble.

class tcren.potts.PottsModel(alpha, h_rec, h_par, g_dist, g_region, g_role, g_class, kernel, coupling=None, beta_matrix=None, coupling_matrix_name=None, coupling_pinned=False, g_orient=<factory>, kernel_se=<factory>, radius=15.0, cutoff=5.0, joint=False, far=False, alphabet='ACDEFGHIKLMNPQRSTVWY', regions=<factory>, roles=<factory>, classes=<factory>, dbin=0.5, n_structures=0, n_sites=0, n_contacts=0, pseudo_loglik=nan, notes='')

Bases: object

Fitted parameters of the coupled contact-map model, plus everything needed to re-apply them.

Variables:

• alpha (float) – Intercept — the log-odds of a contact at the mean of every block.

• h_rec (list[float]) – Receptor field, one entry per AA. Positive = engages an available partner more often than average; negative = declines it.

• h_par (list[float]) – Partner field, same convention and alphabet.

• coupling (list[list[float]] | None) – Free 20x20 J(receptor aa, partner aa), double-centred, or None when coupling_matrix_name is set.

• beta_matrix (float | None) – Scale on the named potential when J is not free.

• coupling_matrix_name (str | None) – Bundled potential the coupling is fixed to (tcren2, mj, mj1996, keskin, tcren, tcren2_em), or None for a free J or for a pin on a matrix that is not one of the bundled five.

• coupling_pinned (bool) – Whether coupling holds a pinned matrix rather than a free J. A pin on an unnamed matrix — one derived in this session rather than shipped — has to store the matrix itself to round-trip, so it lands in coupling and this flag is what tells the two apart. beta_matrix is its single fitted scale either way.

• g_dist (list[float]) – Cα-distance-bin coefficients, bin b covering [b*DBIN, (b+1)*DBIN) Å.

• g_orient (list[float]) – Orientation-bin coefficients over N_ORIENT bins of tcren.geometry.orientation_cosine() — does the residue turn toward the partner. Empty in every model fitted before the block existed, and empty means the block contributes nothing, so those models round-trip and score unchanged.

• g_region (list[float]) – Receptor-region coefficients over REGIONS.

• g_role (list[float]) – Partner-role coefficients over ROLES.

• g_class (list[float]) – Partner-class coefficients over CLASSES.

• kernel (list[float]) – Coupling coefficients in the order of kernel_names().

• kernel_se (list[float]) – Cluster-robust standard errors for kernel, clustered on the structure.

• radius (float) – Availability radius, Å (Cα–Cα). A pair beyond it is not a site.

• cutoff (float) – Contact definition, Å (closest heavy atom).

• joint (bool) – Whether the model was fitted over both partners at once.

• far (bool) – Whether the long-offset exclusion class is fitted. A receptor residue’s grip on the peptide is a contiguous run of two or three positions — K(+0,+1) and K(+0,+2) are the fitted cooperativity that builds it — and nothing in the model says the run stops. This class is what says so, and it is off in every shipped model, so a model written before it round-trips unchanged.

• alphabet (str) – The residue alphabet h and J are indexed by.

• n_contacts (n_structures / n_sites /) – What the fit was estimated on.

• pseudo_loglik (float) – Unpenalized weighted pseudo-log-likelihood at the penalized optimum.

• notes (str) – Free text, e.g. the structure set.

Parameters:

• alpha (float)

• h_rec (list[float])

• h_par (list[float])

• g_dist (list[float])

• g_region (list[float])

• g_role (list[float])

• g_class (list[float])

• kernel (list[float])

• coupling (list[list[float]] | None)

• beta_matrix (float | None)

• coupling_matrix_name (str | None)

• coupling_pinned (bool)

• g_orient (list[float])

• kernel_se (list[float])

• radius (float)

• cutoff (float)

• joint (bool)

• far (bool)

• alphabet (str)

• regions (list[str])

• roles (list[str])

• classes (list[str])

• dbin (float)

• n_structures (int)

• n_sites (int)

• n_contacts (int)

• pseudo_loglik (float)

• notes (str)

alpha: float

h_rec: list[float]

h_par: list[float]

g_dist: list[float]

g_region: list[float]

g_role: list[float]

g_class: list[float]

kernel: list[float]

coupling: list[list[float]] | None = None

beta_matrix: float | None = None

coupling_matrix_name: str | None = None

coupling_pinned: bool = False

g_orient: list[float]

kernel_se: list[float]

radius: float = 15.0

cutoff: float = 5.0

joint: bool = False

far: bool = False

alphabet: str = 'ACDEFGHIKLMNPQRSTVWY'

regions: list[str]

roles: list[str]

classes: list[str]

dbin: float = 0.5

n_structures: int = 0

n_sites: int = 0

n_contacts: int = 0

pseudo_loglik: float = nan

notes: str = ''

coupling_array()

J as a 20x20 array, whether it was fitted freely or fixed to a matrix.

Return type:

ndarray

n_parameters()

Stored fitted coefficients, including redundant/unobserved levels.

Not identifiable degrees of freedom; a pinned matrix contributes its single scale.

Return type:

int

to_json(path)

Write the model as JSON (round-trips through from_json()).

Parameters:

path (str | Path)

Return type:

None

classmethod from_json(path)

Read a model written by to_json().

Parameters:

path (str | Path)

Return type:

PottsModel

classmethod bundled(name='potts_tcr_peptide')

Load a model shipped under tcren/data (default: the TCR:peptide interface).

Parameters:

name (str)

Return type:

PottsModel

tcren.potts.centred_potential(name, alphabet=('A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y'), *, centre=True)

-M onto alphabet, double-centred by default — positive means more contact.

A Potential is signed negative is favourable; a coupling here is a log-odds where positive is more likely, hence the sign flip.

The two settings answer two different questions and are not interchangeable.

centre=True (default) removes every one-body term, so a matrix used as a fixed coupling contributes nothing to the single-residue marginals — the fields carry those, refitted freely, and competing matrices are compared on pair structure alone. This is the setting for ranking potentials against each other.

centre=False keeps the raw matrix. Use it to reproduce a referenced contact-map score: the peptide-referenced energy of tcren.ddg.reference_delta() is a difference of one-body sums, so double-centring re-injects a burial-scaled composition term n_i·c(a) — the position’s contact count times the potential’s partner-residue column mean — and the identity fails. Pinned uncentred, the model’s referenced energy equals the potential’s own up to the fitted scale.

Parameters:

• name – One of tcren2, tcren2_em, tcren, mj, mj1996, keskin; or a Potential — including one derived in this session, which is what tcren.potts.selfcon.fit_selfconsistent() pins on each step and is the reason this is not a name-only lookup.

• alphabet (tuple[str, ...]) – Residue order of the returned array.

• centre (bool) – Double-centre the matrix. See above for which setting a task wants.

Returns:

A (len(alphabet), len(alphabet)) array; cells absent from the potential are 0.

Return type:

ndarray

tcren.potts.kernel_names(joint, far=False)

Names of the coupling classes, in the order tcren.potts.edges() returns them.

Parameters:

• joint (bool)

• far (bool)

Return type:

list[str]

tcren.potts.available_pairs(structure, partner='peptide', *, receptor='tcr', radius=15.0, cutoff=5.0, probe=None)

Every receptor:partner residue pair inside radius, with whether it contacted.

Three arms, and they are three separate Hamiltonians, not one field read three ways: the receptor’s grip on the peptide, its grip on the groove, and the groove’s grip on the peptide. The last is the presentation term, and it is the one an activation read-out cannot do without — a peptide whose anchors are destroyed never presents, so it never activates, whatever the receptor would have done with it.

Parameters:

• structure (Structure) – A chain-typed structure. Whenever the MHC is on either side it must also carry groove regions — run tcren.mhc.annotate_mhc() (or annotate_mhc_batch over a set) first, or every MHC residue is dropped for want of a groove region.

• partner (str) – "peptide" or "mhc".

• receptor (str) – "tcr" (default) or "mhc". receptor="mhc", partner="peptide" is the peptide:MHC arm: the groove residue takes the receptor slot, the peptide residue the partner slot, and region.rec is a groove region rather than a CDR loop — so a model fitted here carries MHC_RECEPTOR_REGIONS, not the TCR loop set.

• radius (float) – Availability radius, Å, on the Cα–Cα distance.

• cutoff (float) – Contact definition, Å, on the closest heavy-atom distance.

• probe (float | None) – Radius the heavy-atom layer is measured out to. None (default) uses radius, so every available pair carries a real d_heavy. Measured on 3UTS the default reaches 0 nulls over 171 pairs where a cutoff-wide probe left 154 (90.1 %), at no measurable cost, and sigma is identical either way — the contact call keys on cutoff and never on the probe. A narrower probe is still accepted and is clamped up to cutoff; it makes d_heavy null for every pair that missed, which tcren.potts.selfcon.observation_field() reads as a confident non-contact. tcren.potts.em needs the layer to reach past cutoff in both directions, because a latent band that only reaches inwards can turn contacts off and never on.

Returns:

pdb.id, the two residue types and positions, the receptor chain and region, the partner role, the partner class, both distances and sigma. pos.par is a within-region offset on both partners (the peptide chain carries one region starting at 0, so there it is the plain 0-based peptide position), because the coupling kernel is defined on sequence offsets inside one loop or one helix.

Return type:

One row per site

Example

>>> from tcren.potts import available_pairs
>>> pairs = available_pairs(structure)
>>> pairs["sigma"].mean()

tcren.potts.site_codes(sites, model=None, *, radius=15.0, dbin=0.5, regions=None, orient=None)

Integer code arrays per one-body block, plus the block sizes and the annotated frame.

Bin edges are global (base 0), never derived from the frame’s own minimum: a parameter vector fitted on one structure set has to index the same bins on every other, or every score is read against the wrong coefficients.

Returns:

(codes, sizes, frame) — codes is one integer array per block in the order receptor residue, partner residue, distance bin, receptor region, partner role, partner class, and — when the model carries a fitted g_orient and the frame carries cos_cb — orientation bin; sizes the number of levels of each; frame the input with loop, pchain and a global sid added.

Parameters:

• sites (DataFrame)

• model (PottsModel | None)

• radius (float)

• dbin (float)

• regions (tuple[str, ...] | None)

• orient (bool | None)

tcren.potts.eta(codes, model)

The one-body log-odds of a contact at every site, under model.

Parameters:

model (PottsModel)

Return type:

ndarray

tcren.potts.edges(q, joint=False, far=False)

One (n, 2) array of site-index pairs per coupling class, in kernel_names() order.

Parameters:

• q (DataFrame) – The frame tcren.potts.site_codes() returns, carrying sid, loop and pchain.

• joint (bool) – Include the cross-class family (both partners in one model).

• far (bool) – Include the long-offset exclusion family, _edges_far().

Return type:

list[ndarray]

tcren.potts.neighbour_counts(edge_lists, sigma, n_sites)

n_k(a) — contacting neighbours of each site in each class, both directions of each edge.

These are the covariates of the pseudolikelihood: the conditional P(sigma_a = 1 | sigma_{-a}) is logistic in eta_a + sum_k K_k n_k(a).

Parameters:

• edge_lists (list[ndarray])

• sigma (ndarray)

• n_sites (int)

Return type:

ndarray

tcren.potts.bucket_edges(edge_lists, starts, kv)

Group every non-zero-coefficient edge by structure. Returns (E, offsets).

Site indices are contiguous per structure, so one searchsorted assigns each edge to its structure. E has columns (sid_a, sid_b, class).

Parameters:

• edge_lists (list[ndarray])

• starts (ndarray)

• kv (ndarray)

tcren.potts.colour(n_sites, ea, eb)

Greedy graph colouring of the coupling graph, largest-degree first.

Sites sharing a colour have no edge between them, so they are conditionally independent given the rest and a whole colour class updates in one vectorised Gibbs step. The property is asserted against the real edge lists wherever it is relied on, never argued.

Parameters:

• n_sites (int)

• ea (ndarray)

• eb (ndarray)

Return type:

list[ndarray]

tcren.potts.coupling_matrix(n_sites, ea, eb, ec, kv)

Dense symmetric A with zero diagonal: A[a,b] is the coefficient of the class linking a and b.

Parameters:

• n_sites (int)

• kv (ndarray)

Return type:

ndarray

tcren.potts.design(codes, sizes, counts, fixed_coupling=None)

Sparse row-level design [intercept | one-hot blocks | J | neighbour counts].

fixed_coupling replaces the free 400-cell J block with a single column carrying the named potential’s value at that cell, so its coefficient is the scale beta_matrix.

Returns:

(X, block_slices, k0) with k0 the first coupling-coefficient column.

Parameters:

• counts (ndarray)

• fixed_coupling (ndarray | None)

tcren.potts.irls(X, y, w, ridge, *, free=(0,), tol=1e-10, maxit=200)

Penalised IRLS for a weighted Bernoulli GLM on a sparse design. Convex; Newton steps.

Returns (coefficients, penalised Hessian). free names columns exempt from the ridge.

Parameters:

• y (ndarray)

• w (ndarray)

• ridge (float)

• tol (float)

• maxit (int)

tcren.potts.gauge(b, slices, free_coupling)

Project to the zero-sum (Ising) gauge, in place on a copy.

Every one-body block is centred and J double-centred, the displaced level absorbed by the intercept, so eta is unchanged while as much as possible sits in the fields and as little as necessary in the couplings (Cocco et al., arXiv:1703.01222). A fixed-matrix coupling is a scale on an already double-centred matrix and is never re-centred.

Parameters:

• b (ndarray)

• free_coupling (bool)

tcren.potts.cluster_se(X, resid, gid, n_groups, H)

Sandwich standard errors clustered on the structure.

A structure contributes hundreds of correlated sites, so the model-based (X'WX)^-1 understates every standard error.

Parameters:

• resid (ndarray)

• gid (ndarray)

• n_groups (int)

Return type:

ndarray

tcren.potts.fit_potts(sites, *, radius=15.0, cutoff=5.0, couplings=True, coupling_matrix=None, pin_centred=True, weights=None, ridge=1.0, joint=None, far=False, orient=False, regions=None, notes='')

Fit the model to a table of available pairs.

Parameters:

• sites (DataFrame) – Rows from tcren.potts.available_pairs(), one or both partners concatenated.

• radius (float) – Availability radius, Å. Sites beyond it are dropped.

• cutoff (float) – Contact definition, Å — recorded on the model; sigma is already computed.

• couplings (bool) – Fit the sigma–sigma kernel. False gives the factorised model, whose partition function is then exact.

• coupling_matrix (str | object | None) – Fix J to one scale on this potential — a bundled name (tcren2, mj, …) or a Potential object, which is how a matrix derived in this session gets pinned instead of fitting 400 free cells. Competing matrices fitted this way carry identical parameter counts, so their pseudo-log-likelihoods compare directly.

• pin_centred (bool) – Double-centre coupling_matrix before pinning (default). False pins the raw matrix, which is what makes the model’s peptide-referenced energy reduce to the potential’s own; see tcren.potts.centred_potential() for which setting a task wants. Ignored when coupling_matrix is None.

• weights (dict[str, float] | None) – Per-structure weights, e.g. from tcren.potential.balanced_weights(), to down-weight redundancy. Missing ids get 1.

• ridge (float) – ℓ2 penalty on every coefficient but the intercept.

• regions (tuple[str, ...] | None) – Receptor-region level set for the g_region block. None (default) is the TCR loop set. The peptide:MHC arm passes tcren.potts.MHC_RECEPTOR_REGIONS, because there the receptor is the groove and every TCR loop level would be unused while every groove region collapsed into other.

• far (bool) – Include the long-offset exclusion class, X(same receptor residue, |dj| > 2). Off by default, so every shipped model is unaffected.

• orient (bool) – Fit the g_orient block — does the receptor residue turn toward the partner. Needs a cos_cb column on sites, which available_pairs() emits. Held at a fixed Cα distance the cosine separates contact from non-contact at ROC-AUC 0.69-0.84 over 4-12 Å, so this is information the distance block does not carry. Off by default, so every shipped model is unaffected.

• joint (bool | None) – Include the cross-class coupling family. Defaults to whether both partner classes are present.

• notes (str) – Free text stored on the model.

Returns:

A PottsModel in the zero-sum gauge.

Return type:

PottsModel

Example

>>> from tcren.potts import available_pairs, fit_potts
>>> pairs = pl.concat([available_pairs(s) for s in structures])
>>> model = fit_potts(pairs, notes="Native2026")
>>> model.to_json("potts.json")

tcren.potts.kernel_table(model)

The coupling coefficients with their cluster-robust standard errors, strongest first.

A coefficient is the log-odds added to a site’s field per contacting neighbour of that class. On TCR:peptide crystals every axial class comes out positive and every off-axis class negative: a made contact recruits its own sequence neighbours onto the same partner residue and suppresses the diagonal one.

Parameters:

model (PottsModel)

Return type:

DataFrame

tcren.potts.gibbs(eta, A, colours, rng, *, chains=64, burn=100, draws=100, thin=3, observer=None, init=None)

Block Gibbs at beta = 1.

Parameters:

• eta (ndarray) – One-body log-odds per site.

• A (ndarray) – Symmetric coupling matrix with zero diagonal.

• colours – Colour classes from tcren.potts.colour(). A subset is accepted only with explicit init; omitted sites remain clamped at that start.

• rng – A numpy.random.Generator.

• chains (int) – Independent chains run in parallel (vectorised, not looped).

• burn (int) – Sweeps discarded before the first draw.

• draws (int) – Draws kept per chain.

• thin (int) – Sweeps between kept draws.

• observer – Optional callable invoked on each kept draw with the (chains, n) matrix of configurations. Use it to accumulate a statistic of whole configurations without materialising every draw — the matrix is reused between calls, so copy what you keep.

• init – Optional binary start, shape (n,) broadcast to every chain, or (chains, n) retaining separate chain configurations for warm starts. None (default) draws the start from the factorised model, which is what every existing caller gets and leaves the random stream untouched. Passing the observed map starts the chain at the structure’s own contact set; the stationary distribution is the same either way, so this changes how fast the chain forgets its start and nothing else.

Returns:

(occupancy, totals) — the mean of sigma per site over all kept draws, and the contact total of every kept draw.

tcren.potts.ais_log_z(eta, A, colours, rng, *, particles=64, steps=256)

log Z by annealed importance sampling. Returns (log_z, effective_sample_size).

The beta = 0 reference is the factorised model, sampled exactly, so no burn-in is needed and the estimator is unbiased in Z. An effective sample size close to particles means the annealing schedule is long enough; a small one means it is not.

Parameters:

• eta (ndarray)

• A (ndarray)

• particles (int)

• steps (int)

Return type:

tuple[float, float]

tcren.potts.exact_log_z(eta, A)

log Z by brute force over all 2^n configurations. For testing; n <= 22.

Parameters:

• eta (ndarray)

• A (ndarray)

Return type:

float

tcren.potts.factorised_log_z(eta)

log Z of the uncoupled model — exact, in closed form.

Parameters:

eta (ndarray)

Return type:

float

tcren.potts.energy(sigma, eta, A)

E(sigma) = -eta . sigma - 1/2 sigma' A sigma. Lower is more favourable.

Parameters:

• sigma (ndarray)

• eta (ndarray)

• A (ndarray)

Return type:

float

tcren.potts.tilt_mean(totals, mu)

<N> under the tilted model E_mu = E - mu N, by reweighting draws taken at mu = 0.

Exact in the tilted family, because every tilt considered here depends on sigma only through N(sigma): the importance weight is exp(mu N) and needs no new sampling.

Parameters:

• totals (ndarray)

• mu (float)

Return type:

float

tcren.potts.mu_star(totals, n_obs, *, span=40.0, tol=1e-08)

The chemical potential at which the model’s mean contact count matches n_obs.

<N>_mu is non-decreasing in mu (its derivative is the tilted variance), so a bisection is exact. Returns nan when n_obs lies outside the sampled support, where reweighting would be extrapolating rather than reweighting.

Parameters:

• totals (ndarray)

• n_obs (float)

• span (float)

• tol (float)

Return type:

float

tcren.potts.count_free_energy(totals)

(N, F(N)) with F(N) = -log p(N) from the sampled contact-count histogram.

Related to log Z(mu) by log-sum-exp over counts; a Legendre transform only under a saddlepoint approximation. This is the same object the protein-folding literature draws along the native-contact coordinate. Only the sampled range is returned: the histogram cannot see the unbound basin of a docked pose, which is what delta_f_empty() is for.

Parameters:

totals (ndarray)

Return type:

tuple[ndarray, ndarray]

tcren.potts.delta_f_empty(log_z)

log[P(N >= 1) / P(N = 0)] = log(Z - 1) — exact, since E(empty) = 0.

The whole-interface reading of the two-state contrast that eta_a is for a single site: the free energy of making any contact at all against making none.

Parameters:

log_z (float)

Return type:

float

tcren.potts.delta_f_threshold(totals, x)

log[P(N >= x) / P(N < x)] under the model. Z cancels, so this needs no AIS.

Parameters:

• totals (ndarray)

• x (int)

Return type:

float

tcren.potts.bound_unbound(sites, model, *, threshold=None, chains=64, burn=100, draws=100, thin=3, particles=64, steps=256, seed=0, workers=None)

The whole-interface two-state free energy, in its three readings.

A single site has two states, sigma_a = 0 and 1, and eta_a is the free-energy difference between them. The same contrast for the whole interface needs a macrostate, and the contact count N(sigma) is the coordinate that defines one. Because every tilt in N is an exponential family, one Gibbs pass serves all three readings.

Columns, one row per structure:

neg_energy

\(-E(\sigma_{\mathrm{obs}})\), higher is more native-like. The \(\Pi\) block of tcren.reliability.s_score(), emitted here so one call supplies both the free-energy readings and the energy term the reliability score needs.

df_empty

log[P(N >= 1) / P(N = 0)] = log(Z - 1), exact, from the AIS log Z. E(empty) = 0, so the empty configuration needs no separate estimate. This is the capacity of the interface: what it would gain by making any contact at all.

df_threshold

log[P(N >= x) / P(N < x)] at threshold, from the sampled histogram. Z cancels, so this needs no AIS — but it is only finite while x sits inside the sampled range.

mu_star

the tilt at which <N>_mu equals the observed count: how much extra chemical potential the model needs to explain the map it was shown. Positive means the structure made more contacts than its fields and couplings warrant. nan outside the sampled support.

n_lo, n_hi, n_mean, n_var

the sampled contact-count range, mean and variance, so every nan above is auditable.

Note

df_empty and df_threshold are not two estimates of one number. The unbound basin of a docked pose is astronomically improbable — the model is conditioned on an available set that already holds the receptor against the peptide — so no sampler reaches N = 0 and only the log Z route gives it. The threshold reading is a local contrast inside the bound basin.

Parameters:

• sites (DataFrame) – Rows from tcren.potts.available_pairs().

• model (PottsModel) – A fitted PottsModel.

• threshold (int | None) – x for df_threshold. None disables that column.

• chains (int) – Gibbs settings; see tcren.potts.gibbs().

• burn (int) – Gibbs settings; see tcren.potts.gibbs().

• draws (int) – Gibbs settings; see tcren.potts.gibbs().

• thin (int) – Gibbs settings; see tcren.potts.gibbs().

• particles (int) – AIS settings; see tcren.potts.ais_log_z().

• steps (int) – AIS settings; see tcren.potts.ais_log_z().

• seed (int) – Seed for both samplers.

• workers (int | None) – Processes to split the structures across. None takes every core, 1 runs serially. The result is identical either way — each structure’s numbers depend on (seed, pdb.id) alone, never on how the work was split.

Returns:

One row per structure.

Return type:

DataFrame

tcren.potts.count_profile(sites, model, *, chains=64, burn=100, draws=100, thin=3, seed=0)

Pooled F(N) = -log p(N) over every structure in sites, plus the observed counts.

The free-energy profile along the contact-count coordinate. Use it to see whether the model’s contact-count landscape has a barrier – if it does not, a threshold reading of the two-state contrast has nothing to key on and mu_star is the meaningful statistic.

Parameters:

• sites (DataFrame)

• model (PottsModel)

• chains (int)

• burn (int)

• draws (int)

• thin (int)

• seed (int)

Return type:

DataFrame

tcren.potts.score_sites(sites, model, *, particles=64, steps=256, seed=0, workers=None)

Energy, log Z and likelihoods for every structure in a table of available pairs.

Columns, one row per structure:

n_sites / n_contacts

available pairs and the contacts among them.

energy / neg_energy

E(sigma_obs), the Hamiltonian of the observed contact map, and its negation. Lower energy is more favourable, so neg_energy is the one that runs the same way as every other score here: higher is more native-like. It is the \(\Pi\) block of tcren.reliability.s_score().

log_z / log_z0

the coupled partition function by AIS, and the factorised one in closed form. log_z0 is the free energy of the available set before any contact is observed.

log_lik / psi

log P(sigma_obs) and the same divided by n_sites, so it compares across interfaces of different size.

pseudo_log_lik / psi_pseudo

the exact, MCMC-free cross-check.

ais_ess

effective sample size of the AIS weights, out of particles. Close to particles means the annealing schedule was long enough.

Parameters:

• sites (DataFrame) – Rows from tcren.potts.available_pairs(), one or many structures.

• model (PottsModel) – A fitted PottsModel.

• particles (int) – AIS particles per structure.

• steps (int) – AIS annealing steps.

• seed (int) – Seed for the sampler.

• workers (int | None) – Processes to split the structures across. None takes every core, 1 runs serially. The result is identical either way — each structure’s numbers depend on (seed, pdb.id) alone, never on how the work was split.

Return type:

DataFrame

tcren.potts.contact_probabilities(sites, model, *, chains=64, burn=100, draws=100, thin=3, seed=0, workers=None)

Per-site contact probability under the model, beside the observed indicator.

Three probabilities per site, because the difference between them is the couplings:

p_independent

1 / (1 + exp(-eta)) — the one-body model alone, ignoring every other site.

p_model

the marginal <sigma_a> of the full coupled model, by block Gibbs. This is the number to use.

p_conditional

P(sigma_a = 1 | the observed rest) — the pseudolikelihood conditional, which reads “given everything else this structure actually did, how likely was this contact?”.

Parameters:

• sites (DataFrame) – Rows from tcren.potts.available_pairs().

• model (PottsModel) – A fitted PottsModel.

• chains (int) – Gibbs settings; see tcren.potts.gibbs().

• burn (int) – Gibbs settings; see tcren.potts.gibbs().

• draws (int) – Gibbs settings; see tcren.potts.gibbs().

• thin (int) – Gibbs settings; see tcren.potts.gibbs().

• seed (int) – Seed for the sampler.

• workers (int | None) – Processes to split the structures across. None (default) uses every core, 1 runs serially. Bit-identical either way; see tcren.potts.score_sites().

Returns:

The input frame with the three probabilities appended, one row per site.

Return type:

DataFrame

tcren.potts.map_agreement(sites, *, column='p_model', observed='sigma')

How well a predicted contact map agrees with the one a structure actually made.

Two numbers, because they answer different questions and no single one answers both:

js

Placement. Jensen–Shannon similarity between the two maps read as distributions over sites, \(q_a = \sigma_a / \sum \sigma\) and \(r_a = p_a / \sum p\):

\[\mathrm{js} = 1 - \tfrac12 D_{\mathrm{KL}}(q \| m) - \tfrac12 D_{\mathrm{KL}}(r \| m), \qquad m = \tfrac12 (q + r)\]

in bits, so it lies in \([0, 1]\): 1 is identical placement, 0 disjoint support. It is scale-free in the level by construction — normalising divides the total out — which is why ratio is reported beside it and not folded in.

\(\mathrm{js} = 1\) requires the model to put all its mass on exactly the contacting cells in equal proportion. This also holds for \(p_a=c\sigma_a\), with any \(0<c<1\); it does not require a deterministic map. The level ratio distinguishes these proportional maps.

ratio

Level. \(\sum_a p_a / \sum_a \sigma_a\). One is the count reproduced.

auc is reported too, as the ordering read alone.

Not R². On a binary outcome 1 - SS_res/SS_tot is the Brier skill score — they agree to three decimals on every map we have checked — so reporting both is reporting one number twice. It is also dominated by the many easy zeros and scales with prevalence, so two maps with different site counts are not comparable under it. That is exactly the comparison this function exists for, which is why the placement read-out here is normalised and symmetric instead.

Parameters:

• sites (DataFrame) – A scored site table — from contact_probabilities(), or tcren potts contacts. Grouped by pdb.id; rows with a null observed are read as non-contacts.

• column (str) – The probability column to score.

• observed (str) – The observed indicator column.

Returns:

sites, contacts, expected, ratio, js, auc.

Return type:

One row per pdb.id

tcren.potts.contact_map(sites, model, *, by='loop', chains=64, burn=100, draws=100, thin=3, seed=0, workers=None)

Predicted contact frequency, closed from per-pair probabilities onto a coarser grid.

contact_probabilities() gives p_model, the coupled model’s marginal for a single receptor-residue : peptide-residue pair. Two coarser readings are what an experiment actually measures:

by="loop"

one row per (structure, CDR loop, peptide position) — the contact-frequency map, the grid a molecular-dynamics trajectory reports as the fraction of frames in which any residue of that loop touches that peptide position.

by="position"

one row per (structure, peptide position) — peptide residue importance: how engaged the model expects that position to be, before any residue identity is scored.

by="pair"

the ungrouped table, exactly contact_probabilities().

The returned p_any uses an independence approximation to the coupled marginals. Under independence the contact count is Poisson-binomial and the event “at least one” is:

\[P(N \ge 1) \;=\; 1 - \prod_j (1 - p_j)\]

where \(p_j\) is the model marginal of pair \(j\) in the group and \(N\) the number of contacts the group makes. Coupled marginals alone do not determine this union probability; an exact coupled union requires joint configuration samples. It is accumulated in \(\log(1 - p)\) so a twelve-residue loop does not underflow, and a pair at \(p_j = 1\) forces the group to 1 exactly rather than to nan.

These are contact frequencies: dimensionless, in \([0, 1]\), higher meaning more often in contact. They are not energies and carry no \(k_\mathrm{B}T\), so nothing here belongs in an energy block — this is the diagnostic and importance side of the model.

Parameters:

• sites (DataFrame) – Rows from tcren.potts.available_pairs().

• model (PottsModel) – A fitted PottsModel.

• by (str) – "loop" (default), "position" or "pair".

• chains (int) – Gibbs settings; see tcren.potts.gibbs().

• burn (int) – Gibbs settings; see tcren.potts.gibbs().

• draws (int) – Gibbs settings; see tcren.potts.gibbs().

• thin (int) – Gibbs settings; see tcren.potts.gibbs().

• seed (int) – Seed for the sampler.

• workers (int | None) – Processes to split the structures across; None uses every core, 1 runs serially. Bit-identical either way.

Returns:

For "pair", contact_probabilities()’ frame unchanged. Otherwise the grouping columns plus p_any (the independence approximation above), p_expected (the expected number of contacts in the group, \(\sum_j p_j\)), n_pairs (available pairs in the group), n_observed (how many of them this structure made) and observed (1 if any did, else 0 — the indicator the prediction is scored against).

Return type:

DataFrame

Example

>>> from tcren.potts import available_pairs, contact_map, PottsModel
>>> m = contact_map(available_pairs(structure), PottsModel.bundled())

tcren.potts.peptide_free_energy(sites, model, *, coupled=False, marginals=None, chains=64, burn=100, draws=100, thin=3, seed=0, workers=None)

Historical log-partition substitution readout, not a calibrated binding free energy.

contact_map() reads how engaged a position is expected to be before any residue identity is scored. This reads what happens when the identity changes. The partner residue enters the one-body field twice — through the partner propensity \(h^{\mathrm{par}}\) and through the pair term \(J\) — so substituting position \(i\) shifts \(\eta\) at every available pair carrying that position, and the uncoupled log-partition changes with it:

\[\Phi^{\mathrm{Potts}}(x) \;=\; \log Z_0\big(\eta(x)\big) \;=\; \sum_a \log\!\big(1 + e^{\eta_a(x)}\big), \qquad \Delta\Phi_i(a) \;=\; \Phi^{\mathrm{Potts}}(x_{i \to a}) \;-\; \tfrac{1}{20}\sum_b \Phi^{\mathrm{Potts}}(x_{i \to b})\]

Higher is more favourable: \(\log Z_0\) is the interface’s capacity to make contacts at all, so a residue that raises it engages more. The reference is the equimolar one — the mean over the twenty residues at that position, not the residue the structure happens to carry — which is the null a positional-scanning library actually holds the other positions at.

log Z0 is dimensionless, in nats. The uncoupled statistical free energy is -log Z0, with the opposite sign; neither is a calibrated binding free energy. These fitted-model readouts are historical, excluded from final fit-free scoring.

Two readings, from the same fields:

coupled=False (default)

\(\log Z_0\) for the coupling-free model, which is exact and closed form — no sampling.

coupled=True

linear response about the observed sequence. Since \(\partial \log Z / \partial \eta_a = \langle\sigma_a\rangle\), the coupled log-partition changes to first order as \(\Delta \log Z \approx \sum_a p_a \Delta\eta_a\) with \(p_a\) the marginal of contact_probabilities() — one Gibbs pass, then a dot product per cell. This is not an exact finite-substitution coupled log-partition difference. The retained log_z0 column then stores the frozen-marginal linear surrogate, not an estimate of the coupled normalizer.

Only aa.par changes: the backbone, the Cα distances, the receptor residues and the partner roles are the structure’s own and are held fixed, which is the same fixed-backbone approximation every threading score in the package makes.

Parameters:

• sites (DataFrame) – Rows from tcren.potts.available_pairs().

• model (PottsModel) – A fitted PottsModel.

• coupled (bool) – Take the linear-response path against the coupled marginals.

• marginals (DataFrame | None) – A frame from contact_probabilities() to reuse when coupled=True. None computes it.

• chains (int) – Gibbs settings, used only when coupled=True.

• burn (int) – Gibbs settings, used only when coupled=True.

• draws (int) – Gibbs settings, used only when coupled=True.

• thin (int) – Gibbs settings, used only when coupled=True.

• seed (int) – Seed for the sampler.

• workers (int | None) – Processes to split the structures across when coupled=True.

Returns:

One row per (pdb.id, pos.par, aa.par) with log_z0 (the whole-interface \(\log Z_0\) under that substitution), dF (its equimolar-referenced effect), n_pairs (available pairs carrying that position) and is_observed (1 for the residue the structure carries). dF sums to zero over the twenty residues at every position.

Return type:

DataFrame

Example

>>> from tcren.potts import available_pairs, peptide_free_energy, PottsModel
>>> peptide_free_energy(available_pairs(structure), PottsModel.bundled())

tcren.potts.connected_correlations(sites, model, *, chains=64, burn=100, draws=100, thin=3, seed=0)

The generative test: two-point correlations in the data against in maps sampled from it.

For each coupling class, pooled over every edge in it and every structure,

\[c_\Delta = \langle \sigma_a \sigma_b \rangle - \langle \sigma_a \rangle \langle \sigma_b \rangle\]

computed once from the observed contact maps and once from maps drawn from the model. This is the criterion Boltzmann-machine DCA trains to, and a pseudolikelihood fit never sees it — the conditionals it maximises are one-site quantities — so agreement here is a real test rather than a restatement of the objective.

Returns:

class, n_edges, c_data, c_model.

Return type:

One row per coupling class

Parameters:

• sites (DataFrame)

• model (PottsModel)

• chains (int)

• burn (int)

• draws (int)

• thin (int)

• seed (int)

tcren.potts.sample_maps(sites, model, *, chains=64, burn=100, draws=100, thin=3, seed=0)

Contact totals of maps drawn from the model, one row per structure per draw.

The generative check: a model that reproduces a real interface must reproduce the spread of its contact count, not only the mean.

Parameters:

• sites (DataFrame)

• model (PottsModel)

• chains (int)

• burn (int)

• draws (int)

• thin (int)

• seed (int)

Return type:

DataFrame

tcren.potts.score_structure(structure, model=None, *, partner='peptide', **kwargs)

Convenience: enumerate one structure’s available pairs and score them.

Parameters:

• structure – A chain-typed structure (MHC-annotated as well, for partner="mhc").

• model (PottsModel | None) – A fitted model; defaults to the bundled TCR:peptide one.

• partner (str) – "peptide" or "mhc".

• **kwargs – Passed to score_sites().

Return type:

dict

tcren.potts.energy_ensemble(sites, model, potential=None, *, beta=1.0, anchor='observed', chains=64, burn=100, draws=100, thin=3, particles=64, steps=256, seed=0, workers=None)

Ensemble read-outs of the interaction energy, one row per structure.

The mean is the number tcren.potts.contact_probabilities() already determines; the variance, the free energy and the entropy are not, and are what a joint model buys over a set of marginals.

column	units	what it is
phi_mean	potential units	\(\langle\Phi\rangle\)
phi_mean_from_marginals	potential units	\(\sum_a p_a \varphi_a\); equal to phi_mean up to Monte-Carlo error
phi_var	(potential units)²	\(\operatorname{Var}\Phi\)
phi_free	potential units	\(-\beta^{-1}\log\langle e^{-\beta\Phi}\rangle\)
map_entropy	nats	\(\langle E\rangle + \log Z\)
hamming_obs	count	mean distance from the observed map
n_undefined	count	sites the potential does not define

Parameters:

• sites (DataFrame) – Rows from tcren.potts.available_pairs().

• model (PottsModel) – A fitted PottsModel.

• potential (Potential | None) – The interaction matrix \(\varphi\) is read from. Default TCRen2.

• beta (float) – Inverse temperature in phi_free. beta -> 0 recovers the mean.

• anchor (str) – "observed" starts every chain at the structure’s own contact map — the ensemble of maps that “do not differ much from \(\delta_{ij}\)” — and hamming_obs reports how far it travels. "random" starts from the factorised model. The stationary distribution is the same; only the burn-in differs.

• chains (int) – Gibbs settings; see tcren.potts.gibbs().

• burn (int) – Gibbs settings; see tcren.potts.gibbs().

• draws (int) – Gibbs settings; see tcren.potts.gibbs().

• thin (int) – Gibbs settings; see tcren.potts.gibbs().

• particles (int) – Annealed-importance settings for log Z; see tcren.potts.ais_log_z().

• steps (int) – Annealed-importance settings for log Z; see tcren.potts.ais_log_z().

• seed (int) – Seed for the sampler.

• workers (int | None) – Processes to split the structures across. Bit-identical either way.

Return type:

DataFrame

Note

anchor="observed" is not a Hamming-ball constraint and deliberately so. An unconstrained chain started at the observed map is still sampling \(P(\sigma)\); restricting it to a ball would sample something else and would need its own justification. hamming_obs is the diagnostic such a constraint would have been reached for.

tcren.potts.map_ensemble(sites, model, potential=None, *, by='position', top_k=20, chains=64, burn=100, draws=100, thin=3, seed=0, workers=None)

The plausible contact configurations, each with its probability.

One row per (structure, rank), ranked by how often the configuration was drawn: which maps the interface is choosing between, how much probability the most likely one carries, whether the crystal’s own map is among them (is_observed), and the mean interaction energy conditional on each.

A configuration has to be coarse-grained to carry a probability at all, and this is measured, not assumed. On 1ao7’s TCR:peptide interface (174 available sites, 6,400 draws) the raw site-level map recurs never — 6,400 distinct configurations out of 6,400 draws, every one at \(\hat p = 1.6 imes10^{-4}\), because the state space is \(2^{174}\). Coarse-grained to which peptide positions are engaged at all, the same draws collapse to 41 distinct patterns with the most likely at \(\hat p = 0.3219\). So:

by="position" (default)

which peptide positions the receptor engages. 9 groups on a 9-mer; resolvable.

by="loop"

which (CDR loop, peptide position) cells engage — the granularity the model’s molecular-dynamics validation works at. 39 groups on 1ao7, 3,308 distinct patterns, top \(\hat p = 0.0216\): at the edge of what 6,400 draws resolve.

by="pair"

the raw sites. Only meaningful for an interface small enough to enumerate, where tcren.potts.exact_log_z() gives the exact distribution instead.

pattern is the bit string over that structure’s groups in engaged-label order; engaged names the groups that are on.

Parameters:

• sites (DataFrame) – Rows from tcren.potts.available_pairs().

• model (PottsModel) – A fitted PottsModel.

• potential (Potential | None) – Read phi from. Default TCRen2.

• by (str) – Coarse-graining, one of "position", "loop", "pair".

• top_k (int) – Configurations kept per structure, most probable first.

• chains (int) – as energy_ensemble().

• burn (int) – as energy_ensemble().

• draws (int) – as energy_ensemble().

• thin (int) – as energy_ensemble().

• seed (int) – as energy_ensemble().

• workers (int | None) – as energy_ensemble().

Return type:

DataFrame

tcren.potts.threaded_ensemble(sites, model, candidates, potential=None, *, beta=1.0, coupled=True, chains=32, burn=100, candidate_burn=10, draws=50, thin=2, seed=0)

Ensemble read-outs with the contact map re-equilibrated for each candidate sequence.

The shipped soft-contact path — ddg(weights=...), cpl.predict_response(tcr_weights=...) — takes one marginal vector computed on the template and applies it to every candidate, giving \(\sum_a p_a(\text{template})\,\varphi_a(x)\). That is not the ensemble mean: a candidate residue enters \(\eta_a\) through \(h^{\mathrm{par}}\) and \(J\), so the contact ensemble itself moves under substitution. This computes \(\sum_a p_a(x)\,\varphi_a(x)\) and the non-linear read-outs beside it.

Parameters:

• sites (DataFrame) – Rows from tcren.potts.available_pairs().

• model (PottsModel) – A fitted PottsModel.

• candidates (DataFrame) – Frame with pdb.id and peptide. Every candidate must match the structure’s own par.len.

• potential (Potential | None) – Read phi from. Default TCRen2.

• beta (float) – Inverse temperature in phi_free.

• coupled (bool) – True re-runs Gibbs per candidate. False holds the marginals at the template and only re-gathers \(\varphi\), which is the shipped approximation and is here so the two can be measured against each other on identical rows.

• chains (int) – Gibbs settings. burn applies to the first candidate of a structure; later candidates warm-start from the previous one and use candidate_burn.

• burn (int) – Gibbs settings. burn applies to the first candidate of a structure; later candidates warm-start from the previous one and use candidate_burn.

• draws (int) – Gibbs settings. burn applies to the first candidate of a structure; later candidates warm-start from the previous one and use candidate_burn.

• thin (int) – Gibbs settings. burn applies to the first candidate of a structure; later candidates warm-start from the previous one and use candidate_burn.

• candidate_burn (int) – Sweeps between candidates once the chain is warm.

• seed (int) – Seed for the sampler.

Returns:

One row per (pdb.id, peptide) with phi_mean, phi_var, phi_free, n_contacts_model and n_sites.

Return type:

DataFrame

tcren.potts.site_energies(frame, potential)

Per-site interaction energy \(\varphi_a\) = potential(aa.rec, aa.par).

A pair the potential does not define — TCRen2 carries no from = C row, so it is 380 cells and not 400 — contributes zero, matching tcren.energetics.scoring._interface_energy()’s nansum. The count of such sites is returned rather than swallowed, because a silent zero on a large fraction of the interface is the failure mode this whole path exists to avoid.

Returns:

(phi, undefined) — the energies, and the per-site boolean mask of the sites the potential does not define, so a caller can count them per structure rather than in bulk.

Parameters:

• frame (DataFrame)

• potential (Potential)

Return type:

tuple[ndarray, int]

tcren.potts.observation_field(frame, *, cutoff=5.0, margin=0.5, clip=40.0)

Per-site log-odds that the crystal’s distance read contributes to \(\eta\).

(cutoff - d_heavy) / margin, clipped. Positive inside the cutoff, negative outside, and steep in proportion to how finely the margin says the distance was resolved.

A null d_heavy is a confident non-contact: tcren.potts.available_pairs() measures the heavy-atom layer only out to probe, so a null means the pair missed by more than that, which is the reading — not missing data. Pass probe = cutoff + 3 * margin and the boundary sits where the field is already below −3.

Parameters:

• frame (DataFrame) – Sites carrying d_heavy (and sigma when margin <= 0).

• cutoff (float) – Contact definition, Å — the model’s own cutoff.

• margin (float) – Resolution of the distance call, Å. 0 is the hard crystal map.

• clip (float) – Saturation of the returned log-odds.

Returns:

One float per row of frame.

Return type:

ndarray

Example

>>> import numpy as np, polars as pl
>>> f = pl.DataFrame({"d_heavy": [4.0, 5.0, 6.0, None]})
>>> np.round(observation_field(f, margin=0.5), 3)
array([  2.,   0.,  -2., -40.])

tcren.potts.contact_posterior(sites, model, *, margin=0.5, clip=40.0, chains=64, burn=100, draws=100, thin=3, seed=0, workers=None)

The E-step: \(P(\sigma_a = 1 \mid \text{the crystal's distances}, \theta)\) per site.

tcren.potts.contact_probabilities() computed against the observation-shifted field, so the couplings still propagate a resampled site to its neighbours while the confidently read sites hold their value. Same sampler, same worker split, same blake2b(seed, pdb.id) streams — so the result is a function of the structure and not of its position in the table.

At margin <= 0 no sampling happens at all: the posterior is a point mass on the observed map, and p_post is sigma exactly. That exactness is what makes the reduction to TCRen2 a test rather than an approximation.

Parameters:

• sites (DataFrame) – Rows from tcren.potts.available_pairs(), ideally with probe = cutoff + 3 * margin.

• model (PottsModel) – The current PottsModel.

• margin (float) – Resolution of the distance call, Å; see observation_field().

• clip (float) – Saturation of the observation log-odds.

• chains (int) – Gibbs settings; see tcren.potts.gibbs().

• burn (int) – Gibbs settings; see tcren.potts.gibbs().

• draws (int) – Gibbs settings; see tcren.potts.gibbs().

• thin (int) – Gibbs settings; see tcren.potts.gibbs().

• seed (int) – Seed for the sampler.

• workers (int | None) – Processes to split structures across. None uses every core.

Returns:

The prepared site frame with eta (post-observation), eta_prior (the model’s own) and p_post appended, one row per site.

Return type:

DataFrame

tcren.potts.weighted_contacts(post, *, weight_col='p_post')

Rename contact_posterior()’s output into the schema derive_tcren() reads.

The two halves of this loop were written years apart and name the same thing differently — aa.rec/aa.par against residue.aa.from/residue.aa.to. This is the whole bridge.

Every available pair is emitted, not only the contacts: a site at p_post near zero contributes near-zero mass and dropping it would be the same tally with more code. That the universe is the 15 Å available set rather than the observed contacts is exactly what lets a near-miss vote.

Parameters:

• post (DataFrame) – Output of contact_posterior().

• weight_col (str) – Which probability column becomes the fractional count.

Returns:

pdb.id, residue.aa.from, residue.aa.to and the weight column.

Return type:

DataFrame

tcren.potts.fit_selfconsistent(sites, *, potential=None, margin=0.5, damping=0.5, steps=15, tol=0.005, pseudocount=1, pin_centred=False, weights=None, ridge=1.0, chains=64, burn=100, draws=100, thin=3, seed=0, workers=None, notes='', **fit_kwargs)

Iterate potential \(\leftrightarrow\) contact model to a fixed point.

Each step recomputes the contacts under the current model (E-step), re-derives the potential from those contacts weighted by their probability, damps it against the previous matrix, and refits the model with the new matrix pinned as its pair term.

Two things are estimated each step, not one. The 20x20 energy matrix is re-derived from the probability-weighted counts; the coupling kernel — the 18 classes over sequence offsets that carry the “crosses”, a made contact pulling its own sequence neighbours onto the same partner residue — is refitted by pseudolikelihood, jointly with every field and with the scale on the matrix. They are not independent: the E-step’s posterior for a borderline site depends on whether its neighbours contacted, so the couplings shape the counts the matrix is derived from. Both land in the trace.

The pair term is pinned, not free. J carries one coefficient — the scale beta_matrix on the current matrix — instead of 400. That is what closes the loop: a free J would absorb the pair structure and leave the derived matrix with no route back into \(P(\sigma)\). It also makes the fixed point well posed and beta_matrix reportable, as the inverse temperature relating TCRen units to contact log-odds.

pin_centred=False by default here, against fit_potts’s own default: double-centring strips the one-body terms, which is right for ranking matrices against each other and wrong for a loop whose output is a matrix meant to be scored raw.

Parameters:

• sites (DataFrame) – Rows from tcren.potts.available_pairs() with probe past the cutoff.

• potential (Potential | None) – Starting matrix. None (default) is TCRen2, the shipped deterministic estimate — the loop’s \(e^{(0)}\).

• margin (float) – Resolution of the distance call, Å. 0 returns after one step having reproduced the starting derivation.

• damping (float) – Mixing weight \(\lambda\) on the newly derived matrix. Fixed-point iteration on a statistic of its own output overshoots undamped.

• steps (int) – Maximum iterations.

• tol (float) – Convergence on max_abs_de, in the potential’s own (natural-log) units. The default sits above the sampler’s noise floor deliberately — see the note below.

• pseudocount (int) – Passed to derive_tcren(); 1 is the shipped value, and dropping it sends 17 of 380 cells to infinity.

• pin_centred (bool) – Double-centre the pinned matrix. See above.

• weights (dict[str, float] | None) – Per-structure redundancy weights, applied to both the derivation and the fit — e.g. tcren.potential.balanced_weights(), which is what ships TCRen2.

• ridge (float) – Ridge on the pseudolikelihood fit.

• chains (int) – Gibbs settings for the E-step.

• burn (int) – Gibbs settings for the E-step.

• draws (int) – Gibbs settings for the E-step.

• thin (int) – Gibbs settings for the E-step.

• seed (int) – Base seed; step t samples at seed + t so the steps are not correlated draws.

• workers (int | None) – Processes for the E-step.

• notes (str) – Free text recorded on the returned model.

• **fit_kwargs – Forwarded to fit_potts() (radius, cutoff, joint, …).

Returns:

(potential, model, trace) — the converged matrix, the model pinned to it, and one trace row per step carrying step, max_abs_de, r_vs_initial, beta_matrix, rms_de, contact_mass, mass_from_uncertain, pseudo_loglik, n_sites, converged, and one column per coupling class of tcren.potts.kernel_names() — the loop re-estimates the matrix and the couplings together, so the trace carries both.

Raises:

ValueError – If steps is below 1, or a derived matrix’s alphabet stops matching.

Return type:

tuple[Potential, PottsModel, DataFrame]

Note

``max_abs_de`` has a Monte-Carlo floor and ``tol`` has to sit above it, the same way tcren.potts.fit_potts_em()’s max_shift does. The E-step is sampled, so once the matrix stops moving systematically the trace oscillates at the sampler’s own noise rather than reaching zero. Measured on Native2026 at --balance both, 359 structures, the default Gibbs settings (chains=64, draws=100, thin=3): max_abs_de falls 0.3499 -> 0.0033 over ten steps and then sits between 0.0033 and 0.0038 for the next five. Hence a default tol of 5e-3. Lower it only together with more draws or chains, which is what lowers the floor; lowering tol alone buys a run that cannot converge.

The floor scales as the sampler says it should, so raising ``draws`` cannot rescue an arbitrarily small ``tol``. Measured on the same set at draws=400, 25 steps: the plateau mean falls from 0.00391 (draws=100, steps 10-15) to 0.00217 (draws=400, steps 13-25), a ratio of 0.553 against the 0.500 that quartering the Monte-Carlo variance predicts. Reaching 1e-3 on max_abs_de would take roughly 1,800 draws, and it would measure nothing: over those same steps r_vs_initial moves by 6.2e-5 and beta_matrix by 2.1e-4, so the matrix has long stopped. rms_de is the statistic to watch when what you want to know is whether the matrix is still moving – averaging over the 380 cells rather than maximising, its floor is lower by about sqrt(380).

tcren.potts.latent_mask(frame, *, cutoff=5.0, margin=0.5)

Boolean mask of the sites whose contact call is inside its own resolution.

|d_heavy - cutoff| <= margin. A pair at 4.9 Å and one at 5.1 Å are recorded as contact and non-contact, and nothing in the coordinates justifies the difference.

The measurement probe must cover cutoff + margin. Tables emitted by available_pairs carry probe and a narrower measurement is rejected. Counts alone cannot diagnose truncation: a complete small table may contain only contacts in the band. Legacy tables without probe provenance cannot be certified here and must be regenerated for publication.

Note

A second class of genuinely-latent site — a residue with no resolved side-chain atoms, so d_heavy is a backbone-only lower bound — is not detected here, because the pair table does not carry atom counts. Such a site reads as a confident non-contact and the fit believes it.

Parameters:

• frame (DataFrame)

• cutoff (float)

• margin (float)

Return type:

ndarray

tcren.potts.impute(sites, model, *, margin=0.5, draws=5, chains=1, burn=50, thin=5, seed=0)

The E-step: draws completions of the contact map, stacked into one frame.

Every site outside the latent set keeps the value the structure gave it; the latent ones are resampled from their conditional under model. Returns (stacked frame, latent mask, mean imputed sigma); the frame’s pdb.id carries a #m draw suffix.

Parameters:

• sites (DataFrame)

• model (PottsModel)

• margin (float)

• draws (int)

• chains (int)

• burn (int)

• thin (int)

• seed (int)

Return type:

tuple[DataFrame, ndarray, float]

tcren.potts.fit_potts_em(sites, *, margin=0.5, draws=5, steps=10, init=None, chains=1, burn=50, thin=5, tol=0.001, seed=0, ridge=1.0, notes='', **fit_kwargs)

Gibbs-imputation/pseudolikelihood hybrid over borderline contact calls.

Parameters:

• sites (DataFrame) – Rows from tcren.potts.available_pairs().

• margin (float) – Half-width of the latent band around cutoff, Å. 0 makes every site observed, in which case this reduces to tcren.potts.fit_potts() exactly.

• draws (int) – Completions per E-step. Each enters the M-step at weight 1/draws.

• steps (int) – Maximum EM iterations.

• init (PottsModel | None) – Starting parameters. None runs the plain pseudolikelihood fit first, which is both the natural warm start and the model the result is compared against.

• chains (int) – Gibbs settings for the conditional sampler.

• burn (int) – Gibbs settings for the conditional sampler.

• thin (int) – Gibbs settings for the conditional sampler.

• tol (float) – Stop when no parameter moves by more than this between iterations.

• seed (int) – Seed for the sampler.

• ridge (float) – Passed to tcren.potts.fit_potts(). A weights entry is per-structure redundancy weighting keyed on the unsuffixed pdb.id; it is multiplied by the 1/draws each completion enters at.

• notes (str) – Passed to tcren.potts.fit_potts(). A weights entry is per-structure redundancy weighting keyed on the unsuffixed pdb.id; it is multiplied by the 1/draws each completion enters at.

• fit_kwargs – Passed to tcren.potts.fit_potts(). A weights entry is per-structure redundancy weighting keyed on the unsuffixed pdb.id; it is multiplied by the 1/draws each completion enters at.

Returns:

step, n_latent, mean_imputed_sigma, pseudo_loglik, max_shift, max_shift_vs_init and converged.

Return type:

(model, trace). The trace has one row per iteration

Warning

Monte-Carlo EM with a fixed draws does not converge to a point — it random-walks in a ball around the optimum whose radius is set by the Monte-Carlo noise, so max_shift has a floor that tol may sit below. Read max_shift_vs_init beside it: a run that is jittering has a flat one, a run that is still moving does not. Raise draws to lower the floor. Everything is seeded, so a run is exactly reproducible regardless.

Note

pseudo_loglik on the returned model is the M-step’s value on completed data and is not comparable to the pseudolikelihood fit’s value on observed data — different response. Compare the two fits through consistency(), not through their objectives.

tcren.potts.consistency(reference, other)

How far a refit moved, per parameter block, in units of the reference’s own scale.

The yardstick is the reference fit’s cluster-robust standard error where it has one — the coupling kernel, the only block tcren.potts.fit_potts() stores kernel_se for. For every other block the yardstick is that block’s own standard deviation in the reference, which answers “did this move by a lot compared with how much it varies at all”.

Returns one row per block: n, max_abs, mean_abs, scale, max_scaled and n_beyond_1 — how many coefficients moved by more than one unit of the yardstick.

Parameters:

• reference (PottsModel)

• other (PottsModel)

Return type:

DataFrame

Constants

The alphabet, the region and role vocabularies, the coupling offsets and the fitting defaults. automodule does not reach re-exported data members, so they are named here explicitly.

tcren.potts.AA = ('A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y')

Built-in immutable sequence.

If no argument is given, the constructor returns an empty tuple. If iterable is specified the tuple is initialized from iterable’s items.

If the argument is a tuple, the return value is the same object.

tcren.potts.REGIONS = ('TRA:CDR1', 'TRA:CDR2', 'TRA:CDR3', 'TRA:FR1', 'TRA:FR2', 'TRA:FR3', 'TRA:FR4', 'TRB:CDR1', 'TRB:CDR2', 'TRB:CDR3', 'TRB:FR1', 'TRB:FR2', 'TRB:FR3', 'TRB:FR4', 'other')

Built-in immutable sequence.

If no argument is given, the constructor returns an empty tuple. If iterable is specified the tuple is initialized from iterable’s items.

If the argument is a tuple, the return value is the same object.

tcren.potts.ROLES = ('tcr_facing', 'anchor', 'HELIX_A1', 'HELIX_A2', 'HELIX_B1', 'GROOVE_FLOOR')

Built-in immutable sequence.

If no argument is given, the constructor returns an empty tuple. If iterable is specified the tuple is initialized from iterable’s items.

If the argument is a tuple, the return value is the same object.

tcren.potts.CLASSES = ('peptide', 'mhc')

Built-in immutable sequence.

If no argument is given, the constructor returns an empty tuple. If iterable is specified the tuple is initialized from iterable’s items.

If the argument is a tuple, the return value is the same object.

tcren.potts.CDR_LOOPS = ('TRA:CDR1', 'TRA:CDR2', 'TRA:CDR3', 'TRB:CDR1', 'TRB:CDR2', 'TRB:CDR3')

Built-in immutable sequence.

If no argument is given, the constructor returns an empty tuple. If iterable is specified the tuple is initialized from iterable’s items.

If the argument is a tuple, the return value is the same object.

tcren.potts.GROOVE_REGIONS = ('HELIX_A1', 'HELIX_A2', 'HELIX_B1', 'GROOVE_FLOOR')

Built-in immutable sequence.

If no argument is given, the constructor returns an empty tuple. If iterable is specified the tuple is initialized from iterable’s items.

If the argument is a tuple, the return value is the same object.

tcren.potts.MHC_PARTNER = ('MHCa', 'MHCb')

Built-in immutable sequence.

If no argument is given, the constructor returns an empty tuple. If iterable is specified the tuple is initialized from iterable’s items.

If the argument is a tuple, the return value is the same object.

tcren.potts.OFFSETS = ((0, 1), (0, 2), (1, -2), (1, -1), (1, 0), (1, 1), (1, 2), (2, -2), (2, -1), (2, 0), (2, 1), (2, 2))

Built-in immutable sequence.

If no argument is given, the constructor returns an empty tuple. If iterable is specified the tuple is initialized from iterable’s items.

If the argument is a tuple, the return value is the same object.

tcren.potts.CROSS_DJ = (0, 1, 2)

Built-in immutable sequence.

If no argument is given, the constructor returns an empty tuple. If iterable is specified the tuple is initialized from iterable’s items.

If the argument is a tuple, the return value is the same object.

tcren.potts.DBIN = 0.5

Convert a string or number to a floating point number, if possible.

tcren.potts.DEFAULT_RIDGE = 1.0

Convert a string or number to a floating point number, if possible.