A Sensitivity Analysis Comparison of Three Models for the Dynamics of Germinal Centers

Germinal centers (GCs) are transient anatomical microenvironments where antibody affinity maturation and memory B cells generation takes place. In the past, models of Germinal Center (GC) dynamics have focused on understanding antibody affinity maturation rather than on the main mechanism(s) driving their rise-and-fall dynamics. Here, based on a population dynamics model core, we compare three mechanisms potentially responsible for this GC biphasic behavior dependent on follicular dendritic cell (FDC) maturation, follicular T helper (Tfh) cell maturation, and antigen depletion. Analyzing the kinetics of B and T cells, as well as its parameter sensitivities, we found that only the FDC-maturation-based model could describe realistic GC dynamics, whereas the simple Tfh-maturation and antigen-depletion mechanisms, as implemented here, could not. We also found that in all models the processes directly related to Tfh cell kinetics have the highest impact on GC dynamics. This suggests the existence of some still unknown mechanism(s) tuning GC dynamics by affecting Tfh cell response to proliferation-inducing stimuli.

1. Introduction

Immune responses to a T-cell dependent antigen (Ag) are initiated in secondary lymphoid tissues like lymph nodes and the white pulp of the spleen (1). In the spleen, Ag-activated T helper (Th) and B lymphocytes from peri-arteriolar lymphocyte sheath (PALS) and follicles, respectively, change their chemokine responsiveness so that they are forced to migrate to the boundary between PALS and follicles, greatly facilitating Ag-specific B cell-Th cell encounters (1). There, activated Th and B cells interact with each other in an Ag-dependent way, and mutually induce a co-stimulus-dependent proliferation (2, 3). Within 1–2 days some progeny of those activated Th and B lymphocytes start migrating into adjacent follicles, which are characterized by a network of follicular dendritic cells (FDC) (4). There Th cells differentiate into so-called follicular T helper (Tfh) cells, characterized by the expression of the chemokine receptor CXCR5, the inhibitory molecule PD-1, and the transcription repressor factor Bcl6 (5, 6). Ag-specific B and Tfh cells that migrate to follicles form germinal centers (GC) (7). These are transient anatomical microenvironments with an average life of up to 3 weeks in a murine primary immune response to protein Ags (8, 9) during which intense proliferation, apoptosis, and V(D)J hypermutation of B cells takes place (7, 10). These processes, together with a still incompletely understood selection process, are essential for the affinity maturation of antibodies (Ab) (11–15), and therefore their kinetics within the global GC dynamics must impact the resulting Ab affinity maturation.

Most models of GC dynamics have been developed to investigate how Ab affinity maturation for Ag proceeds [e.g., (16–18)] rather than trying to understand what are the forces driving GC dynamics itself. Thus, in those works the focus is on different GC processes potentially accounting for a given view of Ab affinity maturation, and as a consequence GC dynamics has been often modeled following one or another hypothetical selecting processes. This procedure may be misleading in that it generally imposes, secondarily, a mechanism driving the dynamics of GCs. There is, therefore, a need for a better understanding of the FDC, B cell, and T cell inputs involved, not in affinity maturation (19), but in the GC dynamics itself.

Following that view, we present here an attempt to uncover the main mechanism(s) that drive the dynamics of GCs. This can provide then an unbiased and general framework on which to base implementations designed to investigate the relevance of different, potential affinity maturation mechanisms.

The question of what are the main processes driving the GC dynamics has received little attention in the past (20, 21). Moreover, because that previous work focused in the termination of the GC reaction, the different interactions between FDCs, B cells, and Tfh cells, and the processes known to be set in motion by those interactions in GCs, were highly simplified. However, several mechanistic hypotheses exist in the literature that could explain the transient, rise-and-fall dynamics of GCs. They focus on different intrinsic dynamical aspects, some of which can be observed during the GC reaction, and range from increasing differentiation of FDCs (8, 22) or Tfh cells (23–27) likely due to repetitive interactions with B cells and/or the changing GC cytokine milieu (28), to Ag depletion (21, 29). Here we implement a population dynamics model core of GCs, and based on it we develop and analyze three different models. Each of these models adds to the model core one out of the three different mechanisms, FDC maturation-based, Tfh maturation-based, and Ag depletion-based, described above as being potentially responsible for the biphasic nature of GC dynamics.

The behavior of a relatively complex dynamical model frequently depends more on some parameters than on others. Knowledge of the parameters' impact on the model behavior is important, for instance, in helping to simplify a complex model (30), and/or in helping to uncover the relative significance of the various input parameters in determining the model dynamics. This latter use would reveal those parameters that are potential control parameters, that is, parameters through which the GC dynamics can be regulated. Here we quantified the parameters' impact by the so-called first and second order relative sensitivities (31, 32), or simply sensitivity (impact of changes in parameter value on a model's behavior) and synergy (impact of changes in parameter value on the sensitivity of the model with respect to another parameter). Our results indicate that the processes having the highest impact on the dynamics of GCs are those related to the kinetics of Tfh cells, irrespective of the model. This strongly suggests that the global dynamics of the GC reaction can be tuned by mechanisms impinging on the way Tfh cells respond to activating stimuli. This shifts therefore the focus toward uncovering such direct or indirect mechanisms affecting Tfh cell behavior, which in turn would lead to a deeper understanding of how the GC B cell repertoire changes with time.

2. Models and Methods

2.1. Models

2.1.1. Model Core

The conceptual framework on which we base our modeling explicitly takes into account the main interactions between FDCs, B cells, and Tfh cells which have been experimentally established to occur in GCs (5, 33). This conceptual framework is depicted in Figure 1. Free GC B cells, denoted B, interact with rate c₁ with Ag deposited in the form of immunocomplex bodies or iccosomes on the dendrites of FDCs (4, 34–36). Iccosomes not bound by B cells are denoted A_f. B cells conjugated to FDCs through iccosomes are denoted B_a. During their interacting time, B_a cells are Ag-signaled and unbind, with rate a₁, as stimulated cells, B_e, carrying with them and processing some Ag. B_e cells can subsequently present Ag-derived peptides to Tfh cells, here denoted T. B_e cells interact in an Ag-specific fashion with T cells, with rate c₂, and form a B cell-T cell conjugate denoted T_b. During their conjugation B cells and T cells are assumed to activate each other. Cells in T_b conjugates detach from each other with rate a₂ as activated B and T cells, respectively, B_d and T_d, and they become new free GC B cells (B) and Tfh cells with rates p₁ and p₂, respectively. B and B_e cells die with constant rate d_b and T cells die with constant rate d_t.

FIGURE 1

Figure 1. Diagrammatic representation of the GC model core. This conceptual framework explicitly takes into account the main experimentally established interactions between FDCs, GC B cells, and Tfh cells (5, 33). Unbound GC B cells (B) interact with rate c₁ with free Ag on FDCs (A_f) (4, 34–36). B cells conjugated to FDCs through Ag (B_a) are Ag-signaled. They unbind, with rate a₁, as stimulated cells (B_e) possibly picking some Ag. B_e cells can subsequently interact with unbound Tfh cells (T), forming a B cell-Tfh cell conjugate (T_b), and present Ag-derived peptides. Cells in T_b conjugates activate each other and detach, with rate a₂, as activated B_d and T_d cells, which divide with rates p₁ and p₂, respectively, to become new free B and Tfh cells. B and B_e cells die with constant rate d_b and T cells die with constant rate d_t.

With the above indicated notation, and following a continuum approach, the population dynamics of this model core can be described by the following system of ordinary differential equations (ODE):

$\begin{array}{ll}d{A}_f/dt=-c_{1}B{A}_f+a_1{B}_a& (1.1)\end{array}$

$\begin{array}{ll}dB/dt=-c_{1}B{A}_f+p_1(1+{\alpha }_B)B_d-d_{b}B& (1.2)\end{array}$

$\begin{array}{ll}d{B}_a/dt=c_{1}B{A}_f-a_1{B}_a& (1.3)\end{array}$

$\begin{array}{ll}d{B}_e/dt=-c_2{B}_{e}T+a_1{B}_a-d_b{B}_e& (1.4)\end{array}$

$\begin{array}{ll}dT/dt=-c_2{B}_{e}T+p_2(1+{\alpha }_T)T_d-d_{t}T& (1.5)\end{array}$

$\begin{array}{ll}d{T}_b/dt=c_2{B}_{e}T-a_2{T}_b& (1.6)\end{array}$

$\begin{array}{ll}d{T}_d/dt=a_2{T}_b-p_2{T}_d& (1.7)\end{array}$

$\begin{array}{ll}d{B}_d/dt=a_2{T}_b-p_1{B}_d& (1.8)\end{array}$

where (1 + α_B) and (1 + α_T) are cell-density factors explained below.

All cell variables are expressed as cell units per GC, while A_f is given in units per GC of Ag complexes on FDCs' membrane (see Supplementary Table 1). For simplicity, Ag-containing A_f units are assumed to interact with B cells following on average a 1:1 stoichiometry. Therefore, B_a represents both the number of Ag-bound B cells and the number of B-cell bound Ag units. Moreover, at any given time, the maximum number of B cells that can interact with Ag in a GC equals the number of Ag-units in that GC. B_d and T_d cells can go to n division rounds before becoming B and T cells. Initially, we consider n = 1. This condition, however, will be relaxed later for B_d cells. Depending on the availability of limiting resources specific for each cell type, the average number of daughter cells generated by B_d and T_d cells can be lower than 2ⁿ and 2, respectively. More specifically, following (21, 37), the effective numbers of daughter cells are assumed to be given by cell-density dependent parameters (1 + α_B(t)) and (1 + α_T(t)), where the time-dependent parameters α_B and α_T are defined as:

$\begin{array}{ll}{\alpha }_B=(2^n-1)\times r_B\text{with}{r}_B=\frac{K_b}{K_b+B_{\text{T}}},& (1.9)\end{array}$

and

${\alpha }_T=\frac{K_t}{K_t+T_{\text{T}}},$

where B_T and T_T are, respectively, total B and T cells per GC (for instance, in the model core B_T = B + B_a + B_e + T_b + B_d and T_T = T + T_b + T_d), and K_b and K_t are limited resources-related parameters.

The three different models are introduced next. The variables and parameters of each model are described in Supplementary Tables 1, 2, respectively.

2.1.2. Model 1

This model is defined by the model core equations, except that Equation (1.1) is modified as follows to include Ag consumption by B cells (see Figure 2A):

$\begin{array}{ll}d{A}_f/dt=-c_{1}B{A}_f+\delta a_1{B}_a& (1.1\text{a})\end{array}$

where 0 < δ < 1 is the fraction of Ag on FDCs bound by a B cell that remains after the B cell detaches from a FDC.

FIGURE 2

Figure 2. Diagrammatic representation of the three alternative models of the GC reaction. (A) Model 1, it assumes that Ag is consumed proportionally to the number of B cells binding Ag on FDCs. (B) Model 2, here FDCs are assumed to change from initial phenotype F to phenotype F_m due to interactions with Ag-specific B cells. B cells interacting with Ag on F_m cells (i.e., A_Fm) differentiate into plasma or memory B cells (B_m), which exit the GC (indicated by the label “output”). (C) Model 3, here it is assumed that Tfh cells change phenotype, so that B cells interacting with those mature Tfh cells differentiate into plasma or memory B cells (B_m), exiting the GC (“output”). In both models, 2 and 3, Ag consumption by B cells is assumed not to lead to significant Ag depletion.

2.1.3. Model 2

The conceptual scheme of this model is depicted in Figure 2B. This model assumes that FDCs differentiate from an initial stage, denoted F, into a mature stage, named F_m, with a rate determined by the intensity of their Ag-dependent interactions with B cells. B cells bound to Ag on this mature fraction of FDCs are denoted B_am and are assumed to be induced to differentiate into memory or long-lived plasma B cells (B_m) which exit the GC as output cells. In addition, it is assumed that Ag is not significantly depleted by B-cell Ag uptake during the GC reaction. It is also assumed that F(t) + F_m(t) = F(0), at any time t. Finally, total Ag (free plus bound) is assumed to be partitioned between F and F_m cells proportionally to the relative amount of those cells within the total FDCs. That is, the total amount of Ag (free plus bound) on F cells is assumed equal to the total amount of Ag on F_m cells. Denoting A_f and A_Fm the amount of free Ag on F and F_m cells, respectively, the above assumption amounts to $A_f+B_a=A_f(0)\frac{F}{F(0)}$, and $A_{F_m}+B_{am}=A_f(0)\left(\frac{F(0)-F}{F(0)}\right)$. The corresponding ODE system is as follows:

$\begin{array}{ll}d{A}_f/dt=-m_1{A}_f-c_{1}B{A}_f+a_1{B}_a& (2.1)\end{array}$

$\begin{array}{ll}dB/dt=-c_{1}B{A}_f+p_1(1+{\alpha }_B)B_d-d_{b}B& (2.2)\end{array}$

$\begin{array}{ll}d{B}_a/dt=-m_1{B}_a+c_{1}B{A}_f-a_1{B}_a& (2.3)\end{array}$

$\begin{array}{ll}d{B}_e/dt=-c_2{B}_{e}T+a_1{B}_a-d_b{B}_e& (2.4)\end{array}$

$\begin{array}{ll}dT/dt=-c_2{B}_{e}T+p_2(1+{\alpha }_T)T_d-d_{t}T& (2.5)\end{array}$

$\begin{array}{ll}d{T}_b/dt=c_2{B}_{e}T-a_2{T}_b& (2.6)\end{array}$

$\begin{array}{ll}d{T}_d/dt=a_2{T}_b-p_2{T}_d& (2.7)\end{array}$

$\begin{array}{ll}d{B}_d/dt=a_2{T}_b-p_1{B}_d& (2.8)\end{array}$

$\begin{array}{ll}d{B}_{am}/dt=m_1{B}_a+c_{1}B{(A}_f(0)-A_f-B_a-B_{am})-a_{1m}{B}_{am}& (2.9)\end{array}$

$\begin{array}{ll}d{B}_m/dt=a_{1m}{B}_{am}-d_{bm}{B}_m& (2.10)\end{array}$

$\begin{array}{ll}dF/dt=-m_{1}F& (2.11)\end{array}$

where α_B and α_T are as in the model core, but with B_T = B + B_m + B_a + B_am + B_e + T_b + B_d and T_T = T + T_b + T_d. The rate of FDC differentiation to the mature stage F_m, denoted m₁, is assumed to be a time-dependent parameter proportional to the average amount of B-cell binding events per F cell at any time, that is, $m_1=\frac{{\mu }_1{B}_a}{B_a+A_f}\left(={\mu }_1\frac{F(0)}{A_f(0)}\frac{B_a}{F}\right)$.

2.1.4. Model 3

This model assumes that Tfh cells differentiate into a mature stage (denoted T_m) with a rate depending on the intensity of their Ag-dependent interactions with B cells. T_m cells are assumed to induce differentiation of GC B cells into memory or long-lived plasma B cells. B cells conjugated to T_m cells are denoted T_bm. Cells in T_bm conjugates detach from each other as activated B_dm and T_dm cells and B_dm cells are assumed to divide and differentiate into memory or long-lived plasma B cells (B_m) which exit the GC as output cells, while T_dm cells divide and become T_m cells. As in model 2, it assumes that Ag is not significantly depleted by B-cell Ag uptake during the GC reaction. The conceptual model is depicted in Figure 2C and the corresponding ODE system is as follows:

$\begin{array}{ll}d{A}_f/dt=-c_{1}B{A}_f+a_1{B}_a& (3.1)\end{array}$

$\begin{array}{ll}dB/dt=-d_{b}B-c_{1}B{A}_f+p_1(1+{\alpha }_B)B_d& (3.2)\end{array}$

$\begin{array}{ll}d{B}_m/dt=-d_{bm}{B}_m+p_{1m}(1+{\alpha }_B)B_{dm}& (3.3)\end{array}$

$\begin{array}{ll}d{B}_a/dt=c_{1}B{A}_f-a_1{B}_a& (3.4)\end{array}$

$\begin{array}{ll}d{B}_e/dt=-c_2{B}_e(T+T_m)+a_1{B}_a-d_b{B}_e& (3.5)\end{array}$

$\begin{array}{ll}dT/dt=-c_2{B}_{e}T+p_2(1-m_2)(1+{\alpha }_T)T_d-d_{t}T& (3.6)\end{array}$

$\begin{array}{ll}d{T}_b/dt=c_2{B}_{e}T-a_2{T}_b& (3.7)\end{array}$

$\begin{array}{ll}d{T}_d/dt=a_2{T}_b-p_2{T}_d& (3.8)\end{array}$

$\begin{array}{ll}d{B}_d/dt=a_2{T}_b-p_1{B}_d& (3.9)\end{array}$

$\begin{array}{l}d{T}_m/dt=-c_2{B}_e{T}_m+p_2{m}_2(1+{\alpha }_T)T_d\\ +p_{2m}(1+{\alpha }_T)T_{dm}-d_t{T}_m& (3.10)\end{array}$

$\begin{array}{ll}d{T}_{bm}/dt=c_2{B}_e{T}_m-a_2{T}_{bm}& (3.11)\end{array}$

$\begin{array}{ll}d{T}_{dm}/dt=a_2{T}_{bm}-p_{2m}{T}_{dm}& (3.12)\end{array}$

$\begin{array}{ll}d{B}_{dm}/dt=a_2{T}_{bm}-p_{1m}{B}_{dm}& (3.13)\end{array}$

where α_B and α_T are as in the model core, but with B_T = B + B_a + B_e + T_b + B_d + B_m + T_bm + B_dm, and T_T = T + T_b + T_d + T_m + T_bm + T_dm; and m₂ is the rate of Tfh cell differentiation to the mature stage mTfh. We assume this parameter is proportional to the fraction of T cells engaged in interactions with B cells, that is $m_2=\frac{{\mu }_2{T}_b}{T+T_b+T_d}$.

In all the three models, different subsets of B cells can be identified with either centrocytes (e.g., B and Ba) or centroblasts (e.g., Bd). Thus, these models reflect the diversity of classical GC B cell phenotypes, and hence do not oversimplify the GC reaction.

2.1.5. Parameter Reference Values

For the common core parameters and model-specific parameters a reference set of values was defined, based on experimental estimations for a majority of them, and on biologically reasonable, educated guess in other cases (see Supplementary Table 2). The particular case, in model 1, of the fraction δ of non depleted antigen per B_a cell deserves some detailed account. Recently, using coverslips covered with plasma membrane sheets decorated with antigens, it has been shown that high-affinity antigen-specific B cells spread extensively on the membrane and then contract, unbind and mechanically extract antigen (38). The amount of antigen extracted per B cell was estimated to be 20–80% of the antigen within the area of 50–100 μm² covered by the B cell spread when the antigen density was 50 molecules/μm². This amounts to 500–4,000 molecules extracted per B cell. Assuming that B cells cover, in FDC-B cell interactions, a dendrite area that is at most 1/10 of that covered on a planar membrane, it is expected that a high-affinity B cell in GCs can extract, upon unbinding, about 50–400 antigen molecules. On the other hand, it has been estimated, in anti-OVA immune responses from mice transferred with high-affinity, anti-OVA transgenic B cells, that in GCs up to ~ 15 % of transgenic B cells extract antigen after unbinding (35). Hence, GC B cells extract at most (0.15 × 50)–(0.15 × 400) ≈ 8-60 antigen molecules. Finally, assuming that each FDC is a depot of a minimum of 10³ antigen molecules in form of immunocomplexes, and that there are about 300 FDCs per GC (Supplementary Table 1), this implies that there are at least 3 × 10⁵ antigen molecules per GC. Therefore, based on all the above, it can be estimated that the fraction of antigen depleted per B_a (antigen-bound B cell) is in the range [8, 60]/(3 × 10⁵) = [2.5 × 10⁻⁵, 2 × 10⁻⁴]. That is, the fraction not depleted is at least δ = 0.9998. Relaxing the assumptions of the area covered by B cells in FDC-B cell interactions and of the number of FDCs per GC, and considering instead a 10-fold higher area covered by B cells and a 3-fold lower number of FDCs per GC, makes the fraction of Ag not depleted per B_a cell to be at least δ = 0.994. For model 1 we take the even more conservative value of δ = 0.99.

2.2. Model Analysis

The fact that GCs have a characteristic biphasic time evolution makes it possible to characterize their evolution by the maximum size attained in the growth phase and by the time taken to attain that size. Since this applies also to each of the variables in the three models (with the exception of A_f in model 1 and F in model 2) we took advantage of it and characterized the models' behavior by two types of quantities: (i) the Peak value of each variable (maximum or minimum), P_x (with x being any model variable), and (ii) the critical Time, T_x, that is, the time at which this peak is attained. The cases of variables A_f in model 1, and F in model 2 are special because these variables have a monotonically decreasing sigmoid behavior. A simple way of dealing with such cases is to characterize them by two quantities related to their (main) inflection point: the slope at the inflection point and the time at which it occurs. If x is such a variable, the said slope and time are precisely the peak value and critical time of its derivative dx/dt. Therefore, in the two special cases above mentioned we analyze the peak and critical time of dA_f/dt for model 1, and dF/dt for model 2. Each model behavior thus characterized was analyzed with respect to parameter sensitivity (how much a change in a parameter value affects the outputs—peaks and critical times—of the model) and synergy (the dependency of the sensitivities of one parameter on another parameter). We used two alternative, independent methods to calculate the sensitivities and synergies to ensure that convergence to the same solutions occurs in different approaches (for details see Supplementary Methods). We excluded parameter n from the sensitivity analysis because this is a discrete parameter and hence a conventional sensitivity analysis makes no sense for it. To circumvent this problem we have compared the impact of n on the GC dynamics for n = 1, 3, 5 in the three models. For calculational convenience, the equations of model 2 were solved using an equivalent system with the variable A_ft = A_f(0) − B_a − B_am (which is the total unbound antigen on FDCs) instead of A_f.

2.3. Software and Calculations

In order to ensure consistency in the analysis of the three models, we initially numerically solved the corresponding ODE systems using both the approximate and the analytical methods as well as two different ODE solvers, including the NDSolve function of Mathematica 9 and 11 (Wolfram Research, Inc.) run on both a Mac Pro and a MacBook Pro computer (for both the approximate and the analytical methods), and our own solver coded in C and run on an HP Linux machine (for the analytical method). The intermediary linear systems of first order ODEs with non constant coefficients resulting from the analytical method were numerically solved using also the NDSolve function of Mathematica and our own C solver. Once checked that the different methods and solvers gave coincident results, we resorted to the analytical method using Mathematica.

2.4. Global Sensitivity and Synergy Analysis: Latin Hypercube Sampling (LHS)

In order to perform a global analysis, a relatively wide range of values for each parameter was defined, ranging from half to double its reference value (see Supplementary Table 2). Parameter values were selected within those ranges using the LHS technique (39). This method allows reducing considerably the size of parameter space sampling, and hence the computational expense, while giving results with reasonably precision (40). Briefly, to generate m samples from the parameter space, the range to be explored for each parameter was divided in m intervals of equal size (in the absence of more detailed knowledge, a uniform probability distribution was assumed for parameter values within each parameter range). Within each interval one value was randomly obtained with the RandomReal function included in Mathematica. The m values thus obtained for a given parameter were randomly paired with the m values of a second parameter, and the resulting m pairs of values were randomly paired with the m values of a third parameter, and so on. In this way a collection of m sets of parameter values were obtained. The results reported below were obtained with m = 20 samples for each of the three models.

2.5. Global System Sensitivity and Synergy

We define the system sensitivity in a way that differs from previous authors (41). Thus, rather than defining it as the 2-norm of all individual sensitivities to a given parameter, we define it here as the geometric mean of the individual sensitivities. In contrast to the 2-norm (which biases the estimation toward the highest individual sensitivity) the geometric mean is a more conservative quantity in that the weight of individual values is very moderate. We extend this concept to calculate also the system synergy.

Let $S_P^i=(S_{P_B}^i,S_{P_{B_a}}^i,S_{P_{B_e}}^i,S_{P_{T_b}}^i,S_{P_{B_d}}^i,S_{P_T}^i,S_{P_{T_d}}^i,S_{P_{A_f}}^i)$ be the local, relative Peak sensitivities of model 1 variables to parameter p_i. We define the Peak system sensitivity to parameter p_i, SSi_P, as the geometric mean of the absolute values of $S_P^i$:

$S{S}_P^i={\left({\displaystyle \prod _{k=1}^n}|S_P^i(k)|\right)}^{\frac{1}{n}},$

where n is the dimension of $S_P^i$, that is, the number of variables in model 1. Similarly, the critical Time system sensitivity to parameter p_i is defined as:

$S{S}_T^i={\left({\displaystyle \prod _{k=1}^n}|S_T^i(k)|\right)}^{\frac{1}{n}}.$

And the system sensitivities of models 2 and 3 are defined in the same fashion.

The geometric mean was also used to define the system synergies of the three models:

$S{R}_P^{ij}={\left({\displaystyle \prod _{k=1}^n}|R_P^{ij}(k)|\right)}^{\frac{1}{n}}\text{ and }S{R}_T^{ij}={\left({\displaystyle \prod _{k=1}^n}|R_T^{ij}(k)|\right)}^{\frac{1}{n}}.$

By using the m sets of parameter values obtained, for each model, with the LHS method, we performed a global system sensitivity and synergy analysis.

3. Results

GCs display an initial expansion phase, dominated by B cell proliferation and to a lesser extent T cell growth, followed by a latter contraction phase, the reason for which is still unclear. Several processes probably contribute to it: decrease of availability of unbound Ag epitopes, B cell apoptosis, and egression of B cells differentiating to long-lived plasma cells (42, 43). Three types of cell interactions are required for a GC to develop, namely, Ag-specific B cells interacting with Ag complexes in iccosomes displayed on FDC membranes (44, 45), B cells interacting with FDCs (45–47), and B cells interacting with Tfh cells (48, 49). Our aim here is to understand what are the cell interactions to which the GC dynamics is robust, and what are the interactions to which it is more sensitive, i.e., the interactions that most likely drive the GC dynamics in a tunable way. To this end, as explained above, we have built a model core representing the basic cell interactions described so far to take place in GCs, after which we developed three models by adding specific, experimentally established processes that can determine the rise-and-fall dynamics of GCs. To consider a model as being minimally realistic we require it to be consistent with the following quantitative features that characterize a typical GC reaction in a murine primary immune response to a protein Ag: (1) a rise-and-fall dynamics in which the peak of the GC reaction involve up to 15,000 B lymphocytes (9); (2) a time to the peak of 8–15 days after immunization (9, 50, 51); (3) a life-span of the GC reaction of up to 4 weeks (9, 50, 51); (4) an amount of Tfh cells at the peak of the response of 5–20% of total lymphocytes (50–53).

We first describe results characterizing the dynamical properties of the different models. After that, we present our local and global system sensitivity and synergy analysis for each of the models.

3.1. Impact of B Cell Initial Conditions and Number of B Cell Division Cycles on the Behavior of the Models: Concentration Kinetics, Sensitivity, and Synergy

Recently, a theoretical reassessment of the classical estimation of the number of B cells seeding GCs predicted that, contrary to the dominant belief (15), GCs are seeded by 50–200 B cells (54). This theoretical prediction has been since then confirmed using a direct experimental approach (55). On the other hand, using an experimental system in which B cells are artificially forced to present to Tfh cells levels of membrane peptide-MHCII (p-MHCII) well above those of mid/high affinity B cells, it has been recently estimated that the number of cellular divisions by GC centroblasts range from 1 to 6 before they return to the centrocyte stage, with a majority of cells following three divisions (56). In contrast, under physiologic conditions mid/high affinity B cells follow on average 2 cell divisions (56). In order to explore the relevance of both estimations, we compared the concentration kinetics as well as the parameter sensitivities and synergies obtained in simulations under three different conditions: case 1, B(0) = 10 cells and ${\alpha }_B(t)=(2^n-1)\times r_B$ (see Equation 1.9) with n = 1 (corresponding to a maximum of one cell division by B_d cells); case 2, B(0) = 100 cells and n = 1; and case 3, B(0) = 100 and n = 3 (corresponding to a maximum of three B_d cell divisions). This kind of analysis was performed for each of the three models, and for each of 20 different sets of parameters randomly generated using the LHS method. The models' simulations performed with each parameter set are referred to in the following as in silico experiments.

For every model, in each of the three cases the 20 in silico experiments, corresponding to the 20 different parameter sets, gave variable results for concentration kinetics (with model 3 displaying the largest variability) but similar sensitivity results. However, for every parameter set cases 1 and 2 gave, by and large, nearly indistinguishable results for concentration kinetics and for sensitivities of Peaks and critical Times, and similar or very similar results for synergies of Peaks and critical Times. Representative examples of the results for one particular parameter set are shown in Figures 3–7.

FIGURE 3

Figure 3. Reference dynamics of model 1. (A) Dynamics of individual compartments. The second row shows in detail the kinetics of those variables reaching low values. (B) Global dynamics. The dynamics of total B cells (B_T), total T cells (T_T), and total lymphocytes (B_T + T_T) are shown. (C) Percent of total T cells out of total lymphocytes. Parameter values used: δ = 0.99; c₁ = 17.33 day⁻¹(molec/GC)⁻¹; c₂ = 16.18 day⁻¹(cell/GC)⁻¹; p₁ = 1.97 day⁻¹; d_b = 2.68 day⁻¹; a₁ = 33.11 day⁻¹; a₂ = 19.92 day⁻¹; p₂ = 1.78 day⁻¹; d_t = 0.22 day⁻¹; K_b = 8285.60 cell/GC; K_t = 383.98 cell/GC.

FIGURE 4

Figure 4. Reference dynamics of model 2. (A) Dynamics of individual compartments. The second row shows in detail the kinetics of those variables reaching low values. (B) Global dynamics. The dynamics of total B cells (B_T), total T cells (T_T), and total lymphocytes (B_T + T_T) are shown. (C) Percent of total T cells out of total lymphocytes. Parameter values used: μ₁ = 0.22; c₁ = 19.04 day⁻¹(molec/GC)⁻¹; c₂ = 30.07 day⁻¹(cell/GC)⁻¹; p₁ = 1.17 day⁻¹; d_b = 1.17 day⁻¹; a₁ = 17.90 day⁻¹; a_1m = 11.94 day⁻¹; a₂ = 13.77 day⁻¹; d_bm = 3.11 day⁻¹; p₂ = 1.04 day⁻¹; d_t = 0.15 day⁻¹; K_b = 9782.02 cell/GC; K_t = 598.82 cell/GC.

FIGURE 5

Figure 5. Reference dynamics of model 3. (A) Dynamics of individual compartments. The second row shows in detail the kinetics of those variables reaching low values. (B) Global dynamics. The dynamics of total B cells (B_T), total T cells (T_T), and total lymphocytes (B_T + T_T) are shown. (C) Percent of total T cells out of total lymphocytes. Parameter values used: μ₂ = 0.10; c₁ = 19.13 day⁻¹(molec/GC)⁻¹; c₂ = 16.77 day⁻¹(cell/GC)⁻¹; p₁ = 2.25 day⁻¹; p_1m = 1.16 day⁻¹; d_b = 2.03 day⁻¹; a₁ = 19.43 day⁻¹; a₂ = 14.29 day⁻¹; d_bm = 1.20 day⁻¹; p₂ = 2.70 day⁻¹; p_2m = 1.63 day⁻¹; d_t = 0.39 day⁻¹; K_b = 12632.48 cell/GC; K_t = 315.18 cell/GC.

FIGURE 6

Figure 6. Local parameter sensitivity analysis of the three models. Histograms of the relative sensitivity of the peaks and critical times of all model variables with respect to each parameter are shown. Upper two panels correspond to Model 1, middle panels to Model 2, and lower panels to Model 3. Left side panels (case 1), center panels (case 2), and right side panels (case 3) of the three models show the results, respectively, for B(0) = 10 and n = 1, B(0) = 100 and n = 1, and B(0) = 100 and n = 3 (see text for details). Representative results from a single set of parameter values are shown. Parameter values used were the same as in Figure 3 for Model 1, Figure 4 for Model 2, and Figure 5 for Model 3.

FIGURE 7

Figure 7. Local parameter synergy analysis of the three models. Histograms show representative results of the relative synergy between parameter p₂ and all other parameters with respect to the peaks and critical times of all model variables. Upper two panels correspond to Model 1, middle panels to Model 2, and lower panels to Model 3. Left side panels (case 1), center panels (case 2), and right side panels (case 3) of the three models show the results, respectively, for B(0) = 10 and n = 1, B(0) = 100 and n = 1, and B(0) = 100 and n = 3 (see text for details). Parameter values used were the same as in Figure 3 for Model 1, Figure 4 for Model 2, and Figure 5 for Model 3.

In contrast, for each parameter set cases 2 and 3 gave very different results with respect to concentration kinetics in all the three models (see Figures 3–5 panels a2 vs. panels a3). In addition, although the sensitivities and synergies of Peaks and critical Times were very similar in cases 2 and 3 of model 1, they were at best moderately similar in these two cases of model 3. Moreover, in model 2 only the sensitivities of critical Times were relatively similar in cases 2 and 3 (Figures 6, 7).

In all the three models, when B cells proliferate with a single replication round (n = 1) the GC total T cell population grows to levels comparable to those of total B cells, with slightly delayed kinetics, irrespective of the initial value of B cells (Figures 3–5 panels b1, b2, c1, c2). However, for higher replication rounds (n = 3) the fraction of total T cells decreases considerably (Figures 3–5 panel c3) approaching estimated values in available experimental data (50–53). In view of the impact of n when increased from n = 1 to n = 3 on the fraction levels of total T cells and peak levels of all cell compartments in the three models, as well as on the critical times in model 3, we performed an additional set of in silico experiments in each model with n = 5 for the same 20 sets of parameter values used before. In this new case, in model 3 the peaks and critical times are considerably larger than for n = 3, particularly with respect to total B cells, while the sensitivities and synergies for the other parameters are comparable to those obtained with n = 3. In models 1 and 2 the critical times are not markedly changed, but the peaks increase 1.5- to 2-fold compared to n = 3.

These results indicate, first, that the particular value for the initial condition B(0) does not impinge on the GC dynamics in any of the three models, and second, that under a moderate level of replication rounds (n = 3) the proportion of total T cells is by far more realistic than under a single replication round (n = 1). In summary, this result uncovered the fact that for T cells to be at most 20% of total GC lymphocytes during most of the GC reaction, the maximum proliferation expansion of dividing B cells (n) should be higher than that of dividing T cells. However, if n is too high (for instance, n = 5) the cells' kinetics in all models display an excessive growth for all or most of the 20 in silico experiments, which is contrary to the above mentioned experimental observation, according to which under physiologic conditions dark zone GC B cells can initiate an average of 2 cell divisions (56). This, in addition to the important differences in sensitivities and synergies between case 3, on the one hand, and cases 1 and 2, on the other hand, for models 2 and 3, prompted us to use in the remaining of this work the conditions of case 3, unless otherwise stated.

3.2. Dynamical Properties

The three GC models were analyzed for their equilibrium points and their ability to reasonably reproduce the global kinetics of GC B and T cells. As expected, the three models exhibit a single fixed point at t = ∞ in which all variables, except A_f in models 2 and 3, and F in model 2, go to zero. In model 2, variables F and A_f decrease in a sigmoidal fashion and go to a positive asymptotic value which depends on the parameter values. In model 3, variable A_f attains a minimum and asymptotically returns to its maximum (initial) value. All other cell variables in the three models attain a maximum and it can be shown that then decrease exponentially to zero. The main reason why GCs in models 2 and 3 ultimately decline is because of the implemented mechanism in each model (FDC maturation in model 2 and Tfh maturation in model 3). Each mechanism leads to an accumulation of either mature F_m cells (model 2) or T_m cells (model 3). Conjugates of these cells with B cells induce them to become B_m cells, which exit the GC as memory B cells or long-lived plasma cells (output cells). Therefore, initially most B cells are induced to proliferate and fewer cells to differentiate, but with time less and less B cells are induced to proliferate and more cells to differentiate and egress from the GC. At some point in time the number of F_m or T_m cells has increased to such a level that the number of egressing, differentiated B cells (that is, B_m cells) equals the number of not differentiated B cells, so that the net variation of B cells is zero. After that time, the number of F_m or T_m cells has increased beyond that critical level so that the fraction of B_m cells egressing from GCs is larger than the fraction of remaining B cells. This shrinks the size of the GC B cell population and as a consequence the size of the GC T cell population.

This general behavior of the models indicates that their biphasic dynamics is independent of, and hence is not due to, the added carrying capacity. We confirmed this is the case by making simulations of the three models with α_B = α_T = 1. As expected, we obtained a biphasic behavior, but with very high peaks (>10⁵ total cells) and very large critical times (> 100 days). This result also justifies the need of considering a carrying capacity. Whether this is required for both B and T cells or for only B or T cells is analyzed below in section 3.3.

All models share the same characteristic kinetics of B_e cells: this is always the first cell population to peak, then it decreases to a value near zero within a relatively short period of time. This triggers a cascade of events in the kinetics of the other cell populations. In model 1, B_e and B_a peak when A_f is near zero, while at this time free B cells start to increase exponentially because they have no free Ag to combine with. Nevertheless, irrespective of whether Ag is consumed or not, in all models, B_a is limited by the amount of total Ag, and therefore when B_a is close to its maximum the combined rates of B_e death plus conjugation to T cells (T_b formation) starts to dominate the kinetics of B_e cells causing them to decline. As a consequence, B_e cells engage less T cells (form less T_b) and so free T cells increase slightly, which makes B_e cells to further decline. When B_e cells fall below a certain level, most free T cells (and free T_m cells in model 3) cannot engage with them and, hence, they become visible, with an exponential increase. Consequently, T_b conjugates decrease. In turn, and because of this, B_d and T_d start also to decline, which leads to a decrease in the generation of new free B and T cells. In models 2 and 3, an increasing fraction of B_e cells become output B cells (B_m). However, they follow different pathways in each model, thus contributing differently to the kinetics of the other cell populations.

In model 1, for the twenty different in silico experiments the critical times of total B and T cells (times at which the peaks are attained) range from 9 to 32 days. By and large, the most frequent critical time values were within 10 and 19 days. Total B cells peaked at values ranging from 10, 000 to 36, 000 cells, with the most frequent values ranging between 10, 000 and 18, 000 cells. Last, in all but one in silico experiment, at the peak of the GC reaction total T cells constituted at least 20% of total lymphocytes, with 16 parameter sets leading to values ranging between 23 and 37%. The critical times and peaks correlated poorly, if at all (R² = 0.15), and negatively. In this model, in order to get a total B cell dynamics compatible with experimental observations, Ag should be degraded with a kinetics that implies a 10-fold reduction in about 1–3 weeks. However, as explained in detail in section 2.1.5, this requirement is in disagreement with semi-quantitative observations of the decay of Ag in GCs (57, 58) and of the GC kinetics for a wide range of Ag deposition on the membrane of FDCs (59). Nevertheless, it must be pointed out that to date it has not been possible to quantify in vivo Ag decay by its uptake by GC B cells (19).

In model 2, the obtained critical times of total B and T cells lay between 8 to 27 days, with the most frequent values ranging from 10 to 14 days. Total B cells peaked at values ranging from 10, 000 to 54, 000 cells, with the most frequent values, by and large, in the range of 20, 000 to 40, 000 cells. In general, in this model the peak value of B_m cells is much higher than that of any other lymphocyte population, so that in most of the performed in silico experiments they constitute 25–50% of the total B cell population. Considering that B_m are cells exiting from GCs, the above result for peak values indicates that in model 2 all but three parameter conditions lead to a peak of total resident GC B cells (i.e., excluding egressing B_m cells) ranging from 3, 000 to 15, 000 cells. In addition, in more than half of the performed in silico experiments total T cells are between 15 and 20% at the critical time. The critical times and peaks were uncorrelated in model 2.

In contrast to models 1 and 2, in model 3 the critical time values of total B and T cells spanned a much broader range, varying from 6 to 105 days, with eight parameter sets leading to values between 20 and 40 days and seven parameter sets leading to values between 50 and 105 days. The peak values of total resident B cells (i.e., excluding egressing B_m cells) also ranged broadly, from 2,000 to 59,000 cells, with 12 out of 20 values varying between 25,000 and 59,000 cells. With respect to the fraction of total T cells at the height of the GC reaction, in all of the 20 in silico experiments they constituted at least 21% of total GC lymphocytes, with 18 out of 20 parameter sets leading to values ranging between 25 and 40%. In this model, several parameter sets lead to two types of total B and T cell dynamics that are far from the conventional rise-and-fall GC lymphocyte dynamics: In the first type, after reaching a peak, B and T cells remain very close to that value, declining very slowly for a prolonged period of time, after which they decay considerably faster. In the second type, after reaching a peak, B cells follow one or more small oscillations and then decline very slowly to zero, and T cells reach a peak and then simply decline very slowly to zero. The critical times and peaks were weakly correlated in model 3 (R² = 0.23).

The availability of some experimental data on the kinetics of GC size (53) allows us to compare it with our theoretical estimations of GC size kinetics in cases 1, 2, and 3 of the three models (Figure 8). In that experimental work, GC size was measured as the average area occupied by individual GCs in a spleen section, while in the theoretical estimations we measured GC size as the average number of total GC B cells in 20 different in silico experiments. The theoretical estimations were performed at the same time points as the experimental measurements. As shown in Figures 8A,B, in cases 1 and 2 the early phase of the GC size kinetics of the three models is very similar to the experimental data until day 10 or 12, but then for models 1 and 3 the GC sizes keep substantially growing until day 21. In contrast, for model 2 the GC size reaches a peak at days 14 (case 1) and 12 (case 2). With respect to case 3 (Figure 8C), the GC size kinetics for the three models has again a similar early phase (until day 8), but then for model 1 the GC size keeps increasing faster than in models 2 and 3, and reaches a peak at day 16. In contrast, in model 3 the GC size increases steadily until day 21. In model 2 the GC size reaches a peak at day 12 with kinetics quite comparable to the experimental one (Figure 8D).

FIGURE 8

Figure 8. Comparison of theoretical and experimental estimations of the kinetics of average GC sizes. (A–C) Show the theoretical estimations obtained for three different conditions (cases 1–3, see below) and (D) shows the experimentally obtained kinetics of GC sizes. The theoretical estimations were performed at the same time points as the experimental ones. GC sizes are expressed either as number of total GC B cells (A–C) or as area of GC sections (in μm², D). (A–C) Correspond, respectively, to case 1 (B(0) = 10; n = 1), case 2 (B(0) = 100; n = 1), and case 3 (B(0) = 100; n = 3). Each of these panels shows the results obtained with the three models. Markers and bars in (A–C) correspond, respectively, to the mean and SEM of 20 in silico experiments. Bars in (D) correspond to the SD of the mean (solid circles). Numerical values of experimental data in (D) were obtained from Figure 2C in (53).

3.3. Global Sensitivity Analysis

In order to cover an ample hypervolume in parameter space and, thus, to perform a global sensitivity analysis, we followed a two-step approach. First, we defined a system sensitivity for each parameter as a way to quantify the overall Peak and critical Time sensitivities of the system variables for each parameter (see Supplementary Methods). And second, for each model we generated 20 different parameter sets using the LHS method, with each parameter value being randomly sampled within a fourfold range centered at the corresponding reference value (Supplementary Table 2).

The results obtained with this approach with respect to both the Peaks and the critical Times are displayed, for each of the three models, in Figure 9. There, parameters are ranked from higher to lower values based on the arithmetic mean of their system sensitivity. It is conspicuous that the variation range of system sensitivities in models 1 and 2 is much smaller for the critical Times than for the Peaks, and moderately smaller in model 3. In model 1, parameters can be classified in three groups: (1) p₂ and K_t, with highest sensitivities; (2) a₁, a₂, p₁, d_b, K_b, with medium to low sensitivities; and (3) c₁, c₂, d_t, with very low sensitivities. Similarly, in model 2 parameters classify in three major groups: (1) p₂, a₁, μ₁, K_t, yield highest sensitivities; (2) a_1m, a₂, K_b, p₁, d_b, with low sensitivities; and (3) c₁, c₂, d_t, d_bm, give very low sensitivities. Last, in model 3 parameters also classify in three major groups: (1) p₂, p_2m, yield highest sensitivities; (2) μ₂, K_t, a₂, K_b, p₁, d_bm, a₁, d_b, with medium to low sensitivities; and (3) c₁, c₂, d_t, p_1m, with very low sensitivities.

FIGURE 9

Figure 9. Global system sensitivity. Symbols in each bin correspond to the system sensitivities obtained in 20 simulations performed with different sets of parameter values obtained with the LHS method. Thick horizontal line in each bin represents the arithmetic mean in that bin. Parameter bins are sorted from higher to lower mean. In all simulations B(0) = 100 and n = 3.

Remarkably, parameter p₂ is ranked highest in all models with respect to both the Peaks and the critical Times.

3.4. Global Synergy Analysis

We analyzed the global synergy following the same approach as for the sensitivity. Thus, the parameters' system synergies were calculated as defined in Supplementary Methods for the same 20 parameter sets used in the analysis of the global sensitivities. For each pair of parameters the arithmetic mean of their 20 system synergy values was calculated. The results, described in Supplementary Results and shown in Supplementary Figures 1–3, add to those of the sensitivity analysis, indicating that parameter p₂ has in all models the highest synergy with respect to Peaks and critical Times.

4. Discussion

The dynamics of GC reactions differ widely for Ags of very different complexities. However, in all or most cases GCs are initiated ensuing the same order of processes as follows: (1) a couple of days after immunization with a protein-containing Ag, small foci of B and Th cells form adjacent to primary follicles; (2) afterwards, in a fraction of foci some Ag-activated B and Th cells migrate inside the adjacent primary follicle and initiate a GC (60). However, foci formation requires the synchronous presence of recirculating Ag-specific Th cells in a T-cell zone and of B cells in an adjacent follicle. For purified, relatively small proteins this only happens in a fraction of follicles due to the low frequency of epitope-specific lymphocytes before immunization and the relatively low number of epitopes presented by those Ags. Because in general the bigger and the more complex a proteinaceous Ag is, the larger the number of total and different epitopes it expresses, the probability of generating a foci, and hence of generating a GC, increases substantially with the size and complexity of the immunizing Ag. Moreover, it also increases the probability that new migrating B cells entering an ongoing GC are specific for some of the many present epitopes, and hence these B cells can be recruited, perturbing the dynamics of the ongoing GC. In addition, very large Ags like virus and virus-like particles can either reproduce and disseminate, or simply leak for a while as observed with iccosomes (34), allowing for a more or less repeated stimulation of new immune reactions in distant lymphoid tissue. This could explain why virus and virus-like particles often generate a chronic GC reaction, possibly with at least some long-lived GCs. In contrast, purified, mid-size proteins trigger in mice a typical acute, primary GC reaction. The modeling analysis performed here of the dynamics of an average GC is based on such acute, primary GC reactions.

The three different models of the GC dynamics implemented here included a number of processes involving different stages of FDCs, B cells, and Tfh cells. Often, some model processes determine more strongly than others the dynamics of the system's behavior. We characterized their relative impact by analyzing the sensitivity of every model variable to each parameter. For each of the 20 in silico experiments performed with every model, we also calculated the system sensitivity, a quantity used to summarize the sensitivity of all the system variables to a given parameter. Irrespective of the mechanism assumed to drive the rise-and-fall dynamics of GCs, the analysis of both the global sensitivity and synergy indicates very clearly that the GC dynamics is most sensitive to parameter p₂, and to a lesser extent to parameter K_t. Surprisingly, the GC dynamics in models 2 and 3 is only moderately sensitive to the model specific parameters μ₁ and μ₂. As an overview of the relative impact of all parameters on the GC dynamics in the three different models that were analyzed, we coarsely identified sensitivities and synergies as high, medium, and low (Table 1).

TABLE 1

Table 1. Overview of parameters impact on the GC dynamics.

One could expect that the outputs used to quantify a dynamical system sensitivity to parameters would be robust to parameters that are foreseen to vary during the system's operation, and sensitive to ones that experience little change. The maximal rates of B cell and Tfh cell proliferation, respectively, p₁ and p₂, are likely to vary importantly, increasing or decreasing at least four-fold (median cell division duration between 5 and 20 h) depending on local microenvironmental conditions. The parameters K_b and K_t, that quantify the effect on B cell and Tfh cell growth of general limiting resources, can also vary by local fluctuations of secreted growth factors and B-cell or T-cell specific metabolic regulators (61–63). Strikingly, the GC dynamics in the three models is robust to p₁ and K_b, but is highly sensitive to p₂ and to a lesser extent to K_t. This suggests that there are still undiscovered processes that regulate the GC dynamics by modulating the Tfh cell growth rate. One possibility, worth to be explored, is the possible contribution of the recently discovered follicular T regulatory cells (51, 64, 65) in the modulation of Tfh cell division during the GC dynamics.

With respect to GC dynamics, the analysis of the total B and T cell dynamics in the different models revealed that, except for a brief transient, the ratio of B and T cells seeding GCs [i.e., B(0)/T(0)] has virtually no effect on the kinetics of any individual system variable and hence does not account for the experimentally observed 5–20% of GC T cells, out of total GC lymphocytes, at the height of the reaction (50–53). However, differences between the maximum number of consecutive division rounds that a proliferating B and T cell may undergo can account for the experimentally observed percentage of GC T cells. In particular, when T-cell-activated B cells were allowed to undergo a maximum of 3 consecutive divisions vs. 1 division for activated T cells, the fraction of total T cells decreased substantially during most of the GC reaction in all three GC models (Figures 3–5 panels c1–c3). Nevertheless, a more detailed analysis of the 20 in silico experiments showed that model 2 is the one that best captures the typical GC dynamics in terms of T-cell fractions when compared with models 1 and 3, with up to 11 in silico experiments resulting in T-cell percentages below 20%, whereas the large majority of the in silico experiments with models 1 and 3 resulted in percentages above 23 and 25%, respectively.

But there are still other aspects of GC dynamics that models 1 and 3 struggle to reproduce. A typical GC reaction in a murine primary immune response takes up to 4 weeks with its height at about days 8–15 (9, 50, 51, 53). This typical GC B cell dynamics is hardly reproduced by model 3 because: (1) the peaks and critical times are highly dependent on parameter values, spanning a very broad range in the 20 in silico experiments, and (2) total B cells either grow to high peak levels (>25, 000 cells, and up to 59,000 cells, excluding egressing Bm cells) with high critical times (>20 days, and up to 105 days), or grow to reasonable peak levels (<10, 000 cells) in 1–2 weeks, but then decrease very slowly after following, at least in some in silico experiments, a couple of dampened oscillations. These dynamical characteristics are reflected in the very different GC size kinetics obtained with model 3 compared to an experimentally obtained one (Figure 8). In contrast, model 1 can reproduce in most in silico experiments the typical GC dynamics. However, this requires an FDC-bound B cell to consume, on average, more than 5% of the Ag displayed by that FDC, which as discussed above (see section 2.1.5), is highly unrealistic. Even an Ag consumption of 1% per FDC-bound B cell, as used here, leads to either high peak values (>15, 000 cells) or critical times too high (>15 days), with an average GC size kinetics that differ substantially from an experimental one (Figure 8). In addition, it was recently shown that Ag trapped on FDCs does not remain on the membrane all the time but is rapidly internalized, remaining intact within a non-degradative cycling compartment, and being displayed periodically on the FDC surface where it is accessible to Ag-specific B cells (36). This indicates that Ag immunocomplexes in GCs are partially protected from circulating Abs and from excessive Ag consumption by GC B cells, suggesting that its concentration changes even less than what we estimated above. For model 1, this implies a duration of the GC reaction in the order of months rather than weeks, indicating that this model is unable to satisfactorily explain the GC dynamics. In contrast, model 2 reproduces adequately the typical GC dynamics of total B cells (Figure 8), with seventeen out of twenty in silico experiments having peak values (excluding egressing Bm cells) between 3,000 and 15,000 cells, among which eleven had critical times between 8 and 13 days.

In summary, model 2 is the one that best captures the typical GC dynamics with respect to the three analyzed features: peak values of total B cells, critical times, and the fraction of total T cells, whereas models 1 and 3 fail to adequately reproduce these features. This does not rule out the mechanism of model 3 as an important contributor to the GC dynamics. Rather, it means that in its present form this mechanism by itself does not adequately explain GC dynamics and, therefore, if it were at work in GCs there should be at least one additional process, perhaps still to be uncovered, working in concert with this mechanism. Also, if future experiments disprove the FDC maturation mechanism of model 2 as a driver of the GC dynamics, it would switch the interest toward trying to uncover mechanisms complementary to the Tfh maturation mechanism of model 3 which, all together, could become an adequate description of the GC dynamics. It is worth mentioning that in model 3 we have assumed that all Tfh cell maturation results in T_m cells that are capable of giving the same quality of help to B cells. However, this is likely an oversimplification and may not be so. For instance, recently matured T_m cells may induce B cells to differentiate into memory B cells, while older T_m cells may induce them to differentiate into plasmablasts. Notably, in all three mechanisms, it is a T cell-centered parameter, p₂, that not only has the highest impact among all parameters but also a high impact per se on the dynamics of the Germinal Center. This general result suggests that the intensity of activation of Tfh cell division could be a natural target of therapies aimed at either potentiating or lessening GC reactions.

Data Availability

All datasets generated for this study are included in the manuscript/Supplementary Files.

Author Contributions

JF conceived and designed the research, performed the analysis, and wrote the manuscript. RG and BvH performed the initial phase of the research and critically revised the manuscript. EF designed the research, performed the analysis, and wrote the manuscript.

Conflict of Interest Statement

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.

Acknowledgments

The authors gratefully acknowledge Michal Or-Guil for reading and providing comments on the text. Partially supported by Grant POCTI/36413/1999 (Fundação para Ciência e Tecnologia, Portugal), and the European Union 7th Framework Program under grant agreements REGPOT-2012-2013.1-316265 (BIOCAPS) and PIRSES-GA-2012-317893 (INDOEUROPEAN-MATHDS), and Xunta de Galicia (Grupo de referencia competitiva, ED431C 2016041).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2019.02038/full#supplementary-material

References

1. Reif K, Ekland E, Ohl L, Nakano H, Lipp M, Forster R, et al. Balanced responsiveness to chemoattractants from adjacent zones determines B-cell position. Nature. (2002) 416:94–9. doi: 10.1038/416094a

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Berek C, Berger A, Apel M. Maturation of the immune response in germinal centers. Cell. (1991) 67:1121–9. doi: 10.1016/0092-8674(91)90289-B

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Liu YJ, Oldfield S, Maclennan I. Memory B cells in T cell-dependent antibody responses colonize the splenic marginal zones. Eur J Immunol. (1988) 18:355–62. doi: 10.1002/eji.1830180306

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Szakal AK, Kosco MH, Tew JG. Microanatomy of lymphoid tissue during humoral immune responses: structure function relationships. Annu Rev Immunol. (1989) 7:91–109. doi: 10.1146/annurev.iy.07.040189.000515

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Crotty S. Follicular helper CD4 T cells (TFH). Annu Rev Immunol. (2011) 29:621–63. doi: 10.1146/annurev-immunol-031210-101400

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Choi YS, Yang JA, Crotty S. Dynamic regulation of Bcl6 in follicular helper CD4 T (Tfh) cells. Curr Opin Immunol. (2013) 25:366–72. doi: 10.1016/j.coi.2013.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Kelsoe G. In situ studies of the germinal center reaction. Adv Immunol. (1995) 60:267–88. doi: 10.1016/S0065-2776(08)60587-8

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Camacho SA, Kosco-Vilbois MH, Berek C. The dynamic structure of the germinal center. Immunol Today. (1998) 19:511–4. doi: 10.1016/S0167-5699(98)01327-9

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Wittenbrink N, Weber TS, Klein A, Weiser AA, Zuschratter W, Sibila M, et al. Broad volume distributions indicate nonsynchronized growth and suggest sudden collapses of germinal center B cell populations. J Immunol. (2010) 184:1339–47. doi: 10.4049/jimmunol.0901040

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Calame KL, Lin KI, Tunyaplin C. Regulatory mechanisms that determine the development and function of plasma cells. Annu Rev Immunol. (2003) 21:205–30. doi: 10.1146/annurev.immunol.21.120601.141138

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Longo NS, Lipsky PE. Why do B cells mutate their immunoglobulin receptors? Trends Immunol. (2006) 27:374–80. doi: 10.1016/j.it.2006.06.007

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Moreira JS, Faro J. Re-evaluating the recycling hypothesis in the germinal centre. Immunol Cell Biol. (2006) 84:404–10. doi: 10.1111/j.1440-1711.2006.01443.x

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Anderson SM, Khalil A, Uduman M, Hershberg U, Louzoun Y, Haberman AM, et al. Taking advantage: high-affinity B cells in the germinal center have lower death rates, but similar rates of division, compared to low-affinity cells. J Immunol. (2009) 183:7314–25. doi: 10.4049/jimmunol.0902452

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Or-Guil M, Faro J. A major hindrance in antibody affinity maturation investigation: we never succeeded in falsifying the hypothesis of single-step selection. Front Immunol. (2014) 5:237. doi: 10.3389/fimmu.2014.00237

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Victora GD, Mesin L. Clonal and cellular dynamics in germinal centers. Curr Opin Immunol. (2014) 28C:90–6. doi: 10.1016/j.coi.2014.02.010

CrossRef Full Text | Google Scholar

16. Keşmir C, Boer RJD. A spatial model of germinal center reactions: cellular adhesion based sorting of B cells results in efficient affinity maturation. J Theor Biol. (2003) 222:9–22. doi: 10.1016/S0022-5193(03)00010-9

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Zhang J, Shakhnovich EI. Optimality of mutation and selection in germinal centers. PLoS Comput Biol. (2010) 6:e1000800. doi: 10.1371/journal.pcbi.1000800

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Amitai A, Mesin L, Victora GD, Kardar M, Chakraborty AK. A population dynamics model for clonal diversity in a germinal center. Front Microbiol. (2017) 8:1693. doi: 10.3389/fmicb.2017.01693

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Bannard O, Cyster JG. Germinal centers: programmed for affinity maturation and antibody diversification. Curr Opin Immunol. (2017) 45:21–30. doi: 10.1016/j.coi.2016.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Moreira JS, Faro J. Modelling two possible mechanisms for the regulation of the germinal center dynamics. J Immunol. (2006) 177:3705–10. doi: 10.4049/jimmunol.177.6.3705

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Keşmir C, De Boer RJ. A mathematical model on germinal center kinetics and termination. J Immunol. (1999) 163:2463–9.

PubMed Abstract | Google Scholar

22. Allen CD, Cyster JG. Follicular dendritic cell networks of primary follicles and germinal centers: phenotype and function. Semin Immunol. (2008) 20:14–25. doi: 10.1016/j.smim.2007.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Ebert LM, Horn MP, Lang AB, Moser B. B cells alter the phenotype and function of follicular-homing CXCR5+ T cells. Eur J Immunol. (2004) 34:3562–71. doi: 10.1002/eji.200425478

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Lim HW, Kim CH. Loss of IL-7 receptor alpha on CD4+ T cells defines terminally differentiated B cell-helping effector T cells in a B cell-rich lymphoid tissue. J Immunol. (2007) 179:7448–56. doi: 10.4049/jimmunol.179.11.7448

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Schmitt N, Ueno H. Blood Tfh cells come with colors. Immunity. (2013) 39:629–30. doi: 10.1016/j.immuni.2013.09.011

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Bentebibel SE, Lopez S, Obermoser G, Schmitt N, Mueller C, Harrod C, et al. Induction of ICOS+CXCR3+CXCR5+ TH cells correlates with antibody responses to influenza vaccination. Sci Transl Med. (2013) 5:176. doi: 10.1126/scitranslmed.3005191

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Weinstein JS, Herman EI, Lainez B, Licona-Limón P, Esplugues E, Flavell R, et al. TFH cells progressively differentiate to regulate the germinal center response. Nat Immunol. (2016) 17:1197–205. doi: 10.1038/ni.3554

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Jandl C, King C. Cytokines in the germinal center niche. Antibodies. (2016) 5:5. doi: 10.3390/antib5010005

CrossRef Full Text | Google Scholar

29. Baumjohann D, Preite S, Reboldi A, Ronchi F, Ansel KM, Lanzavecchia A, et al. Persistent antigen and germinal center B cells sustain T follicular helper cell responses and phenotype. Immunity. (2013) 38:596–605. doi: 10.1016/j.immuni.2012.11.020

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Eno L, Beumee JGB, Rabitz H. Sensitivity analysis of experimental data. Appl Math Comput. (1985) 16:153–63. doi: 10.1016/0096-3003(85)90005-0

CrossRef Full Text | Google Scholar

31. Rabitz H, Kramer M, Dacol D. Sensitivity analysis in chemical kinetics. Annu Rev Phys Chem. (1983) 34:419–61. doi: 10.1146/annurev.pc.34.100183.002223

CrossRef Full Text | Google Scholar

32. Vassiliadis VS, Canto EB, Banga JR. Second-order sensitivities of general dynamic systems with application to optimal control problems. Chem Eng Sci. (1999) 54:3851–60. doi: 10.1016/S0009-2509(98)00432-1

CrossRef Full Text | Google Scholar

33. MacLennan IC. Germinal centers. Annu Rev Immunol. (1994) 12:117–39. doi: 10.1146/annurev.iy.12.040194.001001

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Szakal AK, Kosco MH, Tew JG. A novel in vivo follicular dendritic cell-dependent iccosome-mediated mechanism for delivery of antigen to antigen-processing cells. J Immunol. (1988) 140:341–53.

PubMed Abstract | Google Scholar

35. Suzuki K, Grigorova I, Phan T, Kelly L, Cyster J. Visualizing B cell capture of cognate antigen from follicular dendritic cells. J Exp Med. (2009) 206:1485–93. doi: 10.1084/jem.20090209

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Heesters B, Myers R, Carroll M. Follicular dendritic cells: dynamic antigen libraries. Nat Rev Immunol. (2014) 14:495–504. doi: 10.1038/nri3689

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Reshetova P, van Schaik BDC, Klarenbeek PL, Doorenspleet ME, Esveldt REE, Tak PP, et al. Computational model reveals limited correlation between germinal center B-cell subclone abundancy and affinity: implications for repertoire sequencing. Front Immunol. (2017) 8:221. doi: 10.3389/fimmu.2017.00221

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Natkanski E, Lee WY, Mistry B, Casal A, Molloy JE, Tolar P. B cells use mechanical energy to discriminate antigen affinities. Science. (2013) 340:1587–90. doi: 10.1126/science.1237572

PubMed Abstract | CrossRef Full Text | Google Scholar

39. McKay MD, Beckman RJ, Conover WJ. Comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics. (1979) 21:239–45. doi: 10.1080/00401706.1979.10489755

CrossRef Full Text | Google Scholar

40. Blower SM, Dowlatabadi H. Sensitivity and uncertainty analysis of complex models of disease transmission: an HIV model, as an example. Int Stat Rev. (1994) 229–43. doi: 10.2307/1403510

CrossRef Full Text | Google Scholar

41. Bai P, Banks HT, Dediu S, Govan AY, Last M, Lloyd AL, et al. Stochastic and deterministic models for agricultural production networks. Math Biosci Eng. (2007) 4:373–402. doi: 10.3934/mbe.2007.4.373

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Faro J, Or-Guil M. Reassessing germinal centre reaction concepts. In: Molina-París C, Lythe G, editors. Mathematical Models and Immune Cell Biology. New York, NY: Springer (2011). p. 241–58. doi: 10.1007/978-1-4419-7725-0_12

CrossRef Full Text | Google Scholar

43. Wang P, Shih CM, Qi H, Lan YH. A stochastic model of the germinal center integrating local antigen competition, individualistic T-B interactions, and B cell receptor signaling. J Immunol. (2016) 197:1169–82. doi: 10.4049/jimmunol.1600411

PubMed Abstract | CrossRef Full Text | Google Scholar

44. El Shikh ME, El Sayed RM, Sukumar S, Szakal AK, Tew JG. Activation of B cells by antigens on follicular dendritic cells. Trends Immunol. (2010) 31:205–11. doi: 10.1016/j.it.2010.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

45. El Shikh MEM, Pitzalis C. Follicular dendritic cells in health and disease. Front Immunol. (2012) 3:292. doi: 10.3389/fimmu.2012.00292

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Denton AE, Linterman MA. Stromal networking: cellular connections in the germinal centre. Curr Opin Immunol. (2017) 45:103–11. doi: 10.1016/j.coi.2017.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Heesters BA, van der Poel CE, Das A, Carroll MC. Antigen presentation to B cells. Trends Immunol. (2016) 37:844–54. doi: 10.1016/j.it.2016.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Okada T, Miller M, Parker I, Krummel M, Neighbors M, Hartley S, et al. Antigen-engaged B cells undergo chemotaxis toward the T zone and form motile conjugates with helper T cells. PLoS Biol. (2005) 3:e150. doi: 10.1371/journal.pbio.0030150

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Papa I, Saliba D, Ponzoni M, Bustamante S, Canete PF, Gonzalez-Figueroa P, et al. TFH-derived dopamine accelerates productive synapses in germinal centres. Nature. (2017) 547:318–23. doi: 10.1038/nature23013

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Kelsoe G. The germinal center: a crucible for lymphocyte selection. Semin Immunol. (1996) 8:179–84. doi: 10.1006/smim.1996.0022

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Wollenberg I, Agua-Doce A, Hernández A, Almeida C, Oliveira VG, Faro J, et al. Regulation of the germinal center reaction by Foxp3+ follicular regulatory T cells. J Immunol. (2011) 187:4553–60. doi: 10.4049/jimmunol.1101328

PubMed Abstract | CrossRef Full Text | Google Scholar

52. McHeyzer-Williams MG, Davis MM. Antigen-specific development of primary and memory T cells in vivo. Science. (1995) 268:106–11. doi: 10.1126/science.7535476

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Wittenbrink N, Klein A, Weiser AA, Schuchhardt J, Or-Guil M. Is there a typical germinal center? A large-scale immunohistological study on the cellular composition of germinal centers during the hapten-carrier-driven primary immune response in mice. J Immunol. (2011) 187:6185–96. doi: 10.4049/jimmunol.1101440

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Faro J, Or-Guil M. How oligoclonal are germinal centers? A new method for estimating clonal diversity from immunohistological sections. BMC Bioinformatics. (2013) 14(Suppl. 6):S8. doi: 10.1186/1471-2105-14-S6-S8

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Tas JMJ, Mesin L, Pasqual G, Targ S, Jacobsen JT, Mano YM, et al. Visualizing antibody affinity maturation in germinal centers. Science. (2016) 351:1048–54. doi: 10.1126/science.aad3439

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Gitlin AD, Shulman Z, Nussenzweig MC. Clonal selection in the germinal centre by regulated proliferation and hypermutation. Nature. (2014) 509:637–40. doi: 10.1038/nature13300

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Humphrey JH. The fate of antigen and its relationship to the immune response. The complexity of antigens. Antibiot Chemother. (1969) 15:7–23. doi: 10.1159/000386767

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Tew JG, Mandel TE. Prolonged antigen half-life in the lymphoid follicles of specifically immunized mice. Immunology. (1979) 37:69–76.

PubMed Abstract | Google Scholar

59. Vora KA, Ravetch JV, Manser T. Amplified follicular immune complex deposition in mice lacking the Fc receptor gamma-chain does not alter maturation of the B cell response. J Immunol. (1997) 159:2116–24.

Google Scholar

60. Allen CDC, Okada T, Cyster JG. Germinal-center organization and cellular dynamics. Immunity. (2007) 27:190–202. doi: 10.1016/j.immuni.2007.07.009

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Schmiel SE, Yang JA, Jenkins MK, Mueller DL. Cutting edge: adenosine A2a receptor signals inhibit germinal center T follicular helper cell differentiation during the primary response to vaccination. J Immunol. (2017) 198:623–8. doi: 10.4049/jimmunol.1601686

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Bantug GR, Galluzzi L, Kroemer G, Hess C. The spectrum of T cell metabolism in health and disease. Nat Rev Immunol. (2017) 18:19–34. doi: 10.1038/nri.2017.99

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Buck MD, Sowell RT, Kaech SM, Pearce EL. Metabolic instruction of immunity. Cell. (2017) 169:570–86. doi: 10.1016/j.cell.2017.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Chung Y, Tanaka S, Chu F, Nurieva RI, Martinez GJ, Rawal S, et al. Follicular regulatory T cells expressing Foxp3 and Bcl-6 suppress germinal center reactions. Nat Med. (2011) 17:983–8. doi: 10.1038/nm.2426

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Linterman MA, Pierson W, Lee SK, Kallies A, Kawamoto S, Rayner TF, et al. Foxp3+ follicular regulatory T cells control the germinal center response. Nat Med. (2011) 17:975–82. doi: 10.1038/nm.2425

PubMed Abstract | CrossRef Full Text | Google Scholar