To study the molecular mechanisms behind the development of diffuse large B-cell lymphoma, we used an already existing MSM of PC differentiation during the GC reaction (Merino Tejero et al., 2020). At the cellular level, we used an agent-based model of the GC reaction based on an existing model (Meyer-Hermann et al., 2012). The agent-based model runs for 21 days (504 h) at a time resolution of 0.002 h (7.2 s) (Supplementary Table 1). Each cell is represented in a three-dimensional space representing the GC and undergoes different processes depending on their state. B cells at a CB state divide, undergo somatic hypermutations, and differentiate in the DZ, while at a CC state, they undergo selection in the LZ Figure 1. In the LZ, CCs interact and collect Ag from follicular dendritic cells (FDCs). CCs may further interact and compete for T follicular helper cell signals. Both processes are dependent on the affinity of the BcR for the Ag. Selected CCs recycle back to the DZ as CBs that carry internalized Ag. Recycled CBs divide in an asymmetric manner in 72% of the cases where all Ag and transcription factors segregate toward one of the daughter CBs. In the remaining 28% of divisions, Ag and transcription factors are divided symmetrically (Barnett et al., 2012; Thaunat et al., 2012). This was based on a recent study that showed coupled Ag and transcription factor asymmetric division with a large segregation between daughter B cells recapitulate the GC output transition (temporal switch) from MBCs to PCs and a DZ-to-LZ ratio, which is the ratio of CBs to non-apoptotic CCs present in the GC. After division, CBs may exit the GC as an output cell or remain in the GC for further rounds of selection.

At the molecular level, we used an ODE model of a GRN involved in PC differentiation based on an existing model (Martínez et al., 2012). The GRN consists of three transcription factors, i.e., BCL6 (b), BLIMP1 (p), and IRF4 (r). These transcription factors regulate each other, and the concentration levels depend on their dissociation constant (k), transcription (μ), and degradation rates (λ), as shown in Eqs 1–3. Parameter values are shown in Supplementary Table 2. Bcr0 and cd40 parameters were calibrated using the MoonFit graphical user interface (Robert et al., 2018). “Squares” in the equations represent Hill coefficients representing the cooperative binding of transcription factors. The GRN behaves as a bistable system, with one state being the CB- and CC-like state (BCL6 high and BLIMP1 and IRF4 low) and a second state being the PC-like state (BLIMP1 and IRF4 high and BCL6 low). Upstream signals may instruct the system toward one or the other state. Their effect is integrated through BcR and CD40 parameters, as described in Eqs 4, 5. During the CC-FDC interaction in the LZ, CCs receive BcR-mediated signals. The maximum strength of the signal is represented by the bcr₀ parameter (Eq. 4). The transition from CC toward the PC state is reversible when the BcR signal cessates, bringing the B cell to its previous CC state. During CC-T follicular helper cell interaction in the LZ, CCs receive signals through CD40. The strength of the signal is represented by the cd40 parameter (Eq. 5). The transition from the CC to PC state is irreversibly determined by the CD40 signal strength, which is dependent on the affinity of the BcR for the Ag (Supplementary Table 2). d p d t = μ p + σ p k b 2 k b 2 +   b 2 + σ p r 2 k r 2 +   r 2 − λ p p , (1) d b d t = μ b + σ b k p 2 k p 2 + p 2 k b 2 k b 2 + b 2 k r 2 k r 2 + r 2 − ( λ b + B C R ) b , (2)   d r d t = μ r + σ r   r 2 k r 2 + r 2 + C D 40 − λ r r   , (3) B C R = b c r 0 k b 2 k b 2 + b 2 , (4) C D 40 = c d 40   k b 2 k b 2 + b 2 . (5)

The ODE model was adopted, as described in the simulations section, in order to investigate the effect of different genetic alterations in the GC dynamics.

The definition of output cells, either MBCs or PCs, is based on the Ag status and levels of BLIMP1 (Merino Tejero et al., 2020). Recycled CBs carrying Ag are considered Ag-positive (Ag+). During asymmetric division, Ag inheriting CB is Ag+ while the sibling is Ag-negative (Ag−). During symmetric division, Ag is distributed symmetrically and both siblings are Ag−. This was based on a study that recapitulated the predominant features of Ag segregation in vitro by building a mixed model of B-cell mitosis where 25% of the divisions were fully symmetric, and the remaining were asymmetric (Thaunat et al., 2012). Recycled CBs that finished dividing differentiate to PCs when BLIMP1 levels reach the differentiation threshold (≥8.10⁻⁸ M; Figure 2) regardless of their Ag status. Ag+ output cells with BLIMP1 levels below the threshold are considered to be MBCs. Finally, Ag− CBs with BLIMP1 levels below the threshold stay in the GC and recycle back to the LZ as CCs. This definition leads to a transzone migration rate or the DZ-to-LZ ratio close to experimental observations, as described by Victora et al. (2010), and a temporal switch in the output cell production of the MBC and PC dynamics, as described by Weisel et al. (2016).

We investigated the effect of eight different diffuse large B-cell lymphoma models on both molecular and cellular levels. This was performed to show how the MSM can be used to recapitulate the GC dynamics observed in mouse diffuse large B-cell lymphoma models (Zhang et al., 2015; Calado et al., 2010) in terms of B-cell and PC numbers and transcription factor expression levels on day 7 of the GC reaction (day 10 post-immunization (PI); Table 2). Most diffuse large B-cell lymphomas result in alterations, either an increase or decrease, in GC B-cell numbers and/or PCs produced (see Results section). Each diffuse large B-cell lymphoma model simulates one or more genetic alteration(s) (Figure 3). Models were simulated by modifying different ODEs and/or parameters models, as described in Table 1. Six of the models simulated single genetic alterations based on a previous study by Martínez et al. (2012), that is, BCL6 autoregulatory inactivation (BCL ^auto ), BCL6 constitutive expression (BCL^inc), loss of IRF4-mediated BCL6 silencing (IRF^sil), loss of BLIMP1-mediated BCL6 silencing (BLIMP^sil), BLIMP1 inactivation (BLIMP^loss), and constitutive activation of the NF-kB pathway (IRF^inc). The steady-state levels of BLIMP1, BCL6, and IRF4 after CD40 signal were very close or identical to those derived by Martínez et al. (2012) (Supplementary Table 3). One exception was that the BCL^inc model showed a clear difference in the IRF4 and BLIMP1 levels. This is caused by the fact that in our MSM selected CC differentiation into a PC process was around 11 times faster than in the model presented by Martínez et al. The in silico GC B-cell dynamics were compared to data from experimental models carried out by Zhang et al. (2015) and Calado et al. (2010). In our IRF^inc model, the perturbation applied to IRF4 was only a 2-fold increase in IRF4 transcription rates as opposed to the 10-fold increase applied by Martínez et al. (2012). This was carried out to obtain transcription factor levels (data not shown) that resembled both experimental models (Zhang et al., 2015; Calado et al., 2010). Both experimental models showed a similar increase in the levels of BLIMP1 on day 7 of the GC reaction compared to an experimental reference. Nevertheless, the GC B-cell dynamics were different (Table 2). These differences observed between both experiments were due to the NF-kB pathway was activated. In a study by Zhang et al. (2015), constitutive induction of NF-kB-inducing kinase restrained activation of the alternative NF-kB pathway. In a study by Calado et al. (2010), constitutive induction of the NF-kB-inhibiting kinase restrained activation of the canonical NF-kB pathway. While our IRF^inc model did not explicitly include a distinction between alternative and canonical pathways, it was able to reproduce PC dynamics observed by Calado et al. (2010) (See below). Two additional models were created to simulate experimental findings by combining two of the single models, namely, BCL6 constitutive expression in codominance with constitutive activation of the alternative NF-kB pathway (BCL^incIRF^inc), based on experiments conducted by Zhang et al. (2015), and BLIMP1 inactivation in codominance with constitutive activation of the canonical NF-kB pathway (BLIMP^lossIRF^inc), based on experiments conducted by Calado et al. (2010). We repeated all simulations 10 times with different random seeds, which showed that the amount of variability was low. Thus, this was sufficient to observe significant differences between the reference and most of the diffuse large B-cell lymphoma models when looking at the DZ-to-LZ ratio and output cell production (data not shown).

Diffuse large B-cell lymphomas are characterized by spleen and/or lymph node hyperplasia (Zhang et al., 2015; Calado et al., 2010). Furthermore, the same studies showed that most diffuse large B-cell lymphoma genetic mutations result in altered GC dynamics, either through an increase or decrease, in GC B-cell numbers and/or PCs produced. Thus, we wondered about the effect of modeling eight diffuse large B-cell lymphoma models on the GC B-cell dynamics. We performed an example simulation at the scale of one B-cell, in which ODEs and parameters were defined as described in Eqs 1–5 and Supplementary Table 2. The initial state at BCL6 high and BLIMP1 and IRF4 low switched to BLIMP1 and IRF4 high and BCL6 low, especially as a consequence of CD40 signal induction (Figure 2). At the scale of a full GC simulation, we found that the DZ-to-LZ ratio was similar to the affinity-based CD40 signaling simulation (Scenario 2) of one of our previous studies (Merino Tejero et al., 2020) (Figure 4A; Supplementary Table 4). Nevertheless, it was slightly higher than previous experimental observations of transzone migration rates (Victora et al., 2010). Finally, a temporal transition from MBC to PC production was observed in the GC output (Figure 4B).

Thus, in the reference simulation, the DZ-to-LZ ratio was higher than previously observed transzone migration rates of around 2 (Victora et al., 2010) and similar to the affinity-based CD40 signaling simulation (Scenario 2) (Merino Tejero et al., 2020). A latter large-scale study by Wittenbrink et al. (2011) showed that the GC dynamics is highly variable. They performed a volume-to-volume comparison of both GC zones and showed that dividing B cells are present in both zones constantly over time. Thus, the study did not reflect the dividing vs. non-dividing B-cell DZ-to-LZ ratio, and their measurements are not directly comparable to those by Victora et al. (2010). Further experiments would be required to determine the variability in DZ-to-LZ ratios in GC responses to establish whether or not our model produces deviating values. There was a temporal switch of the GC output from mainly MBC production prior to the peak to mainly PC production after the peak. Finally, there was a 42-fold higher PC than MCB production. Studies have shown the number of PCs and MBCs generated spleen and bone marrow during an immune response (Kishi et al., 2010; Weisel et al., 2016; Sugimoto-Ishige et al., 2020). Nevertheless, these numbers represent percentages of observed PCs and MBCs from the total number of splenic or bone marrow cells, which are not translatable to the number of output cells from a single GC and therefore are not directly comparable with our results.

We now compared the effect of several genetic (models defined in Table 1, see Methods) alterations onto the GC dynamics in order to recapitulate the observed GC phenotype in diffuse large B-cell lymphoma. We first aimed to gain insight into the effect of BLIMP1 suppression (BCL^auto and BCL^inc) or loss of function (BLIMP^loss and BLIMP^lossIRF^inc) on the GC B-cell dynamics. While previous studies showed that the GC phenotype in diffuse large B-cell lymphomas characterized by BCL^auto and BCL^inc alterations had unaltered numbers of B cells and PCs produced, they did lead to spleen and lymph node hyperplasia (Zhang et al., 2015). Another study showed that the GC phenotype BLIMP^loss and BLIMP^lossIRF^inc alterations are characterized by an accumulation of B cells, a block of PC differentiation and spleen hyperplasia (Calado et al., 2010).

We analyzed the transcription factor dynamics and B-cell dynamics in the aforementioned models (Table 1). The transcription factor dynamics simulated in single cells showed that for BCL^auto and BCL^inc models, the system evolved to the BCL6 high and BLIMP1 and IRF4 low steady-state unaffected by the CD40 signal (Figures 5A,B). The steady-state level of BCL6 in the BCL^auto model was around four times higher than in the BCL^inc model, as shown by Martínez et al. (2012), while steady-state levels of IRF4 and BLIMP1 were similar in both models which were not observed by Martínez et al. (2012) (see Methods section). Thus, we found that within the timing of PC differentiation during the GC reaction, both BCL^auto and BCL^inc models had a robust BCL6 overexpression that inhibited IRF4 and BLIMP1 expression and blocked PC differentiation. BLIMP^loss and BLIMP^lossIRF^inc models switched to the BCL6 low and IRF4 high steady-state as a result of CD40 signal. Nevertheless, they both failed to express BLIMP1 (Figures 5C,D).

We observed BCL^auto, BCL^inc, BLIMP^loss, and BLIMP^lossIRF^inc models had a significantly increased DZ-to-LZ ratio throughout the GC reaction, which reached up to a 25-fold increase at day 21 (Supplementary Table 4), compared to the reference model (Figure 6A) resulting on an increased abundance of GC B cells. All output cells were MBCs (Supplementary Table 4), and no PCs were produced (Figure 6B). Statistical analysis of MBC (Figure 7A) and PC (Figure 7B) counts showed a 13,000-fold increase and 35,000-fold decrease compared to the reference. Since BLIMP1 levels were virtually 0 in all cases, and no significant difference in the number of MBCs and PCs was observed between the four models. MBC affinity increased up to its maximum as the GC reaction progressed in all models while it remained low in the reference model (Figure 7C).

To summarize, all models lead to an accumulation of CBs in the DZ. No PCs were generated due to the inability of recycled CBs to express BLIMP1. There was a significant increase in the number and affinity of MBCs. The expansion of GC B cells and significant reduction of PCs in BCL^auto and BCL^inc models were in disagreement with observations showing the BLIMP1 suppression leads to unaltered GC B-cell numbers on days 10 and 21 post-immunization and PC numbers on day 10 post-immunization (Table 2) (Zhang et al., 2015). This could be due to a decreased but present BLIMP1 expression in the experimental setup, which was not the case in our model. Finally, BLIMP^loss and BLIMP^lossIRF^inc models were consistent with observations showing BLIMP1 inactivation lead to GC B-cell expansion (Table 2) and a lack of PC production on days 10 and 21 post-immunization (Calado et al., 2010). Nevertheless, when in codominance with constitutive NF-KB pathway activation (BLIMP^lossIRF^inc), on day 21 post-immunization, experimental data showed unaltered GC B-cell numbers as opposed to the expansion of B cells observed in our model (Table 2).

Next, we assessed the effect of constitutive activation of the NF-kB pathway alone (IRF^inc) and in codominance with the BCL6 overexpression (BCL^incIRF^inc) on the overall GC dynamics (Table 1). Previous studies showed that the GC phenotype in diffuse large B-cell lymphomas that are driven by IRF^inc and BCL^incIRF^inc alterations are characterized by unaltered B-cell numbers (Zhang et al., 2015) or slightly decreased B-cell numbers (Calado et al., 2010), in the case of the IRF^inc alteration, on day 10 post-immunization and spleen and/or lymph node hyperplasia (Calado et al., 2010; Zhang et al., 2015). While alterations that led to the IRF^inc model resulted in an increased PC production on day 10 post-immunization according to Zhang et al. (2015), unaltered PC numbers were observed by Calado et al. (2010). Furthermore, alterations that led to the BCL^incIRF^inc model resulted in unaltered PC production on day 10 post-immunization (Zhang et al., 2015).

The transcription factor dynamics simulated in single cells showed that in both IRF^inc and BCL^incIRF^inc models, there was a switch to the BLIMP1 and IRF4 high steady-state prior to the CD40 signal (Figure 8A). Nevertheless, only the BCL^incIRF^inc model had consistently high BCL6 levels due to the effect of the BCL6 overexpression (Figure 8B).

We observed both diffuse large B-cell lymphoma models had a 2-to-3-fold decrease in the DZ-to-LZ ratio on day 21 of the GC reaction (Supplementary Table 4; Figure 9A) and a slight decrease of GC B cells (data not shown) and increase in output cells produced on day 21 of the GC reaction (Supplementary Table 4) compared to the reference model. All output cells were PCs, and no MBCs were produced (Figures 9B,C). Statistical analysis of PC counts on day 7 of the GC reaction (Figure 10A) showed no significant difference compared to the reference model. Interestingly, on day 21 of the GC reaction (Figure 10B), a 1.14-fold increase in PC count compared to the reference model was observed. Furthermore, no significant difference in the PC count between IRF^inc and BCL^incIRF^inc models was observed.

To summarize, IRF^inc and BCL^incIRF^inc models resulted in a decreased DZ-to-LZ ratio and increased production of PCs on day 21 of the GC reaction. No MBCs were produced due to the fast increase of BLIMP1 independent of CD40 signaling and B-cell affinity. On day 7 of the GC reaction, similar PC counts were observed compared to the reference model. Our IRF^inc model was in agreement with observations showing an unaltered number of PCs on day 10 post-immunization and spleen hyperplasia (Calado et al., 2010). Nevertheless, the model also showed unaltered and reduced GC B-cell numbers on days 7 and 18 of the GC reaction, while the previous study showed the opposite trend. (Table 2) (Calado et al., 2010). Despite reproducing the transcription factor and GC B-cell dynamics, our PC results were in disagreement with observations showing an increased number of PCs on day 10 post-immunization (Table 2) (Zhang et al., 2015). This discrepancy suggested that our model was able to recapitulate the effect of constitutive activation of the canonical but not the alternative NF-kB pathway on PC production. Finally, the BCL^incIRF^inc model was in agreement with observations shown by Zhang et al. (2015) (Table 2). Thus, our results suggested that IRF^inc and BCL^incIRF^inc alterations lead to unaltered PC production on day 7 of the GC reaction but could result in increased PC differentiation at later stages of the GC reaction.

Finally, we assessed the effect of loss of IRF4- and BLIMP1-mediated silencing of BCL6 (IRF^sil and BLIMP^sil) on the overall GC dynamics (Figure 11). For IRF^sil and BLIMP^sil alterations, we did not find observations on the GC phenotype in terms of B-cell numbers and PCs produced. Nevertheless, the transcription factor dynamics simulated in single cells showed in both models the system evolved to a BCL6 low, BLIMP1 and IRF4 high steady-state as a result of the CD40 signal transduction. Similar results were observed by Martínez et al. (2012), indicating that redundancy in the GRN makes the system robust against loss of the BCL6 repression by either BLIMP1 or IRF4 but not by both. Steady-state levels of all transcription factors in the IRF^sil model reached similar concentration levels as in the reference.

We observed the IRF^sil model had a similar DZ-to-LZ ratio, GC B-cell numbers (data not shown), and output cells produced at day 21 of the GC reaction (Supplementary Table 4) compared to the reference model (Figure 12A). In the BLIMP^sil model, there was a slight increase in the DZ-to-LZ ratio and GC B-cell numbers (data not shown) as well as a slight decrease in output cells produced at day 21 of the GC reaction (Supplementary Table 4) compared to the reference model (Figures 12B,C). Statistical analysis of the MBC count at day 21 of the GC reaction (Figure 13A) showed no significant difference between the reference and IRF^sil. There was a significant 1.7-fold increase in the BLIMP^sil model compared to IRF^sil and the reference. This could be due to the faster increase of BLIMP1 in BLIMP^sil, indicating a greater repression effect of BLIMP1 over BCL6 than the effect of IRF4 over BCL6. Statistical analysis of the PC count at day 21 of the GC reaction (Figure 13B) showed no significant difference between the reference, IRF^sil, and BLIMP^sil models.