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 3×10 4 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 Rocha 1999, 2000;Tanchot and Rocha 1998). An adult mouse has a homeostatic population of 10 8 naive T cells (Mason 1998) and a TCR diversity estimated around 2× 10 6 (Zarnitsyna et al. 2013), and a healthy adult human has a homeostatic population of 10 11 naive T cells (Hazenberg et al. 2003) distributed in 2.5 × 10 7 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 pathogenderived 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(Stirk et al. , 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.

2 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 {X (t) : t ≥ 0} on the state space X = N ∪ {0}, 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 i ∈ N ∪ {0} by Stirk et al. (2008): Parameters ν ≥ 0, μ > 0, ϕ > 0, and n ≥ 1 were introduced in Stirk et al. (2008), and they have the following meaning: • μ 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 • n 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 n , 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)).
It was shown in Stirk et al. (2008) that, for any value of the parameters, 0 ∈ X is an absorbing state of the stochastic process, the time to extinction (or absorption) from any other state i ∈ X is finite with probability one, and its mean is also finite.
As introduced above, parameter ϕ is a measure both of the affinity and the crossreactivity 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 Fig. 1 Birth and death process (with absorption) for the time evolution of the number of cells in a given naive T cell clonotype 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 (ν 1), the second is the limit of large inter-clonal competition (ν 1), and the third one is the intermediate regime of competition (ν ≈ 1).
Hard niche clonotype In the special case ν 1, 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, ν 1, 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 ν ≈ 1, 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 variables for i ∈ X , 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 X max 0 = 0 and T max 0 = 0, since 0 is an absorbing state in X . In order to study both variables in parallel, we define the time to reach a total clonal size i , given that the initial clonal size equals i, as the auxiliary random variable We point out here that these auxiliary random variables have their own immunological interest, since attaining a given clonal size i 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 T i,i also allows the study of the random variables of interest X max i and T max i . For a particular clonal size i ∈ X , the random variable T i,i might be analysed in a different manner depending on whether i < i (time to contraction to a given clonal size i ), or i > i (time to expansion to a given clonal size i ). We note that in the latter case, reaching i is not certain, so that T i,i is a defective random variable, while reaching i 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 i < i, given the initial clonal size i
For an initial clonal size i with i > i , the random variable T i,i amounts to the time until absorption into i for a birth and death process defined on a single absorbing state i and the class of transient states, {i + 1, i + 2, . . . }, with birth rates {λ j : j ≥ i + 1} and death rates {μ j : j ≥ i + 1}. Since absorption of the underlying process {X (t) : t ≥ 0} 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 T i,i from the equality where ρ i = 1 and ρ n = n k=i +1 λ −1 k μ k , for n ≥ i +1; note that the expected values E[T k i,0 ] corresponding to i = 0 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 ∞ j=n+1 becomes the finite sum N j=n+1 and T i,i 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 i > i, given the initial clonal size i, and analysis of X max i and T max i
For a fixed (and given) clonal size i ∈ N, we introduce the following notation: where 1 {T i,i <∞} is a random variable that takes the values 1 if T i,i < ∞, and 0 otherwise. The quantities defined in Eq.
(2) will allow us to analyse the random variables X max i and T max i , but the distribution of T i,i is defective, since T i,i = ∞ if the clonotype becomes extinct before reaching the size i ; note that 1 − v i,i is the probability of such an extinction event, as 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 as

123
Then, by using forward elimination, we may obtain In evaluating the restricted moments m (k) i,i , we first derive the probabilities v i,i as the values φ i,i (0) from Eq. (4). In particular, it is seen that v 0,i = 0 and we can write (3). In fact, by taking derivatives in Eq. (3) we obtain In order to solve Eq. (5), for a fixed value k ≥ 1, we introduce some notation as follows: Equation (6) is a recursive procedure that allows us to compute the kth order moments, In the special cases k = 0 and 1, it is readily seen that for any value of the parameter ν.
We finally focus on the random variables X max For practical or computational purposes, the above series should be replaced by the finite sum where K q can be selected as the (100q)th percentile of X max i for a probability q ∈ (0, 1) 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 X max i is greater than q, whence the probability 1 − q 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 N max i of one-step transitions i → i + 1 of X occurring in the random interval [0, T max i ], which allows us to record the number of proliferation events to reach the maximum clonal size, for i ∈ X . This descriptor can be seen as a discrete version of the random variable T max i 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 N i,i of one-step transitions i → i + 1 of X to register X (t) = i for the first time.
We first introduce the following notation: and, consequently, the probabilities u i,i can be computed from Eq. (4), and they amount to the probabilities of reaching the clonal size i before extinction, given that the initial clonal size is i.
By taking derivatives on the equality which is similar to Eq. (5) with the term k m i+1,i . This implies that, similarly to the solution of Eq. (6), the solution of Eq. (8) has the form Similarly to Sect. 2.1.2, it can be seen that lim i →∞m In the case of the random variable N i,i , not only every kth order factorial moment can be obtained from Eq. (9), but also its distribution. We can writẽ so that a recursion on j and i can be obtained in order to compute the probability mass function of N i,i . 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.

123
Finally, the factorial moments of the random variable N max i can be approximated in a similar way as described for Eq. (7). In particular, if we denotem max,(k)

Time to reach a given number of proliferation events
In this section, we fix an initial clonal size i and a number D ≥ 0 of proliferation events, and we study the time T D i to reach a total number D of proliferation events (one-step transitions i → i + 1 of X ), with T 0 i = 0. In order to study this descriptor, we consider the augmented process {(X (t), D(t)) : t ≥ 0}, where D(t) denotes the number of division events up to time t, with (X (0), D(0)) = (i , 0). For the initial clonal size i , the augmented process starts from state (i , 0) and is defined on the finite state space States in X 0 ∪ X D (i ) are absorbing states. In particular, states in X 0 represent the extinction of the clonotype when the number D of proliferation events has not been reached, while states in X D (i ) reflect that the number D of proliferation events has been reached. States in X T (i ) are transient states. Figure 2 illustrates the dynamics of the augmented process in the case D = 3.
In a more general setting, we define the random variable Note that the descriptor T D i is equivalent to the random variable T D (i ,0) , and states in X 0 and X D (i ) lead to the values T D (0,d) = ∞ and T D (i,D) = 0, respectively. For states in X T (i ), the random variable T D (i,d) is defective, and P(T D (i,d) = ∞) > 0 is the probability of reaching X 0 before getting to X D (i ). For and we observe that As a result, the Laplace-Stieltjes transforms D (i,d) (s) verify the boundary conditions which yields a finite system of linear equations that can be solved in a recursive manner (Algorithm 1). The boundary conditions for the momentsm

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 ν, n 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, n 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 μ = 1, 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 Algorithm 2 Computation of the momentsm Step 1: k := 0; d := D; i := 1; Step 2: niche cases amount to values of the signalling rate ϕ, the number ν of competitors, and the characteristic size n 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: (a) The hard niche case, which corresponds to a clonotype (or RTE) with no competitors (ν = 0), in Table 1 Table 2 and Fig. 4, which corresponds to moderate average numbers of competitors, ν ∈ {1, 10, 50}, and moderate homeostatic signalling rates, ϕ = 50, (c) The soft niche case in Table 3 and Fig. 5, which reflects large average numbers of competitors, ν ∈ {500, 1000}, and high homeostatic signalling rates, ϕ = 500.
In Tables 1, 2  ] are computed by considering only clonal sizes up to the 99th percentile K 0.99 of X max i , and by means of the truncating procedures given by Eqs. (7) and (10). However, for parameter values (ν, n ) ∈ {(500, 100), (1000, 50), (1000, 100)} in Table 3 , since K 0.99 = 1 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 i = 1. 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 (                The distribution of X max i is analysed in Tables 1, 2 and 3 by computing its (100q)th percentiles K q for q ∈ {0.25, 0.5, 0.75, 0.99}, where the percentile K q is the first value x ≥ 1, such that P(X max i ≤ x) ≥ q. The probability P(X max i ≥ K q ) of reaching the clonal sizes represented by K q can be exactly obtained from Eq. (4). However, these probabilities are replaced in Tables 1, 2 and 3 by their truncated versionsP(X max i ≥ K q ) = P(K 0.99 ≥ X max i ≥ K q ). 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 K q = 1 (for example, ϕ ∈ {0.1, 1.0} and q = 0.25 in Table 1), they provide the total probability mass of X max i 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 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 T i,K q and N i,K q are defective; that is, we have T i,K q = N i,K q = ∞ if the clonal size K q 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 (K q with q ∈ {0.25, 0.5, 0.75, 0.99}), 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 E[X max i ] and E[T max i ] significantly larger than 1, which correspond to larger values of the homeostatic signalling rate ϕ. It is clear that the mean value, E[X max i ], 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 (ϕ ∈ {0.1, 1.0, 5.0}), the most representative sample path is the one which amounts to the immediate extinction of the clonotype from the initial clonal size i = 1, which gives X max i = 1, 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 X max i around some mode, which is reflected in the concentrated values obtained for the percentiles K 0.25 , K 0.5 , K 0.75 and K 0.99 . In particular, the probability mass function of X max i 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 X max i 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 X max i for the hard niche case for signalling rates ϕ ∈ {0.1, 1.0, 5.0, 10.0}. 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 X max i 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 ϕ = 10.0 in Fig. 3), once the clonotype escapes extinction in the shortterm, the probabilities of reaching different clonal sizes, such as i ∈ {5, 10, 15}, 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 ϕ = 0.1 in Fig. 3), the probability mass function of X max i 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.
. 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 σ [T i,K q  ] 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 ν (ν = 1 and n ∈ {1, 50, 100} in Table 2), or with moderate values of ν together with small values of n (ν = 10 and n = 1 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 E[X max i ], represented by the larger mode of the distribution of X max i (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 X max i 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., ν = 1 and n = 1 in Fig. 4), there exists a non-negligible probability of extinction in the short-term. This probability is identified with the probability mass of X max i 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 (ν, n ) = (50, 50) in Table 2, we obtain the following percentiles of X max i : (K 0.25 , K 0.5 , K 0.75 , K 0.99 ) = (1, 1, 1, 2), 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 i = 1. In particular, for percentiles K q = 1, the corresponding time T i,K q to reach the clonal size K q is equal to 0, since the initial clonal size is already i = 1, and the standard deviation equals 0. Similar comments can be made for N i,K q . The number of proliferation events to reach K q = 2 (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 T i,2 , under the assumption that T i,2 < ∞, can be seen as the distribution of the minimum between two independent and exponentially distributed random variables X and Y with means λ −1 1 and μ −1 1 in the case X < Y . As a result, T i,2 is exponentially distributed with mean (λ 1 + μ 1 ) −1 , and 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  ] near 1, percentiles K q accumulated around 1). This outcome can be identified by analysing the probability mass function of X max i for these cases, which is concentrated around 1 (see Fig. 5).
In Figs. 6,7,8 and 9 and  and 8, respectively), i , and parameters ν, ϕ and n . The results in these figures show relatively little dependence of this time to contraction on the competition parameters (ν, n ). The initial i and final i 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.
In Table 4 and Fig. 9, the time T D i to reach a particular number of proliferation events, D, is studied for the soft niche case. In particular, the probabilities P(T D i < ∞)  Fig. 9 for different values of i, D and n . In contrast to the time to contraction, the results of Table 4 and Fig. 9 indicate that T D i (which is a measure of the proliferation capacity of the clonotype) seems to be highly sensitive to n , 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 n , 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 P(T D i < ∞). At the same time, means and standard deviations in Fig. 9 increase in those scenarios. On the other hand, for small values of n , 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 i = 100 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 i = 1), 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, n , 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: 1. 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 ν = 50, 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 X max i . 2. In this second case, the bimodal structure of the probability mass function of X max i (where one mode is 1 or near 1, and the other mode is far away from it) reflects a binary outcome scenario: • One possibility is the extinction of the clonotype in the short-term, with nonnegligible probability. The probability of this event can be identified with the sum of the probabilities around 1 in the probability mass function of X max i .

123
• 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 X max i . For example, in Fig. 4 for ν = 1 and n = 1, 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 P(X max i ≥ 6). 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 X max i in Tables 1, 2 and 3) should be interpreted as the fact that the set of parameters (ϕ, ν, n ) (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.
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 n , 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 X max i . 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 n ), which is reflected in the concentration of the mass probability function of X max i around the second mode.
Other interesting results that have been derived from our study of the stochastic descriptors defined here are: 1. 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  ]. 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. 2. 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 (ν, n ). 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.