From pre-DP, post-DP, SP4, and SP8 Thymocyte Cell Counts to a Dynamical Model of Cortical and Medullary Selection

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.


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)(5)(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).

www.frontiersin.org
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)(23)(24)(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 −/− Nur77 GFP 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.

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 Nur77 GFP 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 Nur77 GFP 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 CD69 low (26), n 2 , the population of post-selection DP thymocytes, that are TCRβ + and CD69 high (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) φ. We further 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): • ∅ φ − → n 1 -flux of DN thymocytes into compartment n 1 , • n 1 ϕ 1 − → n 2 -differentiation from pre-DP (n 1 ) to post-DP (n 2 ) thymocytes induced by TCR signal, • n 1 µ 1 − → ∅ -death by neglect of pre-DP thymocytes due to lack of (or weak) TCR signal, 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: Frontiers in Immunology | T Cell Biology 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: 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 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): 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: 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.

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): • ∅ φ − → n 1 -flux of DN thymocytes into compartment n 1 , • n 1 ϕ 1 − → n 2 -differentiation from pre-DP (n 1 ) to post-DP (n 2 ) thymocytes induced by TCR signal, www.frontiersin.org FIGURE 2 | Thymic development as hypothesized in the second model. We assume there is a flux, φ, that represents the differentiation of DNs into pre-DPs. Pre-DP thymocytes have two fates: further differentiation into the post-DP pool (ϕ 1 ) or death by neglect (µ 1 ). Post-DP thymocytes have three fates: further differentiation into the SP pool as CD4 SPs (ϕ 4 ), or CD8 SPs (ϕ 8 ), or death by apoptosis (µ 2 ). Finally, CD4 or CD8 SP thymocytes have three fates: maturation and exit into the periphery (ξ 4 ) or (ξ 8 ), death by apoptosis (µ 4 ) or (µ 8 ), or proliferation (λ 4 ) or (λ 8 ).
• n 1 µ 1 − → ∅ -death by neglect of pre-DP thymocytes due to lack of (or weak) TCR signal, • n 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: The time evolution of the four populations can be described by the following set of ODEs: 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: 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 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: 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: Frontiers in Immunology | T Cell Biology 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.

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: 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)(28)(29) τ 1 = 60 h = 2.5 days , τ 2 = 16 h = 0.67 days , In order to derive estimates for the model parameters, we have carried out the following steps: 1. We make use of the experimentally determined mature SP thymocyte flux from the medulla to the periphery, which has been estimated to be 1-4 × 10 6 cells per day (14,26,30). This flux corresponds to about 1% of thymocytes leaving the thymus every day (30). Given this flux, which we denote by φ out , n * 3 , and the fact that φ out = ϕ 3 n * 3 , we can obtain an estimate for ϕ 3 . We have chosen φ out to be 2.5 × 10 6 cells per day (14,30). 2. Given τ 3 , ϕ 3 , and the fact that τ 3 = 1 µ 3 +ϕ 3 , we can obtain an estimate for µ 3 . 3. Given τ 1 , n * 1 , and the fact that n * 1 = φ τ 1 , we can obtain an estimate for φ. 4. We also have n * 2 = ϕ 1 τ 2 n * 1 , which, in principle, allows us to estimate ϕ 1 . We make use of linear regression techniques to do so (31,32).
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: We minimize E(a 1 ) with respect to a 1 , that is dE da 1 = 0. Solving for a 1 , we obtain: The bold font highlights the mean and the standard deviation from the individual mice data.

www.frontiersin.org
Given a 1 , we can then estimate ϕ 1 from the equation ϕ 1 = a 1 τ 2 . 5. Given τ 1 , ϕ 1 , and the fact that τ 1 = 1 µ 1 +ϕ 1 , we can obtain an estimate for µ 1 . 6. We are now left with three remaining parameters: ϕ 2 , µ 2 , and λ 3 . Given the experimental constraints on τ 1 , τ 2 , and τ 3 , we assume that the average time to proliferate, 1/λ 3 , cannot be <7 days. Therefore, we consider λ 3 to be constrained in the interval [1/7, τ −1 3 ], with time measured in days. We sample equally spaced values for λ 3 , and for each value, we compute is computed using the linear regression method described above (see Figure 3). In this way, we obtain an estimate for ϕ 3 . We note that the pvalues for a 1 and a 2 are given by 7.57 × 10 −3 and 6.85 × 10 −3 , respectively (both smaller than the significance level α = 0.05). 7. Given τ 2 , ϕ 2 , and the fact that τ 2 = 1 µ 2 +ϕ 2 , we can obtain an estimate for µ 2 . 8. From steps 6 and 7 above, we have generated (a table of) values for ϕ 2 and µ 2 , given a fiducial value for λ 3 in the interval The mice considered in the experimental study are 5.5-17 weeks old, and their thymus is in steady state (26). Thus, we expect that the parameter values can only be accepted if the corresponding system of ODEs attains steady state by 3 weeks. Therefore, 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). That is, we impose for the given parameter set that the mathematically predicted value n i (t = 21 days) belongs to the interval n * i ± σ i , with i = 1, 2, 3, and where n * i is the (experimental) mean number of cells in compartment i, and σ i is the (experimental) standard deviation in compartment i, as given in Table 1.

Proliferation rate
1.3 × 10 5 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: 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: = 42.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: 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 are 1 µ 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, the average residence times in each compartment, as described in Section 3.1, are given by: In order to derive estimates for the model parameters, we have carried out the following steps: 1. Given τ 4 , µ 4 , and the fact that τ 4 = 1 µ 4 +ξ 4 , we can obtain an estimate for ξ 4 . 2. In the same way, given τ 8 , µ 8 , and the fact that τ 8 = 1 µ 8 +ξ 8 , we can obtain an estimate for ξ 8 . 3. We are now left with four remaining parameters: ϕ 4 , ϕ 8 , λ 4 , and λ 8 . We know that ϕ 2 = ϕ 4 + ϕ 8 . We sample ϕ 4 in the interval [0, ϕ 2 ], where ϕ 2 is the mean value of the interval obtained in Section 2.1, and for each fiducial value for ϕ 4 , we compute the corresponding value for ϕ 8 . 4. Given τ 4 , ϕ 4 , and the fact that , we can compute the fraction a 3 = n * 2 n * 4 by linear regression (see Figure 4), and thus obtain an estimate for λ 4 . Note that we will reject values of λ 4 that imply the proliferation time is larger than 7 days (see Section 3.1). 5. In Section 3.1, we obtained an estimate for the mean of λ 3 , and we know that λ 4 n * 4 + λ 8 n * 8 = λ 3 n * 3 . As before, we can compute the fractions a 4 = n *  Figure 4), and thus obtain an estimate for λ 8 . Note that we will reject values of λ 8 that imply the proliferation time is larger than 7 days (see Section 3.1). We note that the p-values for a 3 , a 4 , and a 5 are given by 8.43 × 10 −3 , 3.33 × 10 −7 , and 4.56 × 10 −8 , respectively (smaller than the significance level). 6. From steps 3, 4, and 5 above, we have generated (a table of) values for ϕ 8 , λ 4 , and λ 8 , given a fiducial value for ϕ 4 in the interval [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: These parameters imply the following thymic selection rates:

www.frontiersin.org
We can compute the stringency of thymic selection, defined by the ratio: 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: 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: = 46.3% , 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%.

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: θ = (τ 1 , τ 2 , τ 3 , φ out , a 1 , a 2 ,n 3 ,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, andn i is the experimental mean value of n i . 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: τ 8 , µ 4 , µ 8 , a 3 , a 4 , a 5  , respectively. We have made use of the means of the following parameters of the first model: φ, ϕ 1 , µ 1 , ϕ 2 , µ 2 , λ 3 . We perturb each entry of the vector θ as described above. Therefore, we have 2 7 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 intervals 2 and minimum-maximum interval ranges, are given in Table 2.

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: 2 We define the 95% trimmed interval to be the result of the sensitivity analysis after trimming the lower and upper 2.5% of values.

Frontiers in Immunology | T Cell Biology
February 2014 | Volume 5 | Article 19 | 8 We present in Table 3 the variability of the selection rates.

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)(28)(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 www.frontiersin.org 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: 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 σ = 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: 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 10 3 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, n i (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).  Table 2.

Frontiers in Immunology | T Cell Biology
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 Nur77 GFP 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. 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 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 10 3 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.069N CD4 + T cells and 0.025N CD8 + T cells.