Irene Balelli
email: balelli@math.univ-paris13.fr

Vuk Milišić
email: milisic@math.univ-paris13.fr

Gilles Wainrib
email: gilles.wainrib@ens.fr

Multi-type Galton-Watson processes with affinity-dependent selection applied to antibody affinity maturation

Keywords: Multi-type Galton-Watson process, Germinal center reaction, Affinity-dependent selection, Evolutionary landscapes

We analyze the interactions between division, mutation and selection in a simplified evolutionary model, assuming that the population observed can be classified into fitness levels. The construction of our mathematical framework is motivated by the modeling of antibody affinity maturation of B-cells in Germinal Centers during an immune response. This is a key process in adaptive immunity leading to the production of high affinity antibodies against a presented antigen. Our aim is to understand how the different biological parameters affect the system's functionality. We identify the existence of an optimal value of the selection rate, able to maximize the number of selected B-cells for a given generation.

Introduction

Antibody Affinity Maturation (AAM) takes place in Germinal Centers (GCs), specialized micro-environnements which form in the peripheral lymphoid organs upon infection or immunization [START_REF] Victora | Clonal and cellular dynamics in germinal centers[END_REF][START_REF] De Silva | Dynamics of b cells in germinal centres[END_REF]. GCs are seeded by ten to hundreds distinct B-cells [START_REF] Tas | Visualizing antibody affinity maturation in germinal centers[END_REF], activated after the encounter with an antigen, which initially undergo a phase of intense proliferation [START_REF] De Silva | Dynamics of b cells in germinal centres[END_REF]. Then, AAM is achieved thanks to multiple rounds of division, Somatic Hypermutation (SHM) of the Bcell receptor proteins, and subsequent selection of B-cells with improved ability of antigen-binding [START_REF] Maclennan | B-cell memory and the persistence of antibody responses[END_REF]. B-cells which successfully complete the GC reaction output as memory B-cells or plasma cells [START_REF] Victora | Germinal centers[END_REF][START_REF] De Silva | Dynamics of b cells in germinal centres[END_REF]. Indirect evidence suggests that only B-cells exceeding a certain threshold of antigen-affinity differentiate into plasma cells [START_REF] Phan | High affinity germinal center b cells are actively selected into the plasma cell compartment[END_REF]. The efficiency of GCs is assured by the contribution of other immune molecules, for instance Follicular Dendritic Cells (FDCs) and follicular helper T-cells (Tfh). Nowadays the key dynamics of GCs are well characterized [START_REF] Maclennan | B-cell memory and the persistence of antibody responses[END_REF][START_REF] De Silva | Dynamics of b cells in germinal centres[END_REF][START_REF] Gitlin | Clonal selection in the germinal centre by regulated proliferation and hypermutation[END_REF][START_REF] Tas | Visualizing antibody affinity maturation in germinal centers[END_REF]. Despite this there are still mechanisms which remain unclear, such as the dynamics of clonal competition of B-cells, hence how the selection acts. In recent years a number of mathematical models of the GC reaction has appeared to investigate these questions, such as [START_REF] Meyer-Hermann | An analysis of b cell selection mechanisms in germinal centers[END_REF][START_REF] Wang | A stochastic model of the germinal center integrating local antigen competition, individualistic t-b interactions, and b cell receptor signaling[END_REF],

where agent-based models are developed and analyzed through extensive numerical simulations, or [START_REF] Zhang | Optimality of mutation and selection in germinal centers[END_REF] where the authors establish a coarse-grained model, looking for optimal values of e.g. the selection strength and the initial B-cell fitness maximizing the affinity improvement.

Our aim in this paper is to contribute to the mathematical foundations of adaptive immunity by introducing and studying a simplified evolutionary model inspired by AAM, including division, mutation, affinity-dependent selection and death. We focus on interactions between these mechanisms, identify and analyze the parameters which mostly influence the system functionality, through a rigorous mathematical analysis. This research is motivated by important biotechnological applications. Indeed, the fundamental understanding of the evolutionary mechanisms involved in AAM have been inspiring many methods for the synthetic production of specific antibodies for drugs, vaccines or cancer immunotherapy [START_REF] Ansari | Identification of conformational b-cell epitopes in an antigen from its primary sequence[END_REF][START_REF] Kringelum | Reliable b cell epitope predictions: impacts of method development and improved benchmarking[END_REF][START_REF] Shen | Towards a mathematical foundation of immunology and amino acid chains[END_REF]. This production process involves the selection of high affinity peptides and requires smart methods to generate an appropriate diversity [START_REF] Currin | Synthetic biology for the directed evolution of protein biocatalysts: navigating sequence space intelligently[END_REF]. Beyond biomedical motivations, the study of this learning process has also given rise in recent years to a new class of bio-inspired algorithms [START_REF] Castro | Learning and optimization using the clonal selection principle[END_REF][START_REF] Pang | Clonal selection algorithm for solving permutation optimisation problems: A case study of travelling salesman problem[END_REF][START_REF] Timmis | Theoretical advances in artificial immune systems[END_REF], mainly addressed to solve optimization and learning problems.

We consider a model in which B-cells are classified into N + 1 affinity classes with respect to a presented antigen, N being an integer big enough to opportunely describe the possible fitness levels of a B-cell with respect to a specific antigen [START_REF] Weiser | Affinity maturation of b cells involves not only a few but a whole spectrum of relevant mutations[END_REF][START_REF] Xu | Key mutations stabilize antigenbinding conformation during affinity maturation of a broadly neutralizing influenza antibody lineage[END_REF]. A B-cell is able to increase its fitness thanks to SHMs of its receptors: only about 20% of all mutations are estimated to be affinityaffecting mutations [START_REF] Shannon | Reconciling repertoire shift with affinity maturation: the role of deleterious mutations[END_REF][START_REF] Shlomchik | Clone: a monte-carlo computer simulation of b cell clonal expansion, somatic mutation, and antigen-driven selection[END_REF]. By conveniently defining a transition probability matrix, we can characterize the probability that a B-cell belonging to a given affinity class passes to another one by mutating its receptors thanks to SHMs.

Therefore we define a selection mechanism which acts on B-cells differently depending on their fitness. We mainly focus on a model of positive and negative selection in which B-cells submitted to selection either die or exit the GC as output cells, according to the strength of their affinity with the antigen. Hence, in this case, no recycling mechanism is taken into account. Nevertheless the framework we set is very easy to manipulate: we can define and study other kinds of affinity-dependent selection mechanisms, and eventually include recycling mechanisms, which have been demonstrated to play an important role in AAM [START_REF] Victora | Germinal center dynamics revealed by multiphoton microscopy with a photoactivatable fluorescent reporter[END_REF]. We demonstrate that independently from the transition probability matrix defining the mutational mechanism and the affinity threshold chosen for positive selection, the optimal selection rate maximizing the number of output cells for the t th generation is 1/t, t ∈ N (Proposition 6).

From a mathematical point of view, we study a class of multi-types Galton-Watson (GW) processes (e.g. [START_REF] Harris | The theory of branching processes[END_REF][START_REF] Athreya | Branching processes[END_REF]) in which, by considering dead and selected B-cells as two distinct types, we are able to formalize the evolution of a population submitted to an affinity-dependent selection mechanism. To our knowledge, the problem of affinity-dependent selection in GW processes has not been deeply investigated so far.

In Section 2 we define the main model analyzed in this paper. We give as well some definitions that we will use in next sections. Section 3 contains the main mathematical results. A convenient use of a multi-type GW process allows to study the evolution of both GC and output cells over time. We determine the optimal value of the selection rate which maximizes the expected number of selected B-cells at any given maturation cycle in Section 3.3. We conclude Section 3 with some numerical simulations. In Section 4 we define two possible variants of the model described in previous sections, and provide some mathematical results and numerical simulations as well. This evidences how the mathematical tools used in Section 3 easily apply to define other affinitydependent selection models. Finally, in Section 5 we discuss our modeling assumptions and give possible extensions and limitations of our mathematical model. In order to facilitate the reading of the paper, some technical mathematical demonstrations, as well as some classical results about Galton-Watson theory are reported in the Appendix for interested readers.

Main definitions and modeling assumptions

This section provides the mathematical framework of this article. Let us suppose that given an antigen target cell x, all B-cell traits can be divided in exactly N + 1 distinct affinity classes, named 0 to N . Definition 1 Let x be the antigen target trait. Given a B-cell trait x, we denote by a x (x) the affinity class it belongs to with respect to x, a x (x) ∈ {0, . . . , N }. The maximal affinity corresponds to the first class, 0, and the minimal one to N . Definition 2 Let x be a B-cell trait belonging to the affinity class a x (x) with respect to x. We say that its affinity with x is given by: aff(x, x) = N -a x (x) Of course, this is not the only possible choice of affinity. Typically affinity is represented as a Gaussian function [START_REF] Wang | A stochastic model of the germinal center integrating local antigen competition, individualistic t-b interactions, and b cell receptor signaling[END_REF][START_REF] Meyer-Hermann | An analysis of b cell selection mechanisms in germinal centers[END_REF], having as argument the distance between the B-cell trait and the antigen in the shape space of possible traits.

In our model this distance corresponds to the index of the affinity class the B-cell belongs to (0 being the minimal distance, N the maximal one). Nevertheless the choice of the affinity function does not affect our model.

During the GC reaction B-cells are submitted to random mutations. This implies switches from one affinity class to another with a given probability.

Setting these probability means defining a mutational rule on the state space {0, . . . , N } of affinity classes indices (the formal mathematical definition will be given in Section 3.2).

The main model we study in this paper is represented schematically in Figure 1. It is defined as follows:

Definition 3 The process starts with z 0 ≥ 1 B-cells entering the GC, belonging to some affinity classes in {0, . . . , N }. In case they are all identical, we denote by a 0 the affinity class they belong to, with respect to the antigen target cell x. At each time step, each GC B-cell can eventually undertake three distinct processes: division, mutation and selection. First of all, each GC Bcell can die with a given rate r d ∈ [0, 1]. If not, each B-cell can divide with rate r div ∈ [0, 1]: each daughter cell may have a mutated trait, according to the mutational rule allowed. Hence it eventually belongs to a different affinity class than its mother cell. Clearly, it also happens that a B-cell stays in the GC without dying nor dividing. Finally, with rate r s ∈ [0, 1] each B-cell can be submitted to selection, which is made according to its affinity with x. A threshold a s is fixed: if the B-cell belongs to an affinity class with index greater than a s , the B-cell dies. Otherwise, the B-cell exits the GC pool and reaches the selected pool. Therefore, for any GC B-cell and at any generation, we have:    Probability of cellular apoptosis: P(death)

= r d Probability of cellular division: P(division) = r div Probability of selection challenge: Once the GC reaction is fully established (∼ day 7 after immunization), it is polarized into two compartments, named Dark Zone (DZ) and Light Zone (LZ) respectively. The DZ is characterized by densly packed dividing B-cells, while the LZ is less densely populated and contains FDCs and Tfh cells. The LZ is the preferential zone for selection [START_REF] De Silva | Dynamics of b cells in germinal centres[END_REF]. The transition of B-cells from the DZ to the LZ seems to be determined by a timed cellular program: over a 6 hours period about 50% of DZ B-cells transit to the LZ, where they compete for positive selection signaling [START_REF] Bannard | Germinal center centroblasts transition to a centrocyte phenotype according to a timed program and depend on the dark zone for effective selection[END_REF][START_REF] Victora | Snapshot: the germinal center reaction[END_REF].

P(selection) = r s Ø x x x 1 x 2 Ø x 1 x 1 t t+1 r d 1-r d r div 1-r div r s 1-r s GC GC Sel aff(x 1 ,x )≥ aff s aff(x 1 ,x )< aff s
Through the entire paper one should keep in mind the following main modeling assumptions:

Modeling assumption 1 In our simplified mathematical model we do not take into account any spatial factor and in a single time step a GC B-cell can eventually undergo both division (with mutation) and selection. Hence the time unit has to be chosen big enough to take into account both mechanisms.

Modeling assumption 2 In this paper we are considering discrete-time models.

The symbol t always denote a discrete time step, hence it is an integral value.

We will refer to t as time, generation, or even maturation cycle to further stress the fact that in a single time interval [t, t + 1] each B-cell within the GC population is allowed to perform a complete cycle of division, mutation and selection.

Modeling assumption 3 Throughout the entire paper, when we talk about death rate (respectively division rate or selection rate) we are referring to the probability that each cell has of dying (respectively dividing or being submitted to selection) in a single time step.

Results

In this Section we formalize mathematically the model introduced above. This enables the estimation of various qualitative and quantitative measures of the GC evolution and of the selected pool as well. In Section 3.1 we show that a simple GW process describes the evolution of the size of the GC and determine a condition for its extinction. In order to do this we do not need to know the mutational model. Nevertheless, if we want to understand deeply the whole reaction we need to consider a (N + 3)-type GW process, which we introduce in Section 3.2. Therefore we determine explicitly other quantities, such as the average affinity in the GC and the selected pool, or the evolution of the size of the latter. We conclude this section by numerical simulations (Section 3.4).

Evolution of the GC size

The aim of this section is to estimate the evolution of the GC size and its extinction probability. In order to do so we define a simple GW process, with respect to parameters r d , r div and r s . Indeed, each B-cell submitted to selection exits the GC pool, independently from its affinity with x. Hence we apply some classical results about generating functions and GW processes ([14], Chapter I), which we recall in Appendix A. Proposition 1 gives explicitly the expected size of the GC at time t and conditions for the extinction of the GC.

Definition 4 Let Z (z 0 ) t
, t ≥ 0 be the random variable (rv) describing the GCpopulation size at time t, starting from z 0 ≥ 1 initial B-cells. (Z

(z 0 ) t ) t∈N is a
Markov Chain (MC) -since each cell behaves independently from the others and from previous generations -on {0, 1, 2, . . . }.

If z 0 = 1 and there is no confusion, we denote Z t := Z

(1) t . By Definition 4, Z 1 corresponds to the number of cells in the GC at the first generation, starting from a single seed cell. Thanks to Definition 3 one can claim that Z 1 ∈ {0, 1, 2}, with the following probabilities:

   p 0 := P(Z 1 = 0) = r d + (1 -r d )r s (1 -r div + r div r s ) p 1 := P(Z 1 = 1) = (1 -r d )(1 -r s )(1 -r div + 2r div r s ) p 2 := P(Z 1 = 2) = r div (1 -r d )(1 -r s ) 2 (1)
As far as next generations are concerned, conditioning to Z t = k, i.e. at generation t there are k B-cells in the GC, Z t+1 is distributed as the sum of

k independent copies of Z 1 : P(Z t+1 = k | Z t = k) = P k i=1 Z 1 = k .
Equalities in (1) are derived by identifying the events leading to 0, 1 or 2 offspring in the GC coming from a single clone. Since these events are independent and disjoint, the result follows. For instance there will be 0 new individuals in the GC if either the mother cell dies, or it does not die, does not divide and is submitted to selection, or it does not die, it does divide and both daughter cells are submitted to selection :

P(Z 1 = 0) = P death ∪ (death C ∩ division C ∩ selection) ∪ (death C ∩ division ∩ selection ∩ selection) = P(death) + P(death C )P(division C )P(selection) + P(death C )P(division)P(selection) 2 = r d + (1 -r d )(1 -r div )r s + (1 -r d )r div r 2 s
We have denoted by A C the complement of A. Expressions for p 1 and p 2 are obtained proceeding as before.

Definition 5 Let X be an integer valued rv, p k := P(X = k) for all k ≥ 0. Its probability generating function (pgf) is given by:

F X (s) = +∞ k=0 p k s k
The pgf for Z 1 :

F (s) = p 0 + p 1 s + p 2 s 2 = r d + (1 -r d )r s (1 -r div + r div r s ) + (1 -r d )(1 -r s )(1 -r div + 2r div r s )s + r div (1 -r d )(1 -r s ) 2 s 2 (2)
By using classical results on Galton-Watson processes (see Appendix A), one can prove:

Proposition 1

(i) The expected size of the GC at time t and starting from z 0 initial B-cells is given by:

E(Z (z 0 ) t ) = z 0 ((1 -r d )(1 + r div )(1 -r s )) t (3)
(ii) Denoted by η z 0 the extinction probability of the GC population starting from z 0 initial B-cells, one has:

-if E(Z (1) 1 ) ≤ 1 ⇔ r s ≥ 1 - 1 (1 -r d )(1 + r div )
, then η z 0 = 1: the process is subcritical otherwise η z 0 = η z 0 < 1, η being the smallest fixed point of (2): the process is supercritical

In particular, the initial number of seed cells z 0 does not affect the criticality of the process. Nevertheless, in the supercritical case, increasing the number of seed B-cells at the beginning of the process makes the probability of extinction decrease. More precisely, in the case η < 1, then η z 0 → 0 if z 0 → ∞, but we recall that GCs seem to be typically seeded by few B-cells, varying from ten to hundreds [START_REF] Tas | Visualizing antibody affinity maturation in germinal centers[END_REF]. This section shows that a classical use of a simple GW process enables to understand quantitatively the GC growth. Moreover, Proposition 1 (ii)

gives a condition on the main parameters for the extinction of the GC: if the selection pressure is too high, with probability 1 the GC size goes to 0, independently from the initial number of seed cells. Intuitively, a too high selection pressure prevents those B-cells with bad affinity to improve their fitness undergoing further rounds of mutation and division. Most B-cells will be rapidly submitted to selection, hence either exit the GC as output cells or die by apoptosis if they fail to receive positive selection signals [START_REF] Maclennan | B-cell memory and the persistence of antibody responses[END_REF]. In Figure 2 we plot the extinction probability of a GC initiated from a single seed cell as a function of r s (r d and r div are fixed), in order to stress the presence of a threshold effect of the selection probability over the extinction probability.

The extinction probability of the GC process can give us some further insights on factors which are potentially involved in determining the success or failure of a GC reaction. This simplified mathematical model suggests that if the selection pressure is too high compared to the division rate (c.f. due to Tfh signals in the LZ), the GC will collapse with probability 1, preventing the generation of high affinity antibodies against the presented antigen, hence an efficient immune response.

3.2 Evolution of the size and fitness of GC and selected pools

The GW process defined in the previous Section only describes the size of the GC. Indeed, we are not able to say anything about the average fitness of GC clones, or the expected number of selected B-cells, or their average affinity. Hence, we need to consider a more complex model and take into account the threshold for positive selection a s , and the transition probability matrix characterizing the mutational rule. Indeed, the mutational process is described as a Random Walk (RW) on the state space {0, . . . , N } of affinity classes indices. The mutational rule reflects the edge set associated to the state-space {0, . . . , N }: this is given by a transition probability matrix.

Definition 6 Let (X t ) t≥0 be a RW on the state-space of B-cell traits describing a pure mutational process of a B-cell during the GC reaction. We denote by Q N = (q ij ) 0≤i,j≤N the transition probability matrix over {0, . . . , N } which gives the probability of passing from an affinity class to another during the given mutational model. For all 0 ≤ i, j ≤ N :

q ij = P(a x (X t+1 ) = j | a x (X t ) = i)
From a biological point of view, these probabilities could be obtained e.g. by identifying which key mutations are the most relevant in determining changes in the fitness of a clone to a specific antigen and at which frequency they are produced.

We introduce a multi-type GW Process (see for instance [START_REF] Athreya | Branching processes[END_REF], chapter V).

i = (i 0 , . . . , i N , 0, 0). Let m ij := E[Z (i)
1,j ] the expected number of offspring of type j of a cell of type i in one generation. We collect all m ij in a matrix, M = (m ij ) 0≤i,j≤N +2 .

We have:

E[Z (i) t ] = iM t (4)
Supposing matrix Q N given (Definition 6), describing the probability to switch from one affinity class to another thanks to a single mutation event, one can explicitly derive the elements of M.

Proposition 2 M is a (N + 3) × (N + 3) matrix defined as a block matrix:

M = M 1 M 2 0 2×(N +1) I 2
Where:

-0 2×(N +1) is a 2 × (N + 1) matrix with all entries 0; -I n is the identity matrix of size n;

-M 1 = 2(1 -r d )r div (1 -r s )Q N + (1 -r d )(1 -r div )(1 -r s )I N +1 -M 2 = (m 2,ij
) is a (N + 1) × 2 matrix where for all i ∈ {0, . . . , N }:

-if i ≤ a s : m 2,i1 = (1 -r d )(1 -r div )r s + 2(1 -r d )r div r s a s j=0 q ij , m 2,i2 = r d + 2(1 -r d )r div r s N j=a s +1 q ij -if i > a s : m 2,i1 = 2(1 -r d )r div r s a s j=0 q ij , m 2,i2 = r d + (1 -r d )(1 -r div )r s + 2(1 -r d )r div r s N j=a s +1 q ij
The proof of Proposition 2 is available in Appendix B. It is based on the computation of the probability generating function of Z 1 .

Remark 1 Independently from the given mutational model, the expected number of selected or dead B-cells that each GC B-cell can produce in a single time step is given by α := r d + (1 -r d )(1 + r div )r s . All rows of M 2 sum to α independently from the probability that each clone submitted to selection has of being positive selected, which we recall is 1 if it belongs to the i th affinity class, i ≤ a s , zero otherwise.

Of course in the multi-type context we recover again results from Section 3.1, such as the extinction probability of the GC (detailed in Appendix C).

In order to determine the expected number of selected cells at a given time t, we need to introduce another multi-type GW process.

Definition 8 Let Z (i) t = ( Z (i) t,0 , . . . , Z (i) t,N +2 ), t ≥ 0 be a MC where for all 0 ≤ j ≤ N , Z (i)
t,j describes the number of GC B-cells belonging to the j th -affinity class with respect to x, Z (i) t,N +1 the number of selected B-cells and Z (i) t,N +2 the number of dead B-cells at generation t, when the process is initiated in state i = (i 0 , . . . , i N , 0, 0) and before the selection mechanism is performed for the t th -generation.

Proceeding as we did for Z (i) t , we can determine a matrix M whose elements

are m ij := E[ Z (i) 1,j ] for all i, j ∈ {0, . . . , N + 2}.
Proposition 3 M is a (N + 3) × (N + 3) matrix, which only depends on matrix Q N , r d and r div and can be defined as a block matrix as follows:

M = M 1 M 2 0 2×(N +1) I 2
Where:

-M 1 = 2(1 -r d )r div Q N + (1 -r d )(1 -r div )I N +1 -M 2 = (0 N +1 , r d • 1 N +1 ), where 0 N +1 (resp. 1 N +1 ) is a (N + 1)-column
vector whose elements are all 0 (resp. 1).

One could prove that:

E Z (i) t = iM t-1 M (5)
Proposition 4 Let i be the initial state, |i| its 1-norm (|i|

:= N +2 j=0 i j ).
-The expected size of the GC at time t:

N k=0 (iM t ) k = |i| ((1 -r d )(1 + r div )(1 -r s )) t (6)
-The average affinity in the GC at time t:

N k=0 (N -k)(iM t ) k N k=0 (iM t ) k (7)
-Let S t , t ≥ 1 denotes the random variable describing the number of selected B-cells at time t. By hypothesis S 0 = 0. (S t ) t∈N is a MC on {0, 1, 2, . . . }.

The expected number of selected B-cells at time t, t ≥ 1:

E(S t ) = r s a s k=0 iM t-1 M k ( 8
)
-The expected number of selected B-cells produced until time t:

E t n=0 S n = E Z (i) t N +1 = iM t N +1 (9)
-The average affinity of selected B-cells at time t, t ≥ 1:

a s k=0 (N -k) iM t-1 M k a s k=0 iM t-1 M k ( 10
)
-The average affinity of selected B-cells until time t:

r s t n=1 a s k=0 (N -k) iM n-1 M k (iM t ) N +1 (11)
Proof Equations ( 6) and ( 9) are a direct application of what stated in Equation [START_REF] Kauffman | The nk model of rugged fitness landscapes and its application to maturation of the immune response[END_REF]. Indeed, Equation [START_REF] Kauffman | The nk model of rugged fitness landscapes and its application to maturation of the immune response[END_REF] states that iM t contains the expectation of the number of all types cells at generation t when the process is started in i. Hence the expectation of the size of the GC at the t th generation is given by N k=0 (iM t ) k , since the GC at generation t contains all alive non-selected B-cells, irrespectively from their affinity. Similarly, the expected number of selected B-cells untill time t (9) corresponds to the expectation of the (N + 1) thtype cell, iM t N +1 .

The proof of Equation ( 8) is based on Equation ( 5), which allows to estimate the number of GC B-cells at generation t which are susceptible of being challenged by selection. One can remark that the expected number of selected B-cells at time t is obtained from the expected number of B-cells in GC at time t (before the selection mechanism is performed) having fitness good enough to be positive selected. This is given by a s k=0 iM t-1 M k , thanks to [START_REF] Balelli | Random walks on binary strings applied to the somatic hypermutation of b-cells[END_REF]. The result follows by multiplying this expectation by the probability that each of these B-cells is submitted to mutation, i.e. r s . Finally, results about the average affinity in both the GC and the selected pool (Equations ( 7), ( 10) and ( 11)) are obtained from the previous ones (c.f. ( 6), ( 8) and ( 9)) by multiplying the number of individuals belonging to the same class by their fitness (Definition 2), and dividing by the total number of individuals in the considered pool. The definition of affinity as a function of the affinity classes, determines Equations ( 7), [START_REF] De Silva | Dynamics of b cells in germinal centres[END_REF] and [START_REF] Dighiero | Chronic lymphocytic leukaemia[END_REF]. Indeed, the affinity of the k th -affinity class is given by N -k.

Remark 2

The expected size of the GC at time t can be obtained applying a simple GW process (Section 3.1) and is given by [START_REF] Athreya | Branching processes[END_REF]. It is possible to prove the equality in brackets in Equation ( 6) starting from the (N + 3)-type GW process. The interested reader can address to Appendix D for the detailed proof.

Optimal value of r s maximizing the expected number of selected B-cells at time t

What is the behavior of the expected number of selected B-cells as a function of the model parameters? In particular, is there an optimal value of the selection rate which maximizes this number? In this section we show that, indeed, the answer is positive.

To do so we detail hereafter the computation of E(S t ) (Equation ( 8)), given by Proposition 4.

Let us suppose, for the sake of simplicity, that Q N is diagonalizable:

Q N = RΛ N L , ( 12
)
where

Λ N = diag(λ 0 , . . . , λ N ), and R = (r ij ) (resp. L = (l ij )
) is the transition matrix whose rows (resp. lines) contain the right (resp. left) eigenvectors of

Q N , corresponding to λ 0 , . . . , λ N .
Proposition 5 Let us suppose that at t = 0 there is a single B-cell entering the GC belonging to the i th -affinity class with respect to the target cell. Moreover, let us suppose that Q N = RΛ N L. For all t ∈ N, the expected number of selected B-cells at time t, is:

E(S t ) = r s (1 -r s ) t-1 (1 -r d ) t N =0 (2λ r div + 1 -r div ) t a s k=0 r i l k ,
The proof of Proposition 5 is detailed in Appendix E.

As an immediate consequence of Proposition 5, we can claim:

Proposition 6 For all t * ∈ N fixed, the value r * s := r s (t * ) which maximizes the expected number of selected B-cells at the t * th maturation cycle is:

r * s = 1 t * Proof Since (1 -r d ) t N =0 (2λ r div + 1 -r div ) t a s k=0 r i l k
is a non negative quantity independent from r s , the value of r s which maximizes E(S t * ) is the one that maximizes r s (1 -r s ) t * -1 . The result trivially follows.

This result suggests that the selection rate in GCs is tightly related to the timing of the peak of a GC response, i.e. the timing corresponding to the maximal production of output cells (this timing can be determined e.g. by observing the concentration in blood of produced specific B-cells after infection or vaccination). In particular, following this model, GCs which peak early (e.g. for whom the maximal output cell production is reached in a few days) are possibly characterized by a higher selection pressure than GCs peaking later.

The peak of a typical GC reaction, measured as the average GC volume, has been estimated to be close to day 12 post immunization or a few days before [START_REF] Wollenberg | Regulation of the germinal center reaction by foxp3+ follicular regulatory t cells[END_REF], which is consistent with the observation of plasma cell response peak after immunization, e.g. [START_REF] Mitchell | Polysaccharide-specific b cell responses to vaccination in humans[END_REF]. Moreover, an high selection rate could also prevent a correct and efficient establishment of an immune response (c.f. results about extinction probability -Proposition 1). In addition, from a biological viewpoint, a too demanding selection pressure could avoid the generation of advantageous mutations, hence their fixation.

Numerical simulations

We evaluate numerically results of Proposition 4. The (N + 3)-type GW process allows a deeper understanding of the dynamics of both populations: inside the GC and in the selected pool. Through numerical simulations we emphasize the dependence of the quantities defined in Proposition 4 on parameters involved in the model.

In previous works [START_REF] Balelli | Random walks on binary strings applied to the somatic hypermutation of b-cells[END_REF][START_REF] Balelli | Branching random walks on binary strings for evolutionary processes[END_REF] we have modeled B-cells and antigens as N -length binary strings, hence their traits correspond to elements of {0, 1} N . In this context we have characterized affinity using the Hamming distance between B-cell and antigen representing strings. The idea of using a N -dimensional shape space to represent antibodies traits and their affinity with respect to a specific antigen has already been employed (e.g. [START_REF] Perelson | Theoretical studies of clonal selection: minimal antibody repertoire size and reliability of self-non-self discrimination[END_REF][START_REF] Meyer-Hermann | An analysis of b cell selection mechanisms in germinal centers[END_REF][START_REF] Kauffman | The nk model of rugged fitness landscapes and its application to maturation of the immune response[END_REF]), and N typically varies from 2 to 4. In the interests of simplification, we chose to set N = 2.

Moreover, from a biological viewpoint, this choice means that we classify the amino-acids composing B-cell receptors strings into 2 classes, which could represent amino-acids negatively and positively charged respectively. Charged and polar amino-acids are the most responsible in creating bonds which determine the antigen-antibody interaction [START_REF] Murphy | Janeway's immunobiology[END_REF].

While performing numerical simulations (Sections 3.4 and 4.2) we refer to the following transition probability matrix on {0, . . . , N }: Definition 9 For all i, j ∈ {0, . . . , N }:

q ij = P(a x (X t+1 ) = j | a x (X t ) = i) =    i/N if j = i -1 (N -i)/N if j = i + 1 0 if |j -i| = 1 Q N := (q ij ) 0≤i
,j≤N is a tridiagonal matrix where the main diagonal consists of zeros.

If we model B-cell traits as vertices of the state-space {0, 1} N , this corresponds to a model of simple point mutations (see [START_REF] Balelli | Random walks on binary strings applied to the somatic hypermutation of b-cells[END_REF] for more details and variants of this basic mutational model on binary strings).

Example 1 One can give explicitly the form of matrix M 2 (Proposition 2)

corresponding to the mutational model defined in Definition 9:

M 2 = 0 . . . a s -1 a s a s + 1 a s + 2 . . . N               α r d . . . . . . α r d α -β + β a s N r d + β N -a s N β a s +1 N r d + α -β + β N -(a s +1) N 0 r d + α . . . . . . 0 r d + α              
, where:

-α := (1 -r d )(1 + r div )r s -β := 2(1 -r d )r div r s
Indeed, due to the particular form of matrix Q N one has straightforward:

-

a s j=0 q ij =          1 if i < a s a s /N if i = a s (a s + 1)/N if i = a s + 1 0 if i > a s + 1 - N j=a s +1 q ij =          0 if i < a s (N -a s )/N if i = a s (N -(a s + 1)) /N if i = a s + 1 1 if i > a s + 1
Remark 4 Note that all mathematical results obtained in previous sections are independent from the mutation model defined in Definition 9.

We suppose that at the beginning of the process there is a single B-cell entering the GC belonging to the affinity class a 0 . Of course, the model we set allows to simulate any possible initial condition. Indeed, by fixing the initial vector i, we can decide to start the reaction with more B-cells, in different affinity classes. When it is not stated otherwise, the employed parameter set for simulations is given in Table 1. We perform numerical simulations to better appreciate how the dynamics of the GC and positive selected clones populations are related and evolve depending on model parameters. In the case of a subcritical GC, by model definition selected clones stabilize at a given level once the GC becomes extinct. Hence we conveniently chose a parameter set (Table 1) which implies a supercritical GC (Proposition 1): with great probability the simulated GC goes through explosion, and so the selected population does.

Evolution of the GC population

The evolution of the size of the GC can be studied by using the simple GW process defined in Section 3.1. Equation ( 3), in the case of a single initial B-cell, evidences that the expected number of B-cells within the GC for this model only depends on r d , r div and r s and it is not driven by the initial affinity, nor by the threshold chosen for positive selection a s , nor by the mutational rule.

Equation (3) evidences that, independently from the transition probability matrix defining the mutational mechanism, the GC size at time t increases with r div and decreases for increasing r s and r d . Moreover, the impact of these last two parameters is the same for the growth of the GC. One could expect this behavior since the effect of both the death and the selection on a B-cell is the exit from the GC.

In order to study the evolution of the average affinity within the GC, we need to refer to the (N + 3)-type GW process defined in Section 3.2.

Proposition 7

Let us suppose that Q N = RΛ N L. The average affinity within the GC at time t, starting from a single B-cell belonging to the i th -affinity class with respect to x is given by:

N - N =0 (2λ r div + 1 -r div ) t N k=0 k • r i l k (1 + r div ) t ,
Proof It follows directly from Equations [START_REF] Castro | Learning and optimization using the clonal selection principle[END_REF] and by considering the eigendecomposition of matrix Q. One has to consider the expression of the t th power of matrix M (which can be obtained recursively, see Appendix E): one can Aff(GC) prove that the first N + 1 components of the i th -row of matrix M t are the elements of the i th -row of matrix RD t L, where

at t =15 a 0 =2 a 0 =3 a 0 =4 a 0 =5 (
D = 2(1 -r d )r div (1 -r s )Λ N + (1 -r d )(1 -r div )(1 -r s )I N +1 is a diagonal matrix.
It is obvious from Proposition 7 that this quantity only depends on the initial affinity with the target trait, the transition probability matrix Q N and the division rate r div . The average affinity within the GC does not depend on a s (as one can clearly see in Figure 3

Evolution of the selected pool

The evolution of the number of selected B-cells during time necessarily depends on the evolution of the GC. In particular, let us suppose we are in the supercritical case, i.e. the extinction probability of the GC is strictly smaller than 1. Than, with positive probability, the GC explodes and so does the selected pool. On the other hand, if the GC extinguishes, the number of selected B-cells will stabilize at a constant value, as once a B-cell is selected it can only stay unchanged in the selected pool.

As demonstrated in Section 3.3, there exists an optimal value of the parameter r s which maximizes the expected number of selected B-cells at time for different values of a 0 , depending on r s and a s respectively. There exists an optimal value of r s maximizing the expected number of selected B-cells for a given generation. This value is independent from a 0 and is equal to 1/t as demonstrated in Proposition 6: the red vertical line in (a) corresponds to this value. (c) Comparison between the expected number of selected B-cells until time t given by evaluation of the theoretical formula (Equation ( 9)), and the empirical value obtained as the mean over 4000 simulations. Vertical bars denotes the corresponding estimated standard deviations. Here N = 7 and r s = 0.3. Proposition 8 M + is a (N + 3) 2 matrix, which we can define as a block matrix in the following way:

t.
a 0 =2 a 0 =3 a 0 =4 a 0 =5 (
Ø x x x 1 x 2 x 1 x 1 t t+1 r d 1-r d r div 1-r div r s 1-r s aff(x 1 ,x )≥ aff s GC GC Sel aff(x 1 ,x )< aff s (a) Positive selection Ø x x x 1 x 2 Ø x 1 t t+1 r d 1-r d r div 1-r div
M + = M + 1 M + 2 0 2×(N +1) I 2
Where:

-M + 1 = (m + 1,ij
) is a (N + 1) 2 matrix. For all i ∈ {0, . . . , N }:

-∀ j ≤ a s : m + 1,ij = 2(1 -r d )r div (1 -r s )q ij + (1 -r d )(1 -r div )(1 -r s )δ ij -∀ j > a s : m + 1,ij = 2(1 -r d )r div q ij + (1 -r d )(1 -r div )δ ij
where δ ij is the Kronecker delta.

-M + 2 = (m + 2,ij
) is a (N + 1) × 2 matrix where for all i ∈ {0, . . . , N }, m + 2,i1 = m 2,i1 , and m + 2,i2 = r d . We recall that m 2,i1 is the i th -component of the first column of matrix M 2 , given in Proposition 2.

Negative selection

Definition 13 Let Z - t (i) = (Z - t,0 (i) , . . . , Z - t,N +1 (i)
), t ≥ 0 be a MC where for all 0 ≤ j ≤ N , Z - t,j (i) describes the number of GC B-cells belonging to the j th -affinity class with respect to x and Z - t,N +1 (i) the number of dead B-cells at generation t, when the process is initiated in state i = (i 0 , . . . , i N , 0), and following the evolutionary model described by Definition 11.

Let us denote by M -= (m - ij ) 0≤i,j≤N +1 the matrix containing the expected number of type-j offspring of a type-i cell corresponding to the model defined by Definition 13.

Proposition 9 M -is a (N + 2) 2 matrix, which we can define as a block matrix in the following way:

M -= M - 1 m - 2 0 N +1 1
Where:

-M - 1 = (m - 1,ij
) is a (N + 1) 2 matrix. For all i ∈ {0, . . . , N }: -0 N +1 is a (N + 1) row vector composing of zeros.

-∀ j ≤ a s : m - 1,ij = 2(1 -r d )r div q ij + (1 -r d )(1 -r div )δ ij -∀ j > a s : m - 1,ij = 2(1 -r d )r div (1 -r s )q ij + (1 -r d )(1 -r div )(1 -r s )δ ij -m - 2 is a (N +
We do not prove Propositions 8 and 9, since the proofs are the same as for Proposition 2 (Appendix B).

Results stated in Proposition 4 hold true for these new models, by simply replacing matrix M with M + (resp. M -). Of course, in the case of negative selection, as we do not consider the selected pool, we only refer to ( 6) and ( 7)

quatifying the growth and average affinity of the GC. Matrix M is the same for both models as only selection principles change.

Because of peculiar structures of matrices M + and M -, we are not able to compute explicitly their spectra. Henceforth we can not give an explicit formula for the extinction probability or evaluate the optimal values of the selection rate r s as we did in Sections 3.2 and 3.3.

Nevertheless, by using standard arguments for positive matrices, the greatest eigenvalue of both matrices M + 1 and M - 1 can be bounded, and hence give sufficient conditions for extinction. Indeed, form classical results about multi-type GW processes, the value of the greatest eigenvalue allows to discriminate between subcritical case (i.e. extinction probability equal to 1) and supercritical case (i.e. extinction probability strictly smaller than 1) [START_REF] Athreya | Branching processes[END_REF].

Proposition 10 Let q + (resp. q -) be the extinction probability of the GC for the model corresponding to matrix M + 1 (resp. M - 1 ).

-If r div ≤ r d 1 -r d , then q + = q -= 1. -If r s < 1 - 1 (1 -r d )(1 + r div )
, then q + < 1 and q -< 1.

Proof Since both matrices M + 1 and M - 1 are strictly positive matrices, the Perron Frobenius Theorem insures that the spectral radius is also the greatest eigenvalue. Then the following classical result holds [START_REF] Minc | Nonnegative matrices[END_REF]:

Theorem 1 Let A = (a ij ) be a square nonnegative matrix with spectral radius ρ(A) and let r i (A) denote the sum of the elements along the i th -row of A. Then:

min i r i (A) ≤ ρ(A) ≤ max i r i (A)
Simple calculations provide:

min i r i (M + 1 ) = (1 -r d )(1 + r div ) -r s (1 -r d )   2r div min i a s j=0 q ij + 1 -r div   max i r i (M + 1 ) = (1 -r d )(1 + r div ) -2r s r div (1 -r d ) max i a s j=0 q ij min i r i (M - 1 ) = (1 -r d )(1 + r div ) -r s (1 -r d )   2r div min i N j=a s +1 q ij + 1 -r div   max i r i (M - 1 ) = (1 -r d )(1 + r div ) -2r s r div (1 -r d ) max i N j=a s +1 q ij
The result follows by observing that for all i ∈ {0, . . . , N }, 0 ≤ a s j=0 q ij , N j=a s +1 q ij ≤ 1, and applying Theorem 3. selection is strictly smaller than 1 (Proposition 1). Indeed keeping the same parameters for all models, the size of the GC for the model of positive and negative selection is smaller than the size of GCs corresponding to both models of only positive and only negative selection. Consequently if the GC corresponding to M has a positive probability of explosion, it will be necessarily the same for M + and M -.

Remark 7

The values of both ρ(M + 1 ) and ρ(M - 1 ) depend on a s , varying from a minimum of (1 -r d )(1 + r div )(1 -r s ) and a maximum of (1 -r d )(1 + r div ).

Figure 6 evidences the dependence on a s of the spectral radius of M + 1 and M - 1 , using matrix Q N given by Definition 9 as transition probability matrix.

Remark 7 and Figure 6 evidences that, conversely to the previous case of positive and negative selection, in both cases of exclusively positive (resp.

exclusively negative) selection the parameter a s plays an important role in the GC dynamics, affecting its extinction probability. In particular, keeping unchanged all other parameters, if a s → N (resp. a s → 0), then ρ(M + 1 ) (resp. In Figure 7 we plot the dependence of greatest eigenvalues of both matrices M + and M -with respect to r s . We fix r d = 0.1 and r div = 0.9 as for Figure 2. One can note that with this parameter set and if the threshold for positive selection a s is chosen not "too small" nor "too large" with respect to N , then the greatest eigenvalue for both matrices is always greater than 1 independently from r s , i.e. the extinction probability is always strictly smaller than 1. From a biological viewpoint we expect that a phisiological threshold for positive selection should not be too strict nor to weak. Indeed, a too demanding threshold for positive selection is not optimal since B-cells should have gained an extremely high affinity in order to be positive selected, which would at least require too much time, avoiding a prompt immune response against the invading pathogen. On the other hand, a too weak threshold results in an unchallenging affinity maturation process: almost any B-cell would be positive selected, irrespective from its affinity level with respect to the presented antigen. This could also entail the generation of auto-reactive clones.

ρ(M - 1 )) → (1 -r d )(1 + r div )(1 -r s ), which implies q + (resp. q -) → 1.

Numerical simulations

The evolution of GCs corresponding to matrices M + and M -respectively are complementary. Moreover, in both cases, keeping all parameters fixed one expects a faster expansion if compared to the model of positive and negative selection, since the selection acts only positively (resp. negatively) on good (resp. bad) clones. In particular, the model of negative selection corresponds to the case of 100% of recycling, meaning that all positively selected B-cells stay in the GC for further rounds of mutation, division and selection.

Figure 8 shows the dependence on a s of the GC size and fitness, comparing M + (left column) and M -(right column). Indeed, for these models the GC dynamics depends on the selection threshold, conversely to the previous case of positive and negative selection, and not only on the selection rate. The effects of a s on the GC are perfectly symmetric: it is interesting to observe that when both selection mechanisms are coupled, then a s does not affect the GC dynamics anymore, as shown for instance in Figure 3 (a). Moreover, Figures 8 (c,d) evidence the existence of a value of a s that minimizes (resp. maximizes) the expected average affinity in the GC for M + (resp. M -). In both cases this value is approximately N/2. This certainly depends on the transition probability matrix chosen for the mutational model, which converges to a binomial probability distribution over {0, . . . , N }.

The evolution of the selected pool for the model of positive selection have some important differences if compared to the model described in Section 3.

For instance, it is not easy to identify an optimal value of r s which maximizes the expected number of selected B-cells at time t. Indeed it depends both on a 0 and a s : if a 0 ≤ a s we find curves similar to those plotted in Figure 4 (a), otherwise Figure 9 (a) shows a substantial different behavior. Indeed, if a 0 > a s , choosing a big value for r s does not negatively affect the number of selected B-cells at time t. In this case, for the first time steps no (or a very few) B-cells will be positively selected, since they still need to improve their affinity to the Aff(GC) Aff(GC) respectively, raised to the power of t = 15 (see Figure 6). Note that thanks to Proposition 1 we know that for this parameter choice the expected size of the GC for the model of positive and negative selection corresponds to 15 , which is equivalently λ 15 max for a s = 10 in Figure 8 t. Moreover we observe important oscillations, which are probably due to the mutational model, and to the fact that the total GC size is still small for small t, since the process starts from a single B-cell. Nevertheless, it seems that for t big enough also in this case the value of r * s tends to approach 1/t.

at t =15 a 0 =2 a 0 =3 a 0 =4 a 0 =5 (
at t =15 a 0 =2 a 0 =3 a 0 =4 a 0 =5 (d)
((1 -r d )(1 + r div )(1 -r s ))
Since in the case of negative selection there is no selected pool, one can suppose that at a given time t the process stops and all clones in the GC pool exit the GC as selected clones. Hence it can be interesting to compare the selected pool of the model of positive selection and the GC pool of the model of negative selection at time t. Clearly to make these two compartments comparable, the main parameters of both systems have to be opportunely chosen. In Figure 10 we compare the size and average fitness of the selected pool for M + and the GC for M -at time t = 30. We test different values of the parameter r s . In particular, we observe that increasing r s the GC size for the model of negative selection decreases and its average fitness increases. For the parameter choices we made for these simulations, Figure 10 (a) shows that the size of the GC for M -is comparable to the size of the selected pool for M + at time t = 30 if, keeping all other parameters fixed, r s = 0.15 for M -. Nevertheless, this does not implies a comparable value for the average affinity: the clones of the selected pool for M + have a significantly greater average affinity than those of the GC for M -. In order to increase the average fitness in the GC for the model of negative selection one has to consider greater values for the parameter r s , but this affects the probability of extinction of the process. Aff(S30) , M + , r s =0.1 Aff(S30) , M + , r s =0.12 Aff(S30) , M + , r s =0.15 Aff(GC30) , M -, r s =0.1 Aff(GC30) , M -, r s =0.12 Aff(GC30) , M -, r s =0.15

(b)

Figure 10: (a) Expected number of B-cells which have been selected until time t = 30 for M + compared to the expected size of the GC for M -for different values of r s . (b) Expected corresponding average affinity for the selected pool (case of positive selection) and the GC (case of negative selection). For some choice of the parameter r s , the size of the selected pool for M + , and the GC for M -, are comparable. Nevertheless, the corresponding average affinities are significantly different.

We can expect this discrepancy between the average affinity for the selected pool for M + and the one of the GC for M -. Indeed, in the first case we are looking to all those B-cells which have been positive selected, hence belong at most to the a th s -affinity class. On the contrary in the case of M -, we consider the average affinity of all B-cells which are still alive in the GC at a given time step. Among these clones, if r s < 1, with positive probability there are also individuals with affinity smaller than the one required for escaping negative selection. They remain in the GC because they have not been submitted to selection. These B-cells make the average affinity decrease. Of course r s is not the only parameter affecting the quantities plotted in Figure 10. In particular, one can observe that choosing a greater value for a s also have a significant effect over the growth of both pools, as discussed in Remark 7.

Conclusions and perspectives

In this paper we formalize and analyze a mathematical model describing an evolutionary process with affinity-dependent selection. We use a multi-type GW process, obtaining a discrete-time probabilistic model, which includes division, mutation, death and selection. This is employed in the context of antibody affinity maturation in GCs. We believe that a probabilistic approach is well suited to the study of Darwinian-like processes such as the one taking place in GCs during an immune response. Indeed, this kind of approaches allows to better take into account local inhomogeneities related to the discrete nature of cells and stochastic fluctuations intrinsic to these processes, conversely to more popular deterministic continuum approaches. There, cell concentrations are described by a set of coupled ODEs changing deterministically and continuously during time, which has many computational advantages and has often been employed to model biological systems (e.g. [START_REF] Meyer-Hermann | A theory of germinal center b cell selection, division, and exit[END_REF][START_REF] Iber | A mathematical model for germinal centre kinetics and affinity maturation[END_REF][START_REF] Moreira | Modelling two possible mechanisms for the regulation of the germinal center dynamics[END_REF] for applications to the GC reaction). In the main model developed here, we choose a selection mechanism which acts both positively and negatively on individuals submitted to selection. This choice is motivated by the fact that there are biological evidence supporting both kind of selection mechanisms: positive affinity-based selection by antigen binding as well as selection-dependent apoptosis [START_REF] Victora | Germinal centers[END_REF][START_REF] Inoue | Generation of memory b cells and their reactivation[END_REF].

The simplified mathematical framework proposed here allow to investigate how different kind of B-cell population evolves during the immune response both in the initial explosion phase and in the later relaxation phase of typical GCs.

Indeed, mathematical analysis of the model leads to build matrix M, which contains the expectations of each type (Proposition 2) and enables to describe the average behavior of all components of the process. Moreover, thanks to the spectral decomposition of M we were able to obtain explicitly some formulas giving the expected dynamics of all types. In addition, we exhibited an optimal value of the selection rate maximizing the expected number of selected clones for the t th -generation (Proposition 6). This is one possible choice of the selection mechanism. From a mathematical point of view, matrix M is particularly easy to manipulate, as we can obtain explicitly its spectrum. On the other hand, the positive and negative selection model leads, for example, to a selection threshold that does not have any impact on the evolution of the GC size. From a biological point of view this seems counterintuitive, since we could expect that the GC dynamics is sensible to the minimal fitness required for positive selection. Moreover, this process does not take into account any recycling mechanism, which has been confirmed by experiments [START_REF] Victora | Germinal center dynamics revealed by multiphoton microscopy with a photoactivatable fluorescent reporter[END_REF] and which improves GCs' efficiency. In addition, we considered that only the selection mechanism is affinity dependent, while in the GC reaction other mechanisms, such as the death and proliferation rate, may depend on fitness [START_REF] Gitlin | Clonal selection in the germinal centre by regulated proliferation and hypermutation[END_REF][START_REF] Anderson | Taking advantage: high-affinity b cells in the germinal center have lower death rates, but similar rates of division, compared to low-affinity cells[END_REF]. Of course it is possible to define models with affinity-dependent division and death mechanisms with our formalism. This would clearly lead to a more complicated model, which can be at least studied numerically. From a biological viewpoint there exist many possibilities to improve the models proposed in this paper. First of all it is extremely important to fix the system parameters, which have to be consistent with the real biological process. The choice of N defines the number of affinity level with respect to a given antigen. This value can be interpreted in different ways. On the one hand it can correspond to the number of key mutations observed during the process of Antigen Affinity Maturation, hence be even smaller than 10. On the other hand, each mutational event implies a change in the B-cell affinity, slight

Mathematical tools used in

or not if it is a key mutation. In this case the affinity can be modeled as a continuous function, hence N corresponds to a possible discretization [START_REF] Weiser | Affinity maturation of b cells involves not only a few but a whole spectrum of relevant mutations[END_REF][START_REF] Xu | Key mutations stabilize antigenbinding conformation during affinity maturation of a broadly neutralizing influenza antibody lineage[END_REF].

To this choice corresponds an appropriate choice of the transition probability matrix defining the mutational model over the affinity classes, Q N . In most numerical simulations we set N = 10, which is a sensible value since experimentalists observe that high-affinity B-cells differ in their BCR coding gene by about 9 mutations from germline genes [START_REF] Iber | A mathematical model for germinal centre kinetics and affinity maturation[END_REF][START_REF] Zhang | Optimality of mutation and selection in germinal centers[END_REF]. Nevertheless all mathematical results are independent from this choice and hold true for all N ≥ 1. The selection, division and death rates have also an important impact in the GC and selected pool dynamics: in the simulations we set them in order to be in a case of explosion of the GC hence appreciate the effects of all parameters over the main quantities, but they are not biologically justified. For instance, the typical proliferation rate of a B-cell has been estimated between 2 and 4 per day and in the literature we found B-cell death rates of the order of 0.5-0.8

per day [START_REF] Meyer-Hermann | An analysis of b cell selection mechanisms in germinal centers[END_REF][START_REF] Zhang | Optimality of mutation and selection in germinal centers[END_REF][START_REF] Keşmir | A mathematical model on germinal center kinetics and termination[END_REF]. Hence, if we suppose that a single time step corresponds e.g. to 6 hours, a consistent proliferation rate would be r div 0.75, while the death rate r d should be around 0.175. Since over a 6 hours period about 50% of B-cells transit from the DZ to the LZ, where they compete for positive selection signaling [START_REF] Bannard | Germinal center centroblasts transition to a centrocyte phenotype according to a timed program and depend on the dark zone for effective selection[END_REF][START_REF] Victora | Snapshot: the germinal center reaction[END_REF], we should choose r s ≤ 0.5. It could be further characterized taking into account its tightly relation with the time of GC peak, as highlighted in Section 3.3.

In Section 3.3 we have explicitly determined the optimal value of the selection rate maximizing the production of output cells at time t for the main model of positive and negative selection. It is equal to 1/t independently from all other parameters. Moreover, numerical estimations for the model of positive selection (Section 4.2) suggest that also in this case there exists an optimal value of r s (t), which tends to 1/t at least for t big enough. One has to interpret this result as the ideal optimal strength of the selection pressure to obtain a peak of the GC production of output cells at a given time step. For example, let us suppose again that a time step corresponds to 6 hours. The peak of the GC reaction has been measured to be close to day 12 [START_REF] Wollenberg | Regulation of the germinal center reaction by foxp3+ follicular regulatory t cells[END_REF], i.e. after ∼ 48 maturation cycles in our model: for the kind of models we built and analyzed in this paper, a constant selection pressure r s of 1/48 0.02 assures that the production of plasma and memory B-cells at the GC peak is maximized. Note that with the parameter choice r d = 0.175, r div = 0.75 and r s = 0.02, the extinction probability of the GC is 0.3 z 0 , z 0 being the number of initial seed cells. Since the extinction probability is strictly smaller than 1, such a GC will explode with high probability and will be able to assure an intense and efficient immune response.

The particular form 1/t of the optimal selection rate for the t th generation obtained in Section 3. by an increased frequency of replacements, while a decreased frequency indicates negative selection. According to [START_REF] Yaari | Quantifying selection in high-throughput immunoglobulin sequencing data sets[END_REF], the selection strength seems to vary and also switch from positive to negative in a different way depending on the location, i.e. if we are looking to complementary determining regions (CDRs), which are more significant for functional selection, or to framework regions. This gives stronger motivation to analyse both kind of selection mechanisms, acting both separately and simultaneously, and observe their effects over affinity maturation, as we have done in this paper using our simplified mathematical framework.

In our models the selection pressure is constant. Since the optimal selection rate above depends on time, this suggests to go further in this direction. Moreover, a time-dependent selection pressure would allow to take into account, for instance, the early GC phase in which simple clonal expansion of B-cells with no selection occurs [START_REF] De Silva | Dynamics of b cells in germinal centres[END_REF]. The hypothesis of a selection pressure changing over time can be easily integrated in our model. Indeed let us suppose that a selection rate r s,1 until time t 1 and r s,2 for all t > t 1 are fixed. Starting from the initial condition i the expectations of each type at time t are given by (iM t r s,1 )

if t ≤ t 1 and (iM t 1 r s,1 M t-t 1 r s,2 ) if t > t 1
, where M r s,i is the matrix containing the expectations of each type for an evolutionary process with constant selection rate r s,i , i = 1, 2. In Figure 11 we plot the expected evolution during time of all types considering an increasing selection rate. We evaluate the expectations of all types following a process with positive and negative selection. We set r s = 0 until t = 5, r s = 0.1 from t = 6 to t = 15 and r s = 0.3 for t > 15.

Numerical simulations show that a time dependent selection rate allows initial explosion of the GC, and then progressive extinction, while when parameters are fixed, a GW process gives only rise either to explosion or to extinction, as shown above. The regulation and termination of the GC reaction has not yet been fully understood. In the literature, an increasing differentiation rate of GC B-cells is thought to be a good explanation [START_REF] Moreira | Modelling two possible mechanisms for the regulation of the germinal center dynamics[END_REF], here we show that other Evolution during time of the expected total size of the GC and the selected pool respectively, for the same set of parameters as in Figure 11 (a).

Z 0 Z 1 Z 2 Z 3 Z 4 Z 5 Z 6 Z 7 Z 8 Z 9 Z 10 Z 11 Z 12 Z 13 (
reasons could be of importance as well. Similarly, we can let other parameters vary for fixed time intervals, as well as decide to alternatively switch on and off the mutational mechanism, as already proposed in [START_REF] Perelson | Immunology for physicists[END_REF]. This can be obtained by alternatively use the identity matrix in place of Q N .

Applications of the models presented here to real biological problems and data should be further investigated. We propose here some contexts for which we believe that our kind of modeling approach could be employed to address biologically relevant questions.

Even if it is still extremely hard to have precise experimental information about the evolution of Antibody Affinity Maturation inside GCs, new refined techniques start to be available to measure clonal diversity in GCs. As an example, in [START_REF] Tas | Visualizing antibody affinity maturation in germinal centers[END_REF] the authors combine multiphoton microscopy and sequencing to understand how different clonal diversification patterns can lead to efficient affinity maturation. The models we propose could be used to infer which are reasonable mutational transitional probability matrices and selection mechanisms/pressure to obtain such different scenario and infer if the tendency of GC to go or not through homogenizing selection is solely due to the hazard or if this is dependent on the kind of antigenic challenge and/or some specific characteristics of the host. If this is the case, these results could be particularly relevant e.g. in the context of vaccination design, where we are interested in find new way to improve the quality of the immune response after vaccination challenge.

Another potential interesting application field is the study of some diseases entailing a dysfunction of the immune system, such as in particular Chronic Lymphocytic Leukemia (CLL), derived from antigen-experienced B-cells that differ in the level of mutations in their receptors [START_REF] Chiorazzi | Chronic lymphocytic leukemia[END_REF]. This is the commonest form of leukemia in the Western world [START_REF] Eichhorst | Chronic lymphocytic leukaemia: Esmo clinical practice guidelines for diagnosis, treatment and follow-up[END_REF]. In CLL, leukemia B-cells can mature partially but not completely, are unable to opportunely undergo mutations in GCs, and survive longer than normal cells, crowding out healthy B-cells. Prognosis varies depending on the ability of host B-cells to mutate their antibody gene variable region. Even if major progresses have been made in the identification of molecular and cellular markers predicting the expansion of this disease in patients, the pathology remains incurable [START_REF] Dighiero | Chronic lymphocytic leukaemia[END_REF][START_REF] Eichhorst | Chronic lymphocytic leukaemia: Esmo clinical practice guidelines for diagnosis, treatment and follow-up[END_REF]. Our modeling approach could be employed to understand how an "healthy" mutational matrix is modified in patients affected by CLL, and if other mechanisms could contribute to get the prognosis worse. This could eventually provide suggestions about the causes that lead to CLL, and motivation for further research on possible treatments.

Acknowledgements

This work was supported by the Labex inflamex, ANR project 10-LABX-0017. (ii) The expected size of the GC at time t and starting from z 0 B-cells is given by:

E(Z (z 0 ) t ) = z 0 (E(Zt)) = z 0 (E(Z 1 )) t , ( 13
)
(iii) η is the smallest fixed point of the generating function F , i.e. η is the smallest s s.t.

F (s) = s. (iv) If E(Z 1 ) =: m is finite, then:
if m ≤ 1 then F has only 1 as fixed point and consequently η = 1;

if m > 1 then F as exactly a fixed point on [0, 1[ and then η < 1.

(v) Denoted by ηz 0 the probability of extinction of (Z

(z 0 ) t
), one has:

ηz 0 = η z 0
where η is given by (iii).

Proposition 1 of Section 3.1 follows by applying Theorem 2 and Equation (1).

B Proof of Proposition 2

For all j ∈ {0, . . . , N + 2} the generating function of Z j gives the number of offspring of each type that a type j particle can produce. It is defined as follows:

f (j) (s 0 , . . . , s N +2 ) = k 0 ,...,k N +2 ≥0 p (j) (k 0 , . . . , k N +2 )s k 0 0 • • • s k N +2 N +2 , (14)
0 ≤ sα ≤ 1 for all α ∈ {0, . . . , N + 2}

where p (j) (k 0 , . . . , k N +2 ) is the probability that a type j cell produces k 0 cells of type 0, k 1 of type 1, . . . , k N +2 of type N + 2 for the next generation.

We denote:

-p(k) = (p (0) (k), . . . , p (N +2) (k)), for k = (k 0 , . . . , k N +2 ) ∈ Z N +3 + -f (s) = (f (1) (s), . . . , f (N +1) (s)), for s = (s 0 , . . . , s N +2 ) ∈ C N +3 := [0, 1] N +3
Then the probability generating function of Z 1 is given by:

f (s) = k∈Z N +3 + p(k)s k , s ∈ C N +3 (15)
Again, the generating function of Zt, ft(s), is obtained as the t th -iterate of f , and it holds true that:

f t+r (s) = ft[fr(s)], s ∈ C N +3 . Let m ij := E[Z (i) 1,j
] the expected number of offspring of type j of a cell of type i in one generation. We collect all m ij in a matrix, M = (m ij ) 0≤i,j≤N +2 . We have [START_REF] Athreya | Branching processes[END_REF]:

m ij = ∂f (i) ∂s j (1)
and:

E[Z (i) t,j ] = ∂f (i) t ∂s j (1) (16)
Finally:

E[Z (i) t ] = iM t (17)
One can explicitly derive the elements of matrix M for the process described in Definition 13.

Proposition M is a (N + 3) × (N + 3) matrix defined as a block matrix:

M = M 1 M 2 0 2×(N +1) I 2
Where:

-0 2×(N +1) is a 2 × (N + 1) matrix with all entries 0; -In is the identity matrix of size n;

-M 1 = 2(1 -r d )r div (1 -rs)Q N + (1 -r d )(1 -r div )(1 -rs)I N +1 -M 2 = (m 2,ij
) is a (N + 1) × 2 matrix where for all i ∈ {0, . . . , N }:

-if i ≤ as: m 2,i1 = (1 -r d )(1 -r div )rs + 2(1 -r d )r div rs as j=0 q ij , m 2,i2 = r d + 2(1 -r d )r div rs N j=as+1 q ij -if i > as: m 2,i1 = 2(1 -r d )r div rs as j=0 q ij , m 2,i2 = r d + (1 -r d )(1 -r div )rs + 2(1 -r d )r div rs N j=as+1 q ij
Proof One has to compute all f (i) (s) for i = 0, . . . , N + 2, which depend on r d , r div , rs, as and the elements of Q N . First, the elements of the (N + 2) th and (N + 3) th -lines are obviously determined: all selected (resp. dead) cells remain selected (resp. dead) for next generations, as they can not give rise to any other cell type offspring (we do not take into account here any type of recycling mechanism). Let i ∈ {0, . . . , N } be a fixed index: we evaluate m ij for all j ∈ {0, . . . , N + 2}. The first step is to determine the value of p (i) (k) for k = (k 0 , . . . , k N +2 ) ∈ Z N +3 + . There exists only a few cases in which p (i) (k) = 0, which can be explicitly evaluated:

p (i) (0, . . . , 0, 1) = r d if i ≤ as r d + (1 -r d )(1 -r div )rs otherwise p (i) (0, . . . , 0, 1, 0) = (1 -r d )(1 -r div )rs if i ≤ as 0 otherwise p (i) (0, . . . , 0, 1 i , 0, . . . , 0, 0) = (1 -r d )(1 -r div )(1 -rs)

-p (i) (0, . . . , 0, 2) = (1 -r d )r div r 2 s N j 1 =as+1 q ij 1 N j 2 =as+1
q ij 2 -p (i) (0, . . . , 0, 2, 0) = (1 -r d )r div r 2 s as

j 1 =0 q ij 1 as j 2 =0
q ij 2 Definition 18 We say that (Zt) is singular if each particle has exactly one offspring, which implies that the branching process becomes a simple MC.

Definition 19 Matrix M is said to be strictly positive if it has non-negative entries and there exists a t s.t. M t ij > 0 for all i, j. (Zt) is called positive regular iff M is strictly positive.

Notation 1 Let u, v ∈ R n . We say that u ≤ v if u i ≤ v i for all i ∈ {1, . . . , n}. Moreover, we say that u < v if u ≤ v and u = v.

Theorem 3 Let (Zt) be non singular and strictly positive. Let ρ be the maximal eigenvalue of M. The following three results hold:

1. If ρ < 1 (subcritical case) or ρ = 1 (critical case) then q = 1. Otherwise, if ρ > 1 (supercritical case), then q < 1.

2. lim t→∞ ft(s) = q, for all s ∈ C N +3 .

3. q is the only solution of f (s) = s in C N +3 .

The spectrum of matrix M defined in Definition 2 (and recalled in Appendix B) is obtained as follows:

Proposition 12 Let M be defined as a block matrix as in Proposition 2. Let λ M,i be its i th -eigenvalue. The spectrum of M is given by:

-For all i ∈ {0, . . . , N }, λ M,i = (1 -r d )(1 -rs)(1 + r div (2λ i -1)), where λ i is the i theigenvalue of matrix Q N .

whereas λ M,N +1 = 1 with multiplicity 2.

Proof As M is a block matrix with the lower left block composed of zeros, then Spec(M) = Spec(M 1 ) ∪ Spec(I 2 ). The result follows.

Therefore we obtain the same condition as in Proposition 1 for the extinction probability in the GC:

Proposition 13 Let q be the extinction probability for the process (Zt) defined in Definition 13 and restricted to the first N + 1 components ( i.e. we refer only to matrix M 1 , which defines the expectations of GC B-cells). Therefore:

-if rs ≥ 1 - 1 (1 -r d )(1 + r div )
, then q = 1 -otherwise q < 1 is the smallest fixed point of f (s) in C N +3 .

Proof Q N is a stochastic matrix, therefore its largest eigenvalue is 1. The corresponding eigenvalue of matrix M 1 is: λ M 1 ,1 = (1 -r d )(1 -rs)(1 + r div ). The proposition is proved by

observing that λ M 1 ,1 ≤ 1 ⇔ rs ≥ 1 - 1 (1 -r d )(1 + r div )
and applying Theorem 3 (note that M 1 is positive regular: this is not the case for matrix M).

D Expected size of the GC derived from the multi-type GW process (Section 3.2)

Proposition Let i be the initial state, z 0 := |i| its 1-norm (|i| := N +2 j=0 i j ). The expected size of the GC at time t:

N k=0 (iM t ) k = |i| ((1 -r d )(1 + r div )(1 -rs)) t
Proof For the sake of simplicity, let us suppose that the process starts from a single B-cell belonging to the affinity class a 0 = i with respect to the target trait. We do not need to specify the transition probability matrix used to define the mutational model allowed.

We recall the expression of M t obtained by iteration:

M t =     M t
(A + B) t = t j=0 C j t A t-j B j ( 19
)
For all j, 0 ≤ j ≤ t:

A t-j B j = 2 t-j (1 -r d ) t-j r t-j div (1 -rs) t-j (1 -r d ) j (1 -r div ) j (1 -rs) j Q t-j N = (1 -r d ) t (1 -rs) t (2r div ) t-j (1 -r div ) j Q t-j N Hence: C j t (2r div ) t-j (1 -r div ) j = (1 -r d ) t (1 -rs) t (2r div + 1 -r div ) t = (1 -r d ) t (1 -rs) t (1 + r div ) t , as stated by Equation (3) for z 0 = 1. This result can be easily generalized to the case of z 0 ≥ 1 initial B-cells.

(A + B) t = (1 -r d ) t (

E Proof of Proposition 5

Proposition Let us suppose that at time t = 0 there is a single B-cell entering the GC belonging to the i th -affinity class with respect to the target cell. Moreover, let us suppose that Q N = RΛ N L. For all t ≥ 1, the expected number of selected B-cells at time t, is:

E(St) = rs(1 -rs) t-1 (1 -r d ) t N =0 (2λ r div + 1 -r div ) t as k=0 r i l k ,
Proof Let us suppose, for the sake of simplicity, that Q N is diagonalizable:

Q N = RΛ N L , (20)
We can prove by iteration that:

M t =     M t 1 t-1 k=0 M k 1 M 2 0 2×(N +1) I 2     (21)
It follows from ( 20) and ( 21) that for all t ≥ 1, M t can be written as:

M t =     RD t L R t-1 k=0 D k L M 2 0 2×(N +1) I 2     , ( 22
)
where D = 2(1 -r d )r div (1 -rs)Λ N + (1 -r d )(1 -r div )(1 -rs)I N +1 is a diagonal matrix. We obtain its expression thanks to Proposition 2.

Moreover, by Proposition 3 and Equation [START_REF] Maclennan | B-cell memory and the persistence of antibody responses[END_REF] we have:

M = R DL M 2 0 2×(N +1) I 2 , ( 23
)
where From Equations ( 22) and ( 23):

M t-1 M =     RD t-1 DL RD t-1 L M 2 + R t-2 k=0 D k L M 2 0 2×(N +1) I 2    
Since, by hypothesis, i = (0, . . . , 0, 1, 0, . . . , 0, 0), with the only 1 being at position i, 0 ≤ i ≤ N , then iM t-1 M denotes the i th -row of matrix M t-1 M. Therefore, we are interested in the sum between 0 and as of the elements of the i th -row of matrix M t-1 M, i.e. of the i th -row of matrix RD t-1 DL, since clearly as ≤ N . D t-1 D is a diagonal matrix whose th -diagonal element is given by:

D t-1 D = (2(1 -r d )r div (1 -rs)λ + (1 -r d )(1 -r div )(1 -rs)) t-1 •(2(1 -r d )r div λ + (1 -r d )(1 -r div )) = (1 -rs) t-1 (1 -r d ) t (2λ r div + 1 -r div ) t
The result follows observing that:

RD t-1 DL ik = N =0 D t-1 D r i l k .

Figure 1 :

 1 Figure 1: Schematic representation of model described by Definition 3. Here we denote by aff s := N -a s , the fitness of each B-cell in the affinity class whose index is a s (see Definitions 1 and 2).

Figure 2 :

 2 Figure 2: Numerical estimation of the extinction probability η of the GC with respect to r s for r d = 0.1 and r div = 0.9.

Remark 3

 3 Under certain hypotheses about the mutational model and the GC evolution, one could justify the claim of Proposition 6 by heuristic arguments, without considering the (N +3)-type GW process. This leads to approximately estimate the expected number of selected B-cells at time t (Appendix F).

Figure 4 (

 4 Figure 4 (a) shows the peak of positive selected B-cells at generation t for a certain set of parameters.

Figure 3 :

 3 Figure 3: (a) Dependence of the expected average affinity in the GC on a s at time t = 15, for different values of a 0 . The average affinity in the GC is constant with respect to a s . (b) The evolution during time of the expected average affinity in the GC for different values of a 0 . The average affinity converges through N/2, due to the stationary distribution of Q N , the binomial probability distribution.

  (a)), nor by r s or r d . One can intuitively understand this behavior: independently from their fitness, all B-cells submitted to mutation exit the GC. Moreover, r s and r d impact the GC size, but not its average affinity, as selection and death affect all individuals of the GC independently from their fitness.It can be interesting to observe the evolution of the expected average affinity within the GC during time. Numerical simulations of our model show that the expected average affinity in the GC converges through N/2, independently from the affinity of the first naive B-cell (Figure3 (b)). This depends on the mutational model we choose for these simulations. Indeed, providing that the GC is in a situation of explosion, for t big enough the distribution of GC clones within the affinity classes is governed by the stationary distribution of matrix Q N . Since for Q N given by Definition 9 one can prove that the stationary distribution over {0, . . . , N } is the binomial probability distribution[START_REF] Balelli | Random walks on binary strings applied to the somatic hypermutation of b-cells[END_REF], the average affinity within the GC will quickly stabilizes at a value of N/2. Note that for these simulations we chose N = 10, hence affinities are in the range [0, 10].

Figure 4 (

 4 a) evidences this fact. Moreover, as expected, simulations show that the expected size of selected B-cells at a given time t increases with the threshold a s chosen for positive selection (Figure4 (b)). This is a consequence of Proposition 5: a s determines the number of elements of the sum a s k=0 r i l k .

Figure 4 (Remark 5 4 (Definition 10 (

 45410 Figure 4 (c) underlines the correspondence between theoretical results given by Proposition 4 and numerical values obtained by simulating the evolutionary process described by Definition 3. In particular Figure 4 (c) shows the expected (resp. average) number of selected B-cells produced until time t = 15 depending on the threshold chosen for positive selection, a s .

Figure 4 :

 4 Figure 4: (a-b) Expected number of selected B-cells for the time step t = 15for different values of a 0 , depending on r s and a s respectively. There exists an optimal value of r s maximizing the expected number of selected B-cells for a given generation. This value is independent from a 0 and is equal to 1/t as demonstrated in Proposition 6: the red vertical line in (a) corresponds to this value. (c) Comparison between the expected number of selected B-cells until time t given by evaluation of the theoretical formula (Equation (9)), and the empirical value obtained as the mean over 4000 simulations. Vertical bars denotes the corresponding estimated standard deviations. Here N = 7 and r s = 0.3.

Figure 5 :Definition 11 (

 511 Figure 5: Schematic representations of models described (a) by Definitions 10 and (b) by Definitions 11 of exclusively positive (resp. exclusively negative) selection.

Figure 6 :Remark 6

 66 Figure 6: Dependence of greatest eigenvalues of matrices M + (blue circles) and M -(green squares) respectively on a s for N = 10, r div = 0.9, r d = r s = 0.1. Hence (1 -r d )(1 + r div )(1 -r s ) = 1.539 and (1 -r d )(1 + r div ) = 1.71.

Figure 7 :

 7 Figure 7: Dependence of greatest eigenvalues of matrices M + (blue) and M -(green) respectively on r s for N = 10, r div = 0.9, r d = 0.1, a s = 3.

Figure 8 :

 8 Figure 8: (a,b) Dependence of the expected size of the GC after 15 time steps on a s for different values of a 0 . The thick black line corresponds in both figures to the value of the greatest eigenvalue of matrices M + 1 and M -

1

 1

Figure 9 :

 9 Figure 9: Model of positive selection. (a) Expected number of selected B-cells for the time step t = 15 for different values of a 0 , depending on r s . (b) Estimation of the optimal r * s maximizing the expected number of selected B-cells for a given generation, comparing the model of positive selection for different values of a 0 and the model described in Section 3 (we plot the exact value, r s (t) = 1/t, as obtained by Proposition 6). In (b), for simulations corresponding to the model of positive selection we set a s = 5.

  ], M + , r s =0.1 [S 30 ], M + , r s =0.12 [S 30 ], M + , r s =0.15 [GC 30 ], M -, r s =0.1 [GC 30 ], M -, r s =0.12 [GC 30 ], M -, r s =0.[START_REF] Iber | A mathematical model for germinal centre kinetics and affinity maturation[END_REF]

  Section 3 can be applied to define and study other selection mechanisms. For instance in Section 4 we propose two variants of the model analyzed in Section 3, in which selection acts only positively, resp. only negatively. This Section shows how our mathematical environment can be modified to describe different selection mechanisms, which can be studied at least numerically. Moreover, it gives a deeper insight of the previous model of positive and negative selection, by highlighting the effects of each selection mechanism individually, when they are not coupled.

  3 certainly derives from the simplified structure of the model of positive and negative selection, even if this trend is further confirmed in the model of exclusively positive selection. Nevertheless it should be interesting to test the existence of an inverse relation between the selection rate and the timing of GC peak. Selection pressure can be quantified e.g. through comparative analysis between groups of sequences derived from different germline V(D)J segments, as proposed by the statistical framework for BAyesian estimation of Antigen-driven SELectIoN (BASELINe) [44]. BASE-LINe takes into account both mutation targeting bias and substitution bias and identifies point mutations grouped by location. Moreover it adresses the question of positive versus negative selection: positive selection is identified

Figure 11 :

 11 Figure 11: (a) Evolution during time of the expected value of all types for the model of positive and negative selection, with r s varying during time and N = 10. In particular we set r s = 0 until t = 5, r s = 0.1 from t = 6 to t = 15 and r s = 0.3 from t = 16 to t = 30. Z 13 denotes the total size of the GC (i.e. N k=0 Z k ), and we recall that Z 11 corresponds to selected B-cells and Z 12 to dead B-cells. We set r div = 0.3, r d = 0.005 and z 0 = 100 initial naive B-cells. All initial B-cells belong to a 0 = 5, and the selection threshold is a s = 3. (b) Evolution during time of the expected total size of the GC and the selected pool respectively, for the same set of parameters as in Figure 11 (a).

Definition 16

 16 We denote by η the extinction probability of the process (Zt) t∈N : The pgf of Z (z 0 ) t , t ∈ N, which represents the population size of the t th -generation starting from z 0 ≥ 1 seed cells, is F (z 0 ) t = (Ft) z 0 , Ft being the t th -iterate of F (Equation (2)).

  Therefore we can claim that (iM t ) k corresponds to the k th -component of the i th -rowof matrix M t 1 = (2(1 -r d )r div (1 -rs)Q N + (1 -r d )(1 -r div )(1 -rs)I N +1 ) t , where Q N is a stochastic matrix. Matrices A := 2(1 -r d )r div (1 -rs)Q N and B := (1 -r d )(1 -r div )(1rs)I N +1clearly commute, therefore:

D = 2 ( 1 -

 21 r d )r div Λ N + (1 -r d )(1 -r div )I N +1 is a diagonal matrix.

Definition 7 Let Z

		(i) t = (Z	(i) t,0 , . . . , Z	(i) t,N +2 ), t ≥ 0 be a MC where for all 0 ≤
	j ≤ N , Z	(i) t,j describes the number of GC B-cells belonging to the j th -affinity
	class with respect to x, Z	(i) t,N +1 the number of selected B-cells and Z	(i) t,N +2 the
	number of dead B-cells at generation t, when the process is initiated in state

Table 1 :

 1 Parameter choice for simulations in Sections 3.4 (unless stated otherwise).

	N	r s	r d	r div	a 0	a s
	10	0.1 0.1	0.9	3	3

  1) column vector s.t. for all i ∈ {0, . . . , N } m + i = m 2,i2 , m 2,i2 being the i th -component of the second column of matrix M 2 , given in Proposition 2.

  1 -rs) t

	t
	C j t (2r
	j=0

div ) t-j (1 -r div ) j Q t-j N And consequently:

N k=0 (iM t ) k = N k=0 i (A + B) t k = (1 -r d ) t (1 -rs) t t j=0 C j t (2r div ) t-j (1 -r div ) j N k=0 iQ t-j N k

Since Q N is a stochastic matrix, for all n, Q n N is still a stochastic matrix, i.e. the entries of each row of Q n N sum to 1. Therefore:

N k=0 (iM t ) k = (1 -r d ) t (1 -rs) t t j=0

Appendix

A Few reminders of classical results on GW processes

We recall here some classical results about GW processes we employed to derive Proposition 1 (Section 3.1). For further details the reader can refer to [START_REF] Harris | The theory of branching processes[END_REF]. Definition 14 Let X be an integer valued rv, p k := P(X = k) for all k ≥ 0. Its probability generating function (pgf) is given by:

F X is a convex monotonically increasing function over [0, 1], and F X (1) = 1. If p 0 = 0 and p 0 + p 1 < 1 then F is a strictly increasing function.

Definition 15

Given F , the pgf of a rv X, the iterates of F are given by:

(ii) If X and Y are two integer valued independent rvs, then X + Y is still an integer valued rv and its pgf is given by

q ij 2 -For all j 1 < j 2 ∈ {0, . . . , N }:

p (i) (0, . . . , 0, 2

For all i ≤ as:

If i > as then f (i) (s) is the same except for the first line, which becomes:

The values of each m ij are now obtained by evaluating all partial derivatives of f (i) (s) in 1, keeping in mind that for all i ∈ {0, . . . , N },

C Deriving the extinction probability of the GC from the multi-type GW process (Section 3.2)

Let us recall some results about the extinction probability for multi-type GW processes [START_REF] Athreya | Branching processes[END_REF].

Definition 17

Let q (i) be the probability of eventual extinction of the process, when it starts from a single type i cell. As above bold symbols denote vectors i.e. q := (q (0) , . . . , q (N +2) ) ≥

F Heuristic proof of Proposition 6

Proposition For all t ∈ N the value rs(t) which maximizes the expected number of selected B-cells at the t th maturation cycle is:

Hypothesis 1 Q N converges through its stationary distribution, denoted by m = (m i ), i ∈ {0, . . . , N }.

Hypothesis 2 Zt explodes, where (Zt) t∈N is given by Definition 4.

Let Zt, t ≥ 0 be the random variable describing the GC-population size at time t before the selection mechanism is performed for this generation. For the sake of simplicity, let us

It follows:

Conditioning to Zt = k, Z t+1 is distributed as the sum of k independent copies of Z 1 , which gives:

Thanks to Hypotheses 1 and 2, if t is big enough, there is approximately a proportion of m i elements in the i th -affinity class with respect to x. Therefore, on average at time t there are approximately as i=0 m i E( Zt) B-cells in the GC belonging to an affinity class with index at most equal to as with respect to x, before the selection mechanism is performed for this generation. Each one of these cells can be submitted to selection with probability rs, and in this case it will be positively selected. Hence:

which is maximized at time t ≥ 1 for rs(t) = 1/t.

Remark 8 One observes that the approximation in [START_REF] Murphy | Janeway's immunobiology[END_REF] gives the same value for the optimal rs(t) as in Proposition 6. Nevertheless, it does not allow to describe exactly the behavior of E(St), since it is obtained by approximating the distribution of B-cells in the GC with their stationary distribution.