[1]Hui Liu

1]College of Computer and Information Engineering, Nanjing Tech University,
Nanjing, 211800 , Jiangsu, China

# A unified cross-attention model for predicting antigen binding specificity
to both HLA and TCR molecules

Chenpeng Yu  Xing Fang<EMAIL_ADDRESS>[

###### Abstract

The immune checkpoint inhibitors have demonstrated promising clinical efficacy
across various tumor types, yet the percentage of patients who benefit from
them remains low. The binding affinity between antigens and HLA-I/TCR
molecules plays a critical role in antigen presentation and T-cell activation.
Some computational methods have been developed to predict antigen-HLA or
antigen-TCR binding specificity, but they focus solely on one task at a time.
In this paper, we propose UnifyImmun, a unified cross-attention transformer
model designed to simultaneously predicts the binding of antigens to both HLA
and TCR molecules, thereby providing more comprehensive evaluation of antigen
immunogenicity. We devise a two-phase progressive training strategy that
enables these two tasks to mutually reinforce each other, by compelling the
encoders to extract more expressive features. To further enhance the model
generalizability, we incorporate virtual adversarial training. Compared to
over ten existing methods for predicting antigen-HLA and antigen-TCR binding,
our method demonstrates better performance in both tasks. Notably, on a large-
scale COVID-19 antigen-TCR binding test set, our method improves performance
by at least 9% compared to the current state-of-the-art methods. The
validation experiments on three clinical cohorts confirm that our approach
effectively predicts immunotherapy response and clinical outcomes.
Furthermore, the cross-attention scores reveal the amino acids sites critical
for antigen binding to receptors. In essence, our approach marks a significant
step towards comprehensive evaluation of antigen immunogenicity.

###### keywords:

Cross-attention mechanism, neoantigen, T-cell receptor, Human leukocyte
antigen, Virtual adversarial training, Integrated gradient

##  1 Introduction

The anticancer immune response involves a sequence of intricate biological
events that lead to effective kill of cancer cells. Initially, tumor antigens
are released by cancer cells through specific mechanism, and are captured and
processed by antigen-presenting cells (APCs) [1]. These APCs present the
antigens on their outer surface (antigen presentation). When naive T cells
recognize and bind to the presented antigens, they become activated (T-cell
activation) and subsequently differentiate into effector T-cells, such as
cytotoxic T lymphocytes (CTLs) [2][3][4][5][6]. The effector T-cells migrate
to the tumor site and attack cancer cells [7], ultimately inducing their death
[8]. Additionally, some T-cells differentiate into memory T-cells, which
persist in the body for extended periods and provide long-term immunological
memory [9][10]. These steps are referred to as Cancer-Immunity Cycle [11],
which indeed highlight a delicate balance between the recognition of non-self
antigens and the prevention of autoimmunity. Within this cycle, antigen
presentation and T cell activation stand out as two steps critical to the
success of the anticancer immune response [12][13].

The binding of peptides to class I human leukocyte antigen (HLA) molecules is
a fundamental process in tumor antigen presentation [12]. HLA alleles are
known for their high specificity and polymorphism in the human population
[14], allowing them to selectively bind to a narrow range of peptides, thereby
forming the peptide-HLA complex (pHLA) [15]. Identifying the peptides that are
chosen for binding by specific HLA alleles is a critical step in epitope
selection. Subsequently, the recognition of the presented antigens by T-cell
receptors (TCR) is the prerequisite to activate T cells and trigger a robust
immune response [15]. This process is also highly selective, and only a small
subset of peptides can be recognized by TCRs. This selectivity, known as TCR
binding specificity, arise from the high diversity of TCR repertoire
(estimated to range from 1015 to 1061 possible receptors in humans) [16]. This
diversity is primarily manifested in complementarity determining region 3
(CDR3) [17], which directly interacts with pHLA complex and determines the TCR
binding specificity [18][19]. These binding specificity ensures that only
relevant antigens are presented and activate immune system, maintaining the
delicate balance between effective immune responses and autoimmune reactions.

The HLA polymorphism and TCR repertoire diversity represent evolutionary
traits that enable the human immune system to respond to a wide array of
pathogens at individual level [20][21]. Some experimental assays like mass
spectrometry (MS)-eluted HLA ligands [22] and techniques such as single-cell
TCR sequencing [23] and T-scan [24] have been developed to detect antigen-HLA
and antigen-TCR bindings, respectively. However, these experimental assays are
often time-consuming, technically complex, and costly. To address these
challenges, some computational methods have emerged as viable alternatives for
predicting peptide binders [25]. These methods include TransPHLA[26],
MHCflurry[27], netMHCpan4.0[28], DeepLigand[29], BERTMHC[30] for peptide-MHC
binding prediction, as well as PanPep[31], pMTnet[32], DLpTCR[33], ERGO2[34],
TITAN[35] and ATMTCR [36] for peptide-TCR binding prediction. Although some
methods have achieved impressive prediction accuracy, they all consider the
antigen-HLA and antigen-TCR binding prediction as separate tasks. Existing
methods overlook a crucial fact: the immunogenicity of antigens is profoundly
influenced by their binding affinity to both HLA and TCR molecules, rather
than solely relying on one type of binding.

Distinct from previous studies that considered HLA or TCR binding specificity
alone, we propose a unified model UnifyImmun, which integrates the predictive
tasks of HLA-antigen and TCR-antigen binding to establish a one-stop deep
learning framework for comprehensive evaluation of antigen immunogenicity. As
shown in Figure 1, UnifyImmun comprises three self-attention-based encoders
that receive HLA, antigen, and TCR sequences as inputs, respectively. While
these encoders share a common network structure, they operate with distinct
parameters. They independently map the input sequences into embedding vectors
in latent space. Subsequently, to effectively fuse the information of peptides
and receptors, two cross-attention modules are employed. These modules receive
as inputs the latent representations of HLA-antigen and TCR-antigen pairs,
respectively. The outputs of cross-attention modules are then passed through
fully-connected layers and softmax transformations to generate predictions for
HLA-antigen and TCR-antigen binding, respectively. Once training finished, our
model can be independently applied to three prediction tasks: HLA-antigen
binding, or TCR-antigen binding, or the binding of HLA- antigen-TCR. In
particular, the cross-attention mechanism is instrumental in integrating
receptor-ligand information, thereby facilitating the capture of critical
information for antigen binding to HLA or TCR. Additionally, attention scores
generated by our model offer valuable insights into the key positions and
amino acids within the antigen sequence that are pivotal for effective antigen
presentation and subsequent T cell activation. This aspect of our model
contributes significantly to its interpretability, enhancing our understanding
of the underlying immunological processes.

Given the vast diversity of HLA and TCR repertoires, the sequences currently
available for model training are limited and even biased, posing a significant
challenge of overfitting in the development of prediction models. To overcome
this limitation, we introduced virtual adversarial training as a means to
enhance the model generalizability. Specifically, we apply adversarial
perturbations to the embeddings generated by three sequence encoders,
generating virtual samples that aim to maximize the loss function. By
enforcing the model to resist the adversarial perturbations, thereby its
generalizability can be significantly improved. This adversarial training
approach has demonstrated effective enhancement of the model performance (see
Section 4.8).

Ideally, our model prefers to be trained using HLA-antigen-TCR triplet
samples. However, the availability of such triplets is currently limited,
whereas HLA-antigen and TCR-antigen binding data are relatively abundant. To
efficiently leverage this available data, we propose a two-stage progressive
training strategy (see Method for details). Through performance evaluation on
multiple test datasets, we have demonstrated that the two-stage progressive
training effectively enhances the feature extraction capabilities of the
encoders, thereby improving the predictive performance of HLA-antigen and TCR-
antigen binding specificity.

Our work achieves a significant advancement in the field, offering at least
three notable contributions compared to previous studies:

  * •

By integrating both prediction tasks within a unified model, our approach
enables simultaneous evaluation of antigen presentation efficiency and T-cell
activation potential. This comprehensive assessment provides a more holistic
view of antigen immunogenicity, offering deeper insights into the overall
quality and potential efficacy of antigens in triggering immune responses.

  * •

We employ a unique antigen encoder to derive representations from antigen
sequences, which are then utilized for both HLA-antigen and TCR-antigen
binding predictions. This multitask learning stands in contrast to previous
models that focused exclusively on either one alone. By sharing antigen
presentation across tasks, our model significantly enhances generalization
capabilities by leveraging commonalities in the underlying biology processes
of antigen recognition.

  * •

Our unified framework not only achieves superior performance for HLA-antigen
binding predictions, but also boosts the predictive accuracy of TCR-antigen
binding specificity. This multifaceted performance enhancement represents a
substantial advantage over previous methodologies, highlighting the
superiority of our approach in capturing the intricate relationships between
HLA, TCR, and antigen sequences.

Figure 1: Illustrative diagram of UnifyImmun framework and two-phase training
strategy, as well as frequency distribution w.r.t sequence lengths. (a)
Architecture of UnifyImmun based on cross-attention mechanism. (b-c) Frequency
of antigen sequences and TCR CDR3 sequences included in positive training
samples with respect to lengths.

##  2 Results

###  2.1 Performance evaluation on HLA-antigen binding prediction

To evaluate the proposed model’s predictive performance in HLA-antigen
binding, we conducted performance comparison on a hold-out independent test,
an external test, HPV and neoantigen validation datasets. We compared
UnifyImmun against twelve established methods, including TransPHLA [26],
NetMHCpan_EL [28], NetMHCpan_BA [28], ANN [37], PickPocket [38], SMMPMBEC
[39], SMM [40],NetMHCcons [41], NetMHCstabpan [42] and Consensus [43] , and
two attention-based methods [26] published recently (ACME [44] and
DeepAttentionPan [45]). which can be obtained from http://tools.iedb.org/
main/tools-api/. These comparative methods were downloaded as executable
packages from the IEDB database [46], and run on the same test sets using
their default parameters. The performance metrics including AUROC, accuracy,
MCC, and F1-score. Additionally, we provided additional metrics in the
supplementary file, including precision, recall, AUPR, and specificity, to
offer a more comprehensive evaluation of the model’s performance.

We firstly evaluated the performance of UnifyImmun against other comparative
algorithms on a 10% hold-out test set. As shown in Figure 2(a), UnifyImmun
remarkably outperformed all other methods across all evaluation metrics. In
particular, compared to the second-best method, TransPHLA, UnifyImmun achieved
at least a 5% improvement in both AUROC and AUPR. To provide a visual
representation of the performance differences, Figure 2(b) presented the ROC
curves and precision-recall curves for all comparative methods on the hold-out
test set. These curves further validated the superior performance of our
proposed method. Figure 2(c) visualized the UMAP feature of peptide-HLA pairs,
implying the positive and negative samples separated remarkably in the latent
space. To further assess the model’s ability to prioritize HLA-antigen pairs,
we presented the positive predictive value (PPV) for the top 100, top 1000,
and top 5000 predicted positive samples. As illustrated in Figure 2(d),
UnifyImmun achieved an impressive 100% PPV for the top 100 predictions and
maintained excellent performance above 97% for both the top 1000 and top 5000
predictions. In contrast, the other methods did not demonstrate comparable
ability in prioritizing HLA-antigen pairs.

Figure 2: Performance evaluation on predicting peptide-HLA binding
specificity. (a) Performance comparison to twelve existing methods on
independent and external test dataset, respectively. (b) ROC curves and AUC
values achieved by UnifyImmun and eight comparative methods on independent
test set. (c) UMAP feature visualization of peptide-HLA pairs. (d) Positive
predictive value (PPV) for the top 100, top 1000, and top 5000 predicted
positive samples.

For objective performance evaluation, we conducted tests on an external
dataset provided by TransPHLA. As shown in Figure 2(b), while TransPHLA
exhibited performance advantages on its own dataset, UnifyImmun achieves
nearly identical performance to TransPHLA. Furthermore, UnifyImmun notably
outperformed all other methods across all performance metrics, except for
TransPHLA.

The HPV antigens came from a previous study [47] that identified 278
experimentally verified pHLA bindings derived from the HPV16 proteins E6 and
E7, consisting of peptides ranging from 8 to 11 amino acids in length [48]
[49]. UnifyImmun achieved an impressive accuracy rate of 83.8%. This
significantly surpassed the performance of the current state-of-the-art model
TransPHLA, which exhibited a screening rate of only 68%. For neoantigen
validation, we used the dataset released by TransPHLA. These neoantigens were
collected from non-small-cell lung cancer, melanoma, ovarian cancer and
pancreatic cancer in recent studies . This dataset includes 221 experimentally
verified pHLA bindings. UnifyImmun achieved 94.1% accuracy (208 out of 221),
which is comparable to TransPHLA with 96.4%. Collectively, the performance
comparison experiments on four distinct datasets clearly demonstrated the
superior generalization capability of UnifyImmun compared to the current
state-of-the-art methods.

###  2.2 UnifyImmun boosts the predictive performance of TCR-antigen binding
specificity

Our model can actually predict TCR-antigen binding specificity. To assess its
performance, we conducted a comparative analysis with four current state-of-
the-art methods in the field, including Pan-Peptide, ERGO2, pMTnet and DLpTCR.
For a fair assessment, we executed the executable codes of Pan-Peptide, ERGO2,
and pMTnet using recommended parameters on the same workstation as UnifyImmun.
For DLpTCR, we accessed its web server for evaluation. We initially evaluated
the performance of all methods on a 10% hold-out test set. As shown in Figure
3(a), UnifyImmun remarkably outperformed all other methods on this hold-out
test set. Specifically, UnifyImmun achieved AUROC and AUPR values of 0.938 and
0.936, respectively, highlighting its exceptional predictive ability in TCR-
antigen binding specificity. Among the comparative methods, only ERGO2
exhibited moderate performance, with AUROC and AUPR values of 0.704 and 0.747,
respectively. The remaining methods, Pan-Peptide, pMTnet, and DLpTCR,
performed close to random guessing, indicating significant limitations in
their predictive capacity for TCR-antigen binding specificity.

To further validate the performance of UnifyImmun objectively, we compiled an
external independent test set consisting of over 100,000 TCR-pHLA binding
pairs collected from a number of publications. This comprehensive dataset
allowed us to assess the predictive capabilities of UnifyImmun toward real-
world scenario, beyond the initial hold-out test set. The experimental results
were consistent with our previous findings, further confirming the superior
performance of UnifyImmun. As shown in Figure 3(b), UnifyImmun achieved AUROC
and AUPR values exceeding 0.9, significantly outperforming the second-best
method, ERGO2. The other three methods—PanPep, pMTnet, and DLpTCR—exhibited
even poorer performance. We even observed negative MCC values for PanPep and
DLpTCR, indicating a high degree of disagreement between the model predictions
and the actual outcomes. This observation further underscores the limitations
of these comparative methods, which in turn validated the superior performance
of UnifyImmun in predicting TCR-antigen binding specificity.

Moreover, we computed the positive predictive value (PPV) for the top-ranked
predicted positive TCR-antigen samples on two distinct datasets. Specifically,
we evaluated the PPV for the top 100, top 1000, and top 5000 predictions
(Figure 3(c-d)). UnifyImmun achieved an impressive 100%, 99% and 97% PPV
values for the top 100, top 1000 and top 5000 predictions, respectively. In
contrast, the other methods performed significantly worse in this aspect.
Additionally, we also presented the ROC curves and AUC values for these
methods on both test sets (Figures 3(e-f)). The UMAP feature visualization of
peptide-HLA pairs implied that the positive and negative samples separated
remarkably in the latent space (Figure 2(g-h)). It can be found that
UnifyImmun exhibited superior performance in predicting TCR-antigen binding
specificity compared to the other methods. It is indeed noteworthy that while
the comparative methods may exhibit better performance on smaller datasets,
their performance decreased seriously when applied to large-scale test set.
his suggests they suffered from weak generalization ability and struggle to
adapt to large real-world data scenarios. In contrast, UnifyImmun demonstrated
strong robustness across distinct datasets, offering a more dependable and
precise tool for predicting TCR-antigen binding specificity.

Figure 3: Performance evaluation on predicting peptide-TCR binding
specificity. (a-b) Performance comparison to four methods on hold-out and
external test dataset, respectively. (c-d) Positive predictive value (PPV) for
the top 100, top 1000, and top 5000 predicted samples on hold-out and external
test dataset, respectively. (e-f) ROC curves and AUC values on hold-out and
external test dataset, respectively. (g-h) UMAP feature visualization of
peptide-CDR3 pairs.

###  2.3 Two-phase progressive training effectively improve model performance

Due to the limited number of HLA-antigen-TCR triplet samples for model
training, we devised a two-phase progressive training strategy aimed at
effectively leveraging the HLA-antigen and TCR-antigen binding samples
available. To validate the performance enhancement by two-stage progressive
training strategy, we randomly divided the benchmark dataset into training and
test sets for model training and subsequent assessment. This process was
independently repeated ten times to account for variations introduced by
random data partitioning. Next, We presented the results in boxplots for each
training round. This approach allows us to confidently evaluate the
effectiveness of our methodology while mitigating any potential biases
introduced by random data partition.

As depicted in Figures 4, the model performance was suboptimal in the absence
of alternating training (Round 0). However, as the number of training rounds
increased, we observed a marked and consistent improvement in model
performance until it converged at a high level. Specifically, on the HLA-
antigen hold-out test set, the two-phase alternative training quickly boosted
both AUROC and AUPR values (Figures 4(a-b)). For TCR-antigen binding
prediction, similar trends were observed on both hold-out and external test
datasets (Figures 4(c-f)). To further validate the effectiveness of our
training strategy, we evaluated the model’s predictive power on a COVID-19
test set (for data detail see Section 2.5). As expected, the results showed
that the performance on the COVID-19 test set improved with an increasing
number of training rounds (Figures 4(g-h)). Notably, while the variance in
performance metrics was higher in the early rounds, it decreased progressively
as the number of alternating rounds increased. This indicates that the model’s
performance became less affected by random data partitioning. These findings
confirm that the two-stage progressive training strategy drives the model to
learn more effective feature representations, thereby achieving more stable
and reliable performance.

Figure 4: Two-phase progressive training improved performance for both
peptide-HLA and peptide-TCR binding prediction tasks. (a-b) AUROC and AUPR
values increased with two-phase training rounds on peptide-HLA hold-out test
set. (c-d) AUROC and AUPR values increased with two-phase training rounds on
peptide-TCR external independent test set. (e-f) AUROC and AUPR values
increased with two-phase training rounds on peptide-TCR hold-out test set.
(g-h) AUROC and AUPR values increased with two-phase training rounds on
peptide-TCR COVID-19 test set.

###  2.4 Cross-attentions and integrated gradients reveal important locus

We employed a cross-attention mechanism to integrate information from antigens
with HLA/TCR, aiming to explore whether cross-attention scores reflect the
crucial positions and amino acids in the antigen sequence for its binding
affinity to HLA and TCR receptors. For this purpose, we aggregated the
attention scores for each amino acid at every position across all antigen
sequences. Higher scores indicate strong influence on the binding affinity to
the receptors. To accommodate variations in antigen length, we generated
separate heatmaps based on the cross-attention scores for antigens of lengths
9-14 (Figure 5(a-b) and Figure S3).

For binding to HLA, the attention scores are significantly higher at the
second and final position of the antigen sequence, indicating that these two
loci are particularly important for strong binding affinity. Notably, the Leu
amino acid (L) consistently receives higher attention scores, especially at
these critical positions, emphasizing its significance in antigen-HLA
interactions. When considering TCR-antigen binding, we also found that the Leu
amino acid at the second position [50] stands out as the most influential. To
gain a more comprehensive understanding of the amino acids contributions, we
calculated the cumulative attention scores for each amino acid across all
positions in antigens of specific lengths. This metric reflects the overall
importance of an amino acid in mediating antigen-receptor binding. Figures
5(e-f) graphically illustrated the attention scores for all 20 amino acids in
the context of HLA and TCR binding, respectively. It is evident that Leu
consistently plays a pivotal role in antigen binding to both receptors,
regardless of antigen length.

To further validate the findings from cross-attention heatmaps, we calculated
the Integrated Gradients (IG) values for each amino acid at every position
within the antigen sequences. As shown in Figure 5(e-f), there were strong
correlations between the IG values and the attention score heatmaps. For
instance, the Leu amino acid at the second position within the antigen
sequence emerged as particularly important for the binding affinity to both
HLA and TCR molecules. The consistency between the attention scores and the IG
values reinforced the validity of our approach and highlights the crucial role
of specific amino acids and their positions in mediating antigen-receptor
affinity.

To intuitively illustrate the crucial positions and amino acids involved in
TCR-pHLA binding, we obtained the crystal structure of the TK3 TCR in complex
with HLA-B*3501/HPVG (PDB ID: 3MV7) from the PDB database. We extracted the
HLA allele, antigen (11-mer), and the CDR3$\beta$ chain (11-mer) and predicted
their binding scores using UnifyImmun. The results indicated a high binding
score between HLA and antigen (0.99), as well as a moderate score between the
CDR3$\beta$ chain and antigen (0.62). Upon careful inspection of the three-
dimensional structure (Figure 5(g)), we observed that the 9-th amino acid Phe
(Y) of the antigen is embedded within the HLA binding groove (represented by
blue helices) and stabilized through hydrogen bonds (yellow lines). This
observation illustrated the chemical mechanism underlying the binding of the
antigen to HLA. Furthermore, the attention heatmap for the HLA-antigen-TCR
ternary complex (Figure 5(h)) revealed a significant cross-attention score
between the 8-th amino acid Tyr (Y) of the antigen and the 8-th amino acid Gly
(G) on the CDR3$\beta$ chain. The crystal structure demonstrated the formation
of hydrogen bonds associated with these two amino acids, indicating their
crucial role in the formation of the pTCR complex, despite their distance of
8.77$\mathring{A}$ (distance between two C$\alpha$ atoms of two amino acids)
is slightly beyond the conventional contact threshold of 6$\mathring{A}$. In
summary, the cross-attention mechanism offers opportunity to explore the
global dependencies between TCR-pHLA interactions, thereby enhancing the
interpretability of our model.

Figure 5: Heatmaps based on cross-attention scores and integral gradients.
(a-b) Accumulative contribution (attention score sum across peptide length) of
each amino-acid type of peptide to pHLA and pTCR bindings, respectively. (c-d)
Attention score-based heatmaps of the amino-acid type at each position of
9-mer peptide to pHLA and pTCR bindings, respectively. (e-f) Integral
gradient-based heatmaps of the amino-acid type at each position of 9-mer
peptide to pHLA and pTCR bindings, respectively. (g-h) Attention score-based
heatmap and 3D structure for TCR complex with HLA-B*35:01/HPVG (PDB ID: 3MV7).

###  2.5 UnifyImmun generalizes to COVID-19 antigen-TCR binding specificity

To validate the generalizability of UnifyImmun, we tested its ability to
predict the bindings between virus-derived antigens and TCRs. We collected a
total of 520,000 positive interactions between antigens derived from COVID-19
virus and human TCRs from the ImmuneCODETM database. Meanwhile, we generated
an equal number of negative samples via random shuffle, creating a million-
scale COVID-19 test set. We compared UnifyImmun with several other methods,
including PanPen, ERGO2, DLpTCR, and pMTnet. Due to the low efficiency of
pMTnet and DLpTCR, they were unable to tackle million-scale test set within a
reasonable time frame. Therefore, we randomly selected 100,000 pair as their
test set for evaluation.

The results showed that UnifyImmun achieved AUROC and AUPR values of 0.644 and
0.632, respectively (Figure 6(a)). In contrast, other methods obtained AUROC
and AUPR values only marginally above 0.5, which is close to random guessing.
Furthermore, we computed the PPV values for the top 100, top 1000, and top
5000 predictions made by each method. Our model consistently achieved 90% PPV
values, remarkably outperformed all the comparative methods whose PPV values
were always less than 60%. Overall, this significant difference provide strong
evidence of the robust generalizability of UnifyImmun. Its ability to
accurately predict antigen-TCR bindings, even for external virus-derived
antigens, highlights its potential for facilitating the development of
effective immune-based therapies and vaccines against COVID-19 viruses,
including .

Figure 6: Performance comparison to four existing methods on COVId-19 antigen-
TCR CDR3 binding data set. (a) AUROC and AUPR values on COVID-19 dataset.
(b)PPV of top100, top1000 and top5000 predicted samples on COVID-19 dataset.

###  2.6 Validation of UnifyImmun predictive power on immunotherapy responses
and clinical outcomes

The antigen presentation to cytotoxic T-cells plays a pivotal role in
determining the efficacy of tumor immunotherapy, particularly with immune
checkpoint inhibitors. To evaluate the predictive power of UnifyImmun, we
conducted a comprehensive study on a metastatic melanoma cohort (MM-HLA) [51]
and a urothelial carcinoma cohort (UC) [52]. The MM-HLA cohort included 110
patients, with each individual harboring an average of 919 potential
neoantigens. The UC cohort included 18 patients, with each patient harboring
115 potential neoantigens. For each patient, we obtained the HLA typing, tumor
antigen sequences, and the clinical immunotherapy responses. For the MM-HLA
cohort, we utilized the RECIST criteria to categorize patients into four
distinct groups: complete response (CR), partial response (PR), stable disease
(SD), and progressive disease (PD). We predicted the binding scores for all
possible HLA-peptide pairs using UnifyImmun and visually represented the
predicted scores for each patient group. As shown in Figure 7(a), analysis of
variance (ANOVA) with an $F$-test revealed statistically significant
differences in the HLA-peptide binding affinity between these groups. Notably,
the PD group demonstrated a highly statistical divergence compared to the
other patient groups (CR, PR and SD). This observation underscored the
distinct differences in neoantigen presentation between benefit vs non-benefit
patient groups from immunotherapy. Moreover, the patients in the CR group
exhibited the prevalence of high-affinity neoantigens, while those in the PD
group showed many low-affinity antigens. Furthermore, we stratified the
patients based on response, non-response, and long survival (PFS$<$180 days
but OS$>$2 years), and the violin plots further confirmed that the three
groups exhibited distinct patterns in antigen binding affinity to HLA
molecules (Figure 7(b)). In the UC cohort, the CR group also exhibited a
notable presence of high-affinity antigens (Figure 7(c)), and ANOVA analysis
further showed noteworthy differences between these groups. Notably, while the
differences between the PR and SD groups were not as pronounced, the PD group
exhibited the highly significant divergence from the others. Figure 7(d)
visualized the distinct patterns of antigen presentation between the clinical
response and non-response groups. The findings supported the clinical utility
of UnifyImmun in guiding treatment decisions and ultimately improving patient
outcomes.

To further validate the predicted antigen-TCR pairs, we conducted analysis on
additional advanced melanoma cohort (MM-TCR) that had received immunotherapy
[53]. This cohort comprised 29 patients who underwent TCR-seq and genomic
sequencing. Using the amino acid resulting from a missense mutation as an
anchor, we generated all possible 9-mer peptides harboring this anchor site.
After extracting the CDR3 sequences from TCR-seq data, we created all possible
CDR3-peptide pairs for each patient, yielding a total of 81,851,486 pairs. We
used UnifyImmun to score these CDR3-peptide pairs, and selected the top 5000
highest-scoring peptide-TCR pairs for each patient. Next, we categorized the
patients into CR, PR, SD, and PD groups based on RECIST criteria, and plotted
the boxplots of the predicted scores for each group (Figure 7(e)). ANOVA
analysis revealed significant differences among the groups ($F$-test), with
the CR and PR groups exhibiting notably higher scores than the SD and PD
groups. By stratifying the patients into benefit, non-benefit, and long-term
survival groups, we generated the violin plots (Figure 7(f)) and found that
patients in the long-term survival group exhibited peptide-TCR pairs with
higher binding scores.

Finally, to confirm the correlation between high-confidence HLA- and TCR-
peptide pairs by UnifyImmun and improved clinical outcomes, we conducted
survival analysis on two melanoma cohorts (MM-HLA and MM-TCR). We considered
the top 2% HLA- and TCR-peptide pairs as high-confidence bindings, and
stratified the patients with such bindings into the high-confidence group,
while the remaining patients into the low-confidence group. The results of
survival analysis showed that patients in the high-confidence group exhibited
significantly higher overall survival ratio compared to the low-confidence
group (Figure 7(e-f)). The p-values were 0.0038 and 0.031 for two cohorts,
respectively. These findings suggest a strong association between the HLA- and
TCR-peptide bindings predicted by UnifyImmun and favorable immunotherapy
responses as well as clinical outcomes.

Figure 7: Correlation between UnifyImmun predicted binding scores and
immunotherapy response and clinical outcomes on three cohorts. (a-b) Boxplots
and violin plots of predicted neoantigen-HLA binding scores regarding the
different immunotherapy response groups for MM-HLA cohort, respectively. (c-d)
Boxplots and violin plots of predicted neoantigen-HLA binding scores regarding
the different immunotherapy response groups for UC cohort, respectively. (e-f)
Boxplots and violin plots of predicted antigen-TCR binding scores regarding
the different immunotherapy response groups for MM-TCR cohort, respectively.
(g-h) Survival curves between stratified patient groups with high- and low-
confidence antigen binding specificity on MM-HLA and MM-TCR cohorts.

##  3 Discussion and Conclusion

In this study, we introduce UnifyImmun, a unified cross-attention-based model
designed to simultaneously predict the binding specificity of antigens to both
HLA and TCR molecules. We have devised a two-phase progressive training
strategy by which the two tasks mutually reinforce each other, guiding the
encoders to capture more distinctive features. To bolster the model’s
generalizability, we have incorporated virtual adversarial perturbation into
the framework. When benchmarked against over ten existing methods for antigen-
HLA and antigen-TCR binding prediction, our method consistently outperforms
them in both tasks on hold-out test set and external independent set, as well
as HPV and COVID-19 test set. Additionally, the cross-attention scores
pinpoint the amino acid sites crucial for antigen binding to receptors.

Most importantly, our model integrates the prediction tasks of antigen
presentation and T-cell activation into a unified framework, offering more
comprehensive evaluation of antigen immunogenicity compared to previous models
that only considered individual tasks alone. However, it is important to note
that the antigen-induced immune system activation involves a series of
cascading biological process, with numerous factors influencing the
immunogenicity of antigens. These include endopeptidase preferences for
polypeptide cleavage sites, antigen concentration, stability of pHLA
complexes, and transporter protein efficiency, all of which affect the degree
of immune response activation. Although our model marks a significant
advancement in the holistic assessment of antigen immunogenicity, it still
represents a high-level simplification of the actual immune response process.
This limitation could account for the suboptimal performance of our method on
some test sets.

The currently available TCR CDR3 sequences constitute just a fraction of the
immense human TCR repertoire. Despite our model could capture some antigen
binding patterns from the TCR sequences, its capacity is hindered by the
scarcity of available data. This constraint becomes particularly serious when
confronted with unseen sequences in the test set, resulting in suboptimal
prediction performance. Fortunately, Fortunately, the remarkable progress in
single-cell transcriptome sequencing has led to a significant increase in
scRNA-seq data of T cells, greatly facilitating the acquisition of CDR3
sequences. By leveraging the power of large language models (LLM) for pre-
training, we can extract more expressive and meaningful features from the
wealth of sequences. This, in turn, has the potential to significantly enhance
the predictive capabilities of our model, enabling it to more accurately
assess the immunogenicity of antigens.

##  4 Methods

###  4.1 Dataset

We created a benchmark dataset of HLA-peptide bindings from over ten previous
studies (for more details see Supplementary Table S1). After removing
duplicates and abnormal sequences (such as missing values or asterisk), we
obtained 423,843 pHLA binders, including 139 HLA alleles and 219,916 unique
peptides. The dataset was subsequently split into the training set and hold-
out test set by 9:1 ratio. As a result, the training set contained 381,287
pairs, spanning 139 HLA alleles and 219,916 antigens. The test set contained
42,556 pairs, covering 124 HLA alleles and 34,052 antigens.

Since previous studies have verified that both $\alpha$ and $\beta$ chains of
TCR are crucial for antigen recognition [25, 54], we considered both chains
and treated them as single CDR3 chains. This approach differed from existing
models that only considered the $\beta$ chain. Consequently, we created a TCR-
peptide binding dataset with 109,554 pairs, which was split into training and
hold-out test set by 9:1 ratio. The training set covered 881 distinct antigens
and 93,727 unique CDR3 sequences, while the test set covered 496 distinct
antigens and 10,896 unique CDR3 sequences. In addition, we built an
independent test set composed of HLA-antigen-TCR triplet, which included
63,324 pairs covering 57 HLA alleles, 998 distinct antigens, and 53,439 CDR3
sequences.

To our best knowledge, both the HLA-antigen and TCR-antigen binding datasets
we built are the largest to date. We shuffled the HLA (TCR) and peptide
sequences, and randomly mismatched them to generate negative samples. Although
false negative samples may be generated, the possibility and proportion of
such samples are very low and can be ignored. We generated the same number of
negative samples to positive samples to obtain balanced datasets.

###  4.2 Sequence embedding

The HLA pseudo sequences have a fixed length of 34 amino acids. Each amino
acid is mapped to a 64-dimensional embedding via a character embedding layer.
Since the order of amino acids is critical to the protein structure and
function, the sine and cosine positional encoding is applied to each position.
The amino-acid embedding and positional embedding are summed to obtain the
sequence embedding.

The peptides are padded to maximum length of 15 to handle the variable input
length, and then each amino acid is mapped to a 64-dimensional embedding.
Similarly, the positional encoding is executed to incorporate positional
information of the amino acids. After the padding and embedding steps, each
peptide is represented as a embedding matrix.

All the TCR CDR3 sequences shorter than 34 amino acids are padded to 34, while
a small number of CDR3 sequences exceeding 34 amino acids are truncated. Next,
the similar embedding process is applied to each CDR3 sequence.

###  4.3 Transformer-based encoder

The encoder is based on self-attention mechanism [55], which has shown
exceptional capability in extracting global correlation and dependency
relationships from protein sequences [56, 57, 58]. Self-attention mechanism
learns the attention of all possible amino acid pairs in the input sequence.
The attention scores are computed from the normalized dot product of query
vectors $Q$ and key vectors $K$ followed by a softmax operation. The self-
attention layer output the weighted sum of the value vectors $V$ by the
attention scores. The operations of a self-attention layer written in matrix
form are as follows:

$\mathrm{Attention}(Q,K,V)=\mathrm{softmax}\left(\frac{\mathrm{QK}^{T}}{\sqrt{d_{k}}}\right)V$
(1)

where $d_{k}$ is the dimension of the vectors (chosen as 64). Taking HLA as an
example, the $Q$, $K$, $V$ are all set to its 34$\times$64 embedding matrix.
Subsequently, the output of the self-attention block is passed through
multiple fully-connected layers, by which the dimension of the gyro first
rises and then falls. For peptides and TCR sequences, we adopt the similar
encoder architecture for feature extraction.

It is important to highlight that we introduce the mask mechanism in
calculating the self-attention scores. Specifically, for peptides or CDR3
sequences shorter than their respective maximum lengths, we exclude non-amino-
acid characters from consideration during model training. For this purpose, we
assign zero attention scores corresponding to these characters, so that they
do not influence the computation of attention scores. In our implementation,
the encoder comprises a one-layer one-head self-attention block.

###  4.4 Cross-attention for feature fusion

Cross-attention mechanism can effectively capture the intricate relationships
and dependencies between different sequences [59, 60]. In our work, we
leverage the cross-attention mechanism to extract crucial information
regarding the interactions between peptide and HLA/TCR molecules. For the
integration of HLA and peptide feature, the HLA embedding act as the $K$ and
$V$ matrices, while the peptide embedding serves as the $Q$ matrix.
Subsequently, the $V$ matrix is weighted based on the computed cross-attention
scores between $Q$ and $K$. The cross-attention layer’s output passes through
three fully-connected layers: first expanding to 512 dimensions and then
compressing to 64 for a condensed representation. The similar process is
applied for the fusion of TCR and peptide features, where TCR embedding act as
$K$ and $V$ matrices, and peptide embedding server as $Q$ matrix.

We also employ mask mechanism in calculating the cross-attention scores. The
mask mechanism greatly reduces the computational overhead and accelerates the
model’s convergence rate. In our implementation, we use one-head and one-layer
cross-attention mechanism, as illustrated in Figure 1.

###  4.5 Prediction of binding specificity

To predict the bindings between peptides and HLA (or TCR) molecules, we
flatten the fused matrix of HLA (or TCR)-peptide pairs outputted by the cross-
attention block, resulting in a 2176-dimensional vector (aka 34×64). This
vector then passes through three fully-connected layers with 256, 64, and 2
nodes, respectively, utilizing the ReLU activation function. The final output
is obtained via a Softmax layer. We adopt cross-entropy as the loss function,
and use the Adam optimizer with a learning rate of 1e-3.

The model training is conducted on a CentOS Linux 8.2.2004 (Core) system,
equipped with an Intel(R) Xeon(R) Silver 4210R CPU operating at 2.40GHz, along
with a GeForce RTX 4090 GPU and 128GB of memory. The model is implemented
using PyTorch 2.2.1.

###  4.6 Virtual adversarial training

Virtual adversarial training [61] introduces subtle perturbations within the
vicinity of the sequence embedding space, rather than directly perturbing the
original sequences. The perturbations are oriented towards the direction of
loss gradient ascent and are typically generated under L2 norm constraints.
This training strategy demands that the model not only minimizes the original
loss but also minimizes the adversarial loss, making the model less sensitive
to slight changes in the input. We applied adversarial perturbations to the
embeddings of all three types of sequences, so that the encoder learn to
extract discriminative features. Our ablation experiments have confirmed that
virtual adversarial learning indeed improves model performance, as shown in
Table S1.

###  4.7 Two-stage progressive training

In the first phase, the model was trained exclusively on the HLA-peptide
pairs, keeping the TCR encoder and the TCR-antigen cross-attention module
fixed. This enforced the model to concentrate solely on learning the intricacy
of HLA-antigen interactions. In the second phase, the model were trained
exclusively using the TCR-peptide pairs, with the HLA encoder and HLA-antigen
cross-attention module fixed. These two phases were alternated until the model
performance converged. Note that throughout the alternating training process,
the antigen encoder remained continuously updated and shared between both HLA-
antigen and TCR-antigen binding prediction tasks. By iteratively refining the
antigen encoder parameters, the model learned to capture the essential
information pertinent to both HLA and TCR binding and thereby enhanced its
overall predictive accuracy.

###  4.8 Model ablation experiments

To validate the contributions of different components, we conducted ablation
studies on our model to assess the performance in predicting HLA-antigen and
TCR-antigen bindings. Specifically, we evaluate the attention masking,
positional encoding, and virtual adversarial perturbation individually. The
results obtained by the ablated models for antigen-HLA and antigen-TCR CDR3
binding are outlined in Table 1 and Table 2, respectively.

The findings reveal that the removal of any component leads to a decrease in
performance when compared to UnifyImmun. Notably, the removal of virtual
adversarial perturbation has the most significant impact, resulting in at
least 4% drop in AUROC for both HLA-antigen and antigen-TCR binding
predictions. This strongly indicates that virtual adversarial training
contributes to improving the model’s overall performance and generalization
capabilities. Furthermore, we observed that the absence of the attention
masking leads to increased computational overhead during the training process.
Without the mask, additional computational resources are expended to process
padding sequences, leading to a rise in computational cost.

## References

  *   * [1] Yang, K., Halima, A. & Chan, T. A.  Antigen presentation in cancer—mechanisms and clinical implications for immunotherapy.  _Nature Reviews Clinical Oncology_ 20, 604–623 (2023).
  * [2] Barry, M. & Bleackley, R. C.  Cytotoxic t lymphocytes: all roads lead to death.  _Nature Reviews Immunology_ 2, 401–409 (2002).
  * [3] Raskov, H., Orhan, A., Christensen, J. P. & Gögenur, I.  Cytotoxic cd8+ t cells in cancer and cancer immunotherapy.  _British journal of cancer_ 124, 359–367 (2021).
  * [4] Weigelin, B. _et al._ Cytotoxic t cells are able to efficiently eliminate cancer cells by additive cytotoxicity.  _Nature communications_ 12, 5217 (2021).
  * [5] Rowen, L., Koop, B. F. & Hood, L.  The complete 685-kilobase dna sequence of the human $\beta$ t cell receptor locus.  _Science_ 272, 1755–1762 (1996).
  * [6] Glanville, J. _et al._ Identifying specificity groups in the t cell receptor repertoire.  _Nature_ 547, 94–98 (2017).
  * [7] Fowell, D. J. & Kim, M.  The spatio-temporal control of effector t cell migration.  _Nature Reviews Immunology_ 21, 582–596 (2021).
  * [8] Lim, A. R., Rathmell, W. K. & Rathmell, J. C.  The tumor microenvironment as a metabolic barrier to effector t cells and immunotherapy.  _Elife_ 9, e55185 (2020).
  * [9] Farber, D. L., Yudanin, N. A. & Restifo, N. P.  Human memory t cells: generation, compartmentalization and homeostasis.  _Nature Reviews Immunology_ 14, 24–35 (2014).
  * [10] Liu, Q., Sun, Z. & Chen, L.  Memory t cells: strategies for optimizing tumor immunotherapy.  _Protein & cell_ 11, 549–564 (2020).
  * [11] Chen, D. S. & Mellman, I.  Oncology meets immunology: the cancer-immunity cycle.  _immunity_ 39, 1–10 (2013).
  * [12] Yewdell, J. W. & Bennink, J. R.  Immunodominance in major histocompatibility complex class i–restricted t lymphocyte responses.  _Annual review of immunology_ 17, 51–88 (1999).
  * [13] Dunn, G. P., Old, L. J. & Schreiber, R. D.  The three es of cancer immunoediting.  _Annu. Rev. Immunol._ 22, 329–360 (2004).
  * [14] Trowsdale, J.  Hla genomics in the third millennium.  _Current opinion in Immunology_ 17, 498–504 (2005).
  * [15] Huppa, J. B. _et al._ Tcr–peptide–mhc interactions in situ show accelerated kinetics and increased affinity.  _Nature_ 463, 963–967 (2010).
  * [16] Nikolich-Žugich, J., Slifka, M. K. & Messaoudi, I.  The many important facets of t-cell repertoire diversity.  _Nature Reviews Immunology_ 4, 123–132 (2004).
  * [17] Zhang, S.-Q. _et al._ Direct measurement of t cell receptor affinity and sequence from naïve antiviral t cells.  _Science translational medicine_ 8, 341ra77–341ra77 (2016).
  * [18] Davis, M. M. & Bjorkman, P. J.  T-cell antigen receptor genes and t-cell recognition.  _Nature_ 334, 395–402 (1988).
  * [19] Krogsgaard, M. & Davis, M. M.  How t cells’ see’antigen.  _Nature immunology_ 6, 239–245 (2005).
  * [20] Chowell, D. _et al._ Evolutionary divergence of hla class i genotype impacts efficacy of cancer immunotherapy.  _Nature medicine_ 25, 1715–1720 (2019).
  * [21] Krishna, C., Chowell, D., Gönen, M., Elhanati, Y. & Chan, T. A.  Genetic and environmental determinants of human tcr repertoire diversity.  _Immunity & Ageing_ 17, 1–7 (2020).
  * [22] Purcell, A. W., Ramarathinam, S. H. & Ternette, N.  Mass spectrometry–based identification of mhc-bound peptides for immunopeptidomics.  _Nature protocols_ 14, 1687–1707 (2019).
  * [23] Zhang, S.-Q. _et al._ High-throughput determination of the antigen specificities of t cell receptors in single cells.  _Nature biotechnology_ 36, 1156–1159 (2018).
  * [24] Kula, T. _et al._ T-scan: a genome-wide method for the systematic discovery of t cell epitopes.  _Cell_ 178, 1016–1028 (2019).
  * [25] Hudson, D., Fernandes, R. A., Basham, M., Ogg, G. & Koohy, H.  Can we predict t cell specificity with digital biology and machine learning?  _Nature Reviews Immunology_ 23, 511–521 (2023).
  * [26] Chu, Y. _et al._ A transformer-based model to predict peptide–hla class i binding and optimize mutated peptides for vaccine design.  _Nature Machine Intelligence_ 4, 300–311 (2022).
  * [27] O’Donnell, T., Rubinsteyn, A., Bonsack, M., Riemer, A. & Hammerbacher, J.  Mhcflurry: open-source class i mhc binding affinity prediction.  _Cold Spring Harbor Laboratory_ (2017).
  * [28] Birkir, R., Bruno, A., Sinu, P., Bjoern, P. & Morten, N.  Netmhcpan-4.1 and netmhciipan-4.0: improved predictions of mhc antigen presentation by concurrent motif deconvolution and integration of ms mhc eluted ligand data.  _Nucleic Acids Research_ .
  * [29] Haoyang, Z. & Gifford, D. K.  Deepligand: accurate prediction of mhc class i ligands using peptide embedding.  _Bioinformatics_ i278–i283 (2019).
  * [30] Jun, C., Kadre, B., Karola, R. & Brandon, M.  Bertmhc: improved mhc–peptide class ii interaction prediction with transformer and multiple instance learning.  _Bioinformatics_ (2021).
  * [31] Gao, Y. _et al._ Pan-peptide meta learning for t-cell receptor–antigen binding recognition.  _Nature Machine Intelligence_ 5, 236–249 (2023).
  * [32] Lu, T. _et al._ Deep learning-based prediction of the t cell receptor–antigen binding specificity.  _Nature machine intelligence_ 3, 864–875 (2021).
  * [33] Zhaochun, X. _et al._ Dlptcr: an ensemble deep learning framework for predicting immunogenic peptide recognized by t cell receptor.  _Briefings in Bioinformatics_ (2021).
  * [34] Springer, I., Tickotsky, N. & Louzoun, Y.  Contribution of t cell receptor alpha and beta cdr3, mhc typing, v and j genes to peptide binding prediction.  _Frontiers in Immunology_ 12, 664514 (2021).
  * [35] Weber, A., Born, J. & Martínez, M. R.  Titan: T-cell receptor specificity prediction with bimodal attention networks.  _Bioinformatics_ 37, i237–i244 (2021).
  * [36] Fang, Y., Liu, X. & Liu, H.  Attention-aware contrastive learning for predicting t cell receptor–antigen binding specificity.  _Briefings in Bioinformatics_ 23, bbac378 (2022).
  * [37] Massimo, A. & Morten, N.  Gapped sequence alignment using artificial neural networks: application to the mhc class i system.  _Bioinformatics_ 511\.
  * [38] Zhang, H., Lund, O. & Nielsen, M.  The pickpocket method for predicting binding specificities for receptors based on receptor pocket similarities.  _Bioinformatics_ (2009).
  * [39] Kim, Y., Sidney, J., Pinilla, C., Sette, A. & Peters, B.  Derivation of an amino acid similarity matrix for peptide:mhc binding and its application as a bayesian prior.  _BMC Bioinformatics_ 10 (2009).
  * [40] Peters, B. & Sette, A.  Generating quantitative models describing the sequence specificity of biological processes with the stabilized matrix method.  _BMC Bioinformatics_ 6, 132 – 132 (2005).
  * [41] Karosiene, E., Lundegaard, C., Lund, O. & Nielsen, M.  Netmhccons: a consensus method for the major histocompatibility complex class i predictions.  _Immunogenetics_ 64, 177–186 (2012).
  * [42] Rasmussen _et al._ Pan-specific prediction of peptide-mhc class i complex stability, a correlate of t cell immunogenicity.  _Journal of Immunology_ (2016).
  * [43] Moutaftsi, M. _et al._ A consensus epitope prediction approach identifies the breadth of murine t(cd8+)-cell responses to vaccinia virus.  _Nature Publishing Group_ (2006).
  * [44] Yan, H. _et al._ Acme: pan-specific peptide–mhc class i binding prediction through attention-based deep neural networks.  _Bioinformatics_ 23\.
  * [45] Jin, J. _et al._ Deep learning pan-specific model for interpretable mhc-i peptide binding prediction with improved attention mechanism.  _Proteins_ 89, 866–883 (2021).
  * [46] Vita, R. _et al._ The immune epitope database (iedb): 2018 update.  _Nucleic acids research_ 47, D339–D343 (2019).
  * [47] Bonsack, M. _et al._ Performance evaluation of mhc class-i binding prediction tools based on an experimentally validated mhc–peptide binding data set.  _Cancer immunology research_ 7, 719–736 (2019).
  * [48] Wells, D. K. _et al._ Key parameters of tumor epitope immunogenicity revealed through a consortium approach improve neoantigen prediction.  _Cell_ 183, 818–834 (2020).
  * [49] Wang, G. _et al._ Ineo-epp: a novel t-cell hla class-i immunogenicity or neoantigenic epitope prediction method based on sequence-related amino acid features.  _BioMed research international_ 2020 (2020).
  * [50] Parker, K. C. _et al._ Sequence motifs important for peptide binding to the human mhc class i molecule, hla-a2.  _Journal of immunology (Baltimore, Md.: 1950)_ 149, 3580–3587 (1992).
  * [51] Van Allen, E. M. _et al._ Genomic correlates of response to ctla-4 blockade in metastatic melanoma.  _Science_ 350, 207–211 (2015).
  * [52] Alexandra, S. _et al._ Contribution of systemic and somatic factors to clinical response and resistance to pd-l1 blockade in urothelial cancer: An exploratory multi-omic analysis.  _PLOS Medicine_ 14 (2017).
  * [53] Tumor and microenvironment evolution during immunotherapy with nivolumab.  _Cell_ 171, 934–949.e16 (2017).
  * [54] Zhao, X. _et al._ Tuning t cell receptor sensitivity through catch bond engineering.  _Science_ 376, eabl5282 (2022).
  * [55] Vaswani, A. _et al._ Guyon, I. _et al._ (eds) _Attention is all you need_.  (eds Guyon, I. _et al._) _Advances in Neural Information Processing Systems_ , Vol. 30 (Curran Associates, Inc., 2017).
  * [56] Madani, A. _et al._ Large language models generate functional protein sequences across diverse families.  _Nature Biotechnology_ 41, 1099–1106 (2023).
  * [57] Brandes, N., Ofer, D., Peleg, Y., Rappoport, N. & Linial, M.  Proteinbert: a universal deep-learning model of protein sequence and function.  _Bioinformatics_ 38, 2102–2110 (2022).
  * [58] Rives, A. _et al._ Biological structure and function emerge from scaling unsupervised learning to 250 million protein sequences.  _Proceedings of the National Academy of Sciences_ 118, e2016239118 (2021).
  * [59] Honda, S., Koyama, K. & Kotaro, K.  _Cross attentive antibody-antigen interaction prediction with multi-task learning_ (2020).
  * [60] Dens, C., Laukens, K., Meysman, P. & Bittremieux, W.  A cross-attention transformer encoder for paired sequence data.  _bioRxiv_ 2023–12 (2023).
  * [61] Miyato, T., Dai, A. M. & Goodfellow, I.  _Adversarial training methods for semi-supervised text classification_ (2017).


# Brief Note on Thurston Geometries in 3D Quadratic Curvature Theories

Gökhan Alkaç<EMAIL_ADDRESS>Physics Engineering Department,
Faculty of Engineering, Hacettepe University, 06800, Ankara, Turkey  Deniz
Olgu Devecioğlu<EMAIL_ADDRESS>School of Physics, Huazhong University of
Science and Technology,
Wuhan, Hubei, 430074, China

###### Abstract

We show that Thurston geometries are solutions to a large class of 3D
quadratic curvature theories, where New Massive Gravity, which was studied in
Flores-Alfonso:2021opl , is a special case.

Thurston’s conjecture Thurston  states that a three-manifold with a given
topology has a canonical decomposition into eight manifolds with high
symmetry, the so-called Thurston geometries (TGs) Thurston2 ; Thurston3 ;
Thurston4 . Based on Hamilton’s Ricci flow Hamilton , which provides a way to
deform a given metric in the form of a generalized heat equation, Perelman
gave his celebrated proof Perelman1 ; Perelman2 . Since many of TGs are not
Einstein spaces, one might wonder if a ‘master theory’ that admits all TGs as
solutions can be achieved by a modification of Einstein’s equation.

In a recent letter Flores-Alfonso:2021opl , it was reported that TGs are
solutions to New Massive Gravity (NMG) Bergshoeff:2009hq , where the authors
state that, to the best of their knowledge, NMG is the only 3D gravity theory
with this property111In Flores-Alfonso:2021opl , ghost-free theories were
considered for simplicity Flores-AlfonsoP .. In a previous work Flores-
Alfonso:2020zif , four of the TGs were shown to be solutions of the most
general quadratic theory. In this paper we shall close a gap and show that the
remaining four TGs are also solutions of the most general quadratic theory,
whose action is displayed in (4). NMG is singled out by being the only ghost-
free theory in this class222NMG is obtained after imposing the absence of the
scalar mode, which is always a ghost, and the unitarity of the massive spin-2
mode around a maximally symmetric background spacetime (see Gullu:2010sd for
an analysis on a flat background and Bergshoeff:2009aq ; Tekin:2016vli for
constant curvature backgrounds (dS/AdS)).. However, ghost-freedom is not
essential in defining a flow where TGs are critical points.

In order to demonstrate this, let us present the logic of a uniformization
theorem along the lines of Gegenberg:2003yz . Assuming that the manifold
admits some metric, one needs to devise a mechanism to flow the initial metric
into a highly symmetric one, which can be obtained by the following parabolic
system

$\partial_{t}g_{\mu\nu}=\mathcal{O}(g_{\mu\nu}),$ (1)

where $\mathcal{O}$ is an elliptic operator. The parabolicity of (1) assures
that the flow will end at fixed points as $t\rightarrow\infty$, for which
$\mathcal{O}(g_{\mu\nu})=0$. However, since $\mathcal{O}(g_{\mu\nu})$ should
transform as a tensor, it has to be constructed from curvature tensors, which
makes it impossible to obtain a strictly parabolic system and leads to
nonlinearity due to the inverse metric. Many issues regarding these
complications should be resolved for a full proof. In his proof, Perelman made
use of the operator

$\mathcal{O}(g_{\mu\nu})=-2\left(R_{\mu\nu}+\nabla_{\mu}\nabla_{\nu}f\right),$
(2)

where $f$ is a scalar field. It follows from the field equations of the action

$S[f,g]=\int\mathop{}\\!\mathrm{d}^{3}x\sqrt{|g|}\,e^{-f}\left[R+(\nabla
f)^{2}\right],$ (3)

and the action (3) is an increasing functional provided that
$\partial_{t}f=\partial_{t}\ln\sqrt{|g|}$. Alternatively, one can also try to
use the field equations of quadratic gravity theories as a candidate for such
an operator, which can be defined through an action that is a functional of
the metric tensor only. As long as one finds an operator that vanishes for all
TGs, the ghost-freedom of the theory plays no role in the construction.

Motivated by this, we consider the most general 3D quadratic gravity theory
described by the action

$S[g]=\int\mathop{}\\!\mathrm{d}^{3}x\sqrt{|g|}\left[\sigma(R-2\Lambda)+\alpha
R^{2}+\beta R_{\mu\nu}R^{\mu\nu}\right]\,,$ (4)

where $\sigma=\pm 1$ is introduced to control the sign of the Einstein-Hilbert
term and $(\alpha,\beta)$ are arbitrary constants. Field equations arising
from the action (4) are given by

$\displaystyle 0$ $\displaystyle=\sigma\left(G_{\mu\nu}+\Lambda
g_{\mu\nu}\right)+\alpha\left[2RR_{\mu\nu}-\dfrac{1}{2}g_{\mu\nu}R^{2}+2g_{\mu\nu}\square
R-2\nabla_{\mu}\nabla_{\nu}R\right]$ (5)
$\displaystyle+\beta\left[\dfrac{3}{2}g_{\mu\nu}R_{\rho\sigma}R^{\rho\sigma}-4R_{\mu}^{\
\rho}R_{\nu\rho}+\square R_{\mu\nu}+\dfrac{1}{2}g_{\mu\nu}\square
R-\nabla_{\mu}\nabla_{\nu}R+3RR_{\mu\nu}-g_{\mu\nu}R^{2}\right].$

As we summarize in Tables I and II, TGs are solutions for different choices of
($\sigma,\alpha,\beta$) that are subjected to two constraints:

  1. i.

The length scale in the metric should satisfy $\ell^{2}>0$ in order to
consider metrics with definite signature.

  2. ii.

Denominators in the expressions for $\ell^{2}$ and $\Lambda$ should be non-
zero.

As a simple example, to ensure that the theory admits Lorentzian TGs given in
Table II as solutions, the following choice of parameters is enough

$\sigma=-1,\qquad\beta=\dfrac{1}{m^{2}}>0,\qquad\alpha>-\dfrac{1}{2m^{2}}.$
(6)

When one demands to obtain Lorentzian TGs as solutions to a pure gravity
theory (no scalar mode) with unitary excitations around a maximally symmetric
spacetime, then NMG is the only choice and one has to fix the last parameter
as $\alpha=-\dfrac{3}{8m^{2}}$ Bergshoeff:2009hq ; Bergshoeff:2009aq ;
Tekin:2016vli , for which the scalar mode decouples from the spectrum since it
becomes infinitely heavy ($m_{s}^{2}\rightarrow\infty$), and there remains
only a unitary massive spin-2 graviton with
$m_{g}^{2}=m^{2}+\dfrac{\Lambda}{2}>0$.

For the Euclidean TGs, some possible choices of parameters are given in Table
III, where NMG is again a special case. Note that $m^{2}<0$ is just having an
imaginary mass. Analogous behaviour was also observed in Gegenberg:2002xj ;
Gegenberg:2003yz , where the Euclidean solutions are obtained from the low
energy limit of string theory supplemented by a one-form gauge field. For some
of the solutions, gauge fields become imaginary, which is equivalent to having
the kinetic terms with opposite signs in the action. However, this should not
be considered as a problem since the aim here is to choose the parameters
according to Table I and show that TGs correspond to fixed points.

TGs were also considered in another modification of Einstein’s equations
called Minimal Massive Gravity (MMG) Bergshoeff:2014pca , whose field
equations read

$G_{\mu\nu}+\Lambda\,g_{\mu\nu}+\frac{1}{\mu}C_{\mu\nu}+\frac{\gamma}{\mu^{2}}J_{\mu\nu}=0,$
(7)

where the Cotton tensor $C_{\mu\nu}$ is related to the Schouten tensor
$S_{\mu\nu}$ as

$C_{\mu\nu}=\epsilon_{\mu}^{\ \alpha\beta}\nabla_{\alpha}S_{\beta\nu},\qquad
S_{\mu\nu}=R_{\mu\nu}-\frac{1}{4}g_{\mu\nu}R,$ (8)

and the J-tensor is given by

$J_{\mu\nu}=\epsilon_{\mu}^{\ \alpha\beta}\epsilon_{\nu}^{\
\rho\sigma}S_{\alpha\rho}S_{\beta\sigma}.$ (9)

The field equations (7) do not arise from the variation of an action and are
covariantly conserved only on-shell. In Charyyev:2017uuu , it was shown that
all TGs except Sol geometry are solutions to (7) for certain (non-zero and
finite) choice of parameters. Sol geometry requires the vanishing of the
coefficient of the Cotton term ($\mu\rightarrow\infty$,
$\gamma\rightarrow\infty$ such that
$\frac{\gamma}{\mu^{2}}\rightarrow\text{constant}$), which in general ruins
the physical consistency of the theory. However, such a solution is still
acceptable if $\epsilon^{\mu\rho\sigma}S_{\rho}^{\ \tau}C_{\sigma\tau}=0$,
which is the case for the Sol geometry. Therefore, in addition to NMG, MMG
defines an operator $\mathcal{O}(g_{\mu\nu})$ that vanishes for all TGs. What
field equations of MMG lacks is the existence of an entropy functional along
the flow which admits a gradient formulation (see Kisisel:2008jx for an
example in Cotton flow), which might be crucial in an attempt for a full proof
of a uniformization theorem.

All in all, realising TGs as vacuum solutions to higher curvature
modifications of Einstein’s equations yield many possibilities where the
ghost-free theories are not necessarily privileged. Due to the higher-
derivative operators, it seems impossible to apply results from elliptic
operator theory for rigorous results. However, the relevance of the unitarity
of the theory might be checked through a detailed numerical analysis of the
flow equation as done in Kisisel:2008jx for the Cotton flow.

Geometry | Metric | $\ell^{2}$ | $\Lambda$
---|---|---|---
$\text{E}^{3}$ | $\mathop{}\\!\mathrm{d}s^{2}=\mathop{}\\!\mathrm{d}x^{2}+\mathop{}\\!\mathrm{d}y^{2}+\mathop{}\\!\mathrm{d}z^{2}$ | no length scale | 0
$\text{S}^{3}$ | $\mathop{}\\!\mathrm{d}s^{2}=\ell^{2}\left(\mathop{}\\!\mathrm{d}x^{2}+\sin^{2}x\mathop{}\\!\mathrm{d}y^{2}+\sin^{2}x\sin^{2}y\mathop{}\\!\mathrm{d}z^{2}\right)$ | arbitrary | $\dfrac{1}{\ell^{2}}-\sigma\dfrac{2(3\alpha+\beta)}{\ell^{4}}$
$\text{H}^{3}$ | $\mathop{}\\!\mathrm{d}s^{2}=\ell^{2}\left(\mathop{}\\!\mathrm{d}x^{2}+\cosh^{2}x\mathop{}\\!\mathrm{d}y^{2}+\cosh^{2}x\cosh^{2}y\mathop{}\\!\mathrm{d}z^{2}\right)$ | arbitrary | $-\dfrac{1}{\ell^{2}}-\sigma\dfrac{2(3\alpha+\beta)}{\ell^{4}}$
$\text{E}^{1}\times\text{S}^{2}$ | $\mathop{}\\!\mathrm{d}s^{2}=\ell^{2}\left(\mathop{}\\!\mathrm{d}x^{2}+\mathop{}\\!\mathrm{d}y^{2}+\sin^{2}y\mathop{}\\!\mathrm{d}z^{2}\right)$ | $\dfrac{1}{2\Lambda}$ | $-\dfrac{\sigma}{4(2\alpha+\beta)}$
$\text{E}^{1}\times\text{H}^{2}$ | $\mathop{}\\!\mathrm{d}s^{2}=\ell^{2}\left(\mathop{}\\!\mathrm{d}x^{2}+\mathop{}\\!\mathrm{d}y^{2}+\cosh^{2}y\mathop{}\\!\mathrm{d}z^{2}\right)$ | $-\dfrac{1}{2\Lambda}$ | $-\dfrac{\sigma}{4(2\alpha+\beta)}$
Nil | $\mathop{}\\!\mathrm{d}s^{2}=\dfrac{\ell^{2}}{4}\left[\mathop{}\\!\mathrm{d}x^{2}+\mathop{}\\!\mathrm{d}y^{2}+(x\mathop{}\\!\mathrm{d}y-\mathop{}\\!\mathrm{d}z)^{2}\right]$ | $-\dfrac{1}{2\Lambda}$ | $-\dfrac{\sigma}{8(\alpha+3\beta)}$
SL(2,$\mathbb{R}$) | $\mathop{}\\!\mathrm{d}s^{2}=\ell^{2}\left[\mathop{}\\!\mathrm{d}r^{2}+\sinh^{2}r\cosh^{2}r\mathop{}\\!\mathrm{d}\theta^{2}+(\mathop{}\\!\mathrm{d}\psi+\sinh^{2}r\mathop{}\\!\mathrm{d}\theta)^{2}\right]$ | $-\dfrac{25\alpha+23\beta}{10\Lambda(\alpha+\beta)}$ | $-\sigma\dfrac{25\alpha+23\beta}{200(\alpha+\beta)^{2}}$
Sol | $\mathop{}\\!\mathrm{d}s^{2}=\ell^{2}\left(e^{-2z}\mathop{}\\!\mathrm{d}x^{2}+e^{2z}\mathop{}\\!\mathrm{d}y^{2}+\mathop{}\\!\mathrm{d}z^{2}\right)$ | $-\dfrac{1}{2\Lambda}$ | $-\dfrac{\sigma}{8(\alpha+\beta)}$

Table 1: Euclidean Thurston Geometries as solutions to an arbitrary quadratic
gravity theory. Parameters should satisfy the constraints (i,ii), where
examples of sufficient but not necessary conditions are given in Table III.

Geometry | Metric | $\ell^{2}$ | $\Lambda$
---|---|---|---
Nil | $\mathop{}\\!\mathrm{d}s^{2}=\dfrac{\ell^{2}}{4}\left[\mathop{}\\!\mathrm{d}x^{2}+\mathop{}\\!\mathrm{d}y^{2}-(x\mathop{}\\!\mathrm{d}y-\mathop{}\\!\mathrm{d}z)^{2}\right]$ | $\dfrac{1}{2\Lambda}$ | $-\dfrac{\sigma}{8(\alpha+3\beta)}$
SL(2,$\mathbb{R}$) | $\mathop{}\\!\mathrm{d}s^{2}=\dfrac{\ell^{2}}{4}\left[-(\mathop{}\\!\mathrm{d}\Psi+\cos\Theta\mathop{}\\!\mathrm{d}\Phi)^{2}+\mathop{}\\!\mathrm{d}\Theta^{2}+\sin^{2}\Theta\mathop{}\\!\mathrm{d}\Phi^{2}\right]$ | $\dfrac{25\alpha+23\beta}{10\Lambda(\alpha+\beta)}$ | $-\sigma\dfrac{25\alpha+23\beta}{200(\alpha+\beta)^{2}}$
New Sol | $\mathop{}\\!\mathrm{d}s^{2}=\ell^{2}\left(e^{-2z}\mathop{}\\!\mathrm{d}x^{2}+e^{2z}\mathop{}\\!\mathrm{d}y^{2}-\mathop{}\\!\mathrm{d}z^{2}\right)$ | $\dfrac{1}{2\Lambda}$ | $-\dfrac{\sigma}{8(\alpha+\beta)}$
Lorentz Sol | $\mathop{}\\!\mathrm{d}s^{2}=\ell^{2}\left(2e^{-z}\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}z+e^{2z}\mathop{}\\!\mathrm{d}y^{2}\right)$ | arbitrary |  $\sigma=0$, no effect on the solution
Third Sol | $\mathop{}\\!\mathrm{d}s^{2}=\ell^{2}\left(-e^{2z}\mathop{}\\!\mathrm{d}y^{2}-2\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}y+\mathop{}\\!\mathrm{d}z^{2}\right)$ | $-4\sigma\beta$ | 0
Lorentz-Heisenberg | $\mathop{}\\!\mathrm{d}s^{2}=\dfrac{\ell^{2}}{4}\left(-\mathop{}\\!\mathrm{d}x^{2}+\mathop{}\\!\mathrm{d}y^{2}+(x\mathop{}\\!\mathrm{d}y-\mathop{}\\!\mathrm{d}z)^{2}\right)$ | $\dfrac{1}{2\Lambda}$ | $-\dfrac{\sigma}{8(\alpha+3\beta)}$
AdS | $\mathop{}\\!\mathrm{d}s^{2}=\ell^{2}\left[\mathop{}\\!\mathrm{d}r^{2}+\sinh^{2}r\cosh^{2}r\mathop{}\\!\mathrm{d}\theta^{2}-(\mathop{}\\!\mathrm{d}\psi+\sinh^{2}r\mathop{}\\!\mathrm{d}\theta)^{2}\right]$ | arbitrary | $-\dfrac{1}{\ell^{2}}-\sigma\dfrac{2(3\alpha+\beta)}{\ell^{4}}$

Table 2: Lorentzian Thurston Geometries as solutions to an arbitrary quadratic gravity theory. Parameters should satisfy the constraints (i,ii), where an example of sufficient but not necessary conditions is given in (6). Geometry | Parameters
---|---
$\text{E}^{1}\times\text{S}^{2}$ |  $\qquad\sigma=-1\qquad$ $\beta=\dfrac{1}{m^{2}}>0\qquad$ $\alpha>-\dfrac{1}{2m^{2}}\qquad$
$\sigma=+1\qquad$ $\beta=\dfrac{1}{m^{2}}<0\qquad$
$\,\alpha<-\dfrac{1}{2m^{2}}$
Others |  $\qquad\sigma=-1\qquad$ $\beta=\dfrac{1}{m^{2}}<0\qquad$ $\alpha<-\dfrac{1}{2m^{2}}\qquad$
$\qquad\sigma=+1\qquad$ $\beta=\dfrac{1}{m^{2}}>0\qquad$
$\alpha>-\dfrac{1}{2m^{2}}\qquad$
Table 3: Choice of parameters for Euclidean Thurston Geometries. Note that
these are sufficient but not necessary conditions.

Acknowledgements: D. O. D. was supported in part by the National Natural
Science Foundation of China under Grant No. 11875136 and the Major Program of
the National Natural Science Foundation of China under Grant No. 11690021. We
thank Nihat Sadık Değer and Daniel Flores-Alfonso for useful discussions. We
thank the anonymous referee for his comments, which led us to improve the
presentation of our work.

## References

  * (1) W. P. Thurston, “Three dimensional manifolds, Kleinian groups and hyperbolic geometry,” Bull. Amer. Math. Soc. 6, 357 (1982).
  * (2) W. P. Thurston, “The Geometry and Topology of Three-Manifolds,”
(1980, http://www.msri.org/publications/books/gt3m/)

  * (3) W. P. Thurston, ‘‘Three-Dimensional Geometry and Topology", Vol. 1 (Princeton University Press, Princeton NJ, 1997).
  * (4) P. Scott, ‘‘The Geometry of 3-Manifolds", Bull. London Math. Soc. 15, 401 (1983)
  * (5) Hamilton R.S., ‘‘Three-manifolds with positive Ricci curvature,” J. Diff. Geom. 17, 255 (1982).
  * (6) Perelman G., ‘‘The entropy formula for the Ricci flow and its geometric applications,” [arXiv/math:0211159]
  * (7) Perelman G., ‘‘Ricci flow with surgery on three-manifolds,” [arXiv/math:0303109]
  * (8) D. Flores-Alfonso, C. S. Lopez-Monsalvo and M. Maceda, ‘‘Thurston Geometries in Three-Dimensional New Massive Gravity,’’ Phys. Rev. Lett. 127, no.6, 061102 (2021) [arXiv:2104.00754]
  * (9) E. A. Bergshoeff, O. Hohm and P. K. Townsend, ‘‘Massive Gravity in Three Dimensions,’’ Phys. Rev. Lett. 102, 201301 (2009) [arXiv:0901.1766]
  * (10) D. Flores-Alfonso (private communication)
  * (11) D. Flores-Alfonso, C. S. Lopez-Monsalvo and M. Maceda, ‘‘Contact Geometry in Superconductors and New Massive Gravity,’’ Phys. Lett. B 815, 136143 (2021) [arXiv:2011.13499]
  * (12) I. Gullu, T. C. Sisman and B. Tekin, ‘‘Canonical Structure of Higher Derivative Gravity in 3D,’’ Phys. Rev. D 81, 104017 (2010) [arXiv:1002.3778]
  * (13) E. A. Bergshoeff, O. Hohm and P. K. Townsend, ‘‘More on Massive 3D Gravity,’’ Phys. Rev. D 79, 124042 (2009) [arXiv:0905.1259]
  * (14) B. Tekin, ‘‘Particle Content of Quadratic and $f(R_{\mu\nu\sigma\rho})$ Theories in $(A)dS$,’’ Phys. Rev. D 93, no.10, 101502 (2016) [arXiv:1604.00891]
  * (15) J. Gegenberg and G. Kunstatter, ‘‘Using 3-D stringy gravity to understand the Thurston conjecture,’’ Class. Quant. Grav. 21, 1197-1208 (2004) [arXiv/hep-th:0306279]
  * (16) J. Gegenberg, S. Vaidya and J. F. Vazquez-Poritz, ‘‘Thurston geometries from eleven-dimensions,’’ Class. Quant. Grav. 19, L199-L204 (2002) [arXiv/hep-th:0205276]
  * (17) E. Bergshoeff, O. Hohm, W. Merbis, A. J. Routh and P. K. Townsend, ‘‘Minimal Massive 3D Gravity,’’ Class. Quant. Grav. 31, 145008 (2014) [arXiv:1404.2867]
  * (18) J. Charyyev and N. S. Deger, ‘‘Homogeneous Solutions of Minimal Massive 3D Gravity,’’ Phys. Rev. D 96, no.2, 026024 (2017) [arXiv:1703.06871]
  * (19) A. U. O. Kisisel, O. Sarioglu and B. Tekin, ‘‘Cotton flow,’’ Class. Quant. Grav. 25, 165019 (2008) [arXiv:0803.1603]