Stochastic descriptors to study the fate and potential of naive T cell clonotypes in the periphery.

The population of naive T cells in the periphery is best described by determining both its T cell receptor diversity, or number of clonotypes, and the sizes of its clonal subsets. In this paper, we make use of a previously introduced mathematical model of naive T cell homeostasis, to study the fate and potential of naive T cell clonotypes in the periphery. This is achieved by the introduction of several new stochastic descriptors for a given naive T cell clonotype, such as its maximum clonal size, the time to reach this maximum, the number of proliferation events required to reach this maximum, the rate of contraction of the clonotype during its way to extinction, as well as the time to a given number of proliferation events. Our results show that two fates can be identified for the dynamics of the clonotype: extinction in the short-term if the clonotype experiences too hostile a peripheral environment, or establishment in the periphery in the long-term. In this second case the probability mass function for the maximum clonal size is bimodal, with one mode near one and the other mode far away from it. Our model also indicates that the fate of a recent thymic emigrant (RTE) during its journey in the periphery has a clear stochastic component, where the probability of extinction cannot be neglected, even in a friendly but competitive environment. On the other hand, a greater deterministic behaviour can be expected in the potential size of the clonotype seeded by the RTE in the long-term, once it escapes extinction.

Introduction

T cells are the set of lymphocytes characterised by the expression of a specialised receptor, called the T cell receptor (TCR). T cells can thus, be classified in “families” (or clonotypes) according to the molecular structure of the TCR they display on their membrane (Murphy et al. 2008). On average a T cell expresses identical copies of a given TCR molecule. The number of T cells in the adaptive immune system of an adult human tends to a stationary “homeostatic” distribution (McLean et al. 1997; Freitas and Rocha 2000), specified by its size (the total number of T cells) and TCR diversity (the number of different clonotypes or TCR molecular structures). The size and diversity of the naive T cell pool (those T cells that have not taken part in a previous immune response) is essential to recognise a broad range of pathogens (Freitas and Rocha 1999, 2000; Tanchot and Rocha 1998). An adult mouse has a homeostatic population of naive T cells (Mason 1998) and a TCR diversity estimated around  (Zarnitsyna et al. 2013), and a healthy adult human has a homeostatic population of naive T cells (Hazenberg et al. 2003) distributed in different TCR specificities (Zarnitsyna et al. 2013; Johnson et al. 2014).

Naive T cells originate from T cell precursors (or thymocytes) that survive positive and negative selection in the thymus (Stritesky et al. 2012; Moran and Hogquist 2012; Bird 2009). Thus, naive T cells have not been activated by exposure to pathogen-derived peptide fragments or antigens. During their lifetime, naive T cells constantly circulate in the blood to visit lymph nodes, which we refer to as the “periphery” (Takada and Jameson 2009). The population of naive T cells, which have never taken part in an immune response, is maintained by continuous, yet slow, division in the periphery (or homeostatic proliferation) (Tanchot and Rocha 1998; Troy and Shen 2003; Takada and Jameson 2009). Antigen presenting cells provide the stimulatory signals that induce naive T cell division. The signals are delivered via the interaction between TCRs (on the surface of T cells) and the ligands (also referred to as self-peptides or self-antigens) expressed on the surface of antigen presenting cells. Although the number of naive T cells is approximately constant for much of an individual’s life, as one ages both the size and the TCR diversity of the naive pool decline (Goronzy et al. 2007). This loss of T cell receptor diversity is also observed for B cells (Gibson et al. 2009), and is associated with increased susceptibility to infection and poor health. In the long-term, during the ageing process, it is important to study both the probability of extinction and the mean time to extinction for different T cell clonotypes. This clearly requires a stochastic framework. Some of the homeostatic mechanisms which ensure that a maximum of “antigen space coverage” is achieved with the available numerical diversity of T cell clonotypes (Correia-Neves et al. 2001; Mahajan et al. 2005), have been analysed mathematically before by means of deterministic (De Boer and Perelson 1997) and stochastic (Stirk et al. 2008) models. Yet, the question of how the naive T cell population remains diverse is still not completely understood (De Boer et al. 2012).

During the homeostatic period the organisation of the peripheral naive T cell pool requires that an existing T cell must die if a new T cell is generated in the thymus or by peripheral cell division (Tanchot and Rocha 1998). It is then reasonable to pose a second challenge: to quantify the probability that a given TCR specificity, namely, a single recent thymic emigrant (RTE), can be established in the periphery, while competing with the pre-existing population of naive T cells. Recent thymic emigrants are those naive T cells that have just migrated out of the thymus, after surviving positive and negative selection, to become part of the peripheral T cell population (Fink 2013). Given that double positive thymocytes hardly divide in the thymic cortex and that single positive thymocytes divide at most once in the medulla of the thymus (Sinclair et al. 2013; Stritesky et al. 2013; Sawicka et al. 2014; Yates 2014), one expects a given TCR specificity to be introduced in the periphery by a single recent thymic emigrant. There is evidence to support the fact that RTEs do not have an intrinsic lifespan (Tanchot and Rocha 1998). Given the requirement of naive T cells for continuous TCR ligation to survive in the periphery, it is natural to hypothesise that the TCR cross-reactivity of a given RTE (how many self-peptides can provide a homeostatic stimulus to the TCR under consideration) (Sewell 2012), the number of different clonotypes (competitor clonotypes) that can bind those self-peptides, and the clonal size of each competitor, are key parameters that will determine the fate and potential in the periphery of a recent thymic emigrant.

The aim of this paper is to make use of a previously introduced mathematical model of naive T cell homeostasis (Stirk et al. 2008, 2010), and define new stochastic descriptors that allow us to provide answers to the previous questions. In Refs. Stirk et al. (2008, (2010), a continuous time birth and death Markov model that describes the time evolution of a peripheral naive T cell clonotype was introduced. In this reference, it was shown that extinction takes place with certainty for all parameter values of the model, and the average time to extinction was computed. However, the question of how the affinity and the cross-reactivity of a given T cell clonotype (or RTE) determines (1) its ability to become part of the peripheral naive T cell repertoire, (2) its maximum size if the clone is to get established in the naive peripheral pool, and (3) the timescale to get established in the periphery at its maximum size, was not considered. We study these questions by defining novel stochastic descriptors for a given TCR clonotype or RTE.

The paper is organised as follows: in Sect. 2 we introduce the stochastic birth and death process and define a number of stochastic descriptors (continuous and discrete) to characterise the maximum clonal size (Sect. 2.1) and the time to reach that size, the number of divisions to reach the maximum clonal size (Sect. 2.2), and the time to get to a given number of divisions (Sect. 2.3). Section 3 provides numerical results for the previously introduced descriptors, exploring the three parameter regimes described in Sect. 2. We conclude with a discussion of our results both from a mathematical and immunological perspective in Sect. 4.

Mathematical model

We make use of the mathematical model introduced in Stirk et al. (2008) to describe the time evolution of the number of cells in a given T cell clonotype, that is, the number of naive T cells that express identical T cell receptor molecules. The main assumption of the model is that T cells of a given clonotype compete for signals delivered by molecules (self-peptides also referred to as pMHC) expressed on antigen presenting cells (APCs). These signals induce one round of cell division, that is, a birth event. All cells of the given clonotype can also die and this is a death event in the model.

The underlying mathematical model is a birth and death process on the state space , where the random variable X(t) describes the number of cells of a given T cell clonotype as a function of time t. Its birth and death rates (Fig. 1) are specified for by Stirk et al. (2008):Parameters , , , and were introduced in  Stirk et al. (2008), and they have the following meaning:It was shown in  Stirk et al. (2008) that, for any value of the parameters, is an absorbing state of the stochastic process, the time to extinction (or absorption) from any other state is finite with probability one, and its mean is also finite.

Fig. 1

Birth and death process (with absorption) for the time evolution of the number of cells in a given naive T cell clonotype

is the per (naive) cell death rate for cells of the clonotype under consideration,

is the per (naive) cell rate of homeostatic proliferation due to signals from the self-peptides that can bind to the TCR of the clonotype under consideration,

is the number of other clonotypes (different from the one under consideration) that can compete for the same homeostatic proliferation signals, and

is the characteristic size of T cell clonotypes that compete for homeostatic proliferation signals with the clonotype of interest. We note that these competing clonotypes are not explicitly modelled and by assuming that they all have a characteristic number of naive T cells given by , the time evolution of the clonotype of interest can be reduced to a univariate Markov process (details about this approximation can be found in  Stirk et al. (2008)).

As introduced above, parameter is a measure both of the affinity and the cross-reactivity of the T cell receptor for self-pMHC molecules (Sewell 2012), in the sense that the rate of homeostatic proliferation depends not only on how well a TCR can bind a given pMHC complex, but how many different pMHC complexes can provide survival signals to the T cell clonotype under consideration, characterised by its TCR molecule. Parameter is the number of competitors of the clonotype under consideration (Stirk et al. (2008)). Three regimes can be identified in parameter space that describe different immunological scenarios: the first one is the limit of no inter-clonal competition (), the second is the limit of large inter-clonal competition (), and the third one is the intermediate regime of competition ().

Hard niche clonotype In the special case , the population under study (the number of T cells that belong to a given TCR clonotype) does not compete with any other clonotypes, and thus the previous expressions for the birth and death rates simplify to (Stirk et al. (2008)) Soft niche clonotype The opposite limit, , represents a highly competitive environment, where the population of T cells under study competes with a large number of different clonotypes, and thus the expressions for the birth and death rates become (Stirk et al. (2008)) Intermediate niche clonotype In the case , the TCR clonotype will be referred to as an intermediate niche clonotype, and its birth and death rates will be given by Eq. (1).

These three regimes, hard, soft and intermediate niche clonotypes, will be explored in Sect. 3. Under these regimes, we introduce in this Section several stochastic descriptors which will allow us to study the dynamics of the clonotype under consideration. In particular, in Sect. 2.1 our interest is in the maximum clonal size reached by the clonotype and the time to reach this maximum size, for which we obtain analytical expressions for its probability mass function and its different order moments, respectively. We analyse the number of proliferation events to reach the maximum clonal size in Sect. 2.2, and in Sect. 2.3 we focus on the time to reach a given number of proliferation (or division) events, as a measure of the proliferative capacity (or potential) of the clonotype under study. Finally, the rate of contraction of the clonotype during its way to extinction can be studied by means of the time to contraction to a given clonal size, which is also analysed in Sect. 2.1.

Maximum clonal size and time to reach it

Our interest in this section is in the maximum number of cells belonging to the clonotype under consideration and the time needed to reach this maximum. To begin with, we define the random variablesfor , which amount to the maximum clonal size attained by the clonotype, under the assumption that the initial size equals i, and the time to reach this maximum clonal size, respectively; note that and , since 0 is an absorbing state in . In order to study both variables in parallel, we define the time to reach a total clonal size , given that the initial clonal size equals i, as the auxiliary random variablefor , with . We point out here that these auxiliary random variables have their own immunological interest, since attaining a given clonal size might indicate that the clonotype has become large enough to mount an immune response in a timely fashion, or it has grown too large and might lead to autoimmunity, so that a tolerant immune state has been lost. The analysis of also allows the study of the random variables of interest and .

For a particular clonal size , the random variable might be analysed in a different manner depending on whether (time to contraction to a given clonal size ), or (time to expansion to a given clonal size ). We note that in the latter case, reaching is not certain, so that is a defective random variable, while reaching is certain in the former case, since absorption at state 0 occurs with probability one (Stirk et al. 2008).

Time to contraction to a given clonal size , given the initial clonal size i

For an initial clonal size i with , the random variable amounts to the time until absorption into for a birth and death process defined on a single absorbing state and the class of transient states, , with birth rates and death rates . Since absorption of the underlying process occurs in a finite time almost surely (see  Stirk et al. (2008)), an appeal to  Karlin and McGregor (1957) allows us to describe the expected values of from the equalitywhere and , for ; note that the expected values corresponding to are related to the total extinction of the clonotype. We refer the reader to  Artalejo et al. (2012), where an analogous descriptor is analysed for the spread of an SIS epidemic among a population consisting of N individuals. In this case, the series becomes the finite sum and can be seen as a phase-type random variable (see e.g.,  [Kulkarni (1996), Sect. 6.7] and  [Latouche and Ramaswami (1999), Chapter 2]).

Time to expansion to a given clonal size , given the initial clonal size i, and analysis of and

For a fixed (and given) clonal size , we introduce the following notation:where is a random variable that takes the values 1 if , and 0 otherwise. The quantities defined in Eq. (2) will allow us to analyse the random variables and , but the distribution of is defective, since if the clonotype becomes extinct before reaching the size ; note that is the probability of such an extinction event, aswith in the case .

By a first-step argument, the restricted Laplace-Stieltjes transforms satisfy the following set of linear equations:These equations can be rewritten in multiplicative form aswithThen, by using forward elimination, we may obtainwithIn terms of the functions , we may rewrite assince satisfieswith . This implies thatfrom which it follows that for clonal sizes , since .

In evaluating the restricted moments , we first derive the probabilities as the values from Eq. (4). In particular, it is seen that and we can writewhere and , for . Then, we haveMoments , for clonal sizes , are derived from Eq. (3). In fact, by taking derivatives in Eq. (3) we obtainIn order to solve Eq. (5), for a fixed value , we introduce some notation as follows:Equation (5) is then equivalent to , for , which implies thatwhere . As a result, for , we haveWe also have since , so thatEquation (6) is a recursive procedure that allows us to compute the kth order moments, , in terms of previously computed moments, , since and . In the special cases and 1, it is readily seen thatand the asymptotic value is always finite, with for initial clonal sizes , since , for any value of the parameter .

We finally focus on the random variables and . Note that the distribution of is readily derived from the values , since , and the kth order moment, , of the time to reach the maximum clonal size can be evaluated asFor practical or computational purposes, the above series should be replaced by the finite sumwhere can be selected as the (100q)th percentile of for a probability close enough to 1. Although an analytical study of does not seem to be feasible, we note that by using the finite sum in Eq. (7) we ensure that the probability mass accumulated by is greater than q, whence the probability can be interpreted as a global error measure. Furthermore, this type of truncation procedure has been efficiently used for epidemics (Almaraz et al. 2016) and competition processes (Gómez-Corral and López García 2011, 2012).

Number of proliferation events to reach the maximum clonal size

We are now interested in the number of one-step transitions of occurring in the random interval , which allows us to record the number of proliferation events to reach the maximum clonal size, for . This descriptor can be seen as a discrete version of the random variable in Sect. 2.1, and its analysis can be carried out by means of the number of proliferation events to reach a total clonal size, which is defined as the number of one-step transitions of to register for the first time.

We first introduce the following notation:for initial clonal sizes . It can be easily verified that , and , for . For clonal sizes , it is seen that , since and, consequently, the probabilities can be computed from Eq. (4), and they amount to the probabilities of reaching the clonal size before extinction, given that the initial clonal size is i.

By taking derivatives on the equalitywe obtainwhich is similar to Eq. (5) with the term replaced by . This implies that, similarly to the solution of Eq. (6), the solution of Eq. (8) has the formSimilarly to Sect. 2.1.2, it can be seen that is always finite, with for initial clonal sizes .

In the case of the random variable , not only every kth order factorial moment can be obtained from Eq. (9), but also its distribution. We can writewith for , if , and 0 if , for and , and for . For , we haveso that a recursion on j and i can be obtained in order to compute the probability mass function of . The probability mass function obtained in this way is consistent with the expressions obtained in Eq. (9) for the factorial moments, but the proof is omitted here.

Finally, the factorial moments of the random variable can be approximated in a similar way as described for Eq. (7). In particular, if we denote , we have

Time to reach a given number of proliferation events

In this section, we fix an initial clonal size and a number of proliferation events, and we study the time to reach a total number D of proliferation events (one-step transitions of ), with . In order to study this descriptor, we consider the augmented process , where D(t) denotes the number of division events up to time t, with . For the initial clonal size , the augmented process starts from state and is defined on the finite state spacewhereStates in are absorbing states. In particular, states in represent the extinction of the clonotype when the number D of proliferation events has not been reached, while states in reflect that the number D of proliferation events has been reached. States in are transient states. Figure 2 illustrates the dynamics of the augmented process in the case .

Fig. 2

Transitions between augmented states in the case

In a more general setting, we define the random variablefor states . Note that the descriptor is equivalent to the random variable , and states in and lead to the values and , respectively. For states in , the random variable is defective, and is the probability of reaching before getting to .

For states , we introduce the notationand we observe thatAs a result, the Laplace-Stieltjes transforms verify the boundary conditions if , and if . For states , we havewhich yields a finite system of linear equations that can be solved in a recursive manner (Algorithm 1).

The boundary conditions for the moments are given bywhere the probabilities , for states , are derived from Algorithm 1 by selecting . If we take derivatives in Eq. (11), we may write downfor states . Eq. (12) can be solved (Algorithm 2) by adapting our arguments in Algorithm 1.

Numerical results

In this Section we carry out a set of numerical experiments in order to analyse the potential of a recent thymic emigrant (of a given clonotype) in the periphery to expand to different clonal sizes, or to proliferate for a given number of divisions, as well as the random times of these events, and the rate of contraction of the clonotype under consideration. We study these random variables under several competition and signalling environments, which are specified by the choice of parameters and . We point out here that parameters and should be seen as intrinsically related to the particular TCR expressed by the T cells within the clonotype under consideration, since this TCR determines the number of different self-peptides that the T cell can interact with. On the other hand, is not directly related to the TCR and should be considered an environmental parameter, since it is the characteristic size of the competing clonotypes. We set the per cell death rate , so that the time unit in the process under study is the mean lifetime of a T cell.

We consider throughout this Section homeostatic signalling rates that have been selected for a wide range of parameter values, in view of some preliminary numerical experiments. In particular, our results regarding the hard, the intermediate and the soft niche cases amount to values of the signalling rate , the number of competitors, and the characteristic size of these competitors ranging from friendly scenarios, where the T cell clonotype under consideration can become established in the periphery, to hostile environments, where this establishment is not possible at all. In particular, we analyse:

In Tables 1, 2 and 3, the mean values , and are computed by considering only clonal sizes up to the 99th percentile of , and by means of the truncating procedures given by Eqs. (7) and (10). However, for parameter values in Table 3, values of , and are computed by considering only clonal sizes up to the percentile of , since in these cases which could lead to misleading results. In this way, the quantities obtained are approximations to the true ones, and we focus on the particular case in which the clonotype under study has been recently released to the periphery (RTE), so that . A direct comparison of our results with results obtained by Gillespie simulations suggests a significant level of precision for the truncation procedures proposed by Eqs. (7) and (10). For example, values , and reported in Table 1 for have Gillespie ( simulations) counterparts , and , leading to relative error values, , 0.0150 and 0.0155, respectively.

Table 1

Expected values , and as a function of for the hard niche case (i.e., ), for parameters and

φ	E[Ximax]	E[Timax]	E[Nimax]	q	0.25	0.5	0.75	0.99
0.1	1.08143	0.07834	0.08617	K_q	1	1	1	2
				P~(Ximax≥K_q)	0.99526	0.99526	0.99526	0.08617
				E[T_i,K_q|Ximax≥K_q]	0	0	0	0.90909
				σ[T_i,K_q|Ximax≥K_q]	0	0	0	0.90909
				E[N_i,K_q|Ximax≥K_q]	0	0	0	1
				σ[N_i,K_q|Ximax≥K_q]	0	0	0	0
1.0	1.84694	0.52477	1.07984	K_q	1	2	3	5
				P~(Ximax≥K_q)	0.99351	0.49351	0.24351	0.02292
				E[T_i,K_q|Ximax≥K_q]	0	0.5	1.25	2.79412
				σ[T_i,K_q|Ximax≥K_q]	0	0.5	1.03078	2.03962
				E[N_i,K_q|Ximax≥K_q]	0	1	2.5	6.29412
				σ[N_i,K_q|Ximax≥K_q]	0	0	0.86603	2.66205
5.0	8.40494	12.72341	67.15374	K_q	4	10	12	16
				P~(Ximax≥K_q)	0.74812	0.54931	0.32916	0.00999
				E[T_i,K_q|Ximax≥K_q]	0.72892	9.43564	21.56767	39.98014
				σ[T_i,K_q|Ximax≥K_q]	0.52755	8.08455	19.77093	37.53902
				E[N_i,K_q|Ximax≥K_q]	4.02410	49.43688	113.44516	214.03020
				σ[N_i,K_q|Ximax≥K_q]	1.43975	39.26485	98.40159	191.6886
10.0	21.50413	927.95677	9.29×10^3	K_q	22	24	26	29
				P~(Ximax≥K_q)	0.80227	0.61532	0.25744	0.01306
				E[T_i,K_q|Ximax≥K_q]	209.05940	736.272000	1745.83800	2435.43200
				σ[T_i,K_q|Ximax≥K_q]	206.35540	733.346400	1742.73800	2432.12500
				E[N_i,K_q|Ximax≥K_q]	2.09×10^3	7.37×10^3	1.75×10^4	2.44×10^4
				σ[N_i,K_q|Ximax≥K_q]	2.06×10^3	7.33×10^3	1.74×10^4	2.43×10^4
20.0	48.33950	1.01×10^7	2.02×10^8	K_q	50	51	53	56
				P~(Ximax≥K_q)	0.79967	0.65553	0.24416	0.01118
				E[T_i,K_q|Ximax≥K_q]	3.81×10^6	7.71×10^6	1.88×10^7	2.51×10^7
				σ[T_i,K_q|Ximax≥K_q]	3.81×10^6	7.71×10^6	1.88×10^7	2.51×10^7
				E[N_i,K_q|Ximax≥K_q]	7.62×10^7	1.54×10^8	3.77×10^8	5.03×10^8
				σ[N_i,K_q|Ximax≥K_q]	7.62×10^7	1.54×10^8	3.77×10^8	5.03×10^8
30.0	75.17354	1.48×10^11	4.44×10^12	K_q	77	78	80	83
				P~(Ximax≥K_q)	0.82173	0.67261	0.24499	0.01096
				E[T_i,K_q|Ximax≥K_q]	5.25×10^10	1.10×10^11	2.73×10^11	3.62×10^11
				σ[T_i,K_q|Ximax≥K_q]	5.25×10^10	1.10×10^11	2.73×10^11	3.62×10^11
				E[N_i,K_q|Ximax≥K_q]	1.57×10^12	3.28×10^12	8.19×10^12	1.09×10^13
				σ[N_i,K_q|Ximax≥K_q]	1.57×10^12	3.28×10^12	8.19×10^12	1.09×10^13
40.0	102.01704	2.44×10^15	9.75×10^16	K_q	104	105	107	110
				P~(Ximax≥K_q)	0.83761	0.69029	0.25417	0.01138
				E[T_i,K_q|Ximax≥K_q]	8.06×10^14	1.72×10^15	4.42×10^15	5.93×10^15
				σ[T_i,K_q|Ximax≥K_q]	8.06×10^14	1.72×10^15	4.42×10^15	5.93×10^15
				E[N_i,K_q|Ximax≥K_q]	3.22×10^16	6.88×10^16	1.77×10^17	2.37×10^17
				σ[N_i,K_q|Ximax≥K_q]	3.22×10^16	6.88×10^16	1.77×10^17	2.37×10^17
50.0	128.84185	4.28×10^19	2.14×10^21	K_q	131	132	134	137
				P~(Ximax≥K_q)	0.85166	0.70967	0.26879	0.01220
				E[T_i,K_q|Ximax≥K_q]	1.30×10^19	2.84×10^19	7.60×10^19	1.04×10^20
				σ[T_i,K_q|Ximax≥K_q]	1.30×10^19	2.84×10^19	7.60×10^19	1.04×10^20
				E[N_i,K_q|Ximax≥K_q]	6.52×10^20	1.42×10^21	3.80×10^21	5.19×10^21
				σ[N_i,K_q|Ximax≥K_q]	6.52×10^20	1.42×10^21	3.80×10^21	5.19×10^21
75.0	196.43319	2.06×10^30	1.55×10^32	K_q	199	200	201	205
				P~(Ximax≥K_q)	0.83576	0.67186	0.44024	0.00960
				E[T_i,K_q|Ximax≥K_q]	7.42×10^29	1.58×10^30	2.77×10^30	4.97×10^30
				σ[T_i,K_q|Ximax≥K_q]	7.42×10^29	1.58×10^30	2.77×10^30	4.97×10^30
				E[N_i,K_q|Ximax≥K_q]	5.57×10^31	1.19×10^32	2.07×10^32	3.73×10^32
				σ[N_i,K_q|Ximax≥K_q]	5.57×10^31	1.19×10^32	2.07×10^32	3.73×10^32
100.0	264.12728	1.01×10^41	1.11×10^43	K_q	267	268	269	273
				P~(Ximax≥K_q)	0.81774	0.63708	0.39927	0.00800
				E[T_i,K_q|Ximax≥K_q]	4.59×10^40	9.55×10^40	1.61×10^41	2.68×10^41
				σ[T_i,K_q|Ximax≥K_q]	4.59×10^40	9.55×10^40	1.61×10^41	2.68×10^41
				E[N_i,K_q|Ximax≥K_q]	4.59×10^42	9.55×10^42	1.61×10^43	2.68×10^43
				σ[N_i,K_q|Ximax≥K_q]	4.59×10^42	9.55×10^42	1.61×10^43	2.68×10^43

Table 2

Expected values , and for different choices of in the intermediate niche case for the parameters and the initial clonal size

ν	⟨n⟩	E[Ximax]	E[Timax]	E[Nimax]	q	0.25	0.5	0.75	0.99
1	1	123.25013	1.57×10^18	7.68×10^19	K_q	127	128	129	133
					P~(Ximax≥K_q)	0.78258	0.60668	0.38305	0.00847
					E[T_i,K_q|Ximax≥K_q]	7.22×10^17	1.43×10^18	2.33×10^18	3.83×10^18
					σ[T_i,K_q|Ximax≥K_q]	7.22×10^17	1.43×10^18	2.33×10^18	3.83×10^18
					E[N_i,K_q|Ximax≥K_q]	3.54×10^19	7.00×10^19	1.14×10^20	1.88×10^20
					σ[N_i,K_q|Ximax≥K_q]	3.54×10^19	7.00×10^19	1.14×10^20	1.88×10^20
	50	68.61117	2.16×10^8	5.92×10^9	K_q	71	73	75	79
					P~(Ximax≥K_q)	0.83010	0.60581	0.27023	0.00942
					E[T_i,K_q|Ximax≥K_q]	5.97×10^7	1.86×10^8	3.75×10^8	5.22×10^8
					σ[T_i,K_q|Ximax≥K_q]	5.97×10^7	1.86×10^8	3.75×10^8	5.22×10^8
					E[N_i,K_q|Ximax≥K_q]	1.64×10^9	5.10×10^9	1.03×10^10	1.43×10^10
					σ[N_i,K_q|Ximax≥K_q]	1.64×10^9	5.10×10^9	1.03×10^10	1.43×10^10
	100	57.61605	2.24×10^7	5.15×10^8	K_q	60	62	63	67
					P~(Ximax≥K_q)	0.78631	0.50329	0.32493	0.01130
					E[T_i,K_q|Ximax≥K_q]	8.71×10^6	2.54×10^7	3.59×10^7	5.44×10^7
					σ[T_i,K_q|Ximax≥K_q]	8.71×10^6	2.54×10^7	3.59×10^7	5.44×10^7
					E[N_i,K_q|Ximax≥K_q]	2.01×10^8	5.85×10^8	8.27×10^8	1.25×10^9
					σ[N_i,K_q|Ximax≥K_q]	2.01×10^8	5.85×10^8	8.27×10^8	1.25×10^9
10	1	75.87474	4.57×10^9	1.82×10^11	K_q	92	95	97	102
					P~(Ximax≥K_q)	0.75703	0.58015	0.30540	0.00584
					E[T_i,K_q|Ximax≥K_q]	5.38×10^8	3.50×10^9	8.10×10^9	1.31×10^10
					σ[T_i,K_q|Ximax≥K_q]	5.38×10^8	3.50×10^9	8.10×10^9	1.31×10^10
					E[N_i,K_q|Ximax≥K_q]	2.15×10^10	1.40×10^11	3.23×10^11	5.24×10^11
					σ[N_i,K_q|Ximax≥K_q]	2.15×10^10	1.40×10^11	3.23×10^11	5.24×10^11
	50	1.11079	0.09735	0.11315	K_q	1	1	1	3
					P~(Ximax≥K_q)	0.99869	0.99869	0.99869	0.01028
					E[T_i,K_q|Ximax≥K_q]	0	0	0	1.48326
					σ[T_i,K_q|Ximax≥K_q]	0	0	0	1.14561
					E[N_i,K_q|Ximax≥K_q]	0	0	0	2.10206
					σ[N_i,K_q|Ximax≥K_q]	0	0	0	0.33538
	100	1.04909	0.04937	0.05227	K_q	1	1	1	2
					P~(Ximax≥K_q)	0.99682	0.99682	0.99682	0.05227
					E[T_i,K_q|Ximax≥K_q]	0	0	0	0.94455
					σ[T_i,K_q|Ximax≥K_q]	0	0	0	0.94455
					E[N_i,K_q|Ximax≥K_q]	0	0	0	1
					σ[N_i,K_q|Ximax≥K_q]	0	0	0	0
50	1	2.99647	1.06499	4.53753	K_q	1	2	3	20
					P~(Ximax≥K_q)	0.99229	0.49229	0.32338	0.00242
					E[T_i,K_q|Ximax≥K_q]	0	0.5	1.00676	10.60944
					σ[T_i,K_q|Ximax≥K_q]	0	0.5	0.82202	6.66161
					E[N_i,K_q|Ximax≥K_q]	0	1	2.33784	77.95166
					σ[N_i,K_q|Ximax≥K_q]	0	0	0.67229	43.81203
	50	1.01919	0.01920	0.01959	K_q	1	1	1	2
					P~(Ximax≥K_q)	0.99959	0.99959	0.99959	0.01959
					E[T_i,K_q|Ximax≥K_q]	0	0	0	0.98000
					σ[T_i,K_q|Ximax≥K_q]	0	0	0	0.98000
					E[N_i,K_q|Ximax≥K_q]	0	0	0	1
					σ[N_i,K_q|Ximax≥K_q]	0	0	0	0
	100	1.00990	0.00990	0.01000	K_q	1	1	1	2
					P~(Ximax≥K_q)	0.99990	0.99990	0.99990	0.01000
					E[T_i,K_q|Ximax≥K_q]	0	0	0	0.98990
					σ[T_i,K_q|Ximax≥K_q]	0	0	0	0.98990
					E[N_i,K_q|Ximax≥K_q]	0	0	0	1
					σ[N_i,K_q|Ximax≥K_q]	0	0	0	0

Table 3

Expected values , and for different choices of in the soft niche case, for the parameters and the initial clonal size

ν	⟨n⟩	E[Ximax]	E[Timax]	E[Nimax]	q	0.25	0.5	0.75	0.99
500	1	3.64779	1.36144	9.36164	K_q	1	1	3	44
					P~(Ximax≥K_q)	0.99048	0.99048	0.32293	0.00066
					E[T_i,K_q|Ximax≥K_q]	0	0	1.00166	23.19127
					σ[T_i,K_q|Ximax≥K_q]	0	0	0.81786	14.13329
					E[N_i,K_q|Ximax≥K_q]	0	0	2.33378	363.97348
					σ[N_i,K_q|Ximax≥K_q]	0	0	0.66722	222.16995
	50	1.01882	0.01884	0.01922	K_q	1	1	1	2
					P~(Ximax≥K_q)	0.99961	0.99961	0.99961	0.01922
					E[T_i,K_q|Ximax≥K_q]	0	0	0	0.98039
					σ[T_i,K_q|Ximax≥K_q]	0	0	0	0.98039
					E[N_i,K_q|Ximax≥K_q]	0	0	0	1
					σ[N_i,K_q|Ximax≥K_q]	0	0	0	0
	100	1.00970	0.00971	0.00980	K_q	1	1	1	1
					P~(Ximax≥K_q)	0.99010	0.99010	0.99010	0.99010
					E[T_i,K_q|Ximax≥K_q]	0	0	0	0
					σ[T_i,K_q|Ximax≥K_q]	0	0	0	0
					E[N_i,K_q|Ximax≥K_q]	0	0	0	0
					σ[N_i,K_q|Ximax≥K_q]	0	0	0	0
1000	1	1.54331	0.34871	0.64267	K_q	1	1	2	6
					P~(Ximax≥K_q)	0.99225	0.99225	0.32536	0.00795
					E[T_i,K_q|Ximax≥K_q]	0	0	0.66689	2.77159
					σ[T_i,K_q|Ximax≥K_q]	0	0	0.66689	1.77604
					E[N_i,K_q|Ximax≥K_q]	0	0	1	7.28252
					σ[N_i,K_q|Ximax≥K_q]	0	0	0	2.52032
	50	1.00970	0.00971	0.00980	K_q	1	1	1	1
					P~(Ximax≥K_q)	0.99010	0.99010	0.99010	0.99010
					E[T_i,K_q|Ximax≥K_q]	0	0	0	0
					σ[T_i,K_q|Ximax≥K_q]	0	0	0	0
					E[N_i,K_q|Ximax≥K_q]	0	0	0	0
					σ[N_i,K_q|Ximax≥K_q]	0	0	0	0
	100	1.00493	0.00493	0.00495	K_q	1	1	1	1
					P~(Ximax≥K_q)	0.99503	0.99503	0.99503	0.99503
					E[T_i,K_q|Ximax≥K_q]	0	0	0	0
					σ[T_i,K_q|Ximax≥K_q]	0	0	0	0
					E[N_i,K_q|Ximax≥K_q]	0	0	0	0
					σ[N_i,K_q|Ximax≥K_q]	0	0	0	0

The hard niche case, which corresponds to a clonotype (or RTE) with no competitors (), in Table 1, with numerical results based on the stochastic descriptors , and , for various small-to-moderate homeostatic signalling rates, . The hard niche case is also analysed in Fig. 3, where the probability mass function of is plotted for various homeostatic signalling rates, ,

Fig. 3

Probability mass function of for different values of in the hard niche case (i.e., ), for and the initial clonal size

The intermediate niche case in Table 2 and Fig. 4, which corresponds to moderate average numbers of competitors, , and moderate homeostatic signalling rates, ,

Fig. 4

Probability mass function of as a function , in the intermediate niche case for (from top to bottom) , the parameters and the initial clonal size

The soft niche case in Table 3 and Fig. 5, which reflects large average numbers of competitors, , and high homeostatic signalling rates, .

Fig. 5

Probability mass function of in the soft niche case, for , the parameters , and the initial clonal size

The distribution of is analysed in Tables 1, 2 and 3 by computing its (100q)th percentiles for , where the percentile is the first value , such that . The probability of reaching the clonal sizes represented by can be exactly obtained from Eq. (4). However, these probabilities are replaced in Tables 1, 2 and 3 by their truncated versions . We have done so as these truncated values result in good approximations to the true probabilities, but with the advantage that, for the particular cases yielding (for example, and in Table 1), they provide the total probability mass of considered in the truncations given by Eqs. (7) and (10).

The time and the number of proliferation events to reach the maximum clonal size can be analysed in greater depth by studying the conditional values , , and , where denotes the standard deviation of the random variable X. We note here that we consider conditional values due to the fact that the random variables and are defective; that is, we have if the clonal size is not reached, which occurs with non-zero probability. These values provide information about the time and the number of proliferation events to reach different clonal sizes ( with ), under the assumption that those clonal sizes have been, in fact, reached.

When focusing on the hard niche case (Table 1; Fig. 3), the scenarios that allow the establishment of the RTE in the periphery are given by values of and significantly larger than 1, which correspond to larger values of the homeostatic signalling rate . It is clear that the mean value, , has a clearly increasing behaviour with respect to the signalling rate . Moreover, results corresponding to those hostile environments within the hard niche case, corresponding to low signal rates, represent the pressure against the survival and expansion of the clonotype in the short-term. In those cases (), the most representative sample path is the one which amounts to the immediate extinction of the clonotype from the initial clonal size , which gives , or the occurrence of a few proliferation events (yielding a maximum size of the clonotype near 1), and the subsequent extinction of the clonotype. On the other hand, in those cases where the clonotype survives in the mid-term (higher signalling rates), we observe a high concentration of the probability mass function of around some mode, which is reflected in the concentrated values obtained for the percentiles , , and . In particular, the probability mass function of is unimodal when the clonotype becomes extinct in the short-term with high probability, and bimodal when the clonotype survives and expands in the mid-term, with moderate probabilities.

The bimodal shape of in friendly, yet competitive peripheral environments, within the hard niche case, becomes more apparent when analysing results plotted in Fig. 3. In particular, we plot in Fig. 3 the probability mass function of for the hard niche case for signalling rates . The clonotype under consideration has higher potential for expansion for larger values of the homeostatic signalling rate . However, in these cases, the probability mass function of shows a bimodal shape, with a mode equal to 1 and the other taking different values depending on the particular parameters. This shape, together with the results of Table 1, should be interpreted as a binary outcome scenario: in these cases, the clonotype has moderate or high probabilities of reaching its potential clonal size represented by the second mode, which will only happen if it escapes the extinction at the beginning of its lifetime; on the other hand, there exists, with some significant probability, the chance of the clonotype becoming extinct in the short-term (during its first transitions), which yields the first mode equal to 1. We point out here that this also means that (for example, see the case in Fig. 3), once the clonotype escapes extinction in the short-term, the probabilities of reaching different clonal sizes, such as , are practically equal, since they represent the probability of reaching clonal sizes around the second mode. Finally, when the rate of signal is not enough (see the case in Fig. 3), the probability mass function of is unimodal around 1. This represents the case in which the clonotype will become immediately extinct (or extinct in the short-term after reaching small sizes around 1) from its initial clonal size 1.

The increase of associated with increasing values of the homeostatic signalling rate in Table 1 leads also to increasing mean values for , , and . This is easily explained by noting that, in these cases, the clonotype has a greater potential to reach bigger clonal sizes, so that the time and number of proliferation events to reach those sizes should also be larger. The variability of these random variables increases, as well, which is observed by analysing the values and . Moreover, the coefficient of variance of these random variables seems to slowly increase to one, when increasing . In any case, the large values of and imply that, although in friendly conditions a substantial clonotype maximum size will be reached with high probability, there exists a significant impediment against this event, represented by the large number of proliferation and cell death events taking place until this maximum size is reached.

The intermediate niche case is analysed in Table 2 and Fig. 4. In this case, friendly environments can be identified with small values of ( and in Table 2), or with moderate values of together with small values of ( and in Table 2), while hostile scenarios correspond to greater numbers of competitors and larger clonal sizes for them. Again, the clonotype (or RTE) can become established in the periphery under friendly enough environments, reaching a maximum clonal size near , represented by the larger mode of the distribution of (given by the value around which its percentiles are concentrated). Under too hostile scenarios, the clonotype will not become established and its maximum clonal size remains near 1. Again, this binary outcome scenario can be properly identified by analysing the results in Fig. 4, where the probability mass function of is plotted for different parameter regimes for the intermediate niche case. It is worth noting that even in the case in which the most probable outcome is the establishment of the clonotype (e.g., and in Fig. 4), there exists a non-negligible probability of extinction in the short-term. This probability is identified with the probability mass of around the value 1.

We note here that there exist several scenarios where a certain expected behaviour can be identified. For example, in the particular case in Table 2, we obtain the following percentiles of : , which represent a high concentration of the probability mass function around the value 1. That is, there is a high probability of the clonotype becoming immediately extinct from its initial clonal size . In particular, for percentiles , the corresponding time to reach the clonal size is equal to 0, since the initial clonal size is already , and the standard deviation equals 0. Similar comments can be made for . The number of proliferation events to reach (given that it is reached) is equal to 1, with probability 1, and then its standard deviation is 0. It should be noted that the distribution of the random time , under the assumption that , can be seen as the distribution of the minimum between two independent and exponentially distributed random variables X and Y with means and in the case . As a result, is exponentially distributed with mean , and . Moreover, we note that one should carefully interpret the value of in Tables 1, 2 and 3, since these values have been obtained making use of the truncation procedure described in Eqs. (7) and (10).

Finally, the soft niche case representing a hostile environment given by a large number of competitors and big clonal sizes for them, is analysed in Table 3 and Fig. 5. Results in Table 3 represent scenarios where the most probable outcome is the short-term extinction of the clonotype (values of , and near 1, percentiles accumulated around 1). This outcome can be identified by analysing the probability mass function of for these cases, which is concentrated around 1 (see Fig. 5).

In Figs. 6, 7, 8 and 9 and Table 4, we focus on the auxiliary descriptors , with , and , which have not been analysed previously. In particular, we consider in Figs. 6, 7, and 8 the intermediate niche case and we analyse the time to contraction to a given clonal size . The mean time, , and the standard deviation, , are computed for different values of i (, 100 and 500 for Figs. 6, 7 and 8, respectively), , and parameters , and . The results in these figures show relatively little dependence of this time to contraction on the competition parameters (). The initial i and final clonal sizes, and the signalling rate seem to have a higher impact on the behaviour of this random variable, since larger values of are related to a higher resistance of the clonotype against contraction.

Fig. 6

Expected value plus and minus standard deviation for , in the intermediate niche case, for and various choices of , , and

Fig. 7

Expected value plus and minus standard deviation for , in the intermediate niche case, for and various choices of , , and

Fig. 8

Expected value plus and minus standard deviation for , in the intermediate niche case, for and various choices of , , and

Fig. 9

Conditional expected value plus and minus conditional standard deviation for (from top to bottom) , in the soft niche case, and for various choices of the number D of proliferation events, and the parameters , and

Table 4

Probabilitiy as a function of the initial clonal size i in the soft niche case, for various choices of the number D of proliferation events, and the parameters , and

i	D	⟨n⟩	1	10	100	1000
1	1		0.83263	0.33322	0.04761	0.00497
	5		0.79944	0.05986	<10^-5	<10^-9
	10		0.79876	0.01653	<10^-9	<10^-18
10	1		0.99999	0.98249	0.38600	0.04865
	5		0.99999	0.74816	0.00113	<10^-7
	10		0.99999	0.41791	<10^-6	<10^-16
100	1		1	1	0.99230	0.39263
	5		1	0.99999	0.58648	0.00022
	10		1	0.99999	0.05467	<10^-9

In Table 4 and Fig. 9, the time to reach a particular number of proliferation events, D, is studied for the soft niche case. In particular, the probabilities are given in Table 4, and the conditional mean values and conditional standard deviations are shown in Fig. 9 for different values of i, D and . In contrast to the time to contraction, the results of Table 4 and Fig. 9 indicate that (which is a measure of the proliferation capacity of the clonotype) seems to be highly sensitive to , the characteristic number of naive T cells in a competitor clonotype. When the competition environment is too hostile, which is represented by large values of , the potential of the clonotype in the periphery to reach a certain number, D, of proliferation events is notably reduced, and this is reflected by the decrease of the probabilities . At the same time, means and standard deviations in Fig. 9 increase in those scenarios. On the other hand, for small values of , there exist parameter regimes where the target number of divisions, D, is easily reached (see, for example, the cases associated with the initial clonal size in Table 4).

Discussion

In this paper we have analysed the fate of a given naive T cell clonotype or a recent thymic emigrant (with initial clonal size ), when competing with pre-established peripheral T cell clonotypes for homeostatic proliferation signals. The dynamics of a given clonotype has been based on a mathematical stochastic model originally developed in Stirk et al. (2008). Here, our objective is to study the fate and potential of naive T cell clonotypes in the periphery by analysing, for a given clonotype, its maximum clonal size, the time to reach this maximum, the number of proliferation events required to reach this maximum, the rate of contraction of the clonotype during its way to extinction, as well as the time to a given number of proliferation events.

We have introduced in Sect. 2 several stochastic descriptors in terms of defective and non-defective random variables. The use of Laplace-Stieltjes transforms, probability generating functions, auxiliary random variables and absorbing augmented Markov processes, together with first-step analysis, allow us to obtain, in some cases, the distribution of the random variables of interest, or to exactly compute or approximate the different order moments of these variables.

In Sect. 3 several numerical experiments have been carried out in order to study the descriptors previously defined, under different environmental conditions for the clonotype under consideration. In particular, when a recent thymic emigrant reaches the periphery, it will compete with a small or large number of different clonotypes depending on the TCR it expresses (which in turn determines its cross-reactivity). The parameter that encodes the number of competitors is in the model. The characteristic clonal size of these competitors, , as well as the total rate of stimulatory signal provided by the self-peptides that the given clonotype can recognise, , also affect the dynamics of the clonotype under study, so that three different regimes are defined (hard, soft or intermediate niche). For these three different scenarios, numerical results have been provided for our descriptors, which allow us to follow the survival and extinction dynamics of the clonotype or RTE. Our main results can be summarised as follows:Recent experimental work suggests that whether a recent thymic emigrant survives or not upon its arrival in the periphery, depends on the clonal sizes of the pre-existing and competing naive T cell clonotypes: under lymphopenic conditions, that is, small values of , RTEs seem to be able to proliferate and become established in the periphery, but when the peripheral naive T cell pool is full, RTEs struggle to get incorporated into the population of recirculating lymphocytes (Fink and Hendricks 2011; Houston et al. 2011; Berkley and Fink 2014). This empirical observation is in line with our results: the fate of a recent thymic emigrant during its initial journey in the periphery has a clear stochastic component, where the probability of extinction cannot be neglected, even under favourable conditions for the RTE, and which can be studied in terms of the mass probability function of . On the other hand, a greater deterministic component can be identified in the potential size of the clonotype seeded by the RTE in the long-term, once it escapes extinction (which is more likely for small values of ), which is reflected in the concentration of the mass probability function of around the second mode.

Two fates can be identified for the dynamics of the clonotype: if the clonotype under study does not receive enough stimulatory signal, or if it faces too hostile a competitive environment (encoded by a large number of competitors and/or large clonal sizes of these competitors), the survival probability of the clonotype in the mid-term is low. This fate can be seen by analysing the results in Fig. 4 for , or Fig. 5. If the clonotype faces moderate competitive environments, the probability of the clonotype escaping extinction and surviving in the long-term is significant, and it can be approximated in terms of the probability mass function of .

In this second case, the bimodal structure of the probability mass function of (where one mode is 1 or near 1, and the other mode is far away from it) reflects a binary outcome scenario: For example, in Fig. 4 for and , the probability mass around 1 represents the probability of the clonotype becoming extinct after a few proliferation and death events (then, the maximum clonal size will be 1 or a value near 1). The long term survival probability of the clonotype can be approximately identified with . We note here that the probability mass for states between 6 and 110 is almost null, which reflects the fact that once the clonotype escapes extinction, it will reach maximum clonal sizes above 110 almost surely. Moreover, the high concentration of the probability mass function around the second mode (which can also been identified by analysing the percentiles of in Tables 1, 2 and 3) should be interpreted as the fact that the set of parameters (and, thus, the molecular properties of the corresponding TCR and the clonal sizes of competing T cell clonotypes) highly determines the maximum size of a clonotype in the periphery.

One possibility is the extinction of the clonotype in the short-term, with non-negligible probability. The probability of this event can be identified with the sum of the probabilities around 1 in the probability mass function of .

If such an event does not occur, the clonotype escapes extinction and it reaches its maximum size, which will be a value near the second mode of the probability mass function of .

Other interesting results that have been derived from our study of the stochastic descriptors defined here are:

There seems to exist some threshold value for the parameter when analysing the time and the number of proliferation events to reach the maximum clonal size. This threshold behaviour can be observed, for example, in Table 1 for parameters and (hard niche case). It can be observed that, for these cases, small changes in are, however, associated with large changes in and . These results suggest that different clonotypes with similar maximum clonal sizes can display very different behaviours as to how this maximum size is achieved (average time and number of divisions), depending on the value of , the homeostatic proliferation signalling rate.

Figures 6, 7 and 8 allow us to show that the rate of contraction of the given clonotype on its way to extinction seems to depend much more on , the homeostatic proliferation signalling rate, than on the competition parameters . At the same time, the competition parameters do significantly affect the proliferation rate of the clonotype, which is analysed in Table 4 and Fig. 9, in terms of the time to reach a particular number of proliferation events.