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–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 co-receptors 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.

Materials and Methods

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: n1, the population of pre-selection DP thymocytes (double positive), that are TCRβlow and CD69low (26), n2, the population of post-selection DP thymocytes, that are TCRβ+ and CD69high (26), and n3, the population of mature SP (single positive) thymocytes.

We assume that DN thymocytes differentiate to become pre-selection 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, n1, n2, and n3, are involved in the following selection events in the cortex and the medulla (see Figure 1):

Figure 1

Thymic development as hypothesized in the first model. The flux, ϕ, 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 two fates: further differentiation into the SP pool (φ2) or death by apoptosis (μ2). Finally, SP thymocytes have three fates: maturation and exit into the periphery (φ3), death by apoptosis (μ3), or proliferation (λ3).

– flux of DN thymocytes into compartment n1,

– differentiation from pre-DP (n1) to post-DP (n2) thymocytes induced by TCR signal,

– death by neglect of pre-DP thymocytes due to lack of (or weak) TCR signal,

– differentiation from post-DP (n2) to SP (n3) thymocytes sustained by intermediate TCR signal,

– apoptosis of post-DP (n2) thymocytes due to strong TCR signal,

– exit of SP thymocytes (n3) to the periphery (thymic maturation),

– proliferation of SP thymocytes (n3) in the medulla, and

– apoptosis of SP (n3) thymocytes due to strong 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:

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 . 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 , 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 n1(0), n2(0), n3(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 and , 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: n1, the population of pre-selection DP (double positive) thymocytes, n2, the population of post-selection DP thymocytes, n4, the population of mature CD4+ SP (single positive) thymocytes, and n8, 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, n1, n2, n4, and n8, with n3 = n4 + n8, are involved in the following selection events in the cortex and the medulla (see Figure 2):

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

– flux of DN thymocytes into compartment n1,

– differentiation from pre-DP (n1) to post-DP (n2) thymocytes induced by TCR signal,

– death by neglect of pre-DP thymocytes due to lack of (or weak) TCR signal,

– differentiation from post-DP (n2) to CD4+ SP (n4) sustained by intermediate TCR signal,

– differentiation from post-DP (n2) to CD8+ SP (n8) sustained by intermediate TCR signal,

– apoptosis of post-DP (n2) thymocytes due to strong TCR signal,

– exit of CD4+ SP thymocytes (n4) to the periphery (thymic maturation),

– exit of CD8+ SP thymocytes (n8) to the periphery (thymic maturation),

– proliferation of CD4+ SP thymocytes (n4) in the medulla,

– proliferation of CD8+ SP thymocytes (n8) in the medulla,

– apoptosis of CD4+SP (n4) thymocytes due to strong TCR signal, and

– apoptosis of CD8+SP (n8) thymocytes due to strong TCR signal.

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 and . 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 , 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: where n1(0), n2(0), n4(0), n8(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 , , and , as it is the unique stable steady state.

Results

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: . 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, τ, 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).

Table 1

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

Mouse	n1*(pre-DP) (cells)	n2*(post-DP) (cells)	n3*(SP) (cells)	n4*(SP CD4) (cells)	n8*(SP CD8) (cells)
1	82.58 × 106	9.30 × 106	18.36 × 106	13.85 × 106	4.51 × 106
2	142.19 × 106	19.94 × 106	26.20 × 106	18.73 × 106	7.46 × 106
3	89.00 × 106	5.98 × 106	15.98 × 106	11.88 × 106	4.10 × 106
4	29.32 × 106	2.09 × 106	5.61 × 106	4.40 × 106	1.21 × 106
5	29.32 × 106	2.09 × 106	5.61 × 106	4.40 × 106	1.21 × 106
6	51.26 × 106	5.93 × 106	9.01 × 106	6.85 × 106	2.16 × 106
7	64.48 × 106	6.81 × 106	11.64 × 106	9.03 × 106	2.61 × 106
8	218.94 × 106	15.42 × 106	40.20 × 106	29.46 × 106	10.74 × 106
Mean	88.39 × 106	8.45 × 106	16.57 × 106	12.33 × 106	4.25 × 106
Standard deviation	60.11 × 106	5.89 × 106	11.05 × 106	7.94 × 106	3.12 × 106

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

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

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

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 × 106 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, , and the fact that we can obtain an estimate for φ3. We have chosen ϕout to be 2.5 × 106 cells per day (14, 30).

Given τ3, φ3, and the fact that we can obtain an estimate for μ3.

Given τ1, , and the fact that we can obtain an estimate for ϕ.

We also have which, in principle, allows us to estimate φ1. We make use of linear regression techniques to do so (31, 32).

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

We minimize E(a1) with respect to a1, that is . Solving for a1, we obtain:

Given a1, we can then estimate φ1 from the equation

Given τ1, φ1, and the fact that we can obtain an estimate for μ1.

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 , with time measured in days. We sample equally spaced values for λ3, and for each value, we compute . The ratio 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 p-values for a1 and a2 are given by 7.57 × 10−3 and 6.85 × 10−3, respectively (both smaller than the significance level α = 0.05).

Figure 3

Linear regression plots for the first model.

Given τ2, φ2, and the fact that , we can obtain an estimate for μ2.

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(t = 21 days) belongs to the interval , with i = 1, 2, 3, and where is the (experimental) mean number of cells in compartment i, and σ is the (experimental) standard deviation in compartment i, as given in Table 1.

We obtain the following parameter values: and

These parameters imply the following thymic selection rates:

Death rates

9.7 × 105 cells/h die by neglect in compartment 1 (), 4.8 × 105 cells/h die by negative selection in compartment 2 (), and 6.9 × 104 cells/h die by negative selection in compartment 3 ().

Differentiation rates

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

Proliferation rate

1.3 × 105 cells/h proliferate in compartment 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:

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: . As described in Section 3.1, we also need residence times for each population subset, τ, 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 μ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:

Given τ4, μ4, and the fact that we can obtain an estimate for ξ4.

In the same way, given τ8, μ8, and the fact that we can obtain an estimate for ξ8.

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.

Given τ4, φ4, and the fact that we can compute the fraction 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).

Figure 4

Linear regression plots for the second model.

In Section 3.1, we obtained an estimate for the mean of λ3, and we know that As before, we can compute the fractions and by linear regression (see 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 a3, a4, and a5 are given by 8.43 × 10−3, 3.33 × 10−7, and 4.56 × 10−8, respectively (smaller than the significance level).

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: and

These parameters imply the following thymic selection rates:

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

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

Proliferation rates

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

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:

These probabilities imply that the probability to exit the thymus as a mature CD4 thymocyte (that has already reached the medulla) is given by , which is 45.1%, and the probability to exit as a mature CD8 thymocyte (that has already reached the medulla) is given by , 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: where a1, a2 are the coefficients of the linear regression of and , respectively, and is the experimental mean value of n.

We perturb each entry of the vector θ by adding and subtracting 10% of its value. Therefore, we now have two values for θ, equal to and . Consequently, we have 28 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: where a3, a4, a5 are the coefficients of the linear regression of , , and , 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 sets of θ, with 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 and minimum–maximum interval ranges, are given in Table 2.

Table 2

Means, 95% trimmed and minimum–maximum intervals of the model parameters.

Parameter	Mean value	95% Trimmed interval	Minimum–maximum interval range
ϕ	35.86 × 106 cells/day	(35.65 × 106, 35.07 × 106) cells/day	(28.93, 43.21 × 106) cells/day
φ1	0.139 day−1	(0.138, 0.140) day−1	(0.112, 0.167) day−1
φ2	0.136 day−1	(0.134, 0.139) day−1	(0.041, 0.274) day−1
φ4	0.140 day−1	(0.136, 0.145) day−1	(0.060, 0.264) day−1
φ8	0.134 day−1	(0.129, 0.138) day−1	(0.010, 0.214) day−1
μ1	0.265 day−1	(0.263, 0.267) day−1	(0.196, 0.333) day−1
μ2	1.372 day−1	(1.365, 1.378) day−1	(1.083, 1.618) day−1
μ4	0.040 day−1	n/a	(0.036, 0.044) day−1
μ8	0.110 day−1	n/a	(0.099, 0.121) day−1
λ4	0.181 day−1	(0.179, 0.184) day−1	(0.116, 0.226) day−1
λ8	0.085 day−1	(0.080, 0.090) day−1	(0.078, 0.092) day−1
ξ4	0.231 day−1	(0.230, 0.233) day−1	(0.229, 0.233) day−1
ξ8	0.152 day−1	(0.150, 0.154) day−1	(0.149, 0.155) day−1

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: We present in Table 3 the variability of the selection rates.

Table 3

Selection rate values (initial and after perturbation) and their variability intervals.

Rate	Initial value (%)	After perturbation (%)	Δ Value (%)	Δ Min–max (%)
p1	65.8	65.66	±0.76	±20.55
p2	91.7	90.98	±0.49	±17.28
p3	22.91	24.07	±0.90	±28.22
q3	42.24	41.74	±0.99	±30.53
s4	4.69	8.51	±0.80	±15.16
s8	3.61	8.13	±0.79	±15.03
p4	8.59	8.84	±0.37	±4.91
q4	46.29	40.07	±1.29	±20.49
p8	32.12	31.71	±2.25	±35.37
q8	27.0	24.5	±3.58	±64.81

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 , (ii) the CD4:CD8 ratio in the SP pool is given by , and (iii) the CD4:CD8 ratio of positive selection in the SP pool (differentiation from SP to peripheral early thymic emigrants) is given by . 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 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, n(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).

Figure 5

Time evolution of the thymocyte populations in the second model. The different trajectories correspond to the parameter values and ranges described in Table 2.

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 intra-epithelial 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.

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.