tcren.tex

\documentclass[11pt]{article}

\usepackage[margin=1in]{geometry}
\usepackage{amsmath,amsthm}
\usepackage{mathtools}
\usepackage{booktabs}
\usepackage{graphicx}
% Old-style journal look matching the mhcmatch appendix (OldStandard text+math, small-caps
% headings, run-in paragraph leaders). Requires lualatex -- pinned in .latexmkrc / Makefile.
\usepackage{OldStandard}
\usepackage{unicode-math}
\setmathfont{OldStandard-Math.otf}
\usepackage{titlesec}
\usepackage[hidelinks]{hyperref}

\linespread{1.15}
\setlength{\parindent}{2em}
% allow a large top float (the full-width workflow diagram) near where it is referenced
\renewcommand{\topfraction}{0.9}
\renewcommand{\floatpagefraction}{0.8}
\titleformat{\section}{\normalfont\large\scshape\bfseries}{\thesection.}{0.6em}{}
\titleformat{\subsection}{\normalfont\scshape\bfseries}{\thesubsection.}{0.5em}{}
\titleformat{\paragraph}[runin]{\bfseries\scshape}{}{0pt}{}[.]
\titlespacing*{\paragraph}{\parindent}{0.4em}{0.6em}

\newtheoremstyle{oldplain}{}{}{\itshape}{}{\scshape}{.}{.5em}{}
\newtheoremstyle{oldrem}{}{}{\normalfont}{}{\scshape}{.}{.5em}{}
\theoremstyle{oldplain}
\newtheorem{proposition}{Proposition}
\theoremstyle{oldrem}
\newtheorem{definition}{Definition}
\newtheorem{remark}{Remark}

\DeclareMathOperator*{\argmin}{arg\,min}
\newcommand{\R}{\mathbb{R}}
\newcommand{\Aset}{\mathcal{A}}
\newcommand{\Cset}{\mathcal{C}}
\newcommand{\xref}{x^{0}}

\title{\scshape tcren2 --- a methods reference:\\ contact-potential scoring, ranking, and structure operations}
\author{Mikhail Shugay\textsuperscript{1}\thanks{Correspondence:
  \href{mailto:mikhail.shugay@gmail.com}{mikhail.shugay@gmail.com}}}
\date{%
  \small
  \textsuperscript{1}ISALGO laboratory\\[6pt]
  \normalsize\textit{tcren technical appendix}%
}

\begin{document}
\maketitle

\begin{abstract}
\texttt{tcren} scores a candidate epitope for a fixed TCR--peptide--MHC structure by summing a
residue-level pairwise statistical potential over the observed contact map~\cite{Karnaukhov2024}. The
score depends only on residue identity and the closest-atom contact geometry, so a candidate is scored
\emph{virtually}: its amino acids are threaded onto the native contacts with no atom moved. This
appendix is the implementation-level companion to the method's biophysical synopsis --- it states, as
the code computes them, the contact map, the statistical potentials and their derivation, the per-task
readouts (epitope ranking, percentile rank, mutational $\Delta\Delta G$), the two coordinate-level
operations (backbone-preserving substitution and knowledge-based pose refinement), and the interface
mechanics. Every default parameter is collected in Table~\ref{tab:params}; the formulas are the exact
quantities evaluated.
\end{abstract}

\section{Workflow}\label{sec:workflow}
Figure~\ref{fig:workflow} is the end-to-end pipeline. A single call (\texttt{tcren pipeline}) takes a
structure through: \textbf{(1)~import} (C-gene-trimmed parse of PDB / mmCIF); \textbf{(2)~annotate} ---
TCR loci and CDR/framework regions by germline alignment (\textsc{arda}), MHC allele, class, chain role
and groove positions; \textbf{(3)~superimpose} (optional) onto the canonical database, placing the
complex in a common MHC-groove $\mathrm C_\alpha$ frame so downstream geometry is comparable across
complexes; \textbf{(4)~contacts} --- the $5\,\text{\AA}$ heavy-atom contact map with a per-residue
region markup; and \textbf{(5)~score} --- each interface summed under its potential. Candidate epitopes,
percentile ranks and $\Delta\Delta G$ are read \emph{virtually} off the resulting map; substitution,
refinement and mechanics operate on coordinates. The same operations are exposed as subcommands
(\texttt{annotate}, \texttt{orient}, \texttt{contacts}, \texttt{score}, \texttt{rank}, \texttt{ddg},
\texttt{refine}) and as a Python API.

\begin{figure}[tbp]\centering
\includegraphics[width=0.9\linewidth]{workflow.pdf}
\caption{The \texttt{tcren2} workflow. The single-call pipeline parses, annotates, orients and builds
the $5\,\text{\AA}$ contact map, then scores each interface (TCRen on TCR--peptide, Miyazawa--Jernigan
on the presentation interfaces). Candidate peptides thread onto the native contacts for virtual scoring,
percentile ranking and $\Delta\Delta G$ (no atom moved); the coordinate operations (substitute, DOPE
refine, external FlexPepDock relaxation) and the interface-mechanics off-rate proxy consume the oriented
structure. Node annotations are the code defaults of Table~\ref{tab:params}.}
\label{fig:workflow}
\end{figure}

\section{Contact map}\label{sec:contacts}
For a chain-typed complex \texttt{tcren} builds the inter-chain contact map: two residues on different
chains are in contact when their closest heavy-atom pair lies within the cutoff $r_c=5\,\text{\AA}$,
reduced to \emph{one} pair per residue pair (the closest). Partitioning by chain type gives the three
interfaces of the ternary complex,
\[
  \Cset = \Cset_{\mathrm{TCR:pep}} \,\cup\, \Cset_{\mathrm{TCR:MHC}} \,\cup\, \Cset_{\mathrm{pep:MHC}} .
\]
Each contact carries the two residue identities and, optionally, its heavy-atom-pair count
$n^{\mathrm{atoms}}_{ij}$ (built with \texttt{count\_atoms=True}). The map is the single object every
downstream step consumes; the contact energy is orientation-invariant, while the optional canonical
frame (step~3) is what makes \emph{geometry} comparable across complexes.

\section{Statistical potentials}\label{sec:potentials}
\paragraph{TCRen} The TCRen potential is the Boltzmann inversion of observed TCR--peptide contact
frequencies against an independent-pairing reference~\cite{Karnaukhov2024,Sippl1990}. From the contact
table, with $n_{ab}$ the count of directed contacts between partner residue $a$ and peptide residue $b$
(each incremented by a pseudocount $\kappa$), row/column/grand totals $N_a=\sum_b n_{ab}$,
$N_b=\sum_a n_{ab}$, $N=\sum_{ab}n_{ab}$,
\begin{equation}\label{eq:tcren}
  \mathrm{odds}_{ab}=\frac{n_{ab}\,N}{N_a\,N_b},
  \qquad
  \phi^{\mathrm{TCRen}}_{ab}=-\ln\mathrm{odds}_{ab}
\end{equation}
(lower is more favourable). The grid enumerates the twenty amino acids; cysteine is dropped from the
partner (``from'') axis \emph{after} the log-odds, a data convention of the reference derivation. An
alignment-matrix variant replaces $-\ln$ by $-\log_2(\cdot)/\beta$ over 21 symbols (a gap included),
$\beta=44$; the classic natural-log form above is the shipped default and reproduces
\texttt{TCRen\_potential.csv} exactly. Optional per-structure weights down-weight redundant complexes in
the counts while keeping all data.

\paragraph{Presentation interfaces} The peptide--MHC and TCR--MHC interfaces are scored with the
Miyazawa--Jernigan quasi-chemical contact potential~\cite{Miyazawa1996} (Keskin parameterisation), the
same log-odds construction \eqref{eq:tcren} applied to a generic protein-contact reference. Using an
independent potential per interface --- TCRen where thymic selection shaped the statistics, MJ where it
did not --- is the combination the benchmarks favour.

\section{Scoring, ranking, and \texorpdfstring{$\Delta\Delta G$}{ddG}}\label{sec:scoring}
\paragraph{Interface energy} A candidate peptide $s$ is scored on an interface by threading its residues
onto the contacted peptide positions and summing the potential,
\begin{equation}\label{eq:score}
  \Phi_{XY}(s)=\!\!\sum_{(i,j)\in\Cset_{XY}}\!\! w_{ij}\,\phi\bigl(a_i,\,s_{\pi(j)}\bigr),
\end{equation}
where $\phi$ is the interface's potential, $a_i$ the fixed partner residue, $\pi(j)$ the contact's
peptide position, and the weight is $w_{ij}=1$ (default, \emph{residue} weighting) or
$w_{ij}=n^{\mathrm{atoms}}_{ij}$ (\emph{atomic} weighting, which tracks the LJ$+$Coulomb atom-pair sum
more closely). Pairs absent from the potential (e.g.\ a cysteine partner) contribute nothing. The total
complex energy is
$\Phi=\Phi_{\mathrm{TCR:pep}}(\text{TCRen})+\Phi_{\mathrm{TCR:MHC}}(\text{MJ})+\Phi_{\mathrm{pep:MHC}}(\text{MJ})$.
Scoring is fixed-backbone and fixed-contact-set --- the native pose stands in for every candidate --- so
a whole library ranks from one structure in microseconds per peptide.

\paragraph{Percentile rank} To turn a raw energy into a calibrated rank, the native peptide is scored
against a background of $N_{\mathrm{bg}}=1000$ peptides of the same length --- uniform-random over the
twenty amino acids, or position-shuffled draws from a supplied epitope list --- and its percentile
reported,
\begin{equation}\label{eq:rank}
  \mathrm{rank\%}=\frac{1}{N_{\mathrm{bg}}}\bigl|\{\,b:\Phi(b)\le\Phi(\text{native})\,\}\bigr| .
\end{equation}
A small \texttt{rank\_pct} marks a strong binder relative to chance; the background seed is fixed
(default $0$) for reproducibility.

\paragraph{Mutational $\Delta\Delta G$} A peptide point mutation is scored \emph{virtually} --- no atom
moves, no re-docking --- by re-scoring the mutant on the same contact map,
\begin{equation}\label{eq:ddg}
  \Delta\Delta G=\Phi(\text{native})-\Phi(\text{mutant}),
\end{equation}
so a positive value flags a destabilising mutation. Only contacts touching the mutated position survive
the difference; an interface without the peptide (TCR--MHC) contributes identically zero. This is the
fast triage matrix (\texttt{alanine\_scan}, \texttt{neoantigen\_ddg}); a quantitative,
sampling-dependent $\Delta\Delta G$ is left to the coordinate path below and, ultimately, to full
relaxation.

\paragraph{Poly-alanine reference} Summing \eqref{eq:ddg} over every peptide position gives the
poly-alanine reference difference $\Delta\Phi=\Phi(s)-\Phi(\mathrm A^{L})$ (\texttt{reference\_delta}) ---
the full-peptide alanine scan, which subtracts the pose's identity-independent baseline $\Phi(\mathrm
A^{L})$ (what the geometry scores when every peptide residue is alanine) and keeps the sequence-specific
part. On a \emph{fixed} contact map $\Delta\Phi=\Phi(s)-\text{const}$, so it re-ranks nothing; it matters
only across candidates that each carry their \emph{own} structure (e.g.\ AlphaFold peptide-swap models),
where it normalises out per-pose interface geometry and rescues forced or mis-registered poses whose
geometry corrupts the raw energy. Both $\Phi$ and $\Delta\Phi$ are \emph{relative, within-receptor}
peptide rankings, not binding constants: they order a peptide panel against one fixed TCR:pMHC (e.g.\ the
graded B\texttt{*}27:05 EC$_{50}$ series, Spearman $\rho\approx0.75$--$0.9$) but predict equilibrium
affinity \emph{across} complexes only weakly --- on the ATLAS SPR set raw $\Phi$ and $\Delta\Phi$ track
$\Delta G$, $K_d$, $k_{\mathrm{off}}$ and $k_{\mathrm{on}}$ at $|\rho|\le0.3$; the off-rate, when a
structure predicts anything, comes from the interface mechanics of \S\ref{sec:mechanics}, not the
contact sum.

\section{Backbone-preserving peptide substitution}\label{sec:sub}
\begin{definition}[Substitution]\label{def:sub}
Let the native peptide chain have residues $r_1,\dots,r_L$ with heavy-atom coordinate sets
$X_1,\dots,X_L$, and let $s\in\Aset^L$ be a new sequence over the twenty amino acids. The
substitution $\mathrm{sub}(\,\cdot\,,s)$ returns a chain whose $i$-th residue has identity $s_i$ and
retains exactly the backbone atoms $\{\mathrm N,\mathrm C_\alpha,\mathrm C,\mathrm O\}$ together with
$\mathrm C_\beta$ when $s_i\neq\mathrm{Gly}$, with all coordinates inherited from $X_i$; side-chain
atoms beyond $\mathrm C_\beta$ are discarded.
\end{definition}

\paragraph{Rationale} A residue-level potential reads the contact map through residue
\emph{identity} and the closest-atom (or $\mathrm C_\beta$) distance~\cite{Miyazawa1996}; the
backbone $+\,\mathrm C_\beta$ frame therefore carries everything the score and the refinement energy
of \S\ref{sec:refine} consume, while the native backbone---rigidly constrained between the MHC
groove walls and the TCR loops---is the appropriate scaffold for an equal-length epitope. Rebuilding
the discarded side chains (a backbone-dependent rotamer search~\cite{Shapovalov2011}) and full
physics relaxation are deliberately \emph{not} performed here; substitution only fixes identity and
$\mathrm C_\beta$ geometry, leaving heavier modelling to the refiner or to an external tool. The map
is purely combinatorial---no alignment, no optimisation---so it is exact and order-preserving on the
peptide's sequential indexing.

\section{Potential-guided refinement}\label{sec:refine}
Write $x\in\R^{3M}$ for the $M$ peptide heavy-atom coordinates (the only degrees of freedom) and
$y\in\R^{3M'}$ for the fixed partner (TCR\,$+$\,MHC) heavy atoms. Each heavy atom carries a
\textsc{dope} atom class $c(\cdot)\in\{1,\dots,158\}$---a MODELLER mean-force atom type, i.e.\ a
(residue, atom-name) label---and let $\varphi_{c,c'}:\R_{\ge0}\to\R$ be the \textsc{dope}
distance-dependent potential of mean force for the class pair $(c,c')$~\cite{Shen2006}, tabulated on
$29$ knots at $0.75$--$14.75\,\text{\AA}$ ($0.5\,\text{\AA}$ spacing), linearly interpolated, and set
to $0$ beyond $d_{\max}=14.75\,\text{\AA}$. Crucially $\varphi$ is \emph{not} TCRen or
Miyazawa--Jernigan---those score the final epitope (\S\ref{sec:scoring}), so reusing them to refine
would optimise the pose against the very quantity reported. \textsc{dope} is a single, chain-agnostic,
atom-resolved potential applied uniformly to every peptide$\leftrightarrow$partner pair.
\begin{definition}[Refinement energy]\label{def:energy}
\[
  E(x) \;=\;
  \underbrace{\sum_{p=1}^{M}\sum_{q=1}^{M'} \varphi_{c(p),\,c(q)}\!\bigl(\lVert x_p - y_q\rVert\bigr)}_{\text{\textsc{dope}}}
  \;+\;
  w_{\mathrm r}\sum_{p=1}^{M}\lVert x_p-\xref_p\rVert^2 .
\]
\textsc{dope}'s short-range knots are strongly repulsive (capped near $+10$), so the potential
supplies its own excluded-volume term: there is no separate clash penalty. The second term is a
harmonic restraint to the input pose $\xref$ (default $w_{\mathrm r}=\tfrac12$).
\end{definition}

\paragraph{Boundedness and the role of the restraint} Unlike a residue contact-count energy---whose
trivial rigid-body minimum can be ``no contacts'' (the peptide ejected)---the \textsc{dope} term is
bounded below and tends to $0$ as the peptide is translated away, so a finite rigid-body minimiser
always exists.
\begin{proposition}\label{prop:restraint}
Restrict $x$ to the rigid-body family $x=R\xref+t$, $R\in SO(3)$, $t\in\R^3$, and let $E_{\mathrm
D}(x)=\sum_{p,q}\varphi_{c(p),c(q)}(\lVert x_p-y_q\rVert)$. Then $E_{\mathrm D}$ is bounded
($|E_{\mathrm D}|\le M M'\max|\varphi|$) and $E_{\mathrm D}\to 0$ as $\lVert t\rVert\to\infty$. Its
rigid-body minimiser is the globally best-packed pose, which need not coincide with the input
$\xref$: the peptide may slide within the groove. Adding $w_{\mathrm r}\sum_p\lVert
x_p-\xref_p\rVert^2$ makes $E$ coercive and confines the minimiser to a neighbourhood of $\xref$
whose radius decreases with $w_{\mathrm r}$, yielding a \emph{local} refinement---a native pose
(already near a \textsc{dope} optimum) barely moves, while a clashed or displaced input relaxes
locally.
\end{proposition}

\paragraph{Sampling} $E$ is minimised by simulated-annealing Metropolis Monte
Carlo~\cite{Metropolis1953,Kirkpatrick1983}. At step $t$ a trial pose is proposed by a rigid move of
the peptide: a Gaussian translation $\delta\sim\mathcal N(0,\sigma_t^2 I_3)$ and a rotation about the
peptide centroid by a random axis and angle $\theta\sim\mathcal N(0,\sigma_r^2)$. The trial is
accepted with probability $\min\!\bigl(1,\exp(-\Delta E/T_t)\bigr)$, with the temperature cooled
linearly $T_t:T_0\!\downarrow\!T_1$; the lowest-energy pose visited is returned. The chain is
deterministic given its pseudo-random seed.

\paragraph{Implementation and cost} Each energy evaluation is $O(MM')$ over heavy-atom pairs;
partner atoms are pre-filtered to those within a $12\,\text{\AA}$ shell of the peptide, so $M'$ is
the interface, not the whole complex, and several thousand steps complete in well under a second.
The kernel is C++ (\texttt{pybind11}); the symmetric $158\times158\times29$ \textsc{dope} table is
passed by buffer, an atom's class is its (residue, atom) label, and $\varphi$ is a single buffer
lookup plus a linear interpolation. The bundled table is parsed once from the MODELLER mean-force
libraries (\texttt{atmcls-mf.lib}, \texttt{dist-mf.lib}) as redistributed by the \texttt{pymod/altmod}
project.\footnote{\url{https://github.com/pymodproject/altmod} (\textsc{dope} libraries from MODELLER;
\cite{Shen2006}).} A native pose is a fixed point of the move set up to thermal noise (it barely
moves); a clashed or displaced input relaxes locally.

\section{Interface mechanics}\label{sec:mechanics}
The affinity-adjacent quantity a native structure actually predicts is not the equilibrium free energy
but the dissociation \emph{rate} $k_{\mathrm{off}}$, which sets the bond lifetime a T~cell reads through
kinetic proofreading. \texttt{tcren} models it mechanically (\texttt{tcren.mechanics}): heavy-atom
interface contacts within $8\,\text{\AA}$ become elastic-network springs with inverse-square-distance
stiffness, and the interface is pulled along the TCR$\leftrightarrow$pMHC dissociation axis. The tensile
stiffness $K_{\mathrm{tens}}$ (\texttt{stiffness\_tensor}) tracks experimental $\log k_{\mathrm{off}}$
--- a stiffer interface ruptures more slowly (Bell--Evans) --- while a breakable-spring variant records
the peak resisting force and work past a strain threshold (\texttt{rupture}). These consume the oriented
structure only; no potential or sequence enters.

\section{Parameter reference}\label{sec:params}
Table~\ref{tab:params} collects every default. Values are the code defaults (function signatures / CLI
options); the last column points to the formula or section that uses each.

\begin{table}[ht]\centering\small
\caption{Default parameters of the \texttt{tcren2} computations. One value per cell; parenthesised names
are the CLI options.}
\label{tab:params}
\begin{tabular}{@{}lll@{}}
\toprule
Quantity & Default & Where \\
\midrule
Contact cutoff $r_c$ & $5\,\text{\AA}$ (\texttt{-{}-cutoff}) & \S\ref{sec:contacts} \\
Contact reduction & closest heavy-atom pair, one per residue pair & \S\ref{sec:contacts} \\
Contact weight & \texttt{residue} ($w{=}1$), or \texttt{atomic} ($w{=}n^{\mathrm{atoms}}$) & \eqref{eq:score} \\
TCRen potential & $-\ln\mathrm{odds}$, 20 aa (classic) & \eqref{eq:tcren} \\
Pseudocount $\kappa$ & $1$ (\texttt{-{}-pseudocount}) & \eqref{eq:tcren} \\
Cysteine & dropped from ``from'' axis (classic) & \eqref{eq:tcren} \\
am-variant temperature $\beta$ & $44$ & \S\ref{sec:potentials} \\
Presentation potential & Miyazawa--Jernigan (Keskin) & \S\ref{sec:potentials} \\
Background size $N_{\mathrm{bg}}$ & $1000$ (\texttt{-{}-background}) & \eqref{eq:rank} \\
Background seed & $0$ & \eqref{eq:rank} \\
$\Delta\Delta G$ & $\Phi(\text{nat})-\Phi(\text{mut})$, peptide interfaces & \eqref{eq:ddg} \\
Substitution frame & backbone $\{\mathrm N,\mathrm C_\alpha,\mathrm C,\mathrm O\}+\mathrm C_\beta$ & Def.~\ref{def:sub} \\
\textsc{dope} atom classes & $158$ MODELLER mean-force types & \S\ref{sec:refine} \\
\textsc{dope} knots & $29$ at $0.75$--$14.75\,\text{\AA}$, $0.5\,\text{\AA}$ spacing & \S\ref{sec:refine} \\
\textsc{dope} cutoff $d_{\max}$ & $14.75\,\text{\AA}$ ($\varphi{=}0$ beyond) & \S\ref{sec:refine} \\
Restraint weight $w_{\mathrm r}$ & $0.5$ (\texttt{-{}-restraint}) & Def.~\ref{def:energy} \\
Monte-Carlo steps & $2000$ (\texttt{-{}-steps}) & \S\ref{sec:refine} \\
Partner shell & $12\,\text{\AA}$ around the peptide & \S\ref{sec:refine} \\
Mechanics cutoff & $8\,\text{\AA}$ & \S\ref{sec:mechanics} \\
Spring stiffness & inverse-distance$^2$ & \S\ref{sec:mechanics} \\
Pull direction & tensile (dissociation axis) & \S\ref{sec:mechanics} \\
\bottomrule
\end{tabular}
\end{table}

\section{Assumptions and limitations}\label{sec:limits}
The peptide is treated as a \emph{rigid} body: there is no backbone flexibility and no side-chain
repacking, so substitution leaves heavy atoms beyond $\mathrm C_\beta$ unmodelled and the refiner
cannot resolve rotameric clashes (a backbone-dependent rotamer pass~\cite{Shapovalov2011} is the
natural next step; until then a substituted side chain contributes no \textsc{dope} terms beyond
$\mathrm C_\beta$). \textsc{dope} is an isotropic, distance-only potential---an orientation-dependent
potential would resolve packing better at the cost of per-atom frames---and the restraint weight
$w_{\mathrm r}$ is a tunable that sets the locality radius. For physically realistic, full-atom
relaxation---side-chain optimisation and backbone minimisation under an all-atom force
field---\texttt{tcren} defers to Rosetta FlexPepDock~\cite{Raveh2010,London2011} with the ref2015
energy function~\cite{Alford2017}, invoked as an external process; the C++ kernel here is a fast
knowledge-based pre-refinement, not a substitute for it.

\bibliographystyle{plain}
\bibliography{refs}

\end{document}