Impact Factor 6.429

The 5th most cited journal in Immunology

This article is part of the Research Topic

Immune system modeling and analysis

Original Research ARTICLE

Front. Immunol., 14 February 2014 | https://doi.org/10.3389/fimmu.2014.00019

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

  • 1Department of Applied Mathematics, School of Mathematics, University of Leeds, Leeds, UK
  • 2Department of Laboratory Medicine and Pathology, Center for Immunology, University of Minnesota, Minneapolis, MN, USA

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

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 (2225). In this report, we develop a deterministic mathematical model of T cell development in the thymus. Some of us recently published a report where we used a novel approach (Bim−/−Nur77GFP mice) that allowed us to calculate the number of cells undergoing positive and negative selection (26). Using previously published data on the relative life-span of DP and SP cells, we estimated the hourly rate of both positive and negative selection (26). In this manuscript, we make use of (i) a subset of this experimental data, and (ii) the asymmetric death rates observed for CD4 and CD8 precursor thymocytes (21), to develop two mathematical models that will enable us to estimate selection rates in the cortex and the medulla, and provide a quantitative measure for the stringency of thymic selection. The first model (see Section 2.1) allows the identification of the following parameters: DN thymocyte influx into the cortex, pre-DP and post-DP death rates, and pre-DP and post-DP differentiation rates. Under the assumption of asymmetric death rates for the CD4 and CD8 SP thymocytes (21), we extend the first model to provide estimates for the following medullary rates (see Section 2.2): CD4 and CD8 SP death, proliferation, and maturation rates.

2. Materials and Methods

2.1. A First Model of Thymic Development after the DN Stage

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

This mathematical model makes use of a data set obtained from the analysis of eight C57BL/6 wild type and Bim deficient mice (average age 9 weeks), that express a Nur77GFP transgene to indicate TCR signal strength experience (26). Flow cytometric analysis, as described in that study, used standard markers to define various stages of T cell development in the thymus. The Nur77GFP reporter and Bim deficiency were novel modifications that allowed us to quantify cells that normally would be deleted by strong TCR signaling. In the mathematical model, we consider the following thymocyte populations: 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):

ϕn1 – flux of DN thymocytes into compartment n1,

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

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

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

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

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

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

n3μ3 – apoptosis of SP (n3) thymocytes due to strong TCR signal.

FIGURE 1
www.frontiersin.org

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

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:

dn1dt=ϕφ1n1μ1n1,dn2dt=φ1n1φ2n2μ2n2,dn3dt=φ2n2φ3n3μ3n3+λ3n3.

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:

n1*=ϕφ1+μ1,n2*=n1*φ1φ2+μ2,n3*=n2*φ2φ3+μ3λ3.

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

A=(φ1+μ1)00φ1(φ2+μ2)00φ2(φ3+μ3λ3).

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

β1=(φ1+μ1),β2=(φ2+μ2),β3=(φ3+μ3λ3).

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:

n1(t)=n1*+n1(0)e(φ1+μ1)t,n2(t)=n2*+n1(0)φ1[(φ2+μ2)(φ1+μ1)]e(φ1+μ1)t+n2(0)e(φ2+μ2)t,n3(t)=n3*+n1(0)φ1[(φ2+μ2)(φ1+μ1)]×φ2[(φ3+μ3λ3)(φ1+μ1)]e(φ1+μ1)t+n2(0)φ3[(φ3+μ3λ3)(φ2+μ2)]e(φ2+μ2)t+n3(0)e(φ3+μ3λ3)t,

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 n1(t)n1*,n2(t)n2* and n3(t)n3*, 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: 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):

ϕn1 – flux of DN thymocytes into compartment n1,

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

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

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

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

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

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

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

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

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

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

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

FIGURE 2
www.frontiersin.org

Figure 2. Thymic development as hypothesized in the second model. We assume there is a flux, ϕ, that represents the differentiation of DNs into pre-DPs. Pre-DP thymocytes have two fates: further differentiation into the post-DP pool (φ1) or death by neglect (μ1). Post-DP thymocytes have three fates: further differentiation into the SP pool as CD4 SPs (φ4), or CD8 SPs (φ8), or death by apoptosis (μ2). Finally, CD4 or CD8 SP thymocytes have three fates: maturation and exit into the periphery (ξ4) or (ξ8), death by apoptosis (μ4) or (μ8), or proliferation (λ4) or (λ8).

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:

φ2=φ4+φ8,ξ4n4+ξ8n8=φ3n3,μ4n4+μ8n8=μ3n3,λ4n4+λ8n8=λ3n3.

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

dn1dt=ϕφ1n1μ1n1,dn2dt=φ1n1(φ4+φ8)n2μ2n2,dn4dt=φ4n2ξ4n4μ4n4+λ4n4,dn8dt=φ8n2ξ8n8μ8n8+λ8n8.

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:

n1*=ϕφ1+μ1,n2*=n1*φ1φ4+φ8+μ2,n4*=n2*φ4ξ4+μ4λ4,n8*=n2*φ8ξ8+μ8λ8.

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

B=(φ1+μ1)000φ1(φ4+φ8+μ2)000φ4(ξ4+μ4λ4)00φ80(ξ8+μ8λ8).

B is also the Jacobian at the steady state n*=(n1*,n2*,n4*,n8*), as the system of ODEs is linear. The four eigenvalues of B are given by:

β1=(φ1+μ1),β2=(φ4+φ8+μ2),β3=(ξ4+μ4λ4),β4=(ξ8+μ8λ8).

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:

n1(t)=n1*+n1(0)e(φ1+μ1)t,n2(t)=n2*+n1(0)φ1[(φ4+φ8+μ2)(φ1+μ1)]e(φ1+μ1)t+n2(0)e(φ4+φ8+μ2)t,n4(t)=n4*+n1(0)φ1[(φ4+φ8+μ2)(φ1+μ1)]×φ4+φ8[(ξ4+μ4λ4)(φ1+μ1)]e(φ1+μ1)t+n2(0)ξ4[(ξ4+μ4λ4)(φ4+φ8+μ2)]e(φ4+φ8+μ2)t+n4(0)e(ξ4+μ4λ4)t,n8(t)=n8*+n1(0)φ1[(φ4+φ8+μ2)(φ1+μ1)]×φ4+φ8[(ξ8+μ8λ8)(φ1+μ1)]e(φ1+μ1)t+n2(0)ξ8[(ξ8+μ8λ8)(φ4+φ8+μ2)]e(φ4+φ8+μ2)t+n8(0)e(ξ8+μ8λ8)t,

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 n1(t)n1*,n2(t)n2*, n4(t)n4*, and n8(t)n8*, 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: n1*,n2*,n3*. 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:

τi=1φi+μi,fori{1,2,3}.

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

TABLE 1
www.frontiersin.org

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 (2729)

τ1=60h=2.5days,τ2=16h=0.67days,τ3=96h=4days.

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

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

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

4. We also have n2*=φ1τ2n1*, 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, a1=n2*n1*, and make use of the experimental data to write: n2*,i=an1*,i+εi, for i = 1, 2, …, 8. Thus, the squared error is given by:

E(a1)=i=18(n2*,ia1n1*,i)2.

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

a1=i=18n1*,in2*,ii=18(n1*,i)2.

Given a1, we can then estimate φ1 from the equation φ1=a1τ2.

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

6. We are now left with three remaining parameters: φ2, μ2, and λ3. Given the experimental constraints on τ1, τ2, and τ3, we assume that the average time to proliferate, 1/λ3, cannot be <7 days. Therefore, we consider λ3 to be constrained in the interval [17,τ31], with time measured in days. We sample equally spaced values for λ3, and for each value, we compute φ2=n3*(φ3+μ3λ3)n2*. The ratio a2=n3*n2* 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).

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

8. From steps 6 and 7 above, we have generated (a table of) values for φ2 and μ2, given a fiducial value for λ3 in the interval [17,τ31]. 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 ni(t = 21 days) belongs to the interval ni*±σi, with i = 1, 2, 3, and where ni* 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
www.frontiersin.org

Figure 3. Linear regression plots for the first model.

We obtain the following parameter values:

ϕ=35.350×106cellsday,μ1=0.263day1,μ3=0.099day1,φ1=0.137day1,φ3=0.151day1,

and

μ2[1.295,1.443]day1,φ2[0.050,0.198]day1,λ3[0.143,0.250]day1.

These parameters imply the following thymic selection rates:

3.1.1. Death rates

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

3.1.2. Differentiation rates

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

3.1.3. Proliferation rate

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

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

φ3n3*ϕ=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:

p1=μ1μ1+φ1=65.8%,p2=μ2μ2+φ2=91.7%,p3=μ3μ3+φ3+λ3=22.9%,q3=λ3μ3+φ3+λ3=42.2%.

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: n1*,n2*,n4*,n8*. 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:

τi=1φi+μi,fori{1,2},andτi=1ξi+μi,fori{4,8}.

Recent experimental data provide support for asymmetric death rates in the CD4 and CD8 SP compartments (21). The estimated death rates for CD4 and CD8 SP thymocytes are1μ4 = 0.04 day−1 and μ8 = 0.11 day−1. We also make use of the estimates derived in Section 3.1 for ϕ, μ1, μ2, φ1, and φ2. Finally, the average residence times in each compartment, as described in Section 3.1, are given by:

τ1=60h=2.5days,τ2=16h=0.67days,τ4=96h=4days,τ8=96h=4days.

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

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

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

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

4. Given τ4, φ4, and the fact that n4*=n2*φ4τ41λ4, we can compute the fraction a3=n2*n4* by linear regression (see Figure 4), and thus obtain an estimate for λ4. Note that we will reject values of λ4 that imply the proliferation time is larger than 7 days (see Section 3.1).

5. In Section 3.1, we obtained an estimate for the mean of λ3, and we know that λ4n4*+λ8n8*=λ3n3*. As before, we can compute the fractions a4=n4*n8* and a5=n3*n8* 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).

6. From steps 3, 4, and 5 above, we have generated (a table of) values for φ8, λ4, and λ8, given a fiducial value for φ4 in the interval [0, φ2]. As discussed in Section 3.1, we only accept parameter values which provide thymocyte cell counts at time t = 21 days that are within one standard deviation from the experimentally determined values (see Table 1).

FIGURE 4
www.frontiersin.org

Figure 4. Linear regression plots for the second model.

We obtain the following parameter values:

μ4=0.04day1,μ8=0.11day1,ξ4=0.21day1,ξ8=0.14day1,

and

φ4=0.070day1,φ8=0.054day1,λ4=0.216day1,λ8=0.093day1.

These parameters imply the following thymic selection rates:

3.2.1. Death rates

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

3.2.2. Differentiation rates

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

3.2.3. Proliferation rates

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

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

ξ4n4*+ξ8n8*ϕ=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:

s4=φ4μ2+φ4+φ8=4.7%,s8=φ8μ2+φ4+φ8=3.6%,p2=μ2μ2+φ4+φ8=91.7%.

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:

p4=μ4μ4+ξ4+λ4=8.6%,q4=λ4μ4+ξ4+λ4=46.3%,p8=μ8μ8+ξ8+λ8=32.1%,q8=λ8μ8+ξ8+λ8=27.0%.

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

3.3. Sensitivity Analysis

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

θ=(τ1,τ2,τ3,ϕout,a1,a2,n¯3,n¯1),

where a1, a2 are the coefficients of the linear regression of n2*n1* and n3*n2*, respectively, and n¯i is the experimental mean value of ni.

We perturb each entry of the vector θ by adding and subtracting 10% of its value. Therefore, we now have two values for θi, equal to θi+110θi and θi110θi. 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:

θ=(τ4,τ8,μ4,μ8,a3,a4,a5),

where a3, a4, a5 are the coefficients of the linear regression of n2*n4*, n4*n8*, and n3*n8*, 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 27nφ4 sets of θ, with nφ4, the number of different values considered for φ4 in the interval (0, φ2). Parameter values will only be accepted if they provide a stable solution before t = 21 days.

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

TABLE 2
www.frontiersin.org

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:

Δpi=1(μi+φi)2(φiΔμi+μiΔφi)fori=1,2,
Δp3=1(μ3+φ3+λ3)2[(φ3+λ3)Δμ3+μ3Δφ3+μ3Δλ3],
Δq3=1(μ3+φ3+λ3)2[λ3Δμ3+λ3Δφ3+(μ3+φ3)Δλ3],
Δsi=1(μ2+φi+φj)2[φiΔμ2+φiΔφj+(μ2+φj)Δφi]      for i=4,j=8 or i=8,j=4,
Δpi=1(μi+ξi+λi)2[μiΔξi+μiΔλi+(ξi+λi)Δμi]      for i=4,8,
Δqi=1(μi+ξi+λi)2[λiΔξi+λiΔμi+(ξi+μi)Δλi]        for i=4,8.

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

TABLE 3
www.frontiersin.org

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 (2729), murine thymic export rate (14, 26, 30), and recently reported asymmetric death rates for the CD4 SP and CD8 SP thymocytes (21). Our preliminary results support the unexpectedly high death rate in the post-DP thymocyte population observed in Ref. (21). We note that our approach is unrelated to that of Sinclair et al. both experimentally and mathematically (21). This rate, μ2, has been estimated to be at least an order of magnitude larger than any of the other death rates in the pre-DP, CD4 SP, or CD8 SP pools (see Table 2). In terms of selection rates, our analysis yields the following: pre-selection thymocytes (pre-DPs) have a 65.8% probability of dying by neglect in the cortex, and a 34.2% probability of becoming post-selection thymocytes (post-DPs). At the post-selection stage, post-DPs have a 91.7% probability of dying by negative selection (apoptosis) in the cortex, a 4.7% probability of becoming CD4 SPs, and a 3.6% probability of becoming CD8 SPs. In the medulla, CD4 SPs have an 8.6% probability of dying by negative selection (apoptosis), whereas CD8 SPs have a 32.1% probability of dying by negative selection. CD4 SPs have a 45.1% probability of exiting the thymus and reaching the periphery as mature thymocytes, whereas that probability for CD8 SPs is only 40.9%. Finally, the data supports some level of cellular proliferation in the medulla, with CD4 SPs having a 46.3% probability of proliferation and CD8 SPs a 27% probability.

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

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

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

σ=ξ4n4*+ξ8n8*ϕ=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:

(1p1)×(1p2)×(1p3)=2.19%.

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

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

FIGURE 5
www.frontiersin.org

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.

Acknowledgments

We thank Ed Palmer, Robin Callard, Andy Yates, and Ben Seddon for helpful discussions, and Andy Yates and Ben Seddon for allowing us to make use of their mathematical estimates for μ4 and μ8. Grant Lythe and Carmen Molina-París thank Bill Hlavacek, Ruy Ribeiro, and Alan Perelson (Los Alamos National Laboratory) for their hospitality. This work was supported by a BBSRC Research Development Fellowship BB/G023395/1 (to Carmen Molina-París), an FP7 IRSES INDOEUROPEAN-MATHDS Network PIRSES-GA-2012-317893 (to Grant Lythe and Carmen Molina-París), and National Institutes of Health Grants R37 AI39560 (to Kristin A. Hogquist) and F32 AI100346 (to Gretta L. Stritesky). Maria Sawicka, Joseph Reynolds, Niloufar Abourashchi, Grant Lythe, and Carmen Molina-París acknowledge the hospitality of the Max Planck Institute for the Physics of Complex Systems, Dresden, and the International Centre for Mathematical Sciences, Edinburgh, where part of this work was discussed and presented.

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(x,y)(t)=Prob{X(t)=x,Y(t)=y|X(0)=1,Y(0)=0},

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

p(w,z),(x,y)(Δt)=Prob{X(t+Δt)=w,Y(t+Δt)=z|X(t)=x,Y(t)=y}=λxΔt+o(Δt),ifw=x+1,z=y,μxΔt+o(Δt),ifw=x1,z=y,ξxΔt+o(Δt),ifw=x1,z=y+1,1(λ+μ+ξ)xΔt+o(Δt),ifw=x,z=y,o(Δt),ifw,zotherwise.

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

dp(x,y)(t)dt=λ(x1)p(x1,y)(t)+μ(x+1)p(x+1,y)(t)+ξ(x+1)p(x+1,y1)(t)(λ+μ+ξ)xp(x,y)(t).

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

dmX(t)dt=(λμξ)mX(t),
dmY(t)dt=ξmX(t),
dmXX(t)dt=2(λμξ)mXX(t)+(λμξ)mX(t),
dmXY(t)dt=(λμξ)mXY(t)+ξ[mXX(t)mX(t)],
dmYY(t)dt=ξ[2mXY(t)+mX(t)].

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

mX(t)=e(λμξ)t.

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

mY(t)=ξμ+ξλ1e(λμξ)t.

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

limt+mY(t)=ξμ+ξλ.

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

mYY(t)=ξ2(λμξ)312e2(λμξ)te(λμξ)tξ2(λμξ)2t1λμξe(λμξ)t+ξλμξe(λμξ)tξ2(λμξ)3ξλμξ.

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 → +∞)

σY2=limtmYY(t)mY(t)2=ξ2(μ+ξλ)3+ξμ+ξλξ2(μ+ξλ)2.

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

∂M∂t=μ1(eθ11)∂Mθ1+φ1(eθ1eθ21)∂Mθ1+(μ2+φ8)(eθ21)∂Mθ2+φ4(eθ2eθ41)∂Mθ2+μ4(eθ41)∂Mθ4+λ4(eθ41)∂Mθ4+ξ4(eθ4eθ81)∂Mθ4.

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

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

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

Citation: Sawicka M, Stritesky GL, Reynolds J, Abourashchi N, Lythe G, Molina-París C and Hogquist KA (2014) From pre-DP, post-DP, SP4, and SP8 Thymocyte Cell Counts to a Dynamical Model of Cortical and Medullary Selection. Front. Immunol. 5:19. doi: 10.3389/fimmu.2014.00019

Received: 11 September 2013; Accepted: 15 January 2014;
Published online: 14 February 2014.

Edited by:

Ramit Mehr, Bar-Ilan University, Israel

Reviewed by:

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

Copyright: © 2014 Sawicka, Stritesky, Reynolds, Abourashchi, Lythe, Molina-París and Hogquist. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Carmen Molina-París, Department of Applied Mathematics, School of Mathematics, University of Leeds, Leeds LS2 9JT, UK e-mail: carmen@maths.leeds.ac.uk;
Kristin A. Hogquist, Department of Laboratory Medicine and Pathology, Center for Immunology, University of Minnesota, 2101 6th Street SE, Minneapolis, MN 55414, USA e-mail: hogqu001@umn.edu

Maria Sawicka and Gretta L. Stritesky have contributed equally to this work.