Marco Molari

Rémi Monasson

Simona Cocco
email: simona.cocco@phys.ens.fr

Survival

Survival probability and size of lineages in antibody affinity maturation

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Affinity Maturation (AM) is a biological process through which our immune system generates potent Antibodies (Ab) against newly-encountered pathogens. AM is also at the base of the efficacy of vaccination, in which this process is artificially elicited through the administration of a dose of Antigen (Ag). The biological mechanisms that govern AM are many and complex, and are the object of many excellent reviews [START_REF] Victora | Germinal centers[END_REF][START_REF] De Silva | Dynamics of b cells in germinal centres[END_REF][START_REF] Bannard | Germinal centers: programmed for affinity maturation and antibody diversification[END_REF][START_REF] Mesin | Germinal center B cell dynamics[END_REF][START_REF] Eisen | Affinity enhancement of antibodies: how low-affinity antibodies produced early in immune responses are followed by high-affinity antibodies later and in memory B-cell responses[END_REF][START_REF] Victora | Clonal and cellular dynamics in germinal centers[END_REF][START_REF] Shlomchik | Linking signaling and selection in the germinal center[END_REF][START_REF] Victora | What are the primary limitations in b-cell affinity maturation, and how much affinity maturation can we drive with vaccination? lessons from the antibody response to hiv-1[END_REF]. Simply speaking, AM works by subjecting a population of B-lymphocytes (or B-cells) to iterative cycles of mutations and selection for Ag binding, which generate a Darwinian evolutionary process that progressively increases their affinity for the Ag. This process is schematically depicted in fig. 1. Maturation takes place in Germinal Centers (GCs), microanatomical structures that appear inside of secondary lymphoid organs. They are divided in two areas: the GC Dark Zone (DZ) in which cells divide and mutate, 1and the Light Zone (LZ) in which they undergo selection. Cells iteratively migrate between these two compartments. Selection in the LZ is completed in two steps. In the first step cells try to bind the Ag, exposed on the surface of Follicular Dendritic Cells (FDCs). In the second step they compete to receive a survival signal from T-follicular helper (Tfh) cells, in the absence of which they undergo apoptosis. Tfh cells are able to probe the amount of Ag captured by B-cells, preferentially delivering this survival signal to the cells that were most successful in capturing Ag. Cells that receive the signal either migrate back to the DZ for additional rounds of mutation and selection, or they can differentiate in Plasma or Memory Cells (PCs/MCs). The former are responsible for the production of Abs to fight the infection, while the latter confer long-lasting protection by remaining quiescent until the same Ag is encountered again, in which case they reactivate and produce Abs or enter GCs for further maturation.

In spite of the many recent experimental advancements in the study of AM, several open questions still remain to be answered, which have important implications in vaccine design. For example understanding the role of Ag availability in controlling maturation might lead to optimization of Ag dosage in vaccines [START_REF] Rhodes | Dose finding for new vaccines: the role for immunostimulation/immunodynamic modelling[END_REF][START_REF] Molari | Quantitative modeling of the effect of antigen dosage on B-cell affinity distributions in maturating germinal centers[END_REF][START_REF] Kang | Affinity inequality among serum antibodies that originate in lymphoid germinal centers[END_REF], or understanding the effect of B-cell precursor frequency and affinity might help improving immunogen design [START_REF] Abbott | Precursor frequency and affinity determine b cell competitive fitness in germinal centers, tested with germline-targeting hiv vaccine immunogens[END_REF][START_REF] Abbott | Factors in b cell competition and immunodominance[END_REF]. Given the complexity of this process, computational models represent an invaluable tool to guide our understanding of AM [START_REF] Chakraborty | A perspective on the role of computational models in immunology[END_REF][START_REF] Buchauer | Calculating germinal centre reactions[END_REF]. In this paper we introduce a stochastic model of AM to study the survival probability of B-cell lineages in GCs. Experimental analysis of vaccine-responsive lineages shows signatures of selection in their reconstructed phylogenies [START_REF] Horns | Signatures of selection in the human antibody repertoire: Selective sweeps, competing subclones, and neutral drift[END_REF]. This selection pressure, which is partially controlled by Ag availability [START_REF] Molari | Quantitative modeling of the effect of antigen dosage on B-cell affinity distributions in maturating germinal centers[END_REF], is important to push lineages towards maturation, but at the same time an excessive pressure might be deleterious. Indeed, several maturation models present a phenomenology termed population bottleneck [START_REF] Zhang | Optimality of mutation and selection in germinal centers[END_REF][START_REF] Wang | Manipulating the selection forces during affinity maturation to generate cross-reactive hiv antibodies[END_REF][START_REF] Wang | Optimal sequential immunization can focus antibody responses against diversity loss and distraction[END_REF], in which strong selection pressure causes a decrease in GC population size, potentially leading to extinction. As a consequence of this trade-off optimal maturation is achieved at intermediate levels of selection pressure. While the bottleneck phenomenology has been studied through numerical simulations so far, we think that it is not incompatible with experimental evidence [START_REF] Firl | Capturing change in clonal composition amongst single mouse germinal centers[END_REF], and might represent a limiting regime for AM as later discussed in section V A.

The study of the survival probability of a population through the maturation bottleneck we report below is both analytical and numerical. We start by considering the dependence of a lineage survival probability on the progenitor affinity. Through the use of recursive relations and probability generating functions we are able to evaluate this probability, and also quantify extinction time and progeny size for lineages that go extinct. We then extend our approach to analyze the extinction probability for the full B-cell population, and its dependence on Ag concentration and initial population size. Last of all, we discuss the biological relevance of the bottleneck phenomenology, and how quantifying lineage survival probability might help vaccine design.

FIG. 1. schematic depiction of the germinal center reaction. Inside of a Germinal Center (GC) B-cells undergo iterative cycles of duplication, mutation and selection. Cell duplication and mutation occurs in the GC dark zone, while selection takes place in the light zone. Selection is divided in two steps: cells must first bind the antigen with their B-cell receptors, and then compete to receive a survival signal from T-follicular helper cells. Failure in any of these steps results in apoptosis. Successfully selected cells have some probability of differentiating into Memory or Plasma Cells (MCs/PCs), and exit the GC. Cells that do not differentiate recycle back in the dark zone to start a new maturation cycle.

selection the population increases its average affinity for the Ag over time. In our model each cell in the population is solely characterized by its affinity for the Ag, measured in terms of binding energy and expressed in units of k B T .

The simulation starts when the GC is mature (roughly 1 week after Ag injection [START_REF] De Silva | Dynamics of b cells in germinal centres[END_REF]). The initial population is composed of N i cells whose binding energy is independently extracted from a Gaussian distribution of naive responders, with mean µ i and standard deviation σ i . Cells undergo iterative rounds of duplication, mutation and selection. These steps are schematized in fig. 2.

At the beginning of the round cells duplicate once in the GC DZ. Each daughter cell can then independently either:

• undergo an affinity-affecting mutation with probability p aa , which causes its binding energy to change by some amount ∆ . We assume that ∆ is a random variable, extracted from a Gaussian distribution with mean µ M and standard deviation σ M (see fig.

2 "mutation");

• not mutate or develop silent mutations, with probability p sil . In both cases its affinity is unchanged;

• be hit by a lethal mutation with probability p let , in which case the cell it is removed from the population.

As a result, the total distribution of the changes ∆ is therefore summarized by the kernel

K(∆ ) = p aa 2πσ 2 M exp - (∆ -µ M ) 2 2 σ 2 M + p sil δ(∆ )
(1) where δ(∆ ) is Dirac delta distribution. Notice that, due to lethal mutations, the integral of the kernel K is not normalized to unity but to p aa + p sil = 1 -p let . Parameters are chosen such that only a small fraction of the mutations is beneficial, i.e. decreases the binding energy (cf. appendix A).

After duplication and mutation cells migrate to the LZ where they try to bind the Ag exposed on the surface of FDCs. Failure to do so results in cell death, and only cells that are able to bind the Ag with sufficient affinity survive this step of selection. Similarly to [START_REF] Molari | Quantitative modeling of the effect of antigen dosage on B-cell affinity distributions in maturating germinal centers[END_REF][START_REF] Wang | Manipulating the selection forces during affinity maturation to generate cross-reactive hiv antibodies[END_REF] we consider the survival probability for a cell with binding energy to be given by the following Langmuir isotherm:

P Ag ( ) = Ce - Ce -+ e -Ag , (2)
where Ag is a threshold binding energy and C represents the dimensionless concentration of Ag available for cells to bind. This concentration controls the strength of selection, making successful binding more likely when more Ag is available to bind. In practice it acts by imparting a shift of magnitude log C to the energy threshold. The functional dependence of the selection probability on and C is displayed in fig. 2 ("selection for Ag binding").

In a second selection step cells compete to receive a survival signal from T-follicular helper cells, with the signal being preferentially delivered to cells that bind more Ag. The survival probability for a cell with binding energy is:

P T ( , ¯ ) = Ce - Ce -+ e -¯ , with e -¯ = e - pop (3)
Where the term e - pop represents the average of this quantity over the population and encodes for competition (see fig. 2 "competitive selection for T-cell help"). The surviving cells can then differentiate into plasma or memory cells with total probability p diff . We do not keep track of these differentiated cells in the simulation.

After this step if the population size exceeds the maximum carrying capacity N max cells are randomly removed until this threshold is met. The surviving cells start then the next round of evolution. The values of the model parameters are reported in table I, and discussed in appendix A.

B. Population bottleneck and lineage survival

Similarly to other AM models [START_REF] Zhang | Optimality of mutation and selection in germinal centers[END_REF][START_REF] Wang | Manipulating the selection forces during affinity maturation to generate cross-reactive hiv antibodies[END_REF], for standard parameter values the population initially undergoes a bottleneck state. This is caused by the strong selection pressure initially imposed by Ag-binding selection, which later relaxes if the average population energy reaches values pop ∼ Ag . By controlling the selection pressure (cf. eqs. ( 2) and ( 3)) Ag concentration also impacts the population survival probability.

As an illustration we report in fig. 3 the average evolution of 1000 stochastic simulations for three different values of the concentration C. For all three values the population size initially decreases under the combined effect of the two selection steps (fig. 3A). This decrease lasts for few evolution rounds, and is accompanied by a quick increase in average affinity (fig. 3B). At this point surviving populations are composed of few high-affinity cells, on which the main acting selection force is competitive selection in eq. ( 3). If this selection pressure is not too strong then the population will later expand and mature. Through a mechanism analogous to the one studied in [START_REF] Molari | Quantitative modeling of the effect of antigen dosage on B-cell affinity distributions in maturating germinal centers[END_REF] Ag concentration then controls the maturation speed, as can be seen by comparing the speed of decrease in average binding energy of the population after the bottleneck in the cases C = 10 and C = 3.5 in fig. 3B.

The fraction of surviving simulations as a function of time is shown in fig. 3C. At low concentration (C = 1) the population goes quickly extinct in all simulations. For such small values of Ag concentration competitive selection alone is sufficiently strong to impede population growth. Intermediate concentration values (C = 3.5), on the contrary, are sufficient to sustain population growth. In this case extinction can nevertheless occur close to the bottleneck state, when population size gets transiently small, see fig. 3A; if some cells are able to survive and pass this bottleneck, then the population again grows to full size and continues maturation. Last of all, at high concentration (C = 10), the bottleneck pressure is not sufficient to significantly endanger population survival, and all simulations are able to overcome the lowpopulation state without going extinct, but maturation proceeds more slowly. We will study in detail in the next section the dependence of the survival probability of the population of cells on Ag concentration.

Survival and future expansion is also strongly dependent on the initial distribution of affinities. This effect can be readily observed on lineages originated from a single ancestor, with energy . In fig. 4 we display three examples of lineage evolution in the form of trees in which each node corresponds to a different cell. These lineages differ by the affinity of their progenitor at the root of the tree. The progeny of the lowest-affinity one (red, i ∼ -0.3) goes extinct in few evolution rounds. In the one with intermediate affinity (orange, i ∼ -0.45) only few individuals are able to survive the bottleneck. The high-affinity one (green, i ∼ -1.3) instead expands and eventually takes over the population. To quantify the population survival probability we will first investigate how the survival of single lineages depends on the progenitor affinity.

III. PROBABILITY OF SURVIVAL AND DISTRIBUTION OF EXTINCTION TIMES

In this section we study the probability that a B-cell lineage descending from a single progenitor cell survives through a population bottleneck, in particular how the probability of survival depends on the affinity of the progenitor. We also determine the distribution of extinction times of the lineage. We then make use of these results to evaluate the survival probability for the full population.

A. Case of one lineage

Probability of survival

Let us consider a progenitor cell with binding energy , present in the population at the beginning of evolution t = 0. At each evolution round this cell will divide and its offspring will have some probability of being removed from the population, either due to selection, differentiation or lethal mutations. Since we are interested in studying the bottleneck phase, in which the population is not at its maximum size, we can neglect the enforcement of the carrying capacity constraint. In fig. 5A we report an example of lineage evolution for a progenitor with binding energy = 1. Color indicates the binding energy of each cell, according to the color-scale on top. In this example cells accumulate deleterious mutations until the lineage eventually goes extinct after t = 26 evolution rounds.

We are interested in computing the probability d t ( ) that all of the offspring of a progenitor with binding energy will be extinct by evolution round t, see fig. 5B.

The expression for t = 1 can easily be written as the probability that both daughter cells generated during the duplication phase will be removed by the end of the round. As stated above, this can occur either by lethal mutation, by failing selection or by differentiation. For each FIG. 3. Average evolution of 1000 different stochastic simulations of the model at three different levels of Ag concentration C = 1, 3.5, 10, color-coded according to the legend on the right. A: population size N as a function of evolution round. Shaded area covers one standard deviation for surviving simulations. The minimum population size on the bottleneck depends strongly on Ag concentration B: same as panel A but for the average population binding energy . Notice how for surviving populations the maturation speed depends on Ag concentration. C: Fraction of surviving simulations as a function of time. At low concentration the bottleneck drives all simulations to extinction, while at high concentration the population survives with high probability. FIG. 4. Examples of stochastic lineage evolution through a population bottleneck. We perform a single simulation of our model at Ag concentration C = 5 and consider three different progenitors with different initial affinities (red i ∼ -0.3, orange i ∼ -0.45 and green i ∼ -1.3). We represent their progeny evolution in the form of a tree with each cell corresponding to a node, and encoding affinity in the branch color. The lineage of the red progenitor quickly goes extinct, while the lineage of the orange progenitor survives the bottleneck but only with few individuals. The green progenitor lineage conversely survives the population bottleneck and undergoes great expansion. Notice how fate correlates with the initial progenitor affinity.

daughter cell this probability is more easily expressed as one minus the probability of not being removed:

d 1 ( ) = 1 -d∆ K(∆ ) P S ( + ∆ ) (1 -p diff ) 2 , (4)
where the expression for K(∆) is the one given in eq. ( 1), and P S ( ) is the probability for a cell with binding energy of surviving selection. In the bottleneck state most of the selection pressure is generated by Ag-binding selection (i.e. ¯ t < Ag ). As a first approximation we therefore neglect competitive selection for T-cell help, and consider simply P S ( ) = P Ag ( ) (cf. eq. ( 2)). This introduces two important simplifications. First, the expression of P S does not depend on time. Second, removing the competitive selection decouples the fate of all cells in the population. The probabilities d t ( ) for t > 1 can be evaluated using recursive relations that express the probability of extinction in t rounds as the probability for each daughter cell to either go extinct in one round, or to survive the first round but to have their respective offspring go extinct in the remaining t -1 rounds:

d t ( ) = 1 -d∆ K(∆ ) P S ( + ∆ ) (1 -p diff ) × (1 -d t-1 ( + ∆ )) 2 (5)
In other words, the probability that all of the offspring goes extinct in t rounds is the probability that each of Probability for a progenitor cell with energy that: the lineage goes extinct by evolution round the lineage goes extinct with a total progeny size of cells (not counting the progenitor) a root branch of the lineage goes extinct with cells on the branch Each node in the tree represents a cell, its binding energy encoded using the colorscale on top. In this example cells progressively accumulate deleterious mutations until after 26 evolution rounds the lineage eventually goes extinct. B. Schematic illustration of the quantities analyzed in our theory. On the left we depict a lineage evolution, stemming from a progenitor with binding energy . The probability that such a lineage goes extinct by time t is indicated with dt( ) (in this example t = 5). The quantity qn( ) represents instead the probability that the lineage goes extinct counting a total of n cells not including the progenitor (n = 9 in this example). Finally, bm( ) is the probability that one of the two sub-lineage stemming from the progenitor goes extinct counting m cells (here m = 6 for the left branch and m = 3 for the right one).

the lineages stemming from the two daughter cells generated during the duplication phase of the first round are removed before the end of round t. Since division is symmetric this probability must be the same for each daughter cell, and is the term inside the square brackets. This term is more easily expressed as one minus the probability that the lineage survives (term in the integral). In turn this can be decomposed as the probability that the daughter cell survives mutation (possibly with an energy change of entity ∆ ), selection and differentiation, multiplied by the probability that its offspring does not go extinct by the following t -1 rounds (term on the second line). A more extended explanation for the derivation of this equation and its numerical evaluation is provided in appendix C.

In fig. 6A we plot the behavior of d t ( ) as a function of evolution round t and progenitor binding energy (orange curves, color indicates extinction round t). As expected, the extinction probability is a monotonically increasing function of time and of energy, and reaches an asymptotic value d ∞ ( ) for large t. Our analytical result is in excellent agreement with simulations for the mean extinction probability (blue dots). The asymptotic probability d ∞ ( ) ranges between:

• d ∞ ( → +∞) = 1
. This is easily understood, since high-energy (i.e. low-affinity) cells will not pass the selection step and their progeny will quickly go extinct.

• d ∞ ( → -∞) = min{1, [1 -1/α] 2 }, with α = (1 -p let )(1 -p diff ).
The value of extinction probability d ∞ for very high-affinity cells may be higher than 0 since lethal mutations and differentiation may still drive the lineage to extinction, especially during the first few evolution rounds when the offspring size is still small. The above expression for d ∞ ( → -∞) can be obtained by searching a fixed point to eq. ( 5), and considering that, for → -∞, mutations do not sensibly change the survival probability and can be neglected. The parameter α defined above is then the probability for a high-affinity cell to survive one round and not be removed by lethal mutations or through differentiation. Notice that, if α < 1/2, we have d ∞ ( → +∞) = 1, as is to be expected when on average less than one individual in the offspring will survive. However this case is pathological: in this regime, irrespective of the progenitor energy, the population always goes extinct (d ∞ ( ) = 1).

It is worthwhile to notice that the "infinite-time" extinction probability d ∞ described in our theory represents the probability for a lineage to go extinct in the bottleneck phase, and does not reflect the probability of lineage fixation or extinction when the population is at maximum size. In eqs. ( 4) and ( 5) we indeed neglected the carrying capacity constraint. This is justified when considering the bottleneck state, in which the population is not at its maximum size. However if some lineages are able to survive this state, then the population will eventually grow back to maximum size, and some cells will be randomly removed through enforcement of the carrying capacity. One should therefore extend the theory and add another term to correctly describe the lineage survival probability far away from the bottleneck.

Distribution of extinction times

The probability that a lineage generated by a progenitor with energy goes extinct exactly at round t can A B C FIG. 6. Comparison between stochastic simulations (blue) and theory (orange) for the probability of extinction (A), lineage extinction time (B) and average progeny size at extinction (C) as a function of the progenitor energy i in absence of competitive selection. For each conditions we consider 5000 different stochastic simulations that terminate with extinction at Ag concentration C = 7. A: stochastic extinction probability (blue dots, error bar indicate the standard error of the mean) evaluated as the fraction of simulations that terminate with extinction over the total number of simulations performed. This is compared to the value of dt( ) as described by our theory. B: mean and standard deviation of extinction time (blue) over 5000 simulations terminating in extinction. This is compared to the theoretical prediction (orange) for the mean and standard deviation of this quantity, obtained using the time extinction probability rt( ). C: same as B but for the progeny size. In this case the theoretical predictions are obtained using the generating function theory.

easily be expressed as

r t ( ) = d t ( ) -d t-1 ( ) . (6)
This allows us to evaluate the mean and variance for the extinction time probabilities (see fig. 6B) simply from the first two moments of the distribution:

t = ∞ t=0 t r t ( ), t 2 = ∞ t=0 t 2 r t ( ), (7)
In fig. 6B we compare, in the approximation of no competitive selection, the average extinction time computed from simulations (blue, error bars indicate the standard deviation of extinction times for each progenitor affinity) with theoretical predictions (orange, shaded area covers one standard deviations). We again find a very good match. The average extinction time shows a peak for intermediate affinities, which can be interpreted as follows. Low-affinity progenitors, i.e. having high binding energy have close-to-one probability of extinction, and very often go extinct in the first few rounds. High-affinity cells on the contrary have a small but non-zero probability of extinction, see value of extinction probability in fig. 6A. This is mainly due to affinity-independent terms (differentiation and lethal mutation probabilities) which confer to the lineage a small chance of going extinct during the first few evolution rounds, when the progeny size is still small. For affinities close to = 1 we observe intermediate values for the probability of survival, and maximum value for average extinction time. This behaviour can be better understood when mutations are turned off, in which case equations can be solved exactly, as shown below.

Exactly solvable case of no mutation

We hereafter consider the case of no affinity-affecting mutations, for which the mutation kernel eq. ( 1) reads

K(∆) = (1 -p let ) δ(∆) . (8)
In this case genealogies belong to the class of Galton-Watson trees [START_REF] Watson | On the probability of the extinction of families[END_REF], and the asymptotic survival probability can be derived exactly. This probability is better expressed by considering the quantity

γ( ) = (1 -p let ) P S ( ) (1 -p diff ) (9)
which represents the probability for a daughter cell to remain in the GC and not be removed by either lethal mutations, selection or differentiation. The infinite-time extinction probability d ∞ ( ) can be found by rewriting eq. ( 5) in the limit t → ∞:

d ∞ ( ) = min 1, 1 γ( ) -1 2 . ( 10
)
As expected, lineages will always go extinct if the average number of surviving offspring at division is not greater than one:

d ∞ ( ) = 1 if γ( ) ≤ 1 2 .
In fig. 7A we report the behavior of d ∞ as a function of γ. This function presents a singularity at the critical value γ = 1/2, for which the Galton-Watson process is critical.

Finding an explicit expression for the distribution of extinction times is harder, but results can be obtained for the critical value γ = 1 2 . We find that the extinction time probability behaves asymptotically as a power law, with infinite mean and variance: r t ∼ 4/t 2 for large t. This result, which is a known feature of critical Galton-Watson processes, can be verified by inserting the Ansatz 5), together with the simplified form of the mutation kernel [START_REF] Victora | What are the primary limitations in b-cell affinity maturation, and how much affinity maturation can we drive with vaccination? lessons from the antibody response to hiv-1[END_REF] and the assumption that γ = 1/2. The only admissible solution is α = 4 which, combined with the definition of r t eq. ( 6), proves the statement. In fig. 7B we display the mean and variance of the extinction time distribution as a function of γ. Comparison with fig. 6B shows that the divergence is removed when evolution includes affinity-affecting mutations. Mutations drive lineages away from the critical line, either to high affinities and survival, or to lower affinities and extinction.

d t ∼ 1 -αt -1 + o(t 2 ) in eq. (

B. Case of full population

Building on the results derived above, we now turn to the problem of quantifying the average probability of extinction for the whole population.

As a first approximation we do not consider competitive selection, since most of the selection pressure in a bottleneck is given by Ag-binding selection. Given a total of N i cells in the initial population, having energies { k } k=1...Ni the probability that all cells will be extinct by evolution round t is simply given by the product of extinction probabilities for all cells k d t ( k ). Moreover, since the initial energies are independently extracted from a Gaussian distribution ϕ( ) with mean µ i and standard deviation σ i , the average extinction probability by round t over all possible extractions of the initial population is given by:

P ext (t) = d ϕ( ) d t ( ) Ni (11)
With the help of this formula we evaluate the average survival probability as a function of Ag concentration C and initial population size N i , and compare the prediction with stochastic simulations in which we turn off 9)) in the approximation of no affinity-affecting mutation. In B and C shaded area covers one standard deviation. Notice how extinction times and genealogy sizes diverge at γ = 1 2 .

T-cell selection. The results, reported in fig. 8B andC (blue), match exactly.

In the presence of competitive selection the empirical survival probability evaluated from simulations slightly decreases, compare blue and orange dotted lines in fig. 8B andC. The theory can be extended to account for Tselection in an effective manner. In practice, one needs first to extend the theory to include a time-dependence of the survival probability. At this point competitive selection can be included by introducing an effective coupling between cells in a 'mean field' fashion, by estimating the average evolution of the term ¯ = -log e - pop contained in the expression for the T-selection survival probability eq. ( 3).

Assume that the survival probability P S ( , t) is now time-dependent. The probability of extinction does not depend solely on the number of evolution rounds anymore, but also on the initial time at which the progenitor is considered. We define d t,s ( ) as the probability that a cell, which at the end of round t has binding energy , will have all of its offspring extinct by the end of round s > t. For any value t ≥ 0 we can write as before the probability of extinction in one round:

d t,t+1 ( ) = 1 -d∆ K(∆) P S ( + ∆, t) (1 -p diff ) 2 (12)
And for any pair of rounds s > t ≥ 0, with s -t > 1, the following recursive relation, analogous to eq. ( 5), holds:

d t,s ( ) = 1 -d∆ K(∆) P S ( + ∆, t) (1 -p diff ) × (1 -d t+1,s ( + ∆)) 2 (13)
Finally, similarly to eq. ( 11), the probability that the full population goes extinct by evolution round t is given by:

P ext (t) = d ϕ( ) d 0,t ( ) Ni , (14)
where ϕ( ) is a Gaussian distribution with mean µ i and standard deviation σ i .

At this point we can make explicit the time dependence of the survival probability including selection for T-cell help: P S ( , t) = P Ag ( ) P T ( , ¯ t ) (cf. eq. ( 3)). Given the stochastic nature of our model, the variable ¯ t which quantifies selection pressure is in reality a stochastic variable. We estimate its average evolution using the large-population-size limit described in appendix B, under which the model becomes deterministic. This allows us to numerically evaluate the extinction probabil-ity eq. ( 14). The outcome, however, underestimates the real survival probability (compare red curve and orange dotted line in fig. 8B andC). This mismatch originates mainly from the fact that in the big-size approximation ¯ evolves faster than in stochastic simulations (cf. blue and orange line in fig. 8A-left). In turn, this occurs because the value of ¯ is strongly dependent on the high-affinity tail of the population, whose evolution is influenced by finite-size effects.

This discrepancy can, however, be reduced with a simple finite-size correction. This correction is based on the consideration that the large-size limit of the model approximates the population binding energy histogram with a continuous distribution, encoded in the density function ρ t ( ) (cf. appendix B). At the beginning of evolution this function takes the shape of a normal distribution, corresponding to the initial binding energy distribution of naive responders, with tails extending indefinitely in both directions. As the population is finite in reality, consisting of N i individuals, we do not expect these tails to be populated. The correction procedure consists in removing these tails, by setting the initial distribution equal to zero outside a range delimited by two values

[ -, + ].
These two values are chosen equal to the expected energy of, respectively, the highest and lowest affinity individual in the population. The probability distribution for their binding energies can be expressed as a function of the naive binding energy distribution ϕ( ) (as before a Gaussian with mean µ i and variance σ 2 i ) from which the energy of all cells is extracted. If we call F ( ) = -∞ d ϕ( ) the cumulative distribution function, then these distributions can be expressed as:

ϕ + ( ) = d d [F ( )] Ni (15)
ϕ -( ) = - d d [1 -F ( )] Ni (16)
The values ± simply correspond to the means of these distributions.

Removing the tails to the initial distribution causes an initial slow-down in the evolution of ¯ (cf. green line in fig. 8A-left). This slow-down is eventually lost, but the agreement remains for a time sufficient for most of the stochastic simulations to go extinct (cf. fig. 8A-right) which is the relevant timescale to capture bottleneck survival.

Taking the value of ¯ obtained by combining the bigsize approximation (cf. appendix B) with the cutoff correction described above, and using it to evaluate the population survival probability, we obtain a much better agreement of the theory with simulations (compare orange curve and orange dotted line in fig. 8B andC). The remaining discrepancy are due to the fact that the average evolution of ¯ is still not exactly captured, and the 'mean-field' nature of our approximation, which neglects the feedback of the energies in the population onto ¯ .

IV. LINEAGE SIZE AT EXTINCTION

In this Section we focus on the distribution of sizes of the progeny at extinction. This size strongly depends on the model parameters, such as the energy of the progenitor. An example is displayed in fig. 5A, in which the lineage consists of a total of 82 cells. Like extinction time, this quantity is well-defined only for lineages that go extinct. Populations that are able to pass the bottleneck undergo exponential growth, with a rate that can be calculated from the large-size theory of Appendix B, see [START_REF] Molari | Quantitative modeling of the effect of antigen dosage on B-cell affinity distributions in maturating germinal centers[END_REF].

Recursion equations for the distribution of sizes

Similarly to what was done in the previous Section for the extinction time and probability, we now derive a recursive formula to quantify the total offspring size. We need to keep track of the sum of two random variables representing the numbers of descendants of each daughter cell. The recursion therefore includes a convolution, which is numerically harder to compute but can be handled using probability generating functions. The recursive relations can be expressed in term of these functions, and can be used to evaluate the moments of the probability distribution without having to numerically perform the convolution.

We name q n ( ) the probability that a progenitor with energy generates a total offspring of exactly n cells before extinction, not counting the progenitor itself (see fig. 5B). This probability can be better expressed if we separate the contribution of the two daughter cells. Considering genealogies encoded as binary trees, we call b m ( ) the probability that along the branch corresponding to one of the daughter cells of a progenitor with energy we find a total of m descendants (including the daughter cell itself) before extinction (see fig. 5). The expression for m = 0 is simply given by the probability that the daughter cell is removed before the end of the round:

b 0 ( ) = 1 -d∆ K(∆) P S ( + ∆) (1 -p diff ) = d 1 ( ) (17)
The recursive relation in this case is composed of two equations. The first is a convolution that decomposes the probability of having n descendants as a sum over all possible repartitions of the descendant number along the two branches:

q n ( ) = n m=0 b m ( ) b n-m ( ) (18)
The second expresses the probability to find m descendants along a branch as the probability that the daughter cell survives and has m -1 descendants:

b m ( ) = d∆ K(∆) P S ( + ∆) (1 -p diff ) q m-1 ( + ∆) (19)
We introduce the generating functions Q(z, ) and B(z, ), defined as:

Q(z, ) = ∞ n=0 q n ( ) z n , B(z, ) = ∞ m=0 b m ( ) z m .
(20) In terms of these generating functions equations eqs. ( 18) and [START_REF] Horns | Signatures of selection in the human antibody repertoire: Selective sweeps, competing subclones, and neutral drift[END_REF] become

Q(z, ) = B(z, ) 2 (21) 1 z [B(z, ) -b 0 ( )] = d∆ K(∆) P S ( + ∆) × (1 -p diff ) Q(z, + ∆) (22)
These relations can be used to evaluate the moments of these distributions with two additional considerations. The first is that

∞ n=0 q n ( ) = d ∞ (
). This sum does not converge to one since it only considers lineages that eventually go extinct. For the functions Q and B this translates into:

Q(z = 1, ) = d ∞ ( ), B(z = 1, ) = d ∞ ( ) . (23)
Secondly, the moments of the distributions can be evaluated from the generating functions as:

n k = 1 d ∞ ( ) ∞ n=0 n k q n ( ) = 1 d ∞ ( ) (z∂ z ) k Q(z, )| z=1 (24)
m k = 1 d ∞ ( ) ∞ m=0 m k b m ( ) = 1 d ∞ ( ) (z∂ z ) k B(z, )| z=1 (25)
Applying the operator z∂ z one and two times on eq. ( 21) restitutes the following relations between the first two moments:

n = 2 m , n 2 = 2 m 2 + 2 m 2 (26)
This corresponds simply to the fact that the total number of descendants is the sum of the descendants along the two branches. Applying the same operator on eq. ( 22) gives:

d ∞ ( ) m = d∆ K(∆) P S ( + ∆) (1 -p diff ) × d ∞ ( + ∆) [2 m +∆ + 1] (27) d ∞ ( ) (m -1) 2 = d∆ K(∆) P S ( + ∆) × (1 -p diff ) d ∞ ( + ∆) [2 m 2 +∆ + 2 m 2 +∆ + 1] (28)
These equations can be solved numerically if we express them as fixed-point equations for the functions m and m 2 :

m = 1 d ∞ ( ) d∆ K(∆) P S ( + ∆) (1 -p diff ) × d ∞ ( + ∆) [2 m +∆ + 1] (29)
m 2 = 1 d ∞ ( ) d∆ K(∆) P S ( + ∆) (1 -p diff ) × d ∞ ( + ∆) [2 m 2 +∆ + 2 m 2 +∆ + 4 m +∆ + 1] (30)
The moments for n can then easily be evaluated using eq. ( 26).

In fig. 6C we compare the theoretical prediction for the first two moments (orange line represents the mean and shaded area covers one standard deviation) with the corresponding quantities from stochastic simulations (blue, error bars cover one standard deviation). Once more we find a good match. The peak at intermediate affinities can be explained, as done above for the extinction time, considering the critical nature of this phenomenon at intermediate values of the binding energy. This is done in the next section.

Exactly solvable case of no mutation

Similarly to what done for extinction probability, in the absence of affinity-affecting mutations we can find an explicit expression for the mean and variance of the population size at extinction. By plugging the simplified expression for the mutation kernel eq. ( 8) into eqs. [START_REF] Wang | Manipulating the selection forces during affinity maturation to generate cross-reactive hiv antibodies[END_REF] and ( 22) one obtains the following second degree equation for the generating function B:

z γ( ) B(z, ) 2 -B(z, ) + 1 -γ( ) = 0 (31)
Where as before γ( ) is the probability for a daughter cell not to be removed from the population during the evolution round, cf. eq. ( 9). This equation has two solutions. The correct one can be chosen by considering that B must be a monotonically increasing function of z. This gives:

B(z, ) = 1 -1 -4 z γ( ) (1 -γ( )) 2 z γ( ) (32)
The function Q(z, ) can be evaluated from eq. ( 21), and the mean and variance for the population extinction sizes can be obtained using eq. ( 24). This results in:

n = 2γ( ) 1-2γ( ) if γ( ) < 1/2 2-2γ( ) 2γ( )-1 if γ( ) > 1/2 (33) n 2 -n 2 = n ( n + 1)( n + 2)/2 (34)
This quantity are reported as a function of γ in fig. 7C.

Similarly to what observed for the mean and variance of the extinction time, Both of these quantities diverge for the critical value γ = 1 2 , but this divergence is removed when mutations are considered.

It is interesting to consider the effect of this divergence on the coefficients q n ( ), that represent the probability of a lineage that stems from a progenitor with binding energy to go extinct with a total progeny of n cells (see fig. 5B). From the definition of the generating function Q(z, ) (cf. eq. ( 20)) it follows that the coefficients q n can be obtained by Taylor expansion of this function around z = 0. In turn Q can be easily obtained from eq. ( 32) using the property Q = B 2 (see eq. ( 21)). The expansion results in the following expression for the coefficients when mutations are absent:

q n = (2n + 2)! (n + 1)!(n + 2)! γ n (1 -γ) n+2 n 1 ∼ (1 -γ) 2 √ πn 3 (4γ(1 -γ)) n (35)
In general these probabilities decay exponentially fast as a function of the size n. At the critical value γ = 1/2 however the term 4γ(1 -γ) becomes equal to 1 and the coefficients algebraically decay as n -3/2 . The asymptotic decay of the coefficients q n for n 1 can also be obtained from the behavior of the generating function Q around its singularity at z c = 1/(4γ(1 -γ)). In particular Q ∝ (1 -z/z c ) 1/2 close to the singularity, and therefore q n ∝ n -3/2 z -m c [START_REF] Flajolet | Analytic combinatorics[END_REF], which gives the expected asymptotic behavior described above.

Case of small-effect mutations

Corrections to the asymptotic behavior above arise when small-effect mutations are considered. In particular in eq. ( 8) we substitute the Dirac delta distribution with a peaked Gaussian, and consider a mutation kernel having the form:

K(∆) = (1 -p let ) 1 √ 2πσ exp - ∆ 2 2σ 2
with σ 1 (36) When the standard deviation is small enough one can approximate the integrals in eqs. ( 17) and ( 22) by Taylorexpanding the functions that multiply the mutation kernel around ∆ = 0. This results in the following approximation for eq. ( 22):

B -1 + γ + σ 2 2 γ = zγB 2 + z σ 2 2 γB 2 , ( 37
)
according to the definition of γ (cf. eq. ( 9)), and with inverted commas indicating derivatives with respect to . This equation is analogous to eq. ( 31) with the addition of perturbation terms of the order of σ 2 . We therefore suppose B(z, ) ∼ B 0 (z, ) + σ 2 ∆B(z, ) where B 0 is given by eq. ( 32) and ∆B represents a perturbation to this σ 2 = 0 solution. By plugging this Ansatz in the previous equation we find the following expression for the perturbation:

∆B = (B 0 -1) 2γr 4 2zγ 2 (1 -2γ) + r γ 2 γ + r 2 γ - γ 2 γ , (38)
where r is defined through

r = 1 -4zγ(1 -γ) = z c -z z c . ( 39
)
The critical value of z is given as before by z c = 1/(4γ(1γ)). Moreover, when the survival probability P S is given by eq. ( 2), the derivatives of γ (cf. eq. ( 9)) can be expressed as

γ = -γ(1 -γ) . ( 40
) γ = γ(1 -γ)(1 -2 γ) , (41)
where γ = P Ag ( ). From the expression for ∆B we can derive the perturbation to Q by using eq. ( 21). Keeping only the higher order in σ 2 , we obtain

Q ∼ Q 0 + σ 2 ∆Q , with ∆Q = 2 B 0 ∆B (42)
As stated in the previous section, the behavior of the probabilities q n is strictly related to the behavior of Q around its singularity z c . In particular the magnitude of the perturbation to the coefficients q n introduced by mutations can be derived from the study of ∆Q. If we operate the substitution z = z c (1 -r 2 ) and expand ∆Q in powers of r we obtain:

∆Q = c 4 r -4 + c 3 r -3 + c 2 r -2 + c 1 r -1 + O(1) (43)
with the following expressions for the coefficients:

c 4 = (1 -2γ) 2 (1 -γ) 2 c 3 = -(1 -2γ) 2 (1 -γ) 2 c 2 = -2(1 -γ)(1 -γ)(1 -2γ γ) c 1 = 2(1 -γ)(1 -γ)(3 -2γ -2γ γ) (44)
Based on this expansion we can derive an expression for the perturbation to the coefficients ∆q n = q n -q 0 n caused by weak mutations [START_REF] Flajolet | Analytic combinatorics[END_REF], where q 0 n represents the value in the absence of mutations (cf. section IV 2). We obtain that, for large n,

∆q n = σ 2 [c 4 (n + 1) + c 3 (n/π) 1/2 (2 + 3/4n) + c 2 + c1(πn) -1/2 + O(n -3/2 )](4γ(1 -γ)) n (45)
In fig. 9A we plot the values of the coefficients c as a function of γ. It is interesting to notice that both c 4 and c 3 , controlling the two leading orders, are null for the critical value γ = 1/2. This leaves the next leading order to c 2 which, for this value of γ, is negative. As a result, we expect that the perturbation tends to lower the values of q n at large n for γ ∼ 1/2, and to raise it for γ ≶ 1/2. This means that large-size extinction events are made less probable by the presence of mutations around the critical value, which is consistent with the removal of the divergence observed in fig. 7C when mutations are present.

In fig. 9B,C,D we compare the above prediction for ∆q n with the value provided by numerical simulations for different values of γ, as indicated in each plot. We expect the theoretical prediction to be accurate for large values of n, and as long as the perturbation ∆q n remains small with respect to the unperturbed value q 0 n . However, for increasing values of n the perturbation grows faster than the unperturbed value, eventually invalidating this assumption. As a proxy for an accuracy upper limit we mark with a green dotted line the value of n at which the perturbation ∆q has a magnitude equal to 10% of the unperturbed value q 0 n . Based on eqn. (45) this value scales as σ -4/5 for small σ's.

V. DISCUSSION

In this work we focused on the effects of a bottleneck on a B-cell population in the course of the affinity maturation process. Through a recursive relation that links the probability of bottleneck survival of a cell to the one of its daughter cells we were able to retrieve the dependence of a lineage extinction probability on its progenitor affinity. For lineages that go extinct we also evaluated the mean and variance of extinction time and progeny size, revealing a peak in extinction time corresponding to average affinity progenitors. Lineages stemming from these progenitors spawn in equilibrium between extinction and survival, and persist in this state until mutations drive the lineage either to survival or to extinction. Building on these results we then evaluated the survival probability for the full population as a function of Ag concentration and population size. We also included the effect of competition in an effective manner, using the deterministic model limit combined with a finite-size correction.

The bottleneck phenomenology was included in different maturation models [START_REF] Zhang | Optimality of mutation and selection in germinal centers[END_REF][START_REF] Wang | Manipulating the selection forces during affinity maturation to generate cross-reactive hiv antibodies[END_REF], which considered as optimal the maturation regime in which the B-cell population was subject to a strong enough selection force to grant good affinity enhancement, while at the same time not strong enough to cause population extinction. While the properties of the above models were numerically evaluated using stochastic simulations, we here present exact or approximate derivations for various quantities of interest through the use of recursive relations and probability generating functions. These techniques are often encountered in the context of branching processes, and similar approaches have been used in the study of AM [START_REF] Balelli | Random walks on binary strings applied to the somatic hypermutation of b-cells[END_REF][START_REF] Balelli | Multi-type galtonwatson processes with affinity-dependent selection applied to antibody affinity maturation[END_REF]. In our case we coupled the theoretical analysis with a more realistic model that included both the effect of mutations and selection for Ag binding and competition. We believe our approach provides a better understanding of what controls the lineage survival probability, while at the same time requiring less computa-tional resources when compared to averaging over many stochastic simulations. This allows one to easily explore the effect of changing different model parameters on the survival probability and the lineage size.

In this last section of the paper we discuss the biological relevance of our results, together with some perspectives. We focus on four aspects: the role of population bottleneck as a limiting case in biologically observed maturation, the affinity of the initial B-cell population, how our results could be applied to explain the role of precursor frequency and affinity on the successful colonization of GCs, and finally how the theory might also offer insight in the case of multiple antigens.

A. Population bottleneck in affinity maturation

The population bottleneck phenomenon, intended as a transitory state of low population with a high extinction risk, is a common feature of many maturation models [START_REF] Zhang | Optimality of mutation and selection in germinal centers[END_REF][START_REF] Wang | Manipulating the selection forces during affinity maturation to generate cross-reactive hiv antibodies[END_REF]. However to our knowledge this phenomenon lacks experimental confirmation. This might be due to two main factors: experimental limitations on one hand, and biologically relevant conditions on the other. By nature the experimental observations of germinal centers tends to be destructive: in order to study the cellular composition of a germinal center the animal must often be sacrificed. This allows for the contemporary observation of multiple germinal centers at the same point in time, but not for the repeated observation of the same germinal center at different time points. As a result it is difficult to estimate for example whether chronic germinal center reactions feature long-lived germinal centers, or many short-lived germinal centers that appear and wane [START_REF] Victora | What are the primary limitations in b-cell affinity maturation, and how much affinity maturation can we drive with vaccination? lessons from the antibody response to hiv-1[END_REF]. To our knowledge the only system capable of circumventing these difficulties is the one introduced by Firl and colleagues [START_REF] Firl | Capturing change in clonal composition amongst single mouse germinal centers[END_REF], who developed an intravital imaging system that allows for observation of germinal centers for an extended period of time. Interestingly, they report that in around 40% of their observations the germinal center failed to form, with the population of founder cells initially expanding until around 5 days after antigen administration, and then starting to wane. The authors observe that this might either be due to a competition with a non-imaged clone, or to a failure to establish a productive germinal center, which would be compatible with a "failed bottleneck" scenario. In short, while lacking explicit experimental confirmation, the bottleneck phenomenology is also not ruled out by experimental observations. A second important consideration is that the bottleneck scenario might only represent a particular possible regime in maturation, which might be relatively rare under biologically relevant conditions. To illustrate this we consider the probability for a progenitor with energy progenitor to survive the bottleneck state, as a function A C D FIG. 9. A: values of the power expansion coefficients in eq. ( 43) as functions of the progenitor survival probability γ = PAg. Both c4 and c3 are null at the critical value γ = 1/2. B,C,D: perturbations to the probabilities qn due to weak mutations (cf. eq. ( 36), with σ = 0.05) for three different progenitor affinities, corresponding to γ = 0.4, 0.5, 0.6 (B,C and D respectively). The perturbation ∆qn = qn -q 0 n is evaluated as the difference between the probability of lineage extinction at progeny size n without (q 0 n ) and with (qn) mutations. These were evaluated by performing 10 8 numeric simulations for each of the two conditions. In the plots we report the value ∆qn (4γ(1 -γ)) -n from simulations (blue) and compare it with the theoretical prediction (orange). The value of n at which |∆qn| = 10% |q 0 n | sets an upper limit for the validity of the theory (green).

of the amount of available antigen concentration and of the energy difference progenitor -Ag (cf. fig. 10A). This can be easily evaluated using the t → ∞ limit of recursion eq. ( 5). Here we only consider Ag-binding selection as the major source of selective pressure in the bottleneck and neglect the effect of competition for T-cell help. One can observe a steep transition between a region of almost certain extinction (red) to a region of possible survival (white and blue). The width of the region is essentially controlled by the probability of developing a beneficial mutation, and the magnitude of such mutations. As evident from the form of the selection survival probability eq. ( 2) operating a multiplicative change of antigen concentration by a factor ∆C has the same effect as considering a change of binding energy of magnitude -log ∆C. One could therefore simply include the effect of concentration by rescaling the energy values. As a next step we consider the probability of survival for a full population, as a function of the initial population size N i and the mean µ i and standard deviation σ i of the initial binding energy distribution of cells. In fig. 10C we display the values at which the population survival probability p surv is equal to 50%. These correspond to lines in the µ i / σ i plane, with different colors corresponding to different initial population sizes N i . The lines split the diagram in two regions, of high and low survival probability (p surv ≶ 50%). To visualize the behavior of the system we consider three different cases, of high (case A), inter-mediate (case B) or low (case C) survival probability. In fig. 10D for each of these conditions we display the evolution of the total population size for 100 different stochastic simulations. In the low-survival-probability region, corresponding to case C, in all of the simulations the population quickly goes extinct. Conversely, in the high-survival-probability region (case A), the population only undergoes a partial reduction in size, with all of the simulations successfully surpassing the bottleneck. In the intermediate region (case C) we observe a consistent reduction in size for all of the simulations, with both successes and failures in overcoming the bottleneck. This intermediate region corresponds indeed to the bottleneck phenomenology that we are interested in studying in this paper, but biologically might only represent a subset of all the possible maturation regimes. In particular one might expect maturation to normally occur in the high-survival-probability regime in most of the biologically relevant cases, especially for pathogens that are similar to ones already encountered in the past for which high-affinity precursor might already be available in the initial population. The low-survival-probability regime conversely might never occur. In this case in fact no precursor with sufficient affinity is expected to be present in the initial population, and GCs might fail to form altogether. Finally the intermediate regime, in which the bottleneck phenomenology is observed, might be a limit case and only occur for "hard" antigens, when only few precursor of acceptable affinity exist in the initial population. These limit cases might be of particular interest when developing novel vaccination strategies, for example against mutable pathogens such as HIV, as will be discussed in section V C.

B. Affinity of the initial population

As shown in the previous section, the bottleneck phenomenology occurs only for values of the parameters that place the system in the boundary between the high and low survival probability region. This is fundamentally controlled by the affinity of the initial population: if enough high-affinity precursors are present in the initial population then the population survival is almost certain, while if all cells in the population have low affinity then survival is rare event. To study this phenomenology in our simulations we chose the size and affinity distribution of the initial population so as to be on this boundary region, while at the same time being compatible with experiments. This was done as follows. The initial population size was set equal N i = 2500. This is in agreement with [START_REF] Eisen | Affinity enhancement of antibodies: how low-affinity antibodies produced early in immune responses are followed by high-affinity antibodies later and in memory B-cell responses[END_REF] which reports around 3000 cells per germinal center. For simplicity we consider the binding energy of each cell to be independently extracted from a normal distribution with mean µ i and standard deviation σ i . This simplification overestimates the heterogeneity of the initial population binding energy, which in reality is usually composed of 50 ∼ 200 different clones [START_REF] Tas | Visualizing antibody affinity maturation in germinal centers[END_REF]. The value of σ i = 1.5 was chosen in agreement with [START_REF] Molari | Quantitative modeling of the effect of antigen dosage on B-cell affinity distributions in maturating germinal centers[END_REF], in which a fit of experimental single-cell affinity measurements of naive responder cells results in an estimated σ i ∼ 1.66. Given these conditions the value of µ i (together with antigen concentration C) controls the probability of surviving the bottleneck (cf. fig. 10C). For most of the simulation we choose a value µ i = 4, so that only around 2% of cells from the initial population have affinity higher than the Ag-binding selection threshold for values of C ∼ 5. This makes so that for standard parameter values the model is close to the phase threshold line (i.e. condition B in fig. 10C andD), which allows us to visualize and study the bottleneck phenomenology. Under these assumptions the initial population will be composed in large part by low-affinity cells. In general, while antigen affinity does indeed drive the recruitment of cells in germinal centers, the presence of cells with undetectable affinity has been observed in different experiments [START_REF] Abbott | Factors in b cell competition and immunodominance[END_REF].

For example in the experiments reported in [START_REF] Eyer | The quantitative assessment of the secreted igg repertoire after recall to evaluate the quality of immunizations[END_REF] only a minor fraction (∼ 30%) of all the antibody-secreting cells elicited by immunization was showing detectable affinity for the administered antigen (K d < 500 nM). Moreover, as experimentally shown in [START_REF] Schwickert | A dynamic T cell-limited checkpoint regulates affinity-dependent B cell entry into the germinal center[END_REF], in the absence of competition even cells with very low affinity for the antigen (K d ∼ 8 µM) can colonize GCs. The role of low-affinity cells in maturation has not yet been experimentally elucidated. Their presence has either been attributed to some form of non-specific activation (e.g. bystander activation [START_REF] Inoue | Generation of memory B cells and their reactivation[END_REF][START_REF] Horns | Memory b cell activation, broad anti-influenza antibodies, and bystander activation revealed by single-cell transcriptomics[END_REF]), or to the specific binding of some unknown "dark antigen" [START_REF] Abbott | Factors in b cell competition and immunodominance[END_REF], or again to specific binding of the administered antigen but with a very low (undetectable) affinity. In either case their number usually decreases during maturation, in favor of higher affinity cells.

C. Effect of affinity and precursor frequency on germinal center colonization

The theory we introduced in this paper offers a framework to quantify the probability for a cell lineage to be able to survive the population bottleneck, successfully colonize a GC and undergo maturation. This has strong connections with the field of vaccine design. In this field the development of new and effective vaccination techniques against pathogens such as HIV requires the stimulation of germline clones that are oftentimes rare or have low affinity for the target. To design an effective vaccine it is indeed important to understand the effect of germline affinity and precursor frequency on the successful recruitment of clones in germinal centers. Recent experiments have shown that at physiological levels of affinity and frequency, successful recruitment depends in a non-trivial way on both of these variables [START_REF] Abbott | Precursor frequency and affinity determine b cell competitive fitness in germinal centers, tested with germline-targeting hiv vaccine immunogens[END_REF][START_REF] Abbott | Factors in b cell competition and immunodominance[END_REF][START_REF] Havenar-Daughton | When designing vaccines, consider the starting material: the human b cell repertoire[END_REF]. In this context, our theory could offer a quantitative prediction for the probability of a clonal family to survive selection, as a function of the initial abundance of cells and of their affinity. To illustrate this we consider a group of N identical cells with initial energy clone . Using our theory we can evaluate the probability that the offspring of at least one cell will survive the population bottleneck. This probability is displayed in fig. 10B. In agreement with what observed in experiments [START_REF] Abbott | Precursor frequency and affinity determine b cell competitive fitness in germinal centers, tested with germline-targeting hiv vaccine immunogens[END_REF][START_REF] Abbott | Factors in b cell competition and immunodominance[END_REF], successful survival depends on both N and clone , with survival probability being increased by both higher affinity and higher precursor frequency. A quantitative understanding of this trade-off might guide decisions in the development of vaccine antigens and improve experiment design.

D. Case of multiple antigens

It would be important to extend the present work to the case of multiple antigens. Contrary to the case of simple antigens, where the response is usually focused on a single dominating epitope, in the case of complex or multiple antigens there might be competition between the binding of different epitopes. ble pathogens such as HIV [START_REF] Victora | What are the primary limitations in b-cell affinity maturation, and how much affinity maturation can we drive with vaccination? lessons from the antibody response to hiv-1[END_REF][START_REF] Chakraborty | Rational design of vaccine targets and strategies for HIV: a crossroad of statistical physics, biology, and medicine[END_REF][START_REF] Robert | Induction of broadly neutralizing antibodies in germinal centre simulations[END_REF]. As a natural extension to our theory one could consider that a single cell might possess different affinities for each antigen mutant. These affinities are potentially correlated, depending on the similarity between the different mutants, and so is the effect of mutations. Taking this into account one could then define mutation and selection probabilities in the presence of multiple Ag mutants, and use a similar approach to the one introduced in this paper to evaluate the lineage and population survival probability. Even without such extension however we believe that the single-antigen framework might still be instructive in the study of multiple antigen maturation. To demonstrate this one can extend the model to the case of two antigens, in a similar fashion to what was done in [START_REF] Wang | Manipulating the selection forces during affinity maturation to generate cross-reactive hiv antibodies[END_REF]. In this extension each cell possesses two binding energies 1 and 2 , encoding its affinity for the two antigens. The val-ues of these energies are initially independently extracted from the same normal distribution of "naive" binding energies. The system then evolves similarly to the singleantigen case, with two differences. The first one is that affinity-affecting mutations change the binding energy of the cell for both antigens with two random contributions ∆ 1 and ∆ 2 , independently extracted from the usual mutation kernel (cf. eq. ( 1)). Notice that as a result of this, mutations that increase the affinity for only one of the two antigens are much more common than mutations that increase the affinity for both. The second difference is in the way selection is performed. Following the "meetall" scenario introduced in [START_REF] Wang | Manipulating the selection forces during affinity maturation to generate cross-reactive hiv antibodies[END_REF] we extend the definition an affinity ceiling and beneficial mutations can accumulate indefinitely. However we believe this approximation to be reasonable since for the bottleneck scenario we consider low-affinity cells that could potentially undergo many affinity-improving mutations. Moreover we consider Ag concentration to be constant, while in reality Ag is subject to natural decay and consumption by B-cells. We believe this approximation to be acceptable when studying bottleneck survival, that in most cases is resolved in few evolution rounds. In our analysis we focused on the evolution of the population of B-cells inside a single GC. In the course of a response however many GCs form inside the body. 3 By averaging over the distribution of energies in the initial population our theory indeed quantifies the average bottleneck survival probability of GCs in this ensemble, in the hypothesis that the affinity of cells each initial population is independent. If instead affinities are correlated or if cells are able to migrate between GCs then the average would be harder to compute. Unfortunately the lack of precise experimental quantification of these processes forbids any meaningful modeling so far. Finally, our model does not account for the fact that GC selection might be highly permissive. Indeed, low affinity cells have been shown to be able to reside in GCs for a extendend periods of time, and permissiveness might especially characterize maturation against complex antigens [START_REF] Molari | Quantitative modeling of the effect of antigen dosage on B-cell affinity distributions in maturating germinal centers[END_REF][START_REF] Finney | Germinal center responses to complex antigens[END_REF]. This effect could be accounted for in our model by adding a small affinity-independent probability of surviving selection. This might have a beneficial impact on maturation, since it could allow low-affinity cells to reside longer in GCs and it might provide them a higher chance of developing beneficial mutations. standard quantity, so that the variable t has no dimension. Similarly, also the binding energy is dimensionless, expressed in standard units of k B T . For simplicity following experiments that indicate a cell-cycle time of 12h or longer [START_REF] Allen | Imaging of germinal center selection events during affinity maturation[END_REF] we consider a single division per round.

Other experiments indicate an average of two cell division in the DZ [START_REF] Gitlin | Clonal selection in the germinal centre by regulated proliferation and hypermutation[END_REF]. We point out that, at the expense of simplicity, our theory can also be extended to account for an higher number of cell divisions. In [START_REF] Zhang | Optimality of mutation and selection in germinal centers[END_REF][START_REF] Wang | Optimal sequential immunization can focus antibody responses against diversity loss and distraction[END_REF] mutations occur at a rate of 0.5 per sequence per division, and are silent, lethal or affinity affecting with probabilities of respectively 0.5, 0.3, 0.2. This fixes our effective mutation probabilities to p sil = 0.75, p let = 0.15, p aa = 0.1. To reproduce the fact that most of the mutations are deleterious we pick for simplicity µ M = σ M . This fixes the amount of beneficial mutations to ∼ 16%, which is somewhat higher but still compatible with other models [START_REF] Molari | Quantitative modeling of the effect of antigen dosage on B-cell affinity distributions in maturating germinal centers[END_REF][START_REF] Zhang | Optimality of mutation and selection in germinal centers[END_REF][START_REF] Wang | Manipulating the selection forces during affinity maturation to generate cross-reactive hiv antibodies[END_REF] in which this fraction is set to 5%. Moreover we set µ M = 0.3 so as to set the mean effect of beneficial mutations to ∆ beneficial ∼ -0.15. This value is slightly smaller than ∆ beneficial ∼ -0.53 used in [START_REF] Molari | Quantitative modeling of the effect of antigen dosage on B-cell affinity distributions in maturating germinal centers[END_REF], but this is compensated by the higher rate of beneficial mutations in our model. As described in the main text, the binding energy distribution of the initial population is set to a Gaussian with standard deviation σ i = 1.5, which is compatible with experimental data [START_REF] Molari | Quantitative modeling of the effect of antigen dosage on B-cell affinity distributions in maturating germinal centers[END_REF]. Since evolution of the population is invariant for shifts of the energy space we set Ag = 0. Under this choice of gauge the zero in the energy space is the threshold energy for Ag-binding selection. Moreover we pick µ i = 4 so that the difference µ i -Ag -log C ∼ 2 σ i for the values of Ag concentrations considered in this work (C ∼ 5) and on average only around 2% of cells from the initial population meet this threshold. For simplicity we independently extract the energy of each cell in the initial population from this initial distribution. By doing so we might overestimate the diversity of the initial population. In fact, experiments probing the clonal composition of GCs estimated that early GCs contain around 50 to 200 different clonal families [START_REF] Tas | Visualizing antibody affinity maturation in germinal centers[END_REF]. A more realistic initiation of our GCs would require us to extract the energies of around a hundred different founder cells, and let them duplicate without mutating up to to the full GC capacity. This would generate a less homogeneous initial population than the one we consider in our simulation, but would not otherwise strongly impact our results. Lastly, as in [START_REF] Molari | Quantitative modeling of the effect of antigen dosage on B-cell affinity distributions in maturating germinal centers[END_REF][START_REF] Wang | Manipulating the selection forces during affinity maturation to generate cross-reactive hiv antibodies[END_REF], the probability of differentiation is set to p diff = 0.1.

probability that all of the progeny of the progenitor cell with energy goes extinct by evolution round t is equal to the probability that both of its daughter cell progenies will go extinct before this round. These probabilities can be expressed as one minus their complement: the probability that extinction does not occur by round t.

In turn this is equal to the probability that the daughter cell will mutate, survive selection of the first round, and some of its progeny will also survive the remaining t -1 rounds. Expressing this latter survival probability in terms of its complement completes the recursive relation. Notice that while the recursion is back in time (t → t-1), the term d t-1 refers to the survival probability of the daughter cell, and in this sense the recursion is forward in time. In fact t refers in this case to the number of turns left to survive, and decreases of one unit at every jump to the next generation. Numerical evaluation of eq. ( 5) was performed according to the following scheme. The fact that cells can accumulate mutations makes so that the recursion requires evaluating an integral over all the possible values of the energy. As a consequence in order to evaluate d t ( ) one needs to know the value of d t-1 ( ) at the previous time interval for any value of . In practice, when evaluating the function at a given final time T one can define a matrix d[t, i], where the first index t = 1, . . . , T runs over all possible values of times, and the index i = 0, . . . , I operates a discretization of the energy space 0 , . . . , I , in such a way that d[t, i] = d t ( i ). We indicate with δ = i+1i the value of the discretization step. The first row of the matrix is given simply by the discretized form of eq. ( 4) in the main text:

d[1, i] =   1 - j δ K[j] P S [i + j] (1 -p diff )   2 (C2)
Where K[i] = K(i δ ) and P S [i] = P S ( i ) are the discretized forms of the mutation kernel and selection survival probability. From here rows of the matrix can be iteratively populated using the discretized form of eq. ( 5):

d[t, i] =   1 - j δ K[j] P S [i + j] (1 -p diff ) (1 -d[t -1, i + j]))   2 (C3)
In practice the evaluation of the extinction probability at time T requires of the order of T convolutions between arrays having a size ( I -0 )/δ . Many efficient algorithms exist to perform fast numerical convolution. In our implementation for example, written in Python version 3.9, the evaluation of d t for t = 150 and discretizing the binding energy interval [-50, 50] at a discretization step δ = 0.01 requires around 0.15 seconds on a laptop. The more general case of a time-varying survival probability, described in eqs. ( 12) and ( 13) in the main text, also presents a similar complexity. In this case one is in-terested in evaluating the probability d 0,T ( ) that a cell having energy at time 0 will have its progeny extinct by evolution round T . In order to do so one defines a matrix d[t, i] = d t,T ( i ), with t = 0, . . . , T -1. The row of the matrix for t = T -1 can be populated from the discretized form of eq. ( 12) in the main text. At this point using the discretized form of eq. ( 13) one can populate row t -1 from row t, until reaching the desired value t = 0. As in the previous case, this requires of the order of T convolutions of arrays having the size of the discretized binding energy space, and requires a similar computational time.

FIG. 2 .

 2 FIG.2. schematic description of the processes that make up a simulated evolution round in our model. At the beginning of the round cells duplicate once. Each cell can then independently develop a mutation. Cells that undergo a lethal mutation (p let = 15%) are removed from the population, while cells that develop an affinity-affecting mutation (paa = 10%) receive an additive change ∆ to their binding energy, extracted from the displayed Gaussian distribution. Notice that most of the mutations have a deleterious effect on affinity. Cells undergo then selection for antigen binding and compete to receive T-cell help. For each selection step we display the functional behavior of the probability of surviving as a function of the progenitor affinity and antigen concentration C. In selection for T-cell help competition is obtained by making the threshold energy ¯ depend on the current affinity distribution of the population, as displayed. Cells that are able to survive selection have a probability p diff = 10% of differentiating into memory or plasma cells, and exiting the cycle. Last of all, if the population size exceeds the threshold value Nmax = 2500 cells in excess are randomly removed. The remaining cells will then begin a new evolution round.

FIG. 5 .

 5 FIG.5. A. Example of lineage issued from a progenitor with binding energy i = 1 obtained from a stochastic simulation performed at Ag concentration C = 7 in the approximation of only Ag-binding selection. Each node in the tree represents a cell, its binding energy encoded using the colorscale on top. In this example cells progressively accumulate deleterious mutations until after 26 evolution rounds the lineage eventually goes extinct. B. Schematic illustration of the quantities analyzed in our theory. On the left we depict a lineage evolution, stemming from a progenitor with binding energy . The probability that such a lineage goes extinct by time t is indicated with dt( ) (in this example t = 5). The quantity qn( ) represents instead the probability that the lineage goes extinct counting a total of n cells not including the progenitor (n = 9 in this example). Finally, bm( ) is the probability that one of the two sub-lineage stemming from the progenitor goes extinct counting m cells (here m = 6 for the left branch and m = 3 for the right one).

FIG. 7 .

 7 FIG.7. Value of the extinction probability d∞ (A), average extinction time t (B) and average progeny size n (C) as a function of the survival probability γ( ) (cf. eq. (9)) in the approximation of no affinity-affecting mutation. In B and C shaded area covers one standard deviation. Notice how extinction times and genealogy sizes diverge at γ = 1 2 .

FIG. 8 .

 8 FIG.8. Probability of population survival in a bottleneck condition as a function of initial population size Ni and Ag concentration C. A: Left: comparison between the evolution of ¯ in stochastic simulations (blue, mean and standard deviation over 5000 simulations) and theoretical prediction without (red) and with (orange) finite-size correction. This correction consists in cutting the tail of the initial energy distribution in proximity of the expected value for the highest-affinity individual. The correction improves the prediction for the evolution of ¯ at short times. Right: fraction of surviving simulations as a function of evolution round. With the finite-size correction the value of ¯ is well-approximated during the time it takes for most of the simulations to go extinct. In this example we set C = 7, Ni = 100. B: bottleneck survival probability as a function of antigen concentration C. Comparison between stochastic simulations (dotted line, error bars indicate the standard error of the mean) and the predictions our theory (full lines). Stochastic simulations are reported both without (blue) and with (orange) competitive selection for T-cell help (T-sel). For the theory instead we consider the case without T-sel (blue), with T-sel (red) and with T-sel plus finite-size correction (orange). In the absence of T-sel all cells evolve independently, and the theory and simulations match exactly. The inclusion of T-sel slightly decreases the survival probability in stochastic simulations. Accounting for this contribution by using the infinite-size estimate for the evolution of ¯ overestimates the selection pressure. Adding the finite-size correction results in a much better estimate. In this example we set Ni = 100. C: same as B, but the survival probability is evaluated as a function of the initial size of the population Ni. Here we set C = 7.

FIG. 10 .

 10 FIG. 10. A: bottleneck survival probability of a lineage stemming from a single progenitor cell with energy progenitor, as a function of antigen concentration. B: probability of bottleneck survival for an initial population composed of N identical cells with binding energy clone . C: population survival probability as a function of its initial size Ni and the mean µi and standard deviation σi of the initial binding energy distribution. Colored lines indicate the region of the phase diagram in which the probability of survival is equal to 50%. Different colors distinguish different initial population sizes according to the colorbar on the right. Each line divides the plane in two regions, corresponding to high (psurv > 50%) and low (psurv < 50%) survival probability. We consider three conditions that belong to either one of these two regions (condition A and C), or to the separation line (condition B) to be visualized in the next panel. D: evolution of total population size N for stochastic model simulations realized under three different conditions as indicated in panel C. For each condition we display 100 different stochastic simulation trajectories, with time t measured in evolution rounds. The simulation were performed at an initial population size Ni = 2500 and antigen concentration C = 10, while the mean and standard deviation (µi, σi) for the initial binding energy distribution were varied according to the condition (A: (2.3,1.2), B: (4.3, 1.), C: (5.3, 0.6)). Results displayed in panels A,B and C do not take into account competitive selection, while in panel D competitive selection is included in simulations.

  FIG.11. Joint probability distribution Pt( 1, 2) of the binding energies ( 1, 2) of the B-cell population in the two-antigens case described in the text. The four different plots correspond to four different times (t = 0, 10, 50, 100), measured in evolution rounds. The probability distribution was estimated by averaging over 500 different stochastic simulations. For each plot we also display the marginal probabilities Pt( 1) and Pt( 2). Simulations were performed at a value of antigen concentration C1 = C2 = 5, and considering an initial average binding energy µi = 3 for both antigens.

  parameter

d 2 = [ 1 - 2 = 1 - 2 = 1 -

 212121 t ( ) = P (progeny extinct by round t) = [P (daughter cell's progeny extinct by round t)] P (daughter cell's progeny not extinct by round t)] d∆ P (mutation ∆ ) P (survive round) P (progeny survives next t -1 rounds) d∆ P (mutation ∆ ) P (survive round) (1 -d t-1 ( + ∆ ))

TABLE I .

 I standard values of model parameters. Unless otherwise specified these are the values used in simulations. The choice of their value is discussed in appendix A .

	value	description
	µi, σi 4, 1.5	mean and variance of the population initial binding energy distribution
	Ni 2500	initial population size
	Nmax 2500	maximum carrying capacity
	p sil , p let , paa 0.75, 0.15, 0.1	probabilities of silent, lethal, affinity-affecting mutations
	µM, σM 0.3, 0.3	mean and variance of distribution of affinity-affecting mutations
	Ag 0	Ag-binding selection threshold energy
	C see figure captions Ag concentration
	p diff 0.1	differentiation probability

In the DZ B-cells express high levels of Activation-Induced cytidine Deaminase, an enzyme that increases the natural rate of DNA mutations up to 10 -3 per base-pair per generation[START_REF] Kleinstein | Estimating hypermutation rates from clonal tree data[END_REF][START_REF] Mckean | Generation of antibody diversity in the immune response of BALB/c mice to influenza virus hemagglutinin[END_REF][START_REF] Berek | Mutation drift and repertoire shift in the maturation of the immune response[END_REF].

The latter case seems less likely given the presence of low-affinity responders in immunization experiments[29,

36].

While their number is not known with accuracy, it could range from many tens to few hundreds since spleen sections revealed around 20-50 GCs in mice[START_REF] Wittenbrink | Is there a typical germinal center? a large-scale immunohistological study on the cellular composition of germinal centers during the hapten-carrierdriven primary immune response in mice[END_REF].

Notice that from the order of the operators in the evolution round it follows that t must be evaluated using eq. (B4) not directly on ρt but rather on S Ag M A ρt

Acknowledgments: We are deeply grateful to Jean Baudry, Arup Chakraborty and Klaus Eyer for many useful discussions and interactions.

of the selection survival probabilities as:

, with e -= e -1 + e -2 pop (47)

Where C 1 and C 2 represent the concentrations of the two antigens. As a result of these definitions a cell will survive selection if it is able to bind at least one of the two antigens with good affinity. The combination of these two effects, namely the fact that beneficial mutation will often affect only one of the two components, and selection will allow for survival of cells that have high affinity for at least one of the two antigens, makes so that evolution will split the population in two subsets: cells with high affinity for one antigen or cells with high affinity for the other. To visualize this we simulated the evolution of the system at antigen concentration C 1 = C 2 = 5. We performed 500 different simulations. In fig. 11 we plot the evolution of the average binding energy distribution P t ( 1 , 2 ) of the population, and the two marginal distributions P t ( 1 ), P t ( 2 ) at four different times t = 0, 10, 50, 100. The marginal distributions show how the population divides in these two subsets. These subsets evolve independently, each one improving its affinity for their "target" antigen over time, while the binding energy for the other antigen gradually decreases due to deleterious mutations not being selected against. The evolution of each of these subsets is therefore well captured by the single-antigen framework along the corresponding axis. We point out that this split in the population is a consequence of the two assumptions we made, on the way mutations and selection operate on the two energy component. Evolution can proceed differently if for example the effect of mutations on the two components is strongly correlated, or if selection requires cells to have high affinity for all of the antigens at the same time. 2 If high affinity for both antigens is required then the population will not split in two subsets, but rather evolve along the diagonal of the affinity space 1 = 2 . One could therefore simply study the dynamics along a single component, for which the single antigen case would still offer a useful framework.

E. Model limitations and perspectives

Compared to real AM, our model is obviously simplified in many aspects, for example we do not impose

The values of the parameters are reported in table I, and were chosen based on existing literature.

Mature GCs have a B-cells population consisting of a few thousands cell [START_REF] Tas | Visualizing antibody affinity maturation in germinal centers[END_REF][START_REF] Jacob | In situ studies of the primary immune response to (4-hydroxy-3nitrophenyl) acetyl. III. the kinetics of V region mutation and selection in germinal center B cells[END_REF][START_REF] Mcheyzer-Williams | Antigen-driven B cell differentiation in vivo[END_REF]. We therefore set the initial and maximum size of the population equal to N i = N max = 2500. This is in agreement with [START_REF] Eisen | Affinity enhancement of antibodies: how low-affinity antibodies produced early in immune responses are followed by high-affinity antibodies later and in memory B-cell responses[END_REF] which reports around 3000 cells per GC. However we stress that GCs are heterogeneous in size [START_REF] Wittenbrink | Is there a typical germinal center? a large-scale immunohistological study on the cellular composition of germinal centers during the hapten-carrierdriven primary immune response in mice[END_REF]. Similarly to [START_REF] Molari | Quantitative modeling of the effect of antigen dosage on B-cell affinity distributions in maturating germinal centers[END_REF][START_REF] Wang | Manipulating the selection forces during affinity maturation to generate cross-reactive hiv antibodies[END_REF] we consider the duration of a turn of evolution to be T round ∼ 12h, which is consistent with timing of cell migration [START_REF] Mesin | Germinal center B cell dynamics[END_REF][START_REF] Victora | Germinal center dynamics revealed by multiphoton microscopy with a photoactivatable fluorescent reporter[END_REF]. Time in our model will be rescaled by this Appendix B: Estimation of ¯ evolution Including the effect of competition in our evaluation of the population survival probability requires us to estimate the evolution of ¯ , defined in eq. ( 3). To obtain a tractable approximation, we consider the deterministic limit of big population size. In this limit the population binding energy can be approximated with a continuous distribution, and the state of the system is completely determined by the density function ρ t ( ). This function represents the density of cells having energy at evolution round t, so that its integral is equal to the size of the population, and its normalized version is the population binding energy distribution. Evolution is expressed in terms of operators acting on this function. In particular:

1. Cell duplication corresponds simply to doubling in size:

2. Mutations are represented as the convolution with the mutation kernel K(∆ ) defined in eq. ( 1). Notice that the kernel K is not normalized, to account for the contribution of lethal mutations. It acts on the distribution as:

3. Ag-binding selection is implemented by in the product of the population function with the survival probability eq. ( 2):

4. Similarly, T-cell help selection is given by the product with the survival probability eq. ( 3):

5. Differentiation consists simply in a product involving the differentiation probability:

6. Finally, the carrying capacity constraint is implemented again by a product and is operated only if the size of the population exceeds the maximum limit:

Where again N = d ρ( ) is the current population size.

From these definitions the evolution of the population density function ρ t ( ) can be expressed as:

Combining with the definition for ¯ eq. (B4) provides a way for us to estimate the average evolution of ¯ t . 4 Appendix C: Derivation and numerical evaluation of recursion equation

The core of our theory rests on recursion eq. ( 5), which allows one to evaluate the probability of extinction of the progeny of a cell from the extinction probabilities of its daughter cells. In this appendix we present an intuitive derivation of this equation, and explain how it was numerically evaluated.

As described in the main text, eq. ( 5) operates a recursion on the function d t ( ), defined as the probability that all the progeny of a given progenitor cell with affinity will go extinct by evolution round t. The recursion consists in expressing this probability as a function of the progeny extinction probability for the two daughter cells of the progenitor in the remaining t -1 rounds. These daughter cells might develop a mutation which introduces a change in binding energy of magnitude ∆ , and the recursion will therefore contain the term d t-1 ( + ∆ ). The derivation proceeds in the following manner, schematized in eq. (C1). The