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.

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 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^−/−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.

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: ϕ, μ₁, μ₂, φ₁, and φ₂.

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₁, the population of pre-selection DP thymocytes (double positive), that are TCRβ^low and CD69^low (26), n₂, the population of post-selection DP thymocytes, that are TCRβ⁺ and CD69^high (26), and n₃, 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, n₁, n₂, and n₃, are involved in the following selection events in the cortex and the medulla (see Figure 1):

• $\varnothing \stackrel{\varphi }{\to }{n}_1$ – flux of DN thymocytes into compartment n₁,

• $n_1\stackrel{{\phi }_1}{\to }{n}_2$ – differentiation from pre-DP (n₁) to post-DP (n₂) thymocytes induced by TCR signal,

• $n_1\stackrel{{\mu }_1}{\to }\varnothing$ – death by neglect of pre-DP thymocytes due to lack of (or weak) TCR signal,

• $n_2\stackrel{{\phi }_2}{\to }{n}_3$ – differentiation from post-DP (n₂) to SP (n₃) thymocytes sustained by intermediate TCR signal,

• $n_2\stackrel{{\mu }_2}{\to }\varnothing$ – apoptosis of post-DP (n₂) thymocytes due to strong TCR signal,

• $n_3\stackrel{{\phi }_3}{\to }\mathrm{periphery}$ – exit of SP thymocytes (n₃) to the periphery (thymic maturation),

• $n_3\stackrel{{\lambda }_3}{\to }2n_3$ – proliferation of SP thymocytes (n₃) in the medulla, and

• $n_3\stackrel{{\mu }_3}{\to }\varnothing$ – apoptosis of SP (n₃) thymocytes due to strong TCR signal.

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 (φ₁) or death by neglect (μ₁). Post-DP thymocytes have two fates: further differentiation into the SP pool (φ₂) or death by apoptosis (μ₂). Finally, SP thymocytes have three fates: maturation and exit into the periphery (φ₃), death by apoptosis (μ₃), or proliferation (λ₃).

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{array}{c}\hfill \left\{\begin{array}{c}\frac{d{n}_1}{dt}=\mathrm{\varphi }-{\mathrm{\phi }}_1{n}_1-{\mathrm{\mu }}_1{n}_1,\hfill \\ \frac{d{n}_2}{dt}={\mathrm{\phi }}_1{n}_1-{\mathrm{\phi }}_2{n}_2-{\mathrm{\mu }}_2{n}_2,\hfill \\ \frac{d{n}_3}{dt}={\mathrm{\phi }}_2{n}_2-{\mathrm{\phi }}_3{n}_3-{\mathrm{\mu }}_3{n}_3+{\mathrm{\lambda }}_3{n}_3.\hfill \end{array}\right.\end{array}$ (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{\mathrm{\varphi }}{{\mathrm{\phi }}_1+{\mathrm{\mu }}_1},n_2^{*}=\frac{n_1^{*}{\mathrm{\phi }}_1}{{\mathrm{\phi }}_2+{\mathrm{\mu }}_2},n_3^{*}=\frac{n_2^{*}{\mathrm{\phi }}_2}{{\mathrm{\phi }}_3+{\mathrm{\mu }}_3-{\mathrm{\lambda }}_3}.$ (2)

This unique steady state exists if and only if φ₃ + μ₃ − λ₃ > 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=\left(\begin{array}{ccc}\hfill -\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right)\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill {\mathrm{\phi }}_1\hfill & \hfill -\left({\mathrm{\phi }}_2+{\mathrm{\mu }}_2\right)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\mathrm{\phi }}_2\hfill & \hfill -\left({\mathrm{\phi }}_3+{\mathrm{\mu }}_3-{\mathrm{\lambda }}_3\right)\hfill \\ \hfill \hfill \end{array}\right).$ (3)

A is also the Jacobian matrix at the steady state $n^{*}=\left(n_1^{*},n_2^{*},n_3^{*}\right)$, as the system of ODEs [equation (1)] is linear. The three eigenvalues of A are given by (as the matrix is lower triangular):

${\mathrm{\beta }}_1=-\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right),{\mathrm{\beta }}_2=-\left({\mathrm{\phi }}_2+{\mathrm{\mu }}_2\right),{\mathrm{\beta }}_3=-\left({\mathrm{\phi }}_3+{\mathrm{\mu }}_3-{\mathrm{\lambda }}_3\right).$

Therefore, the steady state [equation (2)] is stable, if and only if, φ₃ + μ₃ − λ₃ > 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{array}{llll}\hfill n_1\left(t\right)& =n_1^{*}+n_1\left(0\right)e^{-\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right)t},& \hfill & \\ \hfill n_2\left(t\right)& =n_2^{*}+\frac{n_1\left(0\right){\mathrm{\phi }}_1}{\left[\left({\mathrm{\phi }}_2+{\mathrm{\mu }}_2\right)-\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right)\right]}{e}^{-\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right)t}& \hfill & \\ \hfill & +n_2\left(0\right)e^{-\left({\mathrm{\phi }}_2+{\mathrm{\mu }}_2\right)t},& \hfill & \\ \hfill n_3\left(t\right)& =n_3^{*}+\frac{n_1\left(0\right){\mathrm{\phi }}_1}{\left[\left({\mathrm{\phi }}_2+{\mathrm{\mu }}_2\right)-\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right)\right]}\\ \hfill & \times \frac{{\mathrm{\phi }}_2}{\left[\left({\mathrm{\phi }}_3+{\mathrm{\mu }}_3-{\mathrm{\lambda }}_3\right)-\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right)\right]}{e}^{-\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right)t}& \hfill & \\ \hfill & +\frac{n_2\left(0\right){\mathrm{\phi }}_3}{\left[\left({\mathrm{\phi }}_3+{\mathrm{\mu }}_3-{\mathrm{\lambda }}_3\right)-\left({\mathrm{\phi }}_2+{\mathrm{\mu }}_2\right)\right]}{e}^{-\left({\mathrm{\phi }}_2+{\mathrm{\mu }}_2\right)t}& \hfill & \\ \hfill & +n_3\left(0\right)e^{-\left({\mathrm{\phi }}_3+{\mathrm{\mu }}_3-{\mathrm{\lambda }}_3\right)t},& \hfill & \end{array}$ (4)

where n₁(0), n₂(0), n₃(0) represent the initial conditions for the thymocyte populations. Note that in the late time limit, that is, if t → +∞ and φ₃ + μ₃ − λ₃ > 0, then $n_1\left(t\right)\to n_1^{*},n_2\left(t\right)\to n_2^{*}$ and $n_3\left(t\right)\to 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₁, the population of pre-selection DP (double positive) thymocytes, n₂, the population of post-selection DP thymocytes, n₄, the population of mature CD4⁺ SP (single positive) thymocytes, and n₈, 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₁, n₂, n₄, and n₈, with n₃ = n₄ + n₈, are involved in the following selection events in the cortex and the medulla (see Figure 2):

• $\varnothing \stackrel{\varphi }{\to }{n}_1$ – flux of DN thymocytes into compartment n₁,

• $n_1\stackrel{{\phi }_1}{\to }{n}_2$ – differentiation from pre-DP (n₁) to post-DP (n₂) thymocytes induced by TCR signal,

• $n_1\stackrel{{\mu }_1}{\to }\varnothing$ – death by neglect of pre-DP thymocytes due to lack of (or weak) TCR signal,

• $n_2\stackrel{{\phi }_4}{\to }{n}_4$ – differentiation from post-DP (n₂) to CD4⁺ SP (n₄) sustained by intermediate TCR signal,

• $n_2\stackrel{{\phi }_8}{\to }{n}_8$ – differentiation from post-DP (n₂) to CD8⁺ SP (n₈) sustained by intermediate TCR signal,

• $n_2\stackrel{{\mu }_2}{\to }\varnothing$ – apoptosis of post-DP (n₂) thymocytes due to strong TCR signal,

• $n_4\stackrel{{\xi }_4}{\to }\mathrm{periphery}$ – exit of CD4⁺ SP thymocytes (n₄) to the periphery (thymic maturation),

• $n_8\stackrel{{\xi }_8}{\to }\mathrm{periphery}$ – exit of CD8⁺ SP thymocytes (n₈) to the periphery (thymic maturation),

• $n_4\stackrel{{\lambda }_4}{\to }2n_4$ – proliferation of CD4⁺ SP thymocytes (n₄) in the medulla,

• $n_8\stackrel{{\lambda }_8}{\to }2n_8$ – proliferation of CD8⁺ SP thymocytes (n₈) in the medulla,

• $n_4\stackrel{{\mu }_4}{\to }\varnothing$ – apoptosis of CD4⁺SP (n₄) thymocytes due to strong TCR signal, and

• $n_8\stackrel{{\mu }_8}{\to }\varnothing$ – apoptosis of CD8⁺SP (n₈) thymocytes due to strong TCR signal.

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 (φ₁) or death by neglect (μ₁). Post-DP thymocytes have three fates: further differentiation into the SP pool as CD4 SPs (φ₄), or CD8 SPs (φ₈), or death by apoptosis (μ₂). Finally, CD4 or CD8 SP thymocytes have three fates: maturation and exit into the periphery (ξ₄) or (ξ₈), death by apoptosis (μ₄) or (μ₈), or proliferation (λ₄) or (λ₈).

We assume that all model parameters are positive, that is, μ₁, μ₂, μ₄, μ₈, φ₁, φ₂, φ₄, φ₈, ξ₄, ξ₈, ϕ, λ₄, λ₈ > 0, and note that the parameters and thymocyte populations of the first and second model are related by the following equations:

$\begin{array}{llll}\hfill & {\mathrm{\phi }}_2={\mathrm{\phi }}_4+{\mathrm{\phi }}_8,{\mathrm{\xi }}_4{n}_4+{\mathrm{\xi }}_8{n}_8={\mathrm{\phi }}_3{n}_3,& \hfill & \\ \hfill & {\mathrm{\mu }}_4{n}_4+{\mathrm{\mu }}_8{n}_8={\mathrm{\mu }}_3{n}_3,{\mathrm{\lambda }}_4{n}_4+{\mathrm{\lambda }}_8{n}_8={\mathrm{\lambda }}_3{n}_3.\end{array}$ (5)

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

$\begin{array}{c}\hfill \left\{\begin{array}{c}\frac{d{n}_1}{dt}=\mathrm{\varphi }-{\mathrm{\phi }}_1{n}_1-{\mathrm{\mu }}_1{n}_1,\hfill \\ \frac{d{n}_2}{dt}={\mathrm{\phi }}_1{n}_1-\left({\mathrm{\phi }}_4+{\mathrm{\phi }}_8\right)n_2-{\mathrm{\mu }}_2{n}_2,\hfill \\ \frac{d{n}_4}{dt}={\mathrm{\phi }}_4{n}_2-{\mathrm{\xi }}_4{n}_4-{\mathrm{\mu }}_4{n}_4+{\mathrm{\lambda }}_4{n}_4,\hfill \\ \frac{d{n}_8}{dt}={\mathrm{\phi }}_8{n}_2-{\mathrm{\xi }}_8{n}_8-{\mathrm{\mu }}_8{n}_8+{\mathrm{\lambda }}_8{n}_8.\hfill \end{array}\right.\end{array}$ (6)

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:

$\begin{array}{llll}\hfill n_1^{*}& =\frac{\mathrm{\varphi }}{{\mathrm{\phi }}_1+{\mathrm{\mu }}_1},n_2^{*}=\frac{n_1^{*}{\mathrm{\phi }}_1}{{\mathrm{\phi }}_4+{\mathrm{\phi }}_8+{\mathrm{\mu }}_2},& \hfill & \\ \hfill n_4^{*}& =\frac{n_2^{*}{\mathrm{\phi }}_4}{{\mathrm{\xi }}_4+{\mathrm{\mu }}_4-{\mathrm{\lambda }}_4},n_8^{*}=\frac{n_2^{*}{\mathrm{\phi }}_8}{{\mathrm{\xi }}_8+{\mathrm{\mu }}_8-{\mathrm{\lambda }}_8}.\end{array}$ (7)

This unique steady state exists if and only if ξ₄ + μ₄ − λ₄ > 0 and ξ₈ + μ₈ − λ₈ > 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=\left(\begin{array}{cccc}\hfill -\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right)\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill {\mathrm{\phi }}_1\hfill & \hfill -\left({\mathrm{\phi }}_4+{\mathrm{\phi }}_8+{\mathrm{\mu }}_2\right)\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\mathrm{\phi }}_4\hfill & \hfill -\left({\mathrm{\xi }}_4+{\mathrm{\mu }}_4-{\mathrm{\lambda }}_4\right)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\mathrm{\phi }}_8\hfill & \hfill 0\hfill & \hfill -\left({\mathrm{\xi }}_8+{\mathrm{\mu }}_8-{\mathrm{\lambda }}_8\right)\hfill \\ \hfill \hfill \end{array}\right).$ (8)

B is also the Jacobian at the steady state $n^{*}=\left(n_1^{*},n_2^{*},n_4^{*},n_8^{*}\right)$, as the system of ODEs is linear. The four eigenvalues of B are given by:

$\begin{array}{llll}\hfill & {\mathrm{\beta }}_1=-\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right),{\mathrm{\beta }}_2=-\left({\mathrm{\phi }}_4+{\mathrm{\phi }}_8+{\mathrm{\mu }}_2\right),& \hfill & \\ \hfill & {\mathrm{\beta }}_3=-\left({\mathrm{\xi }}_4+{\mathrm{\mu }}_4-{\mathrm{\lambda }}_4\right),{\mathrm{\beta }}_4=-\left({\mathrm{\xi }}_8+{\mathrm{\mu }}_8-{\mathrm{\lambda }}_8\right).& \hfill & \end{array}$

Therefore, the steady state equation (6) is stable if and only if ξ₄ + μ₄ − λ₄ > 0 and ξ₈ + μ₈ − λ₈ > 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{array}{llll}\hfill n_1\left(t\right)& =n_1^{*}+n_1\left(0\right)e^{-\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right)t},& \hfill & \\ \hfill n_2\left(t\right)& =n_2^{*}+\frac{n_1\left(0\right){\mathrm{\phi }}_1}{\left[\left({\mathrm{\phi }}_4+{\mathrm{\phi }}_8+{\mathrm{\mu }}_2\right)-\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right)\right]}{e}^{-\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right)t}& \hfill & \\ \hfill & +n_2\left(0\right)e^{-\left({\mathrm{\phi }}_4+{\mathrm{\phi }}_8+{\mathrm{\mu }}_2\right)t},& \hfill & \\ \hfill n_4\left(t\right)& =n_4^{*}+\frac{n_1\left(0\right){\mathrm{\phi }}_1}{\left[\left({\mathrm{\phi }}_4+{\mathrm{\phi }}_8+{\mathrm{\mu }}_2\right)-\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right)\right]}\hfill & \hfill \\ \hfill & \times \frac{{\mathrm{\phi }}_4+{\mathrm{\phi }}_8}{\left[\left({\mathrm{\xi }}_4+{\mathrm{\mu }}_4-{\mathrm{\lambda }}_4\right)-\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right)\right]}{e}^{-\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right)t}\hfill \\ \hfill & +\frac{n_2\left(0\right){\mathrm{\xi }}_4}{\left[\left({\mathrm{\xi }}_4+{\mathrm{\mu }}_4-{\mathrm{\lambda }}_4\right)-\left({\mathrm{\phi }}_4+{\mathrm{\phi }}_8+{\mathrm{\mu }}_2\right)\right]}{e}^{-\left({\mathrm{\phi }}_4+{\mathrm{\phi }}_8+{\mathrm{\mu }}_2\right)t}\hfill \\ \hfill & +n_4\left(0\right)e^{-\left({\mathrm{\xi }}_4+{\mathrm{\mu }}_4-{\mathrm{\lambda }}_4\right)t},\hfill \\ \hfill n_8\left(t\right)& =n_8^{*}+\frac{n_1\left(0\right){\mathrm{\phi }}_1}{\left[\left({\mathrm{\phi }}_4+{\mathrm{\phi }}_8+{\mathrm{\mu }}_2\right)-\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right)\right]}\hfill \\ \hfill & \times \frac{{\mathrm{\phi }}_4+{\mathrm{\phi }}_8}{\left[\left({\mathrm{\xi }}_8+{\mathrm{\mu }}_8-{\mathrm{\lambda }}_8\right)-\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right)\right]}{e}^{-\left({\mathrm{\phi }}_1+{\mathrm{\mu }}_1\right)t}\hfill \\ \hfill & +\frac{n_2\left(0\right){\mathrm{\xi }}_8}{\left[\left({\mathrm{\xi }}_8+{\mathrm{\mu }}_8-{\mathrm{\lambda }}_8\right)-\left({\mathrm{\phi }}_4+{\mathrm{\phi }}_8+{\mathrm{\mu }}_2\right)\right]}{e}^{-\left({\mathrm{\phi }}_4+{\mathrm{\phi }}_8+{\mathrm{\mu }}_2\right)t}\hfill \\ \hfill & +n_8\left(0\right)e^{-\left({\mathrm{\xi }}_8+{\mathrm{\mu }}_8-{\mathrm{\lambda }}_8\right)t},\hfill \end{array}$ (9)

where n₁(0), n₂(0), n₄(0), n₈(0) represent the initial conditions for the thymocyte populations. Note that in the late time limit, that is, if t → +∞ and ξ₄ + μ₄ − λ₄ > 0 and ξ₈ + μ₈ − λ₈ > 0, then $n_1\left(t\right)\to n_1^{*},n_2\left(t\right)\to n_2^{*}$, $n_4\left(t\right)\to n_4^{*}$, and $n_8\left(t\right)\to 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:

${\mathrm{\tau }}_i=\frac{1}{{\mathrm{\phi }}_i+{\mathrm{\mu }}_i},\mathrm{for}i\in \left\{1,2,3\right\}.$

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

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

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

$\begin{array}{llll}\hfill & {\mathrm{\tau }}_1=60\mathrm{h}=2.5\mathrm{days},{\mathrm{\tau }}_2=16\mathrm{h}=0.67\mathrm{days},& \hfill & \\ \hfill & {\mathrm{\tau }}_3=96\mathrm{h}=4\mathrm{days}.& \hfill & \end{array}$

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⁶ 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 ${\mathrm{\varphi }}_{\mathrm{out}}={\mathrm{\phi }}_3{n}_3^{*},$ we can obtain an estimate for φ₃. We have chosen ϕ_out to be 2.5 × 10⁶ cells per day (14, 30).

2. Given τ₃, φ₃, and the fact that ${\mathrm{\tau }}_3=\frac{1}{{\mathrm{\mu }}_3+{\mathrm{\phi }}_3},$ we can obtain an estimate for μ₃.

3. Given τ₁, $n_1^{*}$, and the fact that $n_1^{*}=\mathrm{\varphi }{\mathrm{\tau }}_1,$ we can obtain an estimate for ϕ.

4. We also have $n_2^{*}={\mathrm{\phi }}_1{\mathrm{\tau }}_2{n}_1^{*},$ which, in principle, allows us to estimate φ₁. We make use of linear regression techniques to do so (31, 32).

Let us introduce a₁ by the following equation, $a_1=\frac{n_2^{*}}{n_1^{*}}$, and make use of the experimental data to write: $n_2^{*,i}=a{n}_1^{*,i}+{\mathrm{\epsilon }}^i$, for i = 1, 2, …, 8. Thus, the squared error is given by:

$E\left(a_1\right)={\displaystyle \sum _{i=1}^8}{\left(n_2^{*,i}-a_1{n}_1^{*,i}\right)}^2.$

We minimize E(a₁) with respect to a₁, that is $\frac{dE}{d{a}_1}=0$. Solving for a₁, 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}.$

Given a₁, we can then estimate φ₁ from the equation ${\mathrm{\phi }}_1=\frac{a_1}{{\mathrm{\tau }}_2}.$

5. Given τ₁, φ₁, and the fact that ${\mathrm{\tau }}_1=\frac{1}{{\mathrm{\mu }}_1+{\mathrm{\phi }}_1},$ we can obtain an estimate for μ₁.

6. We are now left with three remaining parameters: φ₂, μ₂, and λ₃. Given the experimental constraints on τ₁, τ₂, and τ₃, we assume that the average time to proliferate, 1/λ₃, cannot be <7 days. Therefore, we consider λ₃ to be constrained in the interval $\left[1∕7,{\mathrm{\tau }}_3^{-1}\right]$, with time measured in days. We sample equally spaced values for λ₃, and for each value, we compute ${\mathrm{\phi }}_2=\frac{n_3^{*}\left({\mathrm{\phi }}_3+{\mathrm{\mu }}_3-{\mathrm{\lambda }}_3\right)}{n_2^{*}}$. The ratio $a_2=\frac{n_3^{*}}{n_2^{*}}$ is computed using the linear regression method described above (see Figure 3). In this way, we obtain an estimate for φ₃. We note that the p-values for a₁ and a₂ are given by 7.57 × 10⁻³ and 6.85 × 10⁻³, respectively (both smaller than the significance level α = 0.05).

7. Given τ₂, φ₂, and the fact that ${\mathrm{\tau }}_2=\frac{1}{{\mathrm{\mu }}_2+{\mathrm{\phi }}_2}$, we can obtain an estimate for μ₂.

8. From steps 6 and 7 above, we have generated (a table of) values for φ₂ and μ₂, given a fiducial value for λ₃ in the interval $\left[1∕7,{\mathrm{\tau }}_3^{-1}\right]$. 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^{*}\pm {\mathrm{\sigma }}_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.

FIGURE 3

Figure 3. Linear regression plots for the first model.

We obtain the following parameter values:

$\begin{array}{llll}\hfill & \mathrm{\varphi }=35.350\times {10}^6\mathrm{cells}∕\mathrm{day},{\mathrm{\mu }}_1=0.263{\mathrm{day}}^{-1},& \hfill & \\ \hfill & {\mathrm{\mu }}_3=0.099{\mathrm{day}}^{-1},{\mathrm{\phi }}_1=0.137{\mathrm{day}}^{-1},{\mathrm{\phi }}_3=0.151{\mathrm{day}}^{-1},& \hfill & \end{array}$

and

$\begin{array}{llll}\hfill & {\mathrm{\mu }}_2\in \left[1.295,1.443\right]{\mathrm{day}}^{-1},{\mathrm{\phi }}_2\in \left[0.050,0.198\right]{\mathrm{day}}^{-1},& \hfill & \\ \hfill & {\mathrm{\lambda }}_3\in \left[0.143,0.250\right]{\mathrm{day}}^{-1}.& \hfill & \end{array}$

These parameters imply the following thymic selection rates:

3.1.1. Death rates

9.7 × 10⁵ cells/h die by neglect in compartment 1 (${\mathrm{\mu }}_1{n}_1^{*}$), 4.8 × 10⁵ cells/h die by negative selection in compartment 2 (${\mathrm{\mu }}_2{n}_2^{*}$), and 6.9 × 10⁴ cells/h die by negative selection in compartment 3 (${\mathrm{\mu }}_3{n}_3^{*}$).

3.1.2. Differentiation rates

5.0 × 10⁵ cells/h are positively selected in compartment 1, that is, become post-DP from pre-DP (${\mathrm{\phi }}_1{n}_1^{*}$), 4.4 × 10⁴ cells/h are positively selected in compartment 2, that is, become SP from post-DP (${\mathrm{\phi }}_2{n}_2^{*}$), and 1.0 × 10⁵ cells/h leave compartment 3 to go to the periphery (${\mathrm{\phi }}_3{n}_3^{*}$).

3.1.3. Proliferation rate

1.3 × 10⁵ cells/h proliferate in compartment 3 (${\mathrm{\lambda }}_3{n}_3^{*}$).

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

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

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{array}{llll}\hfill & p_1=\frac{{\mathrm{\mu }}_1}{{\mathrm{\mu }}_1+{\mathrm{\phi }}_1}=65.8\%,p_2=\frac{{\mathrm{\mu }}_2}{{\mathrm{\mu }}_2+{\mathrm{\phi }}_2}=91.7\%,& \hfill & \\ \hfill & p_3=\frac{{\mathrm{\mu }}_3}{{\mathrm{\mu }}_3+{\mathrm{\phi }}_3+{\mathrm{\lambda }}_3}=22.9\%,q_3=\frac{{\mathrm{\lambda }}_3}{{\mathrm{\mu }}_3+{\mathrm{\phi }}_3+{\mathrm{\lambda }}_3}=42.2\%.& \hfill & \end{array}$

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{array}{llll}\hfill & {\mathrm{\tau }}_i=\frac{1}{{\mathrm{\phi }}_i+{\mathrm{\mu }}_i},\mathrm{for}i\in \left\{1,2\right\},\mathrm{and}& \hfill & \\ \hfill & {\mathrm{\tau }}_i=\frac{1}{{\mathrm{\xi }}_i+{\mathrm{\mu }}_i},\mathrm{for}i\in \left\{4,8\right\}.& \hfill & \end{array}$

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¹μ₄ = 0.04 day⁻¹ and μ₈ = 0.11 day⁻¹. We also make use of the estimates derived in Section 3.1 for ϕ, μ₁, μ₂, φ₁, and φ₂. Finally, the average residence times in each compartment, as described in Section 3.1, are given by:

$\begin{array}{llll}\hfill & {\mathrm{\tau }}_1=60\mathrm{h}=2.5\mathrm{days},{\mathrm{\tau }}_2=16\mathrm{h}=0.67\mathrm{days},& \hfill & \\ \hfill & {\mathrm{\tau }}_4=96\mathrm{h}=4\mathrm{days},{\mathrm{\tau }}_8=96\mathrm{h}=4\mathrm{days}.& \hfill & \end{array}$

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

1. Given τ₄, μ₄, and the fact that ${\mathrm{\tau }}_4=\frac{1}{{\mathrm{\mu }}_4+{\mathrm{\xi }}_4},$ we can obtain an estimate for ξ₄.

2. In the same way, given τ₈, μ₈, and the fact that ${\mathrm{\tau }}_8=\frac{1}{{\mathrm{\mu }}_8+{\mathrm{\xi }}_8},$ we can obtain an estimate for ξ₈.

3. We are now left with four remaining parameters: φ₄, φ₈, λ₄, and λ₈. We know that φ₂ = φ₄ + φ₈. We sample φ₄ in the interval [0, φ₂], where φ₂ is the mean value of the interval obtained in Section 2.1, and for each fiducial value for φ₄, we compute the corresponding value for φ₈.

4. Given τ₄, φ₄, and the fact that $n_4^{*}=\frac{n_2^{*}{\mathrm{\phi }}_4}{{\mathrm{\tau }}_4^{-1}-{\mathrm{\lambda }}_4},$ we can compute the fraction $a_3=\frac{n_2^{*}}{n_4^{*}}$ by linear regression (see Figure 4), and thus obtain an estimate for λ₄. Note that we will reject values of λ₄ 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 λ₃, and we know that ${\mathrm{\lambda }}_4{n}_4^{*}+{\mathrm{\lambda }}_8{n}_8^{*}={\mathrm{\lambda }}_3{n}_3^{*}.$ As before, we can compute the fractions $a_4=\frac{n_4^{*}}{n_8^{*}}$ and $a_5=\frac{n_3^{*}}{n_8^{*}}$ by linear regression (see Figure 4), and thus obtain an estimate for λ₈. Note that we will reject values of λ₈ that imply the proliferation time is larger than 7 days (see Section 3.1). We note that the p-values for a₃, a₄, and a₅ are given by 8.43 × 10⁻³, 3.33 × 10⁻⁷, and 4.56 × 10⁻⁸, respectively (smaller than the significance level).

6. From steps 3, 4, and 5 above, we have generated (a table of) values for φ₈, λ₄, and λ₈, given a fiducial value for φ₄ in the interval [0, φ₂]. 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).

FIGURE 4

Figure 4. Linear regression plots for the second model.

We obtain the following parameter values:

$\begin{array}{llll}\hfill & {\mathrm{\mu }}_4=0.04{\mathrm{day}}^{-1},{\mathrm{\mu }}_8=0.11{\mathrm{day}}^{-1},& \hfill & \\ \hfill & {\mathrm{\xi }}_4=0.21{\mathrm{day}}^{-1},{\mathrm{\xi }}_8=0.14{\mathrm{day}}^{-1},& \hfill & \end{array}$

and

$\begin{array}{llll}\hfill & {\mathrm{\phi }}_4=0.070{\mathrm{day}}^{-1},{\mathrm{\phi }}_8=0.054{\mathrm{day}}^{-1},& \hfill & \\ \hfill & {\mathrm{\lambda }}_4=0.216{\mathrm{day}}^{-1},{\mathrm{\lambda }}_8=0.093{\mathrm{day}}^{-1}.& \hfill & \end{array}$

These parameters imply the following thymic selection rates:

3.2.1. Death rates

2.05 × 10⁴ cells/h die by negative selection in compartment 4 (${\mathrm{\mu }}_4{n}_4^{*}$) and 1.95 × 10⁴ cells/h die by negative selection in compartment 8 (${\mathrm{\mu }}_8{n}_8^{*}$).

3.2.2. Differentiation rates

2.50 × 10⁴ cells/h are CD4 positively selected in compartment 2, that is, become CD4 SP from post-DP (${\mathrm{\phi }}_4{n}_2^{*}$), 1.90 × 10⁴ cells/h are CD8 positively selected in compartment 2, that is, become CD8 SP from post-DP (${\mathrm{\phi }}_8{n}_2^{*}$), 1.08 × 10⁵ cells/h leave compartment 4 to go to the periphery (${\mathrm{\xi }}_4{n}_4^{*}$), and 2.48 × 10⁴ cells/h leave compartment 8 to go to the periphery (${\mathrm{\xi }}_8{n}_8^{*}$).

3.2.3. Proliferation rates

11.10 × 10⁴ cells/h proliferate in compartment 4 (${\mathrm{\lambda }}_4{n}_4^{*}$) and 1.60 × 10⁴ cells/h proliferate in compartment 8 (${\mathrm{\lambda }}_8{n}_8^{*}$).

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

$\frac{{\mathrm{\xi }}_4{n}_4^{*}+{\mathrm{\xi }}_8{n}_8^{*}}{\mathrm{\varphi }}=8.96\%.$

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{array}{llll}\hfill & s_4=\frac{{\mathrm{\phi }}_4}{{\mathrm{\mu }}_2+{\mathrm{\phi }}_4+{\mathrm{\phi }}_8}=4.7\%,s_8=\frac{{\mathrm{\phi }}_8}{{\mathrm{\mu }}_2+{\mathrm{\phi }}_4+{\mathrm{\phi }}_8}=3.6\%,& \hfill & \\ \hfill & p_2=\frac{{\mathrm{\mu }}_2}{{\mathrm{\mu }}_2+{\mathrm{\phi }}_4+{\mathrm{\phi }}_8}=91.7\%.& \hfill & \end{array}$

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:

$\begin{array}{llll}\hfill p_4& =\frac{{\mathrm{\mu }}_4}{{\mathrm{\mu }}_4+{\mathrm{\xi }}_4+{\mathrm{\lambda }}_4}=8.6\%,q_4=\frac{{\mathrm{\lambda }}_4}{{\mathrm{\mu }}_4+{\mathrm{\xi }}_4+{\mathrm{\lambda }}_4}=46.3\%,& \hfill & \\ \hfill p_8& =\frac{{\mathrm{\mu }}_8}{{\mathrm{\mu }}_8+{\mathrm{\xi }}_8+{\mathrm{\lambda }}_8}=32.1\%,q_8=\frac{{\mathrm{\lambda }}_8}{{\mathrm{\mu }}_8+{\mathrm{\xi }}_8+{\mathrm{\lambda }}_8}=27.0\%.& \hfill & \end{array}$

These probabilities imply that the probability to exit the thymus as a mature CD4 thymocyte (that has already reached the medulla) is given by $\frac{100-\left(8.6+46.3\right)}{100}$, which is 45.1%, and the probability to exit as a mature CD8 thymocyte (that has already reached the medulla) is given by $\frac{100-\left(32.1+27.0\right)}{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:

$\mathrm{\theta }=\left({\mathrm{\tau }}_1,{\mathrm{\tau }}_2,{\mathrm{\tau }}_3,{\mathrm{\varphi }}_{\mathrm{out}},a_1,a_2,{\overline{n}}_3,{\overline{n}}_1\right),$

where a₁, a₂ are the coefficients of the linear regression of $\frac{n_2^{*}}{n_1^{*}}$ and $\frac{n_3^{*}}{n_2^{*}}$, respectively, and ${\overline{n}}_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 ${\mathrm{\theta }}_i+\frac{1}{10}{\mathrm{\theta }}_i$ and ${\mathrm{\theta }}_i-\frac{1}{10}{\mathrm{\theta }}_i$. Consequently, we have 2⁸ 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:

$\mathrm{\theta }=\left({\mathrm{\tau }}_4,{\mathrm{\tau }}_8,{\mathrm{\mu }}_4,{\mathrm{\mu }}_8,a_3,a_4,a_5\right),$

where a₃, a₄, a₅ are the coefficients of the linear regression of $\frac{n_2^{*}}{n_4^{*}}$, $\frac{n_4^{*}}{n_8^{*}}$, and $\frac{n_3^{*}}{n_8^{*}}$, respectively. We have made use of the means of the following parameters of the first model: ϕ, φ₁, μ₁, φ₂, μ₂, λ₃. We perturb each entry of the vector θ as described above. Therefore, we have $2^7{n}_{{\mathrm{\phi }}_4}$ sets of θ, with $n_{{\mathrm{\phi }}_4},$ the number of different values considered for φ₄ in the interval (0, φ₂). 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

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

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:

$\begin{array}{ll}\hfill \mathrm{\Delta }{p}_i& =\frac{1}{{\left({\mathrm{\mu }}_i+{\mathrm{\phi }}_i\right)}^2}\left({\mathrm{\phi }}_i\mathrm{\Delta }{\mathrm{\mu }}_i+{\mathrm{\mu }}_i\mathrm{\Delta }{\mathrm{\phi }}_i\right)\mathrm{for}i=1,2,\end{array}$ (10)

$\begin{array}{ll}\hfill \mathrm{\Delta }{p}_3& =\frac{1}{{\left({\mathrm{\mu }}_3+{\mathrm{\phi }}_3+{\mathrm{\lambda }}_3\right)}^2}\left[\left({\mathrm{\phi }}_3+{\mathrm{\lambda }}_3\right)\mathrm{\Delta }{\mathrm{\mu }}_3+{\mathrm{\mu }}_3\mathrm{\Delta }{\mathrm{\phi }}_3+{\mathrm{\mu }}_3\mathrm{\Delta }{\mathrm{\lambda }}_3\right],\end{array}$ (11)

$\begin{array}{ll}\hfill \mathrm{\Delta }{q}_3& =\frac{1}{{\left({\mathrm{\mu }}_3+{\mathrm{\phi }}_3+{\mathrm{\lambda }}_3\right)}^2}\left[{\mathrm{\lambda }}_3\mathrm{\Delta }{\mathrm{\mu }}_3+{\mathrm{\lambda }}_3\mathrm{\Delta }{\mathrm{\phi }}_3+\left({\mathrm{\mu }}_3+{\mathrm{\phi }}_3\right)\mathrm{\Delta }{\mathrm{\lambda }}_3\right],\end{array}$ (12)

$\begin{array}{l}\Delta s_i=\frac{1}{{({\mu }_2+{\phi }_i+{\phi }_j)}^2}[{\phi }_i\Delta {\mu }_2+{\phi }_i\Delta {\phi }_j+({\mu }_2+{\phi }_j)\Delta {\phi }_i]\\ \text{\hspace{1em}for\hspace{1em}}i=4,j=8\text{\hspace{1em}or\hspace{1em}}i=8,j=4,\end{array}$ (13)

$\begin{array}{l}\Delta p_i=\frac{1}{{({\mu }_i+{\xi }_i+{\lambda }_i)}^2}[{\mu }_i\Delta {\xi }_i+{\mu }_i\Delta {\lambda }_i+({\xi }_i+{\lambda }_i)\Delta {\mu }_i]\\ \text{\hspace{1em}for\hspace{1em}}i=4,8,\end{array}$ (14)

$\begin{array}{l}\Delta q_i=\frac{1}{{({\mu }_i+{\xi }_i+{\lambda }_i)}^2}[{\lambda }_i\Delta {\xi }_i+{\lambda }_i\Delta {\mu }_i+({\xi }_i+{\mu }_i)\Delta {\lambda }_i]\\ \text{\hspace{1em}\hspace{0.17em}\hspace{0.17em}for\hspace{1em}}i=4,8.\end{array}$ (15)

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

TABLE 3

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

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, μ₂, 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 $\frac{{\mathrm{\phi }}_4}{{\mathrm{\phi }}_8}\approx 5:4$, (ii) the CD4:CD8 ratio in the SP pool is given by $\frac{n_4^{*}}{n_8^{*}}\approx 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 $\frac{{\mathrm{\xi }}_4{n}_4^{*}}{{\mathrm{\xi }}_8{n}_8^{*}}\approx 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:

$\mathrm{\sigma }=\frac{{\mathrm{\xi }}_4{n}_4^{*}+{\mathrm{\xi }}_8{n}_8^{*}}{\mathrm{\varphi }}=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 Δσ = 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:

$\left(1-p_1\right)\times \left(1-p_2\right)\times \left(1-p_3\right)=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 10³ 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 (ϕ, φ₁, μ₁, μ₂), 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

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

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 μ₄ and μ₈. 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.

Footnotes

1. ^Private communication from Ben Seddon and Andy Yates.
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.

References

1. Anderson G, Lane P, Jenkinson E. Generating intrathymic microenvironments to establish T-cell tolerance. Nat Rev Immunol (2007) 7(12):954–63. doi: 10.1038/nri2187

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

2. Stritesky GL, Jameson SC, Hogquist KA. Selection of self-reactive T cells in the thymus. Annu Rev Immunol (2012) 30:95. doi:10.1146/annurev-immunol-020711-075035

CrossRef Full Text

3. Palmer E. Negative selection: clearing out the bad apples from the T-cell repertoire. Nat Rev Immunol (2003) 3(5):383–91. doi:10.1038/nri1085

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

4. Jameson S, Hogquist K, Bevan M. Positive selection of thymocytes. Annu Rev Immunol (1995) 13(1):93–126. doi:10.1146/annurev.iy.13.040195.000521

CrossRef Full Text

5. Werlen G, Hausmann B, Naeher D, Palmer E. Signaling life and death in the thymus: timing is everything. Science (2003) 299(5614):1859. doi:10.1126/science.1067833

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

6. Petrie H, Zúñiga-Pflücker J. Zoned out: functional mapping of stromal signaling microenvironments in the thymus. Immunology (2007) 25(1):649. doi:10.1146/annurev.immunol.23.021704.115715

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

7. Zúñiga-Pflücker JC. When three negatives made a positive influence in defining four early steps in T cell development. J Immunol (2012) 189(9):4201–2. doi:10.4049/jimmunol.1202553

CrossRef Full Text

8. Takada K, Ohigashi I, Kasai M, Nakase H, Takahama Y. Development and function of cortical thymic epithelial cells. Curr Top Microbiol Immunol (2013) 373:1–17.

9. Singer A, Adoro S, Park J. Lineage fate and intense debate: myths, models and mechanisms of CD4-versus CD8-lineage choice. Nat Rev Immunol (2008) 8(10):788–801. doi:10.1038/nri2416

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

10. Anderson MS, Su MA. Aire and T cell development. Curr Opin Immunol (2011) 23(2):198–206. doi:10.1016/j.coi.2010.11.007

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

11. Klein L. Dead man walking: how thymocytes scan the medulla. Nat Immunol (2009) 10(8):809–11. doi:10.1038/ni0809-809

CrossRef Full Text

12. Faro J, Velasco S, Gonzalez-Fernandez A, Bandeira A. The impact of thymic antigen diversity on the size of the selected T cell repertoire. J Immunol (2004) 172(4):2247.

Pubmed Abstract | Pubmed Full Text

13. Mehr R, Globerson A, Perelson A. Modeling positive and negative selection and differentiation processes in the thymus. J Theor Biol (1995) 175(1):103–26. doi:10.1006/jtbi.1995.0124

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

14. Thomas-Vaslin V, Altes H, de Boer R, Klatzmann D. Comprehensive assessment and mathematical modeling of T cell population dynamics and homeostasis. J Immunol (2008) 180(4):2240.

Pubmed Abstract | Pubmed Full Text

15. Souza-e Silva H, Savino W, Feijóo R, Vasconcelos A. A cellular automata-based mathematical model for thymocyte development. PLoS One (2009) 4(12):e8233. doi:10.1371/journal.pone.0008233

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

16. Efroni S, Harel D, Cohen I. Emergent dynamics of thymocyte development and lineage determination. PLoS Comput Biol (2007) 3(1):e13. doi:10.1371/journal.pcbi.0030013

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

17. Müller V, Bonhoeffer S. Quantitative constraints on the scope of negative selection. Trends Immunol (2003) 24(3):132–5. doi:10.1016/S1471-4906(03)00028-0

CrossRef Full Text

18. Detours V, Mehr R, Perelson A. A quantitative theory of affinity-driven T cell repertoire selection. J Theor Biol (1999) 200(4):389–403. doi:10.1006/jtbi.1999.1003

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

19. Prasad A, Zikherman J, Das J, Roose J, Weiss A, Chakraborty A. Origin of the sharp boundary that discriminates positive and negative selection of thymocytes. Proc Natl Acad Sci U S A (2009) 106(2):528. doi:10.1073/pnas.0805981105

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

20. Ribeiro RM, Perelson AS. Determining thymic output quantitatively: using models to interpret experimental T-cell receptor excision circle (Trec) data. Immunol Rev (2007) 216(1):21–34.

Pubmed Abstract | Pubmed Full Text

21. Sinclair C, Bains I, Yates AJ, Seddon B. Asymmetric thymocyte death underlies the CD4:CD8 T-cell ratio in the adaptive immune system. Proc Natl Acad Sci U S A (2013) 110(31):E2905–14. doi:10.1073/pnas.1304859110

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

22. Surh CD, Sprent J. T-cell apoptosis detected in situ during positive and negative selection in the thymus. Nature (1994) 372(6501):100–3. doi:10.1038/372100a0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

23. Laufer TM, DeKoning J, Markowitz JS, Lo D, Glimcher LH. Unopposed positive selection and autoreactivity in mice expressing class II MHC only on thymic cortex. Nature (1996) 383(6595):81–5. doi:10.1038/383081a0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

24. van Meerwijk JP, Marguerat S, Lees RK, Germain RN, Fowlkes B, MacDonald HR. Quantitative impact of thymic clonal deletion on the T cell repertoire. J Exp Med (1997) 185(3):377–84. doi:10.1084/jem.185.3.377

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

25. Merkenschlager M, Graf D, Lovatt M, Bommhardt U, Zamoyska R, Fisher AG. How many thymocytes audition for selection? J Exp Med (1997) 186(7):1149–58. doi:10.1084/jem.186.7.1149

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

26. Stritesky GL, Xing Y, Erickson JR, Kalekar LA, Wang X, Mueller DL, et al. Murine thymic selection quantified using a unique method to capture deleted T cells. Proc Natl Acad Sci U S A (2013) 110(12):4679–84. doi:10.1073/pnas.1217532110

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

27. Egerton M, Scollay R, Shortman K. Kinetics of mature T-cell development in the thymus. Proc Natl Acad Sci U S A (1990) 87(7):2579–82. doi:10.1073/pnas.87.7.2579

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

28. Saini M, Sinclair C, Marshall D, Tolaini M, Sakaguchi S, Seddon B. Regulation of Zap70 expression during thymocyte development enables temporal separation of CD4 and CD8 repertoire selection at different signaling thresholds. Sci Signal (2010) 3(114):ra23. doi:10.1126/scisignal.2000702

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

29. McCaughtry TM, Wilken MS, Hogquist KA. Thymic emigration revisited. J Exp Med (2007) 204(11):2513–20. doi:10.1084/jem.20070601

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

30. Scollay RG, Butcher EC, Weissman IL. Thymus cell migration: quantitative aspects of cellular traffic from the thymus to the periphery in mice. Eur J Immunol (1980) 10(3):210–8. doi:10.1002/eji.1830100310

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

31. Kenney JF, Keeping ES. Mathematics of Statistics: Part One. Princeton: van Nostrand (1954).

32. Kenney JF, Keeping ES. Mathematics of Statistics: Part Two. Princeton: van Nostrand (1962).

33. Mehr R, Perelson A, Fridkis-Hareli M, Globerson A. Feedback regulation of T cell development in the thymus. J Theor Biol (1996) 181(2):157–67. doi:10.1006/jtbi.1996.0122

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

34. Scollay R, Godfrey D. Thymic emigration: conveyor belts or lucky dips? Immunol Today (1995) 16(6):268–73. doi:10.1016/0167-5699(95)80179-0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

35. Tough DF, Sprent J. Thymic emigration: conveyor belts or lucky dips? Immunol Today (1995) 16(6):273–4. doi:10.1016/0167-5699(95)80180-4

CrossRef Full Text

36. Sim B-C, Aftahi N, Reilly C, Bogen B, Schwartz RH, Gascoigne NR, et al. Thymic skewing of the CD4/CD8 ratio maps with the T-cell receptor α-chain locus. Curr Biol (1998) 8(12):701–S3. doi:10.1016/S0960-9822(98)70276-3

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

37. Boursalian TE, Golob J, Soper DM, Cooper CJ, Fink PJ. Continued maturation of thymic emigrants in the periphery. Nat Immunol (2004) 5(4):418–25. doi:10.1038/ni1049

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

38. Spits H. Development of αβ T cells in the human thymus. Nat Rev Immunol (2002) 2(10):760–72. doi:10.1038/nri913

CrossRef Full Text

39. Marino JH, Tan C, Taylor AA, Bentley C, Van De Wiele CJ, Ranne R, et al. Differential IL-7 responses in developing human thymocytes. Hum Immunol (2010) 71(4):329–33. doi:10.1016/j.humimm.2010.01.009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

40. Molina-París C, Reynolds J, Lythe G, Coles MC. Mathematical model of naive t cell division and survival il-7 thresholds. Front Immunol (2013) 4:434. doi:10.3389/fimmu.2013.00434

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

41. Allen LJS. An Introduction to Stochastic Processes with Applications to Biology. New Jersey: Pearson Education (2003).

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_{\left(x,y\right)}\left(t\right)=\mathrm{Prob}\left\{X\left(t\right)=x,Y\left(t\right)=y|X\left(0\right)=1,Y\left(0\right)=0\right\},$ (A1)

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

$\begin{array}{ccc}\hfill & p_{\left(w,z\right),\left(x,y\right)}\left(\mathrm{\Delta }t\right)=\text{Prob}\{X(t+\Delta t)=w,\hfill & \hfill \\ \hfill & Y(t+\Delta t)=z|X(t)=x,Y(t)=y\}\hfill \\ \hfill & =\left\{\begin{array}{cc}\mathrm{\lambda }x\mathrm{\Delta }t+o\left(\mathrm{\Delta }t\right),\hfill & \mathrm{if}w=x+1,z=y,\hfill \\ \mathrm{\mu }x\mathrm{\Delta }t+o\left(\mathrm{\Delta }t\right),\hfill & \mathrm{if}w=x-1,z=y,\hfill \\ \mathrm{\xi }x\mathrm{\Delta }t+o\left(\mathrm{\Delta }t\right),\hfill & \mathrm{if}w=x-1,z=y+1,\hfill \\ 1-\left(\mathrm{\lambda }+\mathrm{\mu }+\mathrm{\xi }\right)x\mathrm{\Delta }t+o\left(\mathrm{\Delta }t\right),\hfill & \mathrm{if}w=x,z=y,\hfill \\ o\left(\mathrm{\Delta }t\right),\hfill & \mathrm{if}w,z\mathrm{otherwise}.\hfill \\ \hfill \end{array}\right.\hfill \end{array}$ (A2)

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

$\begin{array}{llll}\hfill \frac{\mathrm{d}{p}_{\left(x,y\right)}\left(t\right)}{\mathrm{d}t}& =\mathrm{\lambda }\left(x-1\right)p_{\left(x-1,y\right)}\left(t\right)+\mathrm{\mu }\left(x+1\right)p_{\left(x+1,y\right)}\left(t\right)& \hfill & \\ \hfill & +\mathrm{\xi }\left(x+1\right)p_{\left(x+1,y-1\right)}\left(t\right)-\left(\mathrm{\lambda }+\mathrm{\mu }+\mathrm{\xi }\right)x{p}_{\left(x,y\right)}\left(t\right).\end{array}$ (A3)

Let m_X(t) be the expected number of cells in the compartment under consideration, and m_Y(t) be the expected number of cells which have left the compartment. Similarly, let m_XX(t) be the expectation of the random variable X(t)², m_XY(t) be the expectation of the random variable X(t)Y (t), and m_YY(t) be the expectation of the random variable Y (t)². Then, making use of the probability generating function technique (41), we derive the time evolution for the first two moments of the system:

$\begin{array}{ll}\hfill \frac{\mathrm{d}{m}_X\left(t\right)}{\mathrm{d}t}& =\left(\mathrm{\lambda }-\mathrm{\mu }-\mathrm{\xi }\right)m_X\left(t\right),\end{array}$ (A4)

$\begin{array}{ll}\hfill \frac{\mathrm{d}{m}_Y\left(t\right)}{\mathrm{d}t}& =\mathrm{\xi }{m}_X\left(t\right),\end{array}$ (A5)

$\begin{array}{ll}\hfill \frac{\mathrm{d}{m}_{XX}\left(t\right)}{\mathrm{d}t}& =2\left(\mathrm{\lambda }-\mathrm{\mu }-\mathrm{\xi }\right)m_{XX}\left(t\right)+\left(\mathrm{\lambda }-\mathrm{\mu }-\mathrm{\xi }\right)m_X\left(t\right),\end{array}$ (A6)

$\begin{array}{ll}\hfill \frac{\mathrm{d}{m}_{XY}\left(t\right)}{\mathrm{d}t}& =\left(\mathrm{\lambda }-\mathrm{\mu }-\mathrm{\xi }\right)m_{XY}\left(t\right)+\mathrm{\xi }\left[m_{XX}\left(t\right)-m_X\left(t\right)\right],\end{array}$ (A7)

$\begin{array}{ll}\hfill \frac{\mathrm{d}{m}_{YY}\left(t\right)}{\mathrm{d}t}& =\mathrm{\xi }\left[2m_{XY}\left(t\right)+m_X\left(t\right)\right].\end{array}$ (A8)

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

$m_X\left(t\right)={\mathrm{e}}^{\left(\mathrm{\lambda }-\mathrm{\mu }-\mathrm{\xi }\right)t}.$ (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\left(t\right)=\frac{\mathrm{\xi }}{\mathrm{\mu }+\mathrm{\xi }-\mathrm{\lambda }}\left[1-{\mathrm{e}}^{\left(\mathrm{\lambda }-\mathrm{\mu }-\mathrm{\xi }\right)t}\right].$ (A10)

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

$\underset{t\to +\infty }{\lim }{m}_Y\left(t\right)=\frac{\mathrm{\xi }}{\mathrm{\mu }+\mathrm{\xi }-\mathrm{\lambda }}.$ (A11)

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

$\begin{array}{llll}\hfill m_{YY}\left(t\right)& =\frac{\mathrm{2\lambda }{\mathrm{\xi }}^2}{{\left(\mathrm{\lambda }-\mathrm{\mu }-\mathrm{\xi }\right)}^3}\left[\frac{1}{2}{\mathrm{e}}^{2\left(\mathrm{\lambda }-\mathrm{\mu }-\mathrm{\xi }\right)t}-{\mathrm{e}}^{\left(\mathrm{\lambda }-\mathrm{\mu }-\mathrm{\xi }\right)t}\right]& \hfill & \\ \hfill & -\frac{\mathrm{4\lambda }{\mathrm{\xi }}^2}{{\left(\mathrm{\lambda }-\mathrm{\mu }-\mathrm{\xi }\right)}^2}\left(t-\frac{1}{\mathrm{\lambda }-\mathrm{\mu }-\mathrm{\xi }}\right){\mathrm{e}}^{\left(\mathrm{\lambda }-\mathrm{\mu }-\mathrm{\xi }\right)t}& \hfill & \\ \hfill & +\frac{\mathrm{\xi }}{\mathrm{\lambda }-\mathrm{\mu }-\mathrm{\xi }}{\mathrm{e}}^{\left(\mathrm{\lambda }-\mathrm{\mu }-\mathrm{\xi }\right)t}-\frac{\mathrm{2\lambda }{\mathrm{\xi }}^2}{{\left(\mathrm{\lambda }-\mathrm{\mu }-\mathrm{\xi }\right)}^3}& \hfill & \\ \hfill & -\frac{\mathrm{\xi }}{\mathrm{\lambda }-\mathrm{\mu }-\mathrm{\xi }}.\end{array}$ (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 → +∞)

$\begin{array}{llll}\hfill {\mathrm{\sigma }}_Y^2& =\underset{t\to \infty }{\lim }\left[m_{YY}\left(t\right)-m_Y{\left(t\right)}^2\right]=\frac{\mathrm{2\lambda }{\mathrm{\xi }}^2}{{\left(\mathrm{\mu }+\mathrm{\xi }-\mathrm{\lambda }\right)}^3}+\frac{\mathrm{\xi }}{\mathrm{\mu }+\mathrm{\xi }-\mathrm{\lambda }}& \hfill & \\ \hfill & -\frac{{\mathrm{\xi }}^2}{{\left(\mathrm{\mu }+\mathrm{\xi }-\mathrm{\lambda }\right)}^2}.\end{array}$ (A13)

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{array}{llll}\hfill \frac{\mathrm{\partial M}}{\mathrm{\partial t}}& ={\mathrm{\mu }}_1\left({\mathrm{e}}^{-{\mathrm{\theta }}_1}-1\right)\frac{\mathrm{\partial M}}{\partial {\mathrm{\theta }}_1}+{\mathrm{\phi }}_1\left({\mathrm{e}}^{-{\mathrm{\theta }}_1}{\mathrm{e}}^{{\mathrm{\theta }}_2}-1\right)\frac{\mathrm{\partial M}}{\partial {\mathrm{\theta }}_1}& \hfill & \\ \hfill & +\left({\mathrm{\mu }}_2+{\mathrm{\phi }}_8\right)\left({\mathrm{e}}^{-{\mathrm{\theta }}_2}-1\right)\frac{\mathrm{\partial M}}{\partial {\mathrm{\theta }}_2}+{\mathrm{\phi }}_4\left({\mathrm{e}}^{-{\mathrm{\theta }}_2}{\mathrm{e}}^{{\mathrm{\theta }}_4}-1\right)\frac{\mathrm{\partial M}}{\partial {\mathrm{\theta }}_2}& \hfill & \\ \hfill & +{\mathrm{\mu }}_4\left({\mathrm{e}}^{-{\mathrm{\theta }}_4}-1\right)\frac{\mathrm{\partial M}}{\partial {\mathrm{\theta }}_4}+{\mathrm{\lambda }}_4\left({\mathrm{e}}^{{\mathrm{\theta }}_4}-1\right)\frac{\mathrm{\partial M}}{\partial {\mathrm{\theta }}_4}& \hfill & \\ \hfill & +{\mathrm{\xi }}_4\left({\mathrm{e}}^{-{\mathrm{\theta }}_4}{\mathrm{e}}^{{\mathrm{\theta }}_8}-1\right)\frac{\mathrm{\partial M}}{\partial {\mathrm{\theta }}_4}.\end{array}$ (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 10³ 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.