1 Faculty of Life Sciences, Bar-Ilan University, Ramat-Gan 52900, Israel 2 Pathology and Laboratory Medicine, University of Pennsylvania School of Medicine, 36th and Hamilton Walk, Philadelphia, PA 19104-6082, USA
Correspondence to: R. Mehr; E-mail: mehrra{at}mail.biu.ac.il
Transmitting editor: I. Pecht
![]() |
Abstract |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Keywords: B lymphocyte subsets, bromo-deoxyuridine, mathematical model
![]() |
Introduction |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
In the pre-B stage (mostly B220dullCD43IgM, Hardys fraction D), cell divisions are terminated and the cells proceed to rearrange their light chain genes. In the immature B cell stage (B220dullCD43IgM+, Hardys fraction E), the cells express IgM BCR on their surface and are subject to strong negative selection as ligand binding yields either cell death or rescue via receptor editing (1724). The threshold for negative selection, in terms of receptor avidity, is much lower than that of mature B cell activation (25). Immature B cells are exported from the bone marrow to the periphery as transitional cells (2629), which are subject to further selection (2830), [reviewed in (31)]. Those transitional cells that survive complete their maturation and join the mature naive peripheral B cell pool (B220brightCD43IgM+IgD+, Hardys fraction F). It has been suggested that the extension of selection to the periphery allows those autoreactive B cells, which express receptors for self-antigens not encountered in the bone marrow, to be deleted from the repertoire before they fully mature and become potentially dangerous. The present paper, however, deals only with the population dynamics of developing B cells in the bone marrow. For the sake of our discussion, the term immature B cells refers to cells that are already expressing their IgM antigen receptors on their surface and are undergoing selection or getting ready to emigrate to the periphery.
One factor that delayed the formation of a complete quantitative model of B cell development was the inconsistency in designation of the different developmental stages, which has been recently resolved (1). Additionally, the results of different attempts to quantify the populations of developing B cells in mouse bone marrow vary greatly [see, e.g. (3238)], due to differences between mouse strains, the use of different markers to identify developing B cell subpopulations, variability in experimental conditions or differences in the methods used to calculate the final numbers. For example, many studies use the percentages of cells in different subsets in bone marrow (obtained by flow cytometry) from one or two femura, multiplied by the total cell count (which is usually much less precise than determination of subset composition), then multiplied by the estimated ratio of the number of cells in total bone marrow to the number in one femur. For the latter ratio there is a single estimate of 15.8 (36), which is used in all calculations, even though this ratio probably varies between individual mice and between different mouse strains.
Estimating kinetic parameters has proven equally difficult. The rates of differentiation from one population to the next can only be estimated based on an a priori model of population dynamics, using data from labeling studies and making assumptions about how much of the labeling in each compartment reflects proliferation versus entry of cells from a previous compartment (35). The various estimates of population growth rates are based on different methods of enumerating cycling cells, which do not always give consistent results (33,39). Additionally, population growth rates are usually estimated as the slope of the linear part of the labeling curve (40), which is wrong (41), because the labeling rate is, at all time points, a non-linear function of all rate parameters for each population and no part of it is truly linear. The estimation of cell death rates is the most difficult, because dead cells are rapidly cleared from the bone marrow (4244). Partially blocking programmed cell death can give an estimate of the numbers of cells that would be otherwise lost (45), but it is difficult to deduce the overall population death rates from these numbers.
An assumption inherent in the models used so far (reviewed above) to obtain kinetic parameter estimates is that developing B cell population dynamics are virtually synchronous and unidirectional. Cells are depicted as progressing in a straightforward manner from one developmental stage to the next, as if on a conveyor belt. However, there are several reasons to doubt that B cells develop along a synchronous, irreversible program. (i) Gene rearrangement success or failure and selection for functional heavy chains introduce variability in the time cells take to traverse the early pro-B stage and in the number of cell divisions performed by the cells which have succeeded in rearranging and expressing a functional heavy chain (46). (ii) The number of cell divisions may vary depending on the level of expression of the heavy chain and the components of the surrogate light chain, the affinity between BCR components, the strength of the signal transduced and the availability of microenvironmental support (47). (iii) Variability is also introduced during the rearrangement of light chain genes in the late pre-B stage. (iv) Marked variation in residence times in the pre-B and immature B cell subsets is introduced by the existence of secondary rearrangements. An example of the inter-cell variability is the recent discovery that a single event of receptor editing causes a delay of at least 2 h in developing B cell progression and that the extent of editing is quite high, at least 25% of developing cells (48). Repeated receptor editing events may thus lead to considerable delays. One may ask, then, how many attempts of rearrangement is a cell allowed to make, and how long does it take to go through the process of rearrangement attempt, failure, and initiation of the next rearrangement and testing cycle. The first question is dealt with in other studies [(49) and our work in preparation]; the second is more complex and will be discussed below.
The findings concerning secondary rearrangements raise the following question: is the progression from pre-B to immature B cells a unilateral, irreversible process? All published calculations have based the estimates of differentiation and death rates in these compartments on the assumption that B cell progression from the pre-B to the immature B cell stage is unidirectional and irreversible. However, the identification of cells as pre-B or immature B cells is largely based on surface IgM expression, which may be reduced or even completely abrogated as a cell undergoes receptor editing. Some down-regulation of BCR expression in immature B cells following signaling through their BCR has indeed been observed (50).
The above findings point at the possibility that some of the cells identified as pre-B cells may be re-entrants from the immature B cell compartment, i.e. cells in the process of performing secondary rearrangements. This possibility has not, so far, been included in the analyses used to estimate kinetic parameters of developing B cell populations. It has also been noted that the number of immature B cells exported to the periphery per day accounts for only about a quarter of the calculated production from the immature B compartment (35). Such a discrepancy may be attributed to cell loss; however it can also be partially due to an inaccurate estimate of cell production in the pre-B cell compartment, which ignores the possibility of cells refluxing from the immature into the pre-B compartment.
Below, we present the first quantitative model for the population dynamics of developing B cells. Using this model we attempt to reconcile experimental estimates of the composition and kinetic parameters of B cell subpopulations and examine the assumption that B cell development is a synchronous and unidirectional process. We show that inter-cell variability and phenotypic refluximmature B cells returning to the pre-B stage for secondary rearrangementsprovide a better explanation of experimental observations.
![]() |
Methods |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Lymphocyte suspensions
Briefly, bone marrow cells were obtained from the two hind limbs of donor animals, and prepared by flushing the femurs and tibias with cold medium and aspirating the cell suspension using sequentially smaller bore needles. Erythrocytes are depleted by incubation in either Geys solution or ammonium chlorideTris.
Antibodies to cell-surface antigen and immunofluorescent analyses
The following reagents were purchased from PharMingen (San Diego, CA): phycoerythrin (PE)- and FITC-labeled anti-CD24 (heat stable antigen) (M1/69); PE-labeled anti-CD43 (leukosialin) (S7); allophycocyanin- and PE-labeled anti-CD45R (B220) (RA3-6B2). Biotin-labeled goat anti-mouse IgM, PE-labeled anti-IgD (SBA-1) and streptavidinalkaline phosphatase were purchased from Southern Biotechnology Associates (Birmingham, AL). FITC-labeled anti-bromo-deoxyuridine (BrdU) (B44) was purchased from Becton Dickinson (San Jose, CA).
Continuous BrdU labeling and cytofluorimetric analysis
Mice were injected i.p. with 0.6 mg BrdU (Sigma, St Louis, MO) in 0.2 ml PBS at 12-h intervals for the duration of each experiment. Cells from BrdU-treated mice were then stained for appropriate surface markers (PE, allophycocyanin and biotin-conjugated antibodies followed by Red 670streptavidin) and washed with cold PBS. BrdU incorporation was analyzed according to previously published procedures. Briefly, the cells are permeabilized by dropwise addition of ice-cold 95% ethanol, washed, and then fixed in PBS with 1% paraformaldehyde and 0.05% Tween 20. The cells were then incubated in buffer plus 100 U/ml DNase to partially degrade and denature their chromatin. Finally, the cells were stained with FITC-labeled mAb to BrdU. Cytometric analyses were performed on a FACSCalibur flow cytometer that was calibrated using the manufacturers beads (Becton Dickinson). Data were collected by gating on all nucleated cells. Data were analyzed using CellQuest (Becton Dickinson) and FlowJo (Tree Star, San Carlos, CA) software.
Statistical analysis
The fitting criterion for the case of two interdependent populations, based on the multivariate normality assumption, is minimization of the expression:
where Ykt refers to the set of experimental measurements, fkt refers to the set of simulation results and these are compared for the two subpopulations, indexed by k (with k' = 1 if k = 2 and vice versa), at each time point t for which there is an experimental result (t = 0,...,5 in units of 12 h from the beginning of labeling). The number of data points at time t is given by nt and is the correlation between the data points for the two populations. The likelihood ratio test statistic for this case, used in obtaining the P value for comparison between the two hypotheses, is:
where is the variance of Ykt.
![]() |
Results |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
The results (Fig. 1) suggest that some cells transit through the various stages faster than others. The labeling kinetics in the pro-B compartment are non-linear and slower, indicating a variability in the timing and length of cell cycle, which would be best visualized as a pool of cells with a wide distribution of times in each compartment. The same is true for the pre-B compartment. The labeling kinetics in the immature B cell compartment are more or less linear, implying that cells entering this compartment do not divide and exit the compartment after some characteristic length of time. On the other hand, the labeling kinetics in the pro- and pre-B compartments are clearly non-linear, exhibiting gradual saturation to a maximal value.
A mathematical model of B cell development
In order to extract more information from our experimental data, we formulated a mathematical model of the population dynamics of B cell development, used here to examine the concept of reflux from the immature back to the pre-B cell compartment. The model is based on differential equations representing developing B lymphocyte subsets. This addresses the issue of intercellular variability, by representing each compartment as a pool of cells, where we fix only the probability (per unit time) that a cell would leave this pool by using a population rate of transition. Cell division and death are similarly treated as probabilistic events by using population rates. A schematic representation of the model is presented in Fig. 2.
|
The input of stem cells into the pro-B compartment is denoted by s (for source) in units of cells per time unit. For the sake of convenience we chose the time unit for which parameter values were defined to be 6 h, which is a minimum estimate for the cell cycle time of bone marrow cells (51). All other parameters are rate parameters, defined in units of fractions of cells per time unit. The rates of differentiation between subsets are represented by X, with X representing the compartment the cells are differentiating from. Proliferation (population growth) is assumed to occur only in the late pro-B compartment (with the rate denoted by
o) and the early pre-B compartment (with rate
e).
Proliferation of developing B cells is known to be limited by the finite space and resources (e.g. contact with the stroma, growth factors, nutrients) in the bone marrow (10,32). Hence the term for proliferation in each compartment is multiplied by a logistic growth-limiting factor: (1 Bo/Ko) for pro-B cells and (1 Be/Ke) for pre-B cells, with Bo = Bor + Boc and Be = Bec + Ber. Ko and Ke denote the carrying capacities of the pro- and pre-B compartments, respectively, i.e. the population sizes for which the corresponding population growth rates become zero. While our model is the first mathematical model of the population dynamics of developing B lymphocytes, much work has been done on modeling the population dynamics of developing T lymphocytes, which are similar in many aspects to developing B lymphocytes (5257). For developing T lymphocytes, the logistic model has been found to be the most appropriate and we had no reason to assume that the growth limitations are different for B lymphocytes. Never theless, we have checked other models as well, defining the growth-limiting term in a more general way as (1Bo/Ko)p and found p = 1 to be the best model (data not shown).
Cell death is assumed in our model to occur only in the non-proliferating cell subsets, because proliferation, gene rearrangement and selection occur in distinct stages and cell death usually occurs only as a result of failure in the latter two processes. The corresponding population mortality rates are denoted by µo, µe, µi for Bor, Ber, Bi respectively.
In immature cells undergoing receptor editing, the surface expression levels of IgM may be very low (50), due to internalization of the old receptor, until a new productive rearrangement is achieved and a new receptor is expressed. Hence, for a short time, the editing cells may be indistinguishable from pre-B cells on the basis of sIgM expression alone. Note that this reflux is only phenotypic; the refluxed cells differ from true or primary pre-B cells in the sense that they have already expressed high levels of sIgM in their past. We implement the reflux assumption by allowing a small fraction, r, of immature B cells to return to the pre-B compartment in each simulation time unit (Fig. 5). The term
rBi is subtracted from the equation for Bi and added to the equation for Ber. Thus, our null hypothesis is that there is no reflux (
r = 0) and it is compared to the hypothesis of reflux (
r > 0).
|
The expressions for the steady state, combined with some more biological knowledge, lead to several constraints on our parameter values, which will help us choose parameter ranges for the simulations. First, as mentioned above, for all populations to be positive, all parameters must be non-negative; differentiation rates (the s) must be strictly positive. Second, for the steady state to be stable, the following additional condition (derived in the Appendix) must be met by the parameters:
(o
oc)(
e
ec) > 0(3)
Third, pro-B cells need a continuous contact with the bone marrow stroma in order to survive and develop further, while pre-B cells do not need that much contact with stroma and can thrive as long as the soluble growth factor IL-7 is supplied by stromal cells (3,58). Hence, Ke > Ko, which allows the pre-B compartment to grow to larger cell numbers.
Modeling labeling kinetics of developing B cells
To model the BrdU labeling experiments, we divided all B cell subsets (other than Bor) to two subsets, unlabeled (Boc,Bec,Ber,Bi) and labeled (LBoc,LBec,LBer,LBi). We assume cells in the first subset, Bor, do not divide and hence are not labeled. Fractions AC are non-dividing and hence cannot be labeled (7). As uncommitted stem cells may have been labeled to a minor extent, however, we have tested this point by adding to our simulations the assumption that up to 20% of the input cells may be labeled. This did not result in any significant change to the results presented here (data not shown), because in the early pro-B subset there is much cell death and most of the expansion occurs in later subsets. Hence we neglect this effect henceforth.
Cells in the Boc subset move to the labeled LBoc subset upon dividing and a similar move occurs in the Bec subset (Fig.3). We assume that once a cell is labeled, it remains labeled throughout the experiment. The new equations, describing labeled and unlabeled populations, are given in the Appendix. These equations were integrated using a simple C program. In the simulations, we first run the model without labeling for 100 time units (which is long enough for cell numbers to arrive at their steady-states), then run them with labeling for 10 units (60 h, as in the experiments) and then run them for the remaining time without labeling. The latter part of each run simulates the decay of labeling over time, under the assumption that this decay is only due to cell death or migration, neglecting any possible reduction of the amount of BrdU due to cell divisions. The decay of labeling was not measured in the experiments, hence our simulations predict what the dynamics of decay would be. This suggests an additional way to check our model, beyond fitting to existing data. The total numbers of cells in both compartments and the fractions of cells in each compartment exhibit the same behavior as in the basic model, because we assume that labeling does not change the kinetics of B cell development.
Model parameters
In choosing parameter values for the simulations of our model, we adhered to the following guidelines. (i) The parameters should be in the experimentally observed orders of magnitude, if published information is available. A summary of the published information is given in Table 1. While these estimates (where available) suffer from the problems discussed above and, additionally, are usually not given in units of population rates, so that interpretation of most of this data depends on the model used, they were useful in suggesting the appropriate ranges for some of the parameters. (ii) The steady-state values obtained using these parameters should be in agreement with our experimental observations on both the total numbers and the composition of bone marrow B cells. Any parameter set which did not conform to this criterion was rejected. (iii) The time of arrival to the steady state should be reasonable (about 1 simulated month). (iv) The fraction of labeled cells in each of the three subpopulations should be within the error bars of the experimental observation for at least five out of the six time points. We call the number of time points for which the simulated fraction of labeled cells for a given population the score for that simulation. The reason we do not require that the score will be 6 for all subpopulations is that there is a large variability at t = 0, where theoretically there should have been no labeled cells. All simulations which had a good fit to the experimental data fell within this criterion (see below).
|
Since the parameters governing the behavior of pre-B and immature B cells have no effect on the behavior of pro-B cells, we first conducted an automated search for the values of the six parameters governing the behavior of pro-B cells. These are s, µo, or,
o, Ko and
oc. According to the criteria above, we rejected parameter sets which resulted in total pro-B cell numbers or fraction of labeled pro-B cells outside the experimental error (i.e. we used seven data points to determine six parameter values). Among the acceptable parameter sets, we looked for the best fit to the data on pro-B cells, defined as the parameter set which gave the minimum value of the sum of squared deviations from experimental data points. The search was performed by repeating the simulations for 1020 possible values for each parameter, i.e.
106 simulations, and recording those parameter sets which gave acceptable results according to the criteria above.
Since the parameters governing the behavior of immature B cells may (if r>0) affect the behavior of pre-B cells, we could not perform independent parameter fittings for these two populations; instead, we performed a automated parameter search based on both populations, i.e. a search for the values of
e, Ke,
ec, µe and
er, plus
r, µi and
i. Thus, here we used 14 data points (total pre-B and immature B cell numbers and fraction of labeled pre-B and immature B cells at six time points) to determine the values of eight parameters. The fitting criterion for this case of two interdependent populations and the likelihood ratio test statistic for this case, are given in Methods.
Simulation results
The parameter set which gave the best fit to experimental data is given in Table 2. A simulation with this set of default parameters, run for 240 time units (60 days), arrives at the steady-state value of 4.53 x 107 cells. In terms of percentages, the pro-B cells consist of 9.0% of total cell numbers, the pre-B 60.3% and immature B cells 30.7%. These numbers are within the experimentally observed ranges. The kinetics of all B cell populations in this simulation are presented in Fig. 4. Simulated labeling kinetics obtained with our best-fit set of parameters, including the post-labeling decay of BrdU-labeled cell fractions, are given in Fig.5.
|
|
The ranges of parameter values that give results within the experimental range (i.e. total cell numbers in each population are within the experimental range and fractions of labeled cells are within the experimental range for at least five out of six time points in each population) are also given in Table 2. Note that the ranges are given here for each parameter separately, hence not all values in the range given for one parameter necessarily work with all values in the range given for other parameters. For example, within the ranges given for µe and er, the pairs that give acceptable results vary depending on the value of
r used (Fig. 6A). The same is true for
e and
ec (Fig. 6B). In this case it is clear that only pairs that obey
e >
ec are acceptable, because these are the only sets that obey the stability constraint (
o
oc)(
e
ec) > 0 (equation 11) with
o >
oc. The values of µi and
i are also dependent on
r; pairs are acceptable as long as their sum remains sufficiently small, so that immature B cell numbers do not become too low (Fig. 6C).
|
|
![]() |
Discussion |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Using mice transgenic for autoreactive BCR, receptor editing has been identified as one of the mechanisms of central tolerance (1921). Computational models of lymphocyte antigen receptor gene rearrangement have shown that cells choose the gene segments to be rearranged in a semi-ordered manner (49), thus maximizing the usage of available gene segments and allowing several rearrangement attempts per cell. If a rearrangement is merely non-productive, the cell will either die or perform a secondary rearrangement. On the other hand, if a rearrangement is productive, the cell will express the resulting antigen receptor and await the judgement of negative selection. While expressing its antigen receptor, the cell will be counted among the immature B cells. However, if, due to negative selection signals, the cell is forced to edit its receptor, it is conceivable that during the time period between negative selection and the expression of the edited receptor, the cell will be IgM low to negative (50) and hence be counted again among the pre-B cells. Thus there is a small probability of phenotypic reflux after each rearrangement attempt. Further support for this possibility comes from studies which show that in the pre-B stage, cells can be found with multiple productive light chain gene rearrangements (59).
The probability that a cell will perform such a reflux after a given rearrangement attempt can be calculated as follows. The probability that a rearrangement is productive is 1/3. The probability that the resulting light chain will pair with the cells heavy chain, Ppair, has been estimated as 0.8 [estimates of all probabilities mentioned here are reviewed in (49)]. Out of all expressed receptors, it has been estimated that about 2/3 will be negatively selected. Thus, a maximum of (Ppair/3)*(2/3) = 2Ppair/9 of all cells will revert to the pre-B phenotype and edit their receptors after each rearrangement. With Ppair = 0.8 (49), this corresponds to about 0.18 of the cells at most refluxing to the pre-B phenotype after each rearrangement. Obviously the actual number should be lower, because not all cells edit their receptors upon negative selection; some of them die in the immature stage. The recent finding that a single event of receptor editing causes a delay of at least 2 h in developing B cell progression (48) means that, at most, 0.18 of the cells in a given cohort reflux to the pre-B phenotype after 2 h, which corresponds to a maximum of 0.54 of cells refluxing after 6 h. In our simulations, the fraction of refluxing cells was found to be
0.4 per 6 h, which is lower than the upper bound above, but is of the same order of magnitude. This unexpected agreement lends further support to our conclusions and also points at the possibility that the kinetics we observe result from multiple such delays.
Two new, testable predictions, which stand in contrast to the intuition on which all the published work on estimates of population dynamic rates (cited above) is based, are thus made by our model. The first is that some immature B cells undergoing receptor editing will down-regulate their surface IgM expression and hence will be phenotypically indistinguishable from pre-B cells. The second prediction is that the loss of B cells in the immature stage is lower than previously estimated. Historically, the difference between the rate of appearance of transitional cells in the periphery and the rate of production of immature B cells in the bone marrow has been assumed to be due to cell death at the immature stage. Our results imply that much of this difference may be due to reflux to the pre-B cell stage, where cells are again as likely to fail in the secondary rearrangement as in the primary rearrangement, so that most cell loss is due to failed rearrangements rather than direct BCR-mediated activation of apoptosis.
Phenotypic reflux is only one possible explanation for the data presented here. Other explanations might also be suggested, e.g. there may be several pathways of pre-B to immature B cell differentiation, which proceed in parallel, but at different rates. The evaluation of such explanations would require the examination of yet unresolved phenotypic/functional subsets and further mathematical modeling. The combination of these methods will enable us in the future to examine alternative explanations of the observed labeling kinetics. The main point to bear in mind, however, is that estimates based on the unidirectional, synchronous model are likely to be inaccurate.
![]() |
Acknowledgements |
---|
![]() |
Abbreviation |
---|
![]() |
Appendix: Model equations and analysis |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Equations (4)(8) have only one all-positive (and hence biologically meaningful) steady state, which is stable. This steady state is given by:
where
and:
These steady-state expressions, combined with biological knowledge, lead to several constraints on our parameter values, which helped us choose parameter ranges for the simulations. First, as mentioned in the text, for all populations to be positive, all parameters must be non-negative; differentiation rates (the s) must be strictly positive. Second, for the steady state to be stable, the following additional condition must be met by the parameters:
Since all the parameters are non-negative and all differentiation rates (the s) are positive, this condition is equivalent to (
o
oc)(
e
ec) > 0 (equation 3). Third, as discussed in the model section, Ke>Ko, as the pre-B compartment is less constrained than the pro-B compartment and can grow to larger cell numbers.
When we model labeled and unlabeled cell subsets separately, the equations become:
where Bo = Bor+Boc+LBoc and Be = Ber+Bec+LBec+LBer denote the total numbers of pro-B and pre-B cells, respectively. The total steady-state cell numbers remain unchanged.
|
|
![]() |
References |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|