The GC model was constructed and simulated as a hybrid, agent-based stochastic model in CompuCell3D. CompuCell3D provides a flexible and customizable platform for simulating multi-cellular behaviors and interactions based on the Glazier-Graner-Hogeweg approach (57). The Cellular Potts Model module in CompuCell3D was employed to simulate the physical properties and movements of individual B cells (58), while the molecular network operating in each individual B cell, as depicted in Figure 1 , was simulated by using the Gillespie’s stochastic algorithm implemented in Tellurium conforming to the Antimony notation (59, 60). The CompuCell3D model consists of four files: an XML file, a Potts initialization file (PIF), and two Python script files. The XML file contains various “Plugins” and “Steppables” that define some default Potts model parameter values. The Chemotaxis plugin defines CXCL12 and CXCL13 as the chemoattractants, and the DiffusionSolverFE steppable designates that CXCL12 and CXCL13 are secreted by CRC and FDC, respectively and specifies the parameter values for secretion, diffusion, and decay in the Medium. The PIF file contains the initial coordinates of medium and cells where applicable. The steppable Python file contains the script that defines several steppable classes, including GCR_Steppable, MitosisSteppable, BCell_GRNSteppable, and VisualizationSteppable, and the Tellurium model. The model is initialized in the “start” section of GCR_Steppable, including the generation of CRCs, FDCs, Tfh cells, and seeding B cells. Each B cell is assigned a Tellurium molecular network model named as BcellNetwork.

Variables of each B cell are saved in a plain text file which is updated every 15 Monte Carlo steps (mcs). The file is named in the format of “Generation_mother ID_cell ID.txt” to facilitate lineage tree construction. The CompuCell3D model files containing the parameter values and Python script for analysis of simulated data are available at https://github.com/pulsatility/2024-GCR-Model.git.

The morphological results of a representative simulation of the GC model are shown in Figure 2 . The simulated space dimension is 250x200x11 pixels, which can be considered as 250x200x11 in µm in real space, to represent a slice of the GC to save computational time. The 2-D projection of the instantiated non-B cells on the X-Y plane is indicated in Figure 2A , while along the Z dimension, these cells are distributed randomly in 3 of the 11 layers (results not shown). There are 45 CRCs distributed on the left half of the field, which becomes the future DZ, and 45 FDCs and 36 Tfh cells distributed on the right half of the field, which becomes the future LZ. Tfh cells are located next to FDCs, reflecting the notion that they also express CXCR5 and are thus drawn to FDCs which secrete chemoattractant CXCL13 (61, 62). Not all FDCs are surrounded by Tfh cells, mimicking the situation that the availability of Tfh cells is a limiting factor in the positive selection of LZ B cells (12, 63–65). CRCs and FDCs secrete CXCL12 and CXCL13 respectively, establishing two opposing chemoattractant fields and thus the polarity of the GC ( Figures 2B, C ).

With the layout of residential cells and chemoattractant fields as above, a simulated GC that succeeds in producing B cells of high BCR/antibody affinities is shown as snapshots in Figure 2D and in Supplementary Video S1 . For this simulation, the GC starts with 200 B cells (clones) of intermediate affinity of 5 (indicated by the color of the cells) between the DZ and LZ. Over a period of 40,000 Monte Carlo steps (mcs, where 100 mcs can be regarded approximately as 1 hour in real time), both the number of total GC B cells and the fraction of high-affinity B cells increase, indicating successful GC population growth and affinity maturation. The 2-D trajectories of 3 select B cell lineage branches leading to different cell fates are shown in Figure 2E . These trajectories cycle between the DZ and LZ for multiple rounds. The first trajectory ends with cell death in the DZ due to damaging mutation during mitosis (left panel), the second trajectory also ends with cell death but in the LZ due to death timer (mid panel), and the last trajectory ends with differentiation into a PC in the LZ (right panel).

The above patterns of phenotypical behaviors of B cells are underpinned by the molecular network operating in each B cell responding to extracellular signaling cues in the GC, including chemoattractants CXCL12 and 13, BCR-mediated antigen signaling by engaging FDCs, and CD40 signaling by engaging Tfh cells ( Figure 1 ). Examining the steady-state GC population at 700,000 mcs revealed that the expression or activity levels of the signaling molecules and genes in the 2352 B cells in the DZ and LZ are highly heterogeneous ( Figure 7A ). While there is no pMHCII expression in the DZ, a large fraction of B cells in the LZ expresses high pMHCII levels ( Figure 7B ), and those expressing low pMHCII levels are B cells that either have just entered the LZ or are en route returning to the DZ with pMHCII being downregulated. A small fraction of B cells in the LZ exhibits high AKT levels as they interact with FDCs ( Figure 7C ), which causes downregulation of FOXO1 ( Figure 7D ). Signaling downstream of CD40 signaling are cRel and MYC, which are active in a small fraction of B cells in the LZ ( Figures 7E-G ). MYC activity induces AP4, which remains upregulated for an extended period of time, even after the B cells have returned to the DZ and MYC has been downregulated ( Figure 7H ), to sustain proliferative bursts. The majority of B cells in the DZ express CXCR4 and the small fraction with CXCR4 downregulated are B cells that have exited the cell cycle and are about to migrate to the LZ ( Figure 7I ). There is no difference in CXCR5 expression in B cells between the two zones ( Figure 7A ). A few cells in the LZ exhibit high RelA ( Figure 7A ) and BLIMP1 ( Figure 7J ) expression representing emerging PCs. The correlations between select pairs of the signaling molecules are shown in Supplementary Figure S2 . There is a strong negative correlation between FOXO1 and AKT ( Supplementary Figure S2A ), positive correlations between cRel and CD40 ( Supplementary Figure S2C ), MYC and CD40 ( Supplementary Figure S2D ), MYC and cRel ( Supplementary Figure S2E ), and AP4 and MYC ( Supplementary Figure S2F ) in B cells in the LZ, while CXCR4 is only expressed in high FOXO1-expressing B cells ( Supplementary Figure S2B ).

We next examined the relationship between the signaling molecules and phenotypical behaviors of B cells that are positively selected in the LZ and return to the DZ. Both the peak level ( Figure 8A ) and duration of MYC expression ( Figure 8B ), which is quantified by the area under the curve (AUC), are positively correlated with the BCR affinity of the cells. The correlation indicates that the strength of the affinity-dependent BCR signaling, which downregulates FOXO1 and upregulates pMHCII to enable Tfh-dependent CD40 signaling, is quantitatively transmitted to MYC. Although MYC is only transiently expressed in the positively selected B cells in the LZ, its encoding of affinity is relayed to AP4. By integrating the MYC signal, AP4, which has a longer half-life than MYC, can be activated for a much longer time as B cells migrate back to the DZ. Its peak level is positively correlated with the AUC of MYC in DZ-returning B cells ( Figure 8C ). The 95^th percentile burst size is positively correlated with the mean AP4 peak level in each affinity bin ( Figure 8D ).

Lastly, the trajectory of a representative single B cell lineage branch that successfully makes it to a PC is presented ( Figure 9 ). As the seeding B cell and its progeny cycle between the DZ and LZ ( Figure 9A ), multiple cell divisions (between 3–6 divisions) occur in each proliferative burst ( Figure 9B ) with increasing cell generation ( Figure 9C ). In this particular simulation result, nearly every cell division results in an increase in BCR affinity despite a few occasions when the affinity decreases ( Figure 9D ). Each time the B cell starts the trip to the LZ the countdown of the death timer is initiated, but in this case it never drops below the predefined death threshold of 50 before the cell is rescued by positive selection where the death timer is reset ( Figure 9E ). Every time the B cell moves into the LZ, pMHCII is re-expressed proportionally to the antigen-specific affinity after encounter with FDCs that triggers BCR signaling ( Figure 9F ). BCR signaling transiently activates AKT ( Figure 9G ) which downregulates FOXO1 transiently ( Figure 9H ). Commitment to cell cycle after the B cell is positively selected allows upregulation of CXCR4 by FOXO1, which drives the B cell to return to the DZ, and CXCR4 is downregulated after the B cell exits the last cell cycle of a proliferative burst ( Figure 9I ), which allows the cell to migrate to the LZ due to constitutively expressed CXCR5 (not shown). In the LZ, encounter with Tfh cells triggers transient activation of the CD40-cRel-MYC axis in the presence of downregulation of FOXO1 ( Figures 9J-L ). By integrating the MYC signal, AP4 is upregulated for an extended period of time, which lasts well into the DZ ( Figure 9M ) to sustain the proliferative bursts. The proliferative bursts terminate when AP4 drops below a predefined threshold of 50. When the affinity increases past a predefined threshold of 10, the BCR signaling triggers RelA activation in a switch-like manner ( Figure 9N ), which in turn activates BLIMP1 that drives the B cell to terminally differentiate into a PC ( Figure 9O ). A representative result of a B cell lineage that ends in death in the DZ due to damaging BCR mutations is shown in Supplementary Figures S3A-S3C , where at the time the damaging mutation occurs the affinity drops to “-1” as an indication, and caspase 3 is upregulated to trigger apoptosis. A representative result of a B cell lineage that ends in death in the LZ because of not being positively selected is shown in Supplementary Figures S3D-S3F , where the death timer dips below the threshold of 50 thus triggering cell death.

Emerging GCs are seeded with up to a few hundred B cell clones (68, 75). Starting with 200 B cells, our GC simulation showed that the DZ: LZ ratio of the numbers of B cells oscillates at early time points ( Figure 3A ). This occurs because there are not many B cells at this stage of the GC and these cells tend to migrate in sync between the two zones. As the GC B cell population approaches a steady state, the DZ: LZ ratio converges to a constant value, nearly 2.6:1, in our simulation, which is comparable to the 2.15~2.2:1 ratio observed in both mouse and human GCs (67). The total number of B cells in the simulation can reach about 2500. Given that the modeled GC space contains only 11 pixels (equivalent to 11 µm) along the Z-dimension to keep the simulation time tractable, we expect that when scaling up by 2–5 times to mimic the actual GC thickness of 20–50 µm (76), the peak number of B cells will reach thousands to over 10,000, consistent with the estimated B cell counts in real GCs (68).

Our simulation showed that new B cells are born in the DZ, while cell deaths occur in both the DZ and LZ at an overall rate that eventually matches the birth rate as the GC approaches the steady state. Although the absolute number of deaths in the DZ is more than twice that in the LZ ( Figure 3I ), the percentage death rates are comparable, at about 40–43% per 600 mcs in both zones ( Figure 3J ). These numbers are concordant with the nearly 50% death rate per 6 hours reported for the GCs in Peyer’s patches in mice and GCs in mice immunized with 4-hydroxy-3-nitrophenylacetyl (NP)-conjugated ovalbumin (NP-OVA) or HIV-1 envelope antigen GT1.1 (15). The steady-state death turnover rate of 50% per 6 hours is expected given that GC B cells proliferate with a cell-cycle length of 4–6 hours (10, 24) and in our model the cell volume doubling time of proliferating B cells is parameterized at 500 mcs. The apoptotic deaths in the DZ are believed to occur in the late G1 phase, triggered by AID-induced BCR-damaging mutations during transcription, including stop codons, insertions and deletions in the Ig sequences, such that the synthesized BCR proteins fail to properly fold and be expressed on the cell surface (14, 15).

In the LZ, the default fate of B cells is apoptosis if not positively selected, regardless of their BCR affinity, and the apoptosis will occur when the death timer goes off, which is initiated after the B cells exit the cell cycle in the DZ (15, 77). Our simulations showed that only a small fraction of deaths in the LZ result from not being positively selected after contact with FDCs and Tfh cells, while the majority of the deaths are due to no access to Tfh cells at all ( Figure 3M ). This result recapitulates the current notion that the availability of Tfh cells is the limiting factor, not the competition for antigen, for positive selection (12, 63–65). Interestingly, our simulations also predicted that at the advanced GC stage, a small fraction of cell deaths in the DZ may also result from the lack of access to Tfh when these cells do not have enough time to migrate through a densely populated DZ to reach the LZ ( Figures 3N, O ) and some of these dead cells could be of high-affinity. Whether deaths of such nature occur in vivo and to what extent remain to be tested experimentally using techniques including intravital imaging of apoptotic GC B cells, single-cell cloning of Ig genes from dying B cells, and recombinant antibody expression as in (15).

The dynamics of the GC B cell population is determined primarily by the cell birth rate and death rate. The loss of B cells as a result of GC exit as memory and PCs is expected to be negligible because they only account for <3% of the GC B cell fates (69–71). Our simulation result is consistent with this estimate: PCs emerge at a fraction of only as high as 2% of the B cell population towards the end of the simulated GC ( Figure 4C ). For a GC B cell population to grow or to avoid population collapse, the average proliferation rate has to be higher than the death rate. Absent any limiting factors, the cell population dynamics in the GC may operate as a positive feedback system, producing bistability of two alternative outcomes – expansion or regression – as predicted by previous mathematical models (78). In the expansion mode, as the overall affinities increase, the probability of cell death in the LZ owing to lack of positive selection decreases, thus more B cells will return to the DZ and proliferate with a larger burst size there, which results in a higher birth rate of progeny cells with potentially even higher affinities returning to the LZ, and the cycle repeats leading to GC B cell expansion. In the regression mode, the positive feedback works in the opposite direction, where lower affinities can lead to fewer B cells positively selected in the LZ and smaller burst size in the DZ, which eventually leads to GC regression. This all-or-none type of GC outcome is consistent with the population bottleneck proposition (79) and suggests that there could be an initial affinity threshold condition for those activated B cells that seed a GC, below which a tangible GC is unlikely to emerge and above which a GC will likely emerge and grow. This suggests that it may take B cells of some intermediate affinity to initiate a GC to produce optimal humoral immune outcomes, i.e., high production of high-affinity PCs. A GC starting with B cells of too low affinity will likely abort prematurely as argued above, while a GC starting with high-affinity B cells may grow but only to a small size before some B cells hit the affinity threshold that triggers terminal differentiation to PCs. This may also explain the observation that affinity selection for memory B cells is less stringent, and they are often formed and then exit the GC when their affinities are still at low or intermediate levels, while the PCs are in general of high affinity (71, 80–83). Upon secondary infection or booster immunization whether the recalled memory B cells directly differentiate to PCs or have to go through GCR again can be determined by many factors (84, 85), and it is likely that their BCR affinity may play a role in this regard.

When a GC reaches a certain size, its growth could be restricted by multiple factors, before other GC-shutoff mechanisms such as antigen depletion and antibody feedback kick in. In the present study, we showed that the availability of Tfh cells is a limiting factor, where a higher fraction of B cells in the LZ die because of lack of access to Tfh cells, thus increasing the overall death probability which balances out the increasing birth rate due to improved antigen affinity. That B cells become PCs once their affinities reach a threshold also helps the GC B cell population to reach an equilibrium by curbing further increase in the number of higher-affinity B cells and thus attenuating the positive feedback mechanism described above. Another limiting factor not considered in the present study is the limited nutritional and energetic resources that could restrict B cell proliferation once the GC has grown to a mature size (68).

To grow or maintain a GC the overall damaging BCR mutation-induced cell death probability cannot be higher than 50% for each cell division. Actually it has to be much lower than 50% because not all but only a small fraction of B cells arriving in the LZ are positively selected and return to DZ. In our model, we encoded a DZ death probability of 30% for each of the two daughter cells born from a cell division, which leads to a probability of 49% (0.7*0.7) to double the number of B cells after each division, of 42% (2*0.3*0.7) to keep the number of B cells constant, and of 9% (0.3*0.3) to eliminate the proliferating B cell. The DZ re-entry probability was estimated to be between 10–30% through mathematical modeling and analysis of experimental data (7, 12, 35). These values suggest that absent any DZ death, the average proliferative burst size has to be greater than 1.73–3.32 divisions to grow a GC. When taking into consideration DZ death, which is very significant, on par with LZ death (15), the average burst size has to be much higher. The DZ re-entry probability depends on the BCR antigen affinity, thus it is likely that at the early stage of GC the re-entry probability is low and at the advanced stage it is high. In our model, the positive selection and thus DZ re-entry probability is set to be proportional to pMHCII, such that when the affinity is intermediate at 5, the probability is 50% and when the affinity approaches 10 or higher, the DZ re-entry probability is 100%. However, the overall DZ re-entry fraction is only 36% in our simulation ( Figure 5A ), which can be attributed in part to the inaccessibility to Tfh cells. A re-entry fraction of 36% requires at least a burst size of 1.47 to grow the GC in the absence of DZ death. With a DZ death probability of 30% for each new born B cell, our model has an average burst size of 2.2 ( Figure 6E ), which is in general agreement with the estimated average of 2 divisions (35, 66, 73) or 3 divisions per burst (74) in the literature.

While the average proliferative burst size is 2–3, each burst can vary between 1–6 divisions (35, 66, 73, 74). This range is quantitatively captured in the burst size distribution produced by our simulation ( Figure 6E ). The right-tailed burst size distribution could result intrinsically and in part from the probabilistic damaging mutation-induced cell death after each cell division, which increases the chance of short bursts but limits the highest attainable number of divisions in a proliferative burst even for high-affinity B cells. GC B cells positively selected are guaranteed to divide once, while the number of additional divisions or the burst size is directly proportional to the amount of antigen captured by B cells from FDCs and presented to Tfh cells (66, 86). Our model reproduces this positive association ( Figures 6G, H ). Moreover, because of potential premature termination of a proliferative burst induced by damaging mutation, the affinity appears to be better correlated with the top attainable burst size than the average burst size.

The translation of BCR antigen affinity into burst size is mediated molecularly by two key transcription factors: MYC and one of its target genes, AP4. Because of the transient nature of B cell interactions with FDCs and Tfh cells (65) and the short half-life of MYC (77), MYC is only transiently expressed in a small fraction of LZ B cells (49, 87). The MYC expression level is in direct proportion to the amount of antigen captured and dictates the proliferative burst size (86). While MYC can initiate cell growth and cell cycle by driving LZ B cells into the S phase, its lasting effect on cell proliferation is mediated by AP4, which is induced by MYC in a delayed fashion in positively selected GC B cells and is sustained after the B cells re-enter the DZ (50). Our model recapitulates the spatiotemporal dynamics of MYC by showing transient MYC expression in LZ B cells, its positive correlation with BCR antigen affinity, and sustained AP4 expression.

Our model recapitulates a typical affinity maturation process along with GC growth ( Figure 4A ). The progressive increase in the mean affinity of the GC B cell population is not because all or the majority of the seeding B cell clones improve their affinities uniformly. Rather, in most cases, the mature GC B cell population is dominated by progenies of one or a few of the initial 200 seeding clones ( Figure 4F ; Supplementary Figures S1B, S1D ). This simulation result of terminal clonal dominance is consistent with the premise that GCs mature oligo-clonally (88, 89). More recently, using multiphoton microscopy and sequencing, Tas and colleagues further revealed that a GC can start with tens to hundreds of distinct B cell clones but loses the clonal diversity over time, converging to one or a few parallelly expanding clones (75). The single dominant clone can constitute 10–100% of the final GC B cell population. They further showed that clonal dominance can be achieved through neutral competition, due to stochastic effect, even when all seeding B cells have equal affinity and cannot undergo SHM, a finding that can be explored with our model in the future.

Beltman et al. analyzed time-lapsing imaging data of GCs and revealed that B cells move at a net speed of 0.2–0.3 µm/min toward the LZ to produce a DZ-to-LZ migration time of a few hours (72). In our simulation, the DZ-to-LZ migration time is about 438 ± 205 mcs ( Figure 5B ), consistent with that estimated by Beltman. The LZ-to-DZ migration time in our simulation is 581± 205 mcs ( Figure 5F ) which is about 33% longer than the DZ-to-LZ migration time. The longer time is consistent with the experimental observations (12), and could be attributed to the fact that B cells returning to the DZ have to move against the much heavier incoming traffic of B cells migrating from DZ to LZ. The overall LZ residence time is 1105 ± 244 mcs ( Figure 5J ). These travel times can be tuned by varying the positions and distributions of CRCs and FDCs/Tfh cells in the DZ and LZ respectively, the CXCL12 and 13 gradients, and the chemotaxis strength parameters in the model.

Many mathematical GC models have been developed in the past two decades using a variety of computational approaches, including deterministic, stochastic, agent-based, and hybrid ones. Focusing on various aspects of the GCR and simulating at various biological scales, these models have greatly aided our understanding of this long known phenomenon by investigating possible modes of mechanisms of GC dynamics (34, 35, 72, 90–94), exploring optimal design of vaccination schemes (95–99), and predicting GC-associated disease outcomes (100). While these models focused on the complex processes of affinity maturation and B cell population dynamics, rarely were molecular networks included to drive the B cell behaviors and GC evolution. Quantitative multiscale mathematical models of GC dynamics have been proposed as predictive frameworks to translate basic immunological knowledge to practical challenges (36, 101). In recent years, modeling efforts in this direction have emerged. For example, Merino Tejero and colleagues have developed multiscale GC models integrating molecular and cellular responses (37, 39), by combining the agent-based model developed in (35) and ODE-based gene regulatory network model comprising BCL6, IRF4, and BLIMP1 developed in (102). The early model was used to study the role of affinity-based CD40 signaling and asymmetric B cell division in temporal switch from memory B cell to PC differentiation and DZ-to-LZ ratio. Lately, the model was adapted to examine the oncogenic effects of genetic alteration of the above key transcription factors on GC-originated diffuse large B cell lymphoma (38). More recently, the model was used to explore the relationship between clonal abundance and affinity as well as affinity variability in intraclonal B cells, with an attempt to make sense of repertoire sequencing data (103).

In comparison, the multiscale spatial GC modeling framework we developed here in the CompuCell3D platform further integrates across the molecular, cellular, and tissue scales and offers several improvements. The framework allows the molecular network to drive multiple cellular behaviors, including B cell growth, division, chemotaxis, survival/death, and PC differentiation, which in turn collectively drive GC tissue pattern formation; reciprocally, the cell-to-cell interactions between B cells and FDCs and Tfh cells drive the responses of the molecular signaling network in B cells. Novel cross-scale strengths include cell cycle and FOXO1-dependent CXCR4 expression driving DZ-reentry chemotaxis, MYC and AP4-dependent cell growth and division burst, and RelA and BLIMP1-dependent PC differentiation. As more molecular species are added to the network, additional cross-scale integrations will become available. Because the simulated B cells comprise multiple pixels, the model allows recapitulation of B cell morphology during chemotaxis and volume growth during cell cycle as well as better mimicking of cell-cell interactions. For future iterations of the model, the CompuCell3D platform can readily include paracrine signaling by ILs and other cytokines secreted by B cells, Tfh cells, and FDCs. Last but not least, with the modular plugins and systems biology markup language (SBML) support, the CompuCell3D platform allows a more structured construction of the GC model that will facilitate future model sharing and integration.

As an initial effort to establish the multiscale GC modeling framework on the CompuCell3D platform, we simulated the GC to the mature stage where the B cell population approaches a steady state. While persistent GCs can exist for months or years under certain viral infections (104–106), and in the Peyer’s patch for mucosal immunity (107), in most other infection or immunization scenarios, GCs eventually regress in 3–4 weeks with mechanisms still incompletely understood. Some GC terminations may be because the antigens stored in FDCs are depleted, the signaling nature of Tfh cells and FDCs has changed, or antibodies produced by the departed PCs circulate back into the GC and block antigen presentation (92, 108, 109). To demonstrate that our modeling framework is capable of self-termination, we simulated a scenario of gradual FDC antigen depletion after the GC matures. The simulation clearly showed that the GC eventually regresses, yet preserving many characteristics seen had it persisted, including DZ: LZ B cell ratio, affinity maturation, and clonal dominance ( Supplementary Figure S4 ).

In the current study, we presented the simulations where all GC seeder B cells have the same initial affinity. Randomizing the initial affinity between 3–7 uniformly produced a similar result on GC B cell population growth, affinity maturation, and clonal dominance ( Supplementary Figure S5 ). Future investigations may explore the outcomes of more complex initial affinity compositions. In addition, it will be interesting to use the model to explore alternative hypotheses for mechanisms underpinning the clonal replacement phenomenon observed in long-lived GCs, as triggered by influenza virus or SARS-CoV-2 infection, where the founder B cell clones of high SHM loads and high affinity in mature GCs are gradually replaced by newcomer B cells that are not necessarily antigen-specific (110). In the current model, memory B cells are not included as a cell fate option. Since like PCs, memory B cells only constitute a very small fraction of the GC B cells fates (109), its exclusion is not expected to affect the overall model behavior. Nonetheless, future iterations of the GC model will include the self-termination and memory B cell formation as driven by the IRF4-BCL6-BLIMP1 network (102), and if needed, CSR, which is believed to occur primarily during pre-GC formation (56).

A small caveat of the current model, as currently implemented on the CompuCell3D platform, is that when the cell density is too high such that the DZ is packed, there is a small chance that some dividing B cells vanish in the DZ when their volumes are too small (as represented by the small left tail of the DZ B cell volume distribution in Figure 3H ). This issue can be prevented in future iterations by imposing a limiting resource for cell growth and division to control DZ cell density. B cells are highly mobile in both the DZ and LZ, with a mobility pattern observing persistent random walk (65, 72, 111). While we did not analyze the B cell mobility in this regard, the CompuCell3D platform, which is based on the Cellular Potts Model that follows the Boltzmann law (58), is capable of simulating persistent random walk (112). In our case, the temperature parameter and local CXCL12 and 13 distribution patterns in the LZ and DZ can be optimized, along with reducing the directional chemotactic forces, to help accentuate the random-walk effect. As indicated above, it was previously showed that the persistent random walk of B cells could be an emergent behavior of B cells in a crowded GC environment (113), thus it will be interesting to inspect our GC model in this regard.

In the current model, the two daughter cells split the parent cell’s volume by following a lognormal distribution. While this approach implemented some degree of cell division asymmetry, in future iterations asymmetrical cell divisions can be better implemented based on experimental studies which showed asymmetrical segregation of pMHCII among the two daughter cells, which plays a role in B cell fate decision-making (114), and was implemented in recent models (37, 39). In the current model the length of each cell cycle in a proliferative burst is targeted as a constant. It has been shown that the S phase, which constitutes the major portion of the cell cycles of proliferating B cells, can be shortened by regulating replication fork progression, while the relative order of replication origin activation is preserved (74). The degree of S-phase shortening depends on the interaction strength of B cells with Tfh cells which in turn depends on BCR antigen affinity. Therefore, positively selected high-affinity GC B cells, upon returning to the DZ, will proliferate not only with a larger burst size but also with accelerated cell cycles. This could be a mechanism to compensate for the tendency of longer DZ residence time of high-affinity B cells due to more cell divisions such that they can return to the LZ sooner. Future iterations of the model may consider incorporating affinity-dependent cell cycle shortening.

The molecular network model running in each B cell uses the Gillespie’s stochastic simulation algorithm, with some of the molecular switching actions implemented as hybrid, rule-based events. For simplicity and parsimony, several genes known to participate in GCR and B cell terminal differentiation are not included in the current implementation of the GC model, including BCL6, IRF4, BACH2, and PAX5. BCL6 is upregulated in antigen-engaged B cells in the early stage of GC initiation, before these cells migrate back into the intrafollicular space and cluster in the GC (115). There does not appear to be cyclic BCL6 expression between the DZ and LZ. Since our current model starts with B cells seeding the GC, not including the initial interactions at the T/B border, the absence of BCL6 in the molecular network should not affect the simulation results. Likewise, IRF4, BACH2, and PAX5 do not seem to exhibit cyclic expression patterns during GCR despite that they are involved in either the initiation of GC, AID induction, or terminal differentiation of B cells to PCs or memory cells (116–118). Therefore, by not including them the current model is made simpler. Nonetheless, BCL6, IRF4, BACH2, BLIMP1, and PAX5 are known to form coupled positive or double-negative feedback loops underpinning multistability-based binary decision-making in B cells (119, 120). In future iterations, the hybrid stochastic and rule-based approach of molecular network simulation can be updated by implementing relevant intracellular feedback circuits comprising these transcription factors to enable bi- and multistability-based decision-making on B cell fates. To this end, single-cell RNA sequencing data of GC cells can be integrated in to the molecular network model to parameterize the molecular abundance of the gene transcripts (70).

There are a number of other considerations that could be potentially added in future iterations as well. Once seeded, a GC reaches peak volume in 3–4 days (121, 122), resulting from an initial clonal expansion without SHM and competitive selection (113, 123–125). Future models should consider implementing this process to bring the GC kinetics to more closely align with such rates of B cell accumulation. Paracrine signals mediated by ILs and other cytokines secreted by B cells, Tfh cells, and FDCs may be considered to better recapitulate the tissue-level molecular signaling milieu in the GC. Mobility of Tfh cells, and their intracellular signal transduction and gen regulatory network can be included. The shapes and locations of residential CRCs and FDCs and as a result the spatial patterns of CXCL12 and 13 gradients can be fine-tuned according to immunohistochemistry data, which may ultimately affect the morphology and polarity of GCs.

In conclusion, we have developed a multiscale spatial computational modeling framework for GC simulation in CompuCell3D. The current model is capable of recapitulating GC features that are both qualitatively and to some extent quantitatively consistent with the literature. Given the complexity of GCR, simulations in such a modeling framework can help investigate a range of research questions on this hallmark event of high-affinity antibody production in response to viral infection and vaccination. Upon further extension and refinement, this open-source modeling framework may also help investigate autoimmunity and lymphoma when the GC goes awry due to genetic or environmental disruptions. Lastly, the GC modeling framework may also be utilized towards building digital twins of the human immune system for precision medicine (126).