A mathematical perspective on CD4(+) T cell quorum-sensing.

Journal of Theoretical Biology 347 (2014) 160–175

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

A mathematical perspective on CD4 þ T cell quorum-sensing Joseph Reynolds a, Inês F. Amado b,c,d, Antonio A. Freitas b,c, Grant Lythe a, Carmen Molina-París a,n a

Department of Applied Mathematics, School of Mathematics, University of Leeds, Leeds LS2 9JT, UK Institut Pasteur, Départment d0 Immunologie, Unité de Biologie des Populations Lymphocytaires, Paris, France c CNRS, URA1961, Paris, France d GABBA, ICBAS, Universidade do Porto, Porto, Portugal b

H I G H L I G H T S

We model peripheral CD4 þ T cell population dynamics with an IL-2 quorum-sensing mechanism. We consider naive, IL-2 producing, IL-2 non-producing and regulatory CD4 þ T cells. We study peripheral CD4 þ T cell establishment in the absence of antigen, and its return to homeostasis after an immune challenge. In the presence of specific antigen, a larger number of regulatory T cells is required to maintain homeostasis. We study the behaviour leading to extinction of the IL-2 producing and regulatory T cell populations.

art ic l e i nf o

a b s t r a c t

Article history: Received 19 July 2013 Received in revised form 24 November 2013 Accepted 16 December 2013 Available online 31 December 2013

We analyse a mathematical model of the peripheral CD4 þ T cell population, based on a quorum-sensing mechanism, by which an optimum number of regulatory T cells can be established and maintained. We divide the population of a single T cell receptor specificity into four pools: naive, IL-2 producing, IL-2 non-producing, and regulatory CD4 þ T cells. Proliferation, death and differentiation of cells are introduced as transition probabilities of a stochastic Markov model, with the assumption that the amount of IL-2 available to CD4 þ T cells is proportional to the size of the population of IL-2 producing CD4 þ T cells. We explore the population dynamics both in the absence and in the presence of specific antigen. We study the establishment of the peripheral CD4 þ T cell pool from thymic output in the absence of antigen, and its return to homeostasis after an immune challenge, by steady state analysis of the deterministic approximation. The number of regulatory T cells at steady state is greater in the presence of antigen than in its absence. We also consider the stochastic dynamics of the model after an immune challenge, in particular the behaviour leading to ultimate extinction of the IL-2 producing and regulatory T cell populations. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Quorum-sensing IL-2 Regulatory T cells Stochastic Markov process Extinction

1. Introduction Regulatory T cells, characterised by the expression of the transcription factor FOXP3, make up 5–10% of the peripheral CD4 þ T cell pool in humans and mice (Sakaguchi et al., 1995; Sakaguchi, 2000; Seddon and Mason, 2000; Powrie and Maloy, 2003; Feuerer et al., 2010; Walker and Sansom, 2011; Josefowicz et al., 2012). They are able to regulate disease in adoptive transfer models, and their absence results in the development of autoimmune disorders caused by unconstrained proliferation of CD4 þ T cells (Sakaguchi, 2004). While a sufficient number of regulatory T cells is required to prevent the onset of autoimmune disease (Almeida et al., 2002; Fontenot et al., 2003), regulatory T cells have also been shown to suppress beneficial immune responses against tumours and viral infections (Antony et al., n

Corresponding author. E-mail address: [email protected] (C. Molina-París).

0022-5193/$ - see front matter & 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jtbi.2013.12.019

2005; Belkaid and Rouse, 2005). These experimental findings suggest that there exists an ideal size for the regulatory T cell pool, which allows for both the prevention of autoimmunity and the effective mounting of adaptive immune responses. Mathematical and computational models have been developed to propose potential cellular (or molecular) mechanisms that control the size of the regulatory T cell population. León et al. (2000, 2001, 2003), Carneiro et al. (2007), and Sepúlveda and Carneiro (2011) have developed the cross-regulation model, based on the interactions between populations of regulatory and effector CD4 þ T cells, and built on experimental studies that support the cell–cell contact hypothesis (Tang et al., 2005). That is, that interactions between regulatory and effector T cells occur at conjugation sites present on the surface of antigen presenting cells (APCs). The cross-regulation model predicts that the peripheral CD4 þ T cell repertoire is comprised of two distinct subsets (León et al., 2000, 2001, 2003; Carneiro et al., 2007; Sepúlveda and Carneiro, 2011): small clones of auto-reactive

J. Reynolds et al. / Journal of Theoretical Biology 347 (2014) 160–175

effector T cells, which are highly suppressed by their T cell receptor (TCR) related regulatory T cells, and effector T cell clones, which barely recognise self-antigen, but have the potential to strongly recognise non-self-antigen. Thus, an effective immune response can take place, without autoimmunity. Other approaches include Fouchet and Regoes, who analysed the interaction between antigen and six populations of APCs and regulatory T cells (Fouchet and Regoes, 2008). Kim et al. (2007) introduced a comprehensive model of immune responses that includes populations of CD4 þ T cells, CD8 þ T cells, APCs, virus infected cells and antigen. Feinerman et al. (2010) provided a model in which binding of IL-2 leads to the up-regulation of IL-2Rα (CD25) in both effector and regulatory T cells, favouring their subsequent ability to bind IL-2. Other mathematical models, based on the idea of bystander T cell activation (Burroughs et al., 2011; Kim et al., 2007), predict that an autoimmune state can arise following an immune response. According to the recent quorum-sensing (QS) hypothesis (Almeida et al., 2012), homeostasis of CD4 þ T cell numbers can be achieved by the ability of the population of regulatory T cells to sense the number of cytokine interleukin-2 (IL-2) producing CD4 þ T cells (Almeida et al., 2012). Quorum-sensing, thus, is a mechanism whereby an optimum ratio between CD4 þ effector and regulatory T cells can be established and maintained (Almeida et al., 2005, 2006a, 2012). The mechanism resembles that observed in species of bacteria, in which gene expression is regulated in response to the density of the bacterial population (Diggle et al., 2007; Miller and Bassler, 2001), but is mediated by another cell population, the regulatory T cell subset, and does not involve a change in gene expression. Furthermore, the CD4 þ T cell QS mechanism considered here also controls the dynamics of the regulatory T cell population itself. Here, we use a mathematical model of four interacting CD4 þ T cell subsets to explore the dynamics in the absence and in the presence of specific antigen. In the absence of antigen, we study the establishment of the peripheral CD4 þ T cell pool from thymic output, and its return to homeostasis after an immune challenge. In the presence of antigen, we make use of the QS mathematical model to understand the steady state of the IL-2 producing and regulatory T cell populations, and how to avoid autoimmunity. We divide the CD4 þ T cell pool into four subsets. The first, population 1, is the set of CD4 þ naive T cells, defined as not having yet encountered their specific antigen. The second subset, population 2, is the set of CD4 þ T cells which have encountered their specific antigen. These cells exhibit effector (or helper) function, such as the production of various cytokines to promote immune responses (Banchereau et al., 2012). In particular, we assume that effector cells produce the cytokine interleukin-2 (IL-2) (Wan and Flavell, 2009), and we refer to this population as the effector or the IL-2 producing population. The third subset, population 3, is the set of CD4 þ T cells, assumed to have encountered their specific (or cognate) antigen, but not to be producing IL-2. Lastly, population 4 is the set of CD4 þ regulatory T cells. These cells are assumed to suppress the activity of the effector T cell population (Almeida et al., 2006b), but to depend on the IL-2 produced by the effector T cells for their survival (Malek and Castro, 2010; Cheng et al., 2011) and proliferation. Proliferation, death and differentiation of cells within these populations will be introduced as transition probabilities of a stochastic Markov model. The quorum-sensing hypothesis will be implemented as follows:

the (per cell) death rate of regulatory T cells depends on the

number of IL-2 producing cells (Almeida et al., 2012), and in the absence of antigen, and due to the fact that regulatory T cells constitutively express the IL-2Rα chain, IL-2 induced proliferation only takes place in the regulatory T cell pool (Almeida et al., 2012).

161

The regulatory T cell suppression mechanism (Sakaguchi, 2004; von Boehmer, 2005) will be implemented mathematically as follows:

the rate at which an IL2-producing cell becomes a non-IL2 producing state is proportional to the number of regulatory T cells. The initial stages of an immune response are characterised by the presence of antigen that perturbs the homeostatic equilibrium of the CD4 þ T cell population, leads to proliferation of the antigen specific T cells and thus, to an increase in the number of IL-2 producing T cells. This, in turn, leads to regulatory T cell proliferation. Eventually, CD4 þ T cell homeostasis is re-established due to the suppression of the IL-2 producing CD4 þ T cell subset by the regulatory T cell population. We assume that the level of IL-2 is proportional to the number of IL-2 producing CD4 þ T cells. The structure of the paper is as follows. In Section 2, we introduce the stochastic quorum-sensing model and the transition probabilities that describe the dynamics of the four different CD4 þ T cell subsets. The deterministic approximation in the absence of specific antigen is discussed in Section 2.2. Two scenarios of immunological relevance are considered. In the first, we study how CD4 þ T cells establish themselves in the periphery (see Section 3.1). In the second, we study how the CD4 þ T cell population returns to its homeostatic equilibrium after antigenic challenge (see Section 3.2). In Section 4 we discuss, from a deterministic perspective, how the homeostatic equilibrium of the CD4 þ T cell population will change in the case of specific antigen presentation (see Section 4). In Section 5, we study the stochastic dynamics of the CD4 þ T cell population after antigenic challenge and its return to homeostasis. Given the potentially small number of specific CD4 þ T cells involved, we consider extinction of either regulatory or IL-2 producing T cells. We show that extinction takes place with probability one, and compute the temporal order in which extinction occurs. We also compute the probability of extinction and the expected time to extinction as a function of the parameter values. We finish the paper with a summary of our results and a discussion. Finally, we have included in the Appendix the mathematical details of the extinction problem for the stochastic model.

2. Stochastic model of CD4 þ T cell quorum-sensing The stochastic model will be defined in terms of the transition probabilities of a multi-variate Markov process. The state space of this Markov process, S, is the four-dimensional lattice N N N N, with N ¼ f0; 1; 2; …g. The transition probabilities describe the jumps the process may take on this lattice in a small time interval Δt. Note that, in this immunological context, jumps correspond to cell proliferation, death, differentiation or migration (thymic export). Let fXðtÞg be a four-dimensional Markov process, where each sub-chain fXi ðtÞg, i ¼ 1; 2; 3; 4, is a uni-variate Markov chain that describes one of the four subsets of the CD4 þ T cell population. Let the four-dimensional vector nðtÞ A S be a realisation of the multi-variate Markov process at time t, and assume that the process starts at N A S, that is, ProbfXð0Þ ¼ Ng ¼ 1:

ð1Þ

The state probabilities of the Markov process are denoted by pn ðtÞ and are defined as follows (Allen, 2003): pn ðtÞ ¼ ProbfXðtÞ ¼ n∣Xð0Þ ¼ Ng:

ð2Þ

The transition probabilities of the Markov process are derived from the assumed possible cellular fates (proliferation, death, differentiation or migration) of the four different CD4 þ T cell subsets (Almeida et al., 2012). We describe these cellular events

162


and their associated transition probabilities in the following section (see Fig. 1).

occurring in that subset is an increasing function of the subset size. The homeostatic proliferation rate is denoted by λiH . We encode these assumptions into the transition probabilities as follows:

2.1. Transition probabilities: immunological background

pfn:ni þ 1g;n ðΔtÞ ¼ λiH ni e ni =κi Δt þoðΔtÞ:

The transition probabilities are denoted by pm;n ðΔtÞ and defined as Allen (2003):

Specific antigen/IL-2 induced proliferation: We assume that all T cells, except naive T cells, may proliferate in response to the recognition of specific antigen mediated via TCR:peptide/MHC (T cell:APC) interactions. Furthermore, we assume that these T cell subsets may proliferate in response to other cytokine derived signals, such as IL-2. It has been shown that naive T cells proceed through a small number of divisions before exhibiting markers typical of effector T cells (Jelley-Gibbs et al., 2000). However, we simplify the model by assuming that antigenic stimulation of naive T cells causes them to differentiate into effector T cells, without proliferating first. CD4 þ effector T cells up-regulate the IL-2 receptor α chain (IL-2Rα or CD25) in response to their specific antigen (Feinerman et al., 2010). Therefore, it is reasonable to assume that effector cells become more sensitive to IL-2 following specific antigenic stimulation. We also assume that this is the case for CD4 þ IL-2 non-producing cells. For IL-2 producing and IL-2 non-producing cells, we assume that the IL-2 derived signals are only able to induce proliferation when specific antigen is being presented. Hence, the transition probabilities for IL-2 induced proliferation for these two subsets will depend on the function f(t) introduced above. Regulatory T cells up-regulate the α chain of the IL-2 receptor in larger numbers at the basal (pre-antigen stimulation) level and, unlike IL-2 producing cells and IL-2 nonproducing cells, are critically dependent on IL-2 for their survival and proliferation (Malek et al., 2002). This is the second hypothesis of the QS mechanism. For this reason, for the regulatory T cell subset: (i) we neglect TCR and/or IL-7 induced proliferation, and (ii) we assume that IL-2 induced proliferation does not depend on specific antigen levels encoded by f(t). The specific antigen/IL-2 induced proliferation rate is denoted by λiA for i ¼ 2; 3; 4. The transition probabilities for specific antigen/IL-2 induced proliferation are

pm;n ðΔtÞ ¼ ProbfXðt þ ΔtÞ ¼ m∣XðtÞ ¼ ng:

ð3Þ

We shall adopt the following notation for the transition probabilities in order to simplify the expressions. Let us introduce the state vector n ¼ ðn1 ; n2 ; n3 ; n4 Þ A S. By writing fn : ni 7 kg we refer to the new state vector in which the entry ni is replaced by ni 7k. For example, fn : n1 1g ¼ ðn1 1; n2 ; n3 ; n4 Þ. In a similar manner, if we refer to changes in two entries of n, we shall write, for example fn : n1 1; n2 þ 1g ¼ ðn1 1; n2 þ 1; n3 ; n4 Þ. Specific antigen presentation: We do not consider specific antigen a dynamical variable. However, we introduce a function, f(t), to represent presentation of specific antigen by APCs, where t denotes time. In the case of a successful immune response (which entirely clears antigen), f ðtÞ-0 þ as t- þ1. IL-2 level: We assume that the IL-2 level is proportional to the number of IL-2 producing cells, that is, n2. Thymic export: We assume that both naive T cells and regulatory T cells develop in the thymus and upon thymic maturation migrate to the periphery (Almeida et al., 2005). We shall assume that this occurs at a constant rate, νi (i¼ 1, 4). In a small time interval Δt, the transition probability encoding the influx of recent thymic emigrants into the peripheral CD4 þ T cell pool is pfn:ni þ 1g;n ðΔtÞ ¼ νi Δt þ oðΔtÞ;

i ¼ 1; 4:

ð4Þ

Death: In a small time interval Δt, it is assumed that for each CD4 þ T cell subset, a cell belonging to that subset may die. The probability of a natural death event occurring for all subsets, with the exception of the regulatory subset, is assumed to be proportional to the number of cells in that subset. The expression for the transition probability for such a death event is pfn:ni 1g;n ðΔtÞ ¼ μi ni Δt þoðΔtÞ;

i ¼ 1; 2; 3;

ð5Þ

with μi being the (per cell) death rate for i ¼ 1; 2; 3. The survival of regulatory T cells is assumed to depend on the available level of IL-2, and thus the number of IL-2 producing cells (Malek and Castro, 2010; Cheng et al., 2011). This is one of the hypotheses of the QS mechanism. The transition probability modelling the death of regulatory T cells is pfn:n4 1g;n ðΔtÞ ¼ μ4

κ4

κ 4 þ n2

n4 Δt þoðΔtÞ;

ð6Þ

where μ4 is the (per cell) death rate of regulatory T cells, and κ4 is the number of IL-2 producing cells at which the probability of death of a regulatory T cell is half the value of the probability when there are no IL-2 producing cells. Homeostatic proliferation: We assume that cells from the first three subsets may proliferate in response to homeostatic signals. A number of different signals constitute the full set of signals that a CD4 þ T cell requires under homeostatic conditions. The main signals regulating homeostasis of T cells are believed to be selfpeptide/MHC signals mediated via the T cell receptor (TCR), and the cytokine interleukin-7 (IL-7) (Seddon et al., 2003; Kimura et al., 2012). We assume that each population (naive, IL-2 producing, and IL-2 non-producing) has a different niche for homeostatic signals. Thus, we consider that these signals are limiting in the sense that as a population of T cells of subset i grows larger than a carrying capacity κi, i ¼ 1; 2; 3, the probability of a cell receiving a signal to proliferate decreases. For cell subsets smaller than κi, we assume that the probability of a proliferation event

pfn:ni þ 1g;n ðΔtÞ ¼ λiA f ðtÞ n2 ni Δt þ oðΔtÞ; pfn:n4 þ 1g;n ðΔtÞ ¼ λ4A n2 n4 Δt þoðΔtÞ:

ð7Þ

i ¼ 2; 3;

ð8Þ ð9Þ þ

Differentiation of naive T cells: We assume that naive CD4 T cells may become activated and start exhibiting effector (helper) functions in response to signals resulting from TCR interaction with their specific antigen (Jenkins et al., 2001). This transition probability will then depend on f(t), and is given by pfn:n1 1;n2 þ 1g;n ðΔtÞ ¼ α12 f ðtÞn1 Δt þ oðΔtÞ;

ð10Þ

where α12 is the naive T cell differentiation rate into IL-2 producing T cells. We additionally assume that naive T cells can differentiate directly into IL-2 non-producing cells through homeostatic-driven proliferation. This event is considered to be independent of TCR recognition of specific antigen (Ge et al., 2002). The transition probability for this differentiation event is pfn:n1 1;n3 þ 1g;n ðΔtÞ ¼ α13 n1 Δt þ oðΔtÞ;

ð11Þ

where α13 is the naive T cell differentiation rate into IL-2 nonproducing T cells. Differentiation of IL-2 producing T cells: We assume that IL-2 producing cells may stop producing IL-2 and become memory cells (Pepper and Jenkins, 2011). The transition probability is pfn:n2 1;n3 þ 1g;n ðΔtÞ ¼ α23 n2 Δt þ oðΔtÞ;

ð12Þ

where α23 is the IL-2 producing cell differentiation rate into IL-2 non-producing T cells. Activation of IL-2 non-producing T cells: We assume that IL-2 non-producing T cells can become activated and start producing


IL-2 in secondary responses to specific antigen (Berard and Tough, 2002). This transition probability will then depend on f(t), and is pfn:n2 þ 1;n3 1g;n ðΔtÞ ¼ α32 f ðtÞn3 Δt þ oðΔtÞ;

ð13Þ

where α32 is the IL-2 non-producing cell differentiation rate into IL-2 producing T cells. Suppression of IL-2 producing T cells: We assume that regulatory T cells suppress the ability of IL-2 producing T cells to secrete IL-2 (Sakaguchi, 2004; von Boehmer, 2005). This cellular process is introduced as a contact term: an IL-2 producing T cell will become an IL-2 non-producing T cell, if it receives a suppressive signal from a regulatory T cell. This transition probability is given by pfn:n2 1;n3 þ 1g;n ðΔtÞ ¼ β n2 n4 Δt þ oðΔtÞ;

ð14Þ

with β the suppression parameter. Remaining probabilities: In a small time interval Δt, it is possible that none of the above transitions may occur (Allen, 2003). The probability of no event occurring within a time interval Δt is " 3

pn;n ðΔtÞ ¼ 1 ν1 þ ν4 þ ∑ μi ni þ μ4 i¼1

3

3

i¼1

i¼2

κ4 n κ 4 þ n2 4

dm2 ¼ μ2 m2 þ λ2H m2 e m2 =κ 2 þ λ2A f ðtÞm22 þ α12 f ðtÞm1 dt α23 m2 þ α32 f ðtÞm3 β m2 m4 ;

ð20Þ

dm3 ¼ μ3 m3 þ λ3H m3 e m3 =κ 3 þ λ3A f ðtÞm2 m3 þ α13 m1 dt þ α23 m2 α32 f ðtÞm3 þ β m2 m4 ;

ð21Þ

dm4 ¼ ν4 μ4 m4 þ λ4A sIL 2 m4 ; dt

ð22Þ

where we have introduced the parameter sIL 2 , which represents an external source of IL-2 (which is, of course, independent of the population n2). In the absence of the QS hypothesis, one needs to assume that the death rate of regulatory T cells does not depend on the IL-2 producing cells. Inspection of (22) allows us to infer that the dynamics of the regulatory T cell population does not depend on any other population, in particular the population of IL-2 producing cells, n2. Thus, the steady state that characterises the population of regulatory T cells is given by mn4 ¼

þ ∑ λiH ni e ni =κ i þ ∑ λiA f ðtÞn2 ni þ λ4A n2 n4

163

ν4

ðμ4 λ4A sIL 2 Þ

;

ð23Þ

with the condition μ4 4 λ4A sIL 2 . Given the experimental observations that support the indexation of regulatory T cell numbers to the size of the IL-2 producing CD4 þ T cells (Almeida et al., 2005, 2006b), we conclude that a mathematical description of the CD4 þ T cell population requires a QS mechanism, as the one proposed earlier and analysed here (Almeida et al., 2012).

þ α12 f ðtÞn1 þ α13 n1 þ α23 n2 # þ α32 f ðtÞn3 þ βn2 n4 Δt þ oðΔtÞ: The probability of any other transition occurring within the proposed Markov process in a small time interval Δt is assumed to be oðΔtÞ (Allen, 2003).

3. Deterministic approximation in the absence of specific antigen

2.2. Deterministic approximation of the QS mathematical model Throughout this paper we assume that the CD4 þ T cell population being modelled is the subset of CD4 þ T cells specific to a given antigen. We write down a system of ordinary differential equations (ODEs) which approximate the average behaviour of the stochastic model. That is, we let mi(t) be an approximation to the first moments of the Markov chain fXi ðtÞg for i ¼ 1; 2; 3; 4. One may derive these equations from the stochastic model using the generating function technique and assuming that the kth central moment is zero for k Z2 (Allen, 2003). The system of ODEs is given by the following equations: dm1 ¼ ν1 μ1 m1 þ λ1H m1 e m1 =κ 1 α12 f ðtÞm1 α13 m1 ; dt

ð15Þ

ð16Þ

dm3 ¼ μ3 m3 þ λ3H m3 e m3 =κ3 þ λ3A f ðtÞm2 m3 þ α13 m1 dt þ α23 m2 α32 f ðtÞm3 þ βm2 m4 ;

ð17Þ

ð24Þ

dm2 ¼ μ2 m2 þ λ2H m2 e m2 =κ 2 α23 m2 βm2 m4 ; dt

ð25Þ

dm3 ¼ μ3 m3 þ λ3H m3 e m3 =κ 3 þ α13 m1 þ α23 m2 þ βm2 m4 ; dt

ð26Þ

dm4 κ4 ¼ ν4 μ4 m þλ m m ; dt κ 4 þ m2 4 4A 2 4

ð27Þ

3.1. CD4 þ T cell establishment in the periphery

ð18Þ

2.3. Deterministic approximation of the model without the QS hypothesis If we do not make use of the QS hypothesis, we need to modify the system of ODEs as follows: dm1 ¼ ν1 μ1 m1 þ λ1H m1 e m1 =κ 1 α12 f ðtÞm1 α13 m1 ; dt

dm1 ¼ ν1 μ1 m1 þ λ1H m1 e m1 =κ 1 α13 m1 ; dt

with initial conditions mð0Þ ¼ M, where mðtÞ ¼ ðm1 ðtÞ; m2 ðtÞ; m3 ðtÞ; m4 ðtÞÞ.

dm2 ¼ μ2 m2 þ λ2H m2 e m2 =κ2 þ λ2A f ðtÞm22 þ α12 f ðtÞm1 dt α23 m2 þ α32 f ðtÞm3 βm2 m4 ;

dm4 κ4 ¼ ν4 μ4 m þλ m m : dt κ 4 þ m2 4 4 A 2 4

We start the mathematical analysis of the deterministic model in the case when specific antigen is not being presented by APCs, by setting f ðtÞ ¼ 0. The system of ODEs simplifies to

ð19Þ

The first question that can be studied with the deterministic approximation of the QS model is the establishment in the periphery of the CD4 þ T cell pool from thymic output. That is, the homeostatic distribution of the subpopulations before any antigenic challenge has taken place. We will analyse the steady states of the model described by Eqs. (24)–(27). We begin with the observation that the differential equation for population m1 ðtÞ does not involve any other CD4 þ T cell subset and, thus, can be treated independently of the remaining equations. Let us introduce gðm1 Þ ¼ ν1 ðμ1 þ α13 Þm1 þ λ1H m1 e m1 =κ 1 :

ð28Þ n

Then the steady state solution of (24), m1 , is found by solving gðm1 Þ ¼ 0. Due to the presence of the term e m1 =κ1 , an analytical

164


solution for mn1 cannot be found. However, the existence and stability of such a steady state solution can be proved, as follows. Existence: In the limit m1 - þ 1, we have gðm1 Þ- 1. Furthermore, at m1 ¼ 0, we have gðm1 ¼ 0Þ ¼ ν1 4 0. Since g is a continuous function, by the intermediate value theorem, there exists at least one solution mn1 A ð0; 1Þ, such that gðmn1 Þ ¼ 0. That is, there exists at least one steady state solution for the population of naive T cells, such that mn1 4 0. Uniqueness: Assume that a solution mn1 4 0, such that gðmn1 Þ ¼ 0, exists. We have

λ1H e m1 =κ1 ¼ μ1 þ α13 n

ν1

mn1

:

ð29Þ

Let us introduce h1 ðmÞ ¼ λ1H e m=κ1 and h2 ðmÞ ¼ μ1 þ α13 ν1 =m for m 4 0. The derivatives of h1 and h2 with respect to m are given by 0 0 h1 ðmÞ ¼ ðλ1H =κ 1 Þe m=κ 1 o 0 and h2 ðmÞ ¼ ν1 =m2 4 0, respectively, for m 4 0. The derivatives of h1 and h2 are, respectively, monotonically increasing and decreasing, and continuous functions (for m 4 0). We conclude that (29) has at most one solution. Stability: The steady state mn1 4 0 is asymptotically stable if and only if dgðm1 Þ j n o0: dm1 m1

ð30Þ

Taking the first derivative of g we have dgðm1 Þ m1 m 1 = κ 1 ¼ ðμ1 þ α13 Þ þ λ1H 1 : e dm1 κ1

ð31Þ

We also know that by definition mn1 satisfies the following equation:

ν1 ðμ1 þ α13 Þmn1 þ λ1H mn1 e m1 =κ1 ¼ 0: n

ð32Þ

Rearranging Eq. (32) we have

ν1 þ mn1 ðλ1H e m1 =κ1 μ1 α13 Þ ¼ 0: n

ð33Þ

Since mn1 and ν1 4 0, the term inside the parenthesis must be negative and so the following inequality must hold: e m1 =κ 1 o n

μ1 þ α13 : λ1H

ð34Þ

We now evaluate Eq. (31) at mn1 : mn1 dgðm1 Þ n e m1 =κ 1 : mn1 ¼ μ1 þ α13 þ λ1H 1 dm1 κ1

ð35Þ

If mn1 4 κ 1 , then the RHS of (35) is negative, and mn1 is stable. Now suppose mn1 o κ 1 . Then, using the inequality in Eq. (34), we have mn1 dgðm1 Þ n e m1 =κ1 mn1 ¼ ðμ1 þ α13 Þ þ λ1H 1 dm1 κ1 mn μ1 þ α13 o ðμ1 þ α13 Þ þ λ1H 1 1 ¼

ðμ1 þ α13 Þmn1

κ1

κ1

λ1H

o 0:

Therefore, we conclude that mn1 40 is stable for all parameter values. IL-2 producing T cells: The value mn2 ¼ 0 is a steady state of (25). That is, as there are no source terms for the population of IL-2 producing cells, if m2 ðt ¼ 0Þ ¼ 0, then m2 ðtÞ ¼ 0 for all t Z 0. We do not consider any other steady states for the IL-2 producing T cell population, since we are interested in the peripheral CD4 þ T cell establishment from thymus output, before any T cell activation takes places. Thus, before specific antigen has been introduced, the number of IL-2 producing cells should be zero. IL-2 non-producing T cells: In a similar manner to the naive T cells, the exponential term in Eq. (26) implies that an analytical solution for mn3 cannot be found. However, setting the derivative

equal to zero in Eq. (26), and assuming that mn2 ¼ 0, gives an equation for mn3 , which is of exactly the same form as that for mn1 . The same arguments used to prove the existence and uniqueness of mn1 can be invoked to show that mn3 40 not only exists, but is unique. It is also stable, under the assumption that mn2 ¼ 0. Regulatory T cells: If m2 ðtÞ ¼ 0, then dm4 ¼ ν4 μ4 m4 ; dt

ð36Þ

and the steady state population is mn4 ¼

ν4 : μ4

ð37Þ

It is easy to show that mn4 is a stable solution in this case. Summary: The steady state analysed in this section corresponds to populations of cells that originate in the thymus and are established in the periphery. There are non-zero populations of naive, IL-2 non-producing and regulatory T cells. The steady state population of naive T cells is determined by a balance between cell death, proliferation, loss due to differentiation and an influx of new thymic emigrants. Similarly, the population of IL-2 non-producing T cells is determined by a balance between cell death, proliferation and differentiation of the naive T cell pool into IL-2 non-producing cells, which we have assumed to be due to homeostatic proliferation. The population of regulatory T cells is set by the balance between migration of new regulatory cells from the thymus to the periphery and cell death. If ν4 ¼ 0, that is, there is no further thymic output for the CD4 þ regulatory T cell population, the regulatory T cell pool will eventually become extinct. Since we have assumed that there is no specific antigen presentation, there can be no differentiation from the naive pool into the IL-2 producing one. If we start with no IL-2 producing cells, this population will always be equal to zero. Thus, the above steady state represents a scenario of T cell homeostasis prior to any specific antigen challenge. That is, it represents the steady state CD4 þ T cell population that, from thymic export, gets homeostatically established in the periphery. Naive T cells in the periphery require TCR signals from self-peptide/MHC complexes for their survival. These TCR signals can, on rare occasions, promote activation of naive T cells to exhibit effector functions. Additionally, IL-2 produced by effector cells of a given TCR specificity (or clonotype) may produce a bystander effect, supporting the regulatory T cell populations of other TCR specificities (Feinerman et al., 2010; Burroughs et al., 2011). 3.2. Homeostasis of CD4 þ T cells after a specific antigenic challenge The second question that can be studied with the deterministic approximation of the QS model is the re-establishment of peripheral homeostasis in the CD4 þ T cell pool after an immune challenge, and when antigen has been cleared. In this section, we explore a scenario in which a specific antigen challenge has already taken place. We assume that the specific antigen of the CD4 þ T cell population being modelled is no longer presented by APCs, and that there is no thymic output. That is, we assume that sufficient time has elapsed to allow all T cells of the specificity under consideration to have egressed from the thymus (Freitas and Rocha, 2000; Bains et al., 2009a, 2009b). This approximation can be justified as follows: estimates of thymic output for mice are 2.5 106 T cells per day (CD4 þ and CD8 þ ) (Scollay et al., 1980), and the αβ T cell receptor diversity for mice is of the order of 2 106 (Nikolich-Zugich et al., 2004). This means that, at most, one CD4 þ T cell of a given specificity is incorporated into the peripheral pool per day. The contribution from thymic output can then be neglected when compared with any of the proliferating terms in the post-challenge scenario considered in this section (homeostatic or specific antigen/IL-2 induced).


This post-challenge scenario can be described by initial conditions, such that all T cell subsets have non-zero cell numbers at time t ¼0. We set ν1 ¼ ν4 ¼ 0 and f ðtÞ ¼ 0. Under these assumptions, the deterministic model is dm1 ¼ μ1 m1 þ λ1H m1 e m1 =κ1 α13 m1 ; dt

ð38Þ

dm2 ¼ μ2 m2 þ λ2H m2 e m2 =κ2 α23 m2 β m2 m4 ; dt

ð39Þ

dm3 ¼ μ3 m3 þ λ3H m3 e m3 =κ3 þ α13 m1 þ α23 m2 þ βm2 m4 ; dt

ð40Þ

dm4 κ4 ¼ μ4 m þλ m m ; dt κ 4 þm2 4 4A 2 4

ð41Þ

κ λ ¼ ðμ2 þ α23 Þ exp 4 2κ 2 ⋆ 2H

165

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !! 4μ4 1þ 1 :

κ 4 λ4A

The third solution is sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! ! n λ2H e m2 =κ2 μ2 α23 κ4 4μ4 n n ðm2 ; m4 Þ ¼ 1 ; 1þ ; 2 κ 4 λ4A β

ð49Þ

ð50Þ

⋆

which exists and is stable provided λ2H o λ2H . The bifurcation diagram provided in Fig. 2 illustrates the steady state analysis for

with initial conditions mi ð0Þ a 0, i ¼ 1; 2; 3; 4. Naive T cells: The steady state solution of population m1 ðtÞ, denoted by mn1 , satisfies 0 ¼ μ1 mn1 þ λ1H mn1 e m1 =κ 1 α13 mn1 : n

ð42Þ n

There are two solutions to Eq. (42): the first is given by m1 ¼ 0 and the second by mn1 ¼ κ 1 log ðλ1H =ðμ1 þ α13 ÞÞ. The vanishing solution exists unconditionally and is stable provided λ1H o μ1 þ α13 . This stability condition can be shown as follows. The steady state mn1 ¼ 0 is asymptotically stable if and only if (Edelstein-Keshet, 2005) dgðm1 Þ j n o0; dm1 m1

ð43Þ

with gðm1 Þ ¼ μ1 m1 þ λ1H m1 e

m1 =κ 1

α13 m1 . We have

λ1 dgðm1 Þ ¼ μ1 þ λ1H e m1 =κ 1 H m1 e m1 =κ 1 α13 ; dm1 κ1

ð44Þ

which evaluated at the steady state vanishing solution, mn1 ¼ 0, yields dgðm1 Þ ¼ μ1 þ λ1H α13 : ð45Þ n dm1 m1 ¼ 0 Thus, in order for the vanishing solution to be stable, the condition λ1H o μ1 þ α13 needs to be imposed. On the other hand, the non-zero solution exists and is stable provided λ1H 4 μ1 þ α13 . This condition is required to ensure that the argument of the logarithmic term is positive. For the analysis of the remaining CD4 þ T cell populations, we shall assume λ1H 4 μ1 , but λ1H o μ1 þ α13 , so that mn1 ¼ 0 is a stable steady state. This can be justified as follows: during an immune challenge, most, if not all, naive T cells specific to the antigen presented will differentiate to become IL-2 producing T cells. Under the assumptions of no thymic output and antigen clearance, all naive T cells have differentiated and thus, the naive T cell population vanishes, so that this CD4 þ population subset is unable to reconstitute itself after the challenge. IL-2 producing T cells and regulatory T cells: The steady state solutions for populations m2 and m4 are found by solving the following pair of coupled equations: 0 ¼ μ2 mn2 þ λ2H mn2 e m2 =κ 2 α23 mn2 βmn2 mn4 ; n

0 ¼ μ4

Fig. 1. Representation of the cellular birth events that can increase the number of T cells in any of the four CD4 þ pools. Dashed lines indicate the events that depend on recognition of specific antigen. Not shown are cell death processes that occur in each T cell pool. The death rate of regulatory T cells, n4, depends on the number of IL-2 producing cells, n2 (QS mechanism). Thymocytes are exported to the periphery into naive, n1, and regulatory, n4, T cells (teal arrows). Homeostatic proliferation takes place for naive, n1, IL-2 producing, n2, and IL-2 non-producing, n3, T cells (blue arrows). IL-2 induced proliferation takes place for IL-2 producing, n2, IL-2 nonproducing, n3, and regulatory, n4, T cells (red arrows). The solid red arrow of IL-2 induced proliferation for n4 is a QS mechanism. Differentiation and activation processes are shown in black. Suppression of effector by regulatory T cells is shown in purple. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.)

ð46Þ

κ4 mn þ λ mn mn : κ 4 þ mn2 4 4A 2 4

ð47Þ n

n

n

n

There exist three solutions ðm2 ; m4 Þ with m2 ; m4 Z 0. The first solution is ðmn2 ; mn4 Þ ¼ ð0; 0Þ, which always exists and is stable provided λ2H o μ2 þ α23 . The second solution is λ2H ;0 ; ð48Þ ðmn2 ; mn4 Þ ¼ κ 2 log μ2 þ α23 which exists provided λ2H 4 μ2 þ α23 , and is stable provided μ2 þ α23 o λ2H o λ⋆2H , where

Fig. 2. Bifurcation diagram for the IL-2 producing and regulatory T cell subsets under the assumptions ν1 ¼ ν4 ¼ 0 and f ðtÞ ¼ 0. For proliferation rates of the IL-2 producing T cells ðλ2H Þ greater than λ⋆ 2H , the steady state number of IL-2 producing cells is held constant by a non-zero population of regulatory T cells. Note that in the ⋆ region λ2H 4 λ2H , increasing the proliferation rate of IL-2 producing cells leads to a linear increase in the number of regulatory T cells at steady state. The vertical black dashed lines at μ2 þ α23 and λ⋆ 2H have been drawn to separate the three different stability regions. The red (regulatory) and green (IL-2 producing) dashed lines in the third region ðλ2H 4λ⋆ 2H Þ indicate the unstable steady state solution. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.)

166


Fig. 3. Deterministic trajectories in the absence of thymic output and specific antigen, computed from various initial conditions, and within the stochastic finite state space. Parameter values are described in Table 1. Blue dots represent the steady state solution that corresponds to extinction of both cell types, IL-2 producing and regulatory T cells. (A) The green dot 1 2 represents the stable steady state solution that exists when μ2 þ α23 o λ2H o λ⋆ h . (B) The red dot represents the stable steady state solution that 2H . In this case λ2H ¼ 2:2 10 1 2 exists when λ2H 4 λ⋆ . In this case λ ¼ 8 10 h . (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.) 2 2H H

the populations of IL-2 producing and regulatory T cells. We note that the existence and stability of the three steady states depend on the parameter λ2H , which is the homeostatic proliferation rate of the effector T cell population. This parameter encodes the level of auto-reactivity of the CD4 þ T cell clonotype under considera⋆ tion. As λ2H increases above λ2H , the number of regulatory T cells required to suppress the effector population also increases. IL-2 non-producing T cells: The steady state solution for the population of IL-2 non-producing T cells is found by solving the following equation: 0 ¼ μ3 mn3 þ λ3H mn3 e m3 =κ 3 þ α23 mn2 þ βmn2 mn4 : n

ð51Þ

An analytical solution for mn3 cannot be found. However, if we assume mn2 and mn4 4 0, a solution can be shown to exist, be unique and be unconditionally stable. This is shown in an identical manner to the proof of existence, uniqueness and stability of mn1 provided in Section 3.1. Summary: The steady state, under the assumptions of no specific thymic output and antigen clearance, has no naive T cells. A steady state solution exists with non-zero population sizes of IL-2 producing and regulatory T cells. This describes a tolerant state of the CD4 þ population; by tolerant, we mean a state in which the steady state number of effector (IL-2 producing) T cells is less than the “natural” carrying capacity of this population in the absence of regulatory T cells. Given that effector T cells may promote the activity of some components of the innate immune system (for example, by secretion of interferon gamma), limiting the effector population number is presumably desirable (Banchereau et al., 2012). When the proliferation rate of IL-2 producing T cells is low, μ2 þ α23 o λ2H o λ⋆2H , the steady state number of IL-2 producing T cells is always smaller than the corresponding steady state number ⋆ of IL-2 producing cells in the parameter region λ2H 4 λ2H (see Fig. 2). This suggests that a population of IL-2 producing cells with poor proliferative capacity does not need regulatory T cell suppression to remain tolerant, as defined above. Another interesting feature of our analysis is the independence, both in the steady states and their stability, of the effector T cell subset on the parameter β. Varying β only affects the steady state number of regulatory T cells. This suggests that at equilibrium, regulatory T cells compensate for either a diminished or an enhanced ability to suppress IL-2 producing T cells by increasing or decreasing their population size.

Table 1 Parameter values used to generate numerical results. The two values given for λ2H are chosen to ensure stability of the types of steady solution shown in the middle and right regions of the bifurcation diagram in Fig. 2, respectively. For the middle region the stability interval for λ2H is [μ2 þ α23 ¼ 0:02; λ⋆ 2H 0:023], hence the specific choice λ2H ¼ 0:022 h

1

.

Parameter

Value

μ2

1 10 2

h1

μ3

1 10

2

h1

1 10

2

μ4 λ2H λ3H λ2A

Units

2:2 10

h1

2

=8 10

2

h1

5 10

2

h1 cells 1 h 1

1 10

5

λ3A

1 10

5

cells 1 h 1

λ4A

1 10 5

cells 1 h 1

κ2

6 10

2

cells

κ3

2 103

cells

κ4

1 101

cells

α23

1 10 2

h1

α32

1 10

3

h1

β

1 10 4

cells 1 h 1

⋆

Finally, for values of λ2H 4 λ2H , the third steady state (mn2 4 0 and mn4 4 0) is globally stable. From a deterministic perspective, and unless the initial condition for the population of regulatory T cells is given by m4 ðt ¼ 0Þ ¼ 0, we shall always observe a regulated number of IL-2 producing cells; that is, mn2 is smaller than the carrying capacity of IL-2 producing cells. However, from a stochastic perspective, there exists a non-zero probability of a fluctuation causing the regulatory T cell population to become extinct. In Section 5, we shall explore both the probability of such an event occurring and the average time to extinction. Fig. 3 shows deterministic trajectories for a range of initial conditions and for parameter values as described in Table 1, in the absence of thymic output and specific antigen.

4. Deterministic approximation in the presence of specific antigen A third question that can be studied with the deterministic approximation of the QS model is the dynamics of an immune


response. It is interesting to see that the QS hypothesis is not only a mechanism of CD4 þ T cell homeostasis, but a mechanism that allows immune responses to take place, yet avoids autoimmunity (Almeida et al., 2012). In the previous sections, we have restricted ourselves to the case f(t)¼ 0, that is, to cases where antigen specific for the CD4 þ T cells is not presented. In this section, however, we consider the behaviour of a clonotype of CD4 þ T cells when its specific antigen is presented at a constant level. We set f ðtÞ ¼ c, where c 4 0 is constant. Antigen is presented to T cells by professional APCs, such as dendritic cells or macrophages (Janeway Jr and Travers, 2001). APCs collect antigen in the lymph and tissue via phagocytosis. Inside the APC, antigen is broken down further into peptides, which are subsequently displayed on the surface of the cell as peptide/MHC complexes (MHC class II complexes in the case of CD4 þ T cells). A naive CD4 þ T cell, upon successful binding with a dendritic cell presenting antigen of the correct specificity, will become activated. After activation, the naive T cell will undergo a phenotypic change and elicit effector functions. For CD4 þ T cells, these effector functions serve to promote the functional ability of CD8 þ effector T cells and boost the innate components of the immune system by, for example, releasing pro-inflammatory cytokines. During a successful immune response, the target of T cells (infected cells) is destroyed and cleared from the body. In a typical response, the elimination of infected cells leads to a decrease in the amount of antigen available for uptake and presentation by APCs. We now consider the steady state behaviour of a clonotype of CD4 þ T cells when its specific antigen is presented at a constant level. This implies f ðtÞ ¼ c, where c 4 0 is constant. We further assume that thymic output is zero, and that consequently there are no naive CD4 þ T cells. That is, all naive T cells have been activated due to the presentation of specific antigen. The deterministic

167

equations describing the model, where c has been absorbed into the relevant parameters, are dm2 ¼ μ2 m2 þ λ2H m2 e m2 =κ2 þ λ2A m2 m2 þ α32 m3 α23 m2 βm2 m4 ; dt

ð52Þ dm3 ¼ μ3 m3 þ λ3H m3 e m3 =κ3 þ λ3A m2 m3 þ α23 m2 α32 m3 þ βm2 m4 ; dt

ð53Þ dm4 κ4 ¼ μ4 m þλ m m : dt κ 4 þ m2 4 4 A 2 4

ð54Þ

The non-linear system of ODEs (52)–(54) admits multiple steady state solutions. We focus on non-vanishing steady state solutions for all three populations. Other steady state solutions, which may be stable, can be found (numerically) in the following combinations: mn2 ¼ mn3 ¼ mn4 ¼ 0, and mn2 ; mn3 4 0; mn4 ¼ 0. Although we did not observe any other vanishing steady states during numerical exploration, this may not be an exhaustive list. Let us look for steady state solutions such that mn2 ; mn3 ; mn4 4 0. The steady state mn2 is unique and found directly from Eq. (54): sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! κ4 4μ4 n 1 : ð55Þ 1þ m2 ¼ 2 κ 4 λ4A We can make use of (52), and write down an expression for mn4 , which depends on mn2 and mn3 . This expression is given by

λ2H e m2 =κ2 μ2 α23 λ2A ðmn2 Þ2 þ α32 mn3 þ : β βmn2 n

mn4 ¼

ð56Þ

The above equation provides a unique value of mn4 for a given (but fiducial) steady state value of mn3 . The problem of determining the number of positive steady state solutions can then be reduced to determining how many stable steady state solutions can be found

Fig. 4. The red line gives an example plot of f(m). Black lines give various examples of g(m) for different choices of c2 and c3. If c2 4 0, then we may either have no ðc3 4 0Þ or one ðc3 o 0Þ fixed point solution. If c2 o 0, we may either have no solution or up to three solutions. For any choice of c2 and c3, there is at most one solution which satisfies the necessary condition for asymptotic stability. We conclude that we have at most one stable steady state solution. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.)

168


for m3. Eq. (52) can be written as c2 α23 mn2 þ α32 mn3 ¼ βmn2 mn4 ;

ð57Þ

where c2 ðμ2 ; λ2H ; λ2A ; κ 2 ; μ4 ; λ4A ; κ 4 Þ ¼ μ2 mn2 þ λ2H mn2 e m2 =κ 2 þ λ2A mn2 mn2 : n

ð58Þ We substitute this expression into (53) to find c2 þ c3 mn3 þ λ3H mn3 e m3 =κ 3 ¼ 0 n

ð59Þ

where c3 ðμ3 ; λ3A ; μ4 ; λ4A ; κ 4 Þ ¼ λ3A mn2 μ3 :

ð60Þ

The constants c2 and c3 may be positive or negative real numbers. Let us define f ðmÞ ¼ λ3H me m=κ3 ;

gðmÞ ¼ c2 þ c3 m:

ð61Þ

We know that dm3 ¼ gðm3 Þ f ðm3 Þ: dt

ð62Þ

Stability of mn3 4 0 requires that dgðm3 Þ df ðm3 Þ n o 0: dm3 dm3 m3

Fig. 5. Deterministic trajectories when f ðtÞ ¼ 1 are shown by solid lines. Solutions for f ðtÞ ¼ 0 are shown by dashed lines. Initial conditions are 3 10 IL-2 producing T cells, 3 103 IL-2 non-producing T cells, and 3 102 regulatory T cells. Parameters used are given in Table 1. For f ðtÞ ¼ 1 and at equilibrium, we have a greater number of regulatory T cells and non-IL-2 producing T cells.

ð63Þ

The definition of derivative and the fact that f ðmn3 Þ ¼ gðmn3 Þ allows us to conclude that mn3 4 0 is a stable steady state if gðmn3 þ εÞ of ðmn3 þ εÞ for ε-0 þ . In Fig. 4 we provide an example plot of f(m) (red curve), which for positive λ3H and κ3 has a qualitatively fixed shape. Black lines give various examples of g(m) for different choices of c2 and c3. If c2 4 0, then we may either have no (c3 40) or one (c3 o0) fixed point solution. If c2 o0, we may either have no solution or up to three solutions. For any choice of c2 and c3, there is at most one solution which satisfies the necessary condition for asymptotic stability. We conclude that we have at most one stable steady state solution in the case of constant specific antigen. We conclude that the model, under the assumptions listed at the start of this section (no thymic output, no naive T cells, and constant presentation of antigen), exhibits one possible stable steady state solution. The solution for mn2 is the same as the solution for mn2 in (50), found under the assumption of no specific antigen presentation. We further note that the solution for mn4 is greater than the solution for mn4 when f ðtÞ ¼ 0 [given by (50)]. The QS model predicts that, when antigen specific to the CD4 þ T cell population is presented at a constant level, the number of regulatory T cells at steady state is greater than the number required to maintain homeostasis is the absence of antigen. The quorum sensing hypothesis, as implemented with this model, predicts that the number of IL-2 producing T cells found at steady state (assuming parameters are such that all populations are nonzero) is the same, whether the IL-2 producing population sees specific antigen or not. In the latter case, when the IL-2 producing population is exposed to specific antigen, a greater number of regulatory T cells is required to maintain the steady state. Of further note is that the steady state solution for the IL-2 producing population does not depend on the parameter describing the rate of suppression, β. Rather, varying β has the effect of controlling the number of regulatory T cells at steady state. In Fig. 5 we present numerical solutions of Eqs. (52)–(54) in the case that f ðtÞ ¼ 0 (dashed lines), and in the case that f ðtÞ ¼ 1 (solid lines). Parameters used are given in Table 1. If antigen is presented at a constant level, a greater number of regulatory T cells is observed, whereas the steady state number of IL-2 producing T cells remains the same. Overall, there is a modest increase in the

Fig. 6. Deterministic solutions for parameters given in Table 1, in which λ4A has been reduced by a factor of 10 ðλ4A ¼ 10 6 Þ. The number of IL-2 producing T cells at equilibrium has increased by approximately 3-fold, whereas the number of regulatory T cells has decreased by approximately 2-fold. Decreasing λ4A further (or increasing μ4) may result in the loss of stability with trajectories tending towards infinity. This can be rectified by replacing each parameter λiA ði ¼ 2; 3; 4Þ by a density dependent term in the same manner as done for parameters of the form λjH ðj ¼ 1; 2; 3Þ.

total number of CD4 þ T cells in the case of constant antigen presentation. Autoimmunity: Let us assume now that the CD4 þ T cell clonotype being modelled is specific to self-antigen. Without loss of generality, let us also assume that self-antigen is presented by APCs at a constant rate, so that the dynamics of the clonotype is governed by Eqs. (52)–(54). We discuss the requirements for the QS model to be in a regime in which autoimmunity can occur. Let us use the trajectories in Fig. 5 as an example of a self-reactive clonotype, in which autoimmunity is controlled due to the relatively small number of IL-2 producing T cells. That is, we assume tolerance can be maintained with a non-zero population of IL-2 producing T cells, provided the size of this population is relatively small (compared to the size of the regulatory T cell population). We suppose that a


typical feature of tolerance is an abundance of regulatory T cells relative to the number of IL-2 producing cells ðmn4 4 mn2 Þ. Given this example of a tolerant state, we define an autoimmune state as one in which the number of IL-2 producing T cells is in excess relative to the number of regulatory T cells. As discussed earlier, the steady state number of IL-2 producing T cells is a function of the parameters λ4A , μ4 and κ4. The question of how to increase the steady state number of IL-2 producing T cells is clear: increase κ4, μ4, decrease λ4A , or a combination of these. In Fig. 6, we decrease the value of λ4A by a factor of 10 (set λ4A ¼ 10 7 cells 1 h 1). The number of IL-2 producing T cells increases by approximately 3-fold, whereas the number of regulatory T cells decreases by approximately 2-fold. Varying the parameters in this manner allows the ratio between mn2 and mn4 to approach the autoimmune scenario described above. The findings from this model suggest that in the event that a clonotype is self-reactive, autoimmunity may arise as a consequence of poor proliferation and/or high death rate of regulatory T cells. Of interest in this model is the fact that the rate of suppression is irrelevant for controlling the number of IL-2 producing T cells at steady state. Indeed, if suppression is reduced, the regulatory population expands to compensate.

5. Extinction: a stochastic analysis In this section, we consider the stochastic model under the assumptions that specific antigen is not presented, f ðtÞ ¼ 0, and thymic output is zero, ν1 ¼ ν4 ¼ 0. This corresponds to the scenario analysed in Section 3.1 from a deterministic perspective. We evaluate the probability of regulatory T cell extinction conditioned on the non-extinction of the IL-2 producing T cell population, and the average time one must wait before regulatory T cell extinction takes place (Allen, 2003). We also compute the probability of IL-2 producing T cell extinction conditioned on the non-extinction of the regulatory T cell population, and the average time to IL-2 producing T cell extinction. In Appendix A, we show that the ultimate fate of the CD4 þ T cell subsets is extinction, and present some results for general birth and death processes (Allen, 2003). It is of interest to record which population is extinguished first. We shall study whether the process visits (before extinction takes place) the region in which the number of IL-2 producing T cells is greater than their characteristic size when regulated. We consider the problem of regulatory T cell extinction conditioned on the non-extinction of the IL-2 producing T cell population, and then study the problem of IL-2 producing T cell extinction conditioned on the nonextinction of the regulatory T cell population. 5.1. Numerical results: extinction of regulatory T cells Fig. 7(A) and (B) presents numerical (Gillespie) simulations of the stochastic process for two different sets of parameters, differing in the value of λ2H . Fig. 7(C) and (D) makes use of the method described in Section A.2.1 to compute probabilities of extinction. Finally, Fig. 7(E) and (F) makes use of the matrix method described in Section A.2.2 to compute the expected time to extinction for both sets of parameter values. The regions R and R, introduced in Appendix A, are given by R ¼ fðn2 ; n4 Þ : n2 a 0; n4 ¼ 0g

and

R ¼ fðn2 ; n4 Þ : n2 ¼ 0g:

ð64Þ

Probabilities of regulatory T cell extinction conditioned on the non-extinction of the IL-2 producing T cell population, and the expected time to extinction, are given in Table 2. We note that the value of λ2H ¼ 0:022 corresponds to a parameter set that yields a deterministic stable steady state with no regulatory T cells, and

169

λ2H ¼ 0:080 to a parameter set that yields a deterministic stable steady state with both effector and regulatory T cells present. 5.2. Numerical results: extinction of IL-2 producing T cells Fig. 8(A) and (B) makes use of the method described in Section A.2.1 to compute probabilities of extinction of IL-2 producing T cells, conditioned on the non-extinction of regulatory T cells. We have generated numerical realisations of the stochastic process for two different sets of parameters, as described in Table 1, differing only in the value of λ2H . In Fig. 8(C) and (D), we make use of the matrix method described in Section A.2.2 to compute the expected time to extinction for both sets of parameter values (Table 3). Summary: We consider two sets of parameter values. In the first, there are no regulatory cells in the steady state of the deterministic model of Section 3.2. However, we may choose an initial condition with numbers of regulatory T cells sufficiently large to drive the IL-2 producing population to extinction. In the second case, the steady state of the deterministic model corresponds to a situation where a non-zero population of regulatory cells moderates the size of the IL-2 producing population. Here, by choosing an initial condition with large numbers of IL-2 producing cells, we are able to explore the possibility that this population can avoid regulation.

6. Conclusions The quorum-sensing mechanism allows an IL-2 producing population to be maintained at a regulated size that does not depend on the absence or presence of specific antigen. In the latter case, a larger number of regulatory T cells is required to maintain the homeostatic distribution of the CD4 þ T cell population. Such scenarios include the late part of an immune response, once it has peaked (Almeida et al., 2012). The model presented here does not include competition for IL-2 at the receptor/molecular level. However, implicit in choosing the proliferation of regulatory T cells to be antigen independent is the assumption that regulatory T cells are more sensitive to IL-2 than IL-2 producing cells (Feinerman et al., 2010). We do not assume that this limits the proliferative capacity of IL-2 producing cells, as in our model the limiting factor of proliferation for IL-2 producing cells is the number of cells themselves; regulatory cells play a role through the suppressive term β (Wing et al., 2008; Qureshi et al., 2011). A possible extension of the model is to include a carrying capacity, for the IL-2 induced proliferation terms, which encodes some concept of space or availability of other trophic factors; availability would become limiting for large enough cell numbers. Another similar approach is that of Yates et al. (Callard et al., 2003), where T cells are assumed to induce fratricidal apoptosis due to Fas–FasL interactions. The simplest method of including such a mechanism in the model would be to add non-linear death terms in the ODEs for the populations of antigen-experienced cells (IL-2 non-producing and IL-2 producing T cells). In vivo, and due to the cross-reactivity of T cell clonotypes, IL-2 producing T cells might be suppressed by regulatory T cells of different TCR specificities. Interclonal suppression would limit the degree of proliferation the effector population could undergo. Additionally, regulatory T cells may make use of IL-2 produced by effector T cells of different TCR specificities, a hypothesis not considered in our model, but that should be included in future mathematical models of CD4 þ T cell populations (Müller et al., 2012). A regulatory T cell phenotype can be induced in non-regulatory cells (Chen et al., 2003; Kretschmer et al., 2005; Bettelli et al., 2006) via the transforming growth factor TGF-β, perhaps providing

170


1

Fig. 7. (A) Gillespie realisation of populations 2 and 4 for the parameter set given in Table 1, taking λ2H ¼ 0:022 h . Using this value ensures the deterministic steady state of the form ðnn2 4 0; nn4 ¼ 0Þ is stable. Initial conditions were chosen such that the stochastic realisation begins with 30 IL-2 producing T cells and 300 regulatory T cells and is run 1 for 30 days. (B) Gillespie realisation for the same initial conditions as in (A), with a time interval of two years, and the same parameter set except λ2H ¼ 8 10 2 h . This proliferation rate guarantees the stability of a deterministic steady state of the form ðnn2 4 0; nn4 4 0Þ. As the populations do not exceed 250 effector T cells or 550 regulatory T cells, we impose these values for N2 and N4, respectively. (C) Probability of reaching the region R from a general state ðn2 ; n4 Þ, with a maximum state space given by N 2 ¼ 250 and N 4 ¼ 550. The parameter set is the one used for (A). (D) Probability of reaching the region R from a general state ðn2 ; n4 Þ with a maximum state space N 2 ¼ 250 effector cells and N 4 ¼ 550 regulatory cells with parameters as chosen in (B). For this parameter set the probability of population 4 going extinct before population 2 is 10 12 for most initial conditions. (E) Expected time to absorption in the region defined by R [ ð0; 0Þ with a maximum state space of N 2 ¼ 250 effector T cells and N 4 ¼ 550 regulatory T cells. Parameters have been chosen as in (A) and (C). (F) Expected time to reach R [ ð0; 0Þ with same parameters as chosen in (B) and (D).


171

Table 2 Probabilities and expected times to extinction that correspond to Fig. 7. The table provides the conditional probability of extinction of the regulatory T cell population and the expected time to this extinction event given initial conditions for λ2H ¼ 0:022 h λ2H ¼ 0:022 h

1

1

and λ2H ¼ 0:080 h

1

.

λ2H ¼ 0:080 h

1

IC ðn2 ; n4 Þ

Prob n4 ext

Expected time (days)

IC ðn2 ; n4 Þ

Prob n4 ext

Expected time (years)

(40,100) (80,200) (120,300) (160,400) (200,500)

0.38732 0.12850 0.039915 0.012389 0.0038638

41.794 34.339 31.820 31.243 31.314

(40,100) (80,200) (120,300) (160,400) (200,500)

9.0657e 12 9.0657e 12 9.0657e 12 9.0657e 12 9.0656e 12

1552.2 1552.2 1552.2 1552.2 1551.9

Fig. 8. (Upper plots) Probability of extinction of the IL-2 producing population conditioned on survival of the regulatory population. The parameter set used for the left (right) plot corresponds to the same parameter set used for the left (right) plot in Fig. 7. (Lower plots) Expected time to absorption in the region R [ ð0; 0Þ with parameter sets corresponding to the above plots. Plot (D) looks qualitatively similar to plot (F) in Fig. 7. This is due to the fact that the most probable path of extinction in Fig. 7(F) is by absorption in the region R. Furthermore, the amount of time spent in region R is short relative to the timescales shown in (D).

additional protection against losses in the regulatory T cell population. It is, however, unclear if regulatory T cells can be induced from all non-regulatory phenotypes or just naive phenotypes. These considerations are out of the scope of this paper. Recent experimental studies have provided a molecular basis for the cell-extrinsic function of the co-receptor CTLA-4 (Qureshi et al., 2011). As regulatory T cells constitutively express the co-receptor CTLA-4 on their surface, this mechanism can allow them to capture, via CTLA-4 trans-endocytosis of CD80 and CD86, these ligands from APCs. In this way, regulatory T cells can deplete CD80 and CD86 molecules from the surface of APCs, and in turn suppress the activation of CD4 þ T cells. This mechanism will be explored in future work.

If a tolerant steady state exists and is stable, the model predicts that the equilibrium number of IL-2 producing T cells is the same regardless of whether these cells see specific antigen or not. In the former case, we observe a larger number of regulatory T cells to compensate for the increased turnover of effector cells. The model further predicts that an imbalance in the number of regulatory and effector T cells can lead to large populations of effector cells. However, this is only in the event these cells see specific antigen. Such an imbalance could result from a stochastic fluctuation, or from the immune activity of other T cell clonotypes during an immune response. Exploring potential mechanisms driving this imbalance would require a more comprehensive future model to be developed.

172


Table 3 Conditional probability of IL-2 producing T cell extinction conditioned on the nonextinction of the regulatory T cell population, and the expected time to this extinction. λ2H ¼ 0:022 h

1

λ2H ¼ 0:080 h

IC ðn2 ; n4 Þ

Prob n2 ext

(40,100) (80,200) (120,300) (160,400) (200,500)

0.61268 34.561 0.87150 19.493 0.96008 11.774 0.98761 7.9934 0.99613 5.9972

Expected time (days)

levels are proportional to the size of the IL-2 producing subset. Such an assumption does not take into account any notion of space or competition for IL-2, which may be the case in vivo.

1

Appendix A. Extinction: a stochastic analysis

IC ðn2 ; n4 Þ

Prob n2 ext

Expected time (years)

(40,100) (80,200) (120,300) (160,400) (200,500)

1.0000 1.0000 1.0000 1.0000 1.0000

1552.2 1552.1 1552.1 1552.1 1551.9

A.1. Extinction for a regular birth and death process For a regular uni-variate birth and death process fXðtÞg we have (Allen, 2003) pji ðΔtÞ ¼ ProbfXðt þ ΔtÞ ¼ j∣XðtÞ ¼ ig

ð65Þ

8 λi Δt þ oðΔtÞ > > > > < μ Δt þ oðΔtÞ i ¼ > 1 ðλi þ μi ÞΔt þoðΔtÞ > > > : oðΔtÞ

ð66Þ

if j ¼ i þ 1; if j ¼ j 1; if j ¼ i; otherwise;

where λi 4 0 for i Z1, μi 4 0 for iZ 1, and λ0 ¼ μ0 ¼ 0. For such a process, the state defined by i ¼0 is an absorbing state (Allen, 2003). Karlin and McGregor (1957) provide conditions for guaranteed absorption (or extinction) of a birth and death process of the form above, as well as giving conditions for absorption in a finite time. Let us introduce rn as follows: rn ¼

λ 1 λ 2 ⋯λ n 1 μ2 μ3 ⋯μn

for n Z 1:

ð67Þ

A sufficient condition for guaranteed absorption of the process at state n¼ 0 is that the series 1

1

n¼1

λn r n

∑

Fig. 9. For the first seven days we have set f ðtÞ ¼ 0. Between days 7 and 14, we set f ðtÞ ¼ 1. After day 14 we set f ðtÞ ¼ e 0:027ðt 1424Þ (an exponentially decreasing function which reduces from 1 to 1/100 in approximately seven days). The para1 meter set used for these trajectories is given in Table 1, taking λ2H ¼ 0:08 h . 1 Furthermore, we have increased α32 by a factor of 10 ðα32 ¼ 0:01 h Þ to give a more distinct peak within the second week period. Initial conditions are given by the steady state solution subject to setting f ðtÞ ¼ 0.

In this paper, we have restricted ourselves to examine the model under the restriction of constant antigen presentation. This restriction has been appropriate as we have not modelled the precise details of an immune response, which require a general function f(t). We conclude our discussion with numerical results that describe the time evolution, in the case of time-dependent specific antigen presentation, of the n2, n3 and n4 populations. In Fig. 9, we present example trajectories for non-constant antigen presentation. We first assume f ðtÞ ¼ 0 for seven days. Between days 7 and 14 we let f ðtÞ ¼ 1, and following day 14 we set f ðtÞ ¼ e 0:027ðt 1424Þ . That is, from day 14, f(t) is a decreasing exponential, which decreases from 1 to 1/100 by approximately day 21. Trajectories were computed using the parameter set given in Table 1; at day 0 we let initial conditions be equal to the steady state solution subject to setting f ðtÞ ¼ 0. This figure shows that there is only a modest increase in the number of IL-2 producing T cells. It seems that the main factor allowing growth of the total CD4 þ pool is not a large increase in the amount of available IL-2 (assumed to be proportional to the number of IL-2 producing T cells). We have assumed that recognition of antigen by the nonregulatory populations allows them to up-regulate IL-2Rα. It would appear that it is this extra responsiveness to the available IL-2, which drives growth of the population. A caveat to these observations is that we have made the simplification that IL-2

ð68Þ

diverges (Karlin and McGregor, 1957). A birth and death process with guaranteed absorption will be absorbed in a finite time if the series 1

∑ λn r n

ð69Þ

n¼1

converges (Karlin and McGregor, 1957). Let a birth and death process be defined by the transition rates

λn ¼ λne n=κ ;

ð70Þ

μn ¼ μn for n Z 0;

ð71Þ

with λ; μ and κ 4 0. We prove, in what follows, that for such a process, absorption at the (absorbing) state n ¼0 is certain and occurs in a finite time. Absorption is certain: For the birth and death process defined by Eqs. (70) and (71) we have λn 1 ðn 1Þ! exp ∑ni ¼ 11 i=κ rn ¼ : ð72Þ μn 1 n! Then, it follows that 1

λn rn

¼

μn 1 ; λn exp ∑ni¼ 1 i=κ

and by the ratio test we have μ exp ðn þ 1Þ=κ λn rn lim ¼ lim - þ 1: n- þ 1λn þ 1 r n þ 1 n- þ 1 λ

ð73Þ

ð74Þ

Thus, the series ∑1 n ¼ 1 1=λn r n diverges and hence the process is absorbed with certainty, as we wanted to show. Time to absorption is finite: For this process we can write λ exp ∑ni¼ 1 i=κ λn rn ¼ n ; ð75Þ n1

μ


and again by the ratio test we have

λn þ 1 r n þ 1 λ ¼ 0: ¼ lim lim n- þ 1 n- þ 1μ exp ðn þ 1Þ=κ λn r n

ð76Þ

λ converges and the process is absorbed Thus, the series in a finite time as we wanted to show. Naive T cells: The Markov sub-chain fX1 ðtÞg can be considered as a birth and death process with transition rates of the type given by Eqs. (70) and (71). For the population of naive T cells we have ∑1 n ¼ 1 n rn

λ ¼ λ1H ; μ ¼ μ1 þ α13 ; κ ¼ κ1 :

173

independently of the remaining three populations, the process is guaranteed to reach (in a finite time) region R1 defined as R1 ¼ fð0; n2 ; n3 ; n4 Þ : ni Z 0 for i ¼ 2; 3; 4g:

ð78Þ

Equivalently, we have shown that population fX2 ðtÞg will become extinct in a finite time, independent of the other populations. Thus, the process is guaranteed to reach (in a finite time) region R2 defined as R2 ¼ fðn1 ; 0; n3 ; n4 Þ : ni Z 0 for i ¼ 1; 3; 4g:

ð79Þ

Since both R1 and R2 are reached with certainty in a finite time, it follows that the intersection of these regions, R3, defined as R3 ¼ fð0; 0; n3 ; n4 Þ : ni Z 0 for i ¼ 3; 4g;

Thus, the population of naive T cells will become extinct with certainty in a finite time. IL-2 producing T cells: Borrowing ideas from Iglehart (Iglehart, 1964), we can bound the Markov sub-chain fX2 ðtÞg by a birth and death process (with transition rates of the form given in Eqs. (70) and (71)), which “travels to infinity” faster than fX2 ðtÞg (Iglehart, 1964). For such a process we have

λ ¼ λ2H ; μ ¼ μ2 þ α23 ; κ ¼ κ2 ; which is absorbed at n¼ 0 in a finite time. Since this process bounds the Markov sub-chain fX2 ðtÞg, it follows that the population of effector T cells will also become extinct with certainty in a finite time. Regulatory T cells: If sufficient time has passed for the population of effector T cells to become extinct, that is, n2 ¼ 0, the transition probabilities for the Markov sub-chain fX4 ðtÞg depend only on n4. Indeed, the only non-zero transition probability for the sub-chain fX4 ðtÞg is pn4 1;n4 ðΔtÞ ¼ ProbfX4 ðt þ ΔtÞ ¼ n4 1∣X4 ðtÞ ¼ n4 g ¼ μ4 n4 Δt þ oðΔtÞ:

ð77Þ Such a process is a simple death process and if not already at the unique absorbing state, n4 ¼ 0, before extinction of the effector T cell population, the process will be absorbed in a finite time (Karlin and McGregor, 1957). IL-2 non-producing T cells: Lastly, suppose sufficient time has elapsed for the naive and effector T cell populations to become extinct. Following extinction of these populations, the Markov sub-chain representing the IL-2 non-producing T cell population can be considered a birth and death process of the form (70) and (71) with

λ ¼ λ3H ; μ ¼ μ3 ; κ ¼ κ3 : It, therefore, follows that the IL-2 non-producing T cell population will become extinct in a finite time, given populations fX1 ðtÞg and fX2 ðtÞg have already become extinct. Let us suppose the process starts in a general state ðn1 ; n2 ; n3 ; n4 Þ A S, where ni Z 1 for i ¼ 1; …; 4. By virtue of the fact that population fX1 ðtÞg will become extinct in a finite time,

ð80Þ

is reached with certainty in a finite time. Given that enough time has elapsed for the populations of naive and effector T cells to have become extinct, we have shown that the population of IL-2 nonproducing T cells will become extinct in a finite time. Therefore, if the process is in R3, the region R4, defined as R4 ¼ fð0; 0; 0; n4 Þ : n4 Z 0g;

ð81Þ

is reached with certainty in a finite time. Equivalently, we have shown that following extinction of the population of effector T cells, the population of regulatory T cells will become extinct in a finite time. Thus, the region R5, defined as R5 ¼ fð0; 0; n3 ; 0Þ : n3 Z 0g;

ð82Þ

is reached with certainty in a finite time for a process starting in R3. Finally, considering that R4 and R5 are both reached with certainty in a finite time, the intersection of these regions, namely the absorbing state ð0; 0; 0; 0Þ is reached with certainty in a finite time from region R3. It, thus, follows the entire process becomes extinct in a finite time starting from an arbitrary state ðn1 ; n2 ; n3 ; n4 Þ A S, with ni Z 0 for i ¼ 1; …; 4. A.2. Probability and time to extinction of populations 2 and 4 We have shown that the populations of effector and regulatory T cells will both become extinct in a finite time. In this section we present a numerical method to approximately calculate the probability that the regulatory T cell population will be driven to the absorbing region defined by R ¼ fðn2 ; n4 Þ : n2 a 0; n4 ¼ 0g;

ð83Þ

before the population of effector T cells becomes extinct. In such an event, when the population of regulatory T cells becomes extinct before the effector T cell population, the CD4 þ T cell population is at risk of unregulated proliferation of the effector T cell subset, following a possible increase in their proliferative capacity, which for example could be induced by presentation of specific antigen. In this scenario, the probability of autoimmunity increases (Almeida et al., 2012). In this section, we let fXðtÞg ¼ fðX2 ðtÞ; X4 ðtÞÞg be the bi-variate Markov process representing populations 2 and 4. The transition probabilities for such a process are

pm;n ðΔtÞ ¼ Prob Xðt þ ΔtÞ ¼ m∣XðtÞ ¼ n 8 ðμ2 þ α23 þ βn4 Þn2 Δt þoðΔtÞ > > > > κ > > > μ4 4 n4 Δt þ oðΔtÞ > > > κ 4 þ n2 > > > < λ2H n2 e n2 =κ 2 Δt þ oðΔtÞ ¼ λ4 n2 n4 Δt þ oðΔtÞ > > > A

> > κ4 > > n4 þ λ2H n2 e n2 =κ 2 þ λ4A n2 n4 Δt þ oðΔtÞ > 1 ðμ2 þ α23 þ β n4 Þn2 þ μ4 > > κ 4 þ n2 > > > : oðΔtÞ

if m ¼ ðn2 1; n4 Þ; if m ¼ ðn2 ; n4 1Þ; if m ¼ ðn2 þ 1; n4 Þ; if m ¼ ðn2 ; n4 þ 1Þ; if m ¼ n; otherwise:

174


To calculate the probability that population 4 becomes extinct before population 2, we first impose that the state space of the Markov process fXðtÞg is finite, such that the maximum number of cells for population 2 is N2 and the maximum number of cells for population 4 is N4. Imposing a finite state space implies that the calculated probabilities will be an approximation to the true ones. Since the process typically does not grow to infinity (due to the carrying capacities of the cell subsets), for increasing values of Ni, i¼2, 4, the likelihood of any realisation (for a given parameter set) reaching either Ni, decreases. Thus, we can increase the accuracy of the calculation to arbitrary precision by increasing N2 and N4. This, however, comes at a further computational cost which we

κ4

κ 4 þ n2

if m ¼ n;

ð85Þ

Let Sn be the set of states which can be reached in a single transition from state n. Then Sn is given by Sn ¼ fðn2 1; n4 Þ; ðn2 þ 1; n4 Þ; ðn2 ; n4 1Þ; ðn2 ; n4 þ 1Þg;

ð86Þ

for 0 o n2 o N 2 and 0 o n4 oN 4 . When n2 ¼ N 2 , we have Sn ¼ fðN 2 1; n4 Þ; ðN 2 ; n4 1Þ; ðN 2 ; n4 þ 1Þg. When n4 ¼ N4 , we have Sn ¼ fðn2 1; N 4 Þ; ðn2 þ 1; N4 Þ; ðn2 ; N 4 1Þg, and when n2 ¼ N 2 , n4 ¼ N4 , we have Sn ¼ fðN2 1; N 4 Þ; ðN2 ; N4 1Þg. If the process reaches either the region defined by R or the region R defined by R ¼ fðn2 ; n4 Þ : n2 ¼ 0g;

ð87Þ

the process stops.

P n ∑ ρm;n P m ¼ s:

ð88Þ

m A Sn

= R. We define s as follows: when n A R, s ¼ 1, and s ¼ 0 if n2 We can express Eq. (88) in the form AQ ¼ a;

ð89Þ

where A is a ðN 2 N 4 Þ ðN 2 N 4 Þ matrix that depends on the variables

ρm;n , Q is a ðN 2 N 4 Þ 1 (column) vector that depends on the

variables P n , and a is a ðN 2 N4 Þ 1 (column) vector that depends on the variables s. In order to construct A, Q and a, we first define a bijection ψ from the set K ¼ f0; …; N 2 g f0; …; N 4 g to the set K ¼ f0; …; ðN 2 þ 1ÞðN 4 þ 1Þ 1g, such that an ordering for all states in K can be introduced. We note that the ordering of the states in K defined by the bijection is irrelevant to the calculation, therefore any choice of bijection ψ will suffice. Let A ¼ ðAi;j Þ, Q ¼ ðQk Þ and a ¼ ðak Þ, with i; j; k A K. We then define if i ¼ j;

Ai;j ¼ ρψ 1 ðiÞ;ψ 1 ðjÞ ak ¼ 1

if ψ

1

ð91Þ

m A S0n

where in is the mean time the process waits in state n. This waiting time is given by in ¼ 1=qn;n (Allen, 2003). Note that S0n is the set of all states that can be reached from state n in a single transition and that are not in Ω. For n A Ω we impose τn ¼ 0. We solve the above equation by setting it up as a linear algebra problem as we did before. Define a bijection φ from the set K ¼ f0; …; N 2 g f0; …; N4 g\Ω to the set K ¼ f0; …; ðN 2 þ 1ÞðN 4 þ 1Þ 1 jΩjg. Given i; j A K, let B ¼ ðBi;j Þ, b ¼ ðb0 ; b1 ; …; bðN2 þ 1ÞðN4 þ 1Þ 1 jΩj ÞT and T ¼ ðT 0 ; T 1 ; …; T ðN2 þ 1ÞðN4 þ 1Þ 1 jΩj ÞT , and define Bi;j ¼ 1

if i ¼ j;

Bi;j ¼ ρφ 1 ðiÞ;φ 1 ðjÞ bk ¼ iφ 1 ðkÞ ;

if φ 1 ðiÞ A S0φ 1 ðjÞ ;

T k ¼ τφ 1 ðkÞ ;

for i; j; k A K. The solution for the expected time to absorption into the set Ω is T ¼ B 1 b:

ð92Þ

References

A.2.1. Probability of extinction By a first step analysis argument, for which further details can be found in Allen (2003), we write

Ai;j ¼ 1

One immediate drawback of this method is the computational time it takes to invert the matrix A. As one increases the size of the state space, A increases with the square of N 2 N4 . It, therefore, follows that the computational time to invert A grows like ðN2 N4 Þ2 as N2 and N4 increase.

τn ∑ ρm;n τm ¼ in ;

and ρm;n as follows: qm;n : qn;n

ð90Þ

if m ¼ ðn2 ; n4 1Þ;

ð84Þ

ρm;n ¼

Q ¼ A 1 a:

if m ¼ ðn2 1; n4 Þ;

if m ¼ ðn2 þ 1; n4 Þ; if m ¼ ðn2 ; n4 þ 1Þ; n4 þ λ2H n2 e n2 =κ2 þ λ4A n2 n4

Since the matrix A has been defined so that Ai;i ¼ 1, A is nonsingular and thus, invertible. In this case, the numerical solution is

A.2.2. Expected time to extinction In a similar manner to the calculation of the probability of extinction, we can also compute the expected time to extinction of the populations of effector and regulatory T cells. Let τn be the expected time to reach the set of states Ω ¼ R [ ð0; 0Þ from n2 = Ω. Then, by a first step analysis argument we have (Allen, 2003)

discuss below. Let P n be the probability of reaching region R from an arbitrary state n ¼ ðn2 ; n4 Þ. By definition we have P ðn2 ;0Þ ¼ 1, P ð0;n4 Þ ¼ 0, P ð0;0Þ ¼ 0, and 0 oP ðn2 ;n4 Þ o 1, with n2 ; n4 4 0. We define 8 ðμ þ α23 þ βn4 Þn2 > > > 2 > κ > > > μ4 4 n4 > > > < κ 4 þ n2 qm;n ¼ λ2H n2 e n2 =κ2 > > > λ4A n2 n4 > > > > > > > : ðμ2 þ α23 þ βn4 Þn2 þ μ4

ak ¼ 0 if ψ 1 ðkÞ2 = R; Qk ¼ P ψ 1 ðkÞ :

if ψ 1 ðiÞ A Sψ 1 ðjÞ ;

ðkÞ A R;

Allen, L., 2003. An introduction to stochastic processes with applications to biology. Pearson/Prentice Hall, Upper Saddle River, New Jersey. Almeida, A., Amado, I., Reynolds, J., Berges, J., Lythe, G., Molina-París, C., Freitas, A., 2012. Front. Immunol. 3. Almeida, A., Rocha, B., Freitas, A., Tanchot, C., 2005. Semin. Immunol. 17 (3), 239. Almeida, A., Zaragoza, B., Freitas, A., 2006a. Int. Immunol. 18 (11), 1607. Almeida, A., Zaragoza, B., Freitas, A., 2006b. J. Immunol. 177 (1), 192. Almeida, A.R., Legrand, N., Papiernik, M., Freitas, A.A., 2002. J. Immunol. 169 (9), 4850. Antony, P., Piccirillo, C., Akpinarli, A., Finkelstein, S., Speiss, P., Surman, D., Palmer, D., Chan, C., Klebanoff, C., Overwijk, W., Rosenberg, S., Restifo, N., 2005. J. Immunol. 174 (5), 2591. Bains, I., Antia, R., Callard, R., Yates, A., 2009a. Blood 113 (22), 5480. Bains, I., Thiébaut, R., Yates, A., Callard, R., 2009b. J. Immunol. 183 (7), 4329. Banchereau, J., Pascual, V., O0 Garra, A., 2012. Nat. Immunol. 13 (10), 925. Belkaid, Y., Rouse, B., 2005. Nat. Immunol. 6 (4), 353. Berard, M., Tough, D., 2002. Immunology 106 (2), 127. Bettelli, E., Carrier, Y., Gao, W., Korn, T., Strom, T., Oukka, M., Weiner, H., Kuchroo, V., 2006. Nature 441 (7090), 235. Burroughs, N., Ferreira, M., Oliveira, B., Pinto, A., 2011. Math. Comput. Modell. 53 (7–8), 1389. Callard, R., Stark, J., Yates, A., 2003. Trends Immunol. 24 (7), 370. Carneiro, J., León, K., Caramalho, I., Van Den Dool, C., Gardner, R., Oliveira, V., Bergman, M.-L, Sepúlveda, N., Paixao, T., Faro, J., Demengeot, J., 2007. Immunol. Rev. 216 (1), 48. Chen, W., Jin, W., Hardegen, N., Lei, K., Li, L., Marinos, N., McGrady, G., Wahl, S., 2003. J. Exp. Med. 198 (12), 1875. Cheng, G., Yu, A., Malek, T., 2011. Immunol. Rev. 241 (1), 63. Diggle, S., Griffin, A., Campbell, G., West, S., 2007. Nature 450 (7168), 411. Edelstein-Keshet, L., 2005. Mathematical Models in Biology. SIAM.


Feinerman, O., Jentsch, G., Tkach, K., Coward, J., Hathorn, M., Sneddon, M., Emonet, T., Smith, K., Altan-Bonnet, G., 2010. Mol. Syst. Biol. 6 (1). Feuerer, M., Hill, J.A., Kretschmer, K., von Boehmer, H., Mathis, D., Benoist, C., 2010. Proc. Natl. Acad. Sci. 107 (13), 5919. Fontenot, J.D., Gavin, M.A., Rudensky, A.Y., 2003. Nat. Immunol. 4 (4), 330. Fouchet, D., Regoes, R., 2008. PLoS One 3 (5), e2306. Freitas, A., Rocha, B., 2000. Annu. Rev. Immunol. 18 (1), 83. Ge, Q., Hu, H., Eisen, H., Chen, J., 2002. Microbes Infect. 4 (5), 555. Iglehart, D., 1964. Ann. Math. Stat., 350–361. Janeway Jr., C.A., Travers, P., Immunobiology: Principles of Innate and Adaptive Immunity in the Immune System in Health and Disease, third edition. Current Biology Ltd, London, Garland Publishing, New York, USA, 2001. Jelley-Gibbs, D., Lepak, N., Yen, M., Swain, S., 2000. J. Immunol. 165 (9), 5017. Jenkins, M., Khoruts, A., Ingulli, E., Mueller, D., McSorley, S., Reinhardt, R., Itano, A., Pape, K., 2001. Annu. Rev. Immunol. 19 (1), 23. Josefowicz, S.Z., Lu, L.-F., Rudensky, A.Y., 2012. Annu. Rev. Immunol. 30, 531. Karlin, S., McGregor, J., 1957. Trans. Am. Math. Soc. 86 (2), 366. Kim, P., Lee, P., Levy, D., 2007. J. Theor. Biol. 246 (1), 33. Kimura, M.Y., Pobezinsky, L.A., Guinter, T.I., Thomas, J., Adams, A., Park, J.-H., Tai, X., Singer, A., 2012. Nat. Immunol. Kretschmer, K., Apostolou, I., Hawiger, D., Khazaie, K., Nussenzweig, M., Von Boehmer, H., 2005. Nat. Immunol. 6 (12), 1219. León, K., Lage, A., Carneiro, J., 2003. J. Theor. Biol. 225 (1), 107. León, K., Pérez, R., Lage, A., Carneiro, J., 2000. J. Theor. Biol. 207 (2), 231. León, K., Pérez, R., Lage, A., Carneiro, J., 2001. J. Immunol. 166 (9), 5356. Malek, T., Castro, I., 2010. Immunity 33 (2), 153.

175

Malek, T., Yu, A., Vincek, V., Scibelli, P., Kong, L., 2002. Immunity 17 (2), 167. Miller, M., Bassler, B., 2001. Annu. Rev. Microbiol. 55 (1), 165. Müller, A.J., Filipe-Santos, O., Eberl, G., Aebischer, T., Späth, G.F., Bousso, P., 2012. Immunity 37 (1), 147. Nikolich-Zugich, J., et al., 2004. Nat. Rev. Immunol. 4 (2), 123. Pepper, M., Jenkins, M., 2011. Nat. Immunol. 131 (6), 467. Powrie, F., Maloy, K., 2003. Science 299 (5609), 1030. Qureshi, O.S., Zheng, Y., Nakamura, K., Attridge, K., Manzotti, C., Schmidt, E.M., Baker, J., Jeffery, L.E., Kaur, S., Briggs, Z., et al., 2011. Science 332 (6029), 600. Sakaguchi, S., 2000. Cell 101 (5), 455. Sakaguchi, S., 2004. Annu. Rev. Immunol. 22, 531. Sakaguchi, S., Sakaguchi, N., Asano, M., Itoh, M., Toda, M., 1995. J. Immunol. 155 (3), 1151. Scollay, R.G., Butcher, E.C., Weissman, I.L., 1980. Eur. J. Immunol. 10 (3), 210. Seddon, B., Mason, D., 2000. Immunol. Today 21 (2), 95. Seddon, B., Tomlinson, P., Zamoyska, R., 2003. Nat. Immunol. 4 (7), 680. Sepúlveda, N., Carneiro, J., 2011. In: Mathematical Models and Immune Cell Biology. Springer, pp. 275–303. Tang, Q., Adams, J., Tooley, A., Bi, M., Fife, B., Serra, P., Santamaria, P., Locksley, R., Krummel, M., Bluestone, J., 2005. Nat. Immunol. 7 (1), 83. von Boehmer, H., 2005. Nat. Immunol. 6 (4), 338. Walker, L.S., Sansom, D.M., 2011. Nat. Rev. Immunol. 11 (12), 852. Wan, Y., Flavell, R., 2009. J. Mol. Cell Biol. 1 (1), 20. Wing, K., Onishi, Y., Prieto-Martin, P., Yamaguchi, T., Miyara, M., Fehervari, Z., Nomura, T., Sakaguchi, S., 2008. Science 322 (5899), 271.

Autoimmunity: a T-cell receptor perspective.

Impact of sepsis on CD4 T cell immunity.

A Mathematical Framework for Understanding Four-Dimensional Heterogeneous Differentiation of CD4+ T Cells.

Population mechanics: A mathematical framework to study T cell homeostasis.

Distinct CD4 T-cell effects on primary versus recall CD8 T-cell responses during viral encephalomyelitis.

Aging of the CD4 T Cell Compartment.

STAT5 and CD4 + T Cell Immunity.

Plasticity of human CD4 T cell subsets.

Unbalanced Neonatal CD4(+) T-Cell Immunity.

CD4 T-cell memory generation and maintenance.

Psoriasis associated with idiopathic CD4+ T-cell lymphopenia: a regulatory T-cell defect?

Interleukin-7 signalling defects in naive CD4+ T cells of HIV patients with CD4+ T-cell deficiency on antiretroviral therapy are associated with T-cell activation and senescence.

CD4+ and CD8+ T cell activation are associated with HIV DNA in resting CD4+ T cells.

CD8+ memory T-cell inflation renders compromised CD4+ T-cell-dependent CD8+ T-cell immunity via naïve T-cell anergy.

Dynamic interaction between CD4+ T cells and parasitic helminths: mathematical models of heterogeneity in outcome.

Discrimination of human CD4 T cell clones based on their reactivity with antigen-presenting T cells.

CD4-T-cell antigen receptor complexes on human leukemia T cells.

Helping themselves: optimal virus-specific CD4 T cell responses require help via CD4 T cell licensing of dendritic cells.

Mathematical modeling of oncogenesis control in mature T-cell populations.

CD28 co-stimulation in T-cell homeostasis: a recent perspective.

Cell-composition effects in the analysis of DNA methylation array data: a mathematical perspective.

or cytoplasmic tail.

delta+ T lymphocytes.

Longitudinal requirement for CD4+ T cell help for adenovirus vector-elicited CD8+ T cell responses.