A MULTISCALE MODEL OF THE CD8 T CELL IMMUNE RESPONSE STRUCTURED BY INTRACELLULAR CONTENT

. During the primary CD8 T cell immune response, CD8 T cells undergo proliferation and continuous diﬀerentiation, acquiring cytotoxic abil- ities to address the infection and generate an immune memory. At the end of the response, the remaining CD8 T cells are antigen-speciﬁc memory cells that will respond stronger and faster in case they are presented this very same antigen again. We propose a nonlinear multiscale mathematical model of the CD8 T cell immune response describing dynamics of two inter-connected physical scales. At the intracellular scale, the level of expression of key proteins involved in proliferation, death, and diﬀerentiation of CD8 T cells is modeled by a delay diﬀerential system whose dynamics deﬁne maturation velocities of CD8 T cells. At the population scale, the amount of CD8 T cells is represented by a discrete density and cell fate depends on their intracellular con- tent. We introduce the model, then show essential mathematical properties (existence, uniqueness, positivity) of solutions and analyse their asymptotic behavior based on the behavior of the intracellular regulatory network. We numerically illustrate the model’s ability to qualitatively reproduce both pri- mary and secondary responses, providing a preliminary tool for investigating the generation of long-lived CD8 memory T cells and vaccine design.


(Communicated by Kevin Painter)
Abstract. During the primary CD8 T cell immune response, CD8 T cells undergo proliferation and continuous differentiation, acquiring cytotoxic abilities to address the infection and generate an immune memory. At the end of the response, the remaining CD8 T cells are antigen-specific memory cells that will respond stronger and faster in case they are presented this very same antigen again. We propose a nonlinear multiscale mathematical model of the CD8 T cell immune response describing dynamics of two inter-connected physical scales. At the intracellular scale, the level of expression of key proteins involved in proliferation, death, and differentiation of CD8 T cells is modeled by a delay differential system whose dynamics define maturation velocities of CD8 T cells. At the population scale, the amount of CD8 T cells is represented by a discrete density and cell fate depends on their intracellular content. We introduce the model, then show essential mathematical properties (existence, uniqueness, positivity) of solutions and analyse their asymptotic behavior based on the behavior of the intracellular regulatory network. We numerically illustrate the model's ability to qualitatively reproduce both primary and secondary responses, providing a preliminary tool for investigating the generation of long-lived CD8 memory T cells and vaccine design.
1. Introduction. We aim, in this paper, at providing a novel continuous multiscale model of the primary and secondary CD8 T cell immune responses, able to describe the kinetics of a classical response against an intracellular pathogen [1] at both molecular and cellular scales, but also able to run fast and efficiently while accurately describing the various processes involved in the immune response. This will allow to investigate the effects of early molecular events on long-term predictions and the generation of immune memory. Indeed, an efficient primary specific immune response should lead to the generation of long-lived memory cells that allow a faster and more efficient response against the same pathogen [32] should a secondary infection occur: vaccines are based on the generation of immune memory. The generation of an efficient immune memory strongly depends on molecular signaling during the early stages of the response, it is then necessary to couple molecular and cellular dynamics within a model of the CD8 T cell immune response in order to investigate the development of memory cells and eventually optimize vaccine design.
Before antigen presentation by Antigen-Presenting Cells (APC), CD8 T cells are in a naive state. Following antigen presentation, naive T cells become activated and start proliferating through a phase of clonal expansion [44]. While proliferating, CD8 T cells differentiate into effector cells, which are functionally characterized by the acquisition of cytolytic abilities [6]. Effector T cells actively fight the infection by eliminating infected cells [6,44]. After elimination of the pathogen a contraction phase drastically reduces the size of the T cell population through apoptosis mechanisms [44]. A surviving fraction approximately equal to 5%-10% of the number of T cells at the peak of the response (around days 7-10 post-infection) is usually assumed to survive the contraction phase and results in a long-lived memory cell population [44].
Modeling of the CD8 T cell immune response attracted much attention over the past two decades. Most models that have been proposed mainly investigate the evolution of CD8 T cell counts over time and do not focus on the interactions between molecular and cellular scales. We mention the pioneering work of De Boer et al. [20] that identified proliferation, death, and differentiation rates of a CD8 T cell population throughout a response to LCMV virus in mice, using linear differential equations models. This work led to subsequent improvements [25,35,49] in the description and understanding of the immune response. It is however clear that the immune response is not strictly dependent on pathogen amount, since the peak of the response and the switch to the contraction phase do not correspond exactly to the elimination of the pathogen [4,33,34,52]: it has been observed that even with a brief pathogen encounter, T cells perform a complete response. In agreement with these findings, Antia et al. [4,5] proposed a modified model accounting for an age structure of the effector cell population and modeled a programmed proliferative response after pathogen encounter. Inspired by this model, Terry et al. [54] developed a model of the primary response to an acute infection based on a controlled regulation of most cellular processes that proved to be able to reproduce CD8 T cell kinetics observed in the literature [44]. This model has also recently been used [16] to identify the CD8 T cell differentiation pathway involved in CD8 T cell responses.
Identification of differentiation pathways involved in the CD8 T cell response, based on relevant markers and allowing to functionally define CD8 T cell subsets, has led to several propositions, sometimes contradictory ones. A diversity of effector and memory populations have been identified [6,14,51]. Memory cells have for instance been classified into two main categories: central memory T cells (TCM, CD62L+) and effector memory T cells (TEM, CD62L-) [42,51], and their precursors have been identified by several works from Kaech's group [18,34,46]. Links between the various phenotipically-defined subsets obtained over the years are still subject to debate because of the markers used that, most of the time, do not define neither a genealogy nor specific functions. Recently, Buchholz et al. [12] suggested that naive cells would differentiate into precursors of TCM, then into TEM precursors, before the effector stage, which is not in agreement with the classical representation where effector cells appear before memory cells. A branching process between effector and memory populations is also supported by a single cell transcriptome approach [7]. These controversial results could result from a subset definition which may identify functionally diverse entities rather than developmental stages. Crauste et al [16] showed that phenotypically defined subsets of CD8 T cells can represent functional classes and proposed the use of Ki67 and Bcl2, together with CD44, to identify CD8 T cell differentiation stages. They could identify memory cells with the CD44+Mki67-Bcl2+ phenotype, and additionally identified two stages of effector cells, early CD44+Ki67+Bcl2-and late CD44+Ki67-Bcl2-effector cells, while naive cells are CD44-Ki67-Bcl2+. The linear differentiation pathway consistsing in naive→early effector→late effector→memory cells has been shown to correctly fit several experimental data, contrary to Buscholz et al [12] model, and to extend previous results [4,5,25].
In order to account for molecular signaling leading to CD8 T cell differentiation and eventually to the generation of memory cells, the above-mentioned models that only describe cell population kinetics are not sufficient. Although some works [3,21,22,31] focused on the modeling of the intracellular molecular scale, to describe either apoptosis, survival, differentiation or proliferation, they highlight the role of molecular events independently from cellular and cell population events.
In order to reach a better and more complete description of the immune response, multiscale models have been proposed: they tend to describe a biological process at several -usually two -scales simultaneously (molecular and cellular scales, cell population and organ scales, etc.). Several multiscale models of the immune response have been proposed in recent years [9,13,19,37,38,56], yet some scales of interest are often not explicitly modeled in these works, for instance by considering that dynamics at one scale is at equilibrium.
Boosted by the improvement in computational power during the last decade, computational multiscale models of the immune response have emerged and have proven themselves useful to understand the many interactions between scales occurring during the immune response [2,8,9,13,37,38,43,48,57]. Most of these computational models allow stochasticity to be observed since all "individuals" obey the same rules but are not necessarily in the same medium or subject to the same external stimuli, yielding variability in the outputs of the models. It is however noticeable that these models does not usually describe intracellular molecular regulation and therefore cannot be used to investigate the consequences of molecular events on the cell population dynamics.
Prokopiou et al. [47] and Gao et al. [26] recently developed a hybrid multiscale model of the CD8 T cell immune response, coupling descriptions of both molecular regulatory mechanisms (through a minimal molecular regulatory network) and cellular dynamics. These latter are based on a description of the activation and differentiation processes in the early stages of the immune response (with an explicit description of spatial interactions). Molecular scale dynamics are based on an ordinary differential equation (ODE) model describing the regulatory network dynamics, within every single CD8 T cell. At the cellular scale, described by a cellular Potts model (agent-based model), differentiation and proliferation rates of CD8 T cells are defined as functions of the level of expression of the molecular players. In this multiscale model, physical interactions between cells are taken into account, as well as diffusion of a cytokine (interleukin-2) in the extracellular environment.
This mathematical and computational multiscale model is able to qualitatively reproduce immune cell dynamics and experimental results on the early phases of the immune response [26]. However, the complexity of this model and of similar ones (large number of interactions and processes, several cell types, etc.) comes with a high computational cost and consequently some limitations of its initial scope.
Using structured transport equations [45], Friedman et al. [23,24] provided a multiscale model of the CD4 T cell immune response, where CD4 T cells evolve according to a transport equation structured by the level of expression of two transcription factors whose dynamics depend on an ODE system. The authors showed that the asymptotic behavior of the intracellular model strongly influences the asymptotic behavior of the population model: the density of cells concentrates around Dirac masses located on the asymptotically stable steady states of the intracellular system of ODE. Even though this model deals with different types of cells, it can be viewed as a continuous counterpart of the hybrid model in Prokopiou et al. [47] and Gao et al. [26] regarding its design and its scope.
In this work, we introduce a continuous multiscale model of the CD8 T cell immune response, validate its biological relevance by performing a stability analysis, and show that it is numerically able to qualitatively describe the kinetics of a classical response against an intracellular pathogen. To do so, we will implement, at the cell population scale, the linear differentiation pathway commonly admitted in the literature [4,16,20,25] and mentioned above. At the intracellular scale, we will use a simplified network based on two markers, identified in Crauste et al. [16], the nuclear protein Ki67 and the protein Bcl2. Although too simple to account for all the necessary molecular mechanisms involved in a CD8 T cell immune response, this molecular model will show the feasibility of our approach and illustrate the expected behavior of such a multiscale model.
In Section 2, we will introduce this model, in which CD8 T cells will be represented in a two-dimensional maturity space, where each maturity variable will correspond to one of the two elements of the molecular network. Differentiation stages will then be associated with levels of expression of the maturity variables. At the intracellular scale, evolution of maturity variables will be modeled by nonlinear differential equations, and the fate of each cell within the cell population will depend on its position in the maturity space. We will show, in Section 3, that this model is mathematically well-posed and has the good theoretical properties to simulate an immune response, in particular the asymptotic behavior analysis of the intracellular model and its impact on the population scale will be considered. In Section 4, we will show that the model efficiently describes primary and secondary CD8 T cell immune responses at the molecular and cell population scales, qualitatively, and runs much faster than equivalent agent-based models. This model will prove to be simple enough to be studied mathematically and generalized to larger molecular networks or other immune responses if needed. 2. Multiscale model. This section is devoted to the presentation of the multiscale model and its components. We first introduce the CD8 T cell subsets that will be considered in the model. Then, we present a single cell model which describes the intracellular dynamics of each CD8 T cell. Finally, we describe a population model based on the single cell model, where cells interact, differentiate, die and proliferate depending on their intracellular content.
2.1. Cellular subpopulations. Based on Crauste et al. [16], we consider that five CD8 T cell subpopulations appear throughout an immune response, with different capabilities regarding proliferation, survival and cytotoxic abilities (see Figure 1). Their identification is based on the expression of CD44, Ki67, and Bcl2 (see Section 2.2).
Naive cells are quiescent CD8 T cells which have not encountered APC yet. We assume neither death nor proliferation for these cells, only activation by APC that induces their differentiation in activated cells [1,47]. As a consequence, and for the sake of simplicity, dynamics of naive cells are not explicitly described in our model: the naive cell population represents a reservoir of cells emptied by the APC-induced activation.
Activated cells correspond to the first stage of differentiation of CD8 T cells. They do proliferate, but do not have cytotoxic abilities yet. They represent a transient state between naive and effector cells, and have a relatively small death rate [47].
Effector cells result from the differentiation of activated cells, and are characterized by acquired cytotoxic abilities [1,16,47]. They also proliferate and exhibit a high death rate. They differentiate into late effector cells [16].
Late effector cells still have cytotoxic abilities but have lost their proliferation potential [16]. They are mainly observed after the peak of the response, during the contraction phase, hence they are characterized by a high death rate. They eventually generate the memory cell population [16,36].
Memory cells are long-lived CD8 T cells (TCM [42,51]) that do not proliferate on short time scales and have a very small death rate. They are almost quiescent cells waiting for reactivation through APC stimulation to start a secondary antigen specific response [1,16].
For these subpopulations, linear and branching differentiation schemes have been confronted with biological data of a primary immune response to influenza virus infected mice, and a model selection procedure pointed towards a linear scheme [16]. Following these conclusions, we assume in our model a linear differentiation scheme (see Figure 1).
Besides the CD8 T cell populations, we also consider the pathogen count -whose quantity at time t is denoted P (t) and satisfies (4) (see Section 2.3) -and, implicitly, the dynamics of APC. The pathogen interacts with CD8 T cells via APC when the immune response starts. Following activation of the response, the pathogen enhances CD8 T cell survival and proliferation during the clonal expansion phase, via the release of pro-survival cytokines for instance [42]. Then effector cells eliminate infected cells, contributing to the elimination of the pathogen. We assume that pathogen evolution is characterized by both a natural clearance rate and a cytotoxic cell-induced clearance rate. Proliferation of the pathogen is neglected since we aim at describing an immune response triggered by vaccination, and thus involving weakened, non-proliferative pathogens. A summary of the main variables of the model is given in Table 1. 2.2. Single cell dynamics. We consider here a single CD8 T cell. Except for the naive CD8 T cell population, the other four CD8 T cell subpopulations (activated, early effector, late effector, memory) mainly distinguish themselves by their ability to survive and proliferate, characterized by the expressions of two relevant biological markers, identified in Crauste et al. [16]: Ki67 that is a proliferation marker [50] and Bcl2 that is a survival marker [21]. We justify hereafter the construction of a

Variables Description
N (t) Naive CD8 T cell count at time t P (t) Pathogen count at time t E(t) Early and Late Effector CD8 T cell count at time t µ 1 (t) Level of expression of Ki67 at time t µ 2 (t) Level of expression of Bcl2 at time t ρ(t, µ 1 , µ 2 ) Density of CD8 T cells with maturities (µ 1 , µ 2 ) at time t v 1 Maturation velocity of µ 1 v 2 Maturation velocity of µ 2 molecular regulatory network based on Ki67 and Bcl2 expressions that will allow a multiscale description of CD8 T cell dynamics. Crauste et al. [16] showed that among the Bcl2-effector cell population, one could distinguish Ki67+ early effector cells and Ki67-late effector cells, while memory cells are Ki67-Bcl2+ cells. Both effector and memory cells are CD44+ cells while naive cells are CD44-cells. Since we do not describe naive cell molecular dynamics (see Section 2.1), only CD44+ cells molecular dynamics will be described, and hence it is only necessary to focus on Ki67 and Bcl2 to identify each CD8 T cell subset. It is worth mentioning that activated cells are not described as a developmental stage in Crauste et al [16], yet they represent a transient state between the Ki67-Bcl2+ naive cells and the Ki67+Bcl2-effector cells. For the sake of clarity, we then introduce the activated cells as Ki67+Bcl2+ cells, since after activation by APC naive cells enter a phase of strong proliferation (Ki67+) and survival (Bcl2+). The simplified network regulating the interactions between Ki67 and Bcl2, and involving the stimulation of a single CD8 T cell by an APC in our model, is given in Figure 2.
Due to a lack of information in the literature regarding interactions between Ki67 and Bcl2, we started from the molecular regulatory network developed in Prokopiou et al. [47] to build the model in Figure 2.
We first assumed variations in Ki67 expression would be qualitatively similar to IL-2 variations, since Ki67 characterizes proliferation and IL-2 is the main cytokine driving CD8 T cell proliferation. We then assumed that APC induce the expression of Ki67 via interactions with the TCR of CD8 T cells. Second, IL-2 was shown to be both activated by its own expression through a self-activating loop involving the IL-2 gene and its receptor, and inhibited via a negative feedback loop involving Blimp-1, this inhibition occurring on a slower time scale [41]. We then assume that Ki67 acts both positively and negatively on its own expression, with a slow inhibition and a fast activation.
Interactions of proteins Bim, Bcl2, Bax and caspases described in Ewings et al. [21] show that Bcl2 is a Bax and caspase inhibitor (therefore an apoptosis inhibitor and pro-survival marker), and stimulation by an APC has a positive effect on Bcl2 expression through Erk and Bim pathways. We then assume an APC-mediated activation of Bcl2. Since naive and memory cell populations have a high expression level of Bcl2 [16], we assume that Bcl2 is able to activate its own expression when under no other stimulus. In addition, we assume that exhaustion from extensive proliferation drives CD8 T cells towards apoptosis, and we consider that Ki67 negatively regulates Bcl2 expression.
We introduce the levels of expression, hereafter referred to as "maturity" levels, µ 1 and µ 2 which correspond respectively to Ki67 and Bcl2 expression levels of a CD8 T cell. Dynamics of µ 1 and µ 2 , schematically represented in Figure 2, are Table 2. Definition of the CD8 T cell subsets as functions of maturity levels.

Cell Subset
Domain Definition Activated cells given by the following system, with all parameters positive. We define the respective maturation velocities of µ 1 and µ 2 , by v 1 = dµ 1 /dt and v 2 = dµ 2 /dt. The use of first order Hill functions in (1) illustrates interactions by contact either between CD8 T cells and APC or proteins and their receptors. The effect of exhaustion is represented by the delay term µ 1 (t − τ ) [28,40]. This quantity is then put into a second order Hill function to input a saturation effect and keep its impact low for cells that did not proliferate much (lower values of µ 1 (t − τ )). The term for Bcl2 self-activation is represented by a logistic law, so that a stable steady state for µ 2 exists for r 2 K 2 > k 2 , when µ 1 = 0 and P = 0, corresponding to high expression levels. The negative impact of Ki67 on Bcl2 is described, in the simplest way, by a mass action law. Ki67 and Bcl2 degradation rates are assumed to be constant, given respectively by k 1 and k 2 .
To determine to which state of the differentiation process a cell belongs to (see Section 2.1) based on its intracellular content, we define two thresholds α and β for maturity values, that separate the cell subpopulations as defined in Table 2.
Next we are going to embed this single cell model into a model that describes the dynamics of the CD8 T cell population.

CD8 T cell population model.
We now introduce the model describing the evolution of a population of CD8 T cells responding to a stimulation by APC. Since naive cell dynamics are not described (Section 2.1), the model focuses on the dynamics of activated, early effector, late effector and memory cells, as well as pathogen count.
Activated, early effector, late effector and memory cell counts are represented by a weighted sum of Dirac masses, depending on time t and position (µ 1 , µ 2 ) in the maturity space. Each cell subset is defined by maturity levels (see Table 2). Since the model is deterministic, cells with the same maturity level behave identically, and are then represented by a single Dirac mass.
Let  . Illustration of the location process of newly activated cells at t = t k with n = 4. The pink disk represents the support of the source of activated cells. The density of newly activated cells at time t k is approximated by n 2 piles of cells. The black disc in each square represents the location of the corresponding pile of cells and its surface represents its weight. That weight is proportional to the intersection between the square and the pink disc.
upon stimulation of naive cells by APC. One may note that the appearance of new activated cells will depend both on the number of remaining naive cells, and on the presence of pathogen that could stimulate them via APC.
We discretize the Ω A domain to account for the generation of a heterogeneous population of activated cells following APC stimulation.The newly produced activated cells are located in a regular square of center (x m , y m ) and of area 4σ 2 (see Figure 3). This square is actually split into n 2 smaller squares, denoted squares (i, j), i, j ∈ {1, . . . , n}, of area h 2 , with h = 2σ/n. Each square (i, j) contains one pile of activated cells, located in its center ( . Hence, at time t k , n 2 Dirac masses may appear in Ω A , with different maturity values. These Dirac masses represent the pool of newly activated cells, which is a heterogeneous cell population in terms of µ 1 and µ 2 expression levels.
The amount of cells with maturity µ 1 and µ 2 is then described with a discrete density ρ, defined for t ∈ [0, T ] and µ 1 , µ 2 > 0, by where δ X i,j,k (t) is the Dirac mass accounting for cells whose intracellular content is The weight ω i,j,k (t), (i, j, k) ∈ {1, . . . , n} 2 × {1, . . . , q}, which gives the quantity of cells with maturities (µ i,j,k with P the amount of pathogen, E the amount of effector cells, given by and χ Ω the indicator function of Ω. Biologically, in the case of a primary immune response, before stimulation of naive cells by APC no differentiated CD8 T cell is observable. Hence, initially all the weights are null. Upon activation of n 2 naive cells at time t k , the weights ω i,j,k , i, j ∈ {1, . . . , n}, become positive to describe the presence of new activated cells. The initial location of the pile of cells δ X i,j,k (t) does not depend on k (see above). The initial weights ω i,j,k (t k ), k ∈ {1, . . . , q}, are proportional to the surface of the disc of center (x m , y m ) and of radius σ within the square (i, j) (see Figure 3). The proportional constants α i,j in (2) are given by with The constants α i,j have been normalized so that n i,j=1 α i,j = 1, and therefore, for k ∈ {1, . . . , q}, This indicates that the total quantity of cells activated at time t k is a correct approximation of the number of naive cells which differentiated between times t k−1 and t k .
In System (2), the rate F is defined by and accounts for the proliferation (first term on the right hand side of the equality) and death (second term on the right hand side of the equality) rates of the different subpopulations. We assume that F is bounded with respect to all variables.
The proliferation rate f div is assumed to depend on the intracellular concentration of Ki67 (µ 1 ), and to be an increasing function of µ 1 since Ki67 characterizes proliferation capacity. In addition, this rate is null for memory and late effector cells (µ 1 < α).
Regarding the death rate, we account for three biological ways for immune cells to die. First, CD8 T cells die due to their lack of survival capabilities. The corresponding death rate, f d (µ 2 ), is assumed to be a decreasing function of the pro-survival protein µ 2 . Second, cells die due to cell-cell contact through Fas-Fas ligand (FasL) interactions [10,53,58]. This physical interaction between a CD8 T cell expressing the membrane marker Fas and a neighboring CD8 T cell expressing the membrane receptor FasL is detrimental to the survival of the former. The rate f F as (E) is then assumed to be an increasing function of the number of early-and late-effector cells E (which both express the ligand of Fas, while all immune cells express Fas). This process acts as a negative feedback at the population scale. Natural apoptosis and Fas-induced death act on different pathways [30], hence the latter process doesn't depend on µ 2 . Finally, we consider death by Activation-Induced-Cell-Death (AICD): after activation, cells who encounter an APC again suffer from AICD and are driven towards apoptosis, unless their survival rate is high enough. This rate f AICD (P, µ 2 ) is assumed to be an increasing function of P and a decreasing function of µ 2 [11,15,27].
Evolution of maturity variables at the intracellular scale must also be discussed.
. . , q}, evolve in the (µ 1 , µ 2 )-plane following the set of differential equations The expected behavior of the model regarding the movement of the cells in the maturity space is displayed on Figure 4. Finally, the ODE system regulating the evolution of P and N , the naive CD8 T cell count, are assumed to be with initial conditions P (0) = P 0 > 0 and N (0) = N 0 > 0, and all parameters positive. Death of naive cells is neglected since they are slowly renewing, and consequently their evolution only depends on their stimulation by APC, which we describe with a first order Hill function to model a contact-induced saturation phenomenon. Similarly, the pathogen count decreases due to its natural clearance with a rate k P and through deadly attacks of infected cells by effector cells. Pathogen's proliferation is neglected in this work, given that we model an immune response triggered by vaccination, thus weakened pathogens, but could easily be added to the pathogen dynamics.
In the next section, we focus on basic mathematical properties of System (2)-(4), and we investigate the asymptotic behavior of solutions to System (1).
3. Mathematical analysis of the model.  Let T > 0 be fixed, as well as ∆t ∈ R * + and q ∈ N * such that T = q∆t, and n ∈ N * such that the domain Ω A contains a regular square of area 4σ 2 split into n 2 smaller squares (see Figure 3 and previous section).
Next, we focus on studying the asymptotic behavior of solutions to System (1) to investigate the asymptotic behavior of CD8 T cell maturity.
3.2. Steady states of the intracellular system. We focus on the existence of positive steady states of (1) driving the fate of a single CD8 T cell. The objective is to evaluate the potential asymptotic maturities of CD8 T cells after a primary or a secondary response. Such analytical results have been proven and observed by Friedman et al [23,24] in the case where the cell population is described by a continuous density whose maturation velocities are defined by an ODE system.
It is straightforward that P (t), solution of (4), satisfies P (t) < exp(−k P t) for all t > 0 and consequently tends to zero asymptotically. A steady state (x, y) of System (1) then satisfies From the second equation, we deduce the values of y from the values of x, since either y = 0, or From the first equation, we rule out x = 0 which is always a solution. We have left Let and Q(X) = aX 3 + bX 2 + cX + d, so that (6) is equivalent to Q(x) = 0. We first rule out the case where a = 0. In this case, c > 0 and Q is a second order polynomial whose roots are strictly negative.
To handle the case where Q is a third degree polynomial (a = 0), we use Sturm's theorem [55] to count the strictly positive roots of Q. Sturm's sequence for Q is given by P 0 = aX 3 + bX 2 + cX + d, For clarity purposes, we introduce so that P 2 = AX + B and First, let us state some properties regarding the constants we have just introduced.
Then, we give the following result. Proof. If we denote S(x) = (P 0 (x), P 1 (x), P 2 (x), P 3 (x)) and V (x) the number of sign changes in this sequence, Sturm's theorem states that the number of distinct roots in an interval [δ 1 , δ 2 ], with δ 1 < δ 2 , is given by V (δ 1 ) − V (δ 2 ). In our case, we are looking for positive roots, and therefore we count the sign changes in S(0) and S(∞). We have S(0) = (d, c, B, P 3 ), and the sign of S(x) when x → ∞ is given by the sign of (a, a, A, P 3 ). We try out all the different cases and we rule out the ones where the properties stated in Lemma 3.2 and the fact that b < 0 and d < 0 do not hold.
Eventually, Q has two strictly positive roots if and only if a < 0, c > 0 and P 3 < 0.
We do the same with the cases for 3, 1 or 0 positive steady states.
Let us mention that y, given by (5), is positive if and only if In the end, there can be up to three positive solutions (x, y) to System (5). Now, we investigate the stability of these steady states.
3.3. Local asymptotic stability. We investigate the stability of the positive steady states of System (1). To do so, we look for solutions of the characteristic equation [39] related to the autonomous following system Let define f (x) = γ r1 x θ r1 + x and g(x) = k 1 + γ l1 x 2 θ l1 + x 2 . System (10) can have up to three positive steady states (see Proposition 1). Let's assume throughout this section that (x, y) is a steady state of System (10), with x > 0 and y > 0. Then we have: Proposition 2. Assume System (10) has three positive steady states (x 1 , y 1 ), (x 2 , y 2 ), (x 3 , y 3 ), such that x 1 < x 2 < x 3 . Then (x 1 , y 1 ) and (x 3 , y 3 ) are unstable for all values of τ . Moreover, (x 2 , y 2 ) is stable if and only if τ < τ l , where Proof. We linearize System (10) and we write the corresponding characteristic equation, which we can simplify using (11) to obtain (λ − f (x)x + g (x)xe −λτ ))(λ + r 2 y) = 0, with (x, y) independent of τ .
The second term gives strictly negative eigenvalues, hence it won't affect the stability and we focus only on solving the following equation, Such a characteristic equation has already been studied. Necessary and sufficient conditions for stability have been given by Hayes [29], but won't provide us with an explicit critical value of τ at which a change of stability occurs. For τ = 0, the only solution of (12) is real and has the same sign as f (x) − g (x), which is positive for x 1 and x 3 , and negative for x 2 . We obtain that for τ = 0, (x 2 , y 2 ) is locally asymptotically stable and the other steady states are unstable. We now look for purely complex eigenvalues whose existence is a necessary condition to obtain a change in stability [39]. We substitute λ = iω into (12) to obtain, f (x) = g (x)cos(ωτ ), and ω = g (x)xsin(ωτ ), and necessarily, Regarding (x 1 , y 1 ) and (x 3 , y 3 ), a purely imaginary solution would be required for the branches of eigenvalues to cross the imaginary axis to become stable. But for those steady states, 0 < g (x) < f (x), and therefore, this situation cannot happen.
As an illustration, Figure 5 displays all the eigenvalues of (12) for τ ∈ [0, 30]. Existence of purely imaginary solutions to (12) means that periodic solutions to (1) emerge for τ = τ l . In fact, we numerically observe that this limit cycle even persists past the bifurcation point, with an amplitude that is an increasing function of τ (not shown).
In their work, Friedman et al. [23,24] proved that the density of CD4 T cells was converging towards a sum of Dirac masses located at the intracellular steady states. In our case, we have shown that there could exist a locally asymptotically stable positive steady state for τ small enough and can expect our discrete density to similarly converge towards a sum of Dirac masses. In the case where the delay τ is large enough, this steady state loses its stability. The solution can then either converge towards one of the other stable steady states (of type (0, y) or (x 2 , 0), not detailed here) or be asymptotically located along the limit cycle. We illustrate those properties in the next section with numerical simulations.

4.
Simulating primary and secondary immune responses. As mentioned in the previous section, it is not possible to theoretically access some properties due to the complexity of our model. Also, as is widely the case in modeling, no analytical solution of the model can be obtained. Hence, numerical methods are necessary to simulate the model and allow a deeper, biologically-motivated analysis. We hereafter present the numerical method used to solve our model.
Let T > 0, ∆t ∈ R + , q ∈ N, n ∈ N be fixed, such that T = q∆t, and the domain Ω A contains a regular square of area 4σ 2 split into n 2 smaller squares (see Figure 3).
To numerically solve System (2)-(4), we use an explicit Euler method with time step δt, such that τ /δt ∈ N * and ∆t/δt ∈ N * . The latter property allows to manually handle the discontinuities in (2).
The solution ρ of System (2)-(4) is expressed as a discrete sum of Dirac masses. In order to display a smooth approximation of ρ, we convolute ρ with a rescaled cut-off function ϕ, so that ε .
For our simulations, we chose ϕ(x, y) = 1 I e The computation of ρ ε is useful only to be able to observe the behavior of ρ, the total density of CD8 T cells, as a smooth density. Subpopulation dynamics can be obtained by summing the weights ω i,j,k (t) over the appropriate domains.
Regarding the rate F in (2), we hereafter motivate the choice of the 4 different functions (f div , f d , f F as and f AICD ) that compose it.
We define the division rate f div (µ 1 ) as We use a first order Hill function to account for the biological boundedness of a division rate. This function is shifted to equal zero for µ 1 < α (i.e., Ki67-cells do not proliferate) and is positive only for µ 1 > α (i.e. for Ki67+ cells). Consequently, memory and late effector cells do not proliferate.
The natural death rate f d is given by a second order Hill function, adjusted to have high values for effector and late effector cells and low values for activated and memory cells (using the threshold β as a parameter).
The rate f F as is defined as It is a first order Hill function describing the fact that Fas-induced cell death occurs through cell-cell contacts, and is consequently associated with a saturation effect. Finally, the AICD rate f AICD is described by a first order Hill function in P , weighted by a term that tempers this effect for high values of µ 2 , We hereafter present two simulations (Figures 6 and 8) of System (2)-(4) that show a qualitatively correct behavior of the model at both scales. All parameter values are given in Table 3, and have been calibrated to qualitatively reproduce a typical CD8 T cell immune response. Figure 6 displays the convoluted density ρ ε at different time steps of a primary and a secondary immune response, in the case where System (1) has a stable steady state located in Ω M .
In Figure 6, activated cells, produced by the stimulation of naive cells, appear in the Ω A domain (top right part of the maturity space) and start 'moving' following the linear differentiation scheme from Figure 1 (see Figures 6.(a) to 6.(c)). The intracellular system behaves as expected on Figure 4: the density ρ tends towards a weighted Dirac mass located at the intracellular steady state similarly to the results shown in Friedman et al. [23,24]. This Dirac mass represents the pool of quiescent memory cells generated after the primary response. The population of memory cells is less heterogeneous than the initial naive and activated cell populations. This is expected: the immune response selects cells able to fight the pathogen, and generates a memory cell population characterized by its ability to fight a given pathogen, and not all pathogens.
Then from Figure 6.(d) more APC bearing the same antigen are added to the model to mimic a re-stimulation of the CD8 T cell population. A secondary response starts, memory cells are being activated. On Figures 6.(e) and 6.(f) cells follow once again the linear differentiation scheme and in the end generate memory cells again. During the secondary response, almost all memory cells initially present are activated at the same time since ρ, the CD8 T cell population density, is very close to a Dirac mass (lesser heterogeneity). Actually, due to the concentration on a Dirac mass of the memory cell population, the pool of memory cells may present an absence of heterogeneity which means that all the memory cells will behave the same way during the secondary response. Although a reduction of heterogeneity is expected during the immune response, such a lack of heterogeneity is not really relevant biologically. Figure 7 shows the subpopulation counts obtained during the primary response illustrated in Figures 6.(a) to 6.(c), as well as experimental data from Crauste et al. [16] of CD8 T cell counts measured in mice infected with vaccinia virus. All CD8 T cell counts generated by the model are in the same order of magnitude as experimental counts. The total population count (Figure 7.(f)) matches quite well the data, despite an underestimation of early and late CD8 T cell counts. The peak of the response occurs between days 5-8 post-infection, and reaches around 2.10 5 cells (see Crauste et al. [17]). The contraction phase has the right order of magnitude: it ends between days 10-15 post-infection, and approximately 10% of CD8 T cells at the peak of the response finally account for long lasting memory cells.
Regarding subpopulation counts, naive cells (Figure 7.(a)) rapidly reach nonmeasurable counts, as activated (Figure 7.(b)), early effector (Figure 7.(c)) and late effector (Figure 7.(d)) cells present a burst before going extinct, while memory cells (Figure 7.(e)) are the only CD8 phenotype remaining after the primary response. From Figure 7.(b), curves have been shifted to the right with a three-day delay, since in our model activation of naive cells and generation of activated cells occur instantly after APC contact, while in reality APC do not present themselves to naive cells right after an infection, and a three-day delay is experimentally observed [20,44]. One can notice that dynamics of late effector cells (Figure 7.(d)) are poorly described, compared to other subpopulations, in particular the late effector population goes extinct much faster than what is experimentally observed. This comes from the calibration of molecular dynamics, that do not allow cells to spend enough time in the late effector state. As a consequence, a low number of cells is observed on Figure 7.(f) around day 12 post-infection. Apart from these dynamics, the kinetics of the total population and the subpopulations are in good agreement with experimental data [16].
We then consider the case where the steady state (x 2 , y 2 ) is unstable due to the delay τ , and generates a limit cycle. The convoluted density ρ ε is then displayed on Figure 8 over the course of a primary and a secondary immune response.
On Figure 8, the emergence of a limit cycle in System (1) impacts the behavior of the density ρ of CD8 T cells. CD8 T cells still follow the differentiation scheme from Figure 1, but instead of converging towards a Dirac mass at the end of the primary response (see Figure 6.(c)), ρ becomes located on a closed curve in the memory compartment (see Figures 8.(a) to 8.(c)). Then, the re-introduction of APC   Figure 8. Simulation of ρ ε in the case where (1) has an unstable positive steady state. Legend is similar to Figure 6, except that in the memory phase (Figures (c) and (f )) one clearly observes that memory cells present more heterogeneity thanks to the limit cycle appearing in the (µ 1 , µ 2 )-space. intracellular dynamics then leads to a qualitative CD8 immune response more in agreement with biological observations than when the steady state is stable, by allowing maintenance of a heterogeneous memory cell population at the end of the primary response.
5. Discussion. We introduced a multiscale model of the CD8 T cell immune response to an acute infection, accounting for both intracellular dynamics and CD8 T Table 3. Parameters values used for numerical results displayed in Figures 6 to 8. Between Figures 6 and 8, only values of τ , ∆t, and the time step δt are modified, but the ratio τ /δt is preserved. Between Figures 6 and 7, only the total simulation time T is modified, as Figure 7 only describes CD8 T cell counts for the primary response. N.U. means 'no unit'.

Parameter type Parameter
Value Unit Figure 6    cell behavior at the population scale. Compared to most previous multiscale models of the CD8 T cell immune response, this model is based on a deterministic description of cell and molecular processes using differential equations and a continuous approach. Stochasticity is nevertheless introduced in the description of the activation process: the population of activated naive cells that starts proliferating has heterogeneous intracellular content. The model is hence able to describe complex behaviors as well as primary and secondary immune responses, from a qualitative point of view at this stage.
Using an intracellular network developed in a previous work [47] as well as recently obtained experimental data [16], we proposed a simplified and reduced network describing the dynamics of only two proteins, Ki67 and Bcl2, whose evolution at the single cell level is influenced by APC. These two proteins allowed to sort out CD8 T cells into four subpopulations, activated, early effector, late effector and memory cells [16] that are characterized by specific behaviors and functions (regarding proliferation and survival). By describing cell population dynamics with a discrete cell density, modeled by Dirac masses whose dynamics explicitly depend on the intracellular model, we built a model coupling the phenomena occurring at the population scale and the intracellular scale for CD8 T cells that efficiently describes CD8 T cell dynamics and runs much faster than agent-based models. Contrary to Prokopiou et al. [47], who modeled only the early stages of a primary response, we were able to describe the generation of memory cells (as in Terry et al. [54] and Crauste et al. [17]), but also secondary responses, and our model can be used to investigate the consequences of molecular events on the generation of memory cells.
After showing that our model was well-posed and solutions positive (important for the model to be biologically relevant), the study of the steady states of the intracellular system showed the existence, depending on model parameters, of a stable positive steady state in the maturation space, which can represent long-lived memory cells. In this case, similarly to the work done in Friedman et al. [23,24] for a continuous model, numerical simulations showed that CD8 T cells could all asymptotically concentrate around that intracellular steady state (see Figure 6). The study of the delay system (1) also revealed delay-dependent stability conditions for this steady state: the steady state can become unstable for a delay τ large enough, with the appearance of periodic solutions past the bifurcation point. In this case, the model was able to generate a heterogeneous pool of memory cells (regarding their intracellular content) at the end of the primary response, which is more in agreement with biological knowledge (see Figure 8).
Numerical simulations also showed that the model behaves correctly qualitatively regarding both its ability to account for the linear differentiation scheme we adopted and its description of CD8 T cell subpopulation counts overt time. Obviously, a parameter estimation procedure based on experimental data would surely allow a better description of CD8 T cell counts, yet this has not been performed here because it needs technical developments and would considerably increase the size of this manuscript. Numerical resolution of the model is very fast (less than a minute on a personal laptop), which is a great advantage compared to most multiscale models of the CD8 T cell response that are based on agent-based approaches (this is the case of model in Prokopiou et al. [47], built in our group), especially in order to estimate parameter values by confronting the model to experimental data and later use the model to make predictions.
On the downside, some features of CD8 T cell responses are not described in our model. The stimulation of CD8 T cells by APC, for instance, is partly a spatial process: it depends on the position of CD8 T cells and APC in a spatial domain, since it occurs through contacts (as modeled in Prokopiou et al. [47]). In our model, all cells are uniformly stimulated by APC, not taking into account the stochastic repartition and movement of the different types of cells, yet this is compensated by the generation of a heterogeneous activated cell population at the intracellular scale (the same issue occurs for cytotoxic cells attacking APC, or leading other cells to apoptosis through the Fas-FasL interaction). In the same way, our model neglects parameter variability among cells: all cellular processes (proliferation, death) are characterized by the same parameter values for a given type of cells. This could be taken into account by using mixed effects models, in which parameters can be given a distribution law to take into account individual variability. Regarding these aspects, our model can be interpreted as the mean behavior of a model that would take these processes into account. Also, the model does not take into account other types of cells involved in a CD8 T cell response, such as CD4 helper T cells which secrete cytokines into the environment to help other lymphocytes, or T-reg cells whose role is to prevent autoimmune responses [1]. This is however a feature common to most multiscale models, as it is difficult to reach a high level of description while maintaining tractability (not too long simulation times, capacity to estimate parameter values, etc.) of the models. As briefly mentioned above, the model does not consider the fact that CD8 T cells encounter APC in lymph nodes and then migrate through the blood flow. The way we describe APC/CD8 T cell stimulation and the subsequent CD8 T cell response corresponds to a simplified description of what would be observed, in terms of cell counts, only in one lymphoid organ, the blood. However, such complexifications of the model could be implemented once the model has been quantitatively validated, by using the same approach.
The next step of our study is to estimate parameter values by confronting the model to experimental data, and to perform a sensitivity analysis. The model roughly presents thirty parameters, which have been calibrated for the simulations shown in this paper to qualitatively reproduce a CD8 T cell response. We have access to biological data of the CD8 T cell subpopulation counts in response to vaccinia virus infection in mice (data also used in Crauste et al. [16]). The high number of parameters, even though we can set values for a few ones, requires an appropriate method. We plan to use Kriging methods to build a metamodel of our model, based on a planned exploration of the parameter space. The computation time being very low, we can run a decent number of simulations to build our metamodel. After parameter estimation, the model should provide us with an efficient tool to investigate specific questions, regarding for instance the influence of molecular events on the cell population behavior, or the design of optimal vaccination strategies.
Appendix A. Proof of Theorem 3.1. Let T > 0 be fixed, as well as ∆t ∈ R * + and q ∈ N * such that T = q∆t, and n ∈ N * such that the domain Ω A contains a regular square of area 4σ 2 split into n 2 smaller squares (see Figure 3 and previous section). In addition, let G be the Banach space equipped with the norm In order to prove Theorem 3.1, we successively introduce operators T 0 , T 1 , T 2 and T 3 to define the solutions to System (2)-(4). Let be defined, for t ∈ [0, T ], by with P the unique solution to Let with N the unique solution to , ω i,j,k ), i, j ∈ {1, . . . , n}, k ∈ {1, . . . , q}, whose components are solution of Systems (2) and (3) with the corresponding autonomous functions E, P and N . Finally, let Next, we state and prove the Lipschitz properties of the operators T 0 to T 3 in Lemmas A.1 to A.4, allowing us to obtain the Lipschitz property of the operatorT 4 .
Lemma A.1 The operator T 0 is well defined, continuous and for V,Ṽ ∈ G, Proof. The result is straightforward. E is well defined and the Lipschitz property is obtained directly using rough bounds.
Lemma A.2 The operator T 1 is well defined, continuous, uniformly bounded and Moreover, for E,Ẽ ∈ L ∞ ([0, T ], R + ), Proof. To prove the existence of P , we use a Banach fixed point theorem on the mapping where L is a contraction mapping for T small enough (T < 1/(k P + γ E )) on L ∞ ([0, T ), R), and Equation (14) has a unique solution on [0, T ). Using a standard argument, we construct the solution on R + by iterating the Banach fixed point on [T /2, 3T /2), [T, 2T )... P is well defined and continuous, therefore strictly positive for t small enough. By contradiction, let us assume there exists t 1 such that P (t 1 ) < 0. Then, there exists t 0 < t 1 such that P (t 0 ) = 0 and for all t ∈ (t 0 , t 1 ), P (t) < 0. Then we have which means P (t 1 ) > 0, hence a contradiction. We thus obtain that 0 ≤ P (t) ≤ P (0) for all t ≥ 0 and (16) follows. We now focus on proving (17). Given E, E ∈ L ∞ ([0, T ], R), we have Therefore, Finally, using the Gronwall inequality we obtain directly (17).
The operator T 2 is well defined, continuous, uniformly bounded and Moreover, for P,P ∈ C 0 ([0, T ], R + ), Proof. The proof is similar to the proof of Lemma A.2, except that we can directly use the Cauchy-Lipschitz theorem for existence of N for T small enough.
Lemma A.4 The operator T 3 is well defined, continuous and for U = (E, P, N ), with K(T ) a constant depending only on T (with finite values for T = 0) and model parameters.
Proof. Let (i, j, k) ∈ {1, . . . , n 2 } 2 × {1, . . . , q}. We first handle the definition of and ω i,j,k , then prove the Lipschitz properties. For both matters, the properties are obvious for t ≤ t k , so we consider t > t k . For the existence of µ i,j,k 1 , we apply a Banach fixed point theorem to the mapping L 1 , First, since µ i,j,k 1 is continuous, for t small enough it is positive. Noticing that with C 1 = γ P 1 + γ r1 + k 1 + γ l1 , we apply Gronwall's inequality to obtain, for t small enough, x i,j 0 e −C1t ≤ µ i,j,k 1 (t) ≤ x i,j 0 e C1t . Therefore, µ i,j,k 1 is positive. Then, we notice that d dz and thus we obtain that Hence L 1 is a contraction and µ i,j,k 1 exists and is unique for T small enough. Then, we build µ i,j,k 1 for all times using the same method as the one we used to build P in the proof of Lemma A.2 by iterating the Banach fixed point theorem in intervals of length T ([T /2; 3T /2), [T ; 2T ), . . . ).
Regarding the existence of µ i,j,k 2 , we use the Cauchy-Lipschitz theorem and µ i,j,k 2 is well defined for t > t k small enough. The positivity of µ i,j,k 2 is obtained by noticing that µ i,j,k 2 = 0 is a stationary solution and that µ i,j,k 2 is positive for t small enough. Next, let C 2 = γ P 2 + k 2 + r 2 K 2 . Noticing that we have µ i,j,k 2 (t) ≤ y i,j 0 e C2T and therefore µ 2 is defined and continuously differentiable for all times.
Regarding the definition of ω i,j,k , the Cauchy-Lipschitz theorem gives its existence and uniqueness for t − t k small enough. Since ω i,j,k is continuous for t > t k , it is also positive for t − t k small enough. Next, noticing that dω i,j,k /dt ≤ ||F || ∞ , we have ω i,j,k (t k )e −||F ||∞(t−t k ) ≤ ω i,j,k (t) ≤ ω i,j,k (t k )e ||F ||∞(t−t k ) .
This completes the proof.
Using properties (15) to (19), we obtain We deduce that for T small enough (T < T l ), T 4 is a contraction on the Banach space G . A solution to System (2)-(4) is a fixed point of the operator T 4 . Hence, there exists a unique solution to this system for t ∈ [0, T l ). We can apply the same method for any time interval of length T l . We do so on [T l /2; 3T l /2), [T l ; 2T l ),. . . and the solution exists and is unique for all times.