Model-based assessment of the Role of Uneven Partitioning of Molecular Content on Heterogeneity and Regulation of Differentiation in CD8 T-cell Immune Responses

TL;DR

This study introduces a multiscale mathematical model to investigate how uneven molecular partitioning during CD8 T-cell division influences immune response heterogeneity, cell differentiation, and memory cell formation.

Contribution

The paper presents a novel multiscale model linking molecular partitioning unevenness to immune response dynamics and cell fate decisions in CD8 T-cells.

Findings

Uneven molecular partitioning explains immune response heterogeneity.

Partitioning influences the emergence of distinct cell fates.

Degree of unevenness affects memory cell generation.

Abstract

Activation of naive CD8 T-cells can lead to the generation of multiple effector and memory subsets. Multiple parameters associated with activation conditions are involved in generating this diversity that is associated with heterogeneous molecular contents of activated cells. Although naive cell polarisation upon antigenic stimulation and the resulting asymmetric division are known to be a major source of heterogeneity and cell fate regulation, the consequences of stochastic uneven partitioning of molecular content upon subsequent divisions remain unclear yet. Here we aim at studying the impact of uneven partitioning on molecular-content heterogeneity and then on the immune response dynamics at the cellular level. To do so, we introduce a multiscale mathematical model of the CD8 T-cell immune response in the lymph node. In the model, cells are described as agents evolving and interacting in a 2D environment while a set of differential equations, embedded in each cell, models the regulation of intra and extracellular proteins involved in cell differentiation. Based on the analysis of in silico data at the single cell level, we 1 show that immune response dynamics can be explained by the molecular-content heterogeneity generated by uneven partitioning at cell division. In particular, uneven partitioning acts as a regulator of cell differentiation and induces the emergence of two coexisting sub-populations of cells exhibiting antagonistic fates. We show that the degree of unevenness of molecular partitioning, along all cell divisions, affects the outcome of the immune response and can promote the generation of memory cells.

Tables

Table 1. Table 1: Description of the molecular signalling pathways in Figure 1 and corresponding bibliographic references.

n°	Description	Réf.
1	Activated TCR induces the development of IL2 receptors	[42, 56]
2	Activated TCR induces the synthesis of Tbet	[45]
3	Deactivation of activated IL2 receptors	[42]
4	Activation of IL2 receptors	[42]
5	Activated IL2 receptors induce the development of new IL2 receptors	[41]
6	Activated IL2 receptors inhibit the expression of the IL2 gene (via Blimp1)	[43]
7	Activated IL2 receptors induce the expression the IL2 gene	[41]
8	Tbet enhances the inhibition of the IL2 gene by activated IL2 receptors	[44]
9	Internal IL2 is secreted in extracellular environment	[42]
10	External IL2 binds the non-activated IL2 receptors to activates them	[41]
11	Tbet inhibits the secretion of IL2	[57, 58, 59]
12	TCR activation activates IL2 gene (via Erk)	[41]
13	TCR activation inhibits the activation of Caspases (via Erk, Bim, Bax and Bcl2)	[54]
14	Activated IL2 receptors inhibit the activation of Caspases (via Stat5, BAX et Bcl2)	[55]
15	Tbet induces the expression of FasL	[46]
16	FasL activates Fas through cell contact	[47]
17	Activated Fas induces Caspases activation	[47]
18	Tbet activates Tbx21 and induces the synthesis of Tbet (positive feedback loop)	[49, 50]
19	Eomes induces the expression of IL2 receptors	[14]
20	Activated IL2 receptors induce the expression of Eomes (via Runx3)	[7, 13, 53]
21	Activated TCR inhibits Eomes gene expression (via Akt, mTOR, Tcf1, Foxo1)	[7, 13, 51]
22	Eomes inhibits the activation of Caspases (via Bcl2)	[12]
23	Tbet inhibits the expression of Eomes (via IFN$\gamma$, IL12R)	[60, 61]

Table 2. Table 2 : Main rules applying to APCs and CD8 T-cells in the model. ✔: able, ❍: unable

$\backslash$ Property	Division	Apoptosis	IL2	FAS	FasL
Cell type $\backslash$			secretion	expression	expression
APC	❍	✔	❍	❍	❍
Naive	❍	❍	❍	❍	❍
Pre-activated	❍	✔	✔	❍	❍
Activated	✔	✔	✔	✔	❍
Effector	✔	✔	✔	✔	✔
Memory	❍	✔	✔	✔	✔

Equations

$$\frac{\mathbf{d}}{\mathbf{d}t}[R]=$$

$$\frac{\mathbf{d}}{\mathbf{d}t}[L\bullet R]=$$

$$\frac{\mathbf{d}}{\mathbf{d}t}[Tb]=$$

$$\frac{\mathbf{d}}{\mathbf{d}t}[F{s}^{\ast }]=$$

$$\frac{\mathbf{d}}{\mathbf{d}t}[Cas]=$$

$$\frac{\mathbf{d}}{\mathbf{d}t}[E]=$$

$$\frac{\partial [IL2]}{\partial t}=D{\nabla }^2[IL2]+\left({\lambda }_{R3}\frac{[L\bullet R]}{{\lambda }_{R4}+[L\bullet R]}+{\lambda }_1{f}_{APC}\right)\frac{1}{1+{\lambda }_{T4}[Tb]}-\delta [IL2],$$

$$\displaystyle\Omega=\$$

$$+\underset{\text{contact}}{\underbrace{\underset{\text{neighbours }(\vec{x},\overrightarrow{x^{\ast }})}{\Sigma }{J}_{\tau (\sigma (\vec{x})),\tau (\sigma (\overrightarrow{x^{\ast }}))}(1-{\delta }_{\sigma (\vec{x}),\sigma (\overrightarrow{x^{\ast }})})}},$$

$${\Delta }_{\text{motility}}=v(\sigma (x_s))\left(cos(\theta (\sigma (x_s),t)),sin(\theta (\sigma (x_s),t))\right)\cdot (x_g-x_s),$$

$$D_{cell}=\frac{1}{(\#S)(\#V)}\sum _{\text{simulation}\mathcal{S}}\cdot \sum _{\text{mouse}\mathcal{V}}\cdot \sum _{\text{cellular type}\mathcal{C}}\cdot \sum _{\text{time step}t}|\mathcal{S}(\mathcal{C},t)-\mathcal{V}(\mathcal{C},t)|$$

$$D_{prot}=\frac{1}{(\#S)(\#V)}\sum _{\text{simulation}\mathcal{S}}\cdot \sum _{\text{mouse}\mathcal{V}}\cdot \sum _{\text{protein}\mathcal{P}}\cdot \sum _{\text{time step}t}|\mathcal{S}(\mathcal{P},t)-\mathcal{V}(\mathcal{P},t)|,$$

Full text

**Keywords: ** Multiscale modeling, immune response, asymmetric division, agent-based models, immune memory

1 Introduction

Following acute infection, the activation of naive CD8 T-cells by antigen presenting cells (APCs) triggers the synthesis of proteins controlling cell proliferation and differentiation up to the memory state. While CD8 T-cell population dynamics have been widely described, it is of great interest to better understand the molecular mechanisms driving the CD8 T-cell response. In particular, determining the effects of molecular events on the generation of memory cells is necessary for vaccine design improvement. In vivo and in vitro studies have demonstrated that a single presentation of the antigen to naive CD8 T-cells is sufficient to trigger a complete CD8 T-cell immune response [1, 2, 3, 4, 5]. Then, once initiated, antigen-independent molecular pathways drive a program of CD8 T-cell proliferation and differentiation [6, 7].

The CD8 T-cell immune response occurs through four main phases. First the activation of naive CD8 T-cells in secondary lymphoid organs such as lymph nodes (LN) or spleen by APCs through MHC class I antigenic peptide/T-cell receptor (TCR) binding, surface co-receptor/ligands interactions and soluble cytokines secretion. Once activated, CD8 T-cells proliferate quickly during the expansion phase, which expands the initial population by a factor of $10^{3}$ to $10^{5}$ [6, 8]. Concomitantly, activated cells differentiate into effector cells, able to kill infected cells through cytotoxicity. At the end of the expansion phase, known as the peak of the response, the CD8 T-cell population begins a contraction phase, where most of the responding cells die yet leaving a quiescent population of cells with strong re-activation potential: the memory cells. The memory cell population survives the contraction phase and may remain for years in the organism (memory phase) to ensure faster and stronger host-protection against subsequent infection by the same pathogen.

The responding effector population is composite and two subsets with antagonistic fates have been described [9]: memory precursor effector cells (MPEC) and short-lived effector cells (SLEC), characterised by the expression of two proteins KLRG1 and CD127 (IL-7 receptor). Both MPEC (KLRG1loCD127hi) and SLEC (KLRG1hiCD127lo) express effector features (cytotoxicity, proliferation) but MPEC are capable of differentiation into memory cells while SLEC are destined to die during the contraction phase [9]. Thus CD8 T-cell population dynamics arise from cell phenotypic heterogeneity, itself resulting from molecular-content heterogeneity.

Among the genes, transcription factors and proteins involved in the CD8 T-cell response, some seem to play key roles in the differentiation processes. Transcription factors Tbet and Eomesodermin (Eomes) appear to play critical roles in the acquisition of effector and memory phenotypes. It has been shown that the expression of Tbet induces the development of SLEC and represses the development of MPEC profiles [9, 10, 11]. Eomes is not involved in the SLEC versus MPEC fate choice [12, 13]. However, Eomes is necessary for the development of several properties of memory cells (survival, lymph node homing capacities, responsiveness to second infection [11, 12, 14]). Along the differentiation from effector to memory, the concentration of Tbet in a CD8 T-cell decreases, while the concentration of Eomes increases [11, 15].

Since a unique initial antigenic signal can trigger a complete response, additional mechanisms are necessary to generate the observed molecular-content heterogeneity. Chang et al. [16, 17, 18, 19] showed that TCR binding to MHC-class-I peptide-complex results in polarised segregation of proteins in activated CD8 T-cell: some proteins migrate on the TCR side of the T-cell, other migrate on the opposite side. The subsequent division of the activated CD8 T-cell splits the mother cell perpendicularly to the polarisation axis, such that the daughter cell coming from the TCR side (proximal cell) receives more proteins associated to effector lineage, including Tbet, while the other one (distal cell) receives more proteins associated to memory lineage. Asymmetric division of polarised naive CD8 T-cells appears to be one of the major mechanisms regulating CD8 T-cell fate decision.

Nevertheless, the role of asymmetric division of polarised naive cells in the T-cell differentiation process appears to be controversial [20]. While there are several evidences for asymmetric division of polarised naive CD8 T-cells [21], it remains uncertain how this polarisation quantitatively depends on the affinity of the TCR for the MHC-class-I peptide-complex, the duration of the binding, external chemokines and interactions with homotypic CD8 T-cells [21]. Since the asymmetric partitioning of Tbet has been evidenced in mice CD8 T-cells, it will be considered hereafter.

Less is known about the partitioning of molecular content in the course of subsequent cell divisions. However, several studies support the hypothesis that when a cell divides, a random, uneven partitioning of the molecular content occurs [22, 23, 24, 25, 26, 27, 28, 29]. Partitioning of CFSE dye, a cell staining dye used to track cell proliferation through dye dilution, during lymphocyte proliferation has been mathematically studied by Luzyanina et al. [23, 26]. Based on comparison with in vitro experimental data, these studies suggest that uneven partitioning, which does not result from cell polarisation, occurs at T-cell division.

We emphasize that the asymmetric first division of naive cells, which goes through an active polarisation of the cell, has to be distinguished from the random partitioning of the molecular content during the subsequent divisions of non-polarised cells, hereafter referred to as uneven partitioning [29].

In a recent work (Girel and Crauste [30]), we studied how stochastic uneven molecular partitioning, repeated at each cell division, could regulate the effector versus memory cell-fate decision in a CD8 T-cell lineage. To do so, we analysed an impulsive differential equation describing the concentration of the protein Tbet in a CD8 T-cell subject to divisions, where impulses were associated with uneven partitioning of Tbet. In this work, high and low Tbet concentrations were associated with effector and memory phenotypes, respectively. We concluded that, for a low degree of unevenness of molecular partitioning, a CD8 T-cell expressing a moderate concentration of Tbet can still generate both memory and effector cells. If the concentration of Tbet in this cell is high or low enough, the phenotype of the cell and its progeny becomes irreversible, with low Tbet-expresser and high Tbet-expresser differentiating in memory or effector cells, respectively. Moreover, our study indicates that the increase in cell cycle length throughout the immune response [31, 32] favours irreversible cell differentiation.

Several works (see [33] and the references therein), focused on modeling molecular mechanisms of the immune response coupled to cell population dynamics. Most of these works involve agent-based models.

Gong et al. [34, 35] developed a two-compartment model to study how the number of dentritic cells and the level of MHC-peptides on their membrane influence the size and composition of T-cell populations. Since they did not model any dynamics at the molecular level, they were limited in studying the molecular origins of cell differentiation and heterogeneity.

Prokopiou et al. [36] and Gao et al. [37] designed a multi-scale agent-based model of the early CD8 T-cell immune response (Day 3 to 5.5 post-infection). At the population scale, a discrete population of CD8 T-cells and APCs in a LN is modeled by a cellular Potts model (CPM) [38]. At the molecular scale, the dynamics of a simplified molecular regulatory network (MRN) containing some key molecular factors is modeled by a system of differential equations, embedded in each cell of the population, whose state determines cell phenotype and fate. Cells communicate with each other through cell-cell contact and secretion of the cytokine IL2 such that the environment of a cell affects its molecular profile. Parameter calibration resulted in good agreement with in vivo data of an immune response in murine LN after influenza infection, at both cellular and molecular levels.

The model presented in this article has been developed from the multi-scale agent-based model previously introduced in [36, 37]. Since the authors in [36, 37] focused on early events following CD8 T-cell activation, they did not consider processes leading to the generation of memory cells. We enriched their model in order to study a complete response, from the activation of naive cells to the generation of memory cells. In particular, Eomes has been added to the MRN.

In this paper, we are interested in understanding how, from the activation of naive CD8 T-cells, an antigen-independent regulation of intra-cellular molecular content can drive a complete CD8 T-cell response. We particularly focus on the role of molecular-content heterogeneity among a CD8 T-cell population in the generation of memory cells. We first verify our model’s ability to reproduce in vivo data at both cellular and molecular scales. Then we study, in an in silico CD8 T-cell population, the impact of molecular-content heterogeneity on the emergence of sub-populations, characterised by their expression of proteins Tbet and Eomes. We discuss how uneven distribution of molecular content at cell division affects the cellular dynamics (population size, cell differentiation and death) and suggest that memory cell generation efficiency is maximal for a moderate degree of unevenness. Finally, we show that memory cells generated by our model are able to reproduce some features of a secondary CD8 T-cell immune response. Indeed, when restimulated by antigen in silico they generate more cells at the peak of the response and in the memory phase.

2 Material, Methods and Model

2.1 Data

$4\times 10^{5}$ naive CD8 T-cells from CD45.1+ F5 TCR transgenic mice (B6.SJL-PtprcaPepcb/BoyCrl-Tg(CD2-TcraF5, CD2-TcrbF5)1Kio/Jmar) recognizing the NP68 epitope were transferred intravenously in congenic CD45.2+ C57BL/6 mice (C57BL6/J). The day after recipient mice were inoculated intranasally with $2\times 10^{5}$ PFU (plaque forming units) of a vaccinia virus expressing the NP68 epitope [39]. From day 4 to day 22 post-infection, the spleens of infected animals where harvested and the number of F5 transgenic CD8 responder T-cells was assessed by flow cytometry, based on CD8/CD45.1/CD45.2 expression, to distinguish F5 TCR-transgenic responder (CD45.1+CD45.2-) from host (CD45.1-CD45.2+) CD8 T-cells. Naive (CD44- Mki67- Bcl2+), effector (CD44+ Bcl2-) and memory (CD44+ Mki67- Bcl2+) CD8 T-cells have been identified [40]. All experimental procedures were approved by an animal experimentation ethics committee (CECCAPP; Lyon, France), and accreditations have been obtained from the French government.

OT1 CD8 T cells mRNA expression data time courses come from the ImmGen project (http://www.immgen.org). According to the information provided on ImmGen.org, the in vivo mRNA data (Figure 4) were generated for OT1 T-cells stimulated in similar experimental settings i.e. the response of transferred OT1 TCR-transgenic CD8 T-cells following infection by vesicular stomatitis virus expressing their cognate antigen.

2.2 Molecular regulation and IL2 diffusion

We aim at describing the molecular regulation within each CD8 T-cell during a response to an acute infection, and how the dynamical molecular state of a cell characterises its differentiation stage. We present on Figure 1 the MRN that will be used throughout this manuscript and give a detailed description in Table 1. It contains several key molecular factors involved in CD8 T cell proliferation, differentiation, apoptosis, and cell communication. This is an updated version of the MRN developed in [36, 37] that was limited to the description of differentiation up to the effector stage. To account for differentiation into memory cells, we introduced the protein Eomes and its interactions with the rest of the network as documented in the literature. Indeed, Eomes is involved in the development of essential properties of memory cells such as survival, lymph node homing capacities or faster response to antigenic stimulation [11, 12, 14].

2.2.1 Molecular regulatory network

This MRN is initiated upon antigen presentation to a naive CD8 T-cell, through the engagement of the TCR. Antigenic stimulation triggers the synthesis of interleukine-2 (IL2) by the CD8 T-cell and the production of IL2 receptors (IL2R) on the cell membrane [41]. The synthesised IL2 is then released in the environment and can bind its receptor [42] to form IL2-IL2R complex, hereafter referred to as activated IL2R. Activated IL2 receptors enhance the expression of IL2 receptors [41] as well as IL2 synthesis [41]. In the meantime, activated IL2 receptors, jointly with protein Tbet (see below), inhibit the activation of the IL2 gene through the action of the mediator protein Blimp1 [43, 44].

Antigenic stimulation independently stimulates Tbet synthesis [45], a protein involved in the acquisition of cell cytotoxicity. Indeed, Tbet is known to induce the expression of Fas ligand (FasL) [46], a transmembrane protein that can bind to the transmembrane protein Fas to induce cell apoptosis via the activation of Caspases in the Fas-expressing cell [47]. Caspases are a family of proteins playing essential role in cell apoptosis [48]. There exist several types of Caspases involved in CD8 T-cell apoptosis yet, for the sake of simplicity, we aggregated them in a unique variable $[Cas]$. Moreover, Tbet induces its own synthesis (via the gene Tbx21) [49, 50].

Eomes expression, involved in the acquisition of memory phenotype [12], is first inhibited during the activation phase due to engagement of the TCR (via activation of the Akt/mTOR pathway and inhibition of FOXO1 and TCF7) [7, 13, 51]. Eomes is induced later [11, 52] and its expression is enhanced by the activation of IL2 receptors [7, 13, 53]. Eomes promotes the development of new IL2 receptors on cell membrane [14].

The activation of IL2 receptors, of the TCR and the protein Eomes prevents apoptosis by inhibiting the activation of Caspases, in particular through the mediator protein Bcl2 [12, 54, 55]

2.2.2 Intracellular molecular dynamics

Based on the above-described reactions, and from the equations used in [36, 37], we describe the dynamics of the concentrations of non-activated IL2 receptors ($[R]$), activated IL2 receptors ($[L\bullet R]$), Tbet ($[Tb]$), activated Fas ($[Fs^{*}]$), Caspases ($[Cas]$) and Eomes ($[E]$) in a CD8 T-cell with the following system of equations

[TABLE]

All parameters are positive. Parameters $\lambda$ are associated to induction and inhibition effects, $\mu$ are associated to activation and deactivation of transmembrane proteins and $k$ are degradation and dilution rates. The concentrations of System (1)-(6) are assumed to be null in naive CD8 T-cells, and remain null until TCR engagement.

The effects of the external environment on the intracellular system (1)-(6) are taken into account through five variables. The variable $f_{APC}$ (equations (1), (3), (5) and (6)) is equal to the number of APCs bound to the considered CD8 T-cell and accounts for TCR engagement. The variable $G$ (equations (5) and (6)) is equal to [math] in naïve CD8 T-cells and to $1$ otherwise, i.e. in cells that have already met with an APC. It accounts for the fact that up-regulation of Caspases and Eomes described by parameters $\lambda_{c1}$ and $\lambda_{E4}$ is not active in naive cells. The variable $H$ (equation (4)) accounts for the expression of FasL by effector and memory T-cells and for the activation of Fas through cell contact. Hence, $H$ is equal to 1 in a non-naive considered CD8 T-cell in contact with an effector or a memory CD8 T-cell, and equal to 0 otherwise. The variable $[IL2^{cm}]$ is equal to the concentration of IL2 at the cell membrane, in the extracellular environment. Finally, $[Tb^{cm}]$ is defined as the sum of Tbet concentrations in effector and memory CD8 T-cells in contact with the considered CD8 T-cell and acts as a proxy for the expression of Fas in those cells.

We introduced the variable $[E]$ and the associated equation (6) to the system used in [37] in order to account for the synthesis of protein Eomes and its interactions with other molecular factors. The term $\lambda_{E1}[E]$ in (1) accounts for the up-regulation of IL2 receptors by Eomes. Eomes also limits cell apoptosis by activating Bcl-2 gene, as do IL2 and activated TCR. This communal target explains the multiplicative form of the inhibition of Caspases by Eomes, IL2 and TCR in equation (5). We also introduced the function $G$ in (5) to update the dynamics of Caspases concentration from [36, 37].

The positive feedback loop on Tbet is modeled with an order $n$ Hill function in order to allow bistable behaviour of Tbet. As discussed in the introduction, the concentration of protein Tbet can be associated to the level of differentiation of an effector CD8 T-cell, with high level of Tbet correlating with fully differentiated effector cell, while low Tbet levels are associated to memory precursor effector cells. Proposition 1 below, reproduced from [30], gives necessary and sufficient conditions to allow bistable behaviour of Tbet concentration.

Proposition 1** ([30]).**

Assume $f_{APC}=0$, $n>1$ and $\lambda_{T2}(n-1)^{\frac{n-1}{n}}>nk_{T}\lambda_{T3}$, then equation (3) has exactly three non-negative steady states: $0<[Tb]_{u}<[Tb]_{s}$, such that [math] and $[Tb]_{s}$ are locally asymptotically stable and $[Tb]_{u}$ is unstable.

In the following, we will assume that the conditions $n>1$ and $\lambda_{T2}(n-1)^{\frac{n-1}{n}}>nk_{T}\lambda_{T3}$ are fulfilled (see Section 3.2).

System (1)-(6) is embedded in every CD8 T-cell. Nevertheless, cell-cell contacts, stochastic events (cell cycle length, protein distribution at division) and external concentrations of IL2 affect the evolution of the system such that each CD8 T-cell develops a unique molecular profile based on its own history.

2.2.3 Extracellular IL2 diffusion

The secretion of IL2 by CD8 T-cells and its isotropic diffusion in the extracellular domain (with periodic boundary conditions) are modeled by the following PDE, introduced by Prokopiou et al. [36],

[TABLE]

where $[IL2]$ is the IL2 concentration. CD8 T-cells react to extracellular IL2 through their IL2 receptors by means of the $[IL2^{cm}]$ term, in (1)-(2), defined as the sum of $[IL2]$ at the considered cell membrane.

2.3 Cell differentiation and division

Rules controlling cell division (including protein distribution at the division), apoptosis and differentiation are summarised in Table 2 and detailed hereafter. It must be noted that cells properties result from their molecular profile. For example, the properties observed in vivo in memory cells (survival, low IL2 secretion, low cytotoxicity) are not imposed by model rules but acquired as a consequence of their molecular profile. One exception is cell cycle duration (see 2.3.2).

2.3.1 Differentiation

We designed a set of rules based on the linear, irreversible differentiation scheme from [36, 37], allowing the description of a full CD8 T cell response, from the activation of naive cells up to the generation of memory cells. The differentiation pathway is illustrated in Figure 2.

A naive CD8 T-cell binding an APC becomes pre-activated and maintains the contact with the APC thanks to good adhesion properties (cf. Section 2.4 and Table 3 from Appendix Section A). If the concentration $[L\bullet R]$ of activated IL2 receptors in a pre-activated CD8 T-cell reaches a given threshold $IL2R_{th}$, the pre-activated CD8 T-cell becomes activated, leaves the APC, and starts to proliferate. When an activated CD8 T-cell divides, it gives birth to two CD8 T-cells whose states are determined by their respective concentrations of protein Tbet by comparison with a given threshold $Tbet_{th}$: activated if $[Tb]<Tbet_{th}$, effector otherwise. Finally, if the concentration of protein Eomes is greater than the threshold $Eomes_{th}$, a dividing activated or effector CD8 T-cell will differentiate into memory cell and stop proliferating.

2.3.2 Cell cycle length

Division is considered only for activated and effector CD8 T-cells. The cell cycle length (hours) of a cell preparing its $k$-$th$ division ($k\geq 0$) is chosen, at cell birth, from uniform law $\mathcal{U}_{[c_{k}-4,c_{k}+4]}$ where $c_{k}=6+28k^{2}/(k^{2}+100)$ such that the mean duration of the cycle length increases with the number of divisions and can range from 2 to 32 hours [31, 32]. At the outcome of a division, activated and effector CD8 T-cells immediately enter a new cycle.

2.3.3 Protein distribution between daughter cells

When a CD8 T-cell divides, the molecular content of the mother cell is randomly divided between the two daughter cells. To account for protein distribution between daughter cells at each division and for each protein, let us introduce the parameter $m$, defined as the degree of unevenness. We say that divisions are $m\%$ uneven if at division one daughter cell inherits up to $(50+m/2)\%$ of the mother cell’s content, while the second daughter cell receives the rest, that is at least $(50-m/2)\%$ of the mother cell’s content. Then, the molecular content of each daughter cell evolves according to System (1)-(6) until the next division.

For the sake of clarity, we emphasise that the degree of unevenness $m$ is not the percentage of proteins received by daughters cells at each division but indicates to what extent stochastic molecular partitioning can be uneven. Based on estimation from [26], we consider that divisions are $10\%$ uneven, so that the most uneven partitioning in this case would split $45\%$ and $55\%$ of the mother cell’s proteins in the two daughter cells respectively.

The exact value of each daughter cell molecular content at birth is randomly chosen according to a probabilistic law, as detailed hereafter. Each protein concentration $[i]$ of the six proteins in System (1)-(6) is unevenly distributed among daughter cells: one cell inherits $k_{i}[i]$ and the other $(2-k_{i})[i]$. Coefficients $k_{i},\ i=1,\dots,6$, are different for each protein, each cell, and each division, and are chosen from the probability law $\mathcal{U}_{[1-m/100,1]}$. Unless otherwise indicated, we consider $10\%$ uneven divisions [26], i.e. $k_{i}\in[0.9,1]$ for $i=1,\dots,6$. One may note that $k_{i}\in[0,1]$ so the quantity of molecular material is preserved at each division, given that the volume of each daughter cell is half the volume of the mother. Different degrees of unevenness will be considered in Section 3.3.

One special case of division is the asymmetric division, and its associated unequal repartition of Tbet between daughter cells. To account for polarisation of naive cells by antigenic signalling and the consecutive asymmetric divisions, the first division of a CD8 T-cell following its activation by an APC is characterised by a very specific uneven distribution of protein Tbet only between the two daughter cells: a coefficient $K$ is randomly chosen from the uniform law $\mathcal{U}_{[0.5,1]}$, one of the daughter cells is arbitrarily designated as the proximal daughter and receives a concentration $(2-K)[Tb]$ for protein Tbet while the other one is designated as the distal one and receives a concentration $K[Tb]$ where $[Tb]$ is the Tbet concentration in the mother cell, so that Tbet accumulates in proximal cells [17, 18]. Other proteins concentrations are partitioned according to the previously mentioned rule, see paragraphs above.

2.3.4 Apoptosis

CD8 T-cell apoptosis occurs as soon as Caspases concentration $[Cas]$ reaches the threshold $Caspases_{th}$. APCs are present from the beginning of the simulation and their lifetime is randomly chosen from the uniform law $\mathcal{U}_{[48,96]}$ (hours). APCs’ only role is to activate naive CD8 T-cells, so we do not model any molecular activity within APCs. Dead cells are removed from the domain.

2.4 Spatial modeling and cellular interactions

At the cell population scale, we use a cellular Potts model (CPM), also known as Glazier-Graner-Hegeweg model [38], to describe a population of CD8 T-cells and APCs evolving in a two-dimensional domain. Basically, a CPM is a time-discrete algorithm where cells, or agents, are defined as sets of nodes and move on a lattice, one node at a time, according to probabilistic rules based on the minimisation of the energy of the system, known as the Hamiltonian.

In our model, based on that from [36, 37], the domain is a square lattice of $S=150\times 150$ nodes with periodic boundary conditions. Each node $\vec{x}$ bears an index $\sigma(\vec{x})$. A set of nodes bearing the same index $\sigma$ defines a cell, also denoted by $\sigma$. Finally, each cell $\sigma$ has a type $\tau(\sigma)$ defining its properties. In our case, the different types are: extracellular medium, APC, naive, pre-activated, activated, effector and memory CD8 T-cell. Note that, technically, the extracellular medium is considered as a cell, denoted by $\sigma_{e}$.

Cell (including extracellular medium) size variation and displacement result from the succession of copies of index from nodes to neighbour nodes, based on the minimisation of the Hamiltonian $\Omega$ (see equation (8)), thanks to a simulated annealing algorithm. More precisely, at each iteration, known as Monte Carlo Step (MCS), of the CPM, the following algorithm is executed $N=3\times S$ times:

**Step 1 **

Randomly choose a source node $x_{s}$ and, among its first order neighbours, a target node $x_{g}$.

**Step 2 **

Compute the Hamiltonian $\Omega$, and the putative Hamiltonian $\Omega^{\prime}$ that would be obtained if node $x_{s}$ would copy its index on node $x_{g}$, i.e. if cell $\sigma(x_{s})$ incorporates the node $x_{g}$.

**Step 3 **

Compute $\Delta\Omega=\Omega-\Omega^{\prime}+\Delta_{\text{motility}}$ (see equation (9) below) to evaluate the energy cost of such a copy. If $\Delta\Omega>0$, $x_{s}$ copies its index $\sigma(x_{s})$ on $x_{g}$, i.e. $x_{g}$ is integrated by cell $\sigma(x_{s})$. Else, the copy is accepted with probability $\exp(-\Delta\Omega/T)$, known as Boltzman probability, where parameter $T$ characterises the propensity of the system to evolve.

Note that it is conventional to consider $N=S$ pixel copy attempts per MCS. However, in that case the maximum speed cells can reach is limited to approximatively 0.1 pixel per MCS [62], which eventually defines a finer time resolution than expected for the integration of differential equations. We emphasise that this limitation can be removed by increasing this number (here $N=3\times S$).

The Hamiltonian $\Omega$ is computed using the following formula:

[TABLE]

where $J_{\tau_{1},\tau_{2}}$ accounts for the contact energy between two cells of types $\tau_{1}$ and $\tau_{2}$. Thanks to the term $1-\delta_{\sigma(\vec{x}),\sigma(\vec{x^{*}})}$, two neighbour nodes belonging to the same cell do not generate contact energy. $p_{\sigma}$ and $a_{\sigma}$ are the actual perimeter and area of cell $\sigma$, respectively, whereas $P_{\tau(\sigma)}$ and $A_{\tau(\sigma)}$ are the target perimeter and area, respectively, for a cell of type $\tau(\sigma)$ ; perimeter and area constraints then penalize the configurations where the effective perimeter and area are distant from the target ones. Parameters $\lambda_{area}$ and $\lambda_{pm}$ define the weights of those two constraints. The perimeter constraint has been added to the definition used in [36, 37] in order to avoid potential cell fragmentation.

The energy $\Delta_{\text{motility}}$ is defined by

[TABLE]

where $v(\sigma(x_{s}))$ is the weight associated to the motility energy for the cell $\sigma(x_{s})$ and $\theta(\sigma(x_{s}),t)$ is the privileged angle of direction for the cell $\sigma(x_{s})$ at time $t$, randomly updated along the simulation. The operator “$\cdot$” stands for the dot product. Thus, $\Delta_{\text{motility}}$ is all the more high (and then the copy is all the more probably accepted) that the copy direction $(x_{g}-x_{s})$ aligns with $(cos(\theta(t)),sin(\theta(t)))$.

2.5 Numerical resolution

The initial cell population is composed of 30 naive CD8 T-cells and 3 APCs. A simulation requires 30,000 iterations (MCS) corresponding to 20 days and 20 hours in the real time, that is, 1 MCS represents 1 minute. When a simulation starts, APCs are already present in the LN, ready to activate naive CD8 T-cells. We consider the initial time to be day 4 post-infection (D4 p.i.) since our in vivo data set starts D4 p.i..

We assume that a node of the lattice corresponds to $4\times 4\mu m^{2}$ for biological interpretation. The target cell area is chosen to be 9 nodes ($144\mu$m2) for CD8 T-cells and 140 nodes ($2240\mu$m2) for APCs. The target perimeter for CD8 T-cells is $48\mu m$ in order to favour compact shapes ; there is no constraint on APC perimeter. The simulations have been performed using CC-IN2P3 servers on Compucell3D software [62] with, unless otherwise stated, the parameter values from Tables 3, 4, 5 and 6 from Appendix Section A.

In Section 3.4, we study the ability of our model to simulate a secondary response, also called memory response. Our model has first been calibrated in order to reproduce an in vivo primary response against Listeria monocytogenes (Lm) infection from [63] (see Figure 8). Then, the same parameter values have been used to simulated both a primary and secondary responses. However, secondary response simulations are performed with initially 3 APCs and 30 memory CD8 T-cells (instead of 30 naive CD8 T-cells for the primary response) that are able to bind an APC to become pre-activated, then the differentiation scheme presented in Section 2.3.1 applies. The molecular profile of the initial memory cells is set as the asymptotic molecular profile developed by memory cells at the end of a primary response, as discussed in Section 3.2.

2.6 Model calibration

Parameters of equations (1)-(9) have been calibrated on in vivo data using parameter values from [36, 37]. Since handling big cell populations with an agent-based model implies expensive computation time, we focused on fitting the proportion, rather than the number, of CD8 T-cells in each state of differentiation among the whole cell population. In order to compare in silico and in vivo data at both cellular and molecular scales we minimised the metric $D=D_{cell}+D_{prot}$ where

[TABLE]

and

[TABLE]

with $\#S$ the number of simulations performed with a given set of parameters and $\#V$ the number of mice from which in vivo data have been collected. $\mathcal{S}(\mathcal{C},t)$ (resp. $\mathcal{V}(\mathcal{C},t)$) is the ratio between the number of cells of type $\mathcal{C}$ and the size of the CD8 T-cell population at time $t$ in the simulation $\mathcal{S}$ (resp. the mouse $\mathcal{V}$). $\mathcal{S}(\mathcal{P},t)$ (resp. $\mathcal{V}(\mathcal{P},t)$) is the ratio between the mean concentration of protein (resp. expression of mRNA) $\mathcal{P}$ among the CD8 T-cell population at time $t$ in the simulation $\mathcal{S}$ (resp. the mouse $\mathcal{V}$) and the maximal concentration (resp. expression) observed among all the time steps.

Since pre-activated and activated cellular types are not identified in in vivo data, we gathered pre-activated with naive T-cells and activated with effector T-cells. Then cellular types $\mathcal{C}$ in equation (10) are: naive/pre-activated, activated/effector and memory. In equation (11), quantities $\mathcal{P}$ are the ones for which we have relevant in vivo mRNA expression data at our disposal: IL2 receptors, Tbet and Eomes.

Note that we did not perform a parameter estimation procedure, but a calibration of our model based on experimental data. Evaluation of accuracy and sensitivity of parameter values have been investigated in previous studies [37, 36]. Since we modified the model to account for differentiation in memory cells, a sensitivity analysis of our model to parameter $Eomes_{th}$ is presented in Appendix Section B.

3 Results

3.1 modeling the CD8 T-cell immune response at both cellular and molecular scales

We first briefly illustrate our model’s ability to reproduce in vivo dynamics at both cellular and molecular scales. The evolution of the composition of a CD8 T-cell population from D4 to D22 p.i. is presented on Figure 3.A. In both in vivo and in silico data, naive CD8 T-cells are negligible after D6 p.i.. At the peak of the response, occurring D8 p.i. both in vivo and in silico, more than $94\%$ of the CD8 T-cells are in the activated or effector state, while the memory population emerges during the subsequent contraction phase. As a result of effector cell death and differentiation, memory cells represent the major part of the population on D22 p.i.. Figure 3.B shows the size, in number of cells, of the CD8 T-cell population. The qualitative in vivo dynamics is quite well reproduced: antigen presentation to naive CD8 T-cells triggers clonal expansion, population size reaches a peak D8 p.i. followed by a contraction phase where most cells ($64\%$ and $67\%$ in vivo and in silico respectively) die.

On Figure 4, in silico predictions are compared to the mean IL2 receptors, Tbet and Eomes mRNA expression levels of CD8 T cells activated in vivo. The kinetics of IL2R and Tbet are well reproduced. Indeed, as a result of TCR engagement, IL2R concentration sharply increases and reaches a peak D5 p.i., allowing cells to capture IL2 and get activated. Then IL2R concentration decreases until D8 p.i. and slowly increases from D8 to D15 p.i.. Tbet concentration increases from D4 to D6 p.i. and remains stable until D8 p.i., then decreases until D15 p.i.. Mean Tbet concentration consistently correlates with the size of effector CD8 T-cell population (Figures 3.A and 4.B) and is in agreement with its role in the control of cytotoxicity and cell apoptosis. Regarding Eomes concentration, the in vivo increase between D4 and D8 p.i. is well reproduced by our model, however the increase observed between D8 and D15 p.i. does not match the in vivo data. As cells evolve towards a memory phenotype, in silico Eomes concentration increases and up-regulates the expression of IL2R (Figure 1) to exacerbate the sensitivity of memory cells to IL2. It should be noted that various works support that Eomes expression increases in effector cells progressing toward a memory phenotype [11, 15, 45], contrary to what is observed in the mRNA dataset from Immgen.

3.2 Cellular dynamics arise from cellular heterogeneity

In our model, each cell develops its own molecular profile, resulting in a heterogeneous cell population. Consequently, studying the mean concentration of a given protein among the population, as shown on Figure 4 for example, is not sufficient to understand the molecular dynamics among the CD8 T-cell population.

To study the molecular-content heterogeneity and its role in cellular dynamics, we show in Figure 5 the in silico concentrations of Tbet, Eomes and Caspases in each CD8 T-cell of the population at different times of the response. Cells were ranked according to their Tbet content. D5 to D8 p.i., corresponding to the clonal expansion phase (see Figure 3), concentrations are heterogeneous but uniformly distributed around the mean value. Most of that heterogeneity comes from the conditions of activation and from molecular partitioning at cell division. Yet from D8 to D24 p.i., corresponding to the contraction phase, two sub-populations of cells clearly emerge: one with high concentration of Tbet (centred around $[Tb]_{s}\approx 118$ mol/L) and one with low concentration of Tbet ($\approx 0$ mol/L). The unstable steady state of (3), defined in Proposition 1 and separating the stable equilibria [math] and $[Tb]_{s}$, is given by $[Tb]_{u}\approx 21$ mol/L. Moreover, cells expressing high levels of Tbet express high levels of Caspases and low levels of Eomes, a molecular profile associated with cell death and poor memory potential. On the contrary, cells expressing low levels of Tbet have good survival and memory differentiation properties since they express low levels of Caspases and high levels of Eomes. Progressively, cells with high concentrations of Tbet die (when their concentrations of Caspases reach the threshold $Caspases_{th}\approx 19$ mol/L) and cells with low concentrations of Tbet differentiate into memory cells and stop proliferating (when the concentration of Eomes reaches $Eomes_{th}=16$ mol/L). On D24 p.i. there is no cell with intermediary profile, most of the cells have differentiated into memory cells while a few effector cells with high Tbet concentrations still survive. One can observe that the molecular profiles of memory cells converge to the same state where $[Tb]=0$ mol/L, $[E]\approx 26$ mol/L and $[Cas]\approx 9$ mol/L.

The coexistence of two sub-populations characterised by their concentrations of Tbet explains the population dynamics observed on Figure 3. That is, the contraction of the cytotoxic effector cell population simultaneously with the emergence of a memory cell population with survival properties.

As discussed in the introduction, responding CD8 T-cells can be distinguished between short-lived (SLEC) and memory precursor (MPEC) effector cells based on the expression of two proteins: KLRG1 and CD127 [7, 15]. In this section, we investigated how, in our model, the heterogeneity of Tbet concentrations among a CD8 T-cell population explains the emergence of two sub-populations of CD8 T-cells. The first one, expressing high concentrations of Tbet, could be comparable to SLEC that exhibit properties such as apoptosis and cytotoxicity, a process regulated by Tbet. The second one (memory potential, survival) would be similar to the MPEC population. This is consistent with the litterature, since Tbet is known to favour the development of SLEC, to the detriment of MPEC [9, 10, 11].

3.3 Moderate uneven molecular partitioning favours efficient generation of memory cells

A major source of heterogeneity in our model is the uneven molecular partitioning at cell division determined by the degree of unevenness $m$ (see Section 2.3.3). We compare on Figure 6 the sizes of the CD8 T-cell population at the peak of the response as well as the sizes of the memory population on D25 p.i. for different degrees of unevenness, that is the extent of unevenness of the stochastic molecular partitioning. We do not however modify the degree of unevenness of the asymmetric first division, consecutive to the polarisation of the cell due to APC binding [17, 18], see Section 2.3.3.

First, Figure 6 shows that the size of the CD8 T-cell population at the peak of the response decreases as the degree of unevenness increases. Indeed, the more uneven the molecular partitioning, the sooner CD8 T-cells expressing high levels of Caspases or Eomes appear and then the sooner cells die by apoptosis or differentiate in non-proliferating memory cells.

Second, the relation between the degree of unevenness and the size of the memory population generated at the end of the response is not monotonous: the biggest memory populations are observed when considering a moderate unevenness ($10\%$-$50\%$).

In Section 3.2, the role of Tbet concentration in determining the fate (death or memory differentiation) of an effector CD8 T-cell has been discussed. Additionally, we showed in [30] that the progression of a cell lineage towards death or memory differentiation can be slowed down or reversed by molecular partitioning depending on cell cycle length, initial Tbet concentration and the degree of unevenness. This stressed, on a simplified model, the influence of the degree of unevenness on cell fate choice regulation.

On the opposite, when molecular partitioning is symmetrical ($m=0$) and no further T-cell-APC interactions are assumed, there is no more source of stochasticity and consequently all the CD8 T-cells of the same lineage express the same concentration of Tbet. As a consequence of Proposition 1, this concentration irreversibly converges either to $[Tb]_{s}$ (high Tbet concentration) or to [math] mol/L (low Tbet concentration). This irreversibly leads to apoptosis (high Tbet concentration) or memory differentiation (low Tbet concentration) of the whole cell lineage.

Thus, our result clearly stresses that uneven partitioning allows the maintenance of a CD8 T-cell compartment with undetermined fate for some time, through cell fate reversibility. As long as it is maintained, this compartment is able to produce both effector cells destined to die and memory cells.

We also showed in [30] that the higher the degree of unevenness, the more reversible the cellular fate. Surprisingly, strong unevenness ($65\%-80\%$) results in smaller memory cell populations (Figure 6). In fact, strong unevenness favours the fast emergence of daughter cells with very high or low concentrations of Tbet such that those cell lineages are likely to die or to generate memory cells. In particular, effector cells with high memory potential poorly expand before they differentiate hence this leads to the generation of fewer memory cells.

To discuss the efficiency of memory cell generation, we compare on Figure 6 the number of memory cells generated at the end of the response to the number of cells at the peak of the response, viewed as an indicator of the energetic cost of the response for the organism (red crosses). Figure 6 suggests that the degree of unevenness in molecular partitioning impacts memory generation, with the better ratio (more than $30\%$) obtained when considering 50$\%$ uneven molecular partitioning.

3.4 Memory response

One of the characteristics of memory cells is their capacity to mount more rapid effector response than naive cells and to generate an increased fraction of memory cells [64]. To test whether the memory cells generated by our model exhibit some of these features we compared the in silico primary response with a secondary response of in silico generated memory cells.

Figure 7 shows the in silico memory response (or secondary response), obtained with an initial population of 30 memory T-cells, as described in Section 2.5. This secondary response is compared to the primary immune response starting with 30 naive CD8 T-cells (Section 3.1). The in silico secondary response is characterised by a bigger CD8 T-cell population, at any time of the response. From the primary to the secondary response, there is a small increase in the size of the sub-population of activated and effector cells but the major change is in the size of the memory population. Indeed, the number of memory CD8 T-cells increases much faster during the secondary response such that D29 p.i. the memory population is two times bigger than during the primary response. This can be explained by the fact that memory cells are activated faster than naive cells, thanks to their molecular profile. Indeed, memory cells express higher concentrations of IL2 receptors than naive cells, since it is sustained by the expression of Eomes. Consequently, the threshold $IL2R_{th}$ (see Section 2.3.1) is reached sooner when starting with memory cells than with naive cells. As a result, the concentration of Tbet, up-regulated during APC binding, is lower after the activation of a memory cell than after the activation of a naive cell, and low Tbet level is associated to memory precursor fate and low cytotoxicity.

On Figure 8, we compare in silico primary and secondary responses from our model with in vivo primary and secondary responses against Lm infection from [63]. Since our model has been calibrated to fit the primary response data, we do not aspire to reproduce the quantitative dynamics of the secondary response, but rather to study its qualitative properties. Namely, the secondary response is characterised by a slower and weaker contraction phase, from the peak of the response D7 p.i. to the last time point D29 p.i.. This weaker contraction could be explained by an early production of memory cells that leads to a large population of memory cells, as it is the case in our model (Figure 7).

4 Discussion

Activation of naive CD8 T-cells triggers a primary immune response, characterised by a well orchestrated program of cell proliferation, differentiation, death and migration. It is now well known that the responding CD8 T-cell population is heterogeneous and that a single naive T-cell can generate differently fated cells [65]. However, evaluating how cellular and molecular events contribute to that heterogeneity and identifying its consequences on the outcomes of the immune response remain fundamental questions.

With this in mind, we expanded a hybrid multi-scale model of the CD8 T-cell immune response, where cell behaviour is determined by intracellular molecular dynamics. Model parameters have been calibrated using in vivo data at both cellular and molecular scales. Because of expensive running time, we were led to simulate small cell populations so that we focused on semi-quantitative fitting criteria. After calibration, our model succeeded in reproducing the temporal dynamics of the response regarding the size of the CD8 T-cell population and the proportion of cells in each differentiation stage. Apart from a discordance between in silico and in vivo mean concentration of Eomes on day 15 p.i., our model captured the dynamics of the mean concentration of IL2 receptors, Tbet and Eomes, which play key roles in the differentiation processes.

In addition to reproduce primary responses, our model easily produces secondary responses. Memory cells generated during the in silico primary response succeeded in mounting a stronger secondary response upon antigenic stimulation (Figure 7 and Figure 8). It should be noted that the differences between outcomes of the primary and secondary in silico immune responses only depend, in this work, on the difference between the molecular profile of memory and naive CD8 T-cells and do not take into account a lot of characteristics of the secondary response described in the literature such as: biggest initial CD8 T-cell population [8], shorter cell cycle [66] or sensitivity to inflammatory cytokines, such as IL12 [15].

We discussed how a deterministic description of molecular concentration dynamics combined with stochastic events, such as uneven partitioning of molecular content at division, can regulate the emergence and the maintenance of two sub-populations of CD8 T-cells. Those sub-populations, characterised by their molecular profiles, coexist but express different properties and antagonistic fates, comparable to those of SLEC and MPEC described in the literature [9]. From that observation, we showed that the dynamics observed at the cellular scale (cell differentiation, population size) could be explained by molecular-content heterogeneity among the cell population, which mostly originates from uneven partitioning of molecular content. We did not however consider the effect of stochastic fluctuations of gene expression, known to be an important source of heterogeneity [67]. Interestingly, Huh and Paulsson [25] showed that both stochastic gene expression and stochastic partitioning of molecular content are equally good to explain the heterogeneity observed at cell division and suggested that much of the heterogeneity usually attributed to the former actually results from the latter.

In our model, cell phenotypic heterogeneity, associated with molecular-content heterogeneity, first arises upon asymmetric division of polarised naive cells. This heterogeneity is thereafter continuously regulated throughout the whole response by means of uneven partitioning of molecular-content at each division. This is in agreement with the observations of Lemaître et al. [68] who state that T-cell diversification is a continuous process, spread over the whole response, including the asymmetric first division and late events occurring throughout subsequent divisions. Besides, Lemaître et al. [68] pointed out that cellular heterogeneity, that could result from variations in naive T-cell responsiveness to cytokines or TCR signalling, pre-exists prior to the first division. In this article we did not consider preexisting heterogeneity among the naive T-cell pool, that could be achieved by varying the parameter values of System (1)-(6) associated to each naive T-cell. We can expect that it would confer to each naive T-cell a predisposition to engender a cell lineage oriented toward either apoptosis or memory differentiation. Moreover, the initial heterogeneity among naive T-cells could be conserved through the response, then leading to a heterogeneous pool of memory cells, a feature that is not reproduced by our model [69].

Polarisation of naive cells upon antigenic stimulation has been observed in CD4 T-cells [21, 70] and B-cells [21, 71, 72]. This polarisation results in asymmetric division of naive cells and may induce heterogeneous cell fates [70, 71, 72]. Regarding the subsequent divisions, it can be thought that they are subject to uneven and random partitioning of molecular content since this phenomenon has been reported in many types of cells, including yeast, bacteria and T-cells [22, 23, 24, 25, 26, 27]. However, the contribution of uneven and random partitioning of non-polarised cells in the development of heterogeneous cell fate has not been studied yet. To that end, it would be interesting to extend the approach developed in our study to the differentiation of other lymphocytes, such as B-cells or CD4 T-cells.

In our study, increasing the degree of unevenness of molecular partitioning reduces the expansion size of the whole CD8 T-cell population whereas the size of the sub-population of memory cells is maximal for intermediate degrees of unevenness. As a consequence, the ratio between the number of memory cells generated and the magnitude of the response at its peak, viewed as a measure of memory generation efficiency, is maximised when considering a $50\%$ degree of unevenness. As discussed above, molecular partitioning is not the only regulator of heterogeneity. In this regard, we can believe that our evaluation overestimates the value of this optimal degree of unevenness and rather indicates that generating a moderate heterogeneity all along the immune response leads to efficient memory generation.

In our manuscript, when the degree of unevenness is $m=10\%$, each daughter cell inherit from 45 to 55$\%$ of the mother cell’s molecular content, with uniform probability distribution. The unevenness of molecular partitioning remains difficult to measure experimentally. Based on in vitro experimental data of CFSE dye expression, Luzyanina et al. [26] estimated that the two daughter cells inherit of 42,3$\%$ and 57,7$\%$ of the mother cell’s molecular content, respectively. Rather than considering a uniform probability distribution and a degree of unevenness, we could consider that the molecular partitioning is a binomial phenomenon [24], i.e. each protein has the same probability to be attributed to each daughter cell. Such a discrete distribution can be approximate by a continuous and truncated (to avoid negative values) normal distribution whose variance would characterise the level of unevenness.

Note that, in works dealing with the CD8 T-cell immune response, it is usual to consider that 5 to 10$\%$ of the cells present at the peak of the response survive the contraction phase and differentiate into memory cells [8]. This is consistent with our results only for symmetric divisions or for divisions with high ($65\%$ to $80\%$) degrees of unevenness. However, this hypothesis can be challenged, as pointed out in [40], as for the actual in vivo data presented in Figure 3, D22 p.i. the memory population size is $19.5\%$ of the whole population at the peak of the response, D8 p.i.. This suggests that the amplitude, and possibly the kinetics, of the cellular contraction is not only an inherent feature of the CD8 immune response but also depends on external factors such as inflammatory factors.

In many mathematical models of the CD8 T-cell immune response, as those referenced in [6], cell proliferation and differentiation depend on the amount of pathogen, in the manner of prey-predator models used in ecology. In our model a brief initial antigenic stimulation of naive CD8 T-cells is sufficient to trigger an autonomous program of proliferation and differentiation, as stated in the literature [1, 2, 3]. However, while dispensable, in vivo inflammatory signals can affect the immune response outcome [73]. A motivating perspective is to evaluate the respective contributions of both the autonomous program and extrinsic inflammatory factors to the immune response, so that the latter could be tuned by mastering the inflammatory environment. For example, extending our model by incorporating the inflammatory cytokine IL12, secreted by APCs, could markedly affect the effector/memory cell balance since IL12 is known to respectively promote and repress Tbet and Eomes synthesis [9, 57, 74].

Cell cycle length depends in our model on the number of divisions the cell has undergone. It would be instructive to introduce a molecular control of cell proliferation, since the putative existence of coexisting sub-populations with disparate cycle lengths could considerably impact the cellular dynamics. One could for instance consider the transcription factor Foxo1, known to induce Eomes expression while repressing that of Tbet and inhibiting cell cycle progression [75], suggesting that the TbetloEomeshi memory precursor cells discussed in Section 3.2 might adopt a longer cycle than the Tbethi Eomeslo cells.

In conclusion, our agent-based multiscale model successfully reproduced several aspects of the CD8 T-cell immune response at both molecular and cellular scales. Even though we cannot infer quantitative conclusions from this study, it suggests that uneven partitioning of molecular content at cell division, as a source of heterogeneity, can modulate cell fate decision and act as a regulator of the magnitude of the response and of the size of the memory cell pool. Actually, we did not consider intermediaries, namely DNA transcription and mRNA translation, between gene activation and protein synthesis. Consequently, our molecular model is an amalgam between gene activity and protein synthesis. Therefore, while our argumentation is based on uneven partitioning of the molecular content, it could also stand for the situation where, when a cell divides, the two daughter cells inherit different gene activity levels for each gene. All in all, our study focus on molecular heterogeneity generation upon cell division in general, rather than the specific case of molecular partitioning. It stresses that dynamics observed at the cellular scale — including the initiation of the contraction phase and the origin of memory cells — can be explained by structural molecular-content heterogeneity, that is continuously regulated along the response, as CD8 T-cells divide.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Author Contributions

All co-authors discussed the problem, approach and results. SG, OG and FC designed the model. SG ran simulations and performed analysis. CA and JM conducted the experimental studies. SG wrote the paper. All the authors approved the final version.

Funding

This work was performed within the framework of the LABEX MILYON (ANR-10-LABX-0070) of Université de Lyon, within the program “Investissements d’Avenir” (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR) and has been supported by ANR grant PrediVac ANR-12-RPIB-0011.

Acknowledgments

The authors thank the Centre de Calcul de l’Institut National de Physique Nucléaire et de Physique des Particules de Lyon (CC-IN2P3) for providing computing resources, particularly Pascal Calvat and Yonny Cardenas for their valuable help. We also thank the BioSyL Federation and the LabEx Ecofect (ANR-11-LABX-0048) of the University of Lyon for inspiring scientific events. We acknowledge the contribution of SFR Biosciences (UMS3444/CNRS, US8/Inserm, ENS de Lyon, UCBL) facilities. We acknowledge the contributions of the CELPHEDIA Infrastructure (http://www.celphedia.eu/), especially the center AniRA in Lyon.

Data Availability Statement

The datasets generated for this study can be found in the Open Science Framework repository https://osf.io/xbq9r.

mRNA expression data analysed in this study come from the ImmGen project (http://www.immgen.org)

Appendix A Parameter values

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix B Sensitivity analysis to parameter $Eomes_{th}$

In this paper, we enriched previous models [37, 36] in order to allow differentiation into memory cell. In our model, this differentiation occurs when the concentration of protein Eomes in an activated or effector cell crosses the threshold $Eomes_{th}$ (see Section 2.3.1). Here, we discuss how sensitive our model is to this parameter. Sensitivity to other parameters has been investigated in [37] and [36].

It is no surprise that the number of memory cells at the end of the response (D25 p.i) is highly sensitive to $Eomes_{th}$ value, as shown on Figure 9.A. Indeed, for low values of $Eomes_{th}$, the memory precursor cells identified in Section 3.3 differentiate early while for higher values of $Eomes_{th}$ they can accomplish a few additional rounds of division before their Eomes concentrations are sufficiently high to trigger memory differentiation, leading to bigger memory cell populations. Indeed, Figure 9.B illustrates that the lower the value of $Eomes_{th}$, the sooner the pool of memory CD8 T-cells stops to expand.

We can see on Figure 9.C that the number of cells at the peak of the response (i.e. the maximal expansion size reached by the CD8 T-cell population) is less sensitive to the value $Eomes_{th}$. However, it appears that the peak population size slightly increases with the value of $Eomes_{th}$. This is mainly due to the reason exposed above, i.e. memory precursor effector cells proliferate more before their differentiation into non-dividing memory cells.

It is also worth noting that if $Eomes_{th}$ is too low, CD8 T-cells can differentiate into memory cells even if their molecular content does not match with the expected properties of memory cells (survival and low cytotoxicity). Indeed, if $Eomes_{th}$ is $90\%$ (resp. $95\%$) of its base value (see Table 6), an average of 3.6 (resp. 0.4) memory cells die during the response (mean over 5 simulations, results not shown) while for the base value of $Eomes_{th}$ or upper values there is no memory cell death.

Figure 10 represents the mean concentration of IL2 receptors, Tbet and Eomes among the CD8 T-cell population from D4 to D15p.i., for different values of $Eomes_{th}$. Since Eomes expression is associated with memory phenotype, it is no surprise that Eomes concentration increases more slowly, or even decreases, for low values of $Eomes_{th}$. On the opposite, high values of $Eomes_{th}$, associated with big memory cell populations, lead to decreasing Tbet concentrations.