3 National Marrow Donor Program, Minneapolis, MN, USA

#### **Edited by:**

Michal Or-Guil, Humboldt University Berlin, Germany

#### **Reviewed by:**

Edward John Collins, The University of North Carolina at Chapel Hill, USA Fernando A. Arosa, University of Beira Interior, Portugal Dmitriy M. Chudakov, M.M. Shemyakin and Yu.A. Ovchinnikov Institute of Bioorganic Chemistry of the Russian Academy of Sciences, Russia

#### **\*Correspondence:**

Can Kesmir, Theoretical Biology and Bioinformatics, Utrecht University, Padualaan 8, 3584 CH Utrecht, Netherlands e-mail<EMAIL_ADDRESS>
Different human leukocyte antigen (HLA) haplotypes (i.e., the specific combinations of HLA-A, -B, -DR alleles inherited together from one parent) are observed in different frequencies in human populations. Some haplotypes, like HLA-A1-B8, are very frequent, reaching up to 10% in the Caucasian population, while others are very rare. Numerous studies have identified associations between HLA haplotypes and diseases, and differences in haplotype frequencies can in part be explained by these associations: the stronger the association with a severe (autoimmune) disease, the lower the expected HLA haplotype frequency.The peptide repertoires of the HLA molecules composing a haplotype can also influence the frequency of a haplotype. For example, it would seem advantageous to have HLA molecules with non-overlapping binding specificities within a haplotype, as individuals expressing such an haplotype would present a diverse set of peptides from viruses and pathogenic bacteria on the cell surface. To test this hypothesis, we collect the proteome data from a set of common viruses, and estimate the total ligand repertoire of HLA class I haplotypes (HLA-A-B) using in silico predictions. We compare the size of these repertoires to the HLA haplotype frequencies reported in the National Marrow Donor Program (NMDP). We find that in most HLA-A and HLA-B pairs have fairly distinct binding motifs, and that the observed haplotypes do not contain HLA-A and -B molecules with more distinct binding motifs than random HLA-A and HLA-B pairs. In addition, the population frequency of a haplotype is not correlated to the distinctness of its HLA-A and HLA-B peptide binding motifs. These results suggest that there is a not a strong selection pressure on the haplotype level favoring haplotypes having HLA molecules with distinct binding motifs, which would result the largest possible presented peptide repertoires in the context of infectious diseases.

**Keywords: haplotypes, HLA antigens, selection, genetic, peptide binding, bioinformatics, computational biology**

# **INTRODUCTION**

The human leukocyte antigen (HLA) genes are the most polymorphic coding loci known in humans. The HLA gene cluster is located on the major histocompatibility complex (MHC) on chromosome 6, and contains over 200 genes. The two groups of loci that contain the MHC class I and II genes dictating T cell responses are the most polymorphic. It is widely accepted that this variability is maintained by balancing selection, as individuals that are heterozygous in their HLA class I and II loci seem to have a better outcome in infections diseases [see e.g., for HIV-1 (1)]. In line with this, it has been demonstrated that in several species, and up to a certain extent in humans, females prefer to mate with males having a dissimilar HLA to increase the chance of their offspring to survive infectious diseases (2, 3). We have argued that the heterozygous advantage is on its own not enough to maintain such a large degree of polymorphism, and that the frequency dependent co-evolution with pathogens should also play a major role (4).

On the functional level the HLA polymorphism seems to be much smaller. It is well established that MHC class II molecules have largely overlapping peptide repertoires [see e.g., (5)]. In the last few years we and others showed this to be also true for class I alleles (6–9). This promiscuous peptide binding is not limited to the alleles within each locus, and extends to alleles from different MHC class I loci (7, 9). This property has important evolutionary consequences as heterozygous individuals carrying genetically different HLA molecules that nevertheless have overlapping peptide repertoires should have a diminished heterozygous advantage.

Recombinations occur frequently in the HLA region and they generate novel HLA haplotypes (10). However, the number of haplotypes found in human populations is far lower than the number of alleles observed (11), suggesting that not all the HLA haplotypes have the same chance of becoming established in human populations. Indeed, different haplotypes occur with very different frequencies in different ethnic groups/subpopulations1,2 .

1www.allelefrequencies.net

2www.nmdp.org

The "usefulness" of HLA molecules is expected to play a role in determining the haplotype frequencies: haplotypes carrying HLA molecules that are protective for certain diseases are expected to be present in relatively high frequencies in endemic areas. Indeed, several HLA haplotype disease associations have been established, e.g., in autoimmune diseases (12), squamous cell cervical cancer (13) and recurrent aphthous stomatitis (14). Along these lines, one can speculate that it should be advantageous to have MHC molecules with distinct binding specificities in a haplotype, because such haplotypes would provide an individual with more epitopes eliciting T cell responses during an infectious disease than haplotypes containing HLA molecules with similar binding motifs. This hypothesis is rather impossible to test experimentally, because one needs to determine the total set of peptides presented by a very large set of HLA molecules in cells infected by different viruses. Such experiments (e.g., eluting peptides from HLA molecules expressed by a cell) are typically done at small scales because they are time consuming. Therefore, we here study the peptide repertoires of HLA molecules that are estimated to be combined in a haplotype (i.e., the ones with a strong linkage disequilibrium) using an *in silico* approach. This approach unfortunately suffers of two main limitations. First, the quality of the peptide-HLA predictions differs between the loci. The quality of HLA-A and HLA-B peptide binding predictions is very reliable. For the other loci (HLA-C, HLA-DR, etc.) the predictions are of low quality (15), and are therefore left out of this analysis. Second, having distinct HLA binding motifs might also be important in the context of cancer or autoimmunity, however, it is very complicated to determine the role of HLA haplotypes in these diseases. Therefore, we perform our analysis explicitly for infectious diseases, and focus on the peptides that can be presented from common viruses only. Still, we think that our study is a first step to investigate the role of HLA binding motifs in the evolution of HLA haplotype frequencies.

### **MATERIALS AND METHODS**

#### **NMDP DATA**

All estimated HLA-A-B haplotype frequencies were downloaded from the National Marrow Donor Program (NMDP) database, which was established to develop and maintain a registry of HLA-typed volunteer unrelated donors for patients requiring a hematopoietic stem cell transplant (16) 3 . Four predominant US ethnic and racial groups were included in this data set: European Americans, African Americans, Asians, and Hispanic (17). Haplotype frequencies were estimated separately for each ethnic group using an implementation of the expectation maximization (EM) algorithm (18, 19).

The linkage disequilibrium, *D*, between two alleles in each haplotype was expressed as the difference between the observed and expected haplotype frequency: *D*= *fAB*−*fAfB*, where *fAB* is the observed (estimated) haplotype frequency, *fA* is the allele frequency of the HLA-A molecule in the haplotype and *fB* is the allele frequency of the HLA-B molecule in this haplotype. *D* is easy to calculate, but has the disadvantage of depending on the frequency of the alleles. In order to overcome this drawback, the normalized

3http://bioinformatics.nmdp.org

measure, D0 , was calculated a *D* 0 = *D D*max , where *D*max is the lesser of *fAB* or (1 − *fA*) (1 − *fB*) when *D* < 0 and is the lesser of *fA*(1 − *fB*) or *fB*(1 − *fA*) when *D* > 0. The advantage of this measure of disequilibrium is that it ranges between −1 and 1, regardless of the allelic frequencies in the sample. *D* 0 = 1 indicates complete LD and *D* 0 = 0 corresponds to total absence of LD.

Linkage disequilibrium statistics were calculated for each haplotype to identify the haplotypes that have a significantly positive *D* 0 .

#### **HLA LIGAND PREDICTION**

To be able to perform our analysis for as many as possible HLA molecules, we used NetMHCpan (15) to predict peptide-HLA binding affinity. NetMHCpan assigns to each peptide-HLA pair a predicted IC50 value, indicative of the predicted binding affinity. To assess whether a peptide binds to an HLA molecule depends on the choice of binding threshold, and the optimal threshold has been discussed (20). If one assumes that all HLA molecules use a fixed threshold, one can use the default threshold of 500 nM (21, 22), otherwise a 5000 nM threshold can be used to allow for the comparison of more weakly binding peptides. However, using a fixed threshold to define predicted binders result in large differences in the predicted repertoire sizes between HLA molecules. For instance using a fixed threshold of 500 nM, the HLA repertoire sizes range between 20 and 6574 peptides for the viral set listed in Table S2 in Supplementary Material. As such a variance could introduce large biases in our analysis, we defined the 1% top-ranking peptides as candidate binders for each HLA molecule. This gives each HLA molecule the same ligand repertoire size (i.e., 570 binders per HLA molecule for the viral set listed in Table S2 in Supplementary Material). To check the consistency of our results with respect to these parameters, we repeated every analysis with the fixed threshold of 500 nM. All results presented below were derived using the threshold of 1% to define candidate binders, and remain similar for a fixed threshold, unless mentioned otherwise. To test whether our results change if one were to use a much larger data set, we also generated predicted binders using a much larger set of viruses (see below), in which case each HLA molecule had 60-fold more binders.

#### **VIRAL DATA**

The proteomes of 17 common human viruses were downloaded from the European Bioinformatics Institute website4 (downloads were made in October 2006, listed in Table S2 in Supplementary Material) as the source of potential HLA ligands. To extend this data set, we downloaded another set of proteomes (downloads were made in October 2008)4 from viruses that are known to infect mammalians (*n* = 904). We used the HLA-peptide binding predictors (see above) to screen all possible unique virus-derived 9-mer peptides.

#### **PEPTIDE REPERTOIRE OVERLAP**

We define the peptide repertoire overlap between two HLA alleles in the same HLA-A-B haplotype, *Fp*, as the fraction of overlapping ligands between these two HLA class I molecules among

4www.ebi.ac.uk

all of their ligands: *Fp* = *A* T *B A* S *B* , where *A* and *B* are the ligand sets for HLA-A and -B molecule, respectively. *Fp* = 1 implies that the HLA-A and HLA-B molecules belonging to the same haplotype have the same epitope repertories while *Fp* = 0 indicates completely different peptide repertoires for two HLA molecules.

#### **RESULTS AND DISCUSSION**

#### **DETERMINING HLA-A-B HAPLOTYPES**

The "true" HLA haplotypes can be determined either by molecular haplotyping or family-based segregation studies (23–25). However, both approaches are expensive and laborious, and therefore, statistical methods are typically used to infer haplotypes from datasets covering large population of individuals with known HLA genotypes (26). Several methods have been proposed to infer HLA haplotypes from genotype data, and in recent studies the performance of two most commonly used approaches, EM algorithm based (implemented in Arlequin V3.0), and the Bayesian algorithm based (implemented in PHASEV2.1.1),have been compared (27, 28). Unfortunately, neither of the methods could infer all of the known haplotypes: incorrect haplotypes were estimated in more than 30% of the cases. However, once the sample size increases, the power of these statistical methods is expected to increase tremendously.

National marrow donor program2 provides, to our knowledge, the largest repository of HLA-typed donors. Here use of statistical methods should become more reliable (16): for the HLA-A-B haplotype, the total chromosome counts (2N) for the four major ethnic groups exceeds 2000. On the NMDP webpage3 , the high-resolution allele and haplotype frequencies [estimated by EM method, (18, 19)] are available (17). Focusing on HLA-A-B haplotypes, the most common haplotypes found in US population (separated into four main ethnic groups) are summarized in **Table 1** (adopted from bioinformatics.nmdp.org, December 2007 version). Alternatively Allele frequencies web server1 , provides allele frequencies established in smaller, but probably better defined studies (29).

In the NMDP database 660 possible haplotypes are reported for European Americans. However, many of these haplotypes are bound to be falsely predicted, e.g., due to the limited number of individuals carrying particular combinations of specific HLA molecules. To decrease the amount of wrongly identified haplotypes in our analysis, we apply a rather strict criterion for considering a predicted haplotype as a "true" haplotype: we demand a positive LD value that is significantly different than zero (*p* < 0.01, see Materials and Methods). This criterion decreases the number of haplotypes to 60 for European Americans. These 60 haplotypes are estimated to cover 58% of the population (see **Table 2**), i.e., current statistical methods and data sets (even the large repositories like NMDP) remain rather limited in providing the HLA haplotype diversity of a population. For other ethnical groups, the number of reliable haplotypes drops to 30–40 per ethnic group, even though the number of possible haplotypes was the same or higher (see **Table 2** and footnote text 3). The population coverage in non-European groups was lower, 30–40% of their respective populations, possibly due to the lower number of individuals with known HLA-typing. All together we detected 120 reliable unique haplotypes by summing over these ethnical groups (the

**Table 1 | Occurrences of the three most common HLA-A-B frequency ranked haplotypes in four major ethnic groups in US (adopted from bioinformatics.nmdp.org)**.

EUR, Caucasian; AFA, African; API, Asian; HIS, Hispanic. F stands for population frequency in percentages.

**Table 2 | Numbers of different haplotypes with a significantly positive LD in four major US ethnic groups**.

Population coverage in percentages is given within parenthesis.

full list of haplotypes can be found in Table S1 in Supplementary Material).

#### **PEPTIDE REPERTOIRE OF AN HAPLOTYPE**

Having identified the HLA-A-B haplotypes for the US population, we next estimated the overlaps between peptide repertoires of HLA-A and -B molecules that belong to the same haplotype.

We used an *in silico* approach and predicted the peptide repertoire of all HLA-A and HLA-B alleles that are part of the 120 predicted haplotypes, using the proteomes of common viruses (see Table S2 in Supplementary Material) and HLA-peptide binding predictor NetMHCpan (15, 30) (see Materials and Methods). NetMHCpan is the only prediction system available right now that can reliably predict the peptide binding affinities for the large set of HLA-A and HLA-B molecules we are taking into account in this study. The analysis of the 120 significant haplotypes demands predictions for 39 HLA-A and 63 HLA-B molecules. This predictor assigns an IC50 value to each peptide-HLA pair, which can be used as a predicted binding affinity. Using the widely accepted IC50 value of 500 nM as a threshold to distinguish binders from non-binders, generated a large variation in the predicted repertoire sizes of different HLA molecules (20–6574 peptides for the viral set listed in Table S2 in Supplementary Material), which could strongly bias our results. As the physiologically relevant IC50 values are difficult to estimate for each HLA molecule, we have chosen a simplified approach and define the peptide repertoire as the top 1% peptides with the highest HLA binding affinities for each HLA molecule. This approach removes the potential bias introduced by different repertoire sizes, but ignores the fact that some HLA molecules can be more specific than others, and therefore present much fewer peptides.

The unique haplotypes listed in Table S1 in Supplementary Material have an average peptide overlap of 1.1% (with a standard deviation of 1.8%, and a median of 0.44%) using top 1% threshold. The distribution of the overlaps is given in **Figure 1A**, and varies somewhat among the haplotypes. Out of 120 haplotypes, 36 (30%) have non-overlapping peptide repertoires, at least for the common viruses that we tested (see Table S1 in Supplementary Material). Only for five haplotypes is the repertoire overlap higher than 5%, and HLA-A0101-B1517 with an overlap of 11.8% is the highest. HLA-A0101-B1517 is a rare haplotype and occurs only in Asian Americans with a population frequency of 0.4%. To test whether or not rare haplotypes tend to have larger overlaps than common haplotypes, we weighted the overlaps found in **Figure 1A** with the population frequency of the HLA-A-B haplotypes (**Figure 1B**). Since the weighted overlaps remain very similar to the unweighted overlaps (**Figures 1A,B**), there is no evidence for a trend of rare haplotypes having the largest overlaps. In line with this, the frequency of a haplotype is only weakly correlated with the degree of peptide overlap between the haplotype's HLA-A and HLA-B molecules (*r* = −0.08, *p* = 0.4, Spearman rank correlation). Apparently, there is no selection pressure increasing the frequency of the haplotypes with a small peptide repertoire overlap.

These results were obtained using the top 1% peptides with the highest HLA binding affinities as the set of presented peptides per molecule. Using the set of common viruses listed in Table S2 in Supplementary Material, this threshold results in approximately 570 predicted binders per HLA molecule. To test whether the results presented in **Figure 1A** would be sensitive to the number of peptides, we collected proteomes for a much larger set of mammalian viruses (*n* = 904), which contains approximately three million unique peptides of nine amino acids. Using the same 1% threshold for this larger set, we predict for each HLA molecule the extended peptide repertoire and calculate the overlaps as explained before. The distribution of the overlaps hardly changes despite the fact that we enlarged the presented peptide repertoire 60-fold per molecule (see Figure S1 in Supplementary Material). As the results presented in **Figure 1A** seem fairly insensitive to the number of peptides used, we perform the rest of the analysis on the small data set of common viruses.

To test whether or not HLA binding motifs affect haplotype compositions requires comparison of the peptide overlaps of "true" haplotypes with those of "random" haplotypes. To do this, we reshuffled HLA-A and HLA-B molecules in the 120 "true" haplotypes to calculate an expected peptide repertoire overlap for randomly made haplotypes. Although the HLA molecules in random haplotypes can have overlaps up to 28% in their peptide repertoires (see **Figure 1C** and results not shown), the distribution of the overlaps is not significantly different from the distribution given in **Figure 1A**. The set of random haplotypes was generated 100 times, and in none of these cases was the distribution

significantly different from the one given in **Figure 1A** (using a Kolmogorov–Smirnov test). Finally, we calculated a weighted peptide repertoire overlap for the random haplotypes by assuming that the frequency of a random haplotype is simply the multiplication of the frequency of their HLA-A and HLA-B alleles (i.e., assuming a complete lack of linkage disequilibrium). Again, the distribution of weighted overlap of random haplotypes (**Figure 1D**) is not different from that of the real haplotypes (**Figure 1B**). Taken together, these results suggest that HLA-A and -B in pairs in general have distinct peptide binding preferences, and that a small overlap is not a unique property of the HLA molecules having a strong linkage disequilibrium.

The overlap distribution presented in **Figure 1** seems to be in contradiction with earlier results, which estimated cross loci peptide overlaps of 23–44% (7,9). However, this overlap was estimated at the population level, i.e., these percentages reflect the fraction of the peptide repertoire of an HLA-A molecule which is also expected to be presented by at least one HLA-B molecule in the population (and vice versa). Within an individual having maximally two different HLA-A and HLA-B molecules, the overlaps should remain much lower than the population based overlaps. In addition, one needs to realize that the low overlaps presented in **Figure 1** depend on the threshold used to define the peptide repertoire of an HLA molecule. When we use a higher threshold of 2 or 5% the average overlap increases to 2.3 and 5.9%, respectively (with the standard threshold of 1%, the average overlap was 1.1%, see above). However, the choice of the threshold hardly affects our main result, namely that the HLA-A and -B pairs that are in a linkage disequilibrium do not have more distinct binding motifs than random HLA-A/B pairs (**Figure 1**).

# **CONCLUSION**

We hypothesized that it should be advantageous to have HLA molecules with distinct binding specificities combined in a haplotype, because during a viral infection such haplotypes would give an individual a larger epitope repertoire than haplotypes containing HLA molecules with similar binding motifs during a viral infection. To test this hypothesis, we used the high-resolution data available in the NMDP database on haplotype frequencies, and employed state of the art peptide-HLA binding prediction tools. We find that for all the haplotypes we could reliably identify in the US population, their HLA-A and HLA-B molecules present largely distinct set of peptides (**Figure 1A**). However, this turned out to be a generic property of HLA-A and HLA-B molecules: when we compared random HLA-A and -B pairs we find a very similar distribution of the presented peptide overlaps (**Figure 1C**). Moreover, there is no evidence for selection as there is no correlation between the population frequency of the HLA-A-B haplotypes and the overlap in the peptide repertoires of their HLA-A and HLA-B molecules. Taken together, these results suggest the complementarity of binding motifs is a general property of HLA-A and HLA-B molecules, and that complementarity does not affect the HLA haplotype composition. We were not able to specifically test the effect of complimentary binding motifs in the context of autoimmunity and cancer, as for both cases it remains unclear which set of human proteins should be taken as possible auto antigens. Complimentary binding motifs are expected to increase the

number of potential self antigens, which could increase the risk of autoimmunity. Finally, the frequency of an HLA haplotype is determined by complex interactions with many different factors, one example is the correlation between birth weight and particular haplotypes (31). Our results suggest that complimentary binding motifs of HLA molecules during viral infections play a minor role, if any, compared to these other factors in the evolution of HLA haplotype frequencies.

# **ACKNOWLEDGMENTS**

This study was financially supported by the University of Utrecht through a High Potential grant. The funders had no role in study design, data collection, and analysis, decision to publish or preparation of the manuscript.

# **SUPPLEMENTARY MATERIAL**

The Supplementary Material for this article can be found online at http://www.frontiersin.org/Journal/10.3389/fimmu.2013.00374/ abstract

# **REFERENCES**

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 23 July 2013; paper pending published: 28 August 2013; accepted: 31 October 2013; published online: 14 November 2013.*

*Citation: Rao X, De Boer RJ, van Baarle D, Maiers M and Kesmir C (2013) Complementarity of binding motifs is a general property of HLA-A and HLA-B molecules and does not seem to effect HLA haplotype composition. Front. Immunol. 4:374. doi: 10.3389/fimmu.2013.00374*

*This article was submitted to T Cell Biology, a section of the journal Frontiers in Immunology.*

*Copyright © 2013 Rao, De Boer, van Baarle, Maiers and Kesmir. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Receptor pre-clustering and T cell responses: insights into molecular mechanisms

#### **Mario Castro1\*, Hisse M. van Santen2\*, María Férez 2 , Balbino Alarcón2 , Grant Lythe3 and Carmen Molina-París 3**

1 Grupo de Dinámica No-Lineal and Grupo Interdisciplinar de Sistemas Complejos (GISC), Escuela Técnica Superior de Ingeniería (ICAI), Universidad Pontificia Comillas, Madrid, Spain

2 Departamento de Biología Celular e Inmunología, Centro de Biología Molecular Severo Ochoa, Consejo Superior de Investigaciones Científicas, Universidad Autónoma de Madrid, Madrid, Spain

3 Department of Applied Mathematics, School of Mathematics, University of Leeds, Leeds, UK

#### **Edited by:**

Rob J. De Boer, Utrecht University, Netherlands

#### **Reviewed by:**

Daniel Coombs, University of British Columbia, Canada Barbara Szomolay, University of Warwick, UK

#### **\*Correspondence:**

Mario Castro, Universidad Pontificia Comillas, C/ Alberto Aguilera 25, Madrid E28015, Spain e-mail<EMAIL_ADDRESS>Hisse M. van Santen, Centro Biología Molecular Severo Ochoa, Calle Nicolás Cabrera 1, Campus de Cantoblanco, Madrid 28049, Spain e-mail<EMAIL_ADDRESS>
T cell activation, initiated byT cell receptor (TCR) mediated recognition of pathogen-derived peptides presented by major histocompatibility complex class I or II molecules (pMHC), shows exquisite specificity and sensitivity, even though the TCR–pMHC binding interaction is of low affinity. Recent experimental work suggests thatTCR pre-clustering may be a mechanism via whichT cells can achieve such high sensitivity.The unresolved stoichiometry of theTCR makesTCR–pMHC binding andTCR triggering, an open question. We formulate a mathematical model to characterize the pre-clustering of T cell receptors (TCRs) on the surface of T cells, motivated by the experimentally observed distribution of TCR clusters on the surface of naive and memory T cells. We extend a recently introduced stochastic criterion to compute the timescales of T cell responses, assuming that ligand-induced cross-linked TCR is the minimum signaling unit. We derive an approximate formula for the mean time to signal initiation. Our results show that pre-clustering reduces the mean activation time. However, additional mechanisms favoring the existence of clusters are required to explain the difference between naive and memory T cell responses. We discuss the biological implications of our results, and both the compatibility and complementarity of our approach with other existing mathematical models.

**Keywords: T cell receptor, clustering, stochastic dynamics, signaling, naive T cells, memory T cells**

# **1. INTRODUCTION**

A hallmark of the adaptive immune system is the ability of T cells, making use of the T cell receptors (TCRs) on their surface, to recognize a given agonist peptide–MHC ligand complex (pMHC) with high sensitivity (1). Some aspects of TCR–pMHC molecular interactions that are of current research interest are the frequency of encounters between T cells and the agonist pMHC, how cell– cell interactions determine the activation of lymphocytes (2), how early interactions change the state of the T cell receptor (3), what are the mechanisms of modulation of receptor–ligand interactions at cell–cell interfaces (4), and how protein organization in the cell membrane (for instance, protein islands or lipid rafts) affects the recognition process (5). Some recent experiments have explored the role of dimensionality on T cell activation and have highlighted the significance of the events taking place at the receptor level [see Refs. (1) and (6) for comprehensive reviews].

These open questions have been addressed with the use of mathematical modeling. Different theories can be classified according to the level of description (7). At the individual TCR– pMHC bond level, the kinetic proof-reading model (8) assumes that the TCR needs to undergo a series of consecutive (phosphorylation) steps before being triggered. Also at the TCR level, the optimal dwell time model (9) reconciles the concurrence of different timescales, providing an *optimal* timescale between the

very short times related to the off rate of TCR–pMHC binding, and the long times related to kinetic proof-reading mechanisms. The TCR occupancy model (10) considers the cell as a *counting device* in which multiple TCR–pMHC interactions are required to activate a T cell. In a similar fashion, the serial triggering model (11) proposed that the same pMHC can engage *serially* different TCRs. This model enriches the viewpoint of the TCR occupancy model, by giving greater relevance to the role of the pMHC itself. Finally, the serial encounter model (12) and the confinement time model (13) combine several of the ideas above and provide some appealing explanations by relaxing some restrictions in those models.

While antigen presenting cells (APCs), such as dendritic cells or B cells, present 103–104 times more self-pMHC than antigenic pMHC, self-pMHC ligands by themselves do not usually elicit a T cell response, even though their affinity for TCRαβ is only 10 times lower than the affinity of the antigenic pMHC (14). This illustrates how a small difference in affinity results in high specificity, when there is only a few antigenic pMHC molecules in a background of self-pMHC ligands (15).

The T cell signaling process begins with (extracellular) TCR– pMHC binding, followed by phosphorylation of the intracellular ITAM domains of the TCR–CD3 complex. When a TCR binds a pMHC molecule, the TCRαβ hetero-dimer binds the peptide, while the CD4 or CD8 co-receptor binds the MHC molecule. The binding of the co-receptor activates the tyrosine kinase Lck, which phosphorylates the ITAMs of the CD3 complex. ITAM phosphorylation allows recruitment of intracellular signaling components that mediate downstream signaling events (16).

It has recently been suggested that, contrary to what happens in TCR micro-clusters and the immunological synapse, clustering is not only induced by the ligand but by an *avidity maturation* mechanism (or pre-clustering) (17), allowing the aggregation of chains of TCRs as long as 20 units (around 200 nm long), and referred to as *nano-clusters* (3, 18). Specifically, multimeric TCR–CD3 complexes are activated at low agonistic pMHC concentrations and monomeric TCRs remain unaffected at low ligand concentration. The TCR nano-clusters could enhance T cell sensitivity by the mechanisms proposed in the models of T cell activation (7), as their existence would reduce the time needed for two (or more) receptors to aggregate (by diffusion). This pre-cluster formation could be explained by three different mechanisms (3):

The existence of these nano-clusters does not exclude additional mechanisms of T cell activation, as long as they involve the *cooperation* of receptors when they aggregate. Thus, while models such as kinetic proof-reading [and improvements as described in Ref. (22)] operate at the level of a single receptor, other models might be used in combination with the fact that the pre-cluster distribution of naive and memory T cells is different.

Additionally, the fact that the TCR stoichiometry has not been resolved under physiological conditions, yet, makes it even more difficult to understand, at a molecular level, the dynamics of TCR pre-clustering (23). TCR pre-clustering could be an example of a more general mechanism of membrane-bound molecular preclustering, as clustering prior to cell–cell interaction has also been observed on the surface of APCs (19–21). It is worth mentioning that monomeric TCRs can still be activated at increasing ligand concentrations, thus, conferring the T cell with a capacity to generate a dose-dependent response at very high pMHC doses, when multimeric TCR–CD3 complexes are already saturated (18). Such mechanisms have been previously described for chemotactic bacteria, as a cellular mechanism to control sensitivity (24).

Various mechanisms have already been suggested, at the population, cellular or molecular level, to explain the capacity of T cells to respond, faster and more strongly, to a second antigenic encounter. However, the underlying mechanisms of the observed changes in the sensitivity of the T cell for pMHC ligand-mediated TCR stimulation (25) have not yet been clearly elucidated. Interestingly, the distribution of clusters in naive and memory T cells is different: memory T cells accommodate larger linear TCR clusters than naive ones. This could explain why memory T cells elicit more rapid responses than naive T cells (17) (see **Figure 1** below).

In this paper, we explore the consequences of TCR preclustering in signaling and in distinguishing naive from memory T cell responses. We present some experimentally obtained distributions of TCR clusters for both types of cells (see **Figure 1**), and two complementary theoretical models: (i) a simple model of receptor oligomerization that describes cluster size distributions, and (ii) a generalization of the stochastic T cell response criterion of Ref. (26), to accommodate the hypothesis that the minimum signaling unit is composed of a TCR receptor cluster that is bound by the same cross-linked multivalent ligand. We find that this signaling unit is able to discriminate between agonist and antagonist pMHC ligands (with greater sensitivity than in the monomeric case), and to explain some of the advantages that higher cluster sizes can provide to memory T cells. The model also points at the need to invoke additional cooperativity mechanisms, to explain the experimentally observed role of clustering in T cell responses (27). Finally, this model of ligand-induced TCR cross-linking can be relevant in physiological conditions, according to the defective ribosomal products (DRiP) hypothesis (28,29), which provides a rapid source of peptide precursors to optimize immuno-surveillance of pathogens and tumors (30).

# **2. MATHEMATICAL MODELING OF TCR PRE-CLUSTERING AND T CELL ACTIVATION**

#### **2.1. MODEL 1: T CELL RECEPTOR PRE-CLUSTERING**

The TCR–CD3 complex consists of the pMHC binding TCRαβ hetero-dimer, associated with the hetero-dimers CD3γ and CD3δ, and the homo-dimer CD3ζ ζ . Binding of a stimulating pMHC ligand by the extracellular domain of TCRαβ results in conformational changes in the intracellular part of the CD3 chain, and phosphorylation of the immuno-receptor tyrosinebased activation motifs (ITAMs) in the intracellular domains of the CD3γ , CD3δ, and CD3ζ ζ dimers, which in turn lead to initiation of downstream signaling cascades and T cell activation.

It has long been recognized that the TCR–CD3 complex forms clusters upon ligand binding (31–36). More recently, it has been shown that in the absence of stimulating pMHC ligand, TCR– CD3 complexes are already expressed at the cell surface as a combination of monomeric and oligomeric TCR complexes or TCR nano-clusters (18). Electron microscopy (EM) analysis of immuno-gold-labeled human and murine T cells showed that these nano-clusters consist of up to 20 TCR–CD3 complexes. The exact stoichiometry of the nano-clusters has not been resolved yet.

The integrity of TCR nano-clusters depends on cholesterol present at the cell surface membrane (18). The formation of the clusters depends, at least, on the trans-membrane region of the CD3ζ ζ homo-dimer (17), perhaps due to the capacity of ζ ζ dimers to form dimers of dimers (37). Other possible mechanisms of cluster formation rely on the capacity of the extracellular domain of TCRα to dimerize (38).

CD3-specific mAb 2C11 and 10 nm gold-conjugated protein-A. Cell surface replicas of the labeled T cells were analyzed by transmission electron microscopy and the number and size of the observed gold clusters were recorded. **(A)** TEM image of surface replicas of a memory and a naive OT-1 T cell. The insets to the right show an enlargement of the boxed areas. **(B)** Quantification (mean ± SD) of gold particles in clusters of the indicated

distribution of clusters of three or more gold particles and statistical analysis (2-tailed Student's t-test: \*p < 0.05, \*\*p < 0.01, and \*\*\*p < 0.001). All naive and memory T cells had clusters with gold. However, whereas in naive T cells the maximum gold cluster size shared by all cells was four, this was eight for memory T cells. Also clusters bigger than twenty three particles were present in four out of five memory T cells, and only two out of seven naive T cells.

This body of experimental evidence allows us to conclude that multimeric TCR–CD3 complexes are co-expressed with TCR monomers on the surface of resting T cells.

A simple model of aggregation of TCRαβ units is depicted in the left panel of **Figure 2**. Given a chain of length *n* (with *n* heterodimers linked), in a small time interval 1*t*, with probability *q*+ 1*t*, the chain increases to length *n* + 1, and with probability *q*− 1*t*, the chain decreases to length *n* − 1. Thus, by probability conservation, the probability to remain the same length *n* is 1 − (*q*+ + *q*−)1*t*.

Mathematically, the dynamics of the process can be described by a continuous time Markov chain (39) (or birth and death process, as we assume that polymerization takes place in unit steps). The state space is {1, 2, 3, . . ., *n* − 1, *n*, *n* + 1, . . .}, where the number denotes the number of TCRs in a cluster:

$$1 \stackrel{q-}{\underset{q\_+}{\rightleftharpoons}} 2 \stackrel{q-}{\underset{q\_+}{\rightleftharpoons}} 3 \stackrel{q-}{\underset{q\_+}{\rightleftharpoons}} \cdots \stackrel{q-}{\underset{q\_+}{\rightleftharpoons}} n-1 \stackrel{q-}{\underset{q\_+}{\rightleftharpoons}} n \stackrel{q-}{\underset{q\_+}{\rightleftharpoons}} n+1 \cdots \cdots$$

The forward Kolmogorov equations for the probability of having a cluster of size *n* are given by (40)

$$\begin{aligned} \frac{dp\_n(t)}{dt} &= q\_+ p\_{n-1}(t) + q\_- p\_{n+1}(t) - (q\_+ + q\_-) p\_n(t), \; \forall n \ge 2, \\\frac{dp\_1(t)}{dt} &= q\_- p\_{2(t)} - q\_+ p\_1(t). \end{aligned}$$

The stationary probability distribution is then given by

$$\lim\_{t \to -\infty} p\_n(t) \equiv \pi\_n = \frac{b^{n-1}(1-b)}{(1-b^{N\_{\text{max}}})},$$

$$b < 1, \quad n \in \{1, 2, 3, \dots, N\_{\text{max}}\},\tag{1}$$

with *b* = *q*+ *q*− , and π*n* the probability (in thermodynamic equilibrium) to have a cluster of size *n*. When *b* < 1 (the number of clusters with a given size, *n*, decreases as *n* increases), and taking into account that peripheral T cells have around *N*max ' 3 × 104 receptors, the latter expression can be further simplified to

$$
\pi\_n = b^{n-1}(1-b)\,, \quad b < 1\,, \quad n \in \{1, 2, 3, \ldots\}\,. \tag{2}
$$

#### **2.2. MODEL 2: A BIVALENT MODEL FOR T CELL ACTIVATION**

The TCR–pMHC binding model introduced in Ref. (26) considered monovalent pMHC ligands binding to TCR monomers on the surface of a T cell. Monovalent ligands have been reported to elicit a T cell response (41–43), but only when they are immobilized on a surface (which makes it difficult to assess whether they are truly monovalent or not). Yet, multivalent receptor–ligand interactions are required to elicit T cell responses in both CD4+ and CD8+ T cells. In what follows, and supported by a body of experimental work (3, 24, 44), we adopt the hypothesis that the minimum activating unit is a TCR–pMHC cross-linked dimeric complex (31, 45–47). We make use of the binding model (Model 2) with pMHC dimers (ligands) and dimeric TCRs (receptors), described in the right panel of **Figure 2**.

Gold-labeling experiments support the existence of nanoclusters with more than two TCRs, yet it can be shown (see Section 5.2) that the key parameter of the mathematical model is the fraction of monomeric to multimeric TCR clusters. Thus, without loss of generality, we will assume that all TCR clusters are dimeric.

The biochemical reactions encoded by the right panel of **Figure 2** are as follows:

oligomerization model not mediated by ligand (Model 1). We assume that receptors are able to diffuse and aggregate to an existing cluster. However, we exclude the possibility of clusters with size larger than one to diffuse. Clusters grow one monomeric unit at a time. Right panel: reactions included in the stochastic activation model (Model 2). Ligands in solution are able to

minimum dwell time for a TCR–pMHC complex is the time the complex must remain bound in order to reach a level of ITAM phosphorylation, which will allow TCR triggering. Any binding, which persists for longer than the minimum dwell time is classified as a productive binding [see Refs. (48) and (26) for further details].

• From an immunological perspective, the relevant parameter is the *mean time to signal initiation*, or MTSI (26). Namely, the MTSI is the average time needed for a T cell response according to the criterion that at least *N* TCR dimers should be bivalently engaged to a bivalent ligand (pMHC) for at least a time τ .

Here we assume that *N* is around 10–100. That is, 10–100 TCRs are required for signaling and *NB* = *b* × *NR* is of the order of 104 , with *NB* the total number of clusters on the T cell surface. This means, under the assumptions of Model 2, that at most, there can be *N* = 100 internalization events, as this is the number of triggered TCRs. Thus, in this approximation, the loss of TCR due to internalization after triggering can be safely neglected. Nevertheless, internalization is an important step in early signaling, and a proper mechanistic model to justify the value of τ will require internalization to be considered. This analysis is out of the scope of this article.

We implement these reactions as a Markov process, and solve them numerically using the standard Gillespie algorithm (49), and with the parameters summarized in **Table 1**. We have made use of three different ligands: 4A, 4P, and 4N, which were also used in Ref. (26). For these ligands, that bind the same TCR with different affinities, a simple estimation of the number of cross-linking events required to elicit a T cell response is summarized in **Table 2**.

There is some evidence that,under physiological conditions,the chance of two specific peptides being presented by two MHC molecules in sufficient proximity and long enough to act as a dimer is very small (46). This will make ligand-induced TCR cross-linking a rare event. However, some recent experimental work on the distribution of cognate pMHC molecules on the surface of APCs shows that both for MHC class I (virus infection models), and for MHC class II (antigen uptake via the endocytic route) clusters of cognate pMHC can be detected (19–21).

We also note that ligand concentration is not the only factor that depends on physiological conditions. According to the DRiP hypothesis (28, 29), rapid viral antigen presentation is possible because antigenic peptides originate from defective ribosomal products that have short half-lives. Although this phenomenon affects the time between viral challenge and antigen presentation, we assume it is independent of the subsequent signaling dynamics of T cell activation.

#### **3. RESULTS**

#### **3.1. DISTRIBUTION OF TCR CLUSTERS**

The mathematical model described in Section 2.1, or Model 1, allows us to obtain the value of *b* that best fits the experimental data. We have used a weighted (by the variance) minimum-square regression to fit the experimental distributions to equation (2). This kind of fit minimizes the value of χ 2 . Thus, in **Figure 3**, we show the agreement between theory and experiment, with values: *b*naive = 0.32 and *b*memory = 0.55. The difference between *b*naive

#### **Table 1 | Summary of the parameters used in the stochastic simulations.**

For typical values of the dissociation rate, kd, we find that k2 is about 10–50 times koff. We have assumed k−2 = koff following Ref. (44). When not explicitly shown, we have used the same parameters as in Ref. (26).

aThe cross-linking rate k2 is adapted from Ref. (44) for bivalent receptors.

**Table 2 | Estimated mean number of cross-linking events,** *N* 0 ' *N e***2***k*−**2**τ **, required to elicit a T cell response (SP thymocytes).**

and *b*memory can be explained by the existence of larger (or at least more localized) lipid rafts on the membrane of memory T cells (50, 51). Thus, the rates *q*± could be the effective combination of two mechanisms: one related to the diffusion of receptors on the membrane, and the other related to the aggregation of the receptors at the molecular level. The presence of cholesterol on the membrane changes the diffusion coefficient of the TCR receptors, as receptor diffusion within the raft is inhibited due to protein anchorage (52) and, thus, stabilizes the formation of clusters (a larger value of *b* means that, once two receptors are embedded in the same lipid raft, it is more difficult for them to become separated from each other).

A consequence of Model 1 is that, as the stationary probabilities need to sum up to one, the fraction of clusters of size larger than one is, precisely, *b*. This fraction is 72% higher for memory T cells than for naive T cells: *b*memory/*b*naive = 1.72.

#### **3.2. MEAN TIME TO SIGNAL INITIATION**

In **Figures 4A–D**, we show how the stochastic criterion is able to provide a ligand hierarchy according to their *potency*. Namely,

**FIGURE 4 | Dependence of the mean time to signal initiation (MTSI),T (N,** τ **) to have N cross-linked ligand–receptor complexes bound for at least a dwell time** τ **for different model parameters as shown in every panel**. The results have been obtained by making use of a Gillespie algorithm, after averaging over 100 realizations for each

fitted to equation (2) using a weighted (by the variance) minimum-square regression. The fitted values are bnaive = 0.32 and bmemory = 0.55.

set of the parameters, summarized in**Table 2** (a python code for the stochastic integration is available upon request). Units of time are seconds. All parameters are taken from**Tables 1** and **2** except **(A)** N = 10, τ = 1 and k 2 = 10 × k off; **(B)** N = 10, τ = 1 and k 2 = 40 × k off; **(C)** N = 10, τ = 4 and k 2 = 10 × k off; **(D)** N = 10, τ = 4 and k 2 = 40 × k off.

the most agonistic ligand, 4P, elicits a T cell response in times of the order of a few seconds in all cases. On the contrary, the most antagonistic ligand, 4N, takes extremely large times to do so (in practical terms, this means it does not elicit a T cell response). Thus, TCR clustering can enhance the *potency* of ligands, when compared to the monomeric case (26), as experimentally observed and theoretically shown.

Following a similar approach to that of Ref. (26), we can derive an approximate formula for the mean time to signal initiation (MTSI), *T*(*N*, τ ), for different ranges of ligand concentration, ρ. We write (see **Figure 5A** and Section 5.3 for further details):

$$T(N,\mathbf{r}) \simeq \begin{cases} \mathbf{r} & \text{at high concentration} \\ \mathbf{r} + \left[ \frac{N \exp(2\mathbf{k}\_{-2}\mathbf{r})}{2\rho N\_{\rm B} \mathbf{k}\_{\rm bak}} \right]^{1/2} & \text{at intermediate} \\ \mathbf{r} + \frac{N \exp(2\mathbf{k}\_{-2}\mathbf{r})}{4\rho N\_{\rm B} \mathbf{k}\_{\rm on} (\mathbf{k}\_{2}/\mathbf{k}\_{-2})} & \text{at low concentration} \end{cases} (3)$$

These three regimes correspond to different immunological scenarios. In the case of high concentration of ligand, ligand is in great excess, so that the required number of *signaling units* is reached, almost as soon as the first signaling unit is formed (time of order τ ). At low ligand concentration, the dynamics is limited by the first binding event, as cross-linking occurs in a slower timescale. So, the MTSI has the same functional form as that for the monomeric case (26). Finally, for intermediate ligand concentration, the competition between binding and cross-linking implies a more complicated mathematical relationship. Of greater relevance to the discussion is the nature of the ligand (with different *k*on and *k*off rates), and the number of TCR clusters on the membrane of the T cell (encoded in the parameter *NB* = *b* × *NR*, with *NR*, the average number of TCRs per T cell, see **Table 1**).

An expression for the variance of the time to signal initiation (TSI) cannot be provided in a closed form [as done in Ref. (26)]. However, the fact that the variance decreases as the ligand concentration increases, suggests that the mathematical formula for the variance in the monovalent case can provide an upper bound to the present (dimeric) case.

Using equation (3),we also can deduce the role of pre-clustering in the signaling time, or MTSI. As the number of bivalent clusters is *b* × *NR*, the larger *b* is, the shorter the response time becomes. The model predicts that, for physiological conditions (not too high ligand concentrations), the ratio of the MTSI for naive and memory T cells is inversely proportional to the ratio of their corresponding values of *b*. Namely, memory cells would respond up to 72% faster than naive ones (**Figure 5B**).

# **4. DISCUSSION**

TCR triggering mechanisms are currently under debate [see, for example, Ref. (53) and (7) for recent reviews]. TCR clustering may be invoked as a description of the experimental results (27). The requirement for multivalent engagement of TCRs by pMHC ligands in CD4+ T cells has been widely shown (45, 47, 54, 55). The same requirement was shown in CD8+ T cells by Stone and Stern (56).

In this paper, we have made use of the concept of *mean time to signal initiation* (MTSI or stochastic criterion) as a method to quantify the effect of TCR clustering on the timescales of T cell responses and, thus, to compare the behavior of naive and memory T cells. This criterion has also allowed us to compare the results in Section 3 for dimeric binding with those of Ref. (26) for monomeric binding. The introduction of the cross-linked ligand–receptor complex as the minimum *signaling* unit gives the response greater sensitivity to small differences in ligand affinity.

A recent and novel feature of TCR immunology is the existence of TCR nano-clusters that are pre-formed, independently of ligand (3). This suggests that a simple stoichiometric clustering model (oligomerization of free TCRs diffusing on the T cell membrane) is enough to account for the distribution of TCR nano-clusters. In the case of naive T cells, Model 1 predicts an effective nondimensional parameter, *b* = *q*+/*q*−, that allows us to explain the experimentally observed TCR cluster distributions. The presence of larger lipid rafts on the membrane of memory T cells might provide support for the different values of *b* for naive and memory cells, *b*naive and *b*memory, respectively. It has recently been shown that receptor diffusion within the raft is inhibited due to protein anchorage (52). This reduction in the TCR diffusion coefficient

**FIGURE 5 | (A)** Comparison of the numerical solution of Model 2 (Gillespie algorithm with the parameters summarized in**Table 2**) and the approximate solution [equation (3)] for ligand 4P and the same parameters as in **Figure 4**. **(B)** Comparison of the mean MTSI for naive (red dashed line) and memory (solid black line) T cells.

would increase the time required for the receptor to escape from the raft [in a similar fashion as other escape problems (57)]. This escape time is inversely proportional to the diffusion coefficient itself. A smaller TCR diffusivity, as would be the case for memory T cells, will imply a larger residence time in the raft, which in turn will increase the probability of receptor aggregation in a given TCR cluster. A more detailed model of TCR diffusion and aggregation on the T cell membrane will be the subject of future work.

Equation (3) shows the explicit dependence of the MTSI, *T*(*N*, τ ), on the parameter *NB*, for given values of *N* and τ . *NB* is the average number of dimeric receptor clusters per T cell, so that *NB* = *b* × *NR*, with *NR*, the average number of TCRs per T cell (see **Table 1**). For large ligand concentration, the predicted T cell response time for memory and naive T cells is the same, and is equal to τ . In the case of intermediate concentrations, the MTSI is proportional to √ 1 *b* . Finally, for low ligand concentration the MTSI is proportional to 1 *b* . This implies that, at low ligand concentration, TCR pre-clustering alone can only account for at most 72% of the reduction in the response time between memory and naive T cells. This behavior is illustrated in **Figure 5A**. This difference is not so large as to be able to account for the observed higher responsiveness of memory T cells. Our results, thus, point to the need for additional mechanisms beyond TCR pre-clustering.

A potential candidate to explain the large differences between memory and naive T cell responses is the conformational change of the CD3 complex (58). This conformational change is essential to enable ITAM phosphorylation and, thus, the transfer of the TCR signal from the ecto-domain to the cytoplasmic tail of the TCR (58). Conformational changes in the CD3 complex occur as a result of the αβ hetero-dimer binding to pMHC. These conformational changes allow the subunits of the CD3 complex (the γ and δ hetero-dimers and the ζ ζ homo-dimer) to become accessible to Lck, which can then phosphorylate their cytoplasmic domains at the ITAMs, leading to T cell signaling (59). In this way, the ligand-induced conformational change of the receptors can be propagated to all the receptors in the same cluster, so that larger clusters would benefit from this conformational change as a *cascade* [see, for example, Ref. (60) and references therein]. Thus, differences in the distribution of cluster sizes could, indeed, explain the immunological differences between memory and naive cells.

Other membrane receptors also exhibit pre-clustering and ligand-induced receptor cross-linking. For instance, in the case of the vascular endothelial growth factor receptor (VEGFR), it has been shown (61) that there are two distinct pathways to receptor dimerization: (i) dynamic pre-dimerization (as the one described in Model 1), and (ii) ligand-induced receptor dimerization. The main conclusion in Ref. (61) is that both mechanisms are almost indistinguishable at low ligand concentration. However, the first mechanism is more sensitive to changes in the binding affinity at large ligand concentration. Although the biological system studied in Ref. (61) is different from the T cell receptor considered here, their conclusions might be generalized as both receptors are tyrosine kinases.

Bachmann et al. (62, 63) considered a model of diffusion and ligand-induced TCR clustering. Their model suggests that the existence of large enough clusters greatly inhibits subsequent multimer diffusion, thus, reducing the relevance that this mechanism might have. This inhibition might be experimentally tested by exploiting the differences between naive (small and few clusters) and memory (large and many clusters). It will be interesting to make use of the models introduced in this paper to investigate the different roles of ligand binding and cellular activation (62), and TCR turnover (64).

Finally, the existence of TCR pre-clusters [and the knowledge of their membrane distribution given by π*n*, equation (2)] can be considered in the kinetic-segregation model (65). In this model, diffusion out of close-contact zones would be inhibited by the existence of nano-clusters, thus, enhancing the number of triggered receptors. In a similar way, consecutive receptor phosphorylation events (66) in TCR nano-clusters would also amplify receptor signaling.

#### **5. MATERIALS AND METHODS**

#### **5.1. EXPERIMENTS**

Naive CD8+ OT-1 T cells, which recognize an ovalbumin-derived peptide presented by the MHC class I molecule H-2Kb , were isolated from superficial and mesenteric lymph nodes of OT-1 TCR transgenic mice (67), via depletion of CD19+ B cells, CD4+ helper T cells and CD11b+ macrophages, using antibodies and Dynal magnetic beads (Invitrogen). Memory OT-1 T cells were generated by adoptively transferring 106 naive OT-1 T cells into congenic C57BL/6 Ly5.1 Pep3b mice, which were simultaneously immunized with 107 PFU MVA-OVA (68). After 6 months, resting memory OT-1 T cells were isolated from the spleen and lymph nodes of these mice by antibody-mediated depletion of macrophages, B cells, and CD4+ T cells, followed by separation of the OT-1 memory T cells from host-derived Ly5.1+ CD8+ T cells via fluorescence-activated cell sorting, using a Ly5.1-specific antibody. Labeling of cells with the CD3-specific antibody 2C11 and 10 nm gold-conjugated protein-A, replica generation and analysis were performed as previously described (17).

### **5.2. MODELS OF SIGNALING WITH DIMERIC AND TRIMERIC RECEPTOR CLUSTERS**

In Section 2, we introduced a model in which ligands are bivalent and receptor clusters are dimeric (that is, composed of two monomeric TCRs). This is, of course, a first approximation that neglects the distribution of cluster sizes experimentally observed. Yet, the results of our mathematical study only change in a quantitative way, but not qualitatively, when we include TCR clusters of larger sizes. In this Section, we illustrate this by considering a system in which clusters of size 1, 2, and 3 coexist and the ligands are bivalent. **Table 3** provides the notation introduced to describe the molecular species considered in the model, as well as a graphical representation.

At large initial ligand concentration, under the stochastic criterion, the MTSI tends to τ . On the other hand, at low initial ligand concentration, the number of receptors, compared to the number of ligands, is so large that we can neglect molecular species *x*4, *y*4, *y*6, and *y*7, which involve more than one bivalent ligand. This has also been confirmed experimentally. Given our stochastic T cell response criterion, in this case, the signaling units correspond to molecular species *x*5, *y*5, and *y*7. Molecular species *z*1 and *z*3 do

All the variables correspond to the total number of molecular species (not concentrations). Hence, all the rates in the mathematical model have units of s−1 .

not contribute to the T cell response and will be neglected in what follows. Thus, we only need to consider the dynamics of dimeric and trimeric T cell receptor clusters.

We introduce the total number of *signaling units*, *S*5(*t*) ≡ *x*5(*t*) + *y*5(*t*). The set of ordinary differential equations for the model is given by:

$$\begin{aligned} \dot{\mathbf{x}}\_1 &= -4k\_+\mathbf{x}\_1\mathbf{x}\_2 + k\_{\text{off}}\mathbf{x}\_3, \\ \dot{\mathbf{x}}\_2 &= -4k\_+\mathbf{x}\_1\mathbf{x}\_2 + k\_{\text{off}}\mathbf{x}\_3, \\ \dot{\mathbf{x}}\_3 &= 4k\_+\mathbf{x}\_1\mathbf{x}\_2 - k\_{\text{off}}\mathbf{x}\_3 - k\_2\mathbf{x}\_3 + 2k\_{-2}\mathbf{x}\_5, \\ \dot{\mathbf{x}}\_5 &= k\_2\mathbf{x}\_3 - 2k\_{-2}\mathbf{x}\_5, \\ \dot{\mathbf{y}}\_1 &= -6k\_+\mathbf{y}\_1\mathbf{y}\_2 + k\_{\text{off}}\mathbf{y}\_3, \\ \dot{\mathbf{y}}\_2 &= -6k\_+\mathbf{y}\_1\mathbf{y}\_2 + k\_{\text{off}}\mathbf{y}\_3, \\ \dot{\mathbf{y}}\_3 &= 6k\_+\mathbf{y}\_1\mathbf{y}\_2 - k\_{\text{off}}\mathbf{y}\_3 - 2k\_2\mathbf{y}\_3 + 2k\_{-2}\mathbf{y}\_5, \\ \dot{\mathbf{y}}\_5 &= 2k\_2\mathbf{y}\_3 - 2k\_{-2}\mathbf{y}\_5, \end{aligned}$$

where *k*+ = *k*on/(*VNA*), *V* is the volume of the experiment and *NA* is Avogadro's number.

Given the symmetry of the problem, and in the limit of low initial ligand concentration, we will assume that the ratio of *x*3 to *y*3 is that of the initial ratio of free TCR dimers to free TCR trimers, namely,

$$\frac{\wp\_3}{\varkappa\_3} \simeq \frac{\pi\_3}{\pi\_2} = b \quad \Rightarrow \quad \wp\_3 \simeq b \rtimes\_3,\tag{4}$$

where we have made use of equation (2) to conclude π3 π2 = *b*. Thus, the total number of signaling units, *S*5(*t*), obeys the following differential equation

$$
\dot{S}\_5 = \dot{x}\_5 + \dot{y}\_5 = k\_2 \left( 1 + 2b \right) x\_3 - 2 \begin{array}{c} k\_{-2} \ S\_5 \end{array} . \tag{5}
$$

Finally, in the low ligand concentration limit as above, let us introduce *S*3 ≡ *x*3 + *y*3. It is easy to show that equation (5) reduces to

$$
\dot{S}\_5 = k\_2 \frac{1+2b}{1+b} \text{ S}\_3 - 2 \text{ k}\_{-2} \text{ S}\_5 \text{ }, \tag{6}
$$

which is identical to the differential equation for *x*5 above, but with *S*5,3 replaced by *x*5,3, respectively. This means that, except for a prefactor 1+2*b* 1+*b* [which, for *b* ∈ (0, 1), is between 1 and 3/2], the study of dimeric and trimeric clusters is reduced to the dimeric case.

#### **5.3. A SIMPLE FORMULA FOR THE MTSI**

The basic idea behind the stochastic criterion is to count the cumulative number of events that may contribute to signaling (26). Here, we calculate the mean number of cross-linking events up to time *t*, *C*(*t*), as the integral,

$$C(t) = k\_2 \int\_0^t \varkappa\_3(s) \, ds\,. \tag{7}$$

It is possible to obtain an expression for *x*3(*t*) with the approximation that the product *x*1(*t*)*x*2(*t*) is constant, so that the pair of equations for *x*3(*t*) and *x*5(*t*) can be solved exactly. This yields (69):

$$C(t) = k\_2 \left[ \frac{c\_1}{\lambda\_1} (\lambda\_1 + 2k\_{-2}) (e^{\lambda\_1 t} - 1) \right.$$

$$+ \frac{c\_2}{\lambda\_2} (\lambda\_2 + 2k\_{-2}) (e^{\lambda\_2 t} - 1) + a\_1 t \right],\tag{8}$$

where

$$\begin{split} c\_{1} &= \frac{-4\lambda\_{2}k\_{\mathrm{on}}\rho N\_{B}}{(\lambda\_{2}-\lambda\_{1})(4k\_{\mathrm{on}}\rho k\_{2}+2k\_{\mathrm{off}}k\_{-2}+8k\_{\mathrm{on}}\rho k\_{-2})}, \\ c\_{2} &= \frac{4\lambda\_{1}k\_{\mathrm{on}}\rho N\_{B}}{(\lambda\_{2}-\lambda\_{1})(4k\_{\mathrm{on}}\rho k\_{2}+2k\_{\mathrm{off}}k\_{-2}+8k\_{\mathrm{on}}\rho k\_{-2})}, \\ \lambda\_{1,2} &= \frac{1}{2}\left(-4k\_{\mathrm{on}}\rho - k\_{\mathrm{off}} - k\_{2} - 2k\_{-2}\right.\end{split}$$

$$\begin{split} \pm \left[(4k\_{\mathrm{on}}\rho + k\_{\mathrm{off}} - k\_{2} - 2k\_{-2})^{2} + 4k\_{\mathrm{off}}k\_{2}\right]^{1/2}, \\ a\_{1} &= \frac{8k\_{-2}k\_{\mathrm{on}}\rho N\_{B}}{4k\_{\mathrm{on}}\rho k\_{2} + 2k\_{\mathrm{off}}k\_{-2} + 8k\_{\mathrm{on}}\rho k\_{-2}}, \\ a\_{2} &= \frac{4k\_{2}k\_{\mathrm{on}}\rho N\_{B}}{4k\_{\mathrm{on}}\rho k\_{2} + 2k\_{\mathrm{off}}k\_{-2} + 8k\_{\mathrm{on}}\rho k\_{-2}}, \end{split}$$

and *NB* is the number of dimeric receptors. The MTSI is then given by the solution of the equation *C*(*T*(*N*, τ )–τ ) = *N* exp(2*k*−2τ ).

The expressions in equation (3) are obtained from equation (8) in the appropriate regimes. At low ligand concentration, *C*(*t*) is simply proportional to time: *C*(*t*) ' *k*2*a*1*t*, so that *C*(*T* − τ ) = *k*2*a*1(*T* − τ ) = *N* exp(2*k*−2τ ). When λ1,2τ 1, on the other hand, the first non-zero term in a Taylor expansion of *C*(*t*) in time is quadratic: *C*(*t*) ∝ *t* 2 . This provides the exponent 1/2 in the second line of equation (3).

#### **ACKNOWLEDGMENTS**

We thank Ed Palmer, Wolfgang Schamel, and Thomas Höfer for helpful discussions. We also thank the Max Planck Institute for the Physics of Complex Systems (Dresden) and the International Centre for Mathematical Sciences (Edinburgh),where part of this work was discussed and presented, for their hospitality. This work has been partially supported through Grants No. FIS2009-12964-C05- 03 (Mario Castro, Grant Lythe, Carmen Molina-París), BFU2009- 08009 from the Ministerio de Ciencia e Innovación (Hisse M. van Santen), FP7 PIRSES-GA-2008-230665 and PIRSES-GA-2012- 317893 (Mario Castro, Grant Lythe, and Carmen Molina-París), BBSRC BB/F003811/1 (Grant Lythe and Carmen Molina-París), and BBSRC BB/G023395/1 (Carmen Molina-París).

#### **REFERENCES**

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 30 September 2013; accepted: 15 March 2014; published online: 30 April 2014.*

*Citation: Castro M, van Santen HM, Férez M, Alarcón B, Lythe G and Molina-París C (2014) Receptor pre-clustering and T cell responses: insights into molecular mechanisms. Front. Immunol. 5:132. doi: 10.3389/fimmu.2014.00132*

*This article was submitted to T Cell Biology, a section of the journal Frontiers in Immunology.*

*Copyright © 2014 Castro, van Santen, Férez, Alarcón, Lythe and Molina-París. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Theories and quantification of thymic selection

# **Andrew J.Yates\***

Departments of Systems and Computational Biology, Microbiology and Immunology, Albert Einstein College of Medicine, New York, NY, USA

#### **Edited by:**

Rob J. De Boer, Utrecht University, Netherlands

#### **Reviewed by:**

Anton van der Merwe, University of Oxford, UK Viktor Müller, Hungarian Academy of Sciences and Eötvös Loránd University, Hungary Johannes Textor, Utrecht University, Netherlands

#### **\*Correspondence:**

Andrew J. Yates, Departments of Systems and Computational Biology, Microbiology and Immunology, Albert Einstein College of Medicine, 1300 Morris Park Avenue, New York, NY 10461, USA e-mail<EMAIL_ADDRESS>
# **INTRODUCTION**

Conventional (CD4+ and CD8+) T cells are an integral part of adaptive immune systems in vertebrates. A key stage in their development is the creation of the T cell receptor (TCR) through a stochastic process of gene rearrangement. The resulting preselection TCR repertoire has the potential to recognize a very large array of peptides derived both from self and from foreign organisms, presented on Major Histocompatibility Complex (MHC) molecules on the surfaces of other cells. Much of T cell development occurs in a specialized organ in the chest called the thymus, within which this diverse potential repertoire of TCR is vetted. A process referred to as positive selection removes cells with TCR conformations that are generally non-responsive to self-peptide– MHC ligands (self-pMHC), and negative selection removes cells that are overly reactive to self-pMHC and pose a threat of autoimmune responses. The post-selection repertoire exported from the thymus comprises T cells that are largely non-responsive to self, yet capable of responding with remarkable specificity to foreign peptides.

There is a very extensive literature relating to thymic development and selection [for reviews, see for example Ref. (1–3)], but here we summarize the key ideas briefly (**Figure 1**). Conventional T cells begin life as lymphoid progenitors, which migrate from the bone marrow to the inner, cortical region of the thymus and begin a process of proliferation and maturation. Early in development in the cortex thymocytes are referred to as double negative (DN), lacking expression of the CD4 and CD8 co-receptors that are involved in TCR signaling. The TCR comprises two chains and is formed by a multi-step gene rearrangement process that first generates the TCRβ, γ , and δ chains (a small proportion of cells diverge at this stage to seed the γ δ T cell lineage) and then the TCRα chain at around the transition from the DN to CD4+CD8+

The peripheral T cell repertoire is sculpted from prototypic T cells in the thymus bearing randomly generated T cell receptors (TCR) and by a series of developmental and selection steps that remove cells that are unresponsive or overly reactive to self-peptide–MHC complexes.The challenge of understanding how the kinetics ofT cell development and the statistics of the selection processes combine to provide a diverse but self-tolerant T cell repertoire has invited quantitative modeling approaches, which are reviewed here.

**Keywords: thymic selection, T cells, mathematical modeling, repertoire selection, theoretical biology**

(double positive, DP) stage. TCRαβ cells then migrate among cortical thymic epithelial cells and dendritic cells, auditioning for the ability to recognize self-pMHC. There is evidence that DP cells with non-functional TCR can undergo repeated TCRα rearrangements (4) to re-audition. Positively-selected cortical thymocytes begin negative selection and eventually move to the outer capsule of the thymus, the medulla. There they complete negative selection through interactions with medullary thymic epithelial cells and dendritic cells. TCRαβ thymocytes, which recognize selfpeptides presented on MHC class I or class II below an acceptable threshold of reactivity develop into the CD8 SP (single-positive, CD4−CD8+) or CD4 SP (CD4+CD8−) lineages, respectively, and are eventually exported into the peripheral circulation as naive T cells.

The topic of thymic selection has received substantial attention from the immunological modeling community, perhaps for two main reasons. First, selection has widely been viewed as a well-delineated optimization problem – how to craft a TCR repertoire that covers the space of possible pMHC ligands as widely as possible, while preserving sufficient specificity to discriminate between self and foreign (and between different foreign) peptides? This question naturally invites quantitative arguments. Second, the biology is well-characterized – a relatively small number of cell types and modes of interaction appear to be involved, and large amounts of experimental data are available. These simplify and constrain the construction of models.

Modeling studies have focused on many aspects of thymic selection but many questions and uncertainties remain. What are the rates and efficiencies of passage through the different phases of development and selection, and in what thymic microenvironments do each take place? How do thymocytes integrate signals received from interactions with pMHC to make fate decisions?

What are the relative contributions of the MHC itself and its associated peptide to TCR signaling and fate determination? What influence do each of these have on the post-selection repertoire's diversity and coverage of the pMHC universe, and its ability to discriminate between self and foreign?  How many TCR interactions contribute to a thymocyte's fate decisions? What evolutionary pressures have determined the typical number of MHC alleles we possess? There have been many different theoretical approaches to these questions – from mean-field population dynamic models of progression through developmental stages, to probabilistic models of selection, to explicitly spatial models of migration within the thymus.

This review groups studies of these topics into broadly labeled categories, but in some cases the grouping is arbitrary – many of these questions are related and have been addressed either alone or in combination. The review has a bottom-up structure, beginning with an overview of experimental quantification of selection and modeling of thymocyte population dynamics. It then moves to studies of higher-level properties of the T cell repertoire, such as TCR cross-reactivity, and concludes with the problem of optimal within-individual MHC diversity.

### **THE POPULATION DYNAMICS OF THYMOCYTES**

Basic elements of a quantitative understanding of thymic development are the steady-state population sizes of different developmental stages, the mean times to transit between them and the proportion surviving at each stage, which we refer to as the efficiencies of selection. While some quantities can be experimentally determined, mathematical models have helped us develop a more complete description of the kinetics of selection, both for the thymocyte population as a whole and for the CD4 and CD8 lineages in isolation.

To estimate the parameters of a dynamical system usually involves observing its response to perturbations. One method is to follow cohorts of cells as they progress through development using intra-thymic injection of a dye or radioisotope label (5–8). Arguably this method is less disruptive than cell transfers, but the uptake of marker can be heterogeneous (5,7) and measurements of death rates using injected dyes rather than congenic markers may be confounded by loss of label (9). More recently, methods have included using GFP (green fluorescent protein) expressed during TCR rearrangement, its decaying intensity then a marker for time spent in development (10); inducible TCR signaling can be used to arrest, release, and follow cohorts of cells from the early DP stage (11); and small numbers of labeled thymocytes isolated at different developmental states can be followed after intra-thymic injection (11, 12). The population dynamics have also been exposed by transiently depleting thymocytes and observing the system's return to equilibrium (13). Various experimental systems, with or without associated dynamical models, are in general agreement over several quantitative aspects of thymic development but inconsistencies and uncertainties remain.

#### **SELECTION EFFICIENCIES AND CELL FLUXES**

Thymocytes begin to select against self-pMHC ligands at the DP stage following TCR rearrangement and so we focus on survival, proliferation, and differentiation from this stage onward. The proportion of DP cells that reach maturity (that is, survive both positive and negative selection) is widely agreed to be 5% or less (6, 11, 13–16). Within this pruning process, the general view is that positive selection is the most stringent, with 75–80% of cells failing to progress from the earliest DP stage, suggesting the majority of TCR generated are unable to recognize peptides in conjunction with MHC class I or II to any useful degree (11, 13, 15, 17, 18). Many studies have estimated that between 20 and 50% of positively-selected thymocytes then survive negative selection (11, 17, 19–22), although Itano and Robey (8) estimated a selection efficiency as high as 90% for DP cells into the CD4 SP lineage.

The rate of production of mature CD4 and CD8 cells in the thymi of young adult mice is roughly 1% of total thymocytes or 1 − 3 × 106 cells/day, a figure arrived at by a variety of labeling methods (5–7). Egerton et al. (6) estimated this to be just over 3% of the rate of entry into the DP population, meaning that fueling this trickle of output requires that roughly 30% of all thymocytes enter the DP stage each day. This again illustrates the extent of the filtering of the pre-selection repertoire that appears to be required to produce a functional and self-tolerant population of naive T cells. The thymus gradually involutes and its rate of output declines with age in both mice (23) and in humans (24), indicating that the bulk of the peripheral T cell repertoire is probably generated early in life.

### **THE MAJORITY OF THYMOCYTE DIVISION LIKELY OCCURS PRE-SELECTION**

Labeled nucleotide uptake assays have revealed that substantial proliferation of thymocytes occurs before selection on self-pMHC ligands begins, stopping at or around the time of TCR rearrangement at the late DN/early DP stage (6, 25, 26). However, it is proliferation following TCR rearrangement that is most relevant for understanding how repertoire diversity is generated. Division during selection means a smaller proportion of TCR clonotypes may pass selection than measures of percentage survival suggest (27). The extent of division early in selection is unclear – estimates of the proportion of newly generated DP cells that are dividing have ranged from 11 to 68% (6, 25, 28), and CFSE labeling in *in vitro* thymic organ cultures showed up to 5 divisions from DP onward (29). However, the DP population comprises cells preand post-TCR rearrangement, and there appears to be very little proliferation within the more mature DP population (6, 11, 15). There is a low level of proliferation during or just before the SP stage (11, 13, 25, 30), with CD8 SP more prone to division than CD4 SP (10).

Perhaps the most reliable experimental measure of average levels of proliferation during selection uses T cell receptor excision circles (TRECs). TRECs are circular DNA fragments that are stable remnants of the recombination events that generate the TCR and are shared randomly between daughter cells on division. The mean TREC content per cell is a rough measure of the mean number of divisions that have taken place since the TCR was generated. One caveat is that TREC studies are used most commonly in humans and much of what we discuss here derive from studies in mice. Another is that standard TREC measurements contain no information about the variance of the division number, and may gloss over even quite extreme heterogeneity in division patterns. Nevertheless a study of human infants observed 1–2 divisions on average between TCR rearrangement at the CD3low CD4+CD8+ stage and mature CD4 or CD8 SP; once shortly after TCR rearrangement, and another at the CD8 (but not CD4) SP stage (31). The high TREC content they observed at the early DP stage may reflect multiple rearrangements taking place in order to generate a functional TCRα-chain. In line with these results, the TREC content of naive CD31+CD4+ recent thymic emigrants in human infants is ~0.1–0.9/cell (32), suggesting that up to three divisions take place on average between TCR rearrangement and export to the periphery, although this may include some post-thymic proliferation and so is an upper limit on the extent of intra-thymic division.

#### **TURNOVER RATES AND TRANSIT TIMES**

Experimental estimates of the times taken to transit different developmental stages (immature DP→mature DP→SP→ Export) are variable, particularly within the SP population (6, 10, 12, 25). Possible reasons for these discrepancies include different labeling protocols, different gating strategies defining thymic subpopulations, heterogeneity of cell populations, and differences in the kinetics of MHC class I-restricted and class II-restricted lineages. It has also been unclear whether selection is a "conveyor belt," first-in first-out, or has a more stochastic "lucky dip" nature (25). From a modeling perspective these are two points on a continuum. If an experimentally identifiable developmental stage comprises several shorter, sequential differentiation steps, the variance in the transit time through that stage is lowered with respect to a single-step model of transit. The more obligate steps, the more conveyor-belt-like the system appears.

There is general agreement that the transition from nondividing mature DP to SP takes on average 3–4 days (6, 12, 15, 28), although it has been argued that it takes significantly longer to reach CD8 SP than CD4 SP (33). This transition is dependent on TCR signaling (15, 34). Observing a well-defined delay in the appearance of labeled SP cells, Egerton et al. (6) argued for a first-in-first-out kinetic in the DP population. This suggests DP cells must transit through a number of obligate steps. Subsequent experimental and modeling studies have addressed this, and are

discussed below. The same study estimated a mean SP residence time of ~12 days, comparable to other estimates of the medullary residence time (6, 28). McCaughtry et al. (10) argued that this is an overestimate of the time mature conventional SP T cells take to develop, because the SP population is heterogeneous, also containing Treg, NKT, and γ δ T cells, which turn over more slowly. They estimated SP CD4/CD8 residence times to be 4.4/4.6 days. Saini et al. (33) arrived at similar estimates. As for DP cells, there are may be several developmental stages within the SP population and so it seems unlikely that SP residence times are exponentially distributed.

Stritesky et al. (12) estimated the total rate (cells per unit time) at which cells are negatively selected to be almost six times greater than the rate of positive selection, and found that both processes occur predominantly at the DP stage. Converting these figures into the relative efficiencies of positive and negative selection requires knowledge of how long cells spend in each selecting phase. If indeed positive selection is the more stringent, their result indicates that negative selection must take place over a relatively short timescale within the DP compartment. This is supported by a recent study observing negative selection of DP thymocytes taking place over ~12 h (35).

Interpreting data on transit or residence times can be problematic when both death and differentiation are taking place, as they clearly are at the DP stage(s) of development. If death and differentiation are modeled as independent processes, then at equilibrium transit rates through a compartment are not necessarily the same as turnover rates. If cells are maturing at rate µ and dying at rate δ, the population turns over at rate µ + δ and the expected time a cell spends in that compartment is 1/(µ + δ). However, the mean time that successfully differentiating cells spend in each compartment is shorter because it is conditioned on survival, and is µ/(µ + δ) 2 (if cells are capable of maturing but are simultaneously at risk of dying, those that successfully mature tend to do so early). This difference can be quite substantial, as we see below.

#### **KINETIC MODELS OF THYMIC DEVELOPMENT**

Data from these experimental studies and others have invited the use of population dynamic models to infer the kinetics of development. In the first studies to model thymic development, Mehr and collaborators utilized ordinary differential equation (ODE) models of the flow from DN→early DP→late DP→CD4/CD8 SP (36, 37). They utilized measures of steady-state population sizes and parameters either inferred from data or explored systematically to ask questions about the underlying dynamics. Mehr et al. (36) argued that positive selection likely involves triggering of proliferation as well as rescue from death, and while they were unable to use the steady-state data to make strong statements about the timing of positive versus negative selection, they inferred that most death at the DP stage is due to failure to positively select, consistent with many experimental and subsequent modeling studies.

There is evidence from fetal thymic organ cultures that populations of mature CD4+ T cells resident in the thymus may enrich for the CD4 lineage while reducing thymic output. Mehr et al. (37) used a similar model with these data to propose that the mature resident cells increase survival of developing single-positive CD4 T cells while reducing proliferation or increasing the rate of differentiation of DP cells. They suggest that mature CD4 T cells exert their influence by restricting the number of available pMHC ligands in the thymus, which could simultaneously reduce proliferation of DP cells (lowering thymic output) and decrease the stringency of negative selection (increasing the efficiency of maturation into the mature SP state). Again, these conclusions were reached using data from the thymus at steady-state.

Mehr and collaborators also studied the seeding of the cortical stroma with bone marrow-derived progenitor cells using a combination of modeling and experiment. They showed how migration between niche sites explained the competitive advantage of younger progenitors over older (38, 39), and that reconstitution of the progenitor population following irradiation is limited by damage to stromal niches and incumbent, surviving cells (40).

Thomas-Vaslin et al. (13) studied naive T cell homeostasis from the thymus through to the periphery. They induced systemic depletion of T cells for 7 days through expression of a suicide gene in dividing cells, and followed the kinetics of reconstitution. To interpret these data they developed a multi-compartment ODE model of T cell development, with a finer-grained treatment of transit through the DN, DP, and SP stages. In their model extensive proliferation occurs through the DN to early DP, with the latter population dividing 5 times. Their best-fitting model assumes all cell death (positive and negative selection) takes place at the late DP stage. They estimated 5% of total thymocytes (DN, DP, and SP) or ~3 × 106 are exported as naive SP cells per day, and that 93% of DP thymocytes are lost, in line with existing estimates, and again suggesting that the bulk of negative selection occurs at DP. The mean times spent overall in the early DP (dividing), late DP (selecting), and SP compartments were estimated to be 1.2, 2.7, and 5.8 days respectively.

Sinclair et al. (11) used a different experimental system, with controllable TCR signaling that allowed arrest and release of cells at the early DP stage, and used a multi-compartment ODE model to quantify transit dynamics and selection efficiencies. Rather than simply early or late, they broke the DP stage into a branched developmental progression defined by the expression levels of CD5 and the TCR (33). In their schema, DP1 thymocytes are pre-selection; progression to DP2 requires a positively-selecting TCR signal; DP2 thymocytes consist of class I- and class II-restricted thymocytes in the first 12–48 h of development; and DP3 thymocytes are predominantly MHC class I-restricted cells that can select into CD8SP only. Thus cells destined for CD4SP transit DP1-DP2 only, and CD8SP transit through DP1, DP2, and DP3.

Sinclair et al. (11) estimated that ~75% in DP1 fail to progress to DP2, reflecting failure to positively select and dying of neglect. Overall, 5% of DP cells become CD4SP and ~2% become CD8SP, and so ~94% of DP cells are lost. They also saw relatively low levels of cell death in the SP compartment. These results suggest again that the bulk of negative selection occurs before cells transition to SP. They saw very little proliferation in their system, using a variety of methods, and so did not model cell division. Mean residence times in DP1 and DP2 were 3.5 and 1.4 days, respectively, with the smaller CD8 lineage spending an additional 7 days in DP3. They estimated 23% of all thymocytes at DP and SP enter the DP compartment per day. These selection efficiencies and the net flux agree with other estimates. Accounting for the selection bias

on maturing cells, the model predicts that successful thymocytes spend on average 1.3 days in DP1 + DP2, 4.5 days in DP3. SP4 and SP8 residence times were 5 and 3.7 days, respectively, with very little cell death occurring. Their analysis therefore suggests that CD4SP/CD8SP cells take ~6.3/9.5 days from entry into DP1 to export.

# **MIGRATION WITHIN THE THYMUS AND THE TIMING OF POSITIVE AND NEGATIVE SELECTION**

From the perspective of modelers attempting to connect models of thymocyte dynamics to data, it is important to understand when and where the different phases of development and selection occur. Selection begins in the thymic cortex, where the majority of thymocytes perform undirected random walks (41) encountering pMHC on cortical thymic epithelial cells. Sensitivity to medullary chemokine receptor signals begins to increase immediately following receipt of a positive selection signal and positively-selected cortical thymocytes eventually display rapid, directed motion toward the medulla (41), where they encounter pMHC on medullary thymic epithelial cells and dendritic cells. Negative selection takes place in the medulla (35, 42–44) but also late in migration through the cortex (45) and possibly even throughout development (46). The mapping between these migratory and selecting processes to developmental stages is not clearly defined. Cells undergoing negative selection in the medulla include DP populations (35), indicating that maturation from DP to SP does not coincide precisely with the cortical–medullary transition but further supporting the conclusion that the extensive cell loss at the DP stage comes from failure of both positive and negative selection. Further, antigen-presenting cells in the cortex and medulla appear to differ in their ability to provide positive or negative selection signals, either through differences in pMHC expression or diversity, or levels of co-stimulation (47–51). It seems therefore that negative selection at the DP stage takes place in at least two distinct spatial and TCR-stimulatory environments.

#### **MODELS OF SELECTION WITHIN THE CORTEX AND MEDULLA**

Motivated by this, Faro et al. (52) took a different perspective; rather than partitioning selecting thymocytes into developmental stages, they used a probabilistic model to describe selection within the cortex and the medulla. They aimed to quantify the number of selecting events, the number of selecting APC encounters and pMHC engagements, and the efficiencies of positive and negative selection in each region. Using the experimental estimates of overall selection efficiencies, and one experimental estimate of the efficiency of negative selection in the medulla, they inferred that most thymocyte death occurs by failure to positive select in the cortex, and cells are ~10 times more likely to be deleted (negatively selected) in the medulla than in the cortex. With these efficiencies, through a parameter search, they were able to infer the number of ligands each thymocyte selects on in each spatial compartment. They came to the striking conclusion that for each cortical thymocyte selection takes places on <60 pMHC ligand interactions, likely in order to achieve in their model the required high level of failure to positively select. However, this needs to be reconciled with the ~3-day mean lifetime of cells at DP1, which suggests cells have far more opportunities to positively

select, either through repeated encounters with APC or through repeated rearrangements of the TCRα chain [see Ref. (53) and refs therein], before dying of neglect.

#### **IDENTIFYING THE SOURCE OF THE CD4:CD8 LINEAGE BIAS IN THYMUS**

CD4 SP outnumber CD8 SP by roughly 4:1 in the thymi of many species. Using time courses of development in control mice and those lacking MHC class I or class II, Sinclair et al. (11) estimated the CD4 and CD8 lineage-specific selection efficiencies. In control animals, the highest death rate was at the positively-selected DP2 stage, and was substantially greater for MHC class I-restricted cells. MHC class I- and class II-restricted cells are indistinguishable at DP1 and DP2, but they were able to back-calculate the rates of production of precursors of the two lineages after TCR rearrangement, and found they were comparable. This suggests that the CD4:CD8 asymmetry in the thymus derives in large part from more stringent selection acting on MHC class I-restricted cells and not from any significant asymmetry in the predisposition of randomly generated TCR to recognize MHC class I or class II. Theirs is a model of CD4/CD8 lineage commitment in which the ability of a DP thymocyte to recognize MHC class I or class II dictates whether it will progress to the CD8 or CD4 lineages, respectively (8, 54). This is contrast to a less efficient, selective process in which a thymocyte's decision to downregulate either CD4 or CD8 expression is stochastic and decoupled from MHC preference, such that potentially viable TCR may fail positive selection [see, for example Ref. (55, 56); and Ref. (57) for a discussion of a hybrid mechanism]. Mehr et al. (36) proposed a purely instructive model of selection, in which pre-selection thymocytes are in principle able to recognize both MHC class I or II, and concluded that the most likely explanation of the CD4 bias is a difference in the *per capita* rates of maturation from DP into the two lineages, rather than differences in death rates.

The majority of models discussed here assume that thymocytes undergo screening independently. Mehr et al. (36, 37) implicitly allowed for competition with density-dependent proliferation rates at each developmental stage. However, there is some evidence that the probabilities of maturation can be impacted by competition between thymocytes, both globally and in lineage-specific ways. The efficiency of selection of transgenic TCRs varies with their abundance and with the availability of cognate pMHC (15, 58–60), and the selection of polyclonal MHC class I-restricted thymocytes is more efficient in the absence of MHC class II and vice versa (11). These observations suggest that selection efficiencies may be limited by competition both within and between lineages for access to pMHC or other resources needed for selection, and so may impact on the CD4:CD8 ratio emerging from the thymus. Two studies have used explicitly spatial, agent-based models of thymocyte migration and development to investigate this issue. Souza-e Silva et al. (61) modeled the movement of DN, DP, and CD4 SP and CD8 SP populations and their interactions with thymic epithelial cells (TEC) and chemokine gradients, using a 2D model. The structure of the epithelial networks was derived from histological samples from both mice and infant humans. Parameters were chosen to give agreement with published data regarding the repopulation of the thymus after sublethal irradiation, although a sensitivity analysis was not performed. In

their model the CD4:CD8 ratio emerges as a result of competition for access to TEC and stochastic variation in the duration of signaling, which has been associated with CD4/CD8 lineage commitment (62). Their simulations also reproduce an observed variation in the CD4:CD8 ratio as irradiated thymi reconstitute and, in their model, the degree of competition increases. Efroni et al. (63) also took an agent-based approach and concluded that MHC class I and class II ligands on TECs are limiting. If continued access to pMHC stimulation is required for survival, and class I restricted cells stay conjugated to MHC for longer than MHC class II-restricted cells, exclusion of competitors leads to a higher death rate of cells developing into the CD8 lineage and a skewing of the CD4:CD8 ratio. Such a competitive model is an experimentally testable explanation of the differential death rates observed by Sinclair et al. (11).

# **CHARACTERISTICS OF THE TCR REPERTOIRE**

Various summary statistics can be used to describe T cell populations pre- or post-selection. The *diversity* (or the *repertoire*) usually denotes the total number of distinct TCR sequences or clonotypes. The *cross-reactivity* measures a TCR's capacity for discrimination, and is quoted as either the average number or the proportion of different pMHC that one TCR responds to above some defined functional threshold. *Specificity* is inversely related to cross-reactivity. A mirror quantity is the *precursor frequency*, also referred to as the response frequency – the average proportion of all TCR capable of recognizing one pMHC. Further, selection operates in the context of an individual's own MHC alleles. *MHC restriction* measures the degree to which a given TCR is limited to recognizing peptides presented by one or more self-MHC; and *alloreactivity* is the proportion of TCR that respond to a foreign MHC, which is relevant for transplantation of tissues from one individual to another. In the sections that follow we describe how theoretical models have been used to understand how these quantities are linked and constrained by thymic selection.

#### **TCR CROSS-REACTIVITY**

A diverse TCR repertoire seems to be a requirement for coverage of pMHC shape space. However, the number of theoretically possible pMHC complexes appears to be far greater than any individual's capacity for unique TCR clonotypes (64–66); a simple calculation for just one MHC class I variant, assuming it presents 2% of all possible 9-residue peptides, yields 209 × 0.02 ' 1010 possible pMHC, compared with the roughly 5 × 107 naive CD8 T cells in a mouse. To minimize the probability that any given foreign pMHC will escape detection by the immune system, some degree of TCR cross-reactivity therefore seems beneficial. Mason (64) used a variety of methods and data sources to estimate that one MHC class I-restricted T cell responds to between 106 and 107 nonamer peptides, or one in 103 to 104 pMHC using the theoretical estimate of the potential pMHC diversity; and Ishizuka et al. (65) used peptide libraries to estimate more directly that one CD8 T cell clone responds to roughly 1 in 3 × 104 peptide–MHC class I ligands. On the other hand, the average degree of cross-reactivity seems necessarily constrained from above, to avoid excessive deletion of the repertoire and to preserve specificity for self/non-self discrimination. It therefore seems plausible that evolutionary pressures might have optimized this trade-off and determined the degree to which TCR can respond to multiple pMHC.

# **OPTIMAL LEVELS OF TCR CROSS-REACTIVITY – PROBABILISTIC ARGUMENTS**

Several variants of essentially the same argument predict that the diversity of self-peptides involved in selection is the strongest influence on the optimum level of TCR cross-reactivity (64, 67–70). One version of the argument is as follows. The proportion of the positively-selected T cell repertoire *R*0 that avoids deletion, *f*, decreases with both the number of self antigens *Ns* and the crossreactivity *r*, *f* = (1 − *r*) *Ns* which is approximately exp(−*rNs*) for *r* . 1/*Ns* . A pathogen escapes immune recognition if all *f R*0 surviving unique clonotypes fail to recognize (cross-react with) all *x* epitopes it generates, with probability

$$P\_E = (1 - r)^{f R\_0 x} \simeq \exp(-r f R\_0 x) \tag{1}$$

where again the approximation holds if *r* . 1/(*f R*0*x*). This ignores MHC restriction, but including this refinement yields similar conclusions (67). Using the expression for *f*,

$$R\_0 \simeq -\log(P\_E) \frac{\exp(rN\_s)}{r\infty}.\tag{2}$$

This equation connects the repertoire before negative selection *R*0, the probability of immune escape *PE* and the pre-selection cross-reactivity *r*. *R*0 is relatively insensitive to *PE* but very sensitive to the diversity of self, *Ns* . In this model, then, the strongest determinant of the size of the pre-selection repertoire is the diversity of self antigens, *Ns* , and not the requirement for minimizing the probability that a pathogen escapes detection (67).

The three-way relation expressed by equation (2) can then be used to estimate the optimal cross-reactivity under different evolutionary constraints. Suppose the potential repertoire size *R*0 is relatively conserved and evolution has selected for the smallest *PE* by tuning TCR cross-reactivity; in this case, the optimal cross-reactivity is simply the inverse of the number of distinct self-pMHC involved in selection, *r* = 1/*Ns* . The same value of *r* arises if evolution is assumed to minimize the required repertoire size *R*0, whatever the value of *PE* (67). Thus the more diverse the self-peptides involved in thymic selection, the more specific (less cross-reactive) the TCR needs to be. The same result can be derived in a very general way using extreme-value theory (70), requiring only the assumption that the negative selection threshold in the thymus is equal to the activation threshold in the periphery.

The induction of tolerance in the thymus is likely incomplete and there may be mature lymphocytes that are able to recognize self-peptides not involved in thymic selection. Borghans and De Boer (71) argued that to minimize the probability of these cells mounting a cross-reactive autoimmune response to this "ignored self" while responding to a pathogen demands higher levels of specificity than predicted by the simplest models. In this model, optimal cross-reactivity is then modulated by the potential diversity of the repertoire; the greater the number of possible T cell clonotypes, the lower cross-reactivity is required.

Percus et al. (72) took a different approach to studying optimal cross-reactivity, prompted by the observation that the sizes of the binding sites of the TCR and the B cell receptor (antibodies) are similar, at roughly 15 amino acids. They concluded that this size is large enough to provide discriminatory power but small enough that there is sufficient cross-reactivity for coverage of foreign antigen shape space. Interestingly this result does not arise from the demand for self–non-self discrimination, but rather from the constraint of the observation that the B and T cell repertoires comprise ~107 different receptors. However, this diversity itself may be derived from the self-tolerance arguments described above (64, 67–69). It has since been established that substantially fewer peptide residues are involved in TCR recognition. Burroughs et al. (73) analyzed the proteomes of humans and several microorganisms and showed that even the seven exposed (non-anchor) residues of the nine-mer peptides bound to one MHC class I allele may promote self/non-self discrimination, with <0.5% overlap in these sequences between humans and different microorganisms.

#### **CONVERGENT ESTIMATES OF LEVELS OF NEGATIVE SELECTION**

Several of these studies concluded that at the optimal level of cross-reactivity the probability of negative selection is roughly 63%, making various assumptions regarding the magnitude of parameters and maximizing the probability that the post-selection repertoire mounts a response to a foreign pMHC. However, the probability of negative selection can be derived without any assumptions regarding parameter values. From above, the fraction of the positively-selected repertoire with cross-reactivity *r* that survives deletion on *Ns* self-peptides is *f* = (1 − *r*) *Ns* . The probability that the post-selection repertoire *R* = *fR*0 fails to recognize one given foreign pMHC is given by equation (1) with *x* = 1,

$$P\_E = (1 - r)^{f\_{\mathcal{R}\_0}} = (1 - r)^{R\_0(1 - r)^{N\_\mathcal{S}}}.\tag{3}$$

This is minimized with respect to *r* at *r* = 1 − exp(− 1/*Ns*), exactly (the optimal cross-reactivity *r* ' 1/*Ns* then obtains if *Ns* 1). So if evolution acts on cross-reactivity to minimize the probability of foreign pMHC escaping detection, the fraction of the positively-selected repertoire that survives negative selection is then simply *f* = (1 − *r*) *Ns* = exp(−1) ' 0.37, or ' 63% of positively-selected thymocytes are deleted.

Mason (64) arrived at the same result assuming heuristically that the quantity to be maximized is the "reactivity" of the repertoire, proportional to the number of peptides each T cell can recognize multiplied by the proportion surviving negative selection;

Reactivity ∼ Cross-reactivity

× P(survive negative selection) ∼ *r* × (1 − *r*) *Ns* .

Maximizing this reactivity is equivalent to minimizing the probability of escape in equation (3) when *r* is assumed to be small. There, using the Taylor expansion gives *PE* ' 1 − *rR*0(1 − *r*) *Ns* , and so the probability of responding (1 − *PE* ) is ~*rR*0(1 − *r*) *Ns* , or Mason's reactivity. Since *r* is small, the probability of negative selection is (1 − *r*) *Ns* ' exp(−*rNs*) and so the reactivity is proportional to *r* exp(− *rNs*), which is maximal with respect to *r* when argument of the exponential is −1. Thus again *f* ' 0.37 and the optimal cross-reactivity *r* ' 1/*Ns* .

An essentially identical argument applies to negative selection of B cells (67, 69). This estimate of *f* is remarkably consistent with estimates of levels of negative selection in the thymus from several experimental and population dynamic modeling studies (11, 17, 19–22).

#### **ALTERNATIVE TREATMENTS OF CROSS-REACTIVITY**

These models assume a universal cross-reactivity parameter *r*, but T cells may have the capacity to modulate their activation thresholds in response to their signaling environment (74, 75). Motivated by this, Scherer et al. (76) developed a model in which T cells tune their activation thresholds (and thus their cross-reactivity) to the level of their strongest interaction with self-pMHC during selection. If combined with a deletion mechanism that removes cells with activation thresholds so high as to be judged functionally inert, this model appears to be a more efficient mechanism of thymic selection than the standard clonal deletion model. Scherer et al. showed that the tuning model increases the probability of mounting an immune response to a given pathogen epitope, given a pre-selection repertoire size *R*0, and the number of self-pMHC ligands involved in selection, *Ns* . The improvement offered by the tuning model is most striking for small pre-selection repertoires, *R*0 *Ns* , but disappears for *R*0 *Ns* . The latter inequality likely holds for mice and humans; the potential number of unique TCR sequences exceeds the estimated 103–105 self-peptides able to be presented by a given MHC allele (73, 77, 78). Further, equation (1) predicts that at the optimal cross-reactivity *r* = 1/*Ns* , the probability of one epitope (*x* = 1) escaping recognition is *PE* = exp(−*R*0/*eNs*) where *e* is the base of the natural logarithm. For *PE* < 0.05, expected in humans and mice, requires *R*0 & 10*Ns* . Despite this, Scherer et al. (76) argue that the tuning model is a more parsimonious mechanism of self-tolerance in the thymus than the standard model of deletion based on evolutionarily optimized cross-reactivity.

Finally, many of these arguments assumed thymic selection alone optimizes cross-reactivity, but the requirement for memory T cells to discriminate between different pathogens may impose a further constraint of its own (79, 80).

# **EXPLORING CROSS-REACTIVITY WITH SEQUENCE-BASED MODELS OF THYMIC SELECTION**

A series of related papers by Detours, Perelson, and Mehr (27, 81– 83) used a model of TCR–pMHC interactions to understand at a more mechanistic level how cross-reactivity, alloreactivity, and MHC restriction emerge in the post-selection repertoire. Here we focus on their treatment of TCR cross-reactivity, and return to alloreactivity and MHC restriction in the next section. Their starting point was an established model of protein binding (81, 84). They described the interaction between the variable region of the TCR and its pMHC ligand with strings of digits, and binding strengths between each digit pair were determined by the degree of complementarity between their binary representations (81). MHC and peptide contributed additively to the affinity of the interaction, the quantity assumed to drive selection. Given the number of digits ascribed to the polymorphic MHC residues in contact with the TCR, and the number of digits representing the peptide, selection could be performed on a randomly generated TCR repertoire using

randomly generated peptide–MHC complexes. Affinity thresholds were then adjusted to give stringencies of positive and negative selection similar to those observed experimentally.

To circumvent the computational costs of selection using realistic numbers of peptides and unique pre-selection TCRs, they derived expressions for the mean-field predictions of the model for given parameter sets. This has the advantage of yielding population-level statements, which average over all possible TCR, MHC, and peptide sequences.

Detours and Perelson (82) estimated the precursor frequency, the proportion of naive T cells able to respond to a particular foreign pMHC. Experimental estimates of this quantity lie in the range 10−6–10−4 (85–89). They term this the response frequency, *R*, and found it to be strongly and inversely related to the number of selecting self-pMHC ligands. Since precursor frequency is positively correlated with cross-reactivity (64), this result is in keeping with the theoretical studies discussed above (64, 67–70). It is also consistent with observations that repertoires selected on a restricted range of peptides exhibit higher cross-reactivity than normal (90–92). For *R* to lie in the observed range constrains the number of distinct peptides each MHC can present to be of the order 103–105 , in line with estimates for murine MHC class I (77), MHC class II (78), and human MHC class I (73).

To explore the effect of thymic selection on specificity in more detail, Chao et al. (93) revisited the complementary digit-string model. Again peptide and MHC were assumed to contribute additively to an antigenic distance from the TCR, which was inversely related to affinity or the strength of a selecting signal. They confirmed that negative selection reduced the coverage of peptide space, defined as the proportion of peptides that are recognized on the selecting MHC. This was equivalent to a reduction in the cross-reactivity of the repertoire; it reduced the mean antigenic distance to foreign pMHC complexes.

Chao et al. (93) then used the model to address the question of why the number of pMHC that one T cell is able to respond to varies widely across TCR (94). Their simulations suggested that the degree of cross-reactivity to a foreign peptide was inversely related to the peptide's similarity to self, which can be understood with the following argument. In their model, in the pre-selection repertoire a TCR's affinity for the MHC and peptide portions of its ligand are uncorrelated. Selection introduces an inverse correlation between a TCR's affinity for its selecting MHC and its strongest affinity for self-peptide; to select, a TCR's strongest interaction with self must lie between the positive and negative selecting thresholds. (The narrower the range of affinities defining the selecting region, the stronger this correlation will be.) Selected T cells with high affinity for MHC then have a relatively low affinity for the selfpeptide component and require only weak binding to foreign peptide to be activated (activation in their model is defined to be an interaction above the negative selection threshold). These cells are therefore cross-reactive to foreign peptides. Conversely, selected TCR that bind relatively weakly to MHC have higher affinity to self and require strong binding to foreign peptide for activation, and therefore have more specificity for foreign antigen. Thus it emerges from their model that a TCR's specificity to foreign peptide is positively correlated to its affinity for self-peptide; or

equivalently, a TCR's cross-reactivity is positively correlated with its affinity for MHC.

The effect of negative selection on cross-reactivity can be understood with a similar argument. A TCR with high affinity for MHC will survive negative selection only if it has low affinity to all selfpeptides, which is unlikely. Negative selection therefore enriches for cells with lower affinity for MHC, which from the argument above tend to be less cross-reactive. This reduction in coverage means specificity to foreign peptide must be increased.

Kosmrlj et al. (95) used a more physical, mechanistic approach to understanding how negative selection increases specificity, with the aim of characterizing the properties of the amino acid sequences of specific and cross-reactive TCR. Using the Miyazawa– Jernigan matrix (96) to quantify the interaction energies of pairs of amino acids, they extended the digit-string model to calculate the binding affinities between the peptide and the CDR3 region of the TCR, with a constant contribution from the MHC. (The variable peptide element of the pMHC ligand can be assumed to include the polymorphic MHC residues; thus their model may allow for MHC restriction, although this was not discussed.) Košmrlj et al. (97) presents an analytical treatment of the model.

They observed that TCRs selected against multiple peptides on the same MHC had peptide contact residues enriched in weakly interacting amino acids. In their model this arises by a sort of buffering mechanism – such sequences are able to withstand multiple substitutions in the peptide sequence to which they bind most strongly, and so are more resistant to negative selection than those TCR with strongly binding residues. For these TCR to survive selection requires that the invariant MHC contribution to the binding energy is of moderate strength – contributing sufficiently to favor positive selection but well below the negative selection threshold.

Kosmrlj et al. (95) argue that it is this enrichment for weakly binding TCR driven by negative selection that underlies antigen specificity. Antigen recognition is assumed to occur when a TCR signal exceeds the negative selection threshold made up by several interactions. This requires the peptide to contain several amino acids capable of binding the most strongly to the generally weakly binding TCR contact residues. Each contributes significantly to the total binding energy, and so any mutation to the peptide sequence has a high probability of abrogating recognition. Thus there is a restricted peptide signature or "barcode" required to trigger the TCR. In their model, TCR selected against a single pMHC were enriched slightly for strongly interacting amino acids. For these TCR, they argue, fewer amino acids contribute on average to the binding energy, triggering is more robust to mutations in the peptide sequence, and so the TCR is more cross-reactive. Thus again the argument emerges that the cross-reactivity is inversely related to the diversity of self driving selection. Kosmrlj et al. (98) employed this idea to put forward an explanation of why the population of elite-controllers of HIV infection is enriched for the HLA-B57 allele. Using a predictive peptide binding algorithm they argued that HLA-B\*5701 binds a lower diversity of self-peptides than average. Cytotoxic T cells restricted to this allele are then expected to be more cross-reactive than average and so are more resistant to virus mutations that might otherwise escape CTL control.

Chao et al. (93) and Kosmrlj et al. (95) took different approaches to the problem of how negative selection increases specificity. They came to the common conclusion that the most specific TCR are those with low to intermediate affinity to MHC – high enough to have a reasonable probability of passing positive selection, but low enough to avoid negative selection by allowing headroom for the additional contribution from the peptide component. The greater this headroom, the smaller the proportion of peptides that can trigger activation and so the greater the specificity.

# **THE EMERGENCE OF SPECIFICITY IN AVIDITY-BASED MODELS OF SELECTION**

Van den Berg et al. (99) developed a statistical framework to study the question of how specificity and self-tolerance can derive from a pre-selection repertoire of relatively promiscuous TCR. In their formalism, T cell activation is avidity-based and related to the rate of TCR triggering. Their starting point is that TCRs are degenerate and low affinity, binding weakly to many pMHC. TCR perceive an average signal derived from endogenous selfpMHC, and are triggered only by pMHC with sufficiently high prevalence and affinity to be visible above this background. The authors introduce the concept of an antigen presentation profile (APP), characterizing the abundances of different pMHC on antigen-presenting cells (APC). Positively-selected cells are selected against a given number of APC each with distinct APPs. In their framework, negative selection acts only on ubiquitous peptides presented on all APCs, and decisions are made on the basis of the entire APP of one APC. TCR that are triggered by this constitutive self-background are deleted. This filtering acts to sharpen the boundary between triggering rates, which give low and high activation probabilities, and so specificity can emerge even from a highly degenerate TCR. Interestingly they predict that negative selection does not have to be particularly stringent to generate an acceptably self-tolerant repertoire. Nevertheless in this model the selected repertoire may still be reactive to self-peptides expressed heterogeneously in the thymus, and in particular to peptides expressed at high levels only on certain cell types. Van den Berg and Rand (100) review avidity-based models of ligand discrimination.

# **ALLOREACTIVITY AND MHC RESTRICTION**

A high proportion (1–24%) of peripheral T cells are reactive to peptides presented on a foreign MHC allele (101–103), reflected clinically by acute T cell mediated rejection of grafts from MHCmismatched donors. These promiscuous "allogenic" responses contrast with the low precursor frequency (10−6–10−4 ) in normal immune responses to peptides presented by an individual's own MHC. Allogenic responses are also apparently counter to the notion of MHC restriction. Reconciling these results may tell us great deal about the relative contributions of peptide and MHC binding motifs to the TCR signals driving selection, and how this breakdown influences the coverage and cross-reactivity of the T cell repertoire.

Detours and Perelson (82) used their digit-string model of TCR–pMHC interactions, described above, to show how the probabilities of responsiveness to self and foreign MHC emerge. Mean alloreactivities of 1–2% arose naturally, at the lower end of the range of experimental estimates, and they showed that the alloreactivities of the pre- and post-selection repertoires are similar, as observed experimentally (17, 22). In essence, the modeling supports the hypothesis that the greater degree of alloreactivity than response frequency arises simply because many more pMHC ligands can be generated from one MHC than can be generated from one peptide (104). In other words, each TCR is triggered by ligands in a subset of pMHC shape space; one particular MHC along with its associated diversity of peptides will cover a far greater region of shape space than covered by one peptide and all the self-MHC alleles capable of presenting it; a given MHC will then stimulate far more of the T cell repertoire than will a given peptide.

They found that alloreactivity correlates with the extent of negative selection and inversely to the degree of MHC restriction. It can be seen intuitively how this emerges from their model. If negative selection is weak, positive selection must be correspondingly stringent in order to yield the selection efficiencies observed experimentally (3–5%). Stringent positive selection imposes an imprint of self-MHC on the repertoire – only those TCRs that bind strongly to self-MHC residues survive. The strength of binding to a randomly generated MHC not involved in selection (i.e., a foreign MHC) is then on average lower to that of self-MHC in the post-selection repertoire. This difference increases, and thus alloreactivity decreases, as the required strength of binding to self-MHC increases.

This trade-off between alloreactivity and restriction might be expected as they appear to be in conflict. However, experimental estimates of these two quantities are variable. The conclusions described above were derived analytically from a model capturing the mean-field behavior of the digit-string selection process, but did not deal with the variance in these measures of the repertoire outputs across specific simulations or experimental systems. The final study of the series (83) took a simulation-based approach, explicitly performing repertoire selection on random TCR and pMHC populations. This confirmed the inverse correlation between alloreactivity and MHC restriction and yielded sufficient variability to account for restriction ranging from absolute to partial in different settings.

Overall the digit-string model explored by Detours and colleagues yields remarkable agreement with many observations. Their model of TCR–pMHC binding is highly abstracted, but appears to be a powerful one. In part this might be because the relevant quantities for selection in their model are the minimum and maximum binding affinities that each TCR experiences during exposure to large samples of randomly generated pMHC strings. These two quantities will be drawn from extreme-value distributions, which should be insensitive to the distribution of binding strengths of randomly chosen TCR–pMHC pairs (70, 105). The additivity of the MHC and peptide contributions to the fatedetermining signal is likely the most questionable assumption, as the authors point out. Fate decisions may be based on the sum of several TCR interactions (which means for example that positive selection may occur though proximal binding of multiple low-affinity ligands) and so an avidity-based model may be more appropriate. Another caveat is that the population-average model

assumes that positive selection takes place on at most one MHC allele, which we will also return to.

# **INSIGHTS INTO FATE DETERMINATION MECHANISMS FROM STOCHASTICITY IN SELECTION**

Regulatory T cells (Treg) are a distinct lineage of CD4SP cells thought to lie at the higher end of the spectrum of acceptable self-reactivity and play a crucial role in the control of autoimmunity and tolerance to innocuous antigens. Many experimental studies of Treg development have shown that cells with the same TCR can develop into conventional and regulatory T cells within the same selecting environment [see, for example, Ref. (58, 106)], illustrating again, as represented in so many models, the stochastic nature of selection. There are at least two possible sources of this stochasticity. In a purely selective model precursors with identical TCR might be predisposed to the conventional or Treg fates through natural variation in expression of factors involved in lineage commitment. In a purely instructive model, cells within a clone are uncommitted, and intra-clonal heterogeneity in fate may derive from variation in the experience of each thymocyte during selection – most likely because each encounters a different random sample of self-peptides.

Bains et al. (107) used a probabilistic, instructive model that reflects this view of fate determination driven entirely by antigenic experience during selection, in conjunction with data from Ref. (58) to infer the number of pMHC binding events involved in fate determination. In that study, the numbers of conventional and Treg cells with a transgenically expressed TCR were measured for varying abundances of that TCR's agonist peptide on thymic epithelial cells. Conventional cell numbers declined monotonically with agonist abundance, while Treg increased and then decreased. Thus as agonist abundance increased, it appeared that T cells were initially diverted into the Treg lineage, before the risk of deletion through exposure to agonist dominated. Using this information and a simple graphical argument they were able to infer that fate decisions could not be affinity-driven (that is, made on the basis of a single pMHC interaction) unless TCR sensitivity varies during development, for which there is evidence [see Ref. (107) and references therein]. This model also explains apparently paradoxical observations regarding the effect of partial and full TCR agonists on the efficiency of Treg production (108).

### **THE LIMITS OF NEGATIVE SELECTION**

The potentially very large number of unique self-pMHC prompts the question of whether it is possible to tolerize thymocytes to all self-peptides within the timescale of thymic development. Müller and Bonhoeffer (109) studied this problem. Using constraints from the mouse proteome and the efficiencies of peptide production and binding to MHC, they estimated an upper limit of approximately 5 × 106 possible self-pMHC class I complexes. Notably, this diversity of self is several orders of magnitude lower than figures derived from the simple combinatoric arguments (64, 66) and is more closely aligned with an estimate that ~105 different nine-mers derived from the human proteome are expected to bind to one human MHC class I allele (73). The key quantity in Müller and Bonhoeffer's calculation is the probability *P* that a given selfpMHC is presented by any given APC in sufficient numbers for