We model the stochastic evolution of the distribution of binding energies of a population of B-cells during the affinity maturation (AM) process. A virtual population of B-cells in the GC is subject to iterative rounds of duplication, mutation and selection, see Figure 1 (Wang et al., 2015). Each B-cell in our model is characterized by the binding energy $ϵ$ between its receptor and the Ag; $ϵ$ is measured in units of $k_{B}T$, where $k_B$ is Boltzmann constant and $T$ the organism temperature (This choice of unit is standard in biophysics, and allows one to simply express Boltzmann factors as $e^{-ϵ}$; in practice, $1k_{B}T\simeq {10}^{-24}$ kcal). This energy is related to the dissociation constant $K_d$ between the B-cell receptor and the Ag through ${\displaystyle ϵ=\log K_d}$ is here expressed in Molar units, other choices of units would shift energies by a constant amount. Hence, lower energies correspond to higher affinities (For example, the dissociation constant ${\displaystyle K_d=}$ 1 μM corresponds to the energy $ϵ=-13.8$, and a tenfold decrease in affinity (with $K_d$ varying from 1 μM to 10 μM) corresponds to an increase in binding energy of 2.3). The main objective of our model is to track the evolution of the distribution of binding energies across the B-cell population, $\rho (ϵ)$, during the GC maturation process. Tracking the full distribution is important for later comparison with experimental data, which themselves consist of affinity distributions. We now describe the main ingredients of the model.

Figure 1

Download asset Open asset

Ag dynamics

In the course of AM, the concentration $C$ of Ag varies over time, due both to gradual release from the adjuvant matrix and to decay and consumption (Figure 2A). At time of injection Ag molecules are trapped in the adjuvant matrix, which constitutes an Ag reservoir. Ag is then quickly released at a fast rate $k^{+}$. Due to recycling of Ag from surface of follicular dendritic cells (FDCs) to endosomal compartments (Heesters et al., 2013; Mesin et al., 2016) available Ag decays at a slow rate $k_{\mathrm{\varnothing }}^{-}$, and are consumed by B-cells at a faster rate, $k_B^{-}{N}_B$, proportional to the number $N_B$ of B-cells. As the amount of Ag is depleted, selection of B-cell is more and more stringent, and the GC eventually dies out.

Figure 2

Download asset Open asset

In this paper, we have investigated the relationship between Ag dosage and quality of immunization outcome. Several studies (Victora and Nussenzweig, 2012; Kang et al., 2015; Eisen and Siskind, 1964; Goidl et al., 1968; Nussenzweig and Benacerraf, 1967) report the fact that better affinity maturation is not always favored by higher doses of Ag, but can instead be enhanced by lower doses. Similarly, the strength of a response to a vaccine, usually measured through the count of responding cells, may show a bell-like curve at intermediate dosages, and understanding the mechanisms underlying this behavior and locating the optimal Ag dose are of crucial importance (Rhodes et al., 2019). Our works provide quantitative theoretical and experimental support to these findings. In particular, the stochastic model for affinity maturation we consider here is capable of explaining and accounting for the existence of an intermediate optimal Ag dosage, that results in the highest average affinity of the recalled population. While our model is inspired by previous studies of the evolution of a population of B-cells in a Germinal Center during Affinity Maturation, such as (Wang et al., 2015, it differs in two substantial ways.

First, our model is amenable to detailed mathematical analysis. We show that the stochastic evolution of the distribution of binding energies can be accurately approximated by a deterministic dynamics (see Figure 6), which we resolve exactly. Under constant Ag concentration, the distribution of binding energies behaves as a traveling wave, whose speed and growth rate can be recovered by solving an appropriate eigenvalue Equation 4. The dependence of these two quantities on Ag concentration reveals the role Ag availability plays in controlling the strength of selection, both in the generated data and models. In particular, high Ag dosage results in low selection pressure and no maturation, and conversely too low Ag dosage in high selection pressure and population extinction. Only intermediate Ag concentration and intermediate selection pressure ensures both population survival and successful AM.

Second, we show that a single set of parameters of our model is able to reproduce quantitatively the many distributions of single-cell affinities measured on IgG-SC extracted from mice immunized against TT corresponding to multiple protocols largely varying in Ag dosages and delays between injections. To determine the best parameters, we introduce a maximum-likelihood-based inference method. Our inference method fully exploits the results of the experimental technique, developed in Eyer et al., 2017, giving access not to the average affinity, as titer measurement would, but to the complete affinity distribution of the recalled Ab-SC population. This population information is crucial for accurate inference of the model parameters and for a meaningful validation of the model. Furthermore, the inferred parameters provide insights on the internal processes of affinity maturation, such as on the role of permissiveness, as discussed later. Inference techniques are powerful instruments in this respect, since they help us investigate experimentally unaccessible features of the system through their indirect but measurable effects. Our inference procedure is very flexible and can readily be applied to new datasets, providing ad-hoc estimates of parameters for different antigens or even different organisms.

Our stochastic model for affinity maturation is subject to homogenizing selection (Figure 3), to degrees depending on the presence of a high-affinity precursor (Figure 3D), in agreement with experimental evidence (Tas et al., 2016; Abbott et al., 2018). In addition, the initial choice of founder clones accounts for a large part of the stochasticity in the maturation outcome (Appendix 1—figure 3E to H), Hence, in our model, selective expansion of high-affinity precursors plays an important role in affinity enhancement. Affinity enhancement is also obtained through the accumulation of beneficial mutations. When observing the distribution of beneficial and deleterious mutations in the MC and PC populations (Appendix 1—figure 10), one finds that, even though on average cells accumulate very few mutations during the AM process, selection tends to favor the fixation of beneficial mutations and the disappearance of deleterious ones.

Both the mathematical analysis and the inference procedure are made possible by the fact that our stochastic maturation model is well-approximated by its deterministic counterpart. This is usually not the case when describing the evolution of fitness waves (Neher and Walczak, 2018). In many systems, stochastic fluctuations may play a major role, for example when the evolving population passes through a bottleneck, and transiently has very low size, before increasing again. Fluctuations may also be acquire crucial importance when the evolution lasts so long that the leading edge of the fitness wave has time to exponentially amplify and govern the bulk of the population. Here, experimentally measured quantities, such as the distribution of affinities, are the outcome of an average over multiple GC reactions in the spleen. While single simulated GCs show signs of individuality, see homogenizing selection and the evolution of clonality in Figure 3, the average product of multiple GCs is well-approximated by our deterministic theory. Moreover, stochastic effects are also partially mitigated by the fact that we consider quantities related to the integral over time of the fitness wave evolution, namely, the MC and PC distributions. Hence temporal fluctuations are smoothed out. Another factor contributing to this mitigation is the limited selection we infer. The permissiveness of selection results in a less drastic decrease of the population size, and a reduced sensitivity to fluctuations from the leading edge of the fitness wave.

The role of permissiveness in germinal center selection is still an open question (Bannard and Cyster, 2017; Mesin et al., 2016; Victora and Mouquet, 2018). Through phenomena such as bystander activation (Horns et al., 2019) and stochastic noise, GC selection may also allow intermediate- and low-affinity clones to survive, rather than maturing exclusively via selection of the few best clones (Lau and Brink, 2020; Tas et al., 2016). These phenomena generate a wider diversity than previously appreciated, especially when considering complex Ags displaying different epitopes (Kuraoka et al., 2016). In Finney et al., 2018 for example the authors try to characterize the GC response to complex Ags such as influenza vaccine, as opposed to simple ones such as haptens. While in the latter case a strong homogenizing selection and affinity maturation is observed, for complex Ags response is more polyclonal and a consistent part of the GC population (20–30%) is composed of low-affinity clones. This suggests a more permissive nature of the GC selection, in which even low-affinity clones have a non-zero probability of receiving T-cell help.

To model these effects, we have introduced two parameters, $a$ and $b$, in the competitive selection process involving survival signals from $T$-helper cells, see Equation 2 and Appendix 1—figure 1D. $a$ corresponds to the baseline probability for cells to survive a selection step, while $b$ is equal to the probability for cells to fail selection even if they have high affinity; this could be due for example to the limited availability of T-cell help, which could increase the stochasticity of the selection process (Krishna and Bachman, 2018). The role of the parameters $a$ and $b$ in controlling the population evolution is studied in section appendix (Permissive and stochastic selection: effect of a, b parameters Our maximum likelihood fit of the data yields $a=0.12$ and $b=0.66$. These values imply that the probability that a high-affinity cell to survive the second step of selection is $1-b=34\%$, about two and a half times the probability for a low-affinity cell, given by $a=12\%$. This observation is in support for the permissive and stochastic nature of selection, at least in our experimental conditions. The non-permissive variant of our model with base-line levels $a=b=0$ in Equation 2, referred to as variant (B), offers a much worse fit of the data, even when taking into account the smaller number of parameters of this variant (see appendix 'Possible model variations').

Our inference procedure supports the statistical prevalence of variant (C), with T-cell-based selection only, with respect to (A), which included Ag-binding selection. The fact that Ag-binding selection does not seem to be a limiting step for GC colonization, at least in the range of our experimentally measurable affinities, is compatible with experiments performed in Schwickert et al., 2011, in which it is shown that in absence of high affinity competitors even clones with low affinity (as low as ${\displaystyle K_d\sim 8\mu \mathrm{M}}$ or equivalently $ϵ\sim -11.7$) can colonize GCs. This is also in accordance with the fact that selection in GCs should be relatively permissive (Bannard and Cyster, 2017; Victora and Nussenzweig, 2012) in order not to limit the diversity of the repertoire. Let us emphasize, however, that the difference in the BIC of the two selection models is rather weak and that our conclusion is contingent on the data set collected and analyzed here.

Our experimental setup does not allow us to identify whether the IgG-SCs we observe originate from reactivated MCs or residual PCs generated during previous immunizations. We therefore compared the experimental measurements with a weighted mixture of the MC and PC populations predicted by our model. This mixture, which we call the Ab-SC population, represents the population of cells that respond to antigenic challenge under particular conditions. We introduced the parameters $g_{\text{recall}}$ and $g_{\text{imm}}$, corresponding to the fraction of reactivated MCs in the Ab-SC populations when measurement is performed one day after boost or four days after the second injection, and fit their value on the experimental measurements. The result of our inference procedure indicates that, when the system is probed 1 day after pure TT boost, most of the response consists in reactivated memory cells ($g_{\text{recall}}=56\%$). This is in agreement with experimental observations performed in Eyer et al., 2020, in which the frequency of IgG-SCs increased from 0.6 ± 0.1% to 1.6 ± 0.2% one day after the boost, indicating that around 64% of IgG-SCs were not present before the boost. When the measurement is performed 4 days after the second injection then we predict that the vast majority of responders consist of residual PCs ($g_{\text{imm}}=0\%$, with a confidence interval extending to 6%, Appendix 1—figure 6). This is consistent with experimental data (unpublished), which indicate that the majority of IgG-SCs are still active 28 days after CFA immunization, and will be secreting at +4 days.

Concerning the biological difference between the MC and PC populations, it has been observed that MCs show on average less maturation than PCs (Inoue et al., 2018; Shinnakasu et al., 2016; Shinnakasu and Kurosaki, 2017), a feature that is reproduced in our model (Figure 2D) as a consequence of the temporal switch we introduced (Appendix 1—figure 1B) and might be important in maintaining diversity in the response, especially against highly mutable pathogens, and mitigating original antigenic sin (Suan et al., 2017; Morens et al., 2010). The results of our inference are in agreement with the fact that experimentally we observe a higher affinity of the responders if measurement is performed 4 days after the last injection (scheme 1) rather than 1 day after boost (schemes 2,3). This difference in affinity could also originate from some form of selection acting on the responder population during the first days of the response, which could selectively expand high-affinity clones in the time between Ag challenge and measurement. Including this selection in the model would result in a different estimate of the fractions $g_{\text{recall}}$ and $g_{\text{imm}}$. However, for simplicity and lack of explicit experimental evidence we did not include this selection in the model.

In building the model, we chose to only keep the minimal features that could allow us to understand the existence of an optimal dosage and be able to reproduce experimental observations, while still being mathematically tractable. Among the simplifications, the number of duplications per cell is considered independent of the cell affinity. It has been, however, shown that an affinity correlates with GC dark zone dwelling time and number of divisions (Gitlin et al., 2015). This phenomenon introduces an effective fitness difference, which is in practice qualitatively accounted for by the selection terms in our model. Moreover, we consider the distribution of affinity-affecting mutations $K_{\text{aa}}(\mathrm{\Delta }ϵ)$ to be independent of the clone’s affinity, similarly to Wang et al., 2015; Wang, 2017; Zhang and Shakhnovich, 2010. In reality, independence holds only away from affinity peaks in the Ab sequence space; close to these peaks, affinity-increasing mutations become rare, and it is expected that Ag affinity of clones eventually saturate, while the binding energy can take arbitrarily low values in our model. However, in the regime defined by the values of the parameters inferred on our experimental data, MCs and PCs generated by our stochastic model accumulate on average very few mutations in the course of evolution (appendix, Quantifying beneficial and deleterious mutation events and Appendix 1—figure 10), with the maximum number of beneficial mutations accumulated being compatible with experimental evidence (Gérard et al., 2020) ( appendix Quantifying beneficial and deleterious mutation events). In this regime, mutations account for only a part of the maturation, the rest being achieved through selection of high-affinity founder clones ( appendix Quantifying beneficial and deleterious mutation events). This is in line with the limited maturation observed in our experiments. In cases where the saturation effect may become relevant, other approaches to model the effect of affinity-affecting mutations might be more appropriate, for example the introduction of a ‘shape space’ representation (Shaffer et al., 2016; Wang et al., 2016). The model and results reported here do not include Ab-feedback (Wang et al., 2015), the phenomenon by which GC B-cells not only have to compete amongst themselves for Ag acquisition but also with Abs produced earlier in maturation (Bannard and Cyster, 2017; Mesin et al., 2016; Zhang et al., 2013), which could prevent B-cells from internalizing Ag by binding to it. We did not include Ab-feedback in our model, however preliminary investigations (not shown) suggest that it would not affect the existence of an optimal dosage range. GC lifetimes reported in literature vary considerably, from 1 to 2 weeks for soluble protein boosting to several months or longer for certain infections (Victora and Mouquet, 2018; Mesin et al., 2016). In alum immunizations GC lifetimes of 3–4 weeks have been observed (Takahashi et al., 1998). In our simulations, a long lifetime for GCs is observed and for a high dose of Ag they can have an effective lifetime lasting up to 3 months (Appendix 1—figure 2E). The concentration of Ag is crucial in determining the strength of selection and the lifetime of the GC in our model. In reality, Ag dosage value also controls the initiation of the GC and AM. In particular, one could expect that for very low dosages the GC reaction would not be initiated at all. For simplicity, we avoid including this phenomenon in our model, and GC reaction takes place in our simulations even at very low Ag dosages, with the result that very few, highly affine MC are produced in this regime. To avoid a discontinuity with respect to the case of null Ag dosage, ${\displaystyle D=0\mu g}$, in which we expect the measured B-cell population to originate directly from naive precursors, we perform differentiation at the beginning of the simulation round, before mutations and selection (see appendix 'Model definition and parameters choice'). This generates a core of low-affinity MCs keeping the average affinity of the population close to ${\mu }_{\text{naive}}$, even when few additional high-affinity MCs are added. However, this might be an unnecessary caution, since when looking at the data we observe that even the lowest tested dosage (${\displaystyle D=0.01\mu \mathrm{g}}$ TT, Figure 5) shows the hallmark of maturation when compared to the the case of zero dosage (${\displaystyle D=0\mu \mathrm{g}}$ TT, Figure 5). This signals that in the dosage range considered in our experiments we expect maturation to occur. Furthermore, in our model Ag inputs, for example resulting from a new injection, cannot enter a GC while the maturation process is ongoing. Our choice is partly justified by the observation that injecting an Ag bolus when a GC maturation process is in place mostly results in disruption of the ongoing GC reaction (Victora and Nussenzweig, 2012; Pulendran et al., 1995; Shokat and Goodnow, 1995; Han et al., 1995; Victora et al., 2010). We only model a single ‘average’ GC, whose output is assumed to be representative of the outcome of AM. In reality, MC and PC populations are generated by many parallel GC reactions, which could in principle weakly interact via invasion of clones from one GC to another (Mesin et al., 2016; Victora and Mouquet, 2018). Last of all, to test the robustness of some of our hypothesis we performed the inference procedure under slightly different conditions. In particuar, we considered the effect of increasing the Ag decay rate, of setting $p_{\text{mem}}$ to be a constant and not depend on the number of MCs accumulated during evolution, and also of considering the MC/PC time-switch to be only partial, with a residual production of MCs all along the evolution. We verified that even in these case the model is in good agreement with the data. The results are reported in appendix sect. 6.

As shown above our model for AM is simple enough to be amenable to detailed mathematical analysis and, yet, is able to accurately reproduce the full affinity distributions of Ab-SCs generated during the immunization process. This finding suggests several extensions to the current work. First our model could be used to predict the outcome of more complex immunization protocols than the ones investigated experimentally in this work. In particular, it would be interesting to consider the case of continuous delivery methods (osmotic pumps, repeated injections…) (Tam et al., 2016; Cirelli et al., 2019), through which the Ag concentration can be precisely controlled over time, and make predictions for the optimal delivery process. Secondly, the quantitative fit of the model parameters was made here possible thanks to the maximum-likelihood algorithm we have introduced, which is flexible and robust. Our inference procedure, whose code is made available with the publication (see Materials and methods code and data availability), could be readily applied to to different measurements, as well as to variants of the present models, with extra parameters corresponding to features of the affinity maturation process that are hardly experimentally accessible, such as selection permissiveness. The combination of quantitative modeling with inference appears as a promising tool to understand the mechanisms governing the immune response and to guide the development of strategies to control and direct it.

The mice used herein were part of the study as published elsewhere (Eyer et al., 2020). In short, BALB-C mice were purchased from Janvier Labs (age 6–8 weeks at start, female) and housed in the animal facilities of Institute Pasteur during experimentation. All immunizations were made intraperitoneal. Each condition was replicated three times in the same cohort ($n=3,N=1$); except when explicitly stated otherwise. From each mouse, between 20’−100’000 cells were assayed in an experimental run.

Spleens were harvested at the indicated time points of the immunization schedule. Spleen cell suspensions were recovered following disassociating using a 40 µm cell strainer. Cellular suspensions were pelleted at 300 g for 5 min, and red blood cell lysis was performed for 1 min using BD Pharm Lyse (BD). Cells were washed twice with MACS buffer and re-suspended in 3 ml of MACS buffer. These cells were further processed according to the manufacturer’s protocol using the Pan B Cell Isolation Kit II (Miltenyi) on a MultiMACS Cell24 Separator Plus (Miltenyi, program depletion). Purity of B-cell lineage was usually above 90% (data not shown).

Ag concentration dynamically changes in parallel with the evolution of the GC. Initially an amount $C_{\text{inj}}$ of Ag is administered through injection (Figure 2A). The value of $C_{\text{inj}}$ determines the initial amount of Ag trapped in the adjuvant matrix, setting the initial value of the Ag reservoir concentration (Equation 5, right). For the sake of comparison with experiments $C_{\text{inj}}$ is proportional to the injected Ag dosage $D$ up to a conversion factor $\alpha$, $C_{\text{inj}}=D/\alpha$, inferred through maximum likelihood fit of the data. The available ($C_{\text{av}}$, appearing in the selection probabilities in Equations 1 and 2) and reservoir ($C_{\text{res}}$) concentrations then evolve as described in the main text (Figure 2A) under the action of release, decay and consumption according to the following equations:

During GC formation (${\displaystyle t<T_{\text{GC formation}}=6\mathrm{d}}$) the number of B-cells, $N_t^B$, appearing in the rate of Ag consumption increases exponentially up to the maximal size $N_i^B=2,500$. More details on the concentration evolution can be found in appendix Model definition and parameters choice.

In the deterministic/infinite size approximation, the stochastic processes that model one round of GC maturation can be written as operators acting on the distribution $\rho$ of binding energies $ϵ$.

• Cell duplication is represented by the amplification operator, consisting in a simple multiplication:

• Mutations are encoded by convolution of the distribution of energies with a mutation kernel $K_{\text{eff}}$ that includes the effect of silent, affinity affecting and lethal mutations (see Equation 26 and appendix sect. 3 'Theoretical solution and eigenvalue equation'):

• Selection for Ag binding (Equation 9) and T-cell help (Equation 10) are encoded simply by a product with the respective probabilities (Equations 1 and 2), where for the latter the probability depends on $\overline{ϵ}$ which in turns depends on the distribution of binding energies, making the operator not linear:

• Finally, the carrying capacity (Equation 11) and the differentiation (Equation 12) processes correspond to multiplications:

The distribution of binding energies at round $t$ then evolves through ${\rho }_{t+1}=\mathbf{𝐄}[{\rho }_t]$, where the complete operator is $\mathbf{𝐄}=\mathbf{𝐃}\mathbf{𝐍}{\mathbf{𝐒}}_{\text{T}}{\mathbf{𝐒}}_{\text{Ag}}\mathbf{𝐑}$, we indicate with $\mathbf{𝐑}=\mathbf{𝐌}\mathbf{𝐀}\mathbf{𝐌}\mathbf{𝐀}$ the operator encoding for two rounds of mutations and amplification. The evolution operator features, in order of application, two rounds of amplification and mutation, Ag-binding selection, T-cell help selection, carrying capacity and differentiation. Notice that for variant (C), there is no Ag-binding selection, and ${\mathbf{𝐒}}_{\text{Ag}}$ is replaced with the identity operator.

The growth rate $\varphi$ and the maturation velocity $u$ shown in Figure 4 are characteristic of the travelling wave nature of the distribution of energies $\rho$ at large ‘times’. When Ag-binding selection is irrelevant (as is the case at large times if $u<0$) and the carrying capacity constraint is omitted, the evolution operator simplifies into

In one round of maturation, we expect the travelling distribution ${\rho }^{*}(ϵ)$ to be shifted by $u$ along the energy axis, and to be multiplied by $e^{\varphi }$. Without loss of generality, we may choose ${\rho }^{*}$ such that $\overline{ϵ}=0$ ; any other choice would merely consists in a translation of ${\rho }^{*}$ along the energy axis. Hence, $\mathbf{𝐄}$ is now a linear operator, and ${\rho }^{*}$ satisfies the eigenvalue Equation 4. In other words, the operator $\mathbf{𝚺}(-u)\cdot \mathbf{𝐄}$, where $\mathbf{𝚺}(-u)$ is the shift operator $ϵ\to ϵ+u$, has for largest eigenvalue $e^{\varphi }$ and associated eigenvector ${\rho }^{*}$. As all the entries of $\mathbf{𝚺}(-u)\cdot \mathbf{𝐄}$ are positive, the Perron-Frobenius theorem ensures that $e^{\varphi }$ is the top eigenvalue associated to the unique eigenvector ${\rho }^{*}$ with all its components positive.

In practice, given a guess value for $u$, one can iterate $\mathbf{𝚺}(-u)\cdot \mathbf{𝐄}$ a sufficient number of times to determine its top eigenvector $v(ϵ;u)$, and compute $\overline{ϵ}(u)$ through

The value of $u$ is then tuned until $\overline{ϵ}(u)=0$. For a graphical representation of the resolution procedure and details on the numerical scheme used, see Appendix 1—figure 4 and appendix Theoretical solution and eigenvalue equation.

Nine parameters of the model have been obtained through maximum likelihood fit of the data:

• the conversion factor $\alpha$, which allows for conversion between experimental administered Ag dosage $D$, measured in micrograms, and the dimensionless administered Ag concentration of our model, $C=D/\alpha$.

• the Ag consumption rate per B-cell $k_B^{-}$, which controls the GC lifetime and also the extent of the affinity maturation.

• the mean ${\mu }_{\text{naive}}$ and variance ${\sigma }_{\text{naive}}^2$ of the Gaussian binding energy distribution for the GC seeder clones, elicited directly from the naive population.

• the binding energy threshold ${ϵ}_{\text{Ag}}$ for a B-cell to be able to bind Ag with sufficient affinity to internalize it (cf Equation 1). This parameter does not appear in variant (C), where selection is mediated by T-helper cells only.

• The T-cell selection characteristic coefficients, $a$ and $b$, encoding respectively the baseline probabilities to survive or not survive selection, see Equation 2 and Appendix 1—figure 1.D.

• The weight parameters $g_{\text{recall}}$, $g_{\text{imm}}$, representing the MC fraction in the measured population of IgG-SCs for the two protocols, respectively for schemes 2 and 3 with measurement one day after boost, and scheme 1 with measurement 4 days after second injection.

We use a procedure that maximizes the average likelihood of experimental affinity measurements. To do so, we perform the following steps:

1. For each of the 15 different experimental conditions $𝒮$ (five different dosages in scheme 1, plus seven different dosages in scheme 2, plus four different injection delays in scheme 3, minus the experiment at dosage 10 ${\displaystyle \mu \mathrm{g}}$ TT and 4 weeks delay between injection which is repeated, being present in both schemes 2 and 3) we evaluate the log-likelihood of the experimental measurements through

   • (15) $\ln ℒ(𝒮)=\sum _{s\in 𝒮}\ln {\rho }_{\text{Ab-SC}}({ϵ}_s,𝒮)$

2. where ${\{{ϵ}_s\}}_{s\in 𝒮}$ are the binding energy ($\log K_D$) single-cell measurements performed in condition $𝒮$, and ${\rho }_{\text{Ab-SC}}({ϵ}_s,𝒮)$ is the normalized distribution of binding energies of Ab-SC predicted by the deterministic version of the model, and defined as a weighted sum of the normalized MC and PC distributions with MC fraction $g$: ${\rho }_{\text{Ab-SC}}({ϵ}_s,𝒮)=g{\rho }_{\text{MC}}({ϵ}_s,𝒮)+(1-g){\rho }_{\text{PC}}({ϵ}_s,𝒮)$. This fraction is either $g_{\text{recall}}$ or $g_{\text{imm}}$, depending on the condition $𝒮$ considered (scheme 2,3 or scheme 1). For good comparison with the data the final normalization is done for the part of the distribution inside the experimental sensitivity range $-23.03={ϵ}_{\min }<ϵ<{ϵ}_{\max }=-14.51$. Notice that a measurement equal to ${ϵ}_{\min }$ could in truth originate from any lower value of the energy, a situation not taken into account in the above expression for the log-likelihood. In our dataset, however, only four such measurements are present; they have a very weak influence on the results.

3. Last of all we sum the log-likelihoods over all the conditions $𝒮$ considered for the three different schemes to get the total log-likelihood:

   • (16) $\ln {ℒ}_{tot}=\sum _{𝒮}\ln ℒ(𝒮)$

4. We maximize this global log-likelihood over the space of the nine parameters through the implementation of the parallel tempering algorithm, whose details are specified in appendix section Maximum likelihood fit procedure and Appendix 1—figure 5 and 6. We chose this algorithm because it ensures an effective search of the maximum even in a rugged parameter landscape.

5. Notice that the total log-likelihood $\ln {ℒ}_{tot}$ is sensitive to the number of measurements, which can vary considerably between different conditions. As such when performing the maximization the algorithm favors accuracy over the distributions with the higher number of measurements. Notice also that by evaluating the total likelihood in this manner we neglect the fact that multiple single cell measurements can come from the same stochastic realization of the process and can present some degree of correlation. To validate this inference procedure we generated 10 synthetic datasets using our stochastic model (see Appendix 1—figure 11), with the same number of measurements per scheme as in the experimental dataset. We then inferred, for each synthetic dataset, the values of the parameters, and compared them to their groundtruth. On average all values of the parameters were correctly recovered, see appendix section Validation of inference procedure on artificially generated data and Appendix 1—table 1.

The code containing the implementation of our stochastic and deterministic model is made publicly available in the following repository: https://github.com/mmolari/affinity_maturation; (Molari, 2020; copy archived at https://github.com/elifesciences-publications/affinity_maturation). The repository also includes the experimental dataset, the code to run the inference procedure and the code to reproduce the figures of the main paper (Figures 2–6). Please refer to README.md file for further details.

Mature GCs usually appear 5–7 days after Ag administration. During this time a population of up to hundreds different founder clones colonizes the GC and expands to a total size of a few thousand B-cells. The first mutations in the repertoire are observed around day 6 (Jacob et al., 1993; McHeyzer-Williams et al., 1993). Early GCs are highly polyclonal and contain 50 to 200 clones according to Tas et al., 2016. In agreement with these experimental findings at the time of Ag injection we pick a population of $N_{\text{found}}=100$ founder clones. The affinities of these clones are extracted independently from an initial gaussian distribution whose mean and variance are chosen via the maximum-likelihood procedure described in appendix sect. 4 ('Maximum likelihood fit procedure') and it matches the experimental distribution of germline responders (i.e. splenic IgG-SCs that are observed 1 day after boost of pure Ag, Appendix 1—figure 7 in scheme two and Ag dosage $D=0$). During the time of GC formation the founder clones expand uniformly without mutating. We chose to start our simulation at $T_{\text{GC}}=6$ days after Ag injection. At this point, the GCs are almost fully formed (De Silva and Klein, 2015). The simulation starts with the GC at its maximal size, set to $N_i=N_{\text{max}}=2500$ clones. The maximal size is in agreement with (Eisen, 2014) which reports around 3000 cells per GC, or (Tas et al., 2016) in which GCs are said to contain up to a few thousands B-cells. However, we stress that GCs are heterogeneous in size (Wittenbrink et al., 2011).

From here the model proceeds in evolution rounds. Similarly to Wang et al., 2015 we set the duration of a round to ${\displaystyle T_{\text{turn}}=12\mathrm{h}}$. This number is consistent with timing of cell migration (Victora et al., 2010; Mesin et al., 2016). We neglect the fact that high affinity cells are found to dwell longer in the GC dark zone (Gitlin et al., 2014) undergoing additional divisions. In addition to this the fact that the average cell-cycle time is 12 hr or longer (Allen et al., 2007b) indicates that 12 hr is probably a lower limit for the round duration.

As described in the main text each round consists in cell division with somatic hypermutation, selection for Ag binding, selection for T-cell help and differentiation. In our simulations before starting the first round we perform only once differentiation. This is done in order to recover the good average energy limit at low Ag concentrations. In fact when Ag dosage is small the population quickly goes extinct, while at the same time maturating very fast. Performing differentiation first provides a nucleus of low-affinity germline-like clones whose binding energy controls the average binding energy of the MC population, even if few high-affinity clones are added later. Notice that this does not change the asymptotic behavior of the model, since it would be equivalent to simply changing the order of operations in the round.

Proceeding with the standard turn order then the first operation performed is cell division and somatic hypermutation. During a round we consider cells to divide twice (Mesin et al., 2016). In GC dark zone cells up-regulate their expression of Activation-Induced Cytidine Deaminase. This enzyme increases the DNA mutation rate, inducing mutations in the region coding for the BCR and possibly changing the affinity for the Ag. Mutation rate has been estimated to an average of 10^-3 mutations per base pair per division (McKean et al., 1984; Kleinstein et al., 2003). Similarly to Wang et al., 2015 in which the total binding energy consisted in the sum of contributions from 46 different residues, we consider $N_{\text{res}}=50$ residues to contribute to the binding energy. The probability that upon division at least a mutation occurs in any of the 150 bp coding for these residues can be estimated as $p_{\text{mut}}=1-{(1-{10}^{-3})}^{150}\sim 0.14$. As done in Wang et al., 2015; Zhang and Shakhnovich, 2010; Wang, 2017 at every division and for each daughter cell independently we consider a $p_{\text{sil}}=.5$ probability of developing a silent mutation, in which case the binding energy of the daughter cell remains unchanged, a probability $p_{\text{let}}=.3$ of undergoing a lethal mutation, in which case the cell is removed, and finally a probability $p_{\text{aa}}=.2$ of developing an affinity-affecting mutation. These change the binding energy of the daughter cell by adding a variation $ϵ\to ϵ+\mathrm{\Delta }ϵ$. As done in Wang, 2017 the variation follows a lognormal distribution $K_{\text{aa}}(\mathrm{\Delta }ϵ)$ (Appendix 1—figure 1A)

Appendix 1—figure 1

Download asset Open asset

The parameters of the distribution are chosen so that only 5% of the mutations confer an increase in affinity, while the vast majority causes an affinity decrease. As a result of this process after the two mutations the population size increases almost 4-fold in size (two duplications but some cells are eliminated due to lethal mutations) and the average affinity decreases slightly due to the mainly negative effect of mutations (cf Figure 1 main text, histograms 1 to 2).

After duplication we implement selection. In order to avoid apoptosis cells must bind and internalize a sufficient amount of Ag. The amount of Ag internalized depends both on the affinity of the BCR and on the availability of Ag on the surface of the FDC (Nowosad et al., 2016; Knežević et al., 2018; Batista and Neuberger, 2000). We model this process by expressing the probability of survival of a cell with BCR having binding energy $ϵ$ as described in the main text:

A sketch of this function is shown in Appendix 1—figure 1C. The value of the threshold binding energy has been obtained via maximum likelihood fit of the data, which yields for example ${ϵ}_{\text{Ag}}=-13.59$ for variant A. This selection is not in present variant C of the model.

At the second step of selection, the one leading maturation in our model, B-cells compete to receive a survival signal from T-cells. T-cells in GCs are motile and continuously scan the surface of B-cells, sensing for density of pMHC-II complexes (Shulman et al., 2014). Cells with the highest pMHC-II density receive survival signal preferentially (Depoil et al., 2005; Victora and Nussenzweig, 2012). We again express the probability of survival through

The parameters $a$ and $b$ in Equation 20 represent, respectively, the probability of survival at very high energy and the deficit in probability of survival at very low binding energy. Their effect is better discussed in appendix section (Permissive and stochastic selection: effect of a, b parameters). The formula interpolates smoothly between these two values, as depicted in Appendix 1—figure 1D. The threshold binding energy $\overline{ϵ}$ depends on the population’s binding energy distribution, introducing competition between the cells; the symbol ${⟨\cdot ⟩}_{GC}$ indicates the average of the quantity over the current GC population.

Cells that are able to survive selection can either re-enter the dark zone and start a new round of evolution or differentiate into Ab-producing PCs or quiescent MCs that can be reactivated upon future Ag injection. There is evidence that MC/PC output undergoes a temporal switch: MCs are preferentially produced early in the response (Weisel et al., 2016). Moreover there seems to be an affinity bias in differentiation (Shinnakasu and Kurosaki, 2017). Even though experiments show that affinity plays a role in deciding fate (Shinnakasu et al., 2016) simply by implementing a time-switch in the MC/PC differentiation probability (respectively ${\mu }_{\text{MC}}$, ${\mu }_{\text{PC}}$, Appendix 1—figure 1B) we effectively recover both of these observations. The parameters of these functions (${\tau }_{\text{diff}}=11d$, $\mathrm{\Delta }{\tau }_{\text{diff}}=2d$) are chosen so as to be compatible with Weisel et al., 2016:

Notice that the sum of the two is constant ${\mu }_{\text{MC}}(t)+{\mu }_{\text{PC}}(t)=p_{\text{diff}}=10\%$, compatible with seminal studies (Oprea and Perelson, 1997) that estimated that around 90% of the cells recirculate in the dark zone. In the model we consider for simplicity a complete switch, meaning that for $t\gg {\tau }_{\text{diff}}$ the probability of generating MCs decreases asymptotically to zero. In appendix section (Permissive and stochastic selection: effect of a, b parameters), we discuss the more realistic case of a partial switch, in which there is a residual probability of MC production even for $t\gg {\tau }_{\text{diff}}$.

If new Ag is administered, we consider a new GCR to start. The new GC is colonized partly by new B-cells coming from the naive pool and partly by reactivated MCs (Inoue et al., 2018). We allow only MCs that have already been generated at time of the second injection to colonize the new GC. This is done by picking a set of $N_i=2500$ cells from the naive pool, with binding energies extracted from the same initial Gaussian distribution, and adding to them all the MCs generated up to the time of second injection. The founder clones of the new GC will consist of $N_{found}=100$ cells randomly extracted from this cumulative population. Notice that the probability of extracting a MC from the cumulative population is an increasing function of the number $N_{\text{MC}}$ of MCs extracted at time of injection: $p=N_{\text{MC}}/(N_i+N_{\text{MC}})$. In appendix sect. 12 we discuss instead the case in which the probability of extracting a seeder clone from the memory pool is set to a constant $p=0.3$.

The concentration of Ag evolves as explained in the main text according to the differential Equations 5 and 6 (main text). These equations account for Ag release, decay and consumption. The release rate was evaluated considering a half-life of 17 hr for Ag in CFA (MacLean et al., 2001), which gives a value for the release rate of ${\displaystyle k^{+}=\ln 2/{\tau }_{1/2}\sim 0.98/\mathrm{d}}$. Ag on FDCs can be maintained for a long time, up to a year (Heesters et al., 2014), through a mechanism of endocytosis and recycling of immune complexes (Heesters et al., 2013). To reproduce this long clearance time, we take Ag lifetime to be 8.1 weeks, as measured in popliteal lymph nodes of mice (Tew and Mandel, 1979). This results in a Ag decay rate of ${\displaystyle k_{\varnothing }^{-}=\ln 2/{\tau }_{1/2}\sim 0.012/\mathrm{d}}$. The case of a faster Ag-decay is discussed in appendix sect. 6 ('Possible model variations'). Finally, the consumption rate per B-cell ${\displaystyle k_B^{-}=2.07\times {10}^{-5}/\mathrm{d}}$ (variant C, see main text) is obtained via the maximum likelihood fit procedure described in appendix sect. 4 ('Maximum likelihood fit procedure'). This quantity controls both the GC lifetime and the extent of AM at the end of evolution. For the range of Ag dosages considered simulated GCs have an effective lifetime that vary between 1 or 2 weeks and 3 months (Appendix 1—figure 2E), compatible with lifetimes of real GCs (Victora and Mouquet, 2018). Equations 5 and 6 (main text) are continuous in time. To include them in our discrete timestep model, we perform an update of the values of the reservoir and available concentrations $C_{\text{av}}(t)$, $C_{\text{res}}(t)$ at each round $t=0,1,\mathrm{\dots }$ after selection for T-cell help and before differentiation. The Ag removal rate is given by the cumulative effect of decay and consumption: $k_t^{-}=k_{\mathrm{\varnothing }}^{-}+N_t^B{k}_B^{-}$, and changes at each evolution round due to its dependence on the number of B-cells $N_t^B$ at this stage of the round. The values of the concentrations at the next round $t+1$ are obtained by evolving the corresponding quantities at round $t$ for a time ${\displaystyle T=12\mathrm{h}}$ equivalent to the duration of the round:

Appendix 1—figure 2

Download asset Open asset

At times smaller than the GC formation time ${\displaystyle T_{\text{GC}}=6\mathrm{d}}$ we do not account for GC evolution but we account for the evolution of concentration. This is done as in the previous equations but in this case the total consumption rate is evaluated considering an exponentially increasing number of cells, that starting from one at injection exponentially grows to $N_{\text{max}}$ at the time of GC formation. In particular concentration update at the end of round $t=0,1,\mathrm{\dots },11$ is done considering the following number of B-cells consuming Ag:

In our simulations, GC evolution stops either naturally when Ag depletion leads to population extinction, or when the total simulation time is elapsed and cells are harvested, in which case the simulation is stopped irrespective of the population size and only cells produced up to that point are considered. The total simulation time depends on the immunization scheme considered and is set to match the time elapsed between injection and experimental measurement.

In order to have some insight on the way our stochastic model evolves and to gauge the magnitude of stochastic effects in Appendix 1—figure 2A to D we report the average evolution of 1000 independent GC Reaction (GCR) simulations, performed at an injected Ag dosage of ${\displaystyle D=1\mu \mathrm{g}}$. The average trajectory is reported as a black dashed line and shaded area covers one standard deviation around the mean. Increase of noise in the average binding energy of GC B-cells as a function of time is due to the fact that at each time the average is performed only over the surviving GCs, thus as progressively GCRs end the average becomes more noisy (and also becomes biased to represent only surviving trajectories). We compare the evolution of the stochastic model with the theoretical prediction (orangle lines). As observed in the main text the theoretical prediction performs well at high population size, and loses accuracy for small population sizes (Appendix 1—figure 2A, C). However, since most of the MC and PC population is generated at time of high population expansion the average energy of these two populations is always well estimated (Appendix 1—figure 2B, D). The inaccuracy at small population sizes comes mainly from an overestimation of the selection pressure in our theoretical solution. In fact, the threshold binding energy for T-cell selection $\overline{ϵ}$ (cf appendix Equation 20) is sensitive to the high-affinity tail of the population affinity distribution. As the population size diminishes this tail gets progressively less populated and the value of $\overline{ϵ}$ deviates from its theoretical prediction, evaluated under the infinite-size limit. This slight decrease of selection pressure at the end of the GC lifetime increases slightly the survival time and generates a slight slow-down in maturation. Moreover, in order to estimate the lifetime of GCs in our model we evaluate average and standard deviation of lifetimes of 1000 independent GC simulations for different values of injected Ag dosage. Since in the simulations GCs can also have a more or less long period of small population size prior to extinction (Appendix 1—figure 2A) we also evaluate an ‘effective’ lifetime, considering the GC effectively extinguished when its size reaches 1% of its original size. These lifetimes (Appendix 1—figure 2E) depend on the amount of Ag administered and can vary between few weeks to few months. The same overestimation of the selection pressure at small population sizes leads the theory to slightly underestimate the lifetime of stochastic simulations.

Recent experiments estimated the number of different clones in early GC to be between 50 and 200 (Tas et al., 2016). In our stochastic model, we consider the population of founder clones to be composed of 100 cells. The limited number of founders controls the diversity of the initial population and increases the stochasticity in evolution. However, it does not strongly influence the average outcome. To verify this we compare 1000 stochastic simulations of the standard model (single Ag injection of ${\displaystyle D=1\mu \mathrm{g}}$ of Ag, model scenario C) with a modified version in which the number of founder clones was set equal to the number of cells in the initial population (2500 cells). Results are reported in Appendix 1—figure 3A to D. We observe that limiting the initial population diversity increases stochasticity in evolution, but does not impact much the average evolution trajectory and outcome. This is especially evident when observing MC/PC population evolution (panels B and D). The final average binding energies of these populations are very similar, but the standard deviation around the mean is halved in for the initial population with more founders.

Appendix 1—figure 3

Download asset Open asset

This observation raises the question of how much the outcome of evolution is controlled by the particular initial choice of the founder clones. In Appendix 1—figure 3E to H we quantify this by comparing 1000 stochastic GC evolutions of the standard model (injected Ag dosage ${\displaystyle D=1\mu \mathrm{g}}$, model scenario C), in which the founder population was re-extracted every time, with a modified version in which the founder population was kept the same amongst all stochastic trials. In the latter case, we observe a considerable reduction in stochasticity, indicating that the outcome depends strongly on the initial founder clones choice. This is also in line with the observations made in section 'Maturation induces progressive loss of clonality' and Figure 3D (main text), where we show that the presence of a high-affinity founder clones correlates with a stronger homogenizing selection.

As described in the main text the theoretical solution of the model is obtained by performing the limit of infinite size $N_B\to \mathrm{\infty }$, upon which the evolution of the system becomes deterministic. All the stochastic processes can be implemented via an operator, whose explicit expression is provided in Materials and methods, and which acts on the distribution $\rho (ϵ)$. This function is the product between the binding energy distribution of the population and the population size. In this formalism, the total evolution operator is a combination of the duplication and mutation, selection, normalization and differentiation operators: $\mathbf{𝐄}=\mathbf{𝐃}\mathbf{𝐍}{\mathbf{𝐒}}_{\text{T}}{\mathbf{𝐒}}_{\text{Ag}}\mathbf{𝐑}$, where the operator $\mathbf{𝐑}=\mathbf{𝐌}\mathbf{𝐀}\mathbf{𝐌}\mathbf{𝐀}$ encodes for two rounds of duplication with mutation. The mutation operator $\mathbf{𝐌}$ in particular consists in a convolution of the distribution $\rho$ with a mutation kernel $K_{\text{eff}}(\mathrm{\Delta }ϵ)$ lethal, silent and affinity-affecting nature of the mutations. This kernel is defined as follows:

The first term corresponds to affinity-affecting mutations, whose probability is the product between the probability of developing a mutation and the probability for this mutation to be affinity-affecting $p_{\text{mut}}{p}_{\text{aa}}$. The expression for the affinity-affecting mutation probability $K_{\text{aa}}(\mathrm{\Delta }ϵ)$ is the one in Equation 18. The second term encodes silent mutations, occurring with probability $p_{\text{mut}}{p}_{\text{s}}$, and also absence of mutations, with probability $1-p_{\text{mut}}$. Lethal mutations occurs with probability $p_{\text{mut}}{p}_{\text{l}}$ but since their effect is the removal of cells in the population this contribution is multiplied by zero and is not present in the kernel expression. The contribution of lethal mutation makes so that the normalization of this kernel is not unitary: $\int 𝑑\mathrm{\Delta }ϵK_{\text{eff}}(\mathrm{\Delta }ϵ)=1-p_{\text{mut}}{p}_{\text{l}}$. In the main text, we show how under constant Ag concentration the affinity distribution of the population behaves as a traveling wave whose velocity and growth rate can be found via an eigenvalue equation. This eigenvalue equation reads $\mathbf{𝚺}(-u)\mathbf{𝐄}\rho =e^{\varphi }\rho$, where $\mathbf{𝚺}(\mathrm{\Delta })$ is the translation operator that shifts the distribution of a value $ϵ\to ϵ+\mathrm{\Delta }$ and we consider the restricted evolution operator $\mathbf{𝐄}=\mathbf{𝐃}{\mathbf{𝐒}}_{\text{T}}\mathbf{𝐑}$.

To verify the correctness of our eigenvalue equation theoretical prediction, we solve the eigenvalue problem at a given Ag concentration $C=10$, graphically illustrating the procedure, and show that the theoretical prediction matches the asymptotic behavior of the system. As a first step we set the value of $\overline{ϵ}=0$ in the T-cell help selection operator. This choice constitutes simply a gauge-fixing that removes the translational invariance of our solution, and it also linearizes the evolution operator. Since we do not a priori know the value of the shift $\mathrm{\Delta }$ we solve the eigenvalue equation $\mathbf{𝚺}(-\mathrm{\Delta })\mathbf{𝐄}\rho =e^{\varphi }\rho$ for different values of the shift. In Appendix 1—figure 4 we plot the maximum eigenvalue eigenvectors (distributions in A, color encodes the value of $\mathrm{\Delta }$) and their corresponding log-eigenvalues (C, representing the growth rate) as a function of the shift $\mathrm{\Delta }$. In order for our solution to be consistent it must satisfy the condition $\overline{ϵ}=0$. Therefore we evaluate $\overline{ϵ}$ after performing mutations for all the solutions at varying values of the shift (Appendix 1—figure 4B) and we pick the one (${\mathrm{\Delta }}^{*}$) for which this condition is satisfied as the eigenproblem solution. Upon repeated application of the restricted evolution operator we expect the population to asymptotically grow and shift at the values of ${\varphi }^{*}$ and $u={\mathrm{\Delta }}^{*}$ corresponding to this solution. In Appendix 1—figure 4D, we verify this by plotting the evolution of the normalized affinity distribution of the population, re-shifted on its mean, upon repeated application of the evolution operator, and compare it with the eigenproblem solution (black dashed lines). Moreover, in E and F we compare the instantaneous growth rate and energy shift with the eigenproblem solution predictions. Color in the three plots encodes to the evolution round. We observe that in all three cases the asymptotic prediction is matched.

Appendix 1—figure 4

Download asset Open asset

To numerically solve the eigenvalue equation and obtain the results displayed in Appendix 1—figure 4, operators were implemented as square matrices by discretizing the interval of energy $[10,-20]k_{B}T$ with a discretization step of $0.002k_{B}T$. The maximum eigenvalue and corresponding eigenvector were obtained by repeated application of the restricted evolution and shift operators on an initially normal distribution with zero mean and unit variance. The application was repeated until the L1 distance between the normalized distributions before and after application of the operator was less than 10^-5. This numerical procedure however is computationally expensive, since solving the eigenproblem for a particular value of the Ag concentration $C$ requires the creation and repeated application of large matrices for many values of the shift ${\displaystyle \mathrm{\Delta }}$. A less demanding numerical method to obtain the value of μ and ${\displaystyle \varphi }$ consists in simply simulating the evolution of the population under constant Ag concentration ${\displaystyle C}$, and without carrying capacity constraint and Ag-binding selection, until convergence to the asymptotic travelling-wave state is reached with sufficient precision. This second technique was used to evaluate the values of ${\displaystyle u}$ and ${\displaystyle \varphi }$ displayed in Figure 4B. In this case, the simulation domain was set to the interval $[-100,50]k_{B}T$, with a discretization step of $0.01k_{B}T$. Evolution of the distribution of binding energies was repeated until the following three conditions were met at the same time. First, the L1 distance between the normalized distributions before and after evolution, once they were re-centered around their mean, was less than 5 × 10^-5. Second, the relative change of the growth rate and wave speed between two rounds of evolution was less than 5 × 10^-5. Finally, as a safety check we also require that at the moment of convergence the mean of the distribution is more than five standard deviations away from the boundaries of the simulation domain.

Nine model parameters $\mathbf{𝐩}$ = (${\mu }_{\text{naive}}$, ${\sigma }_{\text{naive}}$, ${ϵ}_{\text{Ag}}$, $k_B^{-}$, $\alpha$, $a$, $b$, $g_{\text{recall}}$, $g_{\text{imm}}$) were inferred via the maximum likelihood procedure described in Materials and methods. This procedure makes use of the full experimental affinity distribution. Here, we give a more detailed description of the procedure and add considerations on the variation range of the parameters.

The stochastic maximization procedure is based on an implementation of the Parallel Tempering technique. This technique, used in the context of molecular dynamics simulations (Sugita and Okamoto, 1999), consists in simulating different copies of a system at different temperatures, and then allowing the copies to exchange their states with adequate probabilities so that low-energy states are correctly sampled at low temperatures. This is particularly advantageous when the energy landscape is rugged, and low-temperature simulations tend to get stuck in local energy minima, while high-temperature simulation explore the whole space without being able to locate the minima precisely. Allowing for state-exchange between different temperature simulations makes so that high-temperature simulations can help low-temperature ones exit local minima in the energy space and converge to the global optimum, Appendix 1—figure 5A.

Appendix 1—figure 5

Download asset Open asset

Here, we consider the space of all possible values of the parameters as our configuration space, and use the log-likelihood (see Materials and methods for definition) as a proxy for minus the energy. Our algorithm can be summarized as follows. A number $N=10$ of copies of the parameter set is initialized, and at each copy is assigned a simulation temperature $T$ logarithmically evenly spaced between 10³ and 10^-3. The maximization procedure consists in 10,000 rounds of iterative parameter changes and temperature exchanges. For each proposed parameter change, the likelihood difference $\mathrm{\Delta }\log ℒ$ between the new and the original set of parameters is evaluated, and the change is accepted with probability $\min \{\exp \{\mathrm{\Delta }\log ℒ/T\},1\}$, where $T$ is the temperature associated to the parameter set. In the exchange phase the difference in log-likelihood $\mathrm{\Delta }\log ℒ$ and in inverse temperature $\mathrm{\Delta }\beta$ is evaluated for any two parameters sets with consecutive temperatures. An exchange of the two parameters sets is then performed with probability $\min \{\exp \{-\mathrm{\Delta }\beta \mathrm{\Delta }\log ℒ\},1\}$. At the end of the last round the value of the parameters ${\mathbf{𝐩}}_{\text{best}}$ that maximized the likelihood is returned. In Appendix 1—figure 5B we report the trajectories in energy space of all the parameter sets for the inference performed using variant C of selection (see main text), in which eight parameters are inferred (all but the threshold binding energy ${ϵ}_{\text{Ag}}$, which is removed in this variant). The energy for each parameter set is evaluated from the difference in log-likelihood with the best value, $\mathrm{\Delta }E=\log ℒ({\mathbf{𝐩}}_{\text{best}})-\log ℒ(\mathbf{𝐩})$. Notice how trajectories at high temperature explore the space by being able to visit low-likelihood (high-energy) zones of the parameter space. Conversely, low-temperature trajectories gradually converge to the value of the parameters that maximizes the likelihood.

The likelihood-maximization procedure is algorithmically described by the following pseudocode:

Where the initial parameters choice was set to ${\displaystyle {\mathbf{p}}_{\mathbf{i}\mathbf{n}\mathbf{i}\mathbf{t}}=({\mu }_{\text{naive}}=-14.6}$, ${\displaystyle {\sigma }_{\text{naive}}=1.6}$, ${ϵ}_{\text{Ag}}=-13.6$, ${\displaystyle k_B^{-}=2\times {10}^{-5}/\mathrm{d}}$, ${\displaystyle \alpha =0.025\mu \mathrm{g}}$, $a=b=0.2$, $g_{\text{recall}}=g_{\text{imm}}=0.5)$. Single parameter changes are generated as a function of noise level $\eta$:

• For ${\mu }_{\text{naive}}$, ${ϵ}_{\text{Ag}}$ variation is performed by adding a random number extracted with uniform probability in the interval $[-10\eta ,+10\eta ]$.

• For positive parameters ${\sigma }_{\text{naive}}$, $\alpha$, $k_B^{-}$ the variation is performed multiplicatively by introducing a percentage change of the parameter uniformly extracted in the interval $[-\eta ,\eta ]$.

• For the fractions $g_{\text{recall}}$, $g_{\text{imm}}$, $a$ and $b$ variations are performed by adding a random number uniformly extracted in the interval $[-\eta ,+\eta ]$. After the addition some constraints are enforced: the resulting number is constrained in the interval [0, 1], and as a consequence of the definition of the survival probability (Equations 19, 20) also $a+b\le 1$ is imposed.

We set the noise level to depend on the temperature of the system considered, so that higher-temperature simulations also have a higher level of noise, allowing them a faster exploration of the parameter space, while low temperature simulation will perform only a fine-tuning of the parameters. This allows for a more precise convergence. In particular we set $\eta$ to be logarithmically evenly spaced between 10^-2 and 10^-1. Furthermore we propose changes for all the parameters for the four higher-temperature simulations, while changes affect only one randomly chosen parameter for the four lower-temperature simulations.

In Appendix 1—figure 6, we display the evolution of the parameter set that reaches the highest log-likelihood during the maximization procedure. In Appendix 1—figure 6A we report the log-likelihood and the temperature as a function of the round of maximization. Notice how high-temperature states correlate with big fluctuations in the log-likelyhood and in the values of parameters. The peak likelihood value is reached at round 9970, when the trajectory was at the lowest temperature. In Appendix 1—figure 6B we report the evolution of the nine parameters during the maximization rounds. Orange shaded area cover one standard deviation of the posterior distribution around the maximum-likelihood estimate (MLE) of the parameters. This was obtained by approximating the posterior with a Gaussian distribution around the ML value, with a quadratic fit of the log-likelihood variation $\mathrm{\Delta }\log ℒ$ generated by small variations (±5%) of single parameters around the MLE (Appendix 1—figure 6C).

Appendix 1—figure 6

Download asset Open asset

For completeness in Appendix 1—figure 7, we report the comparison between the experimental measurements (orange histograms) and the maximum-likelihood fit of the theoretical model (blue curves) for all the different immunization schemes and all different values of the injected dosage $D$ and time delay between injections $\mathrm{\Delta }T$ used for the fit. All these distributions were used for the likelihood evaluation. The normalization of the distributions is done considering only the part inside the experimental sensitivity range. We find good agreement under all different schemes. Moreover we also report the comparison between the theoretical solution and the stochastic simulations of the model (green histograms, average over 1000 stochastic simulations). This comparison shows that for all regimes considered the theoretical solution is a very good approximation to the stochastic model.

Appendix 1—figure 7

Download asset Open asset

This likelihood maximization procedure is very general and can be easily extended to the inference of any set of parameters in our model, or to other experimental datasets.

Parameters $a$ and $b$ in Equations 19, 20 represent respectively the probability for a B-cell to pass or to fail a selection step irrespective from their affinity. Parameter $a$ can be related to the ‘permissiveness’ of selection, quantifying how likely is for a cell to be positively selected even if its affinity is small. Parameter $b$ conversely encodes the stochasticity in selection, by virtue of which even high affinity cells are not granted survival (e.g. if they don’t manage to encounter enough T-cells). For simplicity we define these parameters as constants, but one could immagine that their value may change over time, for example it may be related to the availability of T-cell help. Here we investigate the effect of these parameters on the asymptotic wave-like behavior of the system. This asymptotic behavior is characterized by two quantities: the asymptotic growth rate $\varphi (C)$ and the asymptotic maturation speed $u(C)$, as defined in the main text Equation 3. These quantities are functions of the Ag concentration $C$. In Appendix 1—figure 8 we report how these functions change when the parameters $a$ and $b$ are progressively increased. In these tests, we consider only the effect of one parameter at a time and we set the value of the other to its standard value $a=0.13$, $b=0.66$. By making the selection more permissive and allowing for survival of even low-affinity cells, parameter $a$ has a double effect on the asymptotic behavior: on the one hand it decreases the asymptotic wave velocity and slows down maturation (Appendix 1—figure 8E), and on the other hand it increases the growth rate of the population (Appendix 1—figure 8C). On the contrary, increasing parameter $b$ corresponds to increasing the chance that high-affinity cells are selected out of the population. This both decreases the growth rate (Appendix 1—figure 8D) and also the maturation speed (Appendix 1—figure 8F).

Appendix 1—figure 8

Download asset Open asset

The selection process in affinity maturation has both a purifying and promoting effect. On the one hand it must negatively select clones that accumulated negative mutations, purifying the population from low-fitness individuals. On the other hand, it must also grant the survival and amplification of the clones that developed affinity-improving mutations. These two properties of selection are weakened by parameters $a$ and $b$ in our model, since they respectively grant survival of low-affinity clones and can cause the removal of high-affinity ones. According to our inference procedure parameters $a$ and $b$ together seem to account for 79% of the selection probability, making so that affinity can make the survival probability vary of only about 21%. This slows down maturation considerably, since it removes any deterrent against accumulating deleterious mutations. The fact that the inference procedure indicated a high value for these parameters suggests that selection in GCs may be permissive, at least for complex Ags, as suggested in Murugan et al., 2018.

In order to test the relative importance of different model parameters we performed the inference procedure using three different variants for selection (see main text). Variant A corresponds to the inference of all nine inferred model parameters (${\mu }_{\text{naive}}$, ${\sigma }_{\text{naive}}$, $k_B^{-}$, $\alpha$, $a$, $b$, $g_{\text{recall}}$, $g_{\text{imm}}$, ${ϵ}_{\text{Ag}}$). Variant B corresponds to the case in which stochasticity and permissivity parameters $a$ and $b$ are set to zero and only the remaining seven parameters are inferred. In variant C instead Ag-binding selection is neglected, and all eight parameters with the exclusion of ${ϵ}_{\text{Ag}}$ are inferred. The resulting maximum likelihood estimate (MLE) of the parameters, along with the corresponding value of the likelihood, is reported for the three cases in Appendix 1—figure 9. The result of the inference procedure show that the removal of Ag-binding selection (variant A vs C, 9 vs 8 parameters) causes only a very modest decrease in log-likelihood, while the removal of stochasticity in T-cell help selection (variant A vs B, 9 vs 7 parameters) generates a consistent log-likelihood decrease. As described in the main text, both the Bayesian Information Criterion (BIC) (Schwarz, 1978) and Akaike Information Criterion (AIC) (Akaike, 1974) are in support of removing the former but not the latter, and accept variant C.

Appendix 1—figure 9

Download asset Open asset

Moreover, to show that the model is robust under minor modifications of the hypotheses we consider three minor variations, keeping variant C of selection, and show that they generate similar MLE predictions for the model parameters. First, we test the effect of introducing a residual asymptotic rate of MC/PC production. This case is labelled soft MC/PC timeswitch. For simplicity in the model we introduced a complete time-switch between MC and PC production in GCs, making so that asymptotically only PCs are produced for $t\gg {\tau }_{\text{diff}}$ (cf appendix Equations 21, 22 and Appendix 1—figure 1B). In this modification instead we introduce a residual rate of MC/PC production ${\mu }_{\text{res}}=10\%$ modifying appendix Equations 21, 22 as:

This makes so that the fraction of MCs in the differentiated population interpolates between ∼90% at small times and ∼10% at big times, granting some residual production of MCs at all times. Applying the inference procedure on this more realistic version of the model results in a better final likelihood than the standard (variant C) version. The inferred values of the parameters are on average similar to the ones obtained with the standard version of the model, with the difference of $g_{\text{recall}}$ and $g_{\text{imm}}$. While the inequality $g_{\text{recall}}>g_{\text{imm}}$ still holds, the MLE for these parameter is higher than in the standard case. This can be expected since these parameters control the fraction of MCs in the elicited Ab-SC population, and in this version of the model the MC population contains differentiated cells that would have belonged to the PC compartment in the standard version.

A further modification involves the fraction of seeder clones extracted from the MC population when colonizing a GC. At the second injection some of the seeder clones for the new GC are extracted from the MC population generated following the first injection. The probability of extracting a seeder clone from the MC pool and not from the initial germline distribution depends in the standard version of the model on the number of accumulated memory cells $N_{\text{mem}}$ as $p_{\text{mem}}=N_{\text{mem}}/(N_{\text{mem}}+N_i)$ (cf main text section Stochastic model for affinity maturation). This should account for the fact that intuitively if more MCs were produced in the previous maturation then also more should be recalled. However one could more simply suppose this probability to be constant. We test this case by setting $p_{\text{mem}}=0.3$. This change generates only a very small likelihood decrease with respect to the standard version of the model, while the MLE for all parameters is almost unchanged.

Finally, we also test the effect of increasing the rate of Ag decay, multiplying it by a factor three (case labelled faster Ag decay). This results in a slight increase in the maximum likelihood, and the values of all the model parameters are again compatible with the one of the standard version of the model, with the exception of the Ag consumption rate $k_B^{-}$ which decreases to compensate the faster decay rate.

To quantify the number and impact of mutation events in our simulations we execute 1000 stochastic simulations of a single GC, at an injected Ag dosage of ${\displaystyle D=1\mu \mathrm{g}}$. In each simulation and for each cell we keep track of the number of beneficial and deleterious mutations that each cell accumulates during the course of evolution on the residues we consider ($N_{\text{res}}=50$, see appendix sect. 1 Model definition and parameters choice). Results are reported in Appendix 1—figure 10. In the top row we display the average number of cells for any value of beneficial and deleterious mutations number and at different times: 10, 30 and 50 days after Ag injection. To have an idea of the population size and maturation state at these timepoints one can refer to Appendix 1—figure 2, in which stochastic simulations are performed under the same conditions. In our simulations after the first days of maturation deleterious mutations start to appear (see Appendix 1—figure 10, t = 10). These are the first mutations to appear since they are much more likely than beneficial mutations (95% vs 5%, see Appendix 1—figure 1A). However during the course of evolution these are gradually removed by selection, until eventually beneficial mutations, despite being much rarer, start to dominate (see Appendix 1—figure 10, t = 50).

Appendix 1—figure 10

Download asset Open asset

In the bottom row of Appendix 1—figure 10 we display in the same way the number of mutations in the MC and PC population. The former is composed of cells that differentiate early ( Appendix 1—figure 1B and Equations 21, 22) and therefore bear less mutations than the PC population. However, in both cases the vast majority of cells harbor very few mutations, with the average number of mutations per cell being 0.13 for MCs and 0.54 for PCs. The accumulation of more than 4–5 beneficial mutations in a single cell is very rare. These numbers are compatible with experiments performed in a recent work (Gérard et al., 2020), in which mice were immunized against Tetanus Toxoid following a protocol similar to the one used in our experiments. The analysis of high-affinity binders showed an average of 6 non-synonymous mutations on the antibody heavy-chain variable region $V_H$ and three mutations in the light-chain variable region $V_L$.

Our model neglects saturation of beneficial mutations, that is the phenomenon by which beneficial mutations cannot be accumulated indefinitely but become rarer as the cell approaches maximum possible affinity. This is partly justified here by the observation that, at least for the inferred value of model parameters, few beneficial mutations are found to accumulate in our simulations. Even when considering clones with the highest number of beneficial mutations, the number of mutations accumulated in our model is compatible with experiments (Gérard et al., 2020). As a final remark, notice that even though MCs are not as strongly skewed towards beneficial mutations as PCs, their average affinity is higher than the one of the starting population (cf Appendix 1—figure 2B). This is because amongst the founder clones selection will expand the ones having higher affinity (cf Figure 3), which will then be overly-represented in the MC population. This shows that in our model maturation is achieved only partially by accumulation of beneficial mutations, the rest being obtained through selective expansion of high-affinity precursors, as also shown in section 'Maturation induces progressive loss of clonality' of the main paper, and confirmed by the strong dependence of the maturation outcome on the initial founder clones population (see Appendix 1—figure 3E to H).

To test the reliability of our inference algorithm we generate 10 artificial datasets using our stochastic model and we apply to them the inference procedure. We then compare the inferred values of the parameters with the real value used to generate the data. To verify that the experimental measurements at our disposal are sufficient to infer the model parameters, when generating the data we took into account the size and composition of the experimental dataset. In particular, we generate each artificial dataset as follows. For every experimentally tested conditions (15 conditions in total, consisting in five tested dosages for scheme 1, seven dosages for scheme 2 and 4 delays for scheme 3, with condition ${\displaystyle D=10\mu \mathrm{g}}$ and $\mathrm{\Delta }T$ = 4 weeks shared between scheme 2 and 3) we run as many stochastic simulations as mice tested for that particular condition. Then from every simulation we extract from the responders population a number of cells equal to the one obtained from each mouse. Extraction of cells is done keeping into account the experimental detection range, therefore we exclude cells having affinity lower that ${\displaystyle K_d=500\mathrm{n}\mathrm{M}}$, and set any affinity higher than the high-affinity threshold ${\displaystyle K_d=0.1\mathrm{n}\mathrm{M}}$ equal to the threshold. Stochastic simulations of each scheme are done in model scenario C, using the standard value of the parameters with one exception: to account for the fact that each mouse contains multiple GCs we raised the number of founder clones in each GC to 2500 instead of 100. This amounts to introducing a greater diversity in the initial population, which in turns reduces the stochasticity in the evolution outcome (see also Appendix 1—figure 3A to D and appendix sect. 2 'Stochastic model analysis') and is similar in spirit to averaging between multiple GCs, as it is the case for the experimental measurements of cells extracted from the spleen of mice. This generation procedure was re-executed 10 times, resulting in 10 different datasets. As an example, in Appendix 1—figure 11B we show the binding energy distribution of the 10 generated datasets for condition ${\displaystyle D=0.5\mu \mathrm{g}}$ of Ag in scheme 2 (histograms in gray). The histograms are close to the prediction of the deterministic model (blue curves), but with some deviations due to stochastic sampling. The average binding energy of the population for all considered condition is reported in Appendix 1—figure 11A. Again, the average binding energy of the generated populations (gray crosses) is close to the prediction of the deterministic model (blue dot), and as expected the noise is higher for the conditions where a smaller amount of experimental measurement was performed (compare with number of measurements displayed in Appendix 1—figure 7).

Appendix 1—figure 11

Download asset Open asset

We then run the inference procedure on each artificially generated dataset, using the same setup and initial condition as the ones used to infer model parameters on the real data, under scenario C. The average results of the inference are reported in Appendix 1—table 1. For every model parameter, the average inferred value is close to the real value used to generate the data, demonstrating that the amount of experimental measurements at our disposal are sufficient for a good recovery of the model parameters.