Contact-map Potts 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)}\]

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.

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

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.

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_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().

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.

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

# the free-energy effect of every substitution at every peptide position
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 closes those pairs onto a coarser grid, which is where they become comparable with an experiment. The residues of one CDR loop are distinct pairs with different marginals, so the number of simultaneous contacts is Poisson-binomially distributed and has no closed form; the event “at least one” does,

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

with \(p_j\) the model marginal of pair \(j\) in the group. That is the quantity a molecular -dynamics trajectory reports as the fraction of frames in which a loop touches a peptide position, so --by loop is directly comparable with a measured contact-frequency map. --by position collapses the loops and reads how engaged each peptide residue is expected to be, before any residue identity is scored. Emitted columns are p_any (the frequency above), p_expected (\(\sum_j p_j\), the expected number of contacts), n_pairs, n_observed and observed. The sum is accumulated in \(\log(1 - p)\), so a twelve-residue loop does not underflow and a saturated pair returns 1 exactly rather than nan.

These are contact frequencies — dimensionless, in \([0, 1]\). They are not free energies and carry no \(k_\mathrm{B}T\), so they belong to the diagnostic and importance side of the model rather than to any energy block; score’s neg_energy and peptide_free_energy()’s log Z0 are the quantities with units.

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.

Substituting a residue: the free energy, not the frequency#

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 interface free energy moves 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. A partner position carrying two different residues has no sequence to substitute into and is rejected rather than averaged.

Unlike map’s frequencies this is an energy: \(\log Z_0\) carries \(k_\mathrm{B}T\).

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}}\]

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

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

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, kernel_se=<factory>, radius=15.0, cutoff=5.0, joint=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='')[source]#

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), or None for a free J.
g_dist (list[float]) – Cα-distance-bin coefficients, bin b covering [b*DBIN, (b+1)*DBIN) Å.
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.
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) – Penalised pseudo-log-likelihood at the 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)
kernel_se (list[float])
radius (float)
cutoff (float)
joint (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#

kernel_se: list[float]#

radius: float = 15.0#

cutoff: float = 5.0#

joint: 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()[source]#

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

Return type:: ndarray

n_parameters()[source]#

Free parameters, counting a fixed-matrix coupling as the single scale it is.

Return type:: int

to_json(path)[source]#

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

Parameters:: path (str | Path)
Return type:: None

classmethod from_json(path)[source]#

Read a model written by to_json().

Parameters:: path (str | Path)
Return type:: PottsModel

classmethod bundled(name='potts_tcr_peptide')[source]#

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)[source]#

-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 (str) – One of tcren2, tcren, mj, mj1996, keskin.
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)[source]#

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

Parameters:: joint (bool)
Return type:: list[str]

tcren.potts.available_pairs(structure, partner='peptide', *, receptor='tcr', radius=15.0, cutoff=5.0)[source]#

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.

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)[source]#

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; 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)

tcren.potts.eta(codes, model)[source]#

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)[source]#

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).

Return type:

list[ndarray]

tcren.potts.neighbour_counts(edge_lists, sigma, n_sites)[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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, regions=None, notes='')[source]#

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 | None) – Fix J to one scale on this bundled potential (tcren2, mj, …) 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.
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)[source]#

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)[source]#

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().
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.

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)[source]#

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)[source]#

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)[source]#

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

Parameters:: eta (ndarray)
Return type:: float

tcren.potts.energy(sigma, eta, A)[source]#

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)[source]#

<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)[source]#

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)[source]#

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

The Legendre partner of log Z(mu), and 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)[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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.contact_map(sites, model, *, by='loop', chains=64, burn=100, draws=100, thin=3, seed=0, workers=None)[source]#

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 residues of a loop are distinct pairs with different probabilities, so the number of simultaneous contacts is Poisson-binomially distributed and has no closed form. The event “at least one” does, and it is what a frequency map measures:

\[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. 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 predicted frequency 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)[source]#

Free energy of the interface with each residue threaded through each partner position.

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 interface free energy moves 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.

Unlike contact_map()’s frequencies this is an energy: \(\log Z_0\) carries \(k_\mathrm{B}T\) and belongs in an energy block.

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 free energy moves 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.

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)[source]#

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)[source]#

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)[source]#

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

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.