Synergy and antagonism in the integration of BCR and CD40 signals that control B-cell population expansion

Abstract

In response to infection or vaccination, lymph nodes must select antigen-reactive B-cells while eliminating auto-reactive B-cells. B-cells are instructed via B-cell receptor (BCR), which binds antigen, and CD40 receptor by antigen-recognizing T-cells. How BCR and CD40 signaling are integrated quantitatively to jointly determine B-cell fate decisions remains unclear. Here, we developed a differential-equations-based model of BCR and CD40 signaling networks activating NFκB. The model recapitulates NFκB dynamics upon BCR and CD40 stimulation, and when linked to established cell decision models of cell cycle and survival control, the resulting cell population dynamics. However, upon costimulation, NFκB dynamics were correctly predicted but the predicted potentiated population expansion was not observed experimentally. We found that this discrepancy was due to BCR-induced caspase activity that may trigger apoptosis in founder cells, unless timely NFκB-induced survival gene expression protects them. Iterative model predictions and sequential co-stimulation experiments revealed how complex non-monotonic integration of BCR and CD40 signals controls positive and negative selection of B-cells. Our work suggests a temporal proof-reading mechanism for regulating the stringency of B-cell selection during antibody responses.

Keywords: B-cell Selection; B-cell Signaling; Mathematical Modeling; NFκB; Systems Immunology.

Figure 1. B-cell signaling model recapitulates early NFκB dynamics in response to T-dependent stimulation.

(A) Schematic of BCR-CD40 receptor model to recapitulate T-dependent activation of B-cells. Each of the 37 biochemical reactions capture a process like phosphorylation, dephosphorylation, synthesis, degradation, association, dissociation, or catalysis, and is annotated with a number consistent with the corresponding reaction and parameter value in Table 1. Species with a colored background are in their active form, while those with white background are in their inactive form. Circle with ‘P’ inside indicates phosphorylation. (B) Schematics of existing T-independent (left) and newly integrated T-dependent (right) NFκB signaling modeling frameworks. T-independent stimulation typically only involves a single ligand (e.g., CpG or LPS), while T-dependent stimulation always involve a more complex receptor signaling system of both BCR and CD40 ligands. Low, medium and high dose of CD40 were set to 6 nM, 12 nM, and 30 nM, respectively, to correspond to the three experimental doses we used (1 μg/mL, 3.3 μg/mL, and 10 μg/mL). (C, D) Line graph from model simulations show (C) peak nuclear RelA (in the form of RelA:p50 and RelA:P52 heterodimers) and (D) peak nuclear cRel (in the form of cRel:p50 and cRel:P52 heterodimers) levels in response to increasing α-CD40 doses, where the shading represents the sample standard deviation of 1000 cells. X-axes are plotted on a log-scale to accommodate a wide range of concentrations. (E) Line graphs from model simulations of 1000 virtual cells and (F) matching experiments with 600 K founder B-cells show temporal trajectories of nuclear RelA level at 0, 7, 24, and 48 h following stimulation with low α-CD40 (1 μg/mL), high α-CD40 (10 μg/mL), or costimulation with high α-CD40 and α-BCR (10 μg/mL each). Darker colored lines represent the average nuclear RelA level from 1000 cells and the lighter shading represents the sample standard deviation of the 1000 cells. (G, H) Line graphs of nuclear cRel level in matching stimulation conditions as (E, F). Source data are available online for this figure.

Figure 2. Multi-scale model recapitulates B-cell population dynamics in response to CD40 stimulation.

(A) Workflow of fitting model simulations to experimental B-cell population dynamics following stimulation. Left: schematic of full T-dependent modeling framework. Right: experimental workflow with Cell Trace Far Red (CTFR) dye dilution. (B) Stacked area plots from model simulations of 1000 virtual B-cells show their population dynamics in response to stimulation with (from left to right) no (0 nM), low (6 nM), medium (12 nM), and high (30 nM) dose of α-CD40. Each subsequent generation of proliferating cells is indicated with a darker gray. (C) Stacked area plots from matching experiments of 19196 founder B-cells show their population dynamics in response to no (0 μg/mL), low (1 μg/mL), medium (3.3 μg/mL), or high (10 μg/mL) dose of α-CD40. (D, E) Root mean square deviation (RMSD) is calculated between simulated and experimental data, and is composed of 2 scores: RMSD of (D) relative population size expansion and RMSD of (E) generation composition. An example of RMSD between model and experimental data is shown on the left side of (D) and (E), and a heatmap of the RMSD scores in matching (diagonal) or mismatching (off-diagonal) model-and-experiment pairs is shown on the right side. Lower RMSD scores correspond to better fit. Source data are available online for this figure.

Figure 3. Model predicts synergistic population expansion in response to BCR and CD40 costimulation, but experiment reveals dose-dependent interaction between the stimuli.

(A) Stacked area plot from model simulations of 1000 virtual B-cells show their population dynamics in response to costimulation with high α-BCR (0.25 nM) and low α-CD40 (6 nM). Each subsequent generation of proliferating cells is indicated with a darker gray. (B) Stacked area plot from model simulations of 1000 virtual B-cells show their population dynamics in response to costimulation with high α-BCR (0.25 nM) and high α-CD40 (30 nM). (C) Stacked area plot from matching experiments as (A) of 19196 founder B-cells show their population dynamics in response to high α-BCR (10 μg/mL) and low α-CD40 (1 μg/mL) costimulation. (D) Stacked area plot from matching experiments as (B) of 19196 founder B-cells show their population dynamics in response to high α-BCR (10 μg/mL) and high α-CD40 (10 μg/mL) costimulation. (E) Stacked bar graph shows a breakdown of total RMSD by types in the 2 costimulation conditions compared to the 4 model-and-experiment pairs in Fig. 2B,C which includes no, low, medium, and high dose of α-CD40. (F) Line graph of experimental population expansion index is higher in response to costimulation than without α-BCR. (G) Stacked bar graph of experimental generation composition dynamics in response to low α-CD40 stimulation with or without high α-BCR costimulation. (H) Stacked bar graph of RMSD score between the two experimental conditions in (G) shows the addition of high α-BCR changes both population expansion and generation composition. (I) Line graph of experimental population expansion index is lower in response to costimulation than without α-BCR. (J) Stacked bar graph of experimental generation composition dynamics in response to high α-CD40 stimulation with high α-BCR costimulation. (K) Stacked bar graph of RMSD score between the two experimental conditions in (J) shows the addition of high α-BCR predominantly affects population expansion. (L) Bar graph to compare the percentage of actively dying cells (DRAQ7+ cells, N = 100–200 for each bar) in response to 4 stimulation conditions. (M) Line graph of cell area (N = 100–200) in the 4 experimental conditions over time, where error bar represents the standard error of the mean (SEM). Source data are available online for this figure.

Figure 4. BCR-induced caspase-dependent apoptosis and NFκB signaling saturation explains the dose-dependent interaction in costimulation.

(A) Schematic of updated T-dependent multi-scale B-cell model where activated BCR induces caspase-8 processing and subsequently triggers cell death (reaction highlighted in red). (B) Stacked area plots from model simulation (left) and matching experiment (right) show B-cell population dynamics in response to costimulation with high α-BCR (0.25 nM and 10 μg/mL) and high α-CD40 (30 nM and 10 μg/mL) with the addition of BCR-induced caspase processing. Each subsequent generation of proliferating cells is indicated with a darker gray. (C) Stacked area plots from model simulation (left) and matching experiment (right) show B-cell population dynamics in response to costimulation with high α-BCR (0.25 nM and 10 μg/mL) and low α-CD40 (6 nM and 1 μg/mL) with the addition of BCR-induced caspase processing. (D) Bar graph of total RMSDs of the 2 costimulation conditions after the addition of BCR-induced caspase processing compared with before the addition. (E) Line graphs with sample standard deviation (SD) from model simulations show the model-simulated distribution of nuclear RelA (left) and cRel (right) dynamics in response to stimulation with low α-CD40 (light blue) or costimulation with low α-CD40 and high α-BCR (pink). (F) Line graphs with sample standard deviation (SD) from model simulations show the model-simulated distribution of nuclear RelA (left) and cRel (right) dynamics in response to stimulation with high α-CD40 (dark green) or costimulation with low α-CD40 and high α-BCR (purple). (G) Violin plots with corresponding line graphs (N = 100–200) connecting the means from multi-channel fluorescence microscopy (left axis) and triangles (N = 600 K) from Western blot fold-change quantification (right axis) show the experimental distribution of nuclear RelA and cRel dynamics in response to matching condition as (E). Western blot fold-change values were min-max normalized to be on the same scale as microscopy fluorescence. (H) Violin plots with corresponding line graphs (N = 100–200) connecting the means from multi-channel fluorescence microscopy (left axis) and triangles (N = 600 K) from Western blot fold-change quantification (right axis) show the experimental distribution of nuclear RelA and cRel dynamics in response to matching condition as (F). Source data are available online for this figure.

Figure 5. BCR-induced apoptosis prevents BCR stimulation from promoting population growth.

(A) Stacked area plots from model simulations of 1000 virtual B-cells (top) and matching experiments of 19196 founder B-cells (bottom) show their population dynamics in response to costimulation with high (30 nM and 10 μg/mL) α-CD40 and no (0 nM and 0 μg/mL), low (0.005 nM and 1 μg/mL), or high (0.25 nM and 10 μg/mL) dose of α-BCR under the impact of AICD. Each subsequent generation of proliferating cells is indicated with a darker gray. (B) Stacked bar graph from model simulations (top) and experiments (bottom) show a breakdown of live B-cells by generation numbers at 96 h post-stimulation-onset. (C, D) Model-simulated Kaplan–Meier survival curve in response to (C) α-CD40 dose and (D) α-BCR dose shows distinct pattern regarding monotonicity. (E) Bar graph from model simulations show percentage live B-cells at 24 h in response to increasing α-CD40 and α-BCR doses. (F–K) Heatmaps from model simulations of 1000 virtual B-cells in response to 25 single- or co-stimulation scenarios (with 5 doses of α-CD40: 0, 6, 12, 18, and 30 nM, and 5 doses of α-BCR: 0, 0.0005, 0.005, 0.05, and 0.25 nM, combinatorially) show the percentage of survived B-cells at 24 h under (F) no AICD and (I) with AICD, the percentage of proliferative B-cells by 84 h out of those that survived (G) without AICD and (J) with AICD, and the relative population size at 96 h (normalized to founder cell population size) (H) without AICD and (K) with AICD, where white contour lines represent 0.5-, 1.0-, 1.5-, and 2.0-fold changes. (L) Heatmap shows the differences between (F) and (I). (M) Heatmap shows the differences between (G) and (J). (N) Heatmap shows the differences between (H) and (K). In (F–N), 25 simulated doses are plotted as colored circles in a scatterplot, whereas the space in between doses is interpolated with a locally estimated scatterplot smoothing (LOESS) curve. Source data are available online for this figure.

Figure 6. Sequential BCR-CD40 simulation reveals a limited window of opportunity to acquire CD40 signal.

(A) In vitro experimental workflow where primary B cells are sequentially stimulated with pulsing α-BCR, followed by α-CD40 stimulation 1, 3, 5, or 8 h later. (B) Line graph from model simulations of 1000 virtual B-cells show their population expansion in response to sequential costimulation with high α-BCR (0.25 nM) and low CD40 (6 nM) (left), high α-BCR (0.25 nM) and high α-CD40 (30 nM) (middle), and low α-BCR (0.005 nM) and high α-CD40 (30 nM) (right), colored by the gap between α-BCR and α-CD40 stimulation. Each thick colored line represents the average population expansion from 1000 cells, and the shading represents the population standard deviation from the 8 simulations, each with 125 founder cells. (C) Line graph from matching experiments of 19196 founder B-cells show their population expansion at 24, 36, 48, 72, and 96 h in response to sequential costimulation with high α-BCR (10 μg/mL) and low α-CD40 (1 μg/mL) (left), high α-BCR (10 μg/mL) and high α-CD40 (10 μg/mL) (middle), and low α-BCR (1 μg/mL) and high α-CD40 (10 μg/mL) (right), colored by the gap between BCR and CD40 stimulation. The lines are smoothed using Excel’s “smoothed line” function which used a Catmull-Rom spline. Source data are available online for this figure.

Figure 7. BCR-induced anti-apoptotic BclXL protects cells from dying from AICD, as a form of paradoxical signaling.

(A) Schematic of paradoxical BCR signaling that promotes proliferation, survival, and death of B-cells through cMyc, Bcl-xL, and caspase-8, respectively. (B) Line plots of Bcl-xL activity (left axis) colored by caspase-8 level (color bar) in 50 model-simulated single-cells show the correlation between Bcl-xL consumption and caspase-8 activity. The thick line overlaid on top is a Kaplan–Meier survival curve (right axis). The pink and green vertical dashed lines represent the timing of BCR and CD40 stimulation, respectively. From left to right, the 3 conditions are: high α-CD40 costimulation, and sequential α-BCR and α-CD40 stimulation with a 1 h and 8 h gap. When a cell dies, the line continues (and becomes pink due to high caspase-8 level). (C) Line plots of Bcl-xL activity colored by RelA (left) and cRel (right) activity in 50 model-simulated single-cells demonstrate the correlation between NFκB activation and BclXL upregulation in B-cells costimulated sequentially with an 8 h gap. When a cell dies, the line discontinues. (D) Violin plot of peak RelA, cRel, and Bcl-xL activity in 2000 model-simulated B-cells in response to sequential costimulated with an 8 h gap show the differences between cells that died within the first 12 h and those that survived. Statistical significance is evaluated using Mann–Whitney U tests, with p-values of 0.0019, 0.0077, and <1e−18, correspondingly. Mann–Whitney U test is appropriate because the groups in comparison are not normally distributed but are independent and similar in shapes. (E–J) Heatmaps from model simulations of 1000 virtual B-cells in response to 25 single or costimulation scenarios (with 5 doses of α-CD40: 0, 6, 12, 18, and 30 nM, and 5 doses of α-BCR: 0, 0.0005, 0.005, 0.05, and 0.25 nM, combinatorially) show the percentage of survived cell at 24 h after stimulation onset (left) and the percentage of survived cells that proliferated by 84 h (right) in (E) 1 h sequential costimulation or (F) 8 h sequential costimulation. (G) Heatmap highlights the difference between (E) and (F). (H, I) Heatmap of relative population size at 96 h between (H) 1 h sequential stimulation and (I) 8 h sequential stimulation shows the biggest difference (J) in the upper right and lower right corners, where white contour lines represent 0.5-, 1.0-, and 1.5-fold changes. In (E–J), 25 simulated doses are plotted as colored circles in a scatterplot, whereas the space in between doses is interpolated with a locally estimated scatterplot smoothing (LOESS) curve.

Figure EV1. Raw experimental data to test the multi-scale B-cell model.

(A) Immunoblot from experiments with 600 K founder B-cells show nuclear RelA, cRel, and p52 levels at 0, 7, 24, 48, and 72 h after stimulation with low α-CD40 (1 μg/mL), high α-CD40 (10 μg/mL), or costimulation with high α-CD40 and α-BCR (10 μg/mL). (B) Line graph of live B-cell count (fold-change) for each timepoint in (A), to which the samples are adjusted when loading to the gel. The cell count fluctuation is due to cell death, cell division, and technical error when transferring cells. (C) Cell Trace Far Red (CTFR) dye dilution fluorescence histogram for B-cells stimulated with (from left to right) low (1 μg/mL), medium (3.3 μg/mL), and high (10 μg/mL) dose of α-CD40 and costimulation of high α-CD40 and α-BCR (10 μg/mL). There is a baseline shift in CTFR fluorescence by about 2-fold from 24 h to 120 h (dotted line), which we adjusted when deconvolving the cells into each generation. (D) Deconvolution of the time courses in (C) into each generation, where the red line indicates the center of the undivided population of cells, the blue lines indicate individual proliferation peaks, and the green line represents the model sum.

Figure EV2. Multi-scale model needs tuning to recapitulate B-cell population dynamics in response to CD40 stimulation.

(A) Stacked area plots from model simulations of 1000 virtual B-cells (top) and matching experiments with 19196 founder B-cells (bottom) show their population dynamics in response to stimulation with (from left to right) no (0 nM and 0 μg/mL), low (6 nM and 1 μg/mL), medium (12 nM and 3.3 μg/mL), and high (30 nM and 10 μg/mL) dose of α-CD40. Each subsequent generation of proliferating cells is indicated with a darker gray. (B) Heatmap shows RMSD of relative population size expansion (top) and generation composition (bottom) in matching (diagonal) or mismatching (off-diagonal) model-and-experiment pairs. Some model doses (medium and low) are more deviated from their matching than mismatching experimental doses (high and medium, respectively), indicating a subpar fit. (C) Bar graph from local sensitivity analysis of parameters in the cell cycle module shows their standard deviations in time to first division (Tdiv0) and time to later divisions (Tdiv1+). Local sensitivity analysis is achieved by repetitive simulations that independently scaling each parameter in the cell cycle module by 0.2, 0.33, 0.4, 0.5, 0.66, 1.0, 1.5, 2.0, 2.5, 3.0, or 5.0-fold. 2 out of 55 parameters stand out as the best candidates for tuning Tdiv0 and Tdiv1+: retinoblastoma (Rb) decay rate and cyclin B (CycB) synthesis rate, respectively. These parameters were tuned in order to achieve a later and more dose-responsive Tdiv0, shorter Tdiv1+, and smaller divider percentage. (D) Box plots from model simulations of 300 virtual B-cells show the mean Tdiv0 increases for all doses after parameter tuning. (E) Box plots from model simulations of 300 virtual B-cells show the mean Tdiv1+ decreases for all doses after parameter tuning. (F) Pie charts from model simulations of 300 virtual B-cells show the percentage of dividing cells (colored slices) out of all founder cells decreases for all doses, while maintaining CD40 dose-responsiveness. Gray slices are the non-dividing founder cells that either die or survive without division.

Figure EV3. Batch-normalized, background subtracted multi-channel microscopy captures single-cell NFκB dynamics.

(A–D) Multi-channel fluorescence microscopy images of live triple-reporter B-cells at 63X/NA1.4 under oil immersion. Panels from top to bottom show mTFP1-cRel (blue) and mVenus-RelA (yellow) cellular localization, with H2B-mCherry (pink) as a marker distinguishing nuclear from cytoplasmic compartment, all compared from left to right at 0 h baseline, 7 h, 24 h, 48 h, and 72 h post-stimulation with (A) high α-CD40 (10 μg/mL), (B) high α-BCR (10 μg/mL) and high α-CD40 (10 μg/mL), (C) low α-CD40 (1 μg/mL), and (D) high α-BCR (10 μg/mL) and low α-CD40 (1 μg/mL).

Figure EV4. Model-simulated cytoplasmic BclXL level recapitulates experimental results.

(A) Immunoblot from experiments with 600 K founder B-cells show cytoplasmic Bcl-xL and β-tubulin levels in response to stimulation with (from left to right) high (10 μg/mL) dose of α-CD40, high α-CD40 and high α-BCR, high α-BCR, and sequential stimulation of high α-BCR and high α-BCR with a 4 h delay. (B, C) Bar graphs from model simulations (top) and experiments (bottom) show consistent max-normalized quantification of cytoplasmic Bcl-xL level at (B) 7 h and (C) 24 h.