We model B cell proliferation and death during the GCR using a birth-death (BD) process (Renshaw, 1991). AID mutates the gene encoding for the BCR (Muramatsu et al., 2000) and as a consequence, affinity for the antigen changes. The resulting increase (or decrease) in affinity translates to a higher (lower) fitness of the B cell. Regarding the stochastic variation of BCR in affinity space as a form of diffusion, the model resembles a “birth-death-diffusion process” (Adke and Moyal, 1963).

In the first days following immunization, while the GC is still coalescing, B cells proliferate without competition, creating a pool of cells on which AM may operate. Few or no mutations are introduced to the BCR sequence at this early stage. We start from a simple birth/death (BD) process using an agent-based model. Each cell is associated with a birth rate λ and a death rate μ (see Figure 1A). We assume that a GCR starts with M different clones and the system evolves for a period of 6 days, which we denote by T_growth (Jacob et al., 1991; see Figure 1B).

During the growth phase, the probability distribution P_ni(t) of the number of cells n_i that belong to clone i evolves in time according to the master equation (Bailey, 1990):

where (in the absence of interactions): λ_n = nλ and μ_n = nμ and P_0i is the probability of extinction of clone i. The average number of cells 〈n_i〉 in clone i, after time t is given by (Bailey, 1990).

The time dependent extinction probability of a clone is

and the size distribution of a clone lineage is

Both equations are the solution of Equation (1). After T_growth, there is a supply of cells on which AM can work, while some clones disappear. This distribution function is the starting point for the competitive phase of the GCR.

For our parameter choices (see Table 1), which represents a GC development, the average lineage size of a clone at the end of the growth phase (6 days) is 〈n_i(6 days)〉 = 20 cells, the total number of surviving cells is 〈N(6 days)〉 = 1000 cells, while p₀(6 days) ≈ 2/3 corresponding to an average of 50 × 1/3 ≈ 17 surviving clones. This number is lower than the number of surviving clones in Tas et al. (2016) which was 50–200 but as we are interested in the qualitative behavior of the system, we choose a smaller number to facilitate the numerical calculations.

After day 6, B cells survival depends on TfhC signals that are a shared resource. Indeed, it has been shown (Victora et al., 2010; Gitlin et al., 2015) that TfhCs have a role in regulating the duration of cell cycle in B cells during AM and controlling their behavior in the GC. To mimic B cell competition over the limited resource of TfhCs, we used the stochastic logistic growth process (Nåsell, 2001), which constrains the B cell population size. The death rate decreases with the population size, from a basal rate of μ₀, to roughly the birth rate λ₀ for a mature population:

where N is the population capacity. Here n = (n₁, n₂, …, n_M) is the vector of cell number n_i for the M lineages. The competitive phase continues for a period (T_comp), which we take to be 16 days (Tas et al., 2016). The total number of cells in the GC grows gradually until reaching the capacity N, where it remains approximately fixed.

Occasionally, B cells undergo a proliferative burst that is proportional to the amount of presented antigen and thus to the BCR affinity (Victora et al., 2010; Gitlin et al., 2015). B cells move then to the DZ, remain there and divide multiple times (4–6) before going back to the LZ to go through another round of selection (Gitlin et al., 2014, 2015; Tas et al., 2016). We model this process as an increase in the birth rate (see Supplementary Information “Heterozygosity of a Moran process”). Since cell-cycle is modified (shortened) in this process, we take the birth rate of cell i as

where w_i is the affinity of cell i, 〈w〉_Population is the mean affinity of the population and λ₀ is the basal birth rate. Indeed, the average birthrate of B cell clones in a GC, was found to be similar (Anderson et al., 2009) in B cell clones with different affinities. The normalization serves to keep the average population birth rate constant at λ₀. Since the clone birth rate λ_i is related to the clone affinity w_i, we designate this scenario “birth limited selection.”.

During the GCR, cells with poor affinity do not receive a survival signal from T helper cells because they do not display a sufficient amount of peptide-MHC molecules. Previous studies model this process by noting that the probability of a B cell being able to successfully compete with other B cells that have internalized antigen and receive T cell help, grows monotonically with the affinity of its BCR for antigen (Zhang and Shakhnovich, 2010; Wang et al., 2015), with surviving cells proliferating at approximately the same rate (Batista and Neuberger, 1998). Additionally, it was found (Anderson et al., 2009) that on average, B cell clones with different affinities differ in their death rate, where the low affinity clone dies at a higher rate than ones with intermediate affinity. Such a scenario is considered “death limited selection” in our scheme with a death rate μ that depends inversely on the affinity. To study the consequences of such a selection mechanism, we constructed the following model

where α is a constant, μ_i is the death rate of a cell with affinity w_i and μin is the GC-size dependent death rate keeping the population size fixed. Thus, higher affinity is related to a lower death rate.

We also examine a model where the birthrate is normalized over the population and as a result, the average of affinity dependent element of death rate, is constant.

During AM B cells mutate their BCR encoding genes. The effect of a single mutation on fitness in models of Wright-Fisher-like selection is often taken to be small (Park and Krug, 2007; Hallatschek, 2011; Goyal et al., 2012; Tas et al., 2016), which allows analytical treatment of the population dynamics as a diffusion problem. In this spirit, we modeled the effect of mutation as a change in the affinity upon cell division, where one of the daughter cells has the parent affinity and for the other daughter:

where N is a normal distribution with zero mean and standard deviation of 2D, with D akin to an effective diffusion coefficient determining the magnitude of affinity change. Within this model, affinity can increase or decrease with equal probability at every division.

At the end of the growth phase we are left with 17.2 clones on average, consistent with the stochastic simulations (Figure 2C). At this point, the size of remaining lineages has a large variability according to Equation (4). We find that changing the “diffusion coefficient” D has a strong impact on the homogenization rate (Figure 2D). For larger values of D, fewer clones survive to be part of a mature GC (Figures 2C,D). The participation ratio, which is the probability that two randomly chosen B cells belong to the same clone, also suggests rapid loss of diversity for large value of D (Figure S2). Surprisingly, we find that the variability of different GC realizations increases with time (Figures 2A,E). Naturally, at long times diversity is lost and only a few clones are left, and the variation in the fraction of the most dominant clones decreases (Figure 2E). Thus, the highest number of possible outcomes, in clonal variability, occurs at an intermediate time, which for high values of D, happens at day 11 of the competitive phase.

The case of a GCR without mutation was also studied experimentally, in a setting in which multiple clones all having the same BCR seeded the GC and the AID gene was genetically deleted (Tas et al., 2016). Interestingly, even with no changes in affinity, there is a gradual and slow homogenization (Figure 2B, empty circles). To study this scenario, we performed numerical simulations in the absence of mutation (D = 0) and saw a gradual take over by the dominant clone (Figures 2C–E), as seen experimentally. As all clones have the same affinity, clonal loss and homogenization in this case is due to random drift (Renshaw, 1991). To gain intuition regarding the selection and fixation process, we recall known results for a case where the population size is fixed, corresponding to a Wright-Fisher process (Bailey, 1990). When affinity differences between the clones are neglected and a starting group of M clones all occupy the same fraction of the population size, the mean time to fixation of a single clone is given by τ_fixation = 2(M−1)log(M/M−1). With non-uniform initial numbers of clones, the probability of a clone to fix is equal to its initial fraction in the population (Bailey, 1990), which in our model is the probability distribution at the end of the growth phase (Equation 4).

A clone whose affinity is relatively higher than that of the other clones in the GC has a better chance of being selected and becoming dominant (Equation 6). Since all clones had the same initial affinity, during the first few days of the competitive phase the affinity distribution of the population relaxes from a delta function (δ(w−w₀)) (Figure 3A). A GC reaches its capacity only a few days after the beginning of the competitive stage (Figure S3A). Before that, diversity loss continues at the same rate of the growth phase and is D independent (Figure 2C). Beyond a certain threshold, the homogenization rate is independent of the birth-rate (Figure S3B).

At later times, the affinity distribution moves as a traveling wave (Tsimring et al., 1996; Hallatschek, 2011; Figure 3A), as fitter strains at the higher end of the affinity distribution function constitute the moving edge while the cells on the other end die. The velocity of the affinity wave depends on D (Cohen et al., 2005; Figure 3B) and since affinity changes upon cell division, it depends also on λ (Figure S3D). As expected for a traveling wave solution, the average affinity grows linearly with time. During this period in the GCR, since the affinity of all clones change due to the same stochastic process, a clone which after a single mutation has an affinity larger than the mean, is likely to outperform the other clones. Such deviations from the mean affinity, are governed by large jumps, which are related to the value of D.

To study if loss of clonal diversity in a GC is the result of homogenizing selection toward high affinity clones, we computed the correlation between the number of surviving clones in a GC and the average affinity of the most dominant clone at the end of the selection phase (Figure 3C). On day 16, the affinity of the dominant clone is a good proxy for the average affinity in the population. Interestingly, while we observe a weak negative correlation (r = −0.53), many GCs maintained diversity in spite of having high affinity clones.

We can consider the width of the affinity distribution of a GC population to be a proxy for its clonal diversity. It was shown that the ratio of the mean affinity to its standard deviation (STD) grows during AM when the amount of antigen used in the immunization was relatively low (Kang et al., 2015). Indeed, the STD of a stochastic variable grows with time (Schuss, 2009), while the growth of the average affinity is evidence of selection (Desai and Fisher, 2007). When the mean grows faster than the STD it is a sign of strong selection. We estimated this ratio from our simulations. Initially, as the affinity distribution spreads from a delta function and before the GC reaches its capacity, the ratio decays, but following the initial relaxation phase, the mean affinity increases faster than the spread of the distribution (Figure 3D). Thus, our system operates in the strong selection limit as in the experimental system studied in Kang et al. (2015).

To what extent does the initial growth phase determine the later state of the GCR? We define the state of a GC as the vector of proportions of clonal lineages at time t; n(t) = (n₁(t), n₂(t), …, n_M(t))/N_tot(t). The correlation with the initial state of the GC is quantified by

and is observed to decay with time (Figure 4A). The initial fractions of clones change when stochastic increases or decreases in the affinity of cells give relative advantages or disadvantages to particular clones (Equations 6, 9). Thus, for larger values of D, C decays faster. Similarly, the decay rate of correlations is inversely proportional to the basal birth rate (Figure S3C) and to N, since the fixation probability of a species in a population is inversely proportional to population size (Desai and Fisher, 2007) (data not shown). This result raises the question of whether a GC effectively filters the best clones, as the system has a finite probability to be “stuck” in an unfavorable state.

To further explore the relation between clonal competition and affinity we performed numerical simulations where each B cell of the M initial ones had different initial affinity w₀. Following growth, we studied clonal dominance in the competitive phase. Interestingly, while the clone with the highest initial affinity (w₀ = 1.5) had the highest probability of becoming the dominant clone, the clone with w₀ = 1.25 still had a chance of becoming dominant (Figure 4B). This exemplifies the stochastic nature of the selection process. The effect of the initial affinity w₀ in determining the second, third and fourth dominant clone is smaller (Figure 4B).

We addressed the relation between affinity and dominance by estimating the correlation between the average clonal affinity and the fraction occupied by the first to fifth dominant clones. Interestingly, we see that often clones with high affinity compose a small fraction of the GC at the end of the GCR (Figure 4C). We also see that this depends on the value of D, and for a larger value the positive correlation between dominance and affinity is stronger (Figure 4D).

To study the effect of a death-limited model on the progression of the GCR we preform stochastic simulations using an affinity-dependent death rate (Equation 7). The GC population's affinity continues to increase throughout the simulation (Figures 5A,B). We assume that clones with higher affinity have a smaller probability of dying, as they are likely to receive a survival signal from the TfhCs. Thus, in our death-limited model, affinity increase results in decrease of the death rate (Equation 7). Thus, we observe a gradual decrease of the death rate distribution of the cell population (Figure 5C). We found two homogenization regimes (Figure 5D). While the GC has not yet reached its capacity and death rate distribution of the cell population relaxes from a delta function, which was the initial condition (w(t = 0) = δ(w − w₀)), to steady state, homogenization is slow. Indeed, for D = 0 the homogenization rate remains constant. In this case, diversity loss is related to random drift only. At later times, homogenization occurs at a fixed rate, dependent on D (Figure 5D). The exponential relation between affinity and death rate in this death-limited selection model acts to modulate large affinity jumps. Thus, homogenization occurs at a slower rate than that of the birth-limited model we studied in the previous section.

To investigate if the difference between the death and birth limit selection model is due to normalization of the birth-rate (Equation 6), we performed simulations where the death rate of cell i was given by Equation (8). When the average affinity dependent death rate remains μ₀, the homogenization rate increases (Figure S4) with respect to the un-normalized case, but still remains slower than that for the birth-limited model. There are experimental evidences that the average birthrate is constant in the GC, independent of the affinity of B cells (Anderson et al., 2009). However, such is not the case for death limited selection, since no survival signal is given to B cells by T cells when no Ag is captured. This presumably will occur when the affinity is small. Thus, it is likely that a dependence of death rate and affinity (Equation 7) exists in the GC.

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