Neoantigen-specific CD8 T cells with high structural avidity preferentially reside in and eliminate tumors

The success of cancer immunotherapy depends in part on the strength of antigen recognition by T cells. Here, we characterize the T cell receptor (TCR) functional (antigen sensitivity) and structural (monomeric pMHC-TCR off-rates) avidities of 371 CD8 T cell clones specific for neoantigens, tumor-associated antigens (TAAs) or viral antigens isolated from tumors or blood of patients and healthy donors. T cells from tumors exhibit stronger functional and structural avidity than their blood counterparts. Relative to TAA, neoantigen-specific T cells are of higher structural avidity and, consistently, are preferentially detected in tumors. Effective tumor infiltration in mice models is associated with high structural avidity and CXCR3 expression. Based on TCR biophysicochemical properties, we derive and apply an in silico model predicting TCR structural avidity and validate the enrichment in high avidity T cells in patients’ tumors. These observations indicate a direct relationship between neoantigen recognition, T cell functionality and tumor infiltration. These results delineate a rational approach to identify potent T cells for personalized cancer immunotherapy.

, while remaining amino acids in the model have negative weights as they are more frequent in low avidity structures, with the exception of G which is more frequent in the CDR3β of high avidity TCRs (Supplementary Fig. 10 and Supplementary Table S5), but more accessible to the solvent in low avidity TCRs. Supplementary Fig. 10 details amino acids frequencies in the CDR3β independently of their 3D exposure. However, in the structure-based logistic regression, we obtained significantly better models by representing the presence (1) or the absence (0) of the amino acid with a normalized solvent accessibility higher than 30% in the CDR3β of the TCR 3D model (the amino acids that have more chances to contact the peptide). Amino acid frequencies for the solvent exposed residues in CDR3β (see figure below) show some variations in amino acid composition when compared to Supplementary Fig. 10. Considering solvent accessibility of the residues, among others, we observe that G is now more frequent in low avidity structures. The decrease in G is expected since the expected role of these residues is to allow loop rearrangements (Supplementary Material 2) that help strong-interacting residues nearby (such as N, E, I, T and Y) to be displayed in front of the pMHC. As a result, G residues are hindered and exhibit decreased solvent accessibility, and consequently they are themselves exposed with lower frequency.
Three independent cross-validations were carried out: First cross-validation A first cross-validation was carried out by randomly splitting the full set into 10 sets of 38 training TCRs (3 high avidity and 7 low avidity) and 10 testing TCRs (8 high avidity and 30 low avidity). For each of the 10 combinations, we trained the model (logistic regression) from the training set and tested the predictions on the test set. Among the 10 sets, the average accuracy was 0.76±0.01. Despite the fact that the training and test sets were significantly smaller, we still observed a significant correlation, illustrating the robustness of the approach. We note that the different training and test sets include different compositions of neo-antigens, viral, TAA and different alleles, highlighting the broad applicability of the model. P < 0.0001 P < 0.0001 P < 0.0001 P < 0.0001 P < 0.0001 P < 0.0001 P < 0.0001 P < 0.0001 P < 0.0001 P < 0.0001 P < 0.0001 P < 0.0001 P < 0.0001 P < 0.0001 P < 0.0001 P < 0.0001 P < 0.0001

Second cross-validation
For the same 10 sets we re-did the cross-validation, this time including the feature selection step in the process. The following table shows the amino acids (features) that were selected for the 10  As can be seen, RNDGIL and F are the features selected in 6/10 sets and RNDL and F were automatically selected as being the most relevant features in all the 10 sets. In addition, the sign of the coefficient of these features remained the same in all 10 cross-validations (data not shown). These results support the robustness of the model and that the CDR3β amino acid representation used is the most relevant for the data under study. Among the 10 test sets, the average accuracy was 0.76±0.01.

Third cross-validation
A leave-one-epitope-out cross validation was also done, where we trained the model on all TCRs except those recognizing one given epitope, before applying the new model on TCRs recognizing this epitope. For example, we trained the data on a set excluding the viral peptide GLCTLVAML/HLA-A*02:01 and then we applied it on a set that only includes the TCRs recognizing the GLCTLVAML/HLA-A*02:01. Then, we trained the data on the set excluding the neo antigen GRKLFGTHF/HLA-B*27:05 and we applied it on a set that only includes TCRs recognizing GRKLFGTHF/HLA-B*27:05. We did it successively for the 12 epitopes under discussion. The summary of the data is given in the table below: Ten leave-one-epitope-out sets, among the 12, converged to the same features (with the same negative and positive signs) as in the full set: R, N, D, G, I, L and F. The two sets that converge to different variables are the ones leaving the EAAGIGILTV (the pink epitope in Fig. 4b) and the NILDAIAEI epitopes out. However, they still keep 4/7 and 6/7 common features, respectively, with the full set. N and I, highlighted in yellow, are the unique amino acids contributing positively to the prediction in all the 12 sets. When EAAGIGILTV (the most frequent epitope, represented in pink in Fig. 4b) is removed from the training set, we create the smallest and most different training set. Here, we removed 32% of the TCRs used in the full set, all of them with low avidity, making the situation the most different compared to the full set, and therefore the most challenging. Despite this, we still keep the features R, N, G and I in the predictor, with the same positive and negative signs. The overall accuracy leaving this epitope out is 88%, the smallest among all the sets but still satisfactory. The accuracy of the avidity predictions for TCRs recognizing this epitope, even if we did not train on this epitope and they are all of low-avidity, is 67%. The present crossvalidation shows that we have predictive power for new epitopes, even in extreme cases, as the one excluding EAAGIGILTV.
For these epitope-out cross validations, the overall average prediction is 72% for the validation sets (predictions for the validation epitope, the epitope not used in the training). The predictions range from 25% to 100%. Of note, the number of TCRs per epitope range from 1 to 12, with an average of 4 and a median of 3. So, for some cases, one single wrong prediction in the validation sets, drastically decrease the accuracy. All these three cross-validation schemes support the reliability of the predictor for the data from this study. The predictor is relevant for TAA and neoantigen and viral peptides and different HLA alleles but can be the subject of further improvements when more data will be available.

TCR-pMHC structure modeling and correlation with structural avidity
Starting from V and J segment identifiers and from the CDR3 sequences, the full sequence of the constant and variable domains of TCRα and TCRβ were reconstituted based on IMGT/GENE-DB reference sequences 1 . Homology models of the TCR-pMHC complexes were generated using the Rosetta 3.10 program 2 . Template libraries include TCR, TCR-pMHC and pMHC structures retrieved from Protein Data Bank 3 . The Rosetta "TCRmodel" protocol 4 was adapted to our approach and applied to find the respective templates and model TCR. Input α and β chains TCR amino acid sequences were parsed using regular expressions to identify CDR1, CDR2 and CDR3 loop regions and framework regions within the variable regions. Top matching CDR and framework templates were identified from the template library using the BLOSUM62 scoring function 5 . Grafting of CDR loop stem residues onto framework region was performed entailing superposition of residue backbone atoms for CDR N and C terminus onto corresponding framework residues. When one structural template matched exactly the CDR1 and CDR2 loops (germline gene template match) only the CDR3 loop was grafted. The orientation of the modeled Vα and Vβ structure was performed based on Vα/Vβ templates, while the orientation of the TCR relative to the pMHC was performed based on TCR-pMHC templates retrieved from Protein Data Bank 3 and identified using sequence similarity. Side chains and backbones of the TCR-p-MHC models were refined using the fast "relax" protocol in Rosetta 6 . A total of 500 models were produced for each TCR-pMHC. These models were subsequently ranked based on a consensus approach that combines the Rosetta energy function as implemented in Rosetta 4 and the Discrete Optimized Potential Energy as implemented in Modeller 7 . This consensus score corresponded to the sum of the normalized (Z-score) Rosetta and DOPE energies calculated over the peptide residues, as well as the CDRs and MHC residues within 6 Å from the peptide. For each TCR-pMHC, the best model according to the consensus score was selected for CDR loop refinement. The later was performed by creating 100 alternative loop conformations using the kinematic closure loop modeling of Rosetta 8 and subsequent refinement using the fast "relax" protocol. Molecular interactions were analyzed in the top5 ranked models over the 600. The final TCR-pMHC structural model is the one with the highest number of favorable interactions within the top 5 high-score models. In these structures, TCRα is chain A, TCRβ is chain B, peptide is chain C, MHC is chain D. Residue numbers start from 1 for each chain. Molecular graphics and analyses are presented for this structure making use of the UCSF Chimera package 9 . We have monitored the interactions between the residues based on the VDW radii of their atoms. The contacts between two atoms i and j are defined as the sum of their VDW radii minus the distance between them and minus an allowance for potentially hydrogen-bonded pairs overlapij = rVDWi + rVDWj -dij. We use a negative cutoff value of 0.4 Å. Correlation between the mean structural avidity of each pMHC-TCR pair with the number of non-polar, , and the number of polar, , contacts between modeled TCR and pMHC was obtained via the equation: This equation represents a simplification of the binding free energy estimation 10 , where are weighting terms applied on the number of apolar and polar contacts, respectively, and K is added to account for contributions that are not a function of the number of polar and non-polar contacts. The , parameters are fitted by multiple linear regression against the experimental pMHC-TCR T1/2. The parameters were optimized using the 10 complexes (Supplementary Table 3) and values of -62.89 s, 2.647 s and 8.747 s were obtained for K, , respectively. The correlation coefficient R, the leave-one-out correlation coefficient, the standard deviation and the p-value are 0.8679, 0.6928, 24 s and 0.005, respectively. The relevance of this correlation was assessed by a randomization test. The latter consisted in attributing randomly, for each TCR, the T1/2 of another TCR (paying attention that each of the ten T1/2 obtained experimentally was re-attributed only once), before applying the same multiple linear regression. This randomization test was performed 10'000 times. A better correlation with the randomized T1/2 than with the true T1/2 was obtained for only 0.5% of the tests, in agreement with a p-value of 0.005.

Supplementary Fig. 1. Structural avidity of neoantigen-, TAA-and virus-specific CD8 T cells. a.
Structural avidity (mean±SEM) of HLA-A*0201-restricted virus-, TAA-or neoantigen-specific CD8 T cells measured by reversible pMHC multimers (NTAmers). The number of clones is indicated in brackets (n=1 independent experiment per clone). b. Structural avidity of genetically unique TCRrestricted virus-, TAA-or neoantigen-specific CD8 T cells measured by NTAmers (mean±SEM). The number of clones is indicated in brackets (n=1 independent experiment per clone). c. Cumulative structural avidities per classes of antigen (virus, TAAs and neoantigen)-specific CD8 T cells restricted to unique TCR clonotypes described in Supplementary Table 1. The number of clones is indicated in brackets. Exact P-values are provided at 95% confidence interval and using two-sided Mann-Whitney test. Source data are provided as a Source Data file.  Fig. 3. Correlation between antigen sensitivity and structural avidity measurement and cross-validation in native clones or TCR-transduced cells. a. Correlation between antigen sensitivity and structural avidity of antigen-specific CD8 T cells. Medians of antigen specificity and structural avidities were calculated for each specificity, separating clones originating from TILs or PBLs. Spearman's coefficient R 2 and exact P-values are indicated when significant and provided at 95% confidence interval and using two-sided t-test. b-c. Paired analyses of (b) the antigen sensitivity (EC50 (M) measured by IFN-γ ELISpot assay) or (c) the structural avidity (monomeric pMHC-TCR dissociation kinetics, T1/2 (s)) of individual native TAA-and neoantigen-specific CD8 T cells and cells transduced with respective cognate TCRs (n=2 independent experiments). d-e. Associations between the antigen sensitivity (d) and the structural avidity (e) of native clones or TCRtransduced cells. TCRs are described in Supplementary Table 1. Spearman's coefficient R 2 and Pvalues are indicated when significant and provided at 95% confidence interval and using two-sided ttest. Source data are provided as a Source Data file.

. Correlations between antigen sensitivity and structural avidity of neoantigen-, TAA-and virus-specific CD8 T cells and their predicted peptide-MHC binding parameters, dissimilarity-to-self or immunogenicity. a-f. Median values of antigen sensitivity (left)
or structural avidity (right) obtained respectively by IFNγ ELISpot assay and reversible pMHC multimers on virus-(white), TAA-(red) and neoantigen-specific (green) CD8 T cells clones and their association with (a) Dissimilarity-to-self (DisToSelf) 11 PBL TIL 1 10 100 P = 0.008  Supplementary Fig. 10. Amino acids frequencies in high and low structural avidity TCRs. Amino acid frequencies of the CDR3β, excluding the first four and last three residues that may not contact the peptide, in agreement with the biophysical-based approach from Atchley 17 for TCRs with high structural avidity (orange, n=12) or low structural avidity (blue, n=39). High avidity is considered as T1/2 > 60s. Data are presented as Mean ±SD. A significance value of the differences between high and low-avidity (exact p-value) is provided at 95% confidence interval and using two-sided t-test. Source data are provided as a Source Data file.  Fig. 4d-e. The RNA coding for two KIF1BS918F-specific TCRs was transfected into recipient activated T cells followed by staining with the KIF1BS918F-specific multimer. A virus-specific TCR targeting the CMV pp65 peptide TPRVTGGGAM (HLA-B*07) was used as irrelevant control. b. The functional avidity of both KIF1BS918F-specific TCRs was measured using TCR-transfected activated T cells. Shown are the normalized relative frequencies of IFNγ-producing T cells (Mean, n=2 independent experiments) and the EC50 is given for each TCR (color-coded according to panel a). c. The reactivity of KIF1BS918F TCR-transfected T cells against the autologous tumor cell line was assessed by flow cytometry and is shown as the fraction of 4-1BB upregulating CD8 + T cells (Mean±SD). The colorcode of the different clones corresponds to that in panel a. For TCR #1, 4-1BB upregulation was assessed from n=2 wells recovered from an ELISpot and pooled in n=1 tube, for TCR #2, n=2 technical replicates were made and Mean ±SD is shown. Source data are provided as a Source Data file.     Amino acid frequencies in CDR3β (first 4 and last 3 residues removed -in agreement with the biophysical-based clustering analysis).