Emily R Stirk

Grant Lythe

Hugo A Van Den Berg

Carmen Molina-París
email: carmen@maths.leeds.ac.uk

Stochastic competitive exclusion in the maintenance of the naïve T cell repertoire

Keywords: T cell homeostasis, TCR diversity, stochastic bivariate birth and death process, niche overlap, large N expansion

Recognition of antigens by the adaptive immune system relies on a highly diverse T cell receptor repertoire. The mechanism that maintains this diversity is based on competition for survival stimuli; these stimuli depend upon weak recognition of self-antigens by the T cell antigen receptor. We study the dynamics of diversity maintenance as a stochastic competition process between a pair of T cell clonotypes that are similar in terms of the self-antigens they recognise. We formulate a bivariate continuous-time Markov process for the numbers of T cells belonging to the two clonotypes. We prove that the ultimate fate of both clonotypes is extinction and provide a bound on mean extinction times. We focus on the case where the two clonotypes exhibit negligible competition with other T cell clonotypes in the repertoire, since this case provides an upper bound on the mean extinction times. As the two clonotypes become more similar in terms of the self-antigens they recognise, one clonotype quickly becomes extinct in a process resembling classical competitive exclusion. We study the limiting probability distribution for the bivariate process, conditioned on non-extinction of both clonotypes. Finally, we derive deterministic equations for the number of cells belonging to each clonotype as well as a linear Fokker-Planck equation for the fluctuations about the deterministic stable steady state.

Introduction

The T cell repertoire exhibits a large diversity of its antigen receptors (TCRs); a human possesses approximately 10 7 -10 8 TCRs [START_REF] Arstila | A direct estimate of the human αβ T cell receptor diversity[END_REF]. This huge diversity of TCRs is generated by a random genetic recombination of the TCR encoding genes during the T cell maturation process [START_REF] Krangel | Developmental regulation of V(D)J recombination at the TCR α/δ locus[END_REF]. The rationale for the diversity is the need to have a specific receptor ready for any pathogenic challenge, whose associated antigens are unpredictable. T cells awaiting activation by an antigen are called naïve; there are 10 11 naïve T cells [START_REF] Goronzy | Aging and T cell diversity[END_REF], which indicates that a given TCR is usually present on numerous T cells, which together constitute a clone bearing that TCR clonotype.

The number of naïve T cells remains approximately constant throughout an individual's lifetime [START_REF] Freitas | Population biology of lymphocytes: the flight for survival[END_REF]. This homeostasis is driven by interactions with antigens derived from the body's own proteins (self-antigens) displayed on the surface of specialised antigen-presenting cells (APCs; [START_REF] Ernst | The peptide ligands mediating positive selection in the thymus control T cell survival and homeostatic proliferation in the periphery[END_REF][START_REF] Freitas | Peripheral T cell survival[END_REF][START_REF] Goldrath | Selecting and maintaining a diverse T cell repertoire[END_REF]. Naïve T cells infrequently divide, contingent on survival signals from such an APC. The survival signal depends on the TCR and on the self-antigens presented on the APC surface [START_REF] Ferreira | Differential survival of naïve CD4 and CD8 T cells[END_REF].

Whereas individual clones are susceptible to extinction, diversity and optimally dispersed coverage of antigen space is maintained [START_REF] Correia-Neves | The shaping of the T cell repertoire[END_REF][START_REF] Mahajan | Homeostasis of T cell diversity[END_REF]. Instrumental to this mechanism is the relationship between clonal lifetime and antigenic overlap with other clones. We previously investigated this relationship for moderate antigenic overlap, treating clonal size as a stochastic univariate birth and death process [START_REF] Stirk | Stochastic niche structure and diversity maintenance in the T cell repertoire[END_REF]. We here extend the analysis to the case of non-negligible antigenic overlap.

A c c e p t e d m

a n u s c r i p t

|Q i ∩ Q j | |Q i | for all j = i
There exists j such that

|Q i ∩ Q j | ∼ |Q i | Q i Q i Q j
Figure 1: The diagram on the left shows the repertoire analysed in [START_REF] Stirk | Stochastic niche structure and diversity maintenance in the T cell repertoire[END_REF]. The situation on the right is considered in this paper, where

|Q i ∩ Q k | |Q i | and |Q j ∩ Q k | |Q j | for all k = i, j.
Our earlier paper showed that clonotype competition for survival stimuli from APCs presenting self-antigens maximises repertoire diversity by promoting the survival of T cell clonotypes that are most different from each other, in terms of the survival signals that they are able to receive. This analysis relied on mean-field approximations concerning the competition between distinct clones (Fig. 1, left panel). We assumed that competition between individual pairs of clones is small, even though T cells of a given clone compete with many other clones for access to survival stimuli. This is plausible for the vast majority of TCR pairs in the repertoire; however, very similar TCRs do occur [START_REF] Kedzierska | Early establishment of diverse T cell receptor profiles for influenza-specific CD8 + CD62L hi memory T cells[END_REF]Wynn et al., 2008). In this paper we introduce a stochastic model for the case of a pair of clonotypes with similar TCRs that exhibit substantial overlap in the survival signals for which they compete (Fig. 1, right panel). In this case, the influence of the competing clone cannot be treated as that of an average clone forming part of the "general background" of competition; a bivariate Markov process is required to describe the time evolution of the number of T cells belonging to both clonotypes.

The paper is organised as follows. Section 2 introduces the model for the clonal sizes of pairs of "similar" clonotypes, i.e., clonotypes that compete with each other for a significant fraction of the survival stimuli that they are able to receive. In Section 3 we prove that both clonotypes become extinct in finite time for all parameter values of the model and we provide an upper bound on the mean time until extinction occurs. Section 4 considers a special case of the model. We compare the time until extinction of both clonotypes to the average time until one of the T cell clonotypes becomes extinct. We show that this latter time decreases as the commonality between the two clonotypes increases, resembling the ecological principle of classical competitive exclusion [START_REF] Begon | Ecology: Individuals, Populations and Communities. Blackwell's Sci[END_REF]. Section 5 introduces the limiting conditional probability distribution, which describes the stationary behaviour of the process before extinction occurs. Finally, in Section 6 we exploit van Kampen's "large N expansion" [START_REF] Van Kampen | A power series expansion of the master equation[END_REF]Van Kampen, 2007) for this special case to approximate the limiting conditional probability distribution for both clonotypes.

Notation is summarised in Table 1. All other notation is introduced and explained in the relevant sections.

Stochastic model for a pair of competing clonotypes

Naïve T cells compete for a limited supply of survival stimuli from APCs presenting self-antigens; this regulates the diversity of the T cell repertoire [START_REF] Troy | Cutting edge: homeostatic proliferation of peripheral T lymphocytes is regulated by clonal competition[END_REF]. A single APC may present many different peptides at a given instant, and the particular peptides presented will continually change over time [START_REF] Mahajan | Homeostasis of T cell diversity[END_REF]. Let us call the antigens displayed at a given instant on the surface of an APC its antigen-presentation profile (APP; van den [START_REF] Van Den Berg | Quantifying the strength of ligand antagonism in TCR triggering[END_REF][START_REF] Van Den Berg | Antigen presentation on MHC molecules as a diversity filter that enhances immune efficacy[END_REF]. Let Q be the set of all possible APPs and Q i be the subset of APPs capable of providing a survival stimulus to T cells of clonotype i. Previously, we assumed that competition between any pair of clonotypes is small. This means that for any fixed clonotype i,

|Q i ∩ Q j | |Q i | for i = j.
Thus, all clonotypes other than i are treated as part of the average background of competition [START_REF] Stirk | Stochastic niche structure and diversity maintenance in the T cell repertoire[END_REF]. We here consider the case of a pair of T cell clonotypes, i and j, such that

|Q i ∩ Q j | ∼ |Q i | but |Q i ∩ Q k | |Q i | and |Q j ∩ Q k | |Q j |
for all k = i, j (see Fig. 1). This means that the extinction of a clonotype k = i, j from the repertoire does not significantly affect the dynamics of clonotypes i and j and so the influence of all other clonotypes can be treated as part of the general background of competition. However, the number of T cells belonging to both clonotypes i and j must be explicitly accounted for (see Section 2.3). The set of all T cells in the naïve repertoire C q

A c c e p t e d m a n u s c r i p t

The set of all T cells in the naïve repertoire that can receive a survival stimulus from APP q Q

The set of all APPs that may occur on an APC Q i

The set of APPs from which T cells of clonotype i can receive survival stimuli Q ij

The set of APPs from which T cells of clonotype i and T cells of clonotype j can receive survival stimuli

Q i/j
The set of APPs from which T cells of clonotype i can receive survival stimuli but T cells of clonotype j cannot Q ijr

The set of APPs that provide survival stimuli to T cells of clonotypes i and j and T cells of r other distinct clonotypes

Q ir/j
The set of APPs that provide survival stimuli to T cells of clonotype i and T cells of r other distinct clonotypes, none of which is j n i

Number of naïve T cells belonging to clonotype i n q

Number of naïve T cells that are capable of receiving survival stimuli from APP q n iq Whenever q ∈ Q i , we have n q ≥ n i and define

n iq = n q -n i n ijq For q ∈ Q i and q ∈ Q j , define n ijq = n q -n i -n j N C
The total number of clonotypes in the naïve repertoire n Average clone size over the naïve repertoire γ q

The collective stimulus rate from APCs that present APP q

Λ i
The stimulus rate received by a naïve T cell of clonotype i μ i

The per-cell death rate for T cells of clonotype i p 1

The probability that a T cell of clonotype 2 recognises an APP chosen at random from the set Q 1 p 2

The probability that a T cell of clonotype 1 recognises an APP chosen at random from the set

Q 2 p •|ij
The probability that an APP chosen at random from

Q ij will belong to the set Q k of a different clone k selected at random p •|i/j
The probability that an APP chosen at random from Q i/j will belong to the set Q k of a different clone k selected at random ν 12

Mean niche overlap for APPs in Q 12 ν 1

Mean niche overlap for APPs in Q 1/2 ν 2

Mean niche overlap for APPs in Q 2/1 ti The time at which clonotype i first becomes extant in the repertoire ϕ i

This parameter is defined to be equal to

γ|Q i | τ n1,n2
Expected remaining time until extinction of both clonotypes of present sizes n 1 , n 2 τn1,n2

Expected remaining time until one clonotype becomes extinct when the present clonal sizes are n 1 , n 2 ℘ n1,n2

Probability that clonotype 1 will go extinct before clonotype 2 when the present sizes of the clones are n 1 , n 2 TCR: T cell antigen receptor; APC: antigen-presenting cell; APP: antigen-presentation profile; QSD: quasistationary probability distribution; LCD: limiting conditional probability distribution.

Derivation of the birth and death rates

Our aim is to derive the per-cell birth and death rates for a pair of clonotypes (referred to as clonotypes 1 and 2) which overlap substantially in terms of the APPs from which they are able to receive survival stimuli. Let C denote the set of all T cells in the naïve repertoire and C q the subset of T cells capable of receiving a survival stimulus from APP q (see Fig. 2). The number of T cells capable of receiving a survival stimulus from APP q equals |C q | = n q . Let n 1 and n 2 denote the number of T cells belonging to clonotypes 1 and 2, respectively.

We assume that the survival stimuli from APP q are distributed equally among all naïve T cells capable of forming a stimulatory contact with an APC upon encountering APP q. Let γ q be the rate of survival stimuli from all APCs that present APP q; for simplicity we assume γ q ≡ γ, ignoring the APP fluctuations, which do not significantly affect the results (cf. [START_REF] Van Den Berg | A reliable and safe T cell repertoire based on lowaffinity T cell receptors[END_REF]). If Λ 1 is the stimulus rate received by a single T cell of clonotype 1, partitioning gives

Λ 1 = q∈Q1 γ |C q | = q∈Q1 γ n q = q∈Q1 γ n 1 + n 1q , ( 1
)
where n 1q def = n q -n 1 . This rate depends not only on n 1 , but also on n i for i = 1 through the term n 1q ; the ). An edge from an APP q to a T cell clonotype i indicates that APP q can provide a survival signal to T cells of this clonotype.

dynamics of the various clonotypes in the repertoire is coupled due to the competition for APPs.

Let Q 12 def = Q 1 ∩ Q 2 , ( 2
)
be the set of APPs from which T cells of both clonotype 1 and clonotype 2 receive a survival stimulus and let

Q 1/2 def = Q 1 ∩ Q2 , ( 3
)
be the set of APPs from which T cells of clonotype 1 receive a survival stimulus, but T cells of clonotype 2 do not. The set Ā denotes the complement of the set A in Q. By construction we have

Q 1 = Q 12 ∪ Q 1/2 and Q 12 ∩ Q 1/2 = ∅. Hence, we have Λ 1 = q∈Q12 γ n 1 + n 2 + n 12q + q∈Q 1/2 γ n 1 + n 1q , (4)
where n 12q def = n q -n 1 -n 2 . To characterise the strength of competition from the background clonotypes, let Q 12r denote the set of APPs which provide survival stimuli to T cells of clonotype 1 and clonotype 2 as well as to T cells of r distinct clonotypes in the repertoire, other than clonotypes 1 and 2, and let Q 1r/2 be the set of APPs which provide survival stimuli to T cells of clonotype 1 and to r other distinct clonotypes other than clonotype 1, none of which is clonotype 2. We write

Λ 1 = γ +∞ r=0 q∈Q12r 1 n 1 + n 2 + n 12q + q∈Q 1r/2 1 n 1 + n 1q ≈ γ +∞ r=0 |Q 12r | n 1 + n 2 + r n + |Q 1r/2 | n 1 + r n , ( 5
)
where n is the average clonal size and the approximation for the background clones follows the argument introduced in [START_REF] Stirk | Stochastic niche structure and diversity maintenance in the T cell repertoire[END_REF] According to this approximation, Λ 1 depends only on the number of cells of clonotype 1 and clonotype 2 and n , but not explicitly on the number of T cells belonging to other clonotypes. It remains to find suitable expressions for |Q 12r | and |Q 1r/2 |. Let p 1 be the probability that a randomly chosen APP provides a survival stimulus to T cells of clonotype 2, given that it provides a survival stimulus to T cells of clonotype 1:

p 1 def = |Q 12 | |Q 1 | = |Q 1 ∩ Q 2 | |Q 1 | . ( 6
)
To estimate the cardinality of the sets Q 12r and Q 1r/2 , let N C be the total number of clonotypes extant in the peripheral repertoire, let p .|12 denote the probability that a randomly chosen APP from the set Q 12 belongs to the set Q k of a different extant clone k = 1, 2 and let p .|1/2 be the probability that a randomly chosen APP from the set Q 1/2 belongs to the set Q k of a different extant clone k = 1, 2. The numbers of APPs in Q 12r and Q 1r/2 follow binomial distributions:

|Q 12r | = |Q 12 | N C -2 r (p .|12 ) r (1 -p .|12 ) NC -2-r , ( 7
)
and

|Q 1r/2 | = |Q 1/2 | N C -2 r (p .|1/2 ) r (1 -p .|1/2 ) NC -2-r . ( 8
| ν r 1 e -ν1 r! = (1 -p 1 )|Q 1 | ν r 1 e -ν1 r! . ( 10
)
The parameter ν 12 is the mean niche overlap for APPs in the set Q 12 , which is the average number of clonotypes that are competing with clonotypes 1 and 2 for an APP in the set Q 12 , whereas ν 1 is the mean niche overlap for APPs in the set Q 1/2 , which is the average number of clonotypes that are competing with clonotype 1 for an APP in the set Q 1/2 . Hence, these parameters represent the strength of the competition between T cells of clonotype 1 and T cells of other clonotypes extant in the repertoire (other than clonotype 2) for APPs in these sets, Q 12 and Q 1/2 , respectively. The parameter ν 12 differs from the mean niche overlap parameter of the corresponding univariate model ν [START_REF] Stirk | Stochastic niche structure and diversity maintenance in the T cell repertoire[END_REF]. The overlap parameters are related by

ν = ν 12 p 1 + ν 1 (1 -p 1 ) . ( 11
)
Substituting Eqs. ( 9) and (10) into Eq. ( 5) results in

Λ 1 = ϕ 1 p 1 e -ν12 +∞ r=0 ν r 12 r! 1 n 1 + n 2 + r n + (1 -p 1 )e -ν1 +∞ r=0 ν r 1 r! 1 n 1 + r n , ( 12
)
where ϕ 1 def = γ|Q 1 | is a parameter proportional to the number of APPs from which T cells of clonotype 1 receive a survival stimulus. Similarly, the per-cell birth rate for T cells of clonotype 2 is given by

Λ 2 = ϕ 2 p 2 e -ν12 +∞ r=0 ν r 12 r! 1 n 1 + n 2 + r n + (1 -p 2 )e -ν2 +∞ r=0 ν r 2 r! 1 n 2 + r n , ( 13
)
where

ϕ 2 def = γ|Q 2 |, p 2 def = |Q 12 |/|Q 2 | and ν 2 def = (N C -2)p .|2/1
. We assume a constant per-cell death rate of μ 1 for T cells of clonotype 1 and μ 2 for T cells of clonotype 2. The model is now closed in n 1 and n 2 , accounting explicitly for their mutual competition, and treating competition with all other clones as mean field.

Competition process

We model the number of naïve T cells belonging to clonotypes 1 and 2 at time t, which we denote by (n 1 (t), n 2 (t)), as a continuous-time bivariate Markov process [START_REF] Allen | An Introduction to Stochastic Processes with Applications to Biology[END_REF]. We assume that the thymus produces T cells of a particular clonotype i within a short space of time and denote the time at which this "burst" occurs as ti . The initial number of T cells of clonotype 1 produced by the thymus is given by ñ1 = n 1 ( t1 ) and the initial number of T cells of clonotype 2 produced by the thymus is given by ñ2 = n 2 ( t2 ). Without loss of generality, we may assume that t1 ≤ t2 . For t < t2 only one of the clonotypes is extant in the repertoire, in which case the univariate model introduced in [START_REF] Stirk | Stochastic niche structure and diversity maintenance in the T cell repertoire[END_REF] may be applied. We model the number of T cells of clonotypes 1 and 2 as a continuous-time bivariate Markov process {(X (t), Y(t)) : t ≥ t2 } on the state space S = {(n 1 , n 2 ) : n 1 , n 2 = 0, 1, 2, . . .}. The initial state of the process is thus given by (n 1 ( t2 ), n 2 ( t2 )). Transitions are only allowed to adjacent states [START_REF] Allen | An Introduction to Stochastic Processes with Applications to Biology[END_REF], leading to a two-dimensional analogue of the birth and death process, which we refer to as a (bivariate) competition process [START_REF] Reuter | Competition Process[END_REF]. The transition probabilities are defined by

p nm (Δt) = P{X (t + Δt) = m 1 , Y(t + Δt) = m 2 |X (t) = n 1 , Y(t) = n 2 } (14) for n = (n 1 , n 2 ) ∈ S and m = (m 1 , m 2 ) ∈ S.
These probabilities satisfy the following as Δt → 0 + :

p nm (Δt) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ λ (1) n1,n2 Δt + o(Δt) m = (n 1 + 1, n 2 ) , λ (2) n1,n2 Δt + o(Δt) m = (n 1 , n 2 + 1) , μ (1) n1,n2 Δt + o(Δt) m = (n 1 -1, n 2 ) , μ (2) n1,n2 Δt + o(Δt) m = (n 1 , n 2 -1) , 1 -(λ (1) n1,n2 + λ (2) n1,n2 + μ (1) n1,n2 + μ (2) n1,n2 )Δt + o(Δt) m = (n 1 , n 2 ) , o(Δt)
o t h e r w i s e .

(15)

A schematic representation of the competition process is given in Fig. 3. The quantity λ

(1) n1,n2 is the birth rate for T cells of clonotype 1 and is the rate of transition from state (n 1 , n 2 ) to (n 1 + 1, n 2 ). Similarly, the birth rate for T cells of clonotype 2, denoted λ

(2) n1,n2 , is the rate of transition from state (n 1 , n 2 ) to (n 1 , n 2 + 1). The death rate for T cells of clonotype 1 is given by μ

(1) n1,n2 and this is the rate of transition from state (n 1 , n 2 ) to (n 1 -1, n 2 ). The death rate for T cells of clonotype 2, μ (1, 0)

λ (1) 1,0 μ (1) 2,0
(2, 0)

λ (1) 2,0 μ (1) 3,0 (3, 0) . . . μ (2) 0,1 μ (2) 1,1 μ (2) 2,1 μ (2) 3,1 (0, 1) μ (1) 1,1
(1, 1)

λ (1) 1,1 μ (1) 2,1
(2, 1)

λ (1) 2,1 μ (1) 3,1 (3, 1) . . . μ (2) 0,2 λ (2) 0,1 μ (2) 1,2 λ (2) 1,1 μ (2) 2,2 λ (2) 2,1 μ (2) 3,2 λ (2) 3,1 (0, 2) μ (1) 1,2
(1, 2)

λ (1) 1,2 μ (1) 2,2
(2, 2)

λ (1) 2,2 μ (1) 3,2 (3, 2) . . . μ (2) 0,3 λ (2) 0,2 μ (2) 1,3 λ (2) 1,2 μ (2) 2,3 λ (2) 2,2 μ (2) 3,3 λ (2) 3,2 (0, 3) μ (1) 1,3
(1, 3)

λ (1) 1,3 μ (1) 2,3
(2, 3) to (n 1 , n 2 -1). The birth and death rates for the process are defined by λ

λ (1) 2,3 μ (1) 3,3 (3, 3) . . . . . . . . . . . . . . .
(1)

n1,n2 = Λ 1 n 1 , λ (2) n1,n2 = Λ 2 n 2 , μ (1) n1,n2 = μ 1 n 1 , μ (2) n1,n2 = μ 2 n 2 , so that λ (1) n1,n2 = ϕ 1 n 1 p 1 e -ν12 +∞ r=0 ν r 12 r! 1 n 1 + n 2 + r n + (1 -p 1 )e -ν1 +∞ r=0 ν r 1 r! 1 n 1 + r n , ( 16
)
λ (2) n1,n2 = ϕ 2 n 2 p 2 e -ν12 +∞ r=0 ν r 12 r! 1 n 1 + n 2 + r n + (1 -p 2 )e -ν2 +∞ r=0 ν r 2 r! 1 n 2 + r n , ( 17
)
μ (1) n1,n2 = μ 1 n 1 , (18) μ (2) n1,n2 = μ 2 n 2 . ( 19
)
Since μ

(1) 0,j = μ

(2) j,0 = 0 for all j ≥ 0, transitions out of the state space S cannot occur. Because no T cells of either clonotype 1 or 2 will be produced from the thymus after the time t = t2 , we have λ

(1) 0,j = λ

(2) j,0 = 0 for all j ≥ 0. Hence the set of states A = {(n 1 , n 2 ) : n 1 = 0 or n 2 = 0} forms an absorbing set, and the state (n 1 , n 2 ) = (0, 0) is an absorbing state, which corresponds to the extinction of both clonotypes. The following constraint applies:

ϕ 1 p 1 = ϕ 2 p 2 , ( 20
)
where 0 ≤ p 1 ≤ 1 and 0 ≤ p 2 ≤ 1 (since Q 1 ∩Q 2 = Q 12 = Q 21 = Q 2 ∩Q 1 , whence γ|Q 12 | = γp 1 |Q 1 | = p 1 ϕ 1 , and γ|Q 21 | = γp 2 |Q 2 | = p 2 ϕ 2 )
. The model has nine independent parameters, as compared to the four parameters of the model for a single T cell clonotype.

In the limits p 1 = p 2 = 0 and p 1 = p 2 = 1, the bivariate process reduces to the univariate case: for p 1 = p 2 = 0, the two clonotypes have independent dynamics, while for p 1 = p 2 = 1, the pair together behaves as a single clonotype (see Appendix A). The birth rates are bounded, λ

(i) n1,n2 ≤ ϕ i (i = 1, 2), as follows: λ (1) n1,n2 = ϕ 1 n 1 p 1 e -ν12 +∞ r=0 ν r 12 r! 1 n 1 + n 2 + r n + (1 -p 1 )e -ν1 +∞ r=0 ν r 1 r! 1 n 1 + r n ≤ ϕ 1 n 1 p 1 e -ν12 +∞ r=0 ν r 12 r! 1 n 1 + (1 -p 1 )e -ν1 +∞ r=0 ν r 1 r! 1 n 1 = ϕ 1 n 1 p 1 e -ν12 e ν12 n 1 + (1 -p 1 )e -ν1 e ν1 n 1 = ϕ 1 p 1 + ϕ 1 (1 -p 1 ) = ϕ 1 . ( 21
|Q i ∩ Q j | ∼ |Q i | but |Q i ∩ Q k | |Q i | and |Q j ∩ Q k |
|Q j | for all other clonotypes k in the repertoire, averaged over 10000 numerical realisations.

An exact simulation for many clonotypes

We have carried out numerical simulations of the many-clonotype model, without the need for any type of approximation of the birth rates. A numerical run begins with N C T cell clonotypes, each consisting of an initial number of cells, n 0 , and |Q| APPs. The connections between these two sets are defined in terms of an N C × |Q| matrix whose elements are either 0 or 1, and equal to 1 if, and only if, a T cell clonotype receives a survival signal from an APP. During the course of the simulation, the birth rate for each clonotype is calculated at each time point from Eq. ( 1) and the per-cell death rate for T cells of clonotype i is given by μ i . In the simulations presented below, time is expressed in units of γ -1 and μ i = 0.1 for all T cell clonotypes. The first simulation satisfies

|Q i ∩ Q j | |Q i | for all (i, j).
The parameters of the simulation are N C = 15, n 0 = 50, and each clone is connected to 10 APPs (out of |Q| = 135 APPs), such that each clonotype competes with other clonotypes for access to 2 of its 10 APPs, while the other 8 APPs are not shared. The overlap between any two clonotypes consists of at most one APP. The number of T cells of a typical clonotype is plotted in Fig. 4, where the output has been averaged over 10000 individual realisations. The behaviour is consistent with the assumption of the model presented in [START_REF] Stirk | Stochastic niche structure and diversity maintenance in the T cell repertoire[END_REF] that, after a short transient, the mean number of T cells per surviving clonotype remains approximately constant.

The second simulation features a pair of clonotypes (i, j) such that

|Q i ∩Q j | ∼ |Q i | but where |Q i ∩Q k | |Q i | and |Q j ∩Q k |
|Q j | for all other clonotypes k = i, j in the repertoire. The connections between the sets of APPs and T cell clonotypes are defined such that, as before, each clonotype receives survival signals from 10 APPs. However, the pair (i, j) competes for access to 6 of the 10 APPs. Competition between all other clonotype pairs remains as in the previous simulation, and |Q| = 130. The dynamics of clonotypes other than i and j remains as in the previous case (see clonotype k in Fig. 5). The number of T cells of clonotypes i and j was reduced due to the increased competition, as one would expect. The large overlap between between Q i and Q j causes this pair to behave differently from the rest of the repertoire. This example, although constructed by hand, shows that there are circumstances in which the pair-explicit treatment of the present paper is appropriate.

In bility p = 0.025. As time increases, some clonotypes die out, so that the number of surviving clonotypes, N (t), decreases (top panel). This system has the property that, after a short transient, the mean number of T cells per surviving clonotype, n , is approximately constant (lower panel).

Guaranteed population extinction and mean extinction times

In this section we show that for the birth and death rates derived in the previous section, both clonotypes are certain to become extinct within finite time, i.e., the probability of absorption at state (n 1 , n 2 ) = (0, 0) is 1 for all values of the parameters. We also give a bound on the mean time until extinction at (0, 0) occurs from any given initial state (n 1 , n 2 ) ∈ S.

Extinction of both clonotypes occurs with probability one

The only absorbing state of the process is (n 1 , n 2 ) = (0, 0). To show that this state is reached with certainty, we follow the method of [START_REF] Iglehart | Multivariate competition processes[END_REF], bounding the bivariate process by a univariate one that moves to the origin at a slower rate than the bivariate process. We show that absorption at the origin is certain for this univariate process and conclude that the bivariate process reaches (0, 0) with probability 1.

We partition state space S as follows:

S k = {(n 1 , n 2 ) : n 1 + n 2 = k} for k ≥ 0 . ( 22
)
Let

λ k = max (n1,n2)∈Sk {λ (1) n1,n2 + λ (2) n1,n2 } , ( 23
)
μ k = min (n1,n2)∈Sk {μ (1) n1,n2 + μ (2) n1,n2 } , ( 24
)
with λ k =μ k = 0 when k = 0. For the process to pass from S k to S k it must pass through all intervening sets. The rate λ k is the maximum rate for the process to move from S k to S k+1 and the rate μ k is the minimum rate for the process to move from S k to S k-1 . These rates define a univariate birth and death process on the state space {S 0 , S 1 , S 2 , . . .} where S 0 is an absorbing state and S k is now treated as a single state rather than a set of states. The transition structure is as follows:

S 0 μ 1 S 1 λ 1 μ 2 S 2 • • • S k-1 λ k-1 μ k S k λ k μ k+1 S k+1 • • • Let π 1 = 1, π k = λ 1 λ 2 . . . λ k-1 μ 2 μ 3 . . . μ k for k ≥ 2 . ( 25
)
Then by Theorem 3 of [START_REF] Iglehart | Multivariate competition processes[END_REF], a sufficient condition for guaranteed absorption at (0, 0) is that the series

+∞ k=1 1 λ k π k (26)
A c c e p t e d m a n u s c r i p t diverges. To prove this for the rates ( 16)-( 19), we observe that

λ k = max (n1,n2)∈Sk {λ (1) n1,n2 + λ (2) n1,n2 } ≤ ϕ 1 + ϕ 2 , ( 27
)
from Eq. ( 21) and

μ k = min (n1,n2)∈Sk {μ (1) n1,n2 + μ (2) n1,n2 } = min (n1,n2)∈Sk {n 1 μ 1 + n 2 μ 2 } = k min(μ 1 , μ 2 ) . ( 28
) Then +∞ k=1 1 λ k π k = +∞ k=1 μ 2 μ 3 . . . μ k λ 1 λ 2 . . . λ k ≥ +∞ k=1 k![min(μ 1 , μ 2 )] k-1 (ϕ 1 + ϕ 2 ) k . ( 29
)
Define

a k def = k![min(μ 1 , μ 2 )] k-1 (ϕ 1 + ϕ 2 ) k , ( 30
) so that a k+1 a k = (k + 1) min(μ 1 , μ 2 ) (ϕ 1 + ϕ 2 ) → +∞ as k → +∞ . ( 31
)
Hence the series +∞ k=1 a k diverges by the ratio test and therefore +∞ k=1 1 λ k πk also diverges. Thus, absorption at (n 1 , n 2 ) = (0, 0) is certain with probability one for all parameter values of the model; the ultimate fate of both clonotypes is extinction.

A bound on the mean time until extinction

Let τ n1,n2 be the mean time until both clonotypes become extinct when the initial state is (n 1 , n 2 ). Theorem 4 of [START_REF] Iglehart | Multivariate competition processes[END_REF] states that τ n1,n2 < +∞ for all (n 1 , n 2 ) ∈ S \ {(0, 0)} if the series +∞ k=1 π k converges. The rates ( 16)-( 19) satisfy

+∞ k=1 π k = +∞ k=1 λ 1 λ 2 . . . λ k-1 μ 2 μ 3 . . . μ k ≤ +∞ k=1 (ϕ 1 + ϕ 2 ) k-1 k![min(μ 1 , μ 2 )] k-1 . ( 32
) Define b k def = (ϕ 1 + ϕ 2 ) k-1 k![min(μ 1 , μ 2 )] k-1 . ( 33
) Then b k+1 b k = ϕ 1 + ϕ 2 (k + 1) min(μ 1 , μ 2 ) → 0 as k → +∞ , ( 34
)
so that the series +∞ k=1 b k converges by the ratio test. Hence +∞ k=1 π k converges. Therefore the mean time to absorption from all initial states (n 1 , n 2 ) ∈ S \ {(0, 0)} is finite.

The bivariate competition process is bounded by the univariate birth and death process with rates λ k and μ k , in the sense that the univariate process moves towards the absorbing state at a slower rate than the bivariate competition process {(X (t), Y(t)) : t ≥ t2 }. For such a univariate birth and death process, the mean time τ m until absorption from an initial state m is [START_REF] Taylor | An Introduction to Stochastic Modelling[END_REF] where

τ m = +∞ l=1 1 λ l ρ l + m-1 j=1 ρ j +∞ k=j+1 1 λ k ρ k , ( 35)
ρ k = k j=1 (μ j /λ j ).
Hence τ m is an upper bound on the mean time to absorption at state (0, 0) from all initial states (n 1 , n 2 ) ∈ S m for m ≥ 1.

Analysis and results for the special case ν 12

1, ν 1 1, ν 2 1

In the univariate model, two special cases arise: ν 1 and ν 1, which are the "hard niche" and "soft niche" limits. For the bivariate competition process introduced in this paper, there are six such special cases, one of which will be studied in more detail here. We will focus on the case ν 12 , ν 1 , ν 2 1 since this case provides an upper bound on the mean extinction times. Two clones with ν 12 1, ν 1 1, ν 2 1 have TCRs that are different from the TCRs of other clones in the repertoire, although the two clones overlap significantly with each other in terms of the APPs that they are able to receive signals from, that is,

|Q 1 | ∼ |Q 12 | = |Q 21 | ∼ |Q 2 |.

A c c e p t e d m a n u s c r i p t

In this case, as ν 12 1, ν 1 1 and ν 2 1, the first term (r = 0) in the sums of Eqs. ( 16)-( 17) dominates and the birth rates become λ (1) 0,n2 = 0 for n 2 ≥ 0 , (36)

λ (2)
n1,0 = 0 for n 1 ≥ 0 (37)

λ (1) n1,n2 ≈ ϕ 1 1 - p 1 n 2 n 1 + n 2
for n 1 ≥ 1 and n 2 ≥ 0 , (38)

λ (2) n1,n2 ≈ ϕ 2 1 - p 2 n 1 n 1 + n 2 for n 1 ≥ 0 and n 2 ≥ 1 . ( 39
)
We focus on this particular case because the region of parameter space ν 12 1, ν 1 1, ν 2 1 is associated with the longest life-span of pairs of clones in the peripheral repertoire, as competition with T cells of other clonotypes is very small. Therefore, this case provides an "upper bound" for the mean times to extinction of the general bivariate competition model when the parameters ν 12 , ν 1 and ν 2 are allowed to take arbitrary values. Furthermore, the number of parameters in the model is reduced from nine to five (viz. ϕ 1 , ϕ 2 , p 1 , μ 1 , μ 2 ). In the rest of this section we will assume μ 1 = μ = μ 2 and μ = 1, i.e., time is measured in units of μ -1 .

Extinction times in the special case ν 12

1, ν 1 1, ν 2 1

Before extinction at (0, 0), the bivariate competition process enters the absorbing set A = {(n 1 , n 2 ) : n 1 = 0 or n 2 = 0} when one of the T cell clonotypes becomes extinct, because in order for the process to enter the state (0, 0), a transition from either state (0, 1) or state (1, 0) must occur. Let τn1,n2 be the mean time until the process enters the absorbing set A when the initial state of the process is given by (n 1 , n 2 ) ∈ S \ A. One can show by first-step analysis [START_REF] Taylor | An Introduction to Stochastic Modelling[END_REF]) that this quantity satisfies the following two-dimensional difference equation (where we write λ

(1)

n1,n2 + λ (2) n1,n2 + μ (1) n1,n2 + μ (2) n1,n2 = α n1,n2 ): τn1,n2 = λ (1) n1,n2 α n1,n2 τn1+1,n2 + λ (2) n1,n2 α n1,n2 τn1,n2+1 + μ (1) n1,n2 α n1,n2 τn1-1,n2 + μ (2) n1,n2 α n1,n2 τn1,n2-1 + 1 α n1,n2 , (40)
with the boundary conditions τj,0 = τ0,j = 0 ∀j ≥ 1, τN+1,j = τj,N+1 = 0 ∀j ≥ 0 where the state space of the process is truncated to be finite, i.e., S = {(n 1 , n 2 ) : n 1 , n 2 = 0, 1, 2, . . . , N} to allow numerical computation. This set of equations can be written in the form Aτ = b, where A is an N 2 × N 2 matrix and τ , b ∈ R N 2 (see Appendix B for further details). Figs. 7-9 show that the time to reach the absorbing set decreases as p 1 increases. In Fig. 7, τn1,n2 is plotted as a function of p 1 for various initial states with n 1 + n 2 = 10 (the initial total number of cells is fixed). In Fig. 8, τn1,n2 is plotted as a function of p 1 for different initial total number of cells with n 1 = n 2 . These plots show that τn1,n2 depends on the initial state of the process and increases as the initial state lies further away from the absorbing set A = {(n 1 , n 2 ) : n 1 = 0 or n 2 = 0}. In Fig. 9, τn1,n2 is plotted as a function of p 1 for several values of ϕ 1 , with all other parameters and the initial conditions fixed. Fig. 10 shows that the time during which both clones are extant in the repertoire, τn1,n2 , compared to the time until both clones become extinct, τ n1,n2 , decreases as p 1 increases. Therefore, as p 1 → 1 -, one of the two clonotypes quickly becomes extinct by a process that resembles the competitive exclusion principle of classical ecology, which states that two species competing for the same set of resources cannot stably coexist [START_REF] Begon | Ecology: Individuals, Populations and Communities. Blackwell's Sci[END_REF]. The numerical results indicate that the quantity τn1,n2 /τ n1,n2 is maximal when both clones have access to the same number of APPs from which they can receive survival stimuli, i.e., ϕ 1 = ϕ 2 , and decreases as one clone gains a competitive advantage over the other, i.e., ϕ 1 > ϕ 2 or vice versa. Also of interest is the probability ℘ n1,n2 that clonotype 1 becomes extinct before clonotype 2 when the process starts from (n 1 , n 2 ) ∈ S \ A. The probability of clonotype 2 becoming extinct before clonotype 1 is given by 1 -℘ n1,n2 , because extinction is guaranteed. The probability ℘ n1,n2 satisfies the difference equation

℘ n1,n2 = λ (1) n1,n2 α n1,n2 ℘ n1+1,n2 + λ (2) n1,n2 α n1,n2 ℘ n1,n2+1 + μ (1) n1,n2 α n1,n2 ℘ n1-1,n2 + μ (2) n1,n2 α n1,n2 ℘ n1,n2-1 , ( 41
)
with the boundary conditions ℘ 0,j = 1 ∀j ≥ 1, ℘ j,0 = 0 ∀j ≥ 1, ℘ N +1,j = ℘ j,N +1 = 0 ∀j ≥ 0, where we have truncated the state space to allow numerical computation. In Fig. 11 the probability that clonotype 1 becomes extinct before clonotype 2 is plotted as a function of p 1 for different values of ϕ 1 . This probability decreases as ϕ 1 increases, due to an increasing survival stimulus for T cells of clonotype 1, giving a competitive advantage over T cells of clonotype 2. This quantity depends on the initial state, as illustrated in Fig. 12. As p 1 increases, so does the effect of the advantage gained from a higher initial number of cells than its competitor. In the case ϕ 1 = ϕ 2 , p 1 = p 2 , μ 1 = μ 2 where n 1 = n 2 , neither clone has a competitive advantage, as shown by the straight line in Fig. 12. However, changes in the initial conditions such that n 1 = n 2 break this symmetry. 1, ν 1 1, ν 2 1 with ϕ 1 = 5, 8, 10, 20, ϕ 2 = 10, μ 1 = μ 2 = 1 from the initial state n 1 = n 2 = 10. Note that for ϕ 1 = 20, ϕ 2 = 10 we require that p 1 ≤ 0.5 so that the condition p 2 ≤ 1 is satisfied because of the constraint ϕ 1 p 1 = ϕ 2 p 2 . Homeostatic (constant) levels of naïve T cells are maintained by clonal competition for self-antigens presented on the surface of APCs [START_REF] Troy | Cutting edge: homeostatic proliferation of peripheral T lymphocytes is regulated by clonal competition[END_REF]). Yet, clonal sizes exhibit a certain statistical distribution under the influence of the homeostatic drive we have been considering. This probability distribution corresponds to the limiting conditional probability distribution (LCD) of the process. The idea is that, before a clonotype goes extinct, its clonal size wanders stochastically according to a probability distribution which is stationary for most of the time before the absorbing set is reached. The LCD is of interest because it may be expected to reflect the distribution of clonal sizes found at any given instant (the latter being experimentally accessible, in principle).

Let p n1,n2 (t) be the probability there are n 1 T cells of clonotype 1 and n 2 T cells of clonotype 2 at time t:

p n1,n2 (t) = P(X (t) = n 1 , Y(t) = n 2 |X ( t2 ) = ñ1 , Y( t2 ) = ñ2 ) , ( 42
)
with +∞ n1=0 +∞ n2=0 p n1,n2 (t) = 1. These probabilities satisfy the following system of differential equations:

dp n1,n2 (t) dt = λ (1) n1-1,n2 p n1-1,n2 (t) + λ (2) n1,n2-1 p n1,n2-1 (t) + μ (1) n1+1,n2 p n1+1,n2 (t) + μ (2) n1,n2+1 p n1,n2+1 (t) -(λ (1) n1,n2 + λ (2) n1,n2 + μ (1) n1,n2 + μ (2) n1,n2 )p n1,n2 (t) , ( 43
)
for (n 1 , n 2 ) ∈ S. We now condition on the event that extinction of either clonotype has not yet occurred. The absorbing set is A = {(n 1 , n 2 ) : n 1 = 0 or n 2 = 0} and p Ā(t) is the probability that the process is not in set A at time t. For all (n 1 , n 2 ) ∈ S \ A we introduce the probability that the process is in (n 1 , n 2 ) at time t, conditioned on the event that it has not yet entered the absorbing set:

q n1,n2 (t) def = p n1,n2 (t) p Ā(t) . ( 44
)
These probabilities satisfy

+∞ n1=1 +∞ n2=1 q n1,n2 (t) = 1, q n1,n2 (t) ≥ 0 for (n 1 , n 2 ) ∈ S \ A and q n1,n2 (t) = 0 for (n 1 , n 2 ) ∈ A.
We have conditioned on the event that the process has not yet reached the absorbing set, because after one clonotype goes extinct the system can be modelled via the univariate process. From Eq. ( 44) we have

dq n1,n2 (t) dt = 1 p Ā(t) dp n1,n2 (t) dt - p n1,n2 (t) (p Ā(t)) 2 dp Ā(t) dt . ( 45
)
By the law of total probability:

p Ā(t) = 1 - +∞ n2=0 p 0,n2 (t) - +∞ n1=0 p n1,0 (t) + p 0,0 (t) . ( 46
)
Substituting n 1 = 0 in Eq. ( 43), summing over n 2 , then substituting n 2 = 0 in Eq. ( 43) and summing over n 1 gives

d dt +∞ n2=0 p 0,n2 (t) = μ 1 +∞ n2=0 p 1,n2 (t) , ( 47
)
d dt +∞ n1=0 p n1,0 (t) = μ 2 +∞ n1=0 p n1,1 (t) , ( 48
)
and substituting n 1 = n 2 = 0 in Eq. ( 43) yields

d dt p 0,0 (t) = μ 1 p 1,0 (t) + μ 2 p 0,1 (t) , ( 49
) so that dp Ā(t) dt = -μ 1 +∞ n2=1 p 1,n2 (t) -μ 2 +∞ n1=1 p n1,1 (t) . ( 50
)
Hence, q n1,n2 (t) obeys

dq n1,n2 (t) dt = λ (1) n1-1,n2 q n1-1,n2 (t) + λ (2) n1,n2-1 q n1,n2-1 (t) + μ (1) n1+1,n2 q n1+1,n2 (t) + μ (2) n1,n2+1 q n1,n2+1 (t) -(λ (1) n1,n2 + λ (2) n1,n2 + μ (1) n1,n2 + μ (2) n1,n2 )q n1,n2 (t) + μ 1 q n1,n2 (t) +∞ n2=1 q 1,n2 (t) + μ 2 q n1,n2 (t) +∞ n1=1 q n1,1 (t) . ( 51
1, ν 1 1, ν 2 1 with ϕ 1 = ϕ 2 = 10, p 1 = p 2 = 0.5, μ 1 = μ 2 = 1.
A probability distribution q, if it exists, is called a quasi-stationary probability distribution (QSD) of the process if it satisfies

0 = λ (1) n1-1,n2 qn1-1,n2 + λ (2) n1,n2-1 qn1,n2-1 + μ (1) n1+1,n2 qn1+1,n2 + μ (2) n1,n2+1 qn1,n2+1 -(λ (1) n1,n2 + λ (2) n1,n2 + μ (1) n1,n2 + μ (2) n1,n2 )q n1,n2 + μ 1 qn1,n2 +∞ n2=1 q1,n2 + μ 2 qn1,n2 +∞ n1=1 qn1,1 , (52)
where +∞ n1=1 +∞ n2=1 qn1,n2 = 1, qn1,n2 ≥ 0 for (n 1 , n 2 ) ∈ S \ A and qn1,n2 = 0 for (n 1 , n 2 ) ∈ A. The limiting conditional probability distribution (LCD) of the process is the limit of the conditional probability distribution q n1,n2 (t) as t → +∞. If it exists, this distribution must be a QSD. For a process with finite state space, there exists a unique QSD which is also the LCD of the process [START_REF] Darroch | On quasi-stationary distributions in absorbing continuous time finite Markov chains[END_REF]. A method of approximating the LCD will be described in the next section. Figs. 13-14 show the LCD of the process for two sets of the parameters in the special case ν 12 1, ν 1 1, ν 2 1. If the state space of the process is denumerably infinite, a QSD might not exist [START_REF] Clancy | Approximations for the long-term behaviour of an open-population epidemic model[END_REF], or if it does exist, it may not be unique. Then the initial probability distribution of the process may affect which QSD gives the LCD. Questions regarding the existence of a QSD and an LCD when the state space of the process is infinite are outside the scope of this paper and will be the subject of a future publication.

6 The large N expansion for the special case ν 12 1, ν 1 1, ν 2 1

Van Kampen's "large N expansion" [START_REF] Van Kampen | A power series expansion of the master equation[END_REF]van Kampen, 2007) is a method to construct a continuous approximation to the discrete stochastic model. Here, we study the fluctuations about the stable steady state for the special case ν 12 1, ν 1 1, ν 2 1. In the large N approximation, continuous variables x i (t) and ξ i (t) are introduced and n i is written as a mean plus fluctuations:

n i (t) = Ωx i (t) + Ω 1 2 ξ i (t) (53)
(for i = 1, 2) where Ωx i (t) = IE(n i (t)). The relative size of the fluctuations ξ i (t) is posited to be of order Ω 1 2 , and the mean of order Ω, with the parameter Ω representing the (large) size of the system. The set of probabilities p n1,n2 (t) and the forward Kolmogorov equation ( 43

1, ν 1 1, ν 2 1 with ϕ 1 = 5, ϕ 2 = 10, p 1 = 0.5, p 2 = 0.25, μ 1 = μ 2 = 1. Π(ξ 1 , ξ 2 , t) that satisfies ∂Π ∂t -Ω 1 2 dx 1 dt ∂Π ∂ξ 1 -Ω 1 2 dx 2 dt ∂Π ∂ξ 2 = -Ω -1 2 ∂ ∂ξ 1 + 1 2 Ω -1 ∂ 2 ∂ξ 2 1 + . . . λ (1) Ωx1+Ω 1 2 ξ1,Ωx2+Ω 1 2 ξ2 Π + -Ω -1 2 ∂ ∂ξ 2 + 1 2 Ω -1 ∂ 2 ∂ξ 2 2 + . . . λ (2) Ωx1+Ω 1 2 ξ1,Ωx2+Ω 1 2 ξ2 Π + Ω -1 2 ∂ ∂ξ 1 + 1 2 Ω -1 ∂ 2 ∂ξ 2 1 + . . . μ (1) Ωx1+Ω 1 2 ξ1,Ωx2+Ω 1 2 ξ2 Π + Ω -1 2 ∂ ∂ξ 2 + 1 2 Ω -1 ∂ 2 ∂ξ 2 2 + . . . μ (2) Ωx1+Ω 1 2 ξ1,Ωx2+Ω 1 2 ξ2 Π. (54)
In the special case

ν 12 1, ν 1 1, ν 2 1, we have ∂Π ∂t -Ω 1 2 dx 1 dt ∂Π ∂ξ 1 -Ω 1 2 dx 2 dt ∂Π ∂ξ 2 = -Ω -1 2 ∂ ∂ξ 1 + 1 2 Ω -1 ∂ 2 ∂ξ 2 1 + . . . Ω φ1 x 1 + x 2 x 1 + x 2 -p 1 x 2 + Ω -1 2 p 1 (ξ 1 x 2 -ξ 2 x 1 ) x 1 + x 2 Π + -Ω -1 2 ∂ ∂ξ 2 + 1 2 Ω -1 ∂ 2 ∂ξ 2 2 + . . . Ω φ2 x 1 + x 2 x 1 + x 2 -p 2 x 1 + Ω -1 2 p 2 (ξ 2 x 1 -ξ 1 x 2 ) x 1 + x 2 Π + Ω -1 2 ∂ ∂ξ 1 + 1 2 Ω -1 ∂ 2 ∂ξ 2 1 + . . . μ 1 (Ωx 1 + Ω 1 2 ξ 1 )Π + Ω -1 2 ∂ ∂ξ 2 + 1 2 Ω -1 ∂ 2 ∂ξ 2 2 + . . . μ 2 (Ωx 2 + Ω 1 2 ξ 2 )Π . (55)

The deterministic limit

The variables x i are described by deterministic differential equations, obtained by collecting terms of order Ω 1 2 from Eq. ( 55), or directly from the rates (18), ( 19), ( 39) and (39). For the case ν 12 1, ν 1 1, ν 2 1, the deterministic equations are where Ω φ1 = ϕ 1 and Ω φ2 = ϕ 2 . The number of parameters in the model can be reduced by non-dimensionalising the above equations and using the constraint on the parameters given by Eq. ( 20):

dx 1 dt = φ1 1 - p 1 x 2 x 1 + x 2 -μ 1 x 1 , ( 56
)
dx 2 dt = φ2 1 - p 2 x 1 x 1 + x 2 -μ 2 x 2 , (57)
τ = μ 1 t , ( 58
)
x i = φ1 μ 1 u i for i = 1, 2. ( 59
)
The scaled equations are:

du 1 dτ = 1 - p 1 u 2 u 1 + u 2 -u 1 , ( 60
)
du 2 dτ = φ - p 1 u 1 u 1 + u 2 -μu 2 , ( 61
)
where φ = ϕ 2 /ϕ 1 and μ = μ 2 /μ 1 . It can be shown that when ν 12 1, ν 1 1, ν 2 1 there exists a unique steady state (ū 1 , ū2 ) (Appendix C), where 1 -p 1 < ū1 < 1 and 1 μ ( φ -p 1 ) < ū2 < φ μ , which is locally asymptotically stable for all values of the parameters (Appendix D).

If we assume that the per-cell death rates for both clonotypes are equal, i.e., μ 1 = μ 2 then, in terms of the original variables n 1 and n 2 , the unique steady state is given by:

n1 = Ωx 1 = ϕ 1 μ 1 (1 -p 1 + ϕ2 ϕ1 )(1 -p 1 ) 1 -2p 1 + ϕ2 ϕ1 , ( 62
) n2 = Ωx 2 = ϕ 1 μ 1 (1 -p 1 + ϕ2 ϕ1 )( ϕ2 ϕ1 -p 1 ) 1 -2p 1 + ϕ2 ϕ1 . ( 63
)
If ϕ 2 > ϕ 1 , the number of T cells of clonotype 2 at the steady state is greater than the number of T cells of clonotype 1, since T cells of clonotype 2 have access to a larger set of APPs from which they are able to receive survival stimuli. Extinction does not happen in the deterministic approximation, which underscores the need for a stochastic treatment.

If the mean time to extinction is large, the maximum of the corresponding LCD can be approximated by the deterministic steady state. For the parameter values ϕ 1 = ϕ 2 = 10, p 1 = p 2 = 0.5 and μ 1 = μ 2 = 1, the steady state is given by (7 1 2 , 7 1 2 ), which corresponds closely to the maximum of the LCD at state (7, 7) as shown in Fig. 13, while for ϕ 1 = 5, ϕ 2 = 10, p 1 = p 2 = 0.5, μ 1 = μ 2 = 1, the steady state is given by (3 1 8 , 9 3 8 ) and the maximum of the LCD is (3, 9), as shown in Fig. 14.

Fluctuations about the steady state and an approximation to the limiting conditional distribution

We now study the fluctuations about the deterministic stable steady state. This provides a method of approximating the limiting conditional probability distribution of the process.

Collecting terms of order Ω 0 from Eq. ( 55) results in the equation

∂Π ∂t = 1 2 φ1 x 1 + x 2 (x 1 + x 2 -p 1 x 2 ) ∂ 2 Π ∂ξ 2 1 + 1 2 φ2 x 1 + x 2 (x 1 + x 2 -p 2 x 1 ) ∂ 2 Π ∂ξ 2 2 - φ1 p 1 (x 1 + x 2 ) 2 ∂ ∂ξ 1 [(ξ 1 x 2 -ξ 2 x 1 )Π] - φ2 p 2 (x 1 + x 2 ) 2 ∂ ∂ξ 2 [(ξ 2 x 1 -ξ 1 x 2 )Π] + 1 2 μ 1 x 1 ∂ 2 Π ∂ξ 2 1 + μ 1 ∂ ∂ξ 1 (ξ 1 Π) + 1 2 μ 2 x 2 ∂ 2 Π ∂ξ 2 2 + μ 2 ∂ ∂ξ 2 (ξ 2 Π) , ( 64
)
which is a linear bivariate Fokker-Planck equation for the probability distribution of the fluctuations, Π(ξ 1 , ξ 2 , t), the solution of which is an Ornstein-Uhlenbeck process (van Kampen, 2007) and is therefore fully determined by the first and second moments. Multiplying Eq. ( 64) by ξ 1 and integrating over all values of ξ 1 and ξ 2 results in the differential equation

d dt ξ 1 = φ1 p 1 x 2 (x 1 + x 2 ) 2 -μ 1 ξ 1 - φ1 p 1 x 1 (x 1 + x 2 ) 2 ξ 2 .
(65)

A c c e p t e d m a n u s c r i p t

Similarly we obtain

d dt ξ 2 = - φ2 p 2 x 2 (x 1 + x 2 ) 2 ξ 1 + φ2 p 2 x 1 (x 1 + z 2 ) 2 -μ 2 ξ 2 , ( 66
)
d dt ξ 2 1 = 2 φ1 p 1 x 2 (x 1 + x 2 ) 2 -μ 1 ξ 2 1 - 2 φ1 p 1 x 1 (x 1 + x 2 ) 2 ξ 1 ξ 2 + φ1 x 1 + x 2 (x 1 + x 2 -p 1 x 2 ) + μ 1 x 1 , ( 67
)
d dt ξ 2 2 = 2 φ2 p 2 x 1 (x 1 + x 2 ) 2 -μ 2 ξ 2 2 - 2 φ2 p 2 x 2 (x 1 + x 2 ) 2 ξ 1 ξ 2 + φ2 x 1 + x 2 (x 1 + x 2 -p 2 x 1 ) + μ 2 x 2 , ( 68
)
d dt ξ 1 ξ 2 = - φ2 p 2 x 2 (x 1 + x 2 ) 2 ξ 2 1 - φ1 p 1 x 1 (x 1 + x 2 ) 2 ξ 2 2 + φ1 p 1 x 1 + x 2 -μ 1 -μ 2 ξ 1 ξ 2 . ( 69
)
6.2.1 Fluctuations about the steady state in the symmetric case

ϕ 1 = ϕ 2 , p 1 = p 2 , μ 1 = μ 2
We now substitute into Eqs. ( 65)-( 69) the values of x 1 and x 2 at the stable steady state for the symmetric case ϕ 1 = ϕ 2 , p 1 = p 2 and μ 1 = μ 2 , as the algebra is somewhat simpler than in the more general situation. The steady state values are given by

x1 = x2 = φ1 (2 -p 1 ) 2μ 1 . ( 70
)
We integrate Eqs. ( 65)-( 66) with initial conditions ξ 1 (0) = ξ 1 0 , ξ 2 (0) = ξ 2 0 , to find

ξ 1 (t) = ( ξ 1 0 + ξ 2 0 ) 2 e -μ1t + ( ξ 1 0 -ξ 2 0 ) 2 e -2μ 1 2-p 1 t , ( 71
)
ξ 2 (t) = ( ξ 1 0 + ξ 2 0 ) 2 e -μ1t - ( ξ 1 0 -ξ 2 0 ) 2 e -2μ 1 2-p 1 t , ( 72
)
so that the (stable) stationary values for the means of the fluctuations are ξ 1 s = ξ 2 s = 0, as described by Eqs. ( 71)-( 72). The solutions of Eqs. ( 67)-( 69) are given by

ξ 2 1 (t) = φ1 (2 -p 1 )(3p 1 -4) 8μ 1 (p 1 -1) + c 1 e -2μ1t + c 2 e μ 1 (3p 1 -4) 2-p 1 t + c 3 e 4μ 1 (p 1 -1) 2-p 1 t , ( 73
) ξ 2 2 (t) = φ1 (2 -p 1 )(3p 1 -4) 8μ 1 (p 1 -1) + c 1 e -2μ1t -c 2 e μ 1 (3p 1 -4) 2-p 1 t + c 3 e 4μ 1 (p 1 -1) 2-p 1 t , ( 74
) ξ 1 (t)ξ 2 (t) = φ1 p 1 (2 -p 1 ) 8μ 1 (p 1 -1) + c 1 e -2μ1t -c 3 e 4μ 1 (p 1 -1) 2-p 1 t , ( 75
)
where

c 1 = 1 4 [ ξ 2 1 0 + ξ 2 2 0 ] + 1 2 ξ 1 ξ 2 0 - φ1 (2 -p 1 ) 4μ 1 , ( 76
) c 2 = 1 2 [ ξ 2 1 0 -ξ 2 2 0 ] , ( 77
) c 3 = 1 4 [ ξ 2 1 0 + ξ 2 2 0 ] - 1 2 ξ 1 ξ 2 0 + φ1 (2 -p 1 ) 2 8μ 1 (p 1 -1) , ( 78
)
and ξ 2 1 (0) = ξ 2 1 0 , ξ 2 2 (0) = ξ 2 2 0 , and ξ 1 (0)ξ 2 (0) = ξ 1 ξ 2 0 . It follows that the (stable) stationary values of the second moments are given by

ξ 2 1 s = ξ 2 2 s = φ1 (2 -p 1 )(3p 1 -4) 8μ 1 (p 1 -1) ≥ 0 , ( 79
)
and

ξ 1 ξ 2 s = φ1 p 1 (2 -p 1 ) 8μ 1 (p 1 -1) ≤ 0 , ( 80
)
as described by Eqs. ( 73)-(75). Hence, the LCD of the competition process may be approximated by a bivariate normal distribution with mean

Ω(x 1 , x2 ) = (n 1 , n2 ) = ϕ 1 (2 -p 1 ) 2μ 1 , ϕ 1 (2 -p 1 ) 2μ 1 , ( 81
) and covariance matrix Ω ξ 2 1 s ξ 1 ξ 2 s ξ 1 ξ 2 s ξ 2 2 s = 3p1-4 4(p1-1) n1 p1 4(p1-1) n1 p1 4(p1-1) n1 3p1-4 4(p1-1) n1 . ( 82
1, ν 1 1, ν 2 1 with ϕ 1 = ϕ 2 , p 1 = p 2 , μ 1 = μ 2 . 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 0 0.2 0.4 0.6 0.8 1 √ Ω ξ 2 1 s n1 p1
1, ν 1 1, ν 2 1, ϕ 1 = ϕ 2 = 10, μ 1 = μ 2 = 1.
For this to be a good approximation we require that the stationary point (n 1 , n2 ) is stable and that both n1 and n2 are large enough so that it is unlikely that the Ornstein-Uhlenbeck process will reach the absorbing boundary at either n1 = 0 or n2 = 0 (see Appendix E).

The Pearson product moment correlation coefficient [START_REF] Grimmett | Probability and random processes[END_REF] between the number of T cells of the two clonotypes, n 1 and n 2 , at the LCD is given by

ρ def = ξ 1 ξ 2 s ξ 2 1 s ξ 2 2 s = p 1 3p 1 -4 , ( 83
)
from Eq. ( 82). As would be expected for two competing populations, this is always negative and, moreover, ρ → -1 + as the two clonotypes overlap more completely with each other, (i.e., p 1 → 1 -), as shown in Fig. 15. The coefficient of variation (CV) for the number of T cells belonging to a particular clonotype at the LCD is a dimensionless measure of the dispersion of the marginal distribution and is given by its standard deviation, Ω ξ 2 1 s , divided by the mean of the distribution, n1 . The CV for T cells of clonotype 1, at the LCD, can be approximated by

Ω ξ 2 1 s n1 = μ 1 (3p 1 -4) 2ϕ 1 (p 1 -1)(2 -p 1 ) . ( 84
)
As p 1 → 0 + , the CV tends to μ 1 /ϕ 1 , in agreement with earlier results from analysis of the univariate process, and as p 1 → 1 -the CV tends to +∞ indicating the breakdown of the deterministic approximation. In this case, numerical results from Section 4.1 suggest that one clonotype would quickly outcompete the other before the process converged to the LCD, because the mean time to reach the absorbing set is short. The deterministic approximation is poor and the large N expansion is not valid. This is illustrated in Fig. 16, which shows the CV for the number of T cells belonging to clonotype 1 at the numerical LCD (solid line) and the CV derived from the normal approximation defined above (dotted line); see Eq. ( 84).

A c c e p t e d m a n u s c r i p t 7 Discussion

Clonal competition is fundamental to T cell homeostasis and the maintenance of diversity [START_REF] Troy | Cutting edge: homeostatic proliferation of peripheral T lymphocytes is regulated by clonal competition[END_REF][START_REF] Stirk | Stochastic niche structure and diversity maintenance in the T cell repertoire[END_REF]. Numerical Monte Carlo simulations show that two clonotypes whose TCRs are very similar cannot be treated as if they exert a mean-field influence on each other, which motivates and justifies the bivariate analysis which we have pursued here.

We have shown that, as the proportion of APPs providing survival stimuli which are shared by both clonotypes increases, the time until one clonotype becomes extinct decreases -resembling the ecological principle of classical competitive exclusion more closely as p 1 → 1 -. We have defined the LCD to represent the stationary behaviour of the process before extinction occurs and used the large N expansion as a means of approximating the LCD when the mean time until extinction of both clones is long. Together with the results of the previous paper, these results introduce a new hypothesis of TCR diversity maintenance: clones vanish at a rate that is continuously dependent on their similarity to other clones, maximising the evenness with which the TCR repertoire covers the space of possible antigens.

The large N approximation can be carried out (along the same lines as we have described here) in the most general case to derive (i) deterministic equations for the number of cells belonging to each clonotype and (ii) a linear Fokker-Planck equation for the fluctuations about the deterministic steady state, in order to approximate the limiting conditional probability distribution. Conversely, deterministic models of the TCR repertoire (e.g., de Boer andPerelson, 1994, 1997) can be viewed and justified as the large N component of a stochastic treatment, which must always be regarded as the more fundamental perspective. Moreover, as we have shown above, we can characterise conditions for which the deterministic approximation breaks down.

More importantly, the stochastic component contains valuable information which cannot be obtained following a deterministic approach. In particular, the correlation of the time series of clonal numbers n i and n j gives a measure of the functional overlap between the TCR clonotypes i and j. Such data can be obtained by means of repeated RT-PCR on peripheral blood samples (e.g., [START_REF] Bousso | The composition of a primary T cell response is largely determined by the timing of recruitment of individual T cell clones[END_REF]. Thus our method gives the immunologist a new tool to characterise functional similarity of TCR clonotypes and to analyse the dynamics of the TCR repertoire.

TCR diversity has two important aspects: dispersal of TCRs over the space of antigens, and evenness of numbers (i.e., little variation of n i from one clone to the next). The latter type of evenness is measured by the numbers-equivalent diversity of order q ≥ 0, defined as

D(q) def = NC i=1 n i N C n q 1/(1-q)
q = 1 , (85) [START_REF] Jost | Partitioning diversity into independent alpha and beta components[END_REF] where, for example, D(0) = N C , D(1) is the exponential of the Shannon entropy

D(1) def = lim q→1 D(q) = exp - NC i=1 n i N C n log n i N C n , ( 86
)
and D( 2) is Simpson's diversity index [START_REF] Jost | Partitioning diversity into independent alpha and beta components[END_REF]. The D(q) diversity spectrum, obtained by varying q, probes evenness of clonal numbers across the repertoire; this spectrum is closely linked to the ability of the immune system to mount an effective and timely response (e.g., exceedance probabilities for the time required by the system to find and activate a T cell clone of the required quality are an exponential function of D(2); [START_REF] Stirk | Stochastic niche structure and diversity maintenance in the T cell repertoire[END_REF]. Crucially, the other aspect of diversity, i.e., dispersal over antigen space, is not tied to the size or fate of any one particular clone. This dispersal is driven by competition for survival stimuli: the mean life time of any given clone is a sharply decreasing function of the clone's recognition overlap with other clonotypes extant in the repertoire [START_REF] Stirk | Stochastic niche structure and diversity maintenance in the T cell repertoire[END_REF]. For a pair of clonotypes (n i (t), n j (t)) that are perfectly overlapping as regards their TCR recognition profile, it is sufficient, as far as maximisation of coverage diversity is concerned, if the pair as an entity behaves in essentially the same way as a single clonotype does in the mean field limit. Then coverage diversity is affected neither by the anti-correlated fluctuations in n i (t) and n j (t), nor by the identity of the first of the pair to go extinct (or the uncertain timing of this event). The pair-approximation analysis in the present paper bears this out: as p 1 approaches unity, the pair, considered as the combined entity n i (t) + n j (t), has the same stochastic dynamics as a single mean field clone (see Appendix A); on the other hand, as p 1 approaches zero, the pair behaves as a pair of mean field entities that are independent apart from their coupling through the "Brownian background" field (provided by the potentially N C -2 other peripheral T cell clonotypes). There is a gradual transition between these two extremes as p 1 assumes intermediate values.

Figure 2 :

 2 Figure2: Schematic diagram of the survival stimuli model. On the left each circle represents an APP, while each circle on the right represents a T cell clonotype (all the T cells with an identical TCR). An edge from an APP q to a T cell clonotype i indicates that APP q can provide a survival signal to T cells of this clonotype.

  , is the rate of transition from state (n 1 , n 2 )

Figure 3 :

 3 Figure 3: A representation of the bivariate competition process and the transitions between different states.

Figure 4 :Figure 5 :

 45 Figure 4: The number of T cells belonging to a typical clonotype for a repertoire when |Q i ∩ Q j | |Q i | for all pairs of clonotypes averaged over 10000 individual numerical realisations.

Figure 6 :

 6 Figure 6: Time series from a single numerical solution with randomly-assigned recognition of 400 APPs by 100 T cell clonotypes. The parameters used were N C = 100, |Q| = 400, n 0 = 50, μ = 0.1 γ = 1.0, p = 0.025.

Figure 7 :Figure 8 :Figure 9 :

 789 Figure7: The quantity τn1,n2 as a function of p 1 for the special case ν 12 1, ν 1 1, ν 2 1 with ϕ 1 = 5, ϕ 2 = 10, μ 1 = μ 2 = 1 for various initial states with n 1 + n 2 = 10.

Figure 10 :Figure 11 :Figure 12 :

 101112 Figure10: The quantity τn1,n2 /τ n1,n2 as a function of p 1 for the special case ν 12 1, ν 1 1, ν 2 1 with ϕ 1 = 5, 8, 10, 20, ϕ 2 = 10, μ 1 = μ 2 = 1 from the initial state n 1 = n 2 = 10.

Figure 13 :

 13 Figure13: The limiting conditional probability distribution in the special caseν 12 1, ν 1 1, ν 2 1 with ϕ 1 = ϕ 2 = 10, p 1 = p 2 = 0.5, μ 1 = μ 2 = 1.

Figure 14 :

 14 Figure 14: The limiting conditional probability distribution in the special case ν 121, ν 1 1, ν 2 1 with ϕ 1 = 5, ϕ 2 = 10, p 1 = 0.5, p 2 = 0.25, μ 1 = μ 2 = 1.

Figure 15 :

 15 Figure 15: Correlation between clonotype 1 and clonotype 2 for the case ν 121, ν 1 1, ν 2 1 with ϕ 1 = ϕ 2 , p 1 = p 2 , μ 1 = μ 2 .

Figure 16 :

 16 Figure16: The coefficient of variation for the number of T cells belonging to clonotype 1 at the LCD withν 12 1, ν 1 1, ν 2 1, ϕ 1 = ϕ 2 = 10, μ 1 = μ 2 = 1.

ΑFigure 17 :

 17 Figure17: Probability inside the ellipsoid that first touches the positive semi-axes (α) as a function of p for μ = 1 and ϕ = 1 (bottom line), 5, 10, 100, 1000 (top line) (left panel). Probability inside the ellipsoid that first touches the positive semi-axes (α) as a function of p for ϕ = 1 and μ = 1 (top line), 5, 10, 100, 1000 (bottom line) (right panel).

Table 1

 1

	: Notation

A c c e p t e d m a n u s c r i p t

	and
			|Q 1r/2 | |Q 1/2
							)
	Since N C	1 and p .|12	1, p .|1/2	1, a Poisson approximation is warranted. Defining ν 12 = (N C -2)p .|12
	and ν 1 = (N C -2)p .|1/2 , we obtain
			|Q 12r | |Q 12 |	ν r 12 e -ν12 r!	= p 1 |Q 1 |	ν r 12 e -ν12 r!	,	(9)

Acknowledgements

ERS was supported by an EPSRC Doctoral Training Grant and the University of Leeds. CM-P was supported by a BBSRC Research Development Fellowship. We thank Dr. Shev MacNamara for a careful reading of the manuscript and for useful discussions and Dr. Stuart Barber for introducing us to the Mahalanobis distance.

Appendix A. The limits p 1 = p 2 = 0 and p 1 = p 2 = 1

In this section we will show that for a pair of clones with ϕ 1 = ϕ 2 , p 1 = p 2 and μ 1 = μ 2 , there are two limits in which the bivariate competition process can be simplified to a univariate process. First, let p 1 = p 2 = 0. Then the birth rates ( 16)-( 17) become λ (i) n1,n2 = ϕ i n i e -νi +∞ r=0 ν i r!

for i = 1, 2, which are of the form of the birth rates for a single clonotype in the univariate model (cf. [START_REF] Stirk | Stochastic niche structure and diversity maintenance in the T cell repertoire[END_REF]). This is because, in this case, Q 1 ∩Q 2 = ∅, and therefore the two clones have independent dynamics. Now let p 1 = p 2 = 1, which means that the sets of APPs from which each clone receives survival stimuli overlap completely, i.e.,

Then the birth rates ( 16)-( 17) become

so that the pair of clones together behave in the same way as a single clonotype in the univariate model. For intermediate values of p 1 and p 2 (0 < p i < 1; i = 1, 2), the bivariate analysis presented in this paper is required (the effect of varying the parameter p 1 is demonstrated in Figs. 7-12).

Appendix B. Numerical solution of the two-dimensional difference equations

Eq. ( 40) can be rearranged to give

A c c e p t e d m a n u s c r i p t

and

The matrix equation (B-2) can be solved numerically using the MATLAB package [START_REF] Guide | MATLAB[END_REF]. The numerical solution of Eq. ( 41) follows a similar method. (The MATLAB code used to solve this matrix equation is available upon request).

Appendix C. Existence of a unique steady state for the deterministic model in the case ν 12 1, ν 1 1, ν 2 1

Steady states are found by setting du1 dt = du2 dt = 0 in Eqs. ( 60)-( 61):

Since (0, 0) is not a solution, the u 1 nullcline is given by [START_REF] Edelstein-Keshet | Mathematical models in biology[END_REF])

and, similarly, the u 2 nullcline is given by

Steady states occur at the intersection of the nullclines. Since we are concerned with a biological system, we require that u 1 , u 2 ≥ 0. Hence 1 -p 1 < u 1 < 1 and 1 μ ( φ -p 1 ) < u 2 < φ μ . This leads us to define the following functions:

Both f 1 (u 1 ) and f 2 (u 2 ) are decreasing and continuous on their domains with

These functions must intersect and therefore a unique steady state exists. We now prove this. Since f 2 (u 2 ) is continuous and decreasing on its domain, it has a unique inverse, f -1 2 , which is also decreasing:

where f -1 2 (0) = φ/μ and lim u2→+∞ f -1 2 (u 2 ) = ( φ -p 1 )/μ. We now restrict the domain of f -1

2 (1) < 0. Then, by the intermediate value theorem, there exists ū1 ∈ (1 -p 1 , 1) such that h(ū 1 ) = 0 and so f 1 (ū 1 ) = f -1 2 (ū 1 ). This steady state is unique as f 1 and f 2 are bijective on their domains. Thus, in the case ν 12 1, ν 1 1, ν 2 1, there exists a unique steady state, (ū 1 , ū2 ), with 1 -p 1 < ū1 < 1 and 1 μ ( φ -p 1 ) < ū2 < φ μ .

A c c e p t e d m a n u s c r i p t

Appendix D. Stability of the steady state for the deterministic model in the case ν 12 1, ν 1 1, ν 2 1

Here we prove that the steady state (ū 1 , ū2 ) is stable for all values of the parameters. Let

The Jacobian matrix is given by

where u 1 and u 2 take their steady state values, ū1 and ū2 , respectively. From Eq. (C-3) we have

and similarly from Eq. (C-4):

We consider Det(J) in the following three different cases:

Hence, the steady state (ū 1 , ū2 ) is locally asymptotically stable for all parameter values.

Appendix E. Validity of the normal approximation

In Section 6 we have approximated the stochastic bivariate continuous-time Markov process by a bivariate normal distribution with mean the stable steady state of the deterministic system. Here we develop estimates of the population sizes at which one must shift from a deterministic to a stochastic description. The bivariate continuous-time Markov process has state space S = {(n 1 , n 2 ) : n 1 , n 2 = 0, 1, 2, . . .} and the bivariate normal distribution has state space S = {(y 1 , y 2 ) : y 1 , y 2 ∈ R}. In order for the deterministic

A c c e p t e d m a n u s c r i p t

description to be a good approximation, the bivariate normal distribution should have sufficient probability density in the region R + × R + where S is included.

The details are as follows: the probability distribution function (pdf) for a bivariate normal random variable Y is given by

with Σ the covariance matrix, m the mean of the bivariate normal distribution and y an arbitrary point in

is chi-square distributed with 2 degrees of freedom and the pdf given in Eq. (E-1) is constant on ellipsoids in 2-dimensional Euclidean space. The centre of the ellipsoid is m, its shape and orientation are determined by Σ and the constant determines its size.

Our task is thus to find the ellipsoid in 2-dimensional Euclidean space that first touches one of the positive semi-axes. This is equivalent to finding points (y 1 , 0) and/or (0, y 2 ), with y 1 , y 2 ≥ 0 that minimise the Mahalanobis distance to the mean of the normal distribution [START_REF] Timm | Applied Multivariate Analysis[END_REF].

In what follows we restrict ourselves to the case considered in Section 6.2.1. In this case, ϕ 1 = ϕ 2 = ϕ, p 1 = p 2 = p and μ 1 = μ 2 = μ, the mean of the bivariate distribution is given by

with n = ϕ(2-p) 2μ , and the covariance matrix is given by

with σ 2 = 3p-4 4(p-1) n and ρ = p 3p-4 . The Mahalanobis distance from y to m is given by

(E-4)

Given the symmetry of the case under consideration, the first ellipsoid in 2-dimensional Euclidean space to touch one of the positive semi-axes, will in fact, touch both of the axis. This can be shown as follows. Consider y = (y 1 , 0) and compute d M (y, m):

this distance is minimised for y * 1 = (y * 1 , 0) = ((1 -ρ)n, 0) and

The same arguments can be carried out for y = (0, y 2 ) to show that the distance from the mean to the y 2 -axis is minimised for y * 2 = (0, y * 2 , 0) = (0, (1 -ρ)n) and

As the random variable Y-m) t Σ -1 (Y-m) is chi-square distributed with 2 degrees of freedom, we can conclude that

with α ∈ [0, 1], as the cumulative distribution function for a chi-square variable with 2 degrees of freedom is given by F (z, 2) = 1exp(-z/2) [START_REF] Timm | Applied Multivariate Analysis[END_REF]. The parameter α indicates the amount of probability density included in the ellipsoid. The particular value of α will be uniquely determined by the values of the parameters ϕ, μ and p. One can show that in fact

Note that in the limit p → 1 -the deterministic approximation is not valid as σ 2 diverges. We can see in Fig. 17 how the approximation performs depending on the specific values of ϕ, μ and p. For example, for ϕ = 1, μ = 1 and any value of p, the ellipsoid will never contain more than 40% of the total probability (right panel). On the other hand, for ϕ = 1000, μ = 1 and any value of p, the ellipsoid will contain almost 100% of the total probability (left panel, top curve), indicating that the deterministic approximation is very good. This analysis can be extended to the general case introduced in Section 2.