Optimality of Mutation and Selection in Germinal Centers

The population dynamics theory of B cells in a typical germinal center could play an important role in revealing how affinity maturation is achieved. However, the existing models encountered some conflicts with experiments. To resolve these conflicts, we present a coarse-grained model to calculate the B cell population development in affinity maturation, which allows a comprehensive analysis of its parameter space to look for optimal values of mutation rate, selection strength, and initial antibody-antigen binding level that maximize the affinity improvement. With these optimized parameters, the model is compatible with the experimental observations such as the ∼100-fold affinity improvements, the number of mutations, the hypermutation rate, and the “all or none” phenomenon. Moreover, we study the reasons behind the optimal parameters. The optimal mutation rate, in agreement with the hypermutation rate in vivo, results from a tradeoff between accumulating enough beneficial mutations and avoiding too many deleterious or lethal mutations. The optimal selection strength evolves as a balance between the need for affinity improvement and the requirement to pass the population bottleneck. These findings point to the conclusion that germinal centers have been optimized by evolution to generate strong affinity antibodies effectively and rapidly. In addition, we study the enhancement of affinity improvement due to B cell migration between germinal centers. These results could enhance our understanding of the functions of germinal centers.


Introduction
As one of the adaptive immune responses [1,2,3,4], Affinity Maturation (AM) is the procedure in germinal centers (GC) to develop Immunoglobulins (Ig), i.e., antibodies, with increased affinities to a new antigen. Understanding the basic functional and physical principles of GC kinetics is not only important in medical science, but also contributes to the fundamental understanding of molecular evolution [5]. Mathematical models by Perelson and coworkers [6,7] made an important effort to describe the B cell population in a typical GC as a result of dynamic interactions between mutation and selection. This effort, complementary to the studies using Ig sequence data (e.g. [8]), is very important in revealing the functions of germinal centers. However, there are still two puzzles in the models. First, the experimentally observed somatic hypermutation [9] rate was viewed as [6] so high that even B cells expressing antibodies with improved affinities are easily spoiled by the majority of deleterious mutations. To resolve the conflict, it is proposed [6,10] that the mutation could be switched off periodically. Second, even if the mutation rate is tuned periodically, the calculated affinity improvement (up to 15-fold) is still not comparable with the observed improvements (,100-fold [11,12,13,14]). Deem and coworkers [15] used a version of random energy model to describe the alternative rounds of mutations and selection, assuming 100-fold affinity improvement. While this is helpful in sketching the AM procedure, it is interesting to explore the values of selection strength and antigen concentration that lead to the sufficient affinity improvement consistent with experiments. Third, a further analysis indicated that the models did not reproduce the ''all or none'' phenomenon [8], i.e., the fraction of B cells with strong affinity Ig is more likely to be high or low, but less likely to be intermediate.
Since most of the parameters in the models are estimated from experiments with considerable uncertainties, it is possible to reconcile the discrepancies by revising the parameter values. However, the number of parameters in Perelson's models is not small, and it is unclear how to find the best parameter values. We notice that the affinity dependent selection results from antigen binding kinetics, salvation, and recirculation of B cells. If we replace these steps by a phenomenological linear function of affinity to represent the selection, the calculations will be significantly simplified, and it will be possible to explore the parameter space and look for the optimal design. We hypothesize that AM has been optimized in evolution, and expect the B cell population dynamics with these optimized parameters to reproduce realistic AM.

The model
Each GC is believed to start from a few precursor B cells [16]. These precursor B cells first replicate at the perimeter of follicles [17], and the number of B cells in a GC reaches thousands when mutations and affinity dependent selections are turned on [7,18,19,20,21]. Our model is constructed as follows: First, the model describes the stage with mutation and selection, where the total number of germline (initial) B cells in the hundreds of GC in a spleen is ,10 6 . B cells replicate at 3,4divisions/day [22,23,24]. To be specific we use 4divisions/day in our calculations, corresponding to exponential birth rate r = 4ln2 = 2.8/day. Upon every division about one of the two daughter cells are mutated [9], so we find the total mutation rate m total &2:8=day from exp ({m total 0:25)&1{1=2, although in the calculation below we will explore different values of m total to search for the optimal design. It is estimated [25,26] that about 50% mutations are silent, 30% are lethal, and the rest 20% are the affinity-affecting mutations. The lethal mutation rate 0:3m total is an effective death rate. After taking care of the lethal mutations and neglecting silent mutations, we concentrate on the affinity-affecting mutations. Affinity can be described by the Ig-antigen binding free energy X. Define DX as the change of X upon a point mutation. The distribution W (DX ) of affinity change DX upon such single point mutations ( Figure 1) is estimated from the protein interaction (PINT) database [27] (see Text S1 for details). Note that only 4.9% of affinity-affecting mutations (equivalent to 1% of all mutations) improve the affinity, i.e. DX ,0. We assume W (DX ) is independent of the affinity X before the mutation. In support of this assumption, we found no significant correlation between DX and X in the data of PINT database (see Figure S1). Although it might become hard to further improve affinity when it is already very high, our model only requires this assumption to be valid in a range 2-4 kcal/mol from the germline affinity.
Second, selection is considered on the basis of the recent twophoton spectroscopic studies which indicate B cells undergo multiple rounds of mutation and selection [28] through migration within GC [3,4], consistent with the recycling hypothesis [6,7]. As the function of GC is to improve the Ig-antigen affinity, B cells compete for antigens and salvation from apoptosis, while other resources such as space is still not rate limiting during AM. The probability for a B cell to survive in each round depends on the probability to bind an antigen in the round, which in turn depends on the affinity of the B cell's Igs to the antigens. Besides binding free energy X, affinity can be characterized by the association constant K a~e xp ({X =kT) (in units of 1/M). According to a standard Langmuir adsorption isotherm the probability to bind an antigen in a round of selection is K a : C A =(K a : C A z1) where C A is the antigen concentration on the follicular dendritic cells of GCs. If we describe the affinity by the binding free energy X~kT ln (1=K a ), the probability for B cells to survive at each round is exp½{(X {X 0 )=kT for weak affinity, X wX 0 , where X 0~k T ln (C A ); and saturates at 1 for strong affinity, X vX 0 , reflecting the observation of an affinity threshold [29]. Define t as the typical time scale for each recycling round. Selection scales the population size as exp½{(t=t)(X {X 0 )=kT after time t for weak affinity B cells, so death/apoptosis rate is (X {X 0 )=(kTt).
Finally, including the replications, lethal mutations, and selections discussed above, the exponential growth rate of the B cell population B(X) is where the linear term with selection strength b~1=(kTt) reflects the rate of apoptosis of B cells discussed at the end of the preceding paragraph. Defining X Ã and K a Ã as the ''neutral'' affinity, B(X Ã )~0, we find X Ã~X 0 z(r{0:3m total )=b from Eq. (1a). Population decrease is expected for weak affinity X wX Ã , and X Ã is controlled by X 0 or C A , which is in turn determined by antigen density, reflecting an antigen dosage effect on AM.

Deterministic differential equation and its analytic solution
Now we are ready to write the mean-field differential equation for the population of B cells presenting Ig with affinity X at time t: where m:0:2m total is rate of the affinity-affecting mutations. (A similar equation was considered, in a different context and for fixed population size by Tsimring et. al. [30]) As mentioned above, replications, lethal mutations and selections are included in B(X),

Author Summary
The antibodies in our immune system could efficiently improve their abilities in recognizing new antigens. This is done with the help of proliferation, mutation and selection of B cells which carry antibodies, but we have difficulties in developing a quantitative description of this adaptation process which is consistent with the various aspects of experimental observations. Based on the knowledge from experiments, here we present a theoretical model to calculate the numbers of B cells with different antigen recognizing abilities all the time, and look for the best possible design that improves the antigen recognizing ability most efficiently. We find that the best possible design is consistent with the experimental observations, pointing to the conclusion that the immune system has been optimized in evolution. We then study the trade-offs leading to the optimization of the design. The results will not only improve our understanding of the functions in immune system, but also reveal the design principles behind the details. In addition, the study enhances our understanding of the population dynamics in evolution. and the affinity-affecting mutations are characterized by W (DX ) in Figure 1, where DX~X {X 0 is the affinity change.
We can solve the mean field equation (2) when Eq. (1) can be simplified as Eq. (1a), i.e. the effective growth rate (i.e. replication rate minus death rate) depends linearly on the binding free energy in the whole range. Then we rewrite growth rate B(X )~b 0 {bX where b 0 :r{0:3m total zbX 0 . We introduce a Fourier transform of the population: Then in terms of the Fourier transform population, Eq. (2) looks simple: where we get We seek solution of this equation in the form: and get for S: or where C is constant to be determined. The function R and the constant C are determined from the initial condition at t = 0, i.e., the germline distribution of affinities. The general results applicable for different initial conditions will be addressed shortly after; for the moment we adopt a most common initial condition, Gaussian distribution of affinity, or in Fourier space: Now the function R and constant C can be easily determined from the relation: which immediately gives us C~ln (N 0 =2p) and final result for R (and Q): Plugging these results to Eq. (7), we obtain the solution and accordingly , and another Fourier transform pair With help of the Fourier convolution theorem we obtain the final exact solution for the selection-diffusion equation (2) N(X ,t)~g(X ,t)z 1 2p whereÑ N(X ,t)~g(X ,t)e mt~N (X ,t~0) exp½(b 0 {bX )t is the population profile in the limit of m = 0. Due to the linearity of Eq. (2), Eq. (13) applies generally to any initial condition. The total population is The effects of mutations are in the exponent containing m, and this exponent applies universally to initial conditions other than Eq.

Intuitive derivation
An intuitive derivation of (14) helps us understand its physical meaning. For convenience we use discrete values of DX i~{ ih in Figure 1, and write the total mutation rate as a sum of individual mutation rates, m~P i m i where m i is the rate of mutation that change affinity by DX i~{ ih. Then Eq. (14) becomes: À? m i bih e biht to ln N(t). Therefore in the long time limit the beneficial mutations with leads to a super-exponential population growth N(t)* exp½exp (t), in comparison to the exponential population growth without mutation. As the high affinity edge moves toward stronger affinity, the subpopulation at the edge grows faster, and the speed of edge movement becomes faster.
The effects of different mutations in Eq. (15) can be factored out, because a) the rate m i for a mutation of DX i~{ ih to emerge is independent of the affinity distribution N(X ,t), and b) the contribution of a series of mutations fDX i g occurring at ft i g to the subpopulation size at a later time t is P i exp½{bDX i (t{t i ), a product of individual mutation factors. Therefore, to understand Eq. (14), we can focus on effect of only one mutation type with affinity change DX i~{ ih and individual mutation rate m i .
First, in the limit of bR0 in Eq. (1), the benefit of mutations is turned off, and the average number of mutations B cells experience in the interval t is m i t, and the probability to experience j mutations is a standard Poisson distribution Second, for the realistic situation with nonzero b, replication rates of B cells are changed by mutations, therefore every m i t in (m i t) j is replaced by where t9 labels the moment that the mutation occurs. Therefore This matches Eq. (15) and more generally Eq. (14). From Eq. (18), the subpopulation of B-cells which undergo j mutations is with affinity jih stronger than the subpopulation without mutations. If the affinity of the initial population are all the same X in , then N j (t) is the subpopulation with affinity X in zjih.
For i.0, the peak of subpopulation, i.e., largest subpopulation, is at i.e., the peak moves exponentially fast in the long time limit. From Eq. (19) a subpopulation with j mutations grows for two reasons, (a) fed from subpopulations with j-1 mutations and (b) self-replication. At short enough time, bht%1, the subpopulation grows dN j (t)=dt&N j{1 (t)m i mainly for the former reason; and for large enough time, bht&1, the subpopulation growth dN j (t)=dt&N j (t)jbh is mainly contributed by the latter reason. Therefore, the artifact of self-replication is insignificant when time step is much shorter than the characteristic time duration 1/bh.

(13)
Numerical calculation for finite population An artifact in the above derivation is that it allows arbitrarily small N(X ,t), and a small subpopulation N(X ,t)%1 within a bin of strong affinity X can self-replicate rapidly. However, the B cell numbers in GCs are non-negative integers, therefore the expected number of B cells within a bin with N(X ,t)%1 should be zero, and cannot become the seed of a rapid growth. So the above derivations actually describe the population dynamics in the limit of infinitely large population size. It does not take into account the fact [31] that the B cell population size in a GC is no more than 10 4 .
To correct this artifact and calculate the B cell population dynamics numerically for various finite initial population sizes including the discreteness effect [7,30], we do not allow small subpopulation N(X ,t)v1 in an affinity bin to self-replicate in our numerical calculation. Instead, it only represents an accumulative probability for the subpopulation in the bin to emerge. Our calculation is done using discrete time steps. From Eq. (15), the subpopulation which go through fj i g mutations DX~{ih between time t and tzDt is where B(X) is given by Eq. (1), and N 0 (X ,t) is the population distribution excluding the bins with less than one B cells. The population distribution after a time step is The time step Dt is set to 0.01/bh, much smaller the characteristic time duration 1/bh,such that the growth of subpopulations with relatively strong affinities is dominated by mutation influx rather than self replication. We include multiple mutations in one time step. This calculation is rapid to perform even for large population sizes, and allows us to explore the parameter space. The population dynamics here differs from the case of constant population size [30,32,33,34,35] Indeed, the ''neutral'' affinity X Ã , which corresponds to zero population growth rate B(X Ã )~0, does not change over time. This makes it straightforward to find solution Eq. (13)(14). Random drift beyond mean field calculations is important in the case of constant population size [32,33,34,35], especially if population growth rate is small, s%1, at the high affinity edge, because birth rate and death rate are very close. However, the growth rate in units of proliferation rate in our model B(X )=r, which is similar to s, differs by about 1/8 between two nearest affinity bins for a typical value b = 0.7/day/(kcal/mol), so s at high affinity edge typically is not close to zero. Hence the effect of random drift is not overly significant in our model, although a future treatment including stochastic calculation will give more precise results.

Results
We calculate the total population size of B cells N(t) for various initial population sizes numerically (dotted lines in Figure 2). The total number of B cells first decreases because the initial affinity is weaker than the neutral affinity K a Ã . The average affinity is improved continuously (see Figure S2) rather than abruptly, in agreement with experiments [28]. Once the average affinity reaches K a Ã the population begins to increase. The picture of decrease and increase of the B cell population was observed experimentally [36] and theoretically [7], although the experimental data on GC temporal development [36] is too limited to verify the model. The lowest total population corresponds to the neutral affinity K a Ã ; and we call it the population bottleneck, because it is the most challenging moment for the population to survive. The analytical result Eq. (14) (solid line) can describe the population size at the decreasing stage. A smaller initial population leads to a slower growth after the bottleneck (yellow). For a small enough initial population (red), the B cell population is extinct when approaching the bottleneck, and cannot recover thereafter. Therefore, the initial population size should be large enough to ensure some B cells can survive through the population bottleneck. Similarly, for a given initial population size, a weaker initial binding leads to a deeper bottleneck, and takes a longer time to recover. If the initial binding is too weak, the bottleneck will be too deep, and the population will go extinct.
The initial population is set to 10 6 B cells in realistic calculations to include the existence of hundreds of GCs in a spleen [37] and the peak number of thousands of B cells per GC [38]. Different GCs in a spleen might not be perfectly synchronized [39]. If the B cell production rate in a spleen is limited by supply of resources, we conjecture the time that GCs start mutation and selection might vary between day 3 and day 8 or later, so that the population peaks of GCs is smeared. This agrees with the observation [40] that the total population at any moment does not exceed 2:5|10 5 B cells. The population bottleneck within a GC might also be smeared by the continued immigration of B cells from nonfollicular sites [20,41], making the bottleneck less pronounced or harder to observe directly. If all the B cells in a GC die away, the antigens are not exhausted, and it is proposed that more B cells immigrate to the GC [39], probably from GCs which have passed the bottleneck and have many B cells, although the current experiments cannot determine whether there is migration between GCs [39]. The calculations below are first performed with 10 6 initial B cells, which is valid in the limit of fast migration between GCs in a spleen. In this case, the few GCs which by chance pass the bottleneck earlier than others may make significant contributions to the affinity improvement of the whole spleen. Then we study the case in the opposite limit, i.e. the affinity improvement of a typical GC with 3000 initial B cells, assuming no B cell migration between GCs. By comparing these two cases we will quantify the contribution of migration to affinity improvement.
We explore the parameter space to look for the optimum design of GCs that maximizes the affinity improvement (Figure 3) in the limit of fast B cell migration between GCs (10 6 initial B cells). AM terminates probably due to exhaustion of available antigens [42,43] or emigration of B cells [42,44]; and the termination could be described to occur after a certain time scale, or when the B cell population size is big enough-probably comparable to its original size. So if the population size recovers the initial value in less than 14 days, we assume that the total affinity stops to change in the calculation of Figure 3a. The affinity improvement is indicated by the total affinity A, i.e., the sum of Ka for all B cells, such that the situation with very few cells at day 14 should not be regarded as efficient improvement. The improvement of affinity in Figure 3a is calculated for various initial binding levels and mutation rates. The highest affinity improvement is 450-fold. This optimal improvement occurs at mutation rate m total &2:8=day (corresponding to the observed fraction of 50% daughter cells mutated), for the binding level between germline (initial) antibodies and antigens X in~X Ã z1kcal=mol or K a (0)=K a Ã~e xp ({1kcal=mol=kT)0 :18, and selection strength b~0:7=day(kcal=mol) (see Figure S3 and S4). This result agrees with several experimental observations. First, the affinity improvement agrees with the observation of ,100-fold [11,12,13,14] improvement. Second, the theoretically optimal value to provide maximal affinity improvement agrees with the observed in vivo somatic hypermutation rate [9]. Third, the improvement of affinity corresponds to ln(450)kT = 3.6kcal/ mol of free energy improvement. Combining with the typical affinity improvement Y &0.4kcal/mol of an affinity improving mutation, we can estimate that a final B cell contains 3.6/0.4 = 9 mutations in their V regions of Ig genes, in agreement with the observed ,9 mutations per Ig gene [9,45,46,47,48]. From the definition b~1=(kTt), the optimal selection strength b = 0.7/day/ (kcal/mol) corresponds to an optimal time of a recycling round t~1=(kTb)&2:4day, compatible with the earlier model [7]. Figure 3b helps to reveal the design principles of GCs, where the affinity improvement is calculated when the B cell population size recovers the initial value after going through the bottleneck, no matter how long it takes. In general, the affinity improvement is not effective for too strong initial bindings, which results in shallow or no bottlenecks and rapid population recoveries, hastily terminating the AM before accumulating adequate improvements. As the initial binding becomes weaker, the population bottleneck is deeper, and the affinity improvement is more effective, but takes a longer time. The improvement shown in Figure 3b could even exceed 1000-fold for a weak initial binding X in~X Ã z 1:5kcal=mol, although such improvements take much longer than 14 days and would be interrupted in the calculation with fixed AM time as shown in Figure 3a. If the initial binding becomes even weaker, the population bottleneck is so deep that the whole population goes extinct and no longer recovers (grey scales in Figure 3b obtained numerically). Therefore, for a given mutation rate, the most effective affinity improvement occurs when initial binding level is near the critical value (red in Figure 3b), where the population can barely survive through the population bottleneck, i.e. only a few GCs can survive. A somewhat stronger initial binding improves affinities less effectively, but it takes a shorter time for the AM to finish, and faster AM is advantageous. Taking all these into consideration, in the optimal design of GCs, the most commonly appeared initial bindings should be adjusted to be somewhat stronger than the critical value, so that affinities are improved effectively, timely, and safely. This picture agrees with the observed dependence of B cell fate on initial antibody-antigen binding level or antigen density [49][50], where too strong initial bindings (beyond Ka* in our model) do not lead to GC formation, moderate initial bindings (e.g. X in~X Ã z0:5kcal=mol in our model) result in GC response which finishes quickly, while weaker initial affinities (e.g. X in~X Ã z1kcal=mol in our model) result in tempered GCs. To achieve the optimal design, the Ka Ã values might have been adjusted in evolution by tuning the antigen density or modulating the diversity of the germline pool achieved through somatic recombination. When mutation rate is so high that 80% or more daughter cells are mutated, the lethal mutations preclude sustainable replications of B cells, i.e., the population growth rate in Eq. (1) becomes negative for any affinity. This lethal mutagenesis region is shown in black in Figure 3b.
B cell migration between GCs could be beneficial to the AM, as indicated by the calculated affinity improvement of a typical GC with 3000 initial B cells in the limit of no migration (Figure 4). A stronger initial binding X in~X Ã z0:5kcal=mol is needed to ensure a typical GC to survive through the bottleneck; and the optimal affinity improvement becomes 70-fold for the GC, still consistent with the experimental value (of the order ,100-fold [11,12,13,14]). This is achieved in 16 days, only slightly longer than two weeks, therefore separate diagrams like in Figure 3a and b is not necessary. The maximal improvement in this case is ,6-7 times lower than the result in the limit of fast migration, where a factor of ,3 comes from the change in initial binding level, and a factor of ,2 comes from the difference between the average of all GCs and the ''typical'' or median individual GC. In other words, fast B cell migration between GCs could enhance the affinity improvement by a factor of about 6-7 in our model.
Our model is consistent with the ''all or none [8]'' phenomenon observed in experiments [20,51,52], i.e., the fraction of strong affinity cells in a GC is most likely either very high or very low, where the strong affinity B cells are characterized by a certain key mutation [51,52] or a unique piece of Ig gene sequence [20]. Presumably the GCs dominated by strong affinity mutants have gone through the population bottleneck, and this phenomenon is the most pronounced in the case without B cell migration between GCs. Indeed, it is possible for different individual GCs to fall into two categories according to whether they have passed the population bottleneck, and the choice of category for a GC depends on both the initial binding level and stochastic effect. Moreover, if many GCs in a spleen start from the same initial condition but each has a random starting time, then the fractions of strong affinity B cells in the GCs at a given moment is expected to be distributed as in Figure 5, where affinities beyond Ka* are defined as strong. Figure 5 is obtained as follows. During the development of a typical GC (Figure 2 and S2), we can track the fraction F(t) of strong affinity B cells in the GC all the time.
Observing the whole ensemble of GCs at a given moment is equivalent to observing a single GC at many arbitrary moments. Therefore, we transform F(t) into t(F), and the distribution is P(F),dt(F)/dF up to a normalization factor. As we see in Figure 5, the probabilities to observe high or low fraction values are significantly larger than that of intermediate fraction values. Therefore, our model is consistent with the ''all or none'' phenomenon. If B cell diffusions between GCs exist, this phenomenon is somewhat smeared.

Discussion
Since the calculated AM procedure at the optimal set of parameters is consistent with the various aspects of experimental observations, it is likely that evolution has chosen the optimal design. Moreover, exploration of the parameter space, illustrated in Figure 3, 4, S3 and S4, helps us sketch some general design principles of the GCs. First, for the affinity to be improved most effectively, the B cell number should first decrease to reach a population bottleneck and then increase. A weaker initial binding leads to a higher affinity improvement, but a too weak initial binding makes it impossible for the B cell population to go through the bottleneck and recover. Second, the seemingly high mutation rate is actually set to optimize the success rate of AM. On one hand, the optimal mutation rate, in agreement with somatic hypermutation in vivo, is quite high because the improvement of affinity comes only from mutations. On the other hand, if the mutation rate gets even higher, the large number of lethal mutations will spoil the cell replications. Third, we expect that the selection strength b is the optimal value b = 0.7/day/(kcal/mol), which sets an important guide to future simulations of GCs. Indeed, a large enough b ensures adequate preference for stronger affinity B cells in selection; but if b is too large, the initial weak binding B cells vanish quickly before enough beneficial mutations can accumulate.
We obtained mutation effects W (DX ) from the PINT database which presents the data on interactions between all kinds of proteins [27] because we do not have adequate real data of affinity change upon mutations of Ig. The existing analysis of Ig sequence data in the immune response of PhOx and NP [8] is consistent with the distribution obtained from the PINT database. Indeed, 3.2% (for PhOx) or 1.0% (for NP) of the affinity affecting mutations could improve affinity strongly (10-fold), while in our model 4.9% improves affinity and 1.4% improves affinity strongly (5-fold).
The random energy NK model of Deem and coworkers [15] included rugged landscapes without special assumptions on the Ig mutagenesis. On the other hand, using an approximate distribution Figure 4. The improvement of affinity for an isolated GC, i.e., in the limit of no B cell migration between GCs. The improvement of A (sum of Ka) is shown in the same color code as in Figure 3a, assuming AM is terminated when the population recovers the initial size (3000 B cells). The optimal improvement of affinity occurs when about 60% daughter cells are mutated at divisions, and takes 16 days. The grey scale shows the probability for a GC to survive through the bottleneck. doi:10.1371/journal.pcbi.1000800.g004 Figure 5. The distribution of F, fraction of strong affinity B cells, for many GCs which follow similar development patterns but each starts from a random time, is consistent with the ''all or none'' phenomenon. Every GC has initial binding level X in~X Ã z0:5kcal=mol, 50% mutated daughter cells, and selection strength b~0:7=day(kcal=mol). The calculation is terminated when the population in the GC recovers the initial size (3000 B cells), and F reaches 85% at the termination moment. doi:10.1371/journal.pcbi.1000800.g005 of affinity improvements upon point mutations, we explored the values of selection strength, initial binding level and mutation rate that maximize the affinity improvement. Similarly to our work, the Oprea-Perelson (OP) model [7] explores the parameter values to optimize AM. What are the major differences between the OP model and our model? First, the affinity improvement we achieve is higher than in the OP model, which resolves the discrepancy between theory and experiment. If we examine the OP model in our framework, the reason that the OP model did not achieve higher improvements might be that it used a smaller selection strength b than the optimal value. Indeed, the choice to enter a recycle in the OP model is stochastic, so the chance for weak affinity B cells to survive is enhanced by choosing not to enter a cycle and therefore avoiding selection. A larger selection strength b brings a high risk of AM failure, but provides a higher affinity improvement, and this can be achieved if all B cells experience similar number of recycles rounds. The low optimal value of mutation rate in the earlier model might result from the small selection strength. Second, the major simplification we make is to use Eq. (1) to replace the binding kinetics between Ig and antigens, salvation, recycling, the distinction between different B cell phenotypes such as centroblasts and centrocytes, and the separation of dark and light zones, which were considered in details in earlier theoretical studies [7,8,42,43,44,53,54]. We assume that competition for antigen and T-cell help effectively results in the dependence of death rate of B cells on their affinity to an antigen, and we provided an argument that such dependence should be linear in binding free energy of Igantigen interaction. We expect that this simplification captures the key factor determining the AM. Such a coarse-grained description reduces the parameter space greatly, making it possible to search the whole parameter space for optimal design of GC. With this basic picture, we can build more detailed simulations to reproduce the complete procedures of AM. Figure S1 The scatter plot of affinity X and the change of affinity DX from PINT database, which does not show significant correlation.

Supporting Information
Found at: doi:10.1371/journal.pcbi.1000800.s001 (0.05 MB TIF) Figure S2 Affinity distribution in the population N(X,t) at t = 3 (solid), 6 (dashed), 9 (dotted), and 12 (dash-dotted) days, starting from germline affinity X in = X * +1kcal/mol or Ka in /Ka * = 0.18, with initial population N 0 = 10 6 , selection strength b = 0.7/day(kcal/mol) and effective mutation rate m total = 2.8/day. The average affinity improves with time; while population size shrinks then grows. Found at: doi:10.1371/journal.pcbi.1000800.s002 (0.06 MB TIF) Figure S3 Optimization of b and m. The color indicates the improvement of total affinity. Here different initial affinities are tried for each mutation rate and b, and the one which gives largest affinity improvement is chosen. b = 1.2/day/(kcal/mol) is the global optimal selection strength. A minor local peak at b = 1.2/ day/(kcal/mol) might be an artifact due to discrete (rather than continuous) choices of initial affinity values. Found at: doi:10.1371/journal.pcbi.1000800.s003 (0.21 MB TIF) Figure S4 The improvement of total affinity as a function of selection strength and initial affinity. Here mutation rate is chosen as the optimum value, i.e. m = 0.55/day/gene or 50% mutated daughter cells.