negative selection to occur. The probability that a thymocyte specific for this (and only this) self-pMHC escapes negative selection is *P*E in their notation – distinct from the probability of immune escape discussed above – and *P*E = (1 − *P*) *n* , where *n* is the number of unique APC encountered during selection. In this model, *P*E is extremely sensitive to the number of copies *h* of a given self-pMHC that an APC needs to present in order to cause deletion – varying *h* between 15 and 1500 gives values of *P*E between 10−11 and 0.8. Favoring the higher estimates of *h*, Müller and Bonhoeffer (109) concluded that negative selection on the potential diversity of self is likely to be very leaky. Instead, they suggest thymic selection operates on a restricted subset of self-pMHC, a constraint imposed by the number of APCs encountered during selection. This requires that further tolerogenic mechanisms operate in the periphery to prevent autoimmune response to self antigens not encountered in the thymus (53, 70).

To support their argument, Müller and Bonhoeffer (109) reverted to the older model of cross-reactivity and selection to generate another estimate of the number of selecting ligands using the observed efficiency of negative selection. Recall that the probability of thymocyte with cross-reactivity *r* escaping negative selection on *Ns* unique selecting ligands is *P* = (1 − *r*) *Ns* ' *e* −*rNs* . Using the estimate of *r* = 2 × 10−5 (88), and *P* ' 0.33, they obtain *Ns* ' 105 unique selecting self-pMHC, or ~4% of the putative total number of self-pMHC. This estimate is consistent with those of Detours et al. (27). Both studies assume that this cross-reactivity *r* of thymocytes with self-pMHC is equal to the cross-reactivity of mature naive T cells to foreign pMHC. Since negative selection likely acts as a filter to reduce cross-reactivity in the post-selection repertoire (see above), this assumption is moot. But the need to meet the empirical constraint *e* −*rNs* ' 0.33 implies that higher values of *r* would reduce the number of unique selecting ligands *Ns* even further.

A subsequent exchange (110, 111) discussed the assumption that each TCR negatively selects only on a single self-pMHC ligand. Müller and Bonhoeffer (111) argued that in the Bernoulli trial model of cross-reactivity and selection, a 33% probability of survival implies that another third of all thymocytes were reactive to one self-pMHC only, giving some quantitative support to their original model. The discussion also addressed whether *Ns* is constrained by the residence time in the thymus or is a result of restricted presentation of self antigens. Müller and Bonhoeffer (111) favored the latter, presuming that evolution has optimized the thymic residence time for the purposes of efficient selection on a subset of self-peptides. More recently it has been argued that incomplete depletion of self-reactive cells in the thymus may be sufficient for robust self/non-self discrimination in the periphery, if interactions facilitating consensus between T cells are required for the initiation or suppression of immune responses (70).

## **OPTIMALITY OF INDIVIDUAL MHC DIVERSITY – CONSTRAINTS ARISING FROM THYMIC SELECTION**

The polymorphism of the MHC is huge, with hundreds of alleles identified at the HLA-A, HLA-B, and HLA-DR loci in humans (MHC is referred to as HLA in humans but hereon the term MHC is generally used, for simplicity). This diversification is thought not to have occurred by genetic drift but by two non-exclusive mechanisms. Heterozygote advantage (112, 113) suggests that individuals expressing more unique MHC alleles gain fitness by being able to present a larger array of pathogen peptides. Overall the evidence for heterozygote advantage in experimental models of infection is equivocal, though, and it has been argued with a quantitative model that this mechanism alone is insufficient to explain the extent of allelic diversity (114). Another theory is that MHC polymorphism is maintained by frequency-dependent selection under pathogen pressure, in which rare alleles confer protection against pathogen subversion of peptide presentation by commonly expressed alleles (115).

Intriguingly, individuals possess only a small proportion of all MHC alleles. Heterozygous humans possess six at the major HLA-A, HLA-B, and HLA-C loci, which code for MHC class I molecules that present peptides to CD8+ T cells, and six to eight at the HLA-DP, HLA-DQ, and HLA-DR MHC class II loci, which present to CD4+ T cells. A common explanation for this restricted withinindividual diversity is that it derives from the need to generate a broad, functional, and self-tolerant TCR repertoire in the thymus without excessive negative selection (116, 117). The qualitative argument is as follows. If *n* is the number of MHC alleles per person, then increasing *n* both increases the diversity of pathogenderived peptides that can be presented and increases the probability that a thymocyte will be able to obtain positively-selecting signals. On the other hand, higher *n* will also increase the range of self-peptides that can be presented. This will increase the stringency of negative selection, leading to inefficient generation of T cells in the thymus and potential gaps in the repertoire's coverage of peptide space. The observed number of different MHC molecules per individual may result from a trade-off between these demands.

The nature of MHC restriction needs to be considered carefully in these arguments. If restriction is absolute and each TCR recognizes only one MHC allele, increasing the number of alleles per person simply increases the size and diversity of the T cell repertoire with no cost because selection operates on each MHCrestricted subset of the pre-selection repertoire independently. In this case an upper limit to within-host MHC diversity might derive only from the need for APC to display sufficient numbers of peptides in conjunction with each MHC molecule to reliably mediate selection or immune activation. The trade-off evident in the qualitative argument above arises when MHC restriction is not absolute and thymocytes are capable of being positively and/or negatively selected on more than one allele.

Woelfing et al. (118) provide an excellent review of theoretical approaches to understanding intra-individual MHC diversity, but we outline the key results here. Nowak et al. (119) were the first to assess the qualitative trade-off argument using a mathematical model. In their analysis they defined *h* and *f* to be the proportions of T cells capable of being positively and negatively selected, respectively, by a given MHC allele. If an individual expresses *n* distinct MHC alleles, they argue that the proportion of the T cell repertoire surviving selection is

$$(1 - (1 - h)^n)(1 - f)^n.$$

The first term represents positive selection; (1 − *h*) *n* is the probability that a TCR fails to be selected by any MHC. The second term represents negative selection; (1 − *f*) *n* is the probability that a TCR is not negatively selected by any MHC. The proportion of the repertoire surviving is maximized at *n* = (1/*h*)log(1 + *h*/*f*). They argue that *h* ≤ *f*, supported by the experimental and modeling consensus is that positive selection is more stringent than negative selection. This gives *n*~1/*f*. However, using only the assumptions that *hn* 1, or that it is rare for a TCR to be positively selected on more than one MHC allele, and that the proportion of all peptides that can bind to a given MHC is 1, they calculate that *n* = 2/*f* maximizes the probability of a response to a randomly chosen foreign pMHC.

Borghans et al. (120) pointed out that this model contains an inconsistency, which allows for cells that fail to be positively selected on one MHC to be negatively selected by the same MHC. They denoted *p* and *n* to be the unconditional probabilities that one TCR is positively and negatively selected by a given MHC molecule. Then *n* < *p*, because the number of cells that fail negative selection on one MHC is necessarily smaller than the number that audition for it following positive selection on that same MHC. The proportion of the original repertoire that survives is then

$$(1 - n)^M - (1 - p)^M. \tag{4}$$

This model effectively lowers the stringency of negative selection expressed in Nowak et al. (119) and so reduces the cost of increasing the number of MHC alleles. They estimated the probabilities *p* and *n* were 0.01 and 0.005 respectively, using the known efficiencies of positive and negative selection in mice with known numbers of MHC alleles. The optimal value of *M* for these parameter values is far larger than observed allele numbers; conversely, asking what values of *p* and *n* correspond to the observed ranges of *M* being optimal leads to unrealistic levels of positive and negative selection. Their analysis therefore questions the trade-off hypothesis as an explanation of limited MHC diversity.

They suggest alternatives. They estimate that existing typical numbers of MHC alleles together with TCR cross-reactivity may be "good enough" for maximizing the probability of responding to a foreign peptide on self-MHC – in this case the selective pressure for increasing MHC alleles is weak or absent. Alternatively, increased numbers of MHC alleles may increase the risk of autoimmunity through cross-reactivity of T cells responding to antigen that have not been fully tolerized to self. Finally, limited numbers of MHC alleles may allow for sufficient densities of pMHC on the surface of antigen-presenting cells to be able to efficiently select and activate MHC-restricted T cells.

MHC restriction is not absolute in the models described above, although it holds approximately for positive selection when the per-allele positive selection probability *p* is small. However, there is evidence to suggest that MHC restriction is not manifest strongly at the positive selection stage. Zerrahn et al. (22) observed that a relatively large proportion of TCR still positively select when a single type of pMHC was expressed in the thymus. In that study, pre-selection TCRs had approximately a 5% chance of responding to a given class II MHC, independently for different alleles, validating one of the assumptions of these simple probabilistic selection models. On similar lines, Huseby et al. (92) found that the positively-selected repertoire contains TCR with a high degree

of cross-reactivity across MHC alleles, and suggested that MHC restriction emerges as a result of negative selection. Finally, the high degree of alloreactivity suggests that positive selection is only weakly MHC-restricted, and that failure to positive select reflects a generic inability to bind to MHC.

Motivated by this possibility, Woelfing et al. (118) revisited these probabilistic models. They assumed positive selection is highly degenerate with respect to MHC and that even very weak cross-reactivity with any allele is sufficient. Under this assumption, one of the presumed advantages of high MHC diversity is removed. Maximizing the probability of mounting an immune response, they estimated the optimal MHC diversity to be in a physiological range of 3–25.

Van den Berg and Rand (121) used a very different and sophisticated approach to the same optimality problem using a mechanistic, stochastic model of TCR triggering rather than the probabilistic repertoire-based models described above. Considering negative selection only, they concluded that limited individual MHC diversity is beneficial for self–non-self discrimination. The essence of their mathematical argument is that restricting the "diversity of foreign" is the key to increasing the signal-to-noise ratio for aTCR attempting to discriminate a foreign peptide from the background of self. This is achieved with a combination of limiting the number of MHC alleles each TCR can recognize (MHC restriction) and limiting the number of peptides that can be presented from one protein on one MHC allele ("peptide selectivity") to be roughly one. However, the need to ensure that every foreign protein is represented requires multiple MHC alleles, placing a theoretical lower bound on their number. An upper bound comes from the requirement that the density of relevant pMHC ligands must not fall too low on the surface of an APC, similar to the suggestion in Borghans et al. (120) – if a given pMHC is diluted by too many MHC, the relevant TCR will experience fluctuations in signaling that may reduce its discriminatory power. They conclude that of the order 10 MHC alleles is optimal. Notably, as in Ref. (118), this estimate arises without any constraints from positive selection.

# **SUMMARY**

This review has outlined how several relatively simple descriptions of single TCR–pMHC interactions have been used to understand aspects of TCR repertoire development. However, the discussion is necessarily incomplete. In particular, there is an extensive literature exploring the molecular mechanisms by which individual or collections of TCR discriminate between ligands of different affinities [see, for example, Ref. (100, 122–126)], which has direct relevance to thymic selection. It remains unclear how proximal TCR signals derived from multiple and diverse pMHC ligands can drive the emergence of specificity and MHC restriction in the post-selection repertoire, although the models of selection on ensembles of ligands have made steps in this direction (99, 121).

Many of the models discussed here assume that a single interaction above a minimum signaling threshold is sufficient for positive selection. However, there is evidence that repeated or sustained TCR signaling is required during the DP stage for positive selection to occur [see, for example, Ref. (17, 127)]. This may explain findings that positive and negative selection take place concurrently (46).

Overall it is remarkable how much insight into the quantitative aspects of thymic selection has emerged from highly abstracted models. However, there remain a lot of open areas for research, and many of the questions raised in the introduction are still unresolved. Regulatory T cell development in particular has received very little attention from modelers, and already it appears that the simplest extension to the simple probabilistic fixed-threshold model to include a fixed range of affinity or avidity for Treg selection is not sufficient to explain many experimental observations (107). The task of synthesizing and reconciling the huge diversity of experimental data related to thymic development is a daunting one, but the information available is perhaps currently underexploited by theorists.

# **ACKNOWLEDGMENTS**

The author thanks Charles Sinclair and the reviewers for helpful comments. This work was supported by the NIH (R01AI093870).

#### **REFERENCES**

chimeric mice as well as of normal mice respond to foreign antigens in the context of H-2 determinants not encountered on thymic epithelial cells. *Proc Natl Acad Sci U S A* (1980) **77**(12):7390–4. doi:10.1073/pnas.77.12.7390

**Conflict of Interest Statement:** The author declares 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: 10 September 2013; paper pending published: 23 October 2013; accepted: 09 January 2014; published online: 04 February 2014.*

*Citation: Yates AJ (2014) Theories and quantification of thymic selection. Front. Immunol. 5:13. doi: 10.3389/fimmu.2014.00013*

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

*Copyright © 2014 Yates. 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.*

# From pre-DP, post-DP, SP4, and SP8 thymocyte cell counts to a dynamical model of cortical and medullary selection

**Maria Sawicka1† , Gretta L. Stritesky 2† , Joseph Reynolds 1 , Niloufar Abourashchi 1 , Grant Lythe1 , Carmen Molina-París 1\* and Kristin A. Hogquist 2\***

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

2 Department of Laboratory Medicine and Pathology, Center for Immunology, University of Minnesota, Minneapolis, MN, USA

#### **Edited by:**

Ramit Mehr, Bar-Ilan University, Israel

#### **Reviewed by:**

Ramit Mehr, Bar-Ilan University, Israel Ruy Ribeiro, Los Alamos National Laboratory, USA Christian Schönbach, Nazarbayev University, Kazakhstan

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

Carmen Molina-París, Department of Applied Mathematics, School of Mathematics, University of Leeds, Leeds LS2 9JT, UK e-mail<EMAIL_ADDRESS>Kristin A. Hogquist, Department of Laboratory Medicine and Pathology, Center for Immunology, University of Minnesota, 2101 6th Street SE, Minneapolis, MN 55414, USA e-mail<EMAIL_ADDRESS>
†Maria Sawicka and Gretta L. Stritesky have contributed equally to this work.

Cells of the mature αβ T cell repertoire arise from the development in the thymus of bone marrow precursors (thymocytes). αβ T cell maturation is characterized by the expression of thousands of copies of identical αβ T cell receptors and the CD4 and/or CD8 co-receptors on the surface of thymocytes. The maturation stages of a thymocyte are: (1) double negative (DN) (TCR−, CD4− and CD8−), (2) double positive (DP) (TCR+, CD4+ and CD8+), and (3) single positive (SP) (TCR+, CD4+ or CD8+). Thymic antigen presenting cells provide the appropriate micro-architecture for the maturation of thymocytes, which "sense" the signaling environment via their randomly generated TCRs. Thymic development is characterized by (i) an extremely low success rate, and (ii) the selection of a functional and self-tolerant T cell repertoire. In this paper, we combine recent experimental data and mathematical modeling to study the selection events that take place in the thymus after the DN stage. The stable steady state of the model for the pre-DP, post-DP, and SP populations is identified with the experimentally measured cell counts from 5.5- to 17-week-old mice.We make use of residence times in the cortex and the medulla for the different populations, as well as recently reported asymmetric death rates for CD4 and CD8 SP thymocytes. We estimate that 65.8% of pre-DP thymocytes undergo death by neglect. In the post-DP compartment, 91.7% undergo death by negative selection, 4.7% become CD4 SP, and 3.6% become CD8 SP. Death by negative selection in the medulla removes 8.6% of CD4 SP and 32.1% of CD8 SP thymocytes. Approximately 46.3% of CD4 SP and 27% of CD8 SP thymocytes divide before dying or exiting the thymus.

**Keywords: thymocytes, negative selection, positive selection, death by neglect, mathematical model, steady state**

## **1. INTRODUCTION**

T cells are a major component of the adaptive immune system that play a crucial role in protection against a wide variety of pathogens. The T cell receptor (TCR) is generated by somatic recombination and has a vast potential to recognize foreign organisms. T cells do not recognize pathogens directly, but rather through binding pathogen fragments displayed by major histocompatibility complex (MHC) proteins on the surface of antigen presenting cells (APCs). Since MHC molecules are highly polymorphic, useful T cells must be selected for in each individual of the species. These T cells must have lineage specific effector functions that may include direct lysis, production of cytokines, and ability to regulate immune responses. Furthermore, some T cells have the potential to drive dangerous autoimmune responses (1). For all of these reasons, the development of a T cell repertoire is a highly specialized and tightly regulated process (2, 3). It takes place in a dedicated organ, the thymus, where unique properties of the microenvironment ensure the production of functional, yet self-tolerant T cells (4–6).

Multi-potent precursors travel from the bone marrow to the thymus through the blood. When they enter the thymus, the precursors that commit to the T cell lineage [or canonical early T cell progenitors (7)], after a 2 week period on average, transition from the double negative (DN) stage, where they do not express the coreceptors CD4 and CD8, to the double positive (DP) stage, where they express both co-receptors (6). At this stage, the majority of the cells have made productive TCR gene rearrangements and express a fully formed αβ TCR on the cell surface. DP cells are located in the cortex region of the thymus, where they use their TCR to survey self-peptides presented by MHCs on cortical thymic epithelial cells (cTECs) (8). DPs that recognize self-peptide-MHC complexes with low affinity undergo positive selection, whereas those with high affinity are deleted by negative selection (3). Those DP that fail to recognize self-peptide-MHC will undergo apoptosis in a process referred to as death by neglect. The DP cells that are positively selected will then transition to the single positive (SP) stage, where they express either the CD4 or CD8 co-receptor, depending upon their MHC class specificity (9). MHC class specificity also dictates gene expression changes that will ultimately determine the effector functions of that T cell: generally, cytotoxicity for CD8 T cells and cytokine production for CD4 T cells. All positively selected cells, whether MHC class I or class II specific, up-regulate the chemokine receptor CCR7, which facilitates their migration to the medulla, where they undergo further selection events. The medulla contains medullary epithelial cells (mTECs) that express tissue-restricted antigens regulated by the nuclear factor Aire (10).

Exposure to tissue-restricted antigens allows for further deletion of T cells specific for self-antigens they may encounter in the periphery. Finally, those cells that have been positively selected, yet have avoided negative selection, will mature and migrate to the periphery (11).

Previous efforts to develop mathematical models of thymic selection have been based on deterministic approaches or cellular automata simulations. These studies have shown the importance of (i) thymic antigen diversity on the size of the selected T cell repertoire (12), (ii) death rates for the more differentiated thymocyte subsets (13), (iii) thymocyte proliferation and residence times (14), (iv) epithelial networks for thymocyte development and migration (15), (v) thymocyte competition for antigen (16), (vi) self-pMHC complexes expressed on dendritic cells (DCs) (17), (vii) receptor–ligand binding affinity (18), and (viii) a sharp threshold in TCR-ligand binding affinity that defines the boundary between negative and positive selection (19). Recent work by Ribeiro and Perelson (20) supports the need to develop appropriate mathematical models to interpret T cell receptor excision circles (or TREC) data,which are used to quantify thymic export (20). Sinclair et al. in Ref. (21) bring together experimental immunology with mathematical modeling to conclude that CD8 precursor thymocytes are more susceptible to death than CD4 precursors. This asymmetry in the death rates underlies the experimentally observed CD4:CD8 T cell ratio in the periphery.

Previous experimental studies have tried to determine the number of cells going through positive and negative selection in the thymus. However, reports estimating the relative number of cells undergoing negative selection compared to positive selection have been widely variable. Some find that very few cells undergo negative selection; others find that two times more cells undergo negative selection than positive selection (22–25). In this report, we develop a deterministic mathematical model of T cell development in the thymus. Some of us recently published a report where we used a novel approach (Bim−/−Nur77GFP mice) that allowed us to calculate the number of cells undergoing positive and negative selection (26). Using previously published data on the relative life-span of DP and SP cells, we estimated the hourly rate of both positive and negative selection (26). In this manuscript, we make use of (i) a subset of this experimental data, and (ii) the asymmetric death rates observed for CD4 and CD8 precursor thymocytes (21), to develop two mathematical models that will enable us to estimate selection rates in the cortex and the medulla, and provide a quantitative measure for the stringency of thymic selection. The first model (see Section 2.1) allows the identification of the following parameters: DN thymocyte influx into the cortex, pre-DP and post-DP death rates, and pre-DP and post-DP differentiation rates. Under the assumption of asymmetric death rates for the CD4 and CD8 SP thymocytes (21), we extend the first model to provide estimates for the following medullary rates (see Section 2.2): CD4 and CD8 SP death, proliferation, and maturation rates.

### **2. MATERIALS AND METHODS**

#### **2.1. A FIRST MODEL OF THYMIC DEVELOPMENT AFTER THE DN STAGE**

In this section, we introduce a deterministic model of thymocyte development after the DN stage. This first model will be required to calibrate the parameter values of the second model introduced in Section 2.2. In particular, and as described in Section 3.2, the first model allows the identification of parameter values for the following rates: φ, µ1, µ2, ϕ1, and ϕ2.

This mathematical model makes use of a data set obtained from the analysis of eight C57BL/6 wild type and Bim deficient mice (average age 9 weeks), that express a Nur77GFP transgene to indicate TCR signal strength experience (26). Flow cytometric analysis, as described in that study, used standard markers to define various stages of T cell development in the thymus. The Nur77GFP reporter and Bim deficiency were novel modifications that allowed us to quantify cells that normally would be deleted by strong TCR signaling. In the mathematical model, we consider the following thymocyte populations: *n*1, the population of preselection DP thymocytes (double positive), that are TCRβ low and CD69low (26), *n*2, the population of post-selection DP thymocytes, that are TCRβ + and CD69high (26), and *n*3, the population of mature SP (single positive) thymocytes.

We assume that DN thymocytes differentiate to become preselection DP thymocytes with rate (cells/day) φ.Wefurther assume that after the DN stage, thymocyte cell fate is determined by the TCR signal, which a given thymocyte has received. Sinclair et al. used CFSE labeling to show that there is no proliferation at the post-DP stage (see Figure A1 of their manuscript) (21). Stritesky et al. looked at proliferation in the post-DP pool with BrdU labeling, and found no evidence (26). We have, thus, only included proliferation in the SP thymocyte population (21, 26). The three populations, *n*1, *n*2, and *n*3, are involved in the following selection events in the cortex and the medulla (see **Figure 1**):

The time evolution of the three populations can be described by the following set of ordinary differential equations (ODEs), which are based on the selection events described above:

$$\begin{cases} \frac{dn\_1}{dt} = \phi - \varphi\_1 n\_1 - \mu\_1 n\_1 \\ \frac{dn\_2}{dt} = \varphi\_1 n\_1 - \varphi\_2 n\_2 - \mu\_2 n\_2 \\ \frac{dn\_3}{dt} = \varphi\_2 n\_2 - \varphi\_3 n\_3 - \mu\_3 n\_3 + \lambda\_3 n\_3 \end{cases} \tag{1}$$

We are interested in studying the steady state of these populations, as the experimental data correspond to population cell numbers in the three stages (pre-DP, post-DP, and SP) for a steady state thymus (26). The steady state of the system of equations [equation (1)] is given by:

$$n\_1^\* = \frac{\phi}{\varphi\_1 + \mu\_1}, \quad n\_2^\* = \frac{n\_1^\* \varphi\_1}{\varphi\_2 + \mu\_2}, \quad n\_3^\* = \frac{n\_2^\* \varphi\_2}{\varphi\_3 + \mu\_3 - \lambda\_3}. \tag{2}$$

This unique steady state exists if and only if ϕ3 + µ3 − λ3 > 0, so that we have *n* ∗ 3 > 0. In order to study the linear stability of the steady state, we calculate *A*, the Jacobian matrix of equation (1), as follows:

$$A = \begin{pmatrix} -(\varphi\_1 + \mu\_1) & 0 & 0\\ \varphi\_1 & -(\varphi\_2 + \mu\_2) & 0\\ 0 & \varphi\_2 & -(\varphi\_3 + \mu\_3 - \lambda\_3) \end{pmatrix} . \tag{3}$$

*A* is also the Jacobian matrix at the steady state *n* ∗ = (*n* ∗ 1 , *n* ∗ 2 , *n* ∗ 3 ), as the system of ODEs [equation (1)] is linear. The three eigenvalues of *A* are given by (as the matrix is lower triangular):

$$\beta\_1 = -\left(\wp\_1 + \mu\_1\right), \quad \beta\_2 = -\left(\wp\_2 + \mu\_2\right), \quad \beta\_3 = -\left(\wp\_3 + \mu\_3 - \lambda\_3\right).$$

Therefore, the steady state [equation (2)] is stable, if and only if, ϕ3 + µ3 − λ3 > 0, which is also the condition for its existence. We conclude this section with the analytical solution of the system of ODEs [equation (1)], given initial conditions, which provides the time evolution of the three thymocyte populations:

$$\begin{aligned} n\_1(t) &= n\_1^\* + n\_1(0) \, e^{-(\varphi\_1 + \mu\_1)t}, \\ n\_2(t) &= n\_2^\* + \frac{n\_1(0) \, \varphi\_1}{[(\varphi\_2 + \mu\_2) - (\varphi\_1 + \mu\_1)]} \, e^{-(\varphi\_1 + \mu\_1)t} \\ &+ n\_2(0) \, e^{-(\varphi\_2 + \mu\_2)t}, \end{aligned}$$

$$n\_3(t) = n\_3^\* + \frac{n\_1(0)\,\varphi\_1}{[(\varphi\_2 + \mu\_2) - (\varphi\_1 + \mu\_1)]} \tag{4}$$

$$\begin{split} \times \frac{\varphi\_2}{[(\varphi\_3 + \mu\_3 - \lambda\_3) - (\varphi\_1 + \mu\_1)]} \, e^{-(\varphi\_1 + \mu\_1)t} \\ + \frac{n\_2(0)\,\varphi\_3}{[(\varphi\_3 + \mu\_3 - \lambda\_3) - (\varphi\_2 + \mu\_2)]} \, e^{-(\varphi\_2 + \mu\_2)t} \\ + n\_3(0)\,\, e^{-(\varphi\_3 + \mu\_3 - \lambda\_3)t} \,, \end{split}$$

where *n*1(0), *n*2(0), *n*3(0) represent the initial conditions for the thymocyte populations. Note that in the late time limit, that is, if *t*→ +∞ and ϕ3 + µ3 − λ3 > 0, then *n*1(*t*) → *n* ∗ 1 , *n*2(*t*) → *n* ∗ 2 and *n*3(*t*) → *n* ∗ 3 , as it is the unique stable steady state.

#### **2.2. A SECOND MODEL OF THYMIC DEVELOPMENT AFTER THE DN STAGE: CD4 AND CD8 SP THYMOCYTES**

As described in Section 2.1, the first deterministic model will allow us to calibrate some of the parameters of a more comprehensive model, which we now introduce. We subdivide the SP thymocyte population in two classes: CD4 SP and CD8 SP thymocytes. This is an extension of the model introduced in the previous section, and is motivated by the fact that experimentally, SP thymocytes express either the CD4 or the CD8 co-receptor. We now have four different thymocyte populations to consider: *n*1, the population of pre-selection DP (double positive) thymocytes, *n*2, the population of post-selection DP thymocytes, *n*4, the population of mature CD4+ SP (single positive) thymocytes, and *n*8, the population of mature CD8+ SP (single positive) thymocytes.

As described in the previous section, we assume that DN thymocytes differentiate to become pre-selection DP thymocytes with rate (cells/day) φ, and that after the DN stage, thymocyte cell fate is determined by the TCR signal, which a given thymocyte has received. Thus, the four populations, *n*1, *n*2, *n*4, and *n*8, with *n*3 = *n*4 + *n*8, are involved in the following selection events in the cortex and the medulla (see **Figure 2**):

We assume that all model parameters are positive, that is, µ1, µ2, µ4, µ8, ϕ1, ϕ2, ϕ4, ϕ8, ξ 4, ξ 8, φ, λ4, λ8 > 0, and note that the parameters and thymocyte populations of the first and second model are related by the following equations:

$$\begin{aligned} \varphi\_2 &= \varphi\_4 + \varphi\_8 \ , \quad \xi\_4 \ n\_4 + \xi\_8 \ n\_8 = \varphi\_3 \ n\_3 \ ,\\ \mu\_4 \ n\_4 + \mu\_8 \ n\_8 &= \mu\_3 \ n\_3 \ , \quad \lambda\_4 \ n\_4 + \lambda\_8 \ n\_8 = \lambda\_3 \ n\_3 \ . \end{aligned} \tag{5}$$

The time evolution of the four populations can be described by the following set of ODEs:

$$\begin{cases} \frac{dn\_1}{dt} = \phi - \varphi\_1 n\_1 - \mu\_1 n\_1\\ \frac{dn\_2}{dt} = \varphi\_1 n\_1 - (\varphi\_4 + \varphi\_8) n\_2 - \mu\_2 n\_2\\ \frac{dn\_4}{dt} = \varphi\_4 n\_2 - \xi\_4 n\_4 - \mu\_4 n\_4 + \lambda\_4 n\_4\\ \frac{dn\_8}{dt} = \varphi s n\_2 - \xi s n s - \mu s n s + \lambda s n s \end{cases} \tag{6}$$

(µ4) or (µ8), or proliferation (λ4) or (λ8).

We are interested in studying the steady state of these populations, as the experimental data correspond to population cell numbers in the four stages (pre-DP, post-DP, CD4 SP, and CD8 SP) for a steady state thymus (26). The steady state of the system of equations [equation (6)] is given by:

$$n\_1^\* = \frac{\phi}{\varphi\_1 + \mu\_1}, \quad n\_2^\* = \frac{n\_1^\* \varphi\_1}{\varphi\_4 + \varphi\_8 + \mu\_2},$$

$$n\_4^\* = \frac{n\_2^\* \varphi\_4}{\xi\_4 + \mu\_4 - \lambda\_4}, \quad n\_8^\* = \frac{n\_2^\* \varphi\_8}{\xi\_8 + \mu\_8 - \lambda\_8}.\tag{7}$$

This unique steady state exists if and only if ξ 4 + µ4 − λ4 > 0 and ξ 8 + µ8 − λ8 > 0, so that we guarantee *n* ∗ 4 > 0 and *n* ∗ 8 > 0. In order to study the linear stability of the steady state, we calculate *B*, the Jacobian matrix of equation (6), as follows:

$$B = \begin{pmatrix} -(\varphi\_1 + \mu\_1) & 0 & 0 & 0 \\ \varphi\_1 & -(\varphi\_4 + \varphi\_8 + \mu\_2) & 0 & 0 \\ 0 & \varphi\_4 & -(\xi\_4 + \mu\_4 - \lambda\_4) & 0 \\ 0 & \varphi\_8 & 0 & -(\xi\_8 + \mu\_8 - \lambda\_8) \end{pmatrix} \tag{8}$$

*B* is also the Jacobian at the steady state *n* ∗ = (*n* ∗ 1 , *n* ∗ 2 , *n* ∗ 4 , *n* ∗ 8 ), as the system of ODEs is linear. The four eigenvalues of *B* are given by:

$$\begin{aligned} \beta\_1 &= - (\varphi\_1 + \mu\_1), & \beta\_2 &= - (\varphi\_4 + \varphi\_8 + \mu\_2) \end{aligned}$$

$$\begin{aligned} \beta\_3 &= - (\xi\_4 + \mu\_4 - \lambda\_4) \end{aligned} \quad \beta\_4 = - (\xi\_8 + \mu\_8 - \lambda\_8) \dots$$

Therefore, the steady state equation (6) is stable if and only if ξ 4 + µ4 − λ4 > 0 and ξ 8 + µ8 − λ8 > 0, which is also the condition for its existence. We conclude this section with the analytical solution of the system of ODEs equation (6), given initial conditions, which provides the time evolution of the four thymocyte populations:

$$\begin{aligned} n\_1(t) &= n\_1^\* + n\_1(0) \ e^{-(\varphi\_1 + \mu\_1)t}, \\ n\_2(t) &= n\_2^\* + \frac{n\_1(0)\,\varphi\_1}{[(\varphi\_4 + \varphi\_8 + \mu\_2) - (\varphi\_1 + \mu\_1)]} \, e^{-(\varphi\_1 + \mu\_1)t} \\ &+ n\_2(0) \, e^{-(\varphi\_4 + \varphi\_8 + \mu\_2)t}, \end{aligned}$$

$$\begin{split} n\_{4}(t) &= n\_{4}^{\*} + \frac{n\_{1}(0)\,\varphi\_{1}}{[(\varphi\_{4} + \varphi\_{8} + \mu\_{2}) - (\varphi\_{1} + \mu\_{1})]} \\ &\quad \times \frac{\varphi\_{4} + \varphi\_{8}}{[(\xi\_{4} + \mu\_{4} - \lambda\_{4}) - (\varphi\_{1} + \mu\_{1})]}e^{-(\varphi\_{1} + \mu\_{1})t} \\ &\quad + \frac{n\_{2}(0)\,\xi\_{4}}{[(\xi\_{4} + \mu\_{4} - \lambda\_{4}) - (\varphi\_{4} + \varphi\_{8} + \mu\_{2})]}e^{-(\varphi\_{4} + \varphi\_{8} + \mu\_{2})t} \\ &\quad + n\_{4}(0)\,e^{-(\xi\_{5} + \mu\_{4} - \lambda\_{4})t}, \\ n\_{8}(t) &= n\_{8}^{\*} + \frac{n\_{1}(0)\,\varphi\_{1}}{[(\varphi\_{4} + \varphi\_{8} + \mu\_{2}) - (\varphi\_{1} + \mu\_{1})]} \\ &\quad \times \frac{\varphi\_{4} + \varphi\_{8}}{[(\xi\_{8} + \mu\_{8} - \lambda\_{8}) - (\varphi\_{1} + \mu\_{1})]}e^{-(\varphi\_{1} + \mu\_{1})t} \\ &\quad + \frac{n\_{2}(0)\,\xi\_{8}}{[(\xi\_{8} + \mu\_{8} - \lambda\_{8}) - (\varphi\_{4} + \varphi\_{8} + \mu\_{2})]}e^{-(\varphi\_{4} + \varphi\_{8} + \mu\_{2})t} \\ &\quad + n\_{8}(0)\,e^{-(\xi\_{8} + \mu\_{8} - \lambda\_{8})t}, \end{split}$$

where *n*1(0), *n*2(0), *n*4(0), *n*8(0) represent the initial conditions for the thymocyte populations. Note that in the late time limit, that is, if *t*→ +∞ and ξ 4 + µ4 − λ4 > 0 and ξ 8 + µ8 − λ8 > 0, then *n*1(*t*) → *n* ∗ 1 , *n*2(*t*) → *n* ∗ 2 , *n*4(*t*) → *n* ∗ 4 , and *n*8(*t*) → *n* ∗ 8 , as it is the unique stable steady state.

#### **3. RESULTS**

#### **3.1. PARAMETER ESTIMATION FOR THE FIRST MODEL (MEANS)**

In this section, we make use of previously published experimental data (26) that provide thymocyte cell counts for the three subsets considered in the first model: pre-DPs, post-DPs, and SPs. The original experiments have been carried out for two types of mice: wild type mice (Bim+/−) and Bim deficient mice (Bim−/−). In this paper, we will only be considering the wild type experimental results. The data will allow us to provide experimental estimates for the steady state thymocyte cell counts: *n* ∗ 1 , *n* ∗ 2 , *n* ∗ 3 . Note that, in order to estimate rates (with units of inverse time), thymocyte cell counts are not enough. Thus, we will make use of the additional knowledge provided by experimentally determined residence times for each population, τ *i* , with *i* = 1, 2, 3. If we make use of the model (see Section 2.1), the residence time in compartment *i* can be expressed as:

$$
\pi\_i = \frac{1}{\varphi\_i + \mu\_i}, \quad \text{for} \quad i \in \{1, 2, 3\} \;.
$$

The experimental data (see **Table 1**) correspond to the number of cells (thymocytes) at steady state (26), in each of the thymic compartments considered in the mathematical models (see Sections 2.1 and 2.2), and for eight different mice (*j* = 1, 2, . . ., 8).

We have made use of the following average residence times in each compartment (27–29)

$$\begin{aligned} \mathbf{r}\_1 &= 60 \,\mathrm{h} = 2.5 \,\mathrm{days} \,, \quad \mathbf{r}\_2 = 16 \,\mathrm{h} = 0.67 \,\mathrm{days} \,, \\\mathbf{r}\_3 &= 96 \,\mathrm{h} = 4 \,\mathrm{days} \,\mathrm{.} \end{aligned}$$

In order to derive estimates for the model parameters, we have carried out the following steps:

Let us introduce *a*1 by the following equation, *a*1 = *n* ∗ 2 *n* ∗ 1 , and make use of the experimental data to write: *n* ∗,*i* 2 = *a n* ∗,*i* 1 + *i* , for *i* = 1, 2, . . ., 8. Thus, the squared error is given by:

$$E(a\_1) = \sum\_{i=1}^{8} \left( n\_2^{\*,i} - a\_1 \ n\_1^{\*,i} \right)^2$$

.

We minimize *E*(*a*1) with respect to *a*1, that is *dE da*1 = 0. Solving for *a*1, we obtain:

$$a\_1 = \frac{\sum\_{i=1}^{8} n\_1^{\*,i} \ n\_2^{\*,i}}{\sum\_{i=1}^{8} \left(n\_1^{\*,i}\right)^2} \ .$$

**Table 1 | Experimental steady state thymocyte cell counts for the wild type pre-DP, post-DP, CD4 SP, and CD8 SP populations.**

The bold font highlights the mean and the standard deviation from the individual mice data.

Given *a*1, we can then estimate ϕ1 from the equation ϕ1 = *a*1 τ2 .

We obtain the following parameter values:

$$\begin{aligned} \phi &= 35.350 \times 10^6 \text{ cells/day}, \ \mu\_1 = 0.263 \text{ day}^{-1}, \\ \mu\_3 &= 0.099 \text{ day}^{-1}, \ \varphi\_1 = 0.137 \text{ day}^{-1}, \ \varphi\_3 = 0.151 \text{ day}^{-1}, \ \mu\_2 \end{aligned}$$

and

$$\begin{aligned} \mu\_2 &\in \{1.295, 1.443\} \text{ day}^{-1}, \ \varphi\_2 \in \{0.050, 0.198\} \text{ day}^{-1}, \\\ \lambda\_3 &\in \{0.143, 0.250\} \text{ day}^{-1}. \end{aligned}$$

These parameters imply the following thymic selection rates:

#### **3.1.1.  Differentiation rates**

5.0 × 105 cells/h are positively selected in compartment 1, that is, become post-DP from pre-DP (ϕ1 *n* ∗ 1 ), 4.4 × 104 cells/h are positively selected in compartment 2, that is, become SP from post-DP (ϕ2 *n* ∗ 2 ), and 1.0 × 105 cells/h leave compartment 3 to go to the periphery (ϕ3 *n* ∗ 3 ).

#### **3.1.3. Proliferation rate**

1.3 × 105 cells/h proliferate in compartment 3 (λ3 *n* ∗ 3 ).

We have also computed the stringency of thymic selection, which we define as given by the ratio:

$$\frac{\varphi\_3 \ n\_3^\*}{\phi} = 6.799\% \ . $$

Finally,we have computed the (per cell) probability to die, given that the cell is in compartment *i* (*i* = 1, 2, 3), as well as the (per cell) probability to proliferate in the medulla. We have obtained:

$$\begin{aligned} p\_1 &= \frac{\mu\_1}{\mu\_1 + \varphi\_1} = 65.8\%, & p\_2 &= \frac{\mu\_2}{\mu\_2 + \varphi\_2} = 91.7\%, \\ p\_3 &= \frac{\mu\_3}{\mu\_3 + \varphi\_3 + \lambda\_3} = 22.9\%, & q\_3 &= \frac{\lambda\_3}{\mu\_3 + \varphi\_3 + \lambda\_3} = 42.2\% \end{aligned}$$

#### **3.2. PARAMETER ESTIMATION FOR THE SECOND MODEL (MEANS)**

In this section, we make use of previously published experimental data (26) that provide thymocyte cell counts for the four subsets considered: pre-DPs, post-DPs, CD4 SPs, and CD8 SPs. We only make use of the experimental data for the wild type mice. The data will allow us to provide experimental estimates for the steady state thymocyte cell counts: *n* ∗ 1 , *n* ∗ 2 , *n* ∗ 4 , *n* ∗ 8 . As described in Section 3.1, we also need residence times for each population subset, τ *i* , with *i* = 1, 2, 4, 8. If we make use of the model (see Section 2.2), the residence time in compartment *i* can be expressed as:

$$\begin{aligned} \pi\_i &= \frac{1}{\varphi\_i + \mu\_i}, \quad \text{for} \quad i \in \{1, 2\}, \quad \text{and} \\\pi\_i &= \frac{1}{\xi\_i + \mu\_i}, \quad \text{for} \quad i \in \{4, 8\} \end{aligned}$$

Recent experimental data provide support for asymmetric death rates in the CD4 and CD8 SP compartments (21). The estimated death rates for CD4 and CD8 SP thymocytes are1 µ4 = 0.04 day−1 and µ8 = 0.11 day−1 . We also make use of the estimates derived in Section 3.1 for φ, µ1, µ2, ϕ1, and ϕ2. Finally,

1Private communication from Ben Seddon and Andy Yates.

the average residence times in each compartment, as described in Section 3.1, are given by:

$$\begin{aligned} \mathbf{r}\_1 &= 60 \text{ h} = 2.5 \text{ days} \,, \quad \mathbf{r}\_2 = 16 \text{ h} = 0.67 \text{ days} \,, \\ \mathbf{r}\_4 &= 96 \text{ h} = 4 \text{ days} \,, \quad \mathbf{r}\_8 = 96 \text{ h} = 4 \text{ days} \, . \end{aligned}$$

In order to derive estimates for the model parameters, we have carried out the following steps:

[0, ϕ2]. As discussed in Section 3.1, we only accept parameter values which provide thymocyte cell counts at time *t* = 21 days that are within one standard deviation from the experimentally determined values (see **Table 1**).

We obtain the following parameter values:

$$\begin{aligned} \mu\_4 &= 0.04 \text{ day}^{-1}, & \mu\_8 &= 0.11 \text{ day}^{-1}, \\ \xi\_4 &= 0.21 \text{ day}^{-1}, & \xi\_8 &= 0.14 \text{ day}^{-1}, \end{aligned}$$

and

$$\begin{aligned} \varphi\_4 &= 0.070 \text{ day}^{-1}, & \varphi\_8 &= 0.054 \text{ day}^{-1}, \\ \lambda\_4 &= 0.216 \text{ day}^{-1}, & \lambda\_8 &= 0.093 \text{ day}^{-1}. \end{aligned}$$

These parameters imply the following thymic selection rates:

# **3.2.1. Death rates**

2.05 × 104 cells/h die by negative selection in compartment 4 (µ4 *n* ∗ 4 ) and 1.95 × 104 cells/h die by negative selection in compartment 8 (µ8 *n* ∗ 8 ).

## **3.2.2. Differentiation rates**

2.50 × 104 cells/h are CD4 positively selected in compartment 2, that is, become CD4 SP from post-DP (ϕ4 *n* ∗ 2 ), 1.90 × 104 cells/h are CD8 positively selected in compartment 2, that is, become CD8 SP from post-DP (ϕ8 *n* ∗ 2 ), 1.08 × 105 cells/h leave compartment 4 to go to the periphery (ξ4 *n* ∗ 4 ), and 2.48 × 104 cells/h leave compartment 8 to go to the periphery (ξ8 *n* ∗ 8 ).

### **3.2.3. Proliferation rates**

11.10 × 104 cells/h proliferate in compartment 4 (λ4 *n* ∗ 4 ) and 1.60 × 104 cells/h proliferate in compartment 8 (λ8 *n* ∗ 8 ).

We can compute the stringency of thymic selection, defined by the ratio:

$$\frac{\xi\_4 \ n\_4^\* + \xi\_8 \ n\_8^\*}{\phi} = 8.969\phi\ .$$

We can provide an estimate for the cortical positive selection probabilities, that is the (per post-DP cell) probability to become a CD4 SP or a CD8 SP, and the probability to be negatively selected in the cortex. We have obtained:

$$\begin{aligned} s\_4 &= \frac{\varphi\_4}{\mu\_2 + \varphi\_4 + \varphi\_8} = 4.7\%, \quad s\_8 = \frac{\varphi\_8}{\mu\_2 + \varphi\_4 + \varphi\_8} = 3.6\%, \\\ p\_2 &= \frac{\mu\_2}{\mu\_2 + \varphi\_4 + \varphi\_8} = 91.7\% \ . \end{aligned}$$

Finally,we have computed the (per cell) probability to die, given that the cell is in compartment *i*, as well as the (per cell) probability to proliferate in the medulla. We obtain:

$$p\_4 = \frac{\mu\_4}{\mu\_4 + \xi\_4 + \lambda\_4} = 8.6\%, \ q\_4 = \frac{\lambda\_4}{\mu\_4 + \xi\_4 + \lambda\_4} = 46.3\%,$$

$$p\_8 = \frac{\mu\_8}{\mu\_8 + \xi\_8 + \lambda\_8} = 32.1\%, \ q\_8 = \frac{\lambda\_8}{\mu\_8 + \xi\_8 + \lambda\_8} = 27.0\% \text{ .}$$

These probabilities imply that the probability to exit the thymus as a mature CD4 thymocyte (that has already reached the medulla) is given by 100−(8.6+46.3) 100 , which is 45.1%, and the probability to exit as a mature CD8 thymocyte (that has already reached the medulla) is given by 100−(32.1+27.0) 100 , which is 40.9%.

#### **3.3. SENSITIVITY ANALYSIS**

In this section, we explore the sensitivity of the parameters to perturbations in the experimental data. For the first model, the experimental data are given in terms of the following eight quantities:

$$\theta = (\mathfrak{r}\_1, \mathfrak{r}\_2, \mathfrak{r}\_3, \phi\_{\text{out}}, a\_1, a\_2, \vec{n}\_3, \vec{n}\_1) \,,$$

where *a*1, *a*2 are the coefficients of the linear regression of *n* ∗ 2 *n* ∗ 1 and *n* ∗

3 *n* ∗ 2 , respectively, and *n*¯*i* is the experimental mean value of *ni* .

We perturb each entry of the vector θ by adding and subtracting 10% of its value. Therefore, we now have two values for θ*i* , equal to θ*i* + 1 10 θ*i* and θ*i* − 1 10 θ*i* . Consequently, we have 2 8 sets of θ, which will be used to compute the corresponding model parameters as described in Section 3.1. Parameter values will only be accepted if they provide a stable solution before *t* = 21 days.

For the second model, the experimental data is given in terms of the following seven quantities:

$$\theta = (\mathfrak{r}\_{4}, \mathfrak{r}\_{8}, \mu\_{4}, \mu\_{8}, a\_{3}, a\_{4}, a\_{5}),$$

where *a*3, *a*4, *a*5 are the coefficients of the linear regression of *n* ∗ 2 *n* ∗ 4 , *n* ∗ 4 *n* ∗ 8 , and *n* ∗ 3 *n* ∗ 8 , respectively. We have made use of the means of the following parameters of the first model: φ*,* ϕ1, µ1, ϕ2, µ2, λ3.

**Table 2 | Means, 95% trimmed and minimum–maximum intervals of the model parameters.**

We perturb each entry of the vector θ as described above. Therefore, we have 27*n*ϕ4 sets of θ, with *n*ϕ4 , the number of different values considered for ϕ4 in the interval (0, ϕ2). Parameter values will only be accepted if they provide a stable solution before *t* = 21 days.

The results of the sensitivity analysis, with 95% trimmed intervals2 and minimum–maximum interval ranges, are given in **Table 2**.

#### **3.4. VARIABILITY IN THE SELECTION RATES**

The (trimmed and minimum–maximum) intervals derived in Section 3.3 allow us to estimate the variability in the different selection rates discussed in Sections 3.1 and 3.2. For example, given variations in the parameters, the corresponding variations in the selection rates can be shown to be:

$$\Delta p\_i = \frac{1}{\left(\mu\_i + \varphi\_i\right)^2} (\varphi\_i \Delta \mu\_i + \mu\_i \Delta \varphi\_i) \quad \text{for} \quad i = 1, 2, \dots \qquad (10)$$

$$\Delta \alpha\_i = \frac{1}{\left[\left(\mu\_{\alpha\_i} + \mu\_{\alpha\_i}\right) \Delta \mu\_{\alpha\_i} + \mu\_{\alpha\_i} \Delta \varphi\_{\alpha\_i} + \mu\_{\alpha\_i} \Delta \lambda\_{\alpha\_i}\right]}$$

$$
\Delta p\_3 = \frac{1}{\left(\mu\_3 + \varphi\_3 + \lambda\_3\right)^2} \left[ (\varphi\_3 + \lambda\_3)\Delta\mu\_3 + \mu\_3\Delta\varphi\_3 + \mu\_3\Delta\lambda\_3 \right],
\tag{11}
$$

$$
\Delta q\_3 = \frac{1}{\left(\mu\_3 + \varphi\_3 + \lambda\_3\right)^2} \left[\lambda\_3 \Delta \mu\_3 + \lambda\_3 \Delta \varphi\_3 + (\mu\_3 + \varphi\_3) \Delta \lambda\_3\right],
\tag{12}
$$

$$\begin{split} \Delta s\_{i} &= \frac{1}{\left(\mu\_{2} + \varphi\_{i} + \varphi\_{j}\right)^{2}} [\varphi\_{i} \Delta \mu\_{2} + \varphi\_{i} \Delta \varphi\_{j} \\ &+ (\mu\_{2} + \varphi\_{j}) \Delta \varphi\_{i}] \quad \text{for} \quad i = 4, j = 8 \quad \text{or} \quad i = 8, j = 4, \dots \end{split} \tag{13}$$

$$\Delta p\_{i} = \frac{1}{\left(\mu\_{i} + \xi\_{i} + \lambda\_{i}\right)^{2}} [\mu\_{i}\Delta\xi\_{i} + \mu\_{i}\Delta\lambda\_{i}]$$

$$+ (\xi\_{i} + \lambda\_{i})\Delta\mu\_{i}] \quad \text{for} \quad i = 4, 8 \,, \tag{14}$$

2We define the 95% trimmed interval to be the result of the sensitivity analysis after trimming the lower and upper 2.5% of values.

$$
\Delta q\_i = \frac{1}{\left(\mu\_i + \xi\_i + \lambda\_i\right)^2} [\lambda\_i \Delta \xi\_i + \lambda\_i \Delta \mu\_i]
$$

$$
+ (\xi\_i + \mu\_i) \Delta \lambda\_i \left[ \begin{array}{c} \text{for} \\ i = 4, 8 \end{array} \right.\tag{15}
$$

We present in **Table 3** the variability of the selection rates.

# **4. DISCUSSION**

We have brought together experimental data with a mathematical compartment model [similar to other progression models of CD4 and CD8 T cell development (13, 14, 18, 21, 33)] to provide estimates for the selection events that take place in the thymus. We have made use of a range of experimental data: (i) steady state thymocyte cell counts (26), mean residence times in each compartment (27–29), murine thymic export rate (14, 26, 30), and recently reported asymmetric death rates for the CD4 SP and CD8 SP thymocytes (21). Our preliminary results support the unexpectedly high death rate in the post-DP thymocyte population observed in Ref. (21). We note that our approach is unrelated to that of Sinclair et al. both experimentally and mathematically (21). This rate, µ2, has been estimated to be at least an order of magnitude larger than any of the other death rates in the pre-DP, CD4 SP, or CD8 SP pools (see **Table 2**). In terms of selection rates, our analysis yields the following: pre-selection thymocytes (pre-DPs) have a 65.8% probability of dying by neglect in the cortex, and a 34.2% probability of becoming post-selection thymocytes (post-DPs). At the post-selection stage, post-DPs have a 91.7% probability of dying by negative selection (apoptosis) in the cortex, a 4.7% probability of becoming CD4 SPs, and a 3.6% probability of becoming CD8 SPs. In the medulla, CD4 SPs have an 8.6% probability of dying by negative selection (apoptosis), whereas CD8 SPs have a 32.1% probability of dying by negative selection. CD4 SPs have a 45.1% probability of exiting the thymus and reaching the periphery as mature thymocytes, whereas that probability for CD8 SPs is only 40.9%. Finally, the data supports some level of cellular proliferation in the medulla, with CD4 SPs having a 46.3% probability of proliferation and CD8 SPs a 27% probability.

Earlier work by Mehr and collaborators combined experimental and theoretical approaches to estimate thymic selection rates (13, 33), neglected death rates in the medulla, but considered potential feedback from mature T cells. In agreement with these authors, our results indicate that thymocyte death is highest at the post-DP stage. However, as death in the medulla had been neglected, these authors concluded that the CD4:CD8 ratio in SP thymocytes is determined by the differentiation rates. In this paper, we have made use of CD4 and CD8, or medullary, death rates, which allowed us to directly compare cortical (DP) to medullary (SP) death rates. Furthermore, our approach allowed us to conclude that medullary, or SP, death was due to negative selection, as it was rescued by Bim deficiency (26). Sinclair et al. also recently addressed the temporal dynamics of thymic selection using an unrelated approach (both experimentally and mathematically) (21). While their experimental approach did not allow them to distinguish death by negative selection from death by other mechanisms, their overall finding was consistent with ours, that thymocyte death is highest at the post-DP stage.

Further attempts to quantify thymic selection rates making use of mathematical models also include those of Faro et al. (12). The mathematical model developed by Faro and collaborators does not include time dynamics, but describes the relationship between the number of selecting ligands and the probability of selection of a given thymocyte. Thomas-Vaslin et al. (14) obtained estimates of thymic selection rates, using an experimental procedure that temporarily blocks thymic output and a mathematical model in which rates of transit from compartment to compartment depend on the number of cell divisions. Their model can capture the thymic "conveyor belt" (34, 35) scheme, but requires more differential equations and more parameters than equation (6). Despite the differences between their theoretical and experimental models and ours, similar estimates for thymic selection rates are found. For example, we estimate that 1.2 million post-DP become CD4 SP thymocytes per day and 0.5 million post-DP become CD8 SP thymocytes per day; their estimates are 0.9 and 0.2 million, respectively. Finally, we estimate that 2.6 million CD4 SP thymocytes per day and 0.6 million CD8 SP thymocytes per day exit the thymus. Their estimates are 2.4 and 0.5 million, respectively.

Our estimates of how many CD4 and CD8 SP thymocytes survive and exit the thymus reflect the skewed CD4:CD8 SP thymocyte ratio observed in C57BL/6 mice, which is approximately 4:1 (36). This ratio is similar to the reported CD4:CD8 ratio of recent thymic emigrants (37), and raises the question of what accounts for the CD4 bias. While we were able to determine death and differentiation rates for both CD4 and CD8 SP thymocytes (see **Table 2**), our experimental approach did not allow us to determine what fraction of the post-DP pool was MHC class I versus II restricted. Therefore, we could not address the issue of when and how the CD4:CD8 bias becomes established. The approach of Sinclair et al., which used MHC class I and class II deficiency, allowed them to address this question. Their data suggest that the skewed CD4:CD8 ratio reflects asymmetry in post-selection DP death rates, rather than more efficient positive selection of CD4 compared to CD8 thymocytes (21). Yet, the parameter estimation allows us to compare the following different CD4:CD8 ratios (see

Section 3.2): (i) the CD4:CD8 ratio of positive selection in the post-DP pool (differentiation from post-DP to either CD4 SP or CD8 SP) is given by ϕ4 ϕ8 ≈ 5 : 4, (ii) the CD4:CD8 ratio in the SP pool is given by *n* ∗ 4 *n* ∗ 8 ≈ 3 : 1, and (iii) the CD4:CD8 ratio of positive selection in the SP pool (differentiation from SP to peripheral early thymic emigrants) is given by ξ4*n* ∗ 4 ξ8*n* ∗ 8 ≈ 4 : 1. Our observations indicate that the CD4 bias is progressively established, as the thymocytes mature from the post-DP stage until the exit of the SP stage to migrate to the periphery.

Our mathematical analysis has also allowed us to estimate the stringency of thymic selection, defined by:

$$
\sigma = \frac{\xi\_4 \ n\_4^\* + \xi\_8 \ n\_8^\*}{\phi} = 8.96\%,
$$

that is, the ratio between the number of thymocytes per unit time that exit the thymus and the number of thymocytes per unit time that enter the pre-DP stage. The sensitivity analysis described in Section 3.3 allows us to provide a value of 1σ = 0.2%, where we have made use of the minimum–maximum interval ranges (see fourth column of **Table 2**). A different measure of stringency could be based on the probability of a cell surviving the maturation process. In our notation, this would correspond to the following:

$$(1 - p\_1) \times (1 - p\_2) \times (1 - p\_3) = 2.19\% \ . $$

We note that this measure of stringency is the probability of not dying in any of the three compartments considered in the model (pre-DP, post-DP, and SP). As discussed in Appendix A.1, and given that in the SP pool, thymocytes may proliferate, there is a need to consider this special case. Our estimates suggest that a population of 103 pre-DP thymocytes will yield 69 CD4 and 25 CD8 SP thymocytes that leave the medulla to get incorporated into the peripheral naive T cell pool (see details in Appendix A.1).

The sensitivity analysis (see Section 3.3) and the variability of the selection rates derived from it (see Section 3.4) give us the confidence to conclude, that our parameter estimation is robust. We are aware that the experimental data we have made use of [steady state thymocyte cell counts (26)] do not provide the exquisite time resolution described in Ref. (21). However, the supporting mathematical model described in Section 2.2, allows us to obtain the time evolution of the thymocyte populations, once the parameters have been estimated. In **Figure 5**, we plot the time evolution of the total number of cells in each compartment of the mathematical model: pre-DP, post-DP, CD4 SP, and CD8 SP thymocytes. We start with no cells at time zero, *ni*(*t* = 0) = 0 for *i* = 1, 2, 4, 8. Trajectories have been plotted for a period of 6 weeks and have been computed for every permutation of the parameter set presented in **Table 2**. The subset of parameters shared with the simple model (φ, ϕ1, µ1, µ2), were fixed at their mean values. Thus, 548 distinct parameter sets were generated. The system of equations (6) was solved using a fourth order Runge–Kutta method (Python source code).

The approaches introduced in this paper have shed some light on the probabilities and timescales that characterize cellular fate in the thymus after the DN stage. We plan to generalize the mathematical model introduced here, making use of experimental data for the strength of TCR binding in Nur77GFP mice (26), to investigate issues such as the death rate in the post-DP pool and the CD4:CD8 ratio. Our model assumes that all progenitors in a particular pool behave with identical kinetics, i.e., move through the various stages of selection at the same rate. Future model refinements will come from consideration of the heterogeneity of the pools, which are known to include cells that will become iNKT cells, regulatory T cells, and intraepithelial lymphocytes (2). It is also possible that progenitors of the same general class move through the selection process with different kinetics (34). The models introduced here can serve as a first step to study human thymic selection, although comprehensive data on human thymic subsets, their sizes, and residence times are not yet available. It would be of great interest to apply the model to data on thymic subsets and cellularity in children, keeping in mind that residence times of human subsets may differ from murine ones (38). Finally, we note that we have not mentioned the relevance of cytokines, such as IL-7, during thymic development. Some differences have already been described for the role of IL-7R in human versus mouse T cell development (38, 39). We hope in the near future to combine mechanistic mathematical models of IL-7 and IL-7R (40) with the T cell development model introduced here to address these issues.

### **ACKNOWLEDGMENTS**

We thank Ed Palmer, Robin Callard, Andy Yates, and Ben Seddon for helpful discussions, and Andy Yates and Ben Seddon for allowing us to make use of their mathematical estimates for µ4 and µ8. Grant Lythe and Carmen Molina-París thank Bill Hlavacek, Ruy Ribeiro, and Alan Perelson (Los Alamos National Laboratory) for their hospitality. This work was supported by a BBSRC Research Development Fellowship BB/G023395/1 (to Carmen Molina-París), an FP7 IRSES INDOEUROPEAN-MATHDS Network PIRSES-GA-2012-317893 (to Grant Lythe and Carmen Molina-París), and National Institutes of Health Grants R37 AI39560 (to Kristin A. Hogquist) and F32 AI100346 (to Gretta L. Stritesky). Maria Sawicka, Joseph Reynolds, Niloufar Abourashchi, Grant Lythe, and Carmen Molina-París acknowledge the hospitality of 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.

#### **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: 11 September 2013; accepted: 15 January 2014; published online: 14 February 2014.*

*Citation: Sawicka M, Stritesky GL, Reynolds J, Abourashchi N, Lythe G, Molina-París C and Hogquist KA (2014) From pre-DP, post-DP, SP4, and SP8 thymocyte cell counts to a dynamical model of cortical and medullary selection. Front. Immunol. 5:19. doi: 10.3389/fimmu.2014.00019*

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

*Copyright © 2014 Sawicka, Stritesky, Reynolds, Abourashchi, Lythe, Molina-París and Hogquist. 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.*

# **APPENDIX**

#### **A.1. STRINGENCY OF THYMIC SELECTION: A STOCHASTIC MODEL**

In this section, we present the details that allow us to compute the stringency of thymic selection for the mathematical model considered in Section 2.2.

Let us assume that at time *t* = 0, there exists a single T cell in a given compartment. In our case, the compartment can be the pre-DP, post-DP, or SP (either CD4 SP or CD8 SP) stages. The cell, at any time, may die (with rate, µ), divide (with rate, λ), and produce two daughter cells, or leave the compartment (with rate, ξ ) to enter a different compartment. The waiting times for each event are assumed to be exponentially distributed, and daughter cells are assumed to behave identically to the initial single cell. We introduce the bivariate Markov process {*X*(*t*),*Y* (*t*)}*t*≥0, where *X*(*t*) is the number of cells in the compartment at time *t*, and *Y* (*t*) is the number of cells which have left the compartment (to enter a different compartment). Our aim is to calculate the expected number of cells (and variance) that leave the compartment.

State probabilities for the Markov process are defined as follows (41)

$$p\_{(\mathbf{x},\mathbf{y})}(t) = \text{Prob}\{X(t) = \mathbf{x}, Y(t) = \mathbf{y} | X(0) = 1, Y(0) = 0\},\tag{\text{A1}}$$

and transition probabilities for this process are defined as follows (41)

$$\begin{aligned} \mathcal{P}\_{(w,z),(x,\gamma)}(\Delta t) &= \text{Prob}\{X(t+\Delta t) = w, \\ Y(t+\Delta t) &= z \vert X(t) = x, Y(t) = \gamma\} \\ &\quad \begin{cases} \lambda x \Delta t + o(\Delta t), & \text{if } w = x+1, z = \gamma \text{ ,} \\ \mu x \Delta t + o(\Delta t), & \text{if } w = x-1, z = \gamma \text{ ,} \\ \xi x \Delta t + o(\Delta t), & \text{if } w = x-1, z = \gamma+1 \text{ ,} \\ 1 - (\lambda + \mu + \xi) x \Delta t + o(\Delta t), & \text{if } w = x, z = \gamma \text{ ,} \\ o(\Delta t), & \text{if } w, z \text{ otherwise } \end{cases} \end{aligned} \tag{A2}$$

The Kolmogorov (or master) equation for this process is given by (41)

$$\begin{split} \frac{\mathrm{d}p\_{(\mathbf{x},\mathbf{y})}(t)}{\mathrm{d}t} &= \lambda \left(\mathbf{x} - \mathbf{1}\right) p\_{(\mathbf{x}-1,\mathbf{y})}(t) + \mu \left(\mathbf{x} + \mathbf{1}\right) p\_{(\mathbf{x}+1,\mathbf{y})}(t) \\ &+ \xi \left(\mathbf{x} + \mathbf{1}\right) p\_{(\mathbf{x}+1,\mathbf{y}-1)}(t) - \left(\lambda + \mu + \xi\right) \ge p\_{(\mathbf{x},\mathbf{y})}(t) \,. \end{split} \tag{A3}$$

Let *mX*(*t*) be the expected number of cells in the compartment under consideration, and *mY*(*t*) be the expected number of cells which have left the compartment. Similarly, let *mXX*(*t*) be the expectation of the random variable *X*(*t*) 2 , *mXY*(*t*) be the expectation of the random variable *X*(*t*)*Y* (*t*), and *mYY*(*t*) be the expectation of the random variable *Y* (*t*) 2 . Then, making use of the probability generating function technique (41), we derive the time evolution for the first two moments of the system:

$$\frac{\mathrm{d}m\_{\mathrm{X}}(t)}{\mathrm{d}t} = (\lambda - \mu - \xi)m\chi(t)\,,\tag{A4}$$

$$\frac{\mathrm{d}m\_{\mathrm{Y}}(t)}{\mathrm{d}t} = \xi \, m\_{\mathrm{X}}(t) \,, \tag{A5}$$

$$\frac{\mathrm{d}m\_{\mathrm{XX}}(t)}{\mathrm{d}t} = 2(\lambda - \mu - \xi)m\_{\mathrm{XX}}(t) + (\lambda - \mu - \xi)m\_{\mathrm{X}}(t), \quad \text{(A6)}$$

$$\mathrm{d}m\_{\mathrm{XY}}(t) \underset{\longrightarrow}{\longrightarrow} \quad \text{ $\tau$ }\,\mathrm{(\mu\nu\dots\tau)}\,\mathrm{; } \tau\mathrm{\&}\,\mathrm{\,\mu\nu\dots}(t) \quad \text{\(\mu\nu\dots\tau\)}\,\mathrm{; } \tau\mathrm{\&}\,\mathrm{/}$$

$$\begin{cases} \frac{d\mathbf{u}\_{X\Gamma}(t)}{dt} = (\lambda - \mu - \xi)m\_{XY}(t) + \xi \left[m\_{X\Gamma}(t) - m\_X(t)\right], & (\text{A7})\\ \mathbf{d}m\_{Y\Gamma}(t) \end{cases}$$

$$\frac{\mathrm{d}m\_{YY}(t)}{\mathrm{d}t} = \xi [2m\_{XY}(t) + m\_X(t)].\tag{A8}$$

Given that we start with a single cell, the expected number of cells at time *t* is given by

$$m\_{\mathbf{X}}(t) = \mathbf{e}^{(\lambda - \mu - \xi)t} \,. \tag{A9}$$

Under the restriction λ < µ + ξ , the expected number of cells tends to zero as *t*→ +∞. This implies that all cells from the single T cell progenitor either die or leave the compartment for sufficiently large times. The expected number of cells which leave the compartment is given by

$$m\_Y(t) = \frac{\xi}{\mu + \xi - \lambda} \left[ 1 - e^{(\lambda - \mu - \xi)t} \right]. \tag{A10}$$

As *t*→ +∞, the expected number of cells to leave the compartment can be shown to be

$$\lim\_{t \to +\infty} m\_Y(t) = \frac{\xi}{\mu + \xi - \lambda} \,. \tag{A11}$$

We now solve the remaining ODEs equations (A6–A8), to find

$$\begin{split} m\_{YY}(t) &= \frac{2\lambda\xi^2}{\left(\lambda-\mu-\xi\right)^3} \left[ \frac{1}{2} \mathbf{e}^{2(\lambda-\mu-\xi)t} - \mathbf{e}^{(\lambda-\mu-\xi)t} \right] \\ &- \frac{4\lambda\xi^2}{\left(\lambda-\mu-\xi\right)^2} \left(t-\frac{1}{\lambda-\mu-\xi}\right) \mathbf{e}^{(\lambda-\mu-\xi)t} \\ &+ \frac{\xi}{\lambda-\mu-\xi} \mathbf{e}^{(\lambda-\mu-\xi)t} - \frac{2\lambda\xi^2}{\left(\lambda-\mu-\xi\right)^3} \\ &- \frac{\xi}{\lambda-\mu-\xi} .\end{split} \tag{A12}$$

It, therefore, follows that the random variable *Y* (*t*), which represents the number of cells leaving the compartment under consideration, has the following variance (in the limit *t*→ +∞)

$$
\sigma\_Y^2 = \lim\_{t \to \infty} \left[ m\_{YY}(t) - m\_Y(t)^2 \right] = \frac{2\lambda \xi^2}{\left(\mu + \xi - \lambda\right)^3} + \frac{\xi}{\mu + \xi - \lambda}
$$

$$
$$

# **A.2. STRINGENCY OF THYMIC SELECTION IN THE FOUR COMPARTMENT MODEL**

The previous example can easily (but laboriously) be extended to the mathematical model introduced in Section 2.2.

We may evaluate the expected number of, for example, CD4+ T cells produced by a single pre-DP progenitor (or more generally *N* pre-DP progenitors). We only present the time evolution of the moment generating function. The deterministic equations describing the mean and variance of the numbers of cells in each compartment (pre-DP, post-DP, SP CD4, and SP CD8) are left for the reader to derive. When counting the number of CD4+ T cells leaving the thymus, the moment generating function satisfies the following partial differential equation

$$\begin{split} \frac{\partial M}{\partial t} &= \mu\_1 (\mathbf{e}^{-\theta\_1} - 1) \frac{\partial M}{\partial \theta\_1} + \varphi\_1 (\mathbf{e}^{-\theta\_1} \mathbf{e}^{\theta\_2} - 1) \frac{\partial M}{\partial \theta\_1} \\ &+ (\mu\_2 + \varphi\_8) (\mathbf{e}^{-\theta\_2} - 1) \frac{\partial M}{\partial \theta\_2} + \varphi\_4 (\mathbf{e}^{-\theta\_2} \mathbf{e}^{\theta\_4} - 1) \frac{\partial M}{\partial \theta\_2} \\ &+ \mu\_4 (\mathbf{e}^{-\theta\_4} - 1) \frac{\partial M}{\partial \theta\_4} + \lambda\_4 (\mathbf{e}^{\theta\_4} - 1) \frac{\partial M}{\partial \theta\_4} \\ &+ \xi\_4 (\mathbf{e}^{-\theta\_4} \mathbf{e}^{\theta\_8} - 1) \frac{\partial M}{\partial \theta\_4} . \end{split} \tag{A14}$$

The symmetry of the mathematical model implies that an equivalent equation for the number of CD8+ T cells leaving the thymus can be obtained by interchanging the indexes 4 and 8. For our derived parameter set, the previous equation allows us to conclude that the expected number of CD4+ T cells, a single thymocyte in the pre-DP compartment produces, is 0.069 (standard deviation 0.96), whereas the expected number of CD8+ T cells which leave the thymus is 0.025 (standard deviation 0.59).

To put this into perspective, a population of 103 pre-DP thymocytes is expected to produce 69 CD4+ T cells which leave the thymus (standard deviation 66), and 25 CD8+ T cells (standard deviation 41). More generally, a population of *N* pre-DP thymocytes is expected to produce 0.069*N* CD4+ T cells and 0.025*N* CD8+ T cells.

# Asymmetry of cell division in CFSE-based lymphocyte proliferation analysis

#### **Gennady Bocharov 1\*,Tatyana Luzyanina2 , Jovana Cupovic 3 and Burkhard Ludewig3**

1 Institute of Numerical Mathematics, Russian Academy of Sciences, Moscow, Russia

2 Institute of Mathematical Problems in Biology, Russian Academy of Sciences, Pushchino, Russia

3 Institute of Immunobiology, Kantonal Hospital St. Gallen, St. Gallen, Switzerland

#### **Edited by:**

Rob J. De Boer, Utrecht University, Netherlands

#### **Reviewed by:**

Vitaly V. Ganusov, University of Tennessee, USA Ollivier Hyrien, University of Rochester, USA

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

Gennady Bocharov, Institute of Numerical Mathematics, Russian Academy of Sciences, Gubkina Street 8, Moscow 119333, Russia e-mail<EMAIL_ADDRESS>
Flow cytometry-based analysis of lymphocyte division using carboxyfluorescein succinimidyl ester (CFSE) dye dilution permits acquisition of data describing cellular proliferation and differentiation. For example, CFSE histogram data enable quantitative insight into cellular turnover rates by applying mathematical models and parameter estimation techniques. Several mathematical models have been developed using different types of deterministic or stochastic approaches. However, analysis of CFSE proliferation assays is based on the premise that the label is halved in the two daughter cells. Importantly, asymmetry of protein distribution in lymphocyte division is a basic biological feature of cell division with the degree of the asymmetry depending on various factors. Here, we review the recent literature on asymmetric lymphocyte division and CFSE-based lymphocyte proliferation analysis. We suggest that division- and label-structured mathematical models describing CFSEbased cell proliferation should take into account asymmetry and time-lag in cell proliferation. Utilization of improved modeling algorithms will permit straightforward quantification of essential parameters describing the performance of activated lymphocytes.

**Keywords:T cells, CFSE assay, asymmetric division, mathematical modeling**

# **INTRODUCTION**

The ability of the immune system to protect the host organism against live-threatening infections and tumors directly depends on the reactivity of lymphocytes to antigenic stimulation, with a key role of clonal T cell responses (1). The perception of infections as a race between the invading pathogen and immunity suggests that it is the knowledge of the proliferation and death rates of T cells which provides a quantitative basis for assessing the quality of the host immunity (2). For almost 20 years, flow cytometry-based analysis of intracellular fluorescent dye distribution has been used to assess the proliferative performance and differentiation patterns of lymphocytes (3–5). Since the prototype dye for this analysis is CFSE, the assay is commonly referred to as CFSE dilution assay or – more simply – CFSE assay. A quantitative characterization of T cell turnover which can be elaborated from time series of CFSE histograms ranges from"static"measures such as precursor cell frequency or mean generation number, to "dynamic"parameters characterizing the cell cycle progression and apoptosis rates (6). However, estimation of turnover parameters requires formulation of mathematical models of cell growth which can take various forms and differ in their complexity depending on the parameters of interest and the richness of the available data [comprehensively reviewed by De Boer and Perelson (7)]. Importantly, current approaches to the analysis of CFSE proliferation data are based on the assumption that cell division is symmetric, i.e., the fluorescent label is halved in the two daughter cells (3, 5, 7–9). However, a random and uneven partition of mass between the sister cells is considered as an axiom in cell biology since many years (10). Although the detailed knowledge of the intracellular

reactions which affect the turnover and intracellular heterogeneity of CFSE labeled proteins is currently limited (11), it is broadly accepted that CFSE binds indiscriminately to intracellular proteins and the fluorescence intensity of any single cell is roughly proportional to the total number of CFSE molecules bound to proteins within that cell (12). The latter study proposed a method for the analysis of CFSE-labeling experiments which also considered the possibility of an unequal division of CFSE molecules between the daughter cells.

The inequality of the mass (protein) distribution to the daughter cells directly suggests that CFSE labeled proteins are unequally partitioned between daughter cells. Indeed, recent studies describing T cell activation showed that asymmetric cell division can be an inherent part of T cell growth and differentiation (13–16). However, direct experimental evidence for asymmetric partition of CFSE between daughter cells is still missing. Nevertheless, the existing deterministic mathematical frameworks should be amended to facilitate a quantitative analysis of CFSE-based lymphocyte proliferation when asymmetry of cell division associated with unequal partition of CFSE labeled proteins between the two daughter cell results in a poor resolution of divisional clusters in CFSE histograms. Here, we briefly summarize recent findings describing asymmetric lymphocyte division and progress in the analysis of CFSE-based lymphocyte activation. Moreover, cell proliferation is not an instantaneous process and it takes a finite time for a cell to progress from the G1-phase of the cell cycle to the completion of the M-phase. The duration of the continuous progression is called a time-lag and in general, needs to be explicitly parameterized in the model equations. Finally, we suggest that

mathematical models describing CFSE-based lymphocyte proliferation should consider both asymmetry in division and time-lag in proliferation.

# **ASYMMETRIC LYMPHOCYTE DIVISION**

Symmetric or asymmetric cell divisions refer to the mode of cell division which results in two phenotypically identical- or differentdaughter cells, respectively. The phenotypic features could be the cell size, cell surface receptors, intracellular components such as proteins (including those labeled with CFSE), transcription factors, or messenger RNA (17). Hence, these phenotypic differences provide the basis for the functional differences in the daughter cells, i.e., their cell fates.

Following encounter with their antigen displayed in the context of major histocompatibility complex molecules, naïve T lymphocytes go through well-orchestrated series of divisions generating different populations of cells that fulfill immediate effector functions or generate long-lived immunological memory. Two basic models explain the generation of such functionally distinct T cell phenotypes. According to the "one naïve cell – one fate" model, naïve lymphocytes are instructed to generate either effector or memory progeny (18). In this model, instruction of T cells, for example, is achieved through interaction with professional APCs (19). Hence, to preserve the instructing signal(s) received during activation and to maintain equality of the cells throughout division, T cells should divide in a symmetric fashion. The alternative model proposes asymmetric cell division as the mechanism that allows naïve T cells to give rise to two different daughter cells. These are referred to as proximal or distal daughter cell depending on their proximity to the APC. Such asymmetric T cell division represents a process that allows single cells to give rise to two, phenotypically and functionally different daughter cells, and thereby permits diversification of cell populations. In other words, one of the daughter cells inherits the potential to differentiate into full effector cell (proximal daughter), while the second daughter maintains the stemness of the mother cell. This principle feature of asymmetric cell division has also been described in developmental studies examining neurogenesis (20). Likewise, adaptation of adult tissues to changing environmental conditions such as the content of the gut requires rapid adaptation of one cell fraction while other cells maintain their high proliferative potential (21).

The processes involved in activation and differentiation of T cells, for example during infection have to swiftly generate cells with direct effector function to efficiently restrict viral replication (1). At the same time, some T cells should retain their ability to proliferate in order to prevent exhaustion of certain T cell subsets (22) and to facilitate generation of long-lived memory T cells (23). Indeed, Chang et al. (13) demonstrated that division of CD8+ T cells specific for a viral peptide leads to the generation of daughter cells with different characteristics. CFSE-based assays revealed that asymmetry is established already during the first round of division and is dependent on the presence of the cognate antigen (13). Assessment of the protein content in the daughter cells generated during the first cell division showed that asymmetry established during mitosis is preserved throughout cytokinesis. Moreover, proximal and distal daughter cells exhibit different protein expression profiles and functional properties with proximal daughter cells exhibiting higher immediate protective capacity (13). The finding that proximal daughter cells exhibit higher CD8 co-receptor and LFA-1 expression facilitating formation of more frequent and longer lasting interactions with antigen presenting APCs (14) further emphasized that asymmetric division critically determines both T cell phenotype and function.

Asymmetric cell division is not only an important feature of CD8+ T cell activation (13, 14), but also occurs during the activation and differentiation of CD4+ T cells (13, 24) and B cells (25, 26). While naïve CD8+ T cells require only one or only few encounters with APCs to proliferate and differentiate into effector cells, naïve CD4+ T cells depend on multiple encounters in order to differentiate and to exhibit specialized effector functions (27). Hence, it is likely that CD4+ T cells acquire their distinct phenotypes, e.g., Th1, Th2, or Th17, through multiple sequential asymmetric cell divisions. However, recent studies suggest that asymmetric cell division cannot be considered as the only mechanism that leads to the profound heterogeneity of T cell lineages (16). Thus, more research is required to resolve the contribution of sequential asymmetric T cell division to the generation of diverse T cell phenotypes. We suggest that a combination of CFSE-based T cell proliferation analysis with mathematical modeling may help – at least in part – to clarify this issue.

# **CURRENT MATHEMATICAL MODELS FOR CFSE-BASED LYMPHOCYTE PROLIFERATION ANALYSIS**

Several mathematical models have been established for the analysis of CFSE-based proliferation assays (7, 9, 12, 28–32). The existing modeling frameworks can be subdivided on the basis of the major requirements for CFSE histogram data processing into two main categories (**Table 1**). The first group requires a decomposition of the CFSE histograms characterizing the distribution of cells with respect to the fluorescent dye into the distinct generations of cells. The procedure is based on fitting the CFSE histogram with a series of log-normal Gaussian distributions differing in their means and standard deviation and is implemented in commercially available standard software packages. Importantly, the assignment of distinct cell generations to CFSE clusters has remained an empiric process which depends heavily on initial labeling homogeneity, label degradation, cellular auto-fluorescence, and other factors including experimental skills of the researcher (33). As long as the division is symmetric (or almost symmetric) (**Figure 1A**), these factors can be tuned in a proper way to enable resolution of successive generations as distinct CFSE clusters (**Figure 1B**). Under these conditions a range of existing mathematical models can be tuned to estimate the turnover parameters of the stimulated lymphocyte population. The key features of the corresponding families of the models are outlined in **Table 1**, rows one to three. These models describe the population dynamics of cells which differ in the number of completed divisions and ignore the heterogeneity of the cells within a generation with respect to the CFSE content. The immunologically relevant issues that were addressed with the models of this group include regulatory effects of IL-2 on the T cell responses (34, 35), regulation of hematopoietic stem cells cycling (36), and kinetics of mouse erythroid progenitor cell differentiation (37).

**Table 1 | Major features of mathematical models describing CFSE-based proliferation assays.**

1The following notations are used: "A–B state model" refers to the model of cell cycle in (55) in which the intermitotic period is composed of an A-state (major part of G1-phase) and a B-phase (conventional S, G2, and M phases); "G0 model" refers to the view of the cell cycle with two states (47), i.e., resting- (G0) and cyclingstates (G1, S, G2, M). Conceptually, it is equivalent to the A–B state model. "Random birth-death" model refers to a discrete compartmental (generation structured) model of cycling cells (56). "Generation structured" refers to the mathematical model in which the cell population is decomposed into cohorts of cells which differ with respect to the number of completed division cycles; "CFSE-structured" represents the mathematical description of cell population in which the distribution (heterogeneity) of cells with respect to the fluorescence intensity is followed by considering the cell distribution function.

2DDE, delay differential equations; hPDE, hyperbolic partial differential equations; ODE, ordinary differential equations; IE, integral equations.

The second group of models which refer directly to the CFSE histograms seems to be more appropriate when the generational structure of the labeled population cannot be easily resolved (**Table 1**, rows four to five). The initially proposed model describes the evolution of the labeled cells distribution with respect to the CFSE level (38). Although this and similar models proved to be functional in estimating the proliferation- and death rates as functions of the structure variable directly from the histogram data (38, 39), the problem of translating the estimated functions into biologically meaningful parameters still requires the knowledge of the division structure of the lymphocyte population. A major breakthrough in the improvement of the distributed parameter models for the dynamics of heterogeneous CFSE labeled cell populations were recently proposed division- and label-structured mathematical models (32, 40). The major potential of this framework as an analytical tool is based upon the following features: (i) no need for CFSE histogram decomposition, (ii) characterization of cell growth in terms of generation dependent division- and death rates, (iii) an explicit form of the dependence of solution on the turnover parameters.

Another class of recently developed mathematical models which allow a direct fitting of the CFSE histograms is based on branching processes (12, 41). The approach allows for probabilistic characterization of cell activation, proliferation, and death from the CFSE dilution data and does not require the assumption about equality of CFSE division between the two daughter cells.

The first and the second group of models rely on the premise of symmetric cell division. However, tracing proliferation of other cell types such as cancer cells has been reported to be difficult (42) due to poorly resolved peaks of the different cell generations. Since cytokinesis is not perfect, it was suggested that the two daughter cells are unlikely to inherit exactly half of the CFSE fluorescence dye of the mother cell. An increase in the degree of the asymmetry of mass partition between daughter cells and hence disparate distribution of fluorescently labeled proteins should result in a poorer resolution of generational clusters as shown in **Figures 1C,D** for lower asymmetry and in **Figures 1E,F** for higher asymmetry. This in turn will lead to the generations overlap in CFSE histograms thus posing a limit to experimentalists' ability to resolve the individual generations using conventional decomposition methods.

# **MODELING ASYMMETRIC DIVISION OF CFSE LABELED T CELLS**

We have been recently dealing with the analysis of the proliferative performance of monoclonal CD8+ T cells recognizing an H2-Kb-binding epitope derived from the S protein of the mouse hepatitis virus (MHV). Clearance of MHV during acute infection is achieved through the combined action of type I interferons (43) and CD8+ T cells (44). Moreover, CD8+ T cells essentially contribute to control of the virus during persistent infection, for example in the central nervous system (45). We have initiated a project on the generation of avidity-tuned, antigen-specific T cells

**FIGURE 1 | Impact of asymmetry inT cell division impinges on fluorescent protein partition between daughter cells**. **(A,B)** Symmetric cell division with equal distribution of the fluorescent dye between daughter cells **(A)** and modeled time course analysis of T cell proliferation as determined by flow cytometry [**(B)**, solid black lines]. Dashed red lines in **(B)** indicate the evolution of CFSE intensity of the cohorts (generations) of cell which differ in the number of completed divisions with the assumption of symmetric division. **(C,D)** Asymmetric cell division with "low" asymmetry **(C)** and modeled flow cytometric time course analysis of CFSE dilution [**(D)**, solid black lines] that corresponds to an asymmetry 46/54% [**(D)**, dashed red lines describe the CFSE distributions for cell cohorts differing in terms of the completed divisions]. **(E,F)** T cells dividing with "high" asymmetry **(E)** and corresponding model-generated flow cytometric CFSE dilution patterns [**(F)**, solid black lines] with asymmetry values of 42/58% describing the behavior of the T cells in this setting [**(F)**, dashed red lines describe the cell

cohorts corresponding to different generations]. **(G)** Schematic representation of the structure of a mathematical approach which considers the division- and CFSE label-heterogeneity of proliferating cells as well as asymmetry and time of cell division. Some cells from the cohort of cells which completed "i" divisions are activated (αi characterizes the activation rate) and progress through the cell cycle (τi stands for the duration of the progression through S-G2-M phases), resulting to the generation of daughter cells which differ with respect to their CFSE content. Asymmetric mitosis refers to cell division which results into appearance of two phenotypically different daughter cells with a smaller and larger cell mass, respectively. These cells are characterized by an unequal amount of CFSE labeled proteins (m1 and m2 = 1−m1, describe the fractions of CFSE from the mother cell inherited by the two daughter cells). The natural decay of the CFSE fluorescence intensity is taken into account (kx – stands for an exponential decay of CFSE loss).

for adoptive transfer as an option to augment antiviral immune responses during chronic infection. To this end, MHV-specific T cell receptors (TCRs) were cloned and tested in retrogenic systems (46). *In vitro* re-stimulation of the CFSE labeled monoclonal CD8+ T cells showed that CFSE dilution was characterized by broadly varying patterns from highly distinct peaks to poorly resolved generational clusters. We propose that an explicit consideration of the asymmetry in protein partition between the daughter cells facilitates a consistent mathematical description of CFSE histogram time series data (**Figure 1G**). The appropriate mathematical framework should describe the population of CFSE labeled T cells by the distribution of cells with respect to CFSE amount (unit of intensity, UI). The subpopulations differ in terms of completed rounds of division and are further distinguished in resting and proliferating states, with the respective notation and *i* standing for the generation (number of completed divisions), *t* – for time and *x* – for CFSE amount per cell. A conceptual scheme of the modeling approach is shown in **Figure 1G** suggesting that such a model can be naturally formulated as an extension of a generation- and division-structured population balance model with the cell cycle represented according to the G0 model (47) and the division asymmetry explicitly taken into account.

Under conditions of symmetric CD8+ T cell division with the difference of protein partition between the sister cells being equal to zero (i.e., every daughter cell inherits half of the fluorescently labeled proteins of the mother cell), the model should predict clearly distinct generations (**Figure 1B**, dashed red lines). If the division is "weakly" asymmetric, i.e., the protein partition between the sister cells is different, the width of the CFSE distribution of the successive generations should become broader (**Figure 1D**, dashed red lines). Further increase in the degree of the asymmetry would result in a substantial overlap of the distinct cell generations (**Figure 1F**, dashed red lines). Obviously, this type of behavior of T cells – and other cells such as tumor cells needs to be regarded as a cause of a poor resolution of the generations in CFSE histograms (**Figures 1D,F**, solid black lines) thus creating an obstacle on the application of standard CFSE analysis tools.

The fitting of mathematical models for asymmetric cell division as conceptualized in **Figure 1G** to the time series data provides a tool for the estimation of the cell physiology parameters such as: (i) the generation-specific activation and death rates (α*i* , β*i*); (ii) the duration of the division cycle characterized by the time-lag (τ*i*); (iii) the division asymmetry factors (*m*1 + *m*2 = 1), specifying the fraction of proteins which is inherited by the first and the second daughter cells, respectively; and (iv) the natural decay of the CFSE fluorescence intensity of the labeled cells (parameterized as

#### **REFERENCES**

by flow cytometry. *J Immunol Methods* (1994) **171**:131–7. doi:10. 1016/0022-1759(94)90236-4

*kx*). Taken together, asymmetric cell division improves assessment of T cell performance parameters from CFSE-based proliferation assays, even under conditions of poorly separated peaks.

#### **CONCLUDING REMARKS**

It is considered that the regulation of cell expansion and differentiation can occur by modulating the degree of asymmetry of cell divisions (17). It has been clearly shown that T lymphocyte division in response to pathogen exhibits unequal partitioning of proteins that mediate signaling and cell fate determination (13). Hence, asymmetric T lymphocyte division provides an additional mechanism for generating functionally heterogeneous populations of CD8+ T cells both in primary and memory adaptive immune responses (48). Since a precise mechanistic link between the quantitative differences in partitioning of specific proteins between daughter cells and the developmental path of antigenspecific T cells remains to be established (49), mathematical modeling is now a key "instrument" for understanding the regulation of individual cell fates (15, 16, 50).

The addition of asymmetric T cell division to the analysis of CFSE-based proliferation data fills important gaps as it: (i) allows one to estimate the proliferation parameters for asymmetrically dividing cells directly from CFSE histograms with poorly resolved generations peaks and (ii) introduces a quantitative parameter which characterizes the difference in the partition of the fluorescent proteins between daughter cells and can be directly estimated from the same CFSE dilution data. A further question in CFSE analyses open for examination is the interplay between experimental variability, biological variability, and model parsimony. We expect that new mathematical tools for the analysis of a fundamental property of cell division, i.e., the phenotypic identity or differences of the daughter cells known as asymmetry, will be developed and introduced into daily experimental work. Thereby, a better understanding of the diversity and mechanisms underlying activation and homeostasis of T cell responses will be achieved.

#### **ACKNOWLEDGMENTS**

This work has been supported by the Swiss National Science Foundation (130823 and 141918 to Burkhard Ludewig), the Vontobel Foundation (to Burkhard Ludewig), the Russian Foundation of Basic Research (11-01-00117a to Gennady Bocharov and Tatyana Luzyanina), and the Program of the Russian Academy of Sciences "Basic research for Medicine" (to Gennady Bocharov and Tatyana Luzyanina). We thank the reviewers and the editor for their insightful comments and the thorough work on our manuscript.

and in vivo using CFSE-like fluorescent dyes.*J Immunol Methods* (2012) **379**:1–14. doi:10.1016/j.jim. 2012.02.012

compartmental model for computing cell numbers in CFSEbased lymphocyte proliferation assays. *Math Biosci Eng* (2012) **9**: 699–736. doi:10.3934/mbe.2012.9. 699

**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: 03 June 2013; accepted: 19 August 2013; published online: 02 September 2013.*

*Citation: Bocharov G, Luzyanina T, Cupovic J and Ludewig B (2013) Asymmetry of cell division in CFSEbased lymphocyte proliferation analysis. Front. Immunol. 4:264. doi: 10.3389/fimmu.2013.00264*

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

*Copyright © 2013 Bocharov, Luzyanina, Cupovic and Ludewig . 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.*

# Dynamical and mechanistic reconstructive approaches of T lymphocyte dynamics: using visual modeling languages to bridge the gap between immunologists, theoreticians, and programmers

# **Véronique Thomas-Vaslin1,2\*, Adrien Six 1,2, Jean-Gabriel Ganascia3,4 and Hugues Bersini 5**

1 UPMC Univ Paris 06, UMR7211, Immunology, Immunopathology, Immunotherapy (I3),Integrative Immunology, Paris, France

4 CNRS, UMR 7606, ACASA-LIP6, CNRS, Paris, France

5 Université Libre de Bruxelles, IRIDIA-Code, Bruxelles, Belgium

#### **Edited by:**

Rob J. de Boer, Utrecht University, Netherlands

#### **Reviewed by:**

Christopher E. Rudd, University of Cambridge, UK Yoram Louzoun, Bar-Ilan University, Israel Hillel Kugler, Microsoft Research, UK **\*Correspondence:**

Véronique Thomas-Vaslin, UMR 7211, Immunology, Immunopathology, Immunotherapy, Integrative Immunology, UPMC-CNRS, 83 Boulevard de l'Hôpital, Paris 75013, France e-mail: veronique.thomas-vaslin@ upmc.fr

Dynamic modeling of lymphocyte behavior has primarily been based on populations based differential equations or on cellular agents moving in space and interacting each other. The final steps of this modeling effort are expressed in a code written in a programing language. On account of the complete lack of standardization of the different steps to proceed, we have to deplore poor communication and sharing between experimentalists, theoreticians and programmers. The adoption of diagrammatic visual computer language should however greatly help the immunologists to better communicate, to more easily identify the models similarities and facilitate the reuse and extension of existing software models. Since immunologists often conceptualize the dynamical evolution of immune systems in terms of "state-transitions" of biological objects, we promote the use of unified modeling language (UML) state-transition diagram. To demonstrate the feasibility of this approach, we present a UML refactoring of two published models on thymocyte differentiation. Originally built with different modeling strategies, a mathematical ordinary differential equation-based model and a cellular automata model, the two models are now in the same visual formalism and can be compared.

**Keywords: state-transition diagram, computer modeling, cell dynamics, agent-based model, complex system**

The perspective is to encourage immunologists involved into mathematical modeling or software productions, to adopt a visual graphical language, here mainly the unified modeling language (UML) "state-transition" diagram to ease the communication, the reuse and the extension of their models.

# **COMPLEXITY OF THE IMMUNE SYSTEM**

The immune system is a complex biological adaptive, highly diversified, robust and resilient system, characterized by complexity at different levels. Lymphocytes are the central actors of the immune system, in the middle of a multi-scale biological organization, "from molecule to organism". Multi-scale modeling remains a challenge, as for other biological systems (1). Despite recent systems biology initiatives to understand and model the immune system (2), we are still far from having the appropriate tools to understand its dynamics and to easily communicate among various researchers who observe this system at different levels of granularity and attempt through software modeling to answer different questions. Several complementary experimental methods and models have been used to explore lymphocyte dynamics and turnover (3, 4) and to model it in health, aging and diseases (5).

# **DRAWBACKS OF CURRENT DYNAMICS LYMPHOCYTE MODELING AND EVOLUTION**

System dynamics models deal with time, formalized with two distinct concepts as "discrete time," by a succession of time points and intervals, or as "continuous time" (6). Models of lymphocyte population dynamics and turnover (3) have primarily been based on mechanistic reconstruction with continuous time models. The fluxes of cell populations are then described by differential equations. These mathematical models describe for example the thymocyte differentiation and selection (7), until the thymic export (8–10), the homeostasis of CD4/CD8 T cells (11, 12), the CD8 immune response (13, 14), or the Bromodeoxyuridine or deuterium labeling (15) to account for turnover. Up to now, simulations and validation of some of these models reveal interesting T cell dynamics properties: how the system grows, self-maintains as well as the effects of perturbations, i.e. how the system reacts to antigens, collapses and reorganizes. However, integrating the heterogeneity of cell populations, phenotypes, lineages, cell location and interactions, cell differentiation across generations (16) in the different biological, and time scales, is problematic in such a mathematical form, which make these models particularly difficult to handle.

The evolution of homogeneous mathematical model of cell populations (17) toward "spatialized," discrete, and heterogeneous software models (18) has allowed the reproduction and observation of more detailed and thus complex behaviors. For example, this made possible to model lymphocyte dynamics from thymic selection (19, 20) up to quantitative modeling of immune responses, as extensively reviewed (21) with development of agentbased and automata models (22). However, both population-based mathematical model (a top-down approach) and discrete cellbased model (a bottom-up approach) and the various platforms developed have limitations (23). Conversely, the Statecharts language (24) and the visual reactive tools (25) such as biocharts (26) and reactive animation applied to various systems (27) developed by Harel et al. are a powerful way to simulate complex dynamical biological behavior with more didactic representation than equations. Such models have revealed emergent properties during thymic differentiation (28) pancreatic islet organogenesis but also the immune response in lymph node (29).

# **LACK OF INTEROPERABILITY, UNDER-USE OF SOFTWARE MODELS**

Although in Immunology there is more than 20-years tradition of software and mathematical modeling, very few of them have been the object of further exploitation once published and made available (30). Models are often under-used because experimentalists can be reluctant to entertain mathematical formalization and because published models are largely disposable: rapidly forgotten after being published, instead of providing a foundation to build upon. Moreover, the various expressions of these models with different mathematical descriptions, programing languages, software libraries and graphical packages, require much effort in understanding and running the software and prevent interoperability.

# **USING VISUAL LANGUAGE TO COMMUNICATE AND EXECUTE MODELS**

Immunologists often conceptualize the dynamical evolution of their systems in terms of "state-transitions" of biological objects and do it by means of personalized and informal graphical illustrations. Thus, the adoption of a more formal and standard type of state-transition diagram could improve the current situation to not only help biologists to better understand each other but also to facilitate the production and the reading of software code executing these visual transitions, at level of populations or agents (31).

Thus, in this paper, we promote the development and usage of a visual, computational approach more comprehensible than mathematical equations and programing instructions. This should improve our understanding of lymphocyte dynamics, the exchange on this understanding and simplify the implementation of models by non-specialists delivered from the production of executable code or mathematical equations, to concentrate to *in silico* experiments.

#### **LEVEL OF ABSTRACTION AND MULTI-SCALE MODELING**

A model describes a complex system from the "real word" and thus requires abstraction. This abstraction is performed as the immunologist decides about an experimental protocol in order to observe selected objects at different scales and to follow them in time and space. For example, the capacities of the immune system to preserve the homeostasis and to provide rapid adaptation to an antigen and anamnestic responses can be observed at the organism level, through physiological or pathological clinical observations that relate to lower scale levels. At molecular level, the somatic generation of the diversity of an immuno-receptor, as the TCR, allows for a dynamic network of interactions with antigens. At the cell level, this leads to clonal selection, activation and division. At the organ level, the fluidity of the system insures constant tissue redistribution of cells and molecules, cell migration from thymus to spleen and lymph nodes.

Thus, models of lymphocyte population dynamics and turnover consist in reconstructing the components or "entities" of the system across various scales, from molecules to organisms, to determine the relations/interactions through space (varying from micrometer to meters) and "processes" through time (varying between microseconds to years) as explained below. However, the formalization and abstraction of the immune objects as entities undergoing processes, with the help of spatial and dynamic ontologies, respectively defined as SNAP and SPAN (32), as well as cell/molecule interactions (33), is rarely done, maintaining a language-barrier between biologists and theoreticians. In the following, some examples will be given to help the immunologists with the transition between current mathematical models to computer ones and with the terms currently used in modeling.

#### **DEFINE ENTITIES, STATES, LOCATION, INTERACTIONS, GRANULARITY**

The immune "entities" could be described according to the language used by the modeler. A cell exists in one "state": it could be quiescent or in a given phase of the cell cycle or dead. In addition the phenotype and/or a function of a cell define a given state, as CD4 helper T cells. Cells are "located" in various tissues and are in "relation" with other entities. Finally, cells can be considered at various level of "granularity." For example, T lymphocyte populations are"aggregation"of T cells according to criteria of phenotype, structure or function, although heterogeneity still prevails inside these populations at lower granularity. Accordingly, cells can be modeled at population level (with continuous model as ordinary equation) or at cell level according to space (with discrete model as automata or multi-agent system).

#### **DEFINE PROCESSES**

According to ontologies, cells participate to various processes, such as division, activation, differentiation,interaction, clonotype selection, apoptosis or migration. According to the states of the cells, their evolution can thus be modeled as "state-transition" that can be applied to various processes in parallel: for example, a thymocyte can differentiate while migrating in cortex and medulla. Note that processes at other levels like molecular or organ levels can similarly be described and modeled. Finally, all these process will determine the global cell dynamics and turnover.

# **THE UNIFIED MODELING LANGUAGE FOR HIGH-LEVEL MODELING**

"High-level programing languages" are based on abstraction and use of natural language that is easier to understand as compared to "low level programing language," based on codes. Thus, visual modeling language considers biological-object as conceptual abstract-objects that endure processes. The level of abstraction

allowed by these diagrams makes possible to distinguish more easily the"entities"as T cells and the"processes" that occur at different levels such as differentiation, migration and cell cycle. Moreover, such "state-transition diagrams" allow computing parallel pathways at various scales to avoid redundancy that is inherent in the formal description of multi-level, heterogeneous and concurrent systems and to model heterogeneity in a very simplified and economical form (as compared to mathematical equations) (31).We have thus used the well-established Unified Modelling Language (UML – a sofware standard) that still remains approachable to the lab-immunologist, convenient for the theoretician and that can be directly adopted for the high-level graphical depiction (31, 34). The adoption of UML state-transition diagrams that transcends any programing language or computer platform, will allow both experimentalists and theorists to work together at a higher level than writing software code or mathematical equations. This final step is progressively more and more automatized out of the diagrams. Example of basic transformations of mathematical equations into state-transition diagrams including elementary, parallel, independent, or coupled state transition have been given (31), to familiarize biologists with the general approach.

# **REFACTORING FROM LOW LEVEL (CODE, EQUATIONS) TO HIGH-LEVEL (DIAGRAM) MODELING LANGUAGES**

To convince the immunologist of the feasibility of this approach as well as the benefit gained by adopting it, we sketch in the rest of the paper how existing "low level programed" immune models should gain in readability and accessibility by adopting a "highlevel graphical" representation under the form of "state-transition diagrams."We present a "refactoring" of two published models of T cell biology in the thymus. Refactoring consists in restructuring the code or equations of a model to improve its expression, readability and extensibility, without changing its external behavior. One model consists in cell population differentiation modeling with differential equations (continuous model). The other one is a discrete model. Originally it was an automata model consisting of a discrete lattice, where each site (cell) in a given state, follows some rules in space and time that depends on local neighbors (18). It has been refactorized as an agent-based model (ABM), depicting individual cell behavior through thymus differentiation and migration. It would be much too long and redundant to describe in details the behavior of these two models. We do not pretend here to modify at all the results obtained by the running of these models (the readers interested in these results are invited to access the original papers). We have just reshape them into a statetransition diagrammatic form that allows execution of simulations reproducing the original results with similar parameter values.

# **POPULATION-BASED MODEL DESCRIBING THE CONVEYOR-BELT T CELL DIFFERENTIATION IN THYMUS**

The original model (8) is a compartmentalized ordinary differential equation (ODE) model, rather complex to read and manage by immunologists. This model reflects the conceptual "conveyor belt" model of thymic T cell differentiation, schematically represented by immunologists by the continuous ordered transition of cells through the different stages with time (35–37). **Figure 1** represents a biological schema, originally published and the "state-transition" description of the model (in a UML state-transition diagram) as

it is proposed now. Although the original model is composed of 30 differential equations, the whole mathematical description and the code that captures it, can easily be deduced and regenerated from the **Figure 1**. Conversely, the mathematical equations can be automatically generated from the state-transition diagram as previously described (38). In essence, the model is summarized by the input, transition and output from the thymus, by "parallel processes" that occur concomitantly, as differentiation, cell cycle, proliferation and death, and by exit from the thymus. Note that these parallel processes concern various biological levels and time scales. The "differentiation" process represents each stage of the conveyor belt, from double negative (DN) to single positive (SP) cells, with flows into and from a particular stage according to the general equation:

$$d\mathfrak{x}\_{i}/dt = 2\mathfrak{y}\ p\mathfrak{x}\_{i-1} - (\mathfrak{p} + d + \mathfrak{u}(i))\mathfrak{x}\_{i}$$

*p*, *d*, and *u* represent proliferation, death, and differentiation, respectively, and *xi* represents the *i*th stage. The model mainly consists of constant hematopoietic progenitor influx in thymus (Sn); differentiation between thymocyte developmental phenotypes DN and double-positive (DP) cells differentiation into CD4+ or CD8+ cells, then egression of SP stage, either, to the periphery [Us4(i), Us8(i)]; proliferation [Pn(y), Pp(y), and Ps(y)]; positive and negative selection (a4, a8); and natural cell death (Dn, Dp, and Ds). In parallel, the "cell cycle" is represented: the cell switches between quiescence (G0) and cycle with division (S/M). The parameter γ set to 1 represents the cell division into two daughter cells. There is the possibility to induce a perturbation into the system through the specific depletion of T cells entering the S/M cycle phase (*p*), if γ is set to 0. This represents the presence or absence of a pharmacogenetic conditional treatment by ganciclovir that induces apoptosis related to the incorporation of a nucleotide analog during DNA elongation. This rule applies to all cell populations except in late DP quiescent cells as indicated in the schema. The "proliferation" depicts that the daughters of a proliferating cell transit into the next generational compartment, except during treatment (γ = 0) when dividing cells die by apoptosis and are lost from the model. The parameter *u* is an increasing function of generation (G1 to *n*), making cells more likely to differentiate between phenotypic compartments as they progress through the cell generations.

The parallelism in this graphical model largely simplifies the original formulation of ODEs while remaining faithful to it. Hierarchy and compound states are present again clearly reducing the diagram clutter. Other representation of the model depicting the differentiation with linear cell generations is also possible (31).

### **DISCRETE MODEL DESCRIBING THE DIFFERENTIATION ALONG THE MIGRATION OF T CELLS IN THE THYMUS IN A 2-D ENVIRONMENT**

The original model (20) is a discrete-based "cellular automata" computer model. The model depicts the behavior of individual thymocytes that evolve in the 2-D epithelial cell network, guided by chemokines gradients. The current model (**Figure 2**) is now an Agent-Based Model (ABM). Again, the interested reader is referred to the original paper for a detailed understanding of the simulation. Although available for download, the 40 pages of FORTRAN source code are far from easy to understand. After refactoring, the transition rules of any agent (thymocytes) map onto a parallel

state-transition diagram. The conception of the state-transition diagram, as done here, should considerably improve the understanding of the model (evenfor the original programmer) allowing the researchers to progress further with the existing simulator. A complete description of the model includes additional implementation details that are abstractions of the mechanisms behind cell decisions to differentiate. The parallel state-transition diagram represents the different simultaneous transitions taking place in the model and coded as various cellular automata rules: a cell in the model as it differentiates transits in successive states, it may be bound or not to thymic epithelial cells via TCR/MHC,it moves and may be located into one of several anatomical compartments of the thymus. As indicated in the boxes, the gradient of chemokines (k) orients the migration of cells for each specific stage of differentiation, and chemokines are localized in specific areas of the thymus. A T cell sums up the time and number of interactions with the same or different epithelial cells. This sum value determines whether the T cell is positively or negatively selected. DP and SP phenotypes have their own threshold parameters. They cannot differentiate until they cross a threshold number of interactions. If, after a given time, a DP cell has not reached this threshold, it enters apoptosis by neglect. If the time is too long, it is negatively selected. A threshold parameter simulates the phenotype decision. With this "signal-duration" hypothesis, long duration TCR-MHC interactions promote the CD4 phenotype and short duration promotes the CD8 phenotype.

It is important to notice that both refactored models, population-based and ABM can now be compared, are directly executable and can provide simulations of physiology, pathologies, and treatment, while not being the scope of this paper. Moreover, their parameters can be automatically tuned to fit experimental data. Any ABM model can also be run as a population version to save time in simulation (McEwan et al., manuscript in preparation). Moreover, the flexibility of these diagrams allows assembling the parts of the biological puzzle piece after piece and improving the models.

# **PERSPECTIVES**

As shown here, state-transition diagrams can represent highlevel semantics suitable to clarify immunological concepts and to aid communication among interdisciplinary researchers. It can also represent low levels quantitative information suitable for individual-based ABM and population-based ODE modeling. Organization of immune knowledge using a standardized, diagrammatic formal language should greatly improve knowledge integration at multi-scale levels and sharing between experimentalist and theoretician collaborators, rendering their software more readable, scalable, and usable. We are currently working on ways of automatically generating executable code out of these state-transition diagrams. State-transition diagrams supports the extension and interoperability of published models. This will help for dynamic computational modeling of lymphocyte behavior

in health and diseases, and for "*in silico*" experiments to predict and explain the puzzling T cells dynamics and the effect of immunological perturbations.

# **AUTHOR CONTRIBUTIONS**

Véronique Thomas-Vaslin and Hugues Bersini designed the project and wrote the article,Adrien Six and Jean-Gabriel Ganascia participated to discussions.

# **REFERENCES**

(2013). p. 1145–49. doi:10.1007/ 978-1-4419-9863-7\_95

**ACKNOWLEDGMENTS**

The authors thank C. McEwan for computer modeling developments and discussions, B. Zerath, B. Horsmans for early refactoring of models, W. Chaara, and E. Mariotti-Ferrandiz from critical reading of the manuscript. Hugues Bersini is an invited Professor at UPMC. Financial supports: Emergence-UPMC program, PHC Tournesol EGIDE, CNRS PEPS BMI, Investissements d'Avenir program (ANR-11-IDEX-0004-02), LabEx Transimmunom.

from: http://teoriasistemas. mex.tl/imagesnew/6/6/4/5/2/

immature (CD4+8+) thymocytes in normal and genetically modified mice. *J Immunol* (1995) **154**(10):5103–13.

**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: 06 June 2013; paper pending published: 07 July 2013; accepted: 09 September 2013; published online: 01 October 2013.*

*Citation: Thomas-Vaslin V, Six A, Ganascia J-G and Bersini H (2013) Dynamical and mechanistic reconstructive approaches of T lymphocyte dynamics: using visual modeling languages to bridge the gap between immunologists, theoreticians, and programmers. Front. Immunol. 4:300. doi: 10.3389/fimmu.2013.00300*

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

*Copyright © 2013 Thomas-Vaslin, Six, Ganascia and Bersini. 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, providedthe 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.*

# A mechanistic model for naive CD4 T cell homeostasis in healthy adults and children

# **Tharindi Hapuarachchi, Joanna Lewis and Robin E. Callard \***

Institute of Child Health and CoMPLEX, University College London, London, UK

#### **Edited by:**

Miles Davenport, University of New South Wales, Australia

#### **Reviewed by:**

Grant Lythe, University of Leeds, UK Mark Dowling, Walter and Eliza Hall Institute, Australia

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

Robin E. Callard, Institute of Child Health and CoMPLEX, University College London, 30 Guilford Street, London WC1N 1EH, UK e-mail<EMAIL_ADDRESS>
The size and composition of the T lymphocyte compartment is subject to strict homeostatic regulation and is remarkably stable throughout life in spite of variable dynamics in cell production and death during T cell development and immune responses. Homeostasis is achieved by careful orchestration of lymphocyte survival and cell division. New T cells are generated from the thymus and the number of peripheral T cells is regulated by controlling survival and proliferation. How these processes combine is however very complex. Thymic output increases in the first year of life and then decreases but is crucial for establishing repertoire diversity. Proliferation of new naive T cells plays a crucial role for maintaining numbers but at a potential cost toTCR repertoire diversity. A mechanistic twocompartment model of T cell homeostasis is described here that includes specific terms for thymic output, cell proliferation, and cell death of both resting and dividing cells. The model successfully predicts the homeostatic set point for T cells in adults and identifies variables that determine the total number of T cells. It also accurately predicts T cell numbers in children in early life despite rapid changes in thymic output and growth over this period.

**Keywords: naive T cells, homeostasis, CD4 T cells, mathematical modeling, mechanistic modeling, children**

# **INTRODUCTION**

The naive T cell compartment in humans is generated early in development by the thymus and then maintained throughout life by continued export from the thymus and cell division in the periphery. In adult humans, the naive T cell compartment is comprised of roughly 1011 cells circulating between the blood and the peripheral lymphoid organs. It is estimated to comprise at least 108 different T cell receptor specificities (1) providing a broad spectrum of protection in a diverse pathogen environment. The size and composition (T cell receptor diversity) of the naive T cell compartment are subject to strict homeostatic regulation and are remarkably stable throughout adult life despite changing rates of cell production and death during T cell development and immune responses (2, 3). Homeostasis is achieved by control of lymphocyte survival and cell division. Naive T cell survival and peripheral cell division depends on access to the cytokine IL7 (4–7) and TCR signals (8, 9) through contact with self-peptide MHC (spMHC) on dendritic cells (10). In lymphoreplete mice, naive T cells are largely non-cycling (11) whereas homeostatic cell division plays an important role in maintaining naive T cell homeostasis in humans, where cell division is evident in the naive pool (12, 13).

In children, homeostatic control of the T cell compartment may be affected by both the growth of the child with the accompanying increased blood volume (14) and changes in thymic output, which increases to a maximum over the first year of life and then declines to reach an approximately steady level by the age of about 20 years (15). As a result, the CD4 naive T cell count (cells/µl) in children declines over the first 10–20 years of life whereas the total number of naive CD4 cells increases as the child grows (**Figure 1**).

This raises important questions about whether the homeostatic mechanisms themselves change during early life or whether the numbers of naive CD4 T cells observed are determined only by the changes in thymic output and growth.

To date, our understanding of the processes controlling survival and proliferation of T cells has been largely qualitative and detailed quantitative knowledge of how homeostatic responses result in the observed equilibrium of the T cell pool with a given size and composition is lacking. Here, a two-compartment mathematical model of homeostasis is presented incorporating specific terms for thymic export into the naive CD4 compartment, rates of entry into cell division and death (survival) rates for both the resting and dividing cell compartments. In this sense, the model can be considered as mechanistic in comparison to empirical or descriptive models where the parameters have no direct biological meaning. The results illustrate the importance of T cell dynamics for the maintenance of constant naive CD4 T cell numbers in adults and the growth of the T cell compartment in children.

### **MATERIALS AND METHODS**

### **A MODEL OF NAIVE T CELL HOMEOSTASIS**

T cell homeostasis can be described using a two-compartment model of resting and dividing cells with input from the thymus into the resting compartment as shown in **Figure 2** (16, 17). In this model we will consider only naive CD4 T cells assuming no antigenic stimulation and maturation of naive to memory cells. The same model could in principle also be applied to memory cells and CD8 T cells.

This model can be expressed mathematically using two coupled ordinary non-linear differential equations:

An appropriate form for the rate of entry into cell cycle λ is (19)

$$
\lambda = \lambda\_0 \exp\left[-N\left(t\right)/\varepsilon\right],
$$

= θ + 2r*Y* − λ(*X*+*Y* )*X* − δ(*X*+*Y* )*X dt* = λ(*X*+*Y* )*X* − r*Y* − µ(*Y* )