|Home | About | Journals | Submit | Contact Us | Français|
Lyn, a Src-family protein tyrosine kinase expressed in B lymphocytes, contributes to initiation of BCR signaling and is also responsible for feedback inhibition of BCR signaling. Lyn-deficient mice have a decreased number of follicular B cells and also spontaneously develop a lupus-like autoimmunity. We used flow cytometric analysis, BrdU labeling, and our mathematical models of B cell population dynamics, to analyze how Lyn deficiency impacts B-cell maturation and survival. We found that Lyn-deficient transitional 1 (T1) cells develop normally, but T2 cells develop primarily from the T1 subset in the spleen and fail to also develop directly from bone marrow immature B cells. Lyn-deficient T2 cells either mature to the follicular B-cell type at a close to normal rate, or die in this compartment rather than access the T3 anergic subset. The ~40% of wild type follicular cells that were short-lived exited primarily by joining the T3 anergic subset, whereas the ~15% Lyn−/− follicular cells that were not long-lived had a high death rate and died in this compartment rather than entering the T3 subset. We hypothesize that exaggerated BCR signaling resulting from weak interactions with self-antigens is largely responsible for these alterations in Lyn-deficient B cells.
B cells are derived from bone marrow (BM) precursors that transit through a series of differentiation stages associated with immunoglobulin (Ig) gene rearrangement and the assembly of a functional B-cell receptor [1–10]. The expression of surface immunoglobulin allows progression to the immature B-cell stage of development, and if self-reactivity is lacking or is not too great, these cells are exported to the spleen as "transitional" B cells [11–16]. Three transitional subsets of B220+AA4+ cells can be identified in adult murine spleen based on differential sIgM and CD23 expression: AA4+CD23−sIgMhigh (T1), AA4+CD23+sIgMhigh (T2), AA4+CD23+sIgMlow (T3); mature follicular (Fo) B cells are AA4−CD23+sIgMlow .
Lyn is a member of the Src family of protein tyrosine kinases, and is mainly expressed in hematopoietic cells [17, 18]. In B lymphocytes, Lyn acts as both a positive and a negative regulator of BCR signaling pathways . Lyn is activated following BCR ligation, and in its positive role contributes to the initiation of the BCR-signaling cascades by phosphorylating the ITAMs of the Igα/Igβ BCR subunits, and also by phosphorylating the proximal signaling molecule CD19. Lyn is not essential for the initiation of BCR signaling, as this process still occurs in Lyn-deficient cells, suggesting that this function is shared with other Src family of protein tyrosine kinases expressed in B cells, such as Blk and Fyn . In mice deficient in Lyn, Blk, and Fyn, however, a strong block in B cell development is seen early in B cell development, prior to surface immunoglobulin expression, when pre-BCR signaling is required .
In contrast to this redundancy for the positive role of Lyn in BCR signaling, Lyn plays a critical and non-redundant role in the negative regulation of BCR signaling. Upon BCR ligation, Lyn phosphorylates ITIM-bearing co-receptors, including CD22 and FcγRIIB, resulting in the recruitment and activation of the phosphatases SH2 domain-containing phosphatase 1 (SHP1) and SHIP-1. In turn, these phosphatases dephosphorylate signaling components and small molecule signaling intermediates to switch off the activation pathways .
While these biochemical functions of Lyn are well established, how these reactions contribute to the normal functioning of B cells is less well understood. Lyn-deficient mice produce large quantities of anti-nuclear antibodies, which are the defining characteristic of the human autoimmune disease systemic lupus erythematosus [20, 23–27]. In addition, these mice exhibit reduced numbers of mature follicular B cells, a complete absence of marginal zone B cells, and a greater proportion of immature cells with a higher than normal turnover rate [28, 29]. Deficiency of Lyn has an especially strong effect on mature B cells. In mature follicular Lyn−/− B cells, BCR signaling is strongly enhanced [19, 20, 30]. Lyn−/− B-1 cells also exhibit elevated responses to antigen in vitro and spontaneous production of autoantibodies in vivo .
Studies targeted to well-defined B cell subsets are necessary to further our understanding of Lyn signaling abnormalities in autoimmune disease. Here, we used flow cytometric analysis, BrdU labeling, and our mathematical models of B cell development [32–36] to investigate which B cell maturation processes are altered in Lyn−/− mice. In our study, we compared the B cell population dynamics in Lyn−/− mice to those in WT mice. Our mathematical modeling reconfirms the suggestion [34, 37] that in wild-type mice spleen T2 cells develop directly from immature BM B cells as well as from splenic T1 cells. In contrast, we find that Lyn-deficient T2 cells almost all develop directly from the T1 subset in the spleen.
Additionally, we found that transitional and maturing B cells in Lyn−/− mice undergo accelerated death rates in the T2 and mature Fo subsets, with the majority of the dying cells not passing through the anergic T3 stage. The mature Fo subset is about ten-fold smaller than in WT mice but, surprisingly, also undergoes a lower turnover. We hypothesize that the higher death rate and lower rate of passage into the anergic T3 subset in Lyn-deficient mice result, at least in part, from elevated BCR signaling in those B cells that are weakly self-reactive . Thus, the lack of Lyn-mediated attenuation of this signaling leads to cell death rather than anergy .
We performed continuous in vivo BrdU labeling of Lyn−/− and WT mice for 8 days, and sampled transitional and mature subsets on days 2, 3, 4, 6 and 8 (n=4 mice per time point). No differences between WT and Lyn−/− mice were seen in the T1 BrdU labeling kinetics (Figure 1A, p=0.846). However there were significant differences in the T2 labeling kinetics starting at day 4 (Figure 1B, p= 0.001). The fractions of labeled cells in T2 Lyn−/− mice reached 80% by day 8, compared with 50% in WT mice. The small differences between WT and Lyn−/− mice seen in the T3 labeling kinetics from day 4 onwards (Figure 1C) were also significant (p= 0.014), however may or may not be meaningful, as the final fractions at day 8 only differed slightly (30% in Lyn−/− compared with 35% in WT mice). There were also significant differences from day 2 onwards in the mature follicular (Fo) B cell labeling kinetics between WT and Lyn−/− mice (Figure 1D, p= 0.001), where the fractions of labeled cells in Fo mature Lyn−/− mice reached 7%, compared with 4% in WT mice.
In a previous study , we asked whether the T3 subset, which contains most newly formed cells undergoing anergic death, could also include mature B cells destined for elimination. To interrogate this hypothesis and its implications, we applied mathematical models to previously generated in vivo labeling data . Our analyses revealed that the death rate of T3 B cells is far higher than the death rates of all other splenic B cell subpopulations. Furthermore, the model in which the T3 pool includes both newly formed and mature primary B cells destined for apoptotic death gave the best fit to the data, and showed that cell loss at the T3 stage can account for nearly the entire mature B cell turnover. This finding has been incorporated in the models used here.
We fit the BrdU data from WT and Lyn−/− mice to our mathematical model of B cell maturation in the spleen [34, 35], where the numbers of cells in the T1, T2, T3 and mature follicular subsets are represented by the variables T1, T2, T3 and M respectively (Figure 2). The current data did not include BM subpopulations, hence we decided, based on previous studies [34, 37], to include in the model two BM subsets, BM1 and BM2, to represent the cells in the BM that differentiate to T1 and T2, respectively. BM1 cells migrate from the bone marrow and differentiate into the splenic T1 subset with a constant rate denoted by δBM1 (per six hours, which is our simulation time unit), and BM2 migrate from the bone marrow and differentiate into the T2 subset with a constant rate, δBM2. The BM1 and BM2 subpopulations are renewed from previous subsets at the constant rates s1 and s2, and proliferate with rates γ1 and γ2 (Figure 2). T1 differentiate to T2 or T3 at rates δ12 and δ13, respectively. T2 differentiate to T3 or M at rates δ23 and δ2m, respectively; the differentiation of M to T3 is denoted by δm3. The death rates are denoted by μ1, μ2, μ3 and μm for T1, T2, T3 and M, respectively. Based on previous studies, it is assumed that none of the transitional subpopulations are cycling [11, 34, 35].
We fitted our mathematical model of transitional and mature B cell dynamics to the BrdU labeling data, and compared the parameter ranges (i.e. the 95% confidence interval on parameter values, obtained using Akaike's Information Criterion) for the two mouse types. The best fit parameters values are presented in Table 1, and Figure 3 shows the model’s fit to the data with these parameter value sets. Parameter ranges that show no overlaps between the two mouse types – and hence can be considered different – are given in Table 1 in bold. The simulated total cell numbers at each stage of development in the spleen of WT and Lyn−/− mice were within the experimentally observed ranges (Figure 4), for all simulations included in the parameter estimates. The following observations can be made based on the parameter estimates.
Since the labeling in these experiments was lower than in our previous studies [11, 28, 34, 35], which used data from a different lab, we had to assume that the labeling fraction in the BM compartments was less than 100%, that is, not all daughter cells were labeled in each division.
The BM1 and BM2 subpopulations described above represent all developing B cell populations – the proliferating pro-B and pre-B cells, and also the immature subset, part of which is destined to exit the BM and feed the transitional B cell pools. Our simulation results reveal differences in the renewal rates s1 and s2 of the BM1 and BM2 subpopulations, respectively: s1 is lower, and s2 higher (although the s2 ranges are only marginally different), in Lyn-deficient mice relative to the wild type mice. This may imply that, in spite of the similar bone marrow cell numbers in Lyn-deficient and wild type mice, the balance between specific B cell subsets may be altered by Lyn deficiency. However, since the current experiments did not include analysis of bone marrow B cell subsets, we could not investigate this point in further detail. The BM1 proliferation and differentiation rates were found here to be almost the same in WT and Lyn−/− mice; however there were differences in the BM2 kinetic parameters (Table 1). In 96% of the Lyn−/− simulation runs, T2 cells did not develop directly from the BM2 subset (δBM2=0), but rather developed mostly from the T1 subset in the spleen (higher values of δ12). In contrast, in the WT mice, most T2 subset development was from the BM rather than from the T1 subset in the spleen. The BM subset (BM2) did not proliferate in Lyn−/− mice (γ2 =0). Thus, both BM2 proliferation and differentiation are decreased in Lyn−/− mice.
More differences in splenic B cell developmental stages were seen in the mature Fo and T3 subsets (Table 1) than in the less developed compartments. The T3 subpopulation developed mostly or only from the T1 subset in Lyn−/− mice (zero or very small δ23), while in the WT mice the input into the T3 subpopulation came from both the T1 and T2 subsets as well as the mature Fo B cell population. Additionally, in WT mice most splenic B cell death was seen in the T3 compartment, as in our previous study , confirming the assumption that the T3 subset contains mostly anergic cells destined for apoptosis. However in Lyn−/− mice, death was higher than in the WT in the T2 and mature Fo B cell subsets, with most of the death occurring in the mature Fo B cell compartment. Taken together, these results imply that B cells in Lyn−/− mice undergo accelerated death rates in the T2 and Fo compartments, with the majority of the dying cells not passing through the anergic T3 stage. This may be the reason why death in the T3 compartment is lower in Lyn−/− mice than in WT mice: since most cells may die in either the T2 or Fo compartment, we may speculate that in Lyn−/− mice most cells receive negative selection signals strong enough to cause them to die quickly, and only a few receive weaker signals that render them anergic and thus make it to the T3 compartment.
Compared to data we used in our previous study , here the BrdU labeling fractions of the WT mature subset were lower, and the total cell number was higher. Thus, the model could not be fit to the current data as is, and to get reasonable fits we had to include in the simulation a large subset of static, non-dividing, non-dying cells that have the same phenotype and hence are included within the mature subset. These could simply be mature B cells that have not turned over during the 8 days of the experiment. The fraction of these cells (%staticm in Table 1) was higher in Lyn−/− mice than in WT mice.
In a number of studies it has been observed that mice with genetic mutations causing enhanced BCR signaling have substantially reduced numbers of follicular B cells. Examples of this phenomenon include mutations that ablate, Lyn  or SHP1 , and a mutation that enhances CD45 activity . The first two of these molecules act in a feedback inhibitory pathway in which Lyn phosphorylates a tyrosine residue in the cytoplasmic domain of CD22, leading to recruitment of the protein tyrosine phosphatase SHP1, which inhibits BCR signaling at a relatively upstream step. A major function of CD45 is to maintain Src-family tyrosine kinases in a primed state by dephosphorylating a C-terminal negative regulatory tyrosine phosphorylation, and in this way it is thought to regulate the magnitude of signaling by the BCR. The reason for increased BCR signaling caused by these mutations leading to a substantial decrease in the size of the follicular B cell compartment is not known. In the present study, we used in vivo BrdU labeling and flow cytometry to measure the numbers and turnover rates of various immature and mature B cell subpopulations in wild type and Lyn−/− mice. We then applied our most recent mathematical model of B cells development and maturation to these data sets to gain insight into how Lyn deficiency affects B cell maturation and survival. The results obtained from the model indicate that the B cells in Lyn-deficient mice exhibit perturbations at several maturation and survival steps, with an especially dramatic increase in the death rate of follicular B cells. These alterations likely underlie the large decrease in the number of follicular B cells in these mice.
The mathematical model of B cell maturation and survival dynamics that we had previously derived from analysis of wild type mice was found to apply to the data sets developed with Lyn−/− mice, with a few minor alterations. First, we found that it was necessary to assume that the BrdU labeling protocol used in these experiments did not achieve 100% labeling of proliferating B cell precursors, which was consistent with the experimental observations. Second, it was clear that the follicular B cells in both wild type and Lyn−/− mice were comprised of cells that were not uniform in their half-lives, but rather consisted of a minority of cells that were turning over relatively rapidly and a majority of cells that had a turnover that was slow enough it could not be estimated from these experiments. The long-lived follicular B cells included about 85% of follicular B cells in Lyn−/− mice and about 60% of the much greater number of follicular B cells in wild type mice. For the purposes of modeling these data, the very slowly turning over B cells were considered to be "static" cells, or non-dividing, non-dying cells. Please note that we only refer to these cells as "static” because they were not labeled in the 8 days of the experiment; estimating how long they actually live would require longer experiments. Interestingly, we found that such static cells are necessary to explain incomplete cell labeling in some of our other studies  and also in our current work in progress (unpublished data). Moreover, after carefully examining our results in the study with the Allman lab (Figure 2 in  and Figure 3 in ), we see that even there, none of the populations has reached 100% labeling by day 7. Thus, this behaviour seems to be the rule rather than the exception.
When the above two characteristics were incorporated into our mathematical equations of B cell maturation and turnover, the range of maturation rate and death rate parameters needed to model the new wild type data set were similar to the corresponding parameters calculated from earlier data sets . Moreover, it was also possible to obtain compatible sets of maturation rate and death rate parameters with the data from the Lyn−/− mice. Some of those parameters were similar to the wild type parameters, but several of them were strikingly different. Thus, we were able to model successfully the population dynamics of B cell populations from Lyn−/− mice and assess how deficiency of Lyn affects B cell maturation and survival in the spleen.
The earliest B cell maturation state in the spleen, the T1 transitional B cells, labeled with BrdU with very similar kinetics in wild type and Lyn−/− mice, and our mathematical model found that the rate of maturation from bone marrow immature B cells to splenic T1 cells of these two mouse strains was identical (parameter δBM1 in Table 1). These results are consistent with data from mixed bone marrow chimeras, which found that there was not a statistically significant competitive disadvantage of Lyn-deficient B cell precursors accessing the T1 population (Gross et al. manuscript submitted). Those mixed bone marrow chimera experiments did detect a competitive disadvantage of Lyn-deficient B cells in accessing the next maturation state, the T2 stage in the spleen, and moreover the steady state size of this population was decreased by approximately 3-fold in Lyn−/− mice. Our mathematical modeling previously concluded that wild type T2 cells arise both from splenic T1 cells and from bone marrow immature B cells that have reached the T2 phenotype in that location before exiting to the spleen . In contrast, we found here that Lyn-deficient T2 cells only arose from splenic T1 cells with minimal contribution from bone marrow T2-like cells. This difference may in part explain the decreased size of the T2 population in Lyn−/− mice. This finding may be tested by following the kinetics of reconstitution of bone marrow chimeras, where B cells – whether wild type or Lyn−/− – are allowed to develop in irradiated hosts. Following the appearance of the immature populations in a relatively homogeneous wave, as done by Allman et al. , will show how Lyn−/− and wild type B cells compare in terms of reconstitution kinetics, vs. the steady state analysis on which our model is based. The model's predicted kinetics of growth of the modeled B cell populations, simulated with the best fit parameter values for wild type and Lyn−/− cells, are presented in Figure 5. These kinetics clearly show a slower growth of the T2 compartment in Lyn−/− mice compared with its growth in the wild type mice, due to the direct differentiation of immature BM B cells into the T2 compartment, which is almost absent in the Lyn−/− mice.
Once Lyn-deficient cells reached the T2 stage, they were able to mature to the follicular population at a similar but somewhat higher rate in Lyn−/− mice than in wild type mice (δ2m parameter in Table 1). Whereas the rate of maturation of T2 cells was altered by Lyn deficiency to a modest degree, there was a much greater impact on alternative fates of this cell type: in Lyn−/− mice, the T2 cells that didn’t mature all died, and none accessed the T3 anergic population . In contrast, wild type T2 cells that did not mature were most likely to become anergic T3 cells, and only a minor fraction died directly from the T2 population (δ23 and μ2 parameters in Table 1). In order to validate the different developmental fates from T2 compartment predicted by our model, one may isolate T2 cells from wild type vs. Lyn−/− mice, transfer them into wild type recipient mice (distinguished with the Ly5 allelic difference), and measure what fraction of the transferred cells turn into Fo vs. T3 B cells. A more complicated experiment may be possible if a way is found to conditionally delete the Lyn gene at various developmental stages, as then we could follow the developmental fates of the progeny of the Lyn-deleted cells.
The most dramatic change in B cell population dynamics as indicated by our mathematical model was in the fate of those follicular mature B cells that were turning over moderately fast. In wild type mice, these less stable follicular B cells shifted to the anergic T3 phenotype and then eventually died from this population. In contrast in Lyn−/− mice, the less stable follicular B cells died directly. Interestingly, the numbers of follicular B cells in Lyn−/− mice are increased by approximately 10-fold by expression of a B cell-specific Bcl-2 transgene or by knockout of the pro-apoptotic BH3-only factor Bim (Gross et al, manuscript submitted), which is consistent with a high death rate contributing to the decreased numbers of follicular B cells in Lyn−/− mice.
We previously found that the elevated BCR signaling in B cells in Lyn−/− mice is most profound in T3 and follicular B cell populations, and is elevated to a lesser degree in T1 and T2 transitional B cells . Thus, an attractive hypothesis to explain the large increase in the death rate of follicular B cells resulting from Lyn deficiency is that the dying cells are those cells with a low degree of self-reactivity that, combined with strongly elevated BCR signaling, leads to clonal deletion, whereas in wild type mice, this degree of self-reactivity instead leads to clonal anergy and acquisition of the T3 phenotype. This hypothesis may be validated by using Ig transgenic mice with a weak self-reactivity, crossed with Lyn-deficient mice, and examining how the developmental kinetics would differ between Lyn−/− and wild type mice with this BCR transgene.
According to this hypothesis, T2 cells with a moderate degree of self-reactivity die directly and those with a lower degree of self-reactivity mature to the follicular B cell type. In wild type B cells, maturation to the follicular stage is accompanied by multiple developmental changes in signaling that balance enhanced positive signaling with increased activity of the inhibitory Lyn-CD22-SHP1 pathway, resulting in a threshold behavior in which low levels of BCR signaling are dampened out, but the B cell can still respond to foreign antigen in sufficient amounts . Thus, wild type follicular B cells can either adapt to a low level of self-reactivity or can enter the anergic population, depending on the degree of self-reactivity. In Lyn−/− B cells, however, the developmental changes increasing BCR signaling are unopposed by SHP1 and so a low degree of self-reactivity induces excessive BCR signaling leading to clonal deletion. A related hypothesis may explain the low access of Lyn-deficient T2 cells to the T3 anergic population, that is, as the cells begin to transition to the T3 stage, their signaling becomes magnified and this leads to rapid death before the T3 phenotype is fully realized. According to this view, access to the T3 population directly from the T1 population (δ13 in Table 1), rather than from the T2 or Fo subsets, may be possible in Lyn−/− mice because elevated BCR signaling is only manifest in later stages, giving more time for other negative regulatory mechanisms to balance off BCR signaling induced by a low degree of self-reactivity.
It should be noted that Lyn−/− mice have elevated levels of BAFF , a TNF family cytokine which promotes survival of T2, T3 and follicular B cells and may promote maturation of transitional B cells to mature B cells . It may be that the increased death rates of Lyn−/− T2 and follicular B cells would be further increased if BAFF levels were not elevated in these mice. Elevated BAFF levels have previously been shown to support survival of anergic B cells  and thus the increased BAFF levels in Lyn−/− mice may explain the decreased death rate we observed in the T3 population (parameter μ3 in Table 1).
Interestingly, the percent of "static" cells – cells that did not become labelled during the 8 days of the experiments – was higher in Lyn-deficient mice. The low degree of self-reactivity that is sufficient for Lyn-deficient cells to mature may thus not be sufficient to maintain a normal homeostatic division rate in this population. In summary, our previously derived mathematical model of B cell population dynamics was applied to analysis of the alterations in B cell maturation and survival in Lyn-deficient mice and was found to apply well to this substantial perturbation of B cells and to provide insights into the nature of those perturbations. The maturation rate and death rate parameters determined by the model were generally in accord with experimental observations regarding changes in the importance of negative regulation by Lyn as B cells mature and suggested reasonable hypotheses regarding how such changes might impact mechanisms for tolerizing self-reactive B cells in the periphery. These results provide insights into how alterations in BCR signaling affect B cell population dynamics and tolerance induction, although how tolerance breaks down in Lyn−/− mice to give spontaneous lupus-like autoimmunity remains to be elucidated.
Mice 7–12 weeks old were used for most experiments. C57BL/6 wild-type mice were obtained from Jackson Laboratory (Bar Harbor, ME, USA). Lyn−/− mice were used as described , and backcrossed at least fifteen generations onto C57/BL6 background. All animals were housed in a specific pathogen-free facility at UCSF according to University and National Institutes of Health (NIH) guidelines. Animal use was approved by the UCSF Institutional Animal Care and Use Committee.
Fluorophore conjugated antibodies directed against the following molecules were used: CD19 (1D3), CD23 (B3B4), IgM (II/41), BrdU, all from BD Pharmingen; CD93 (AA4.1) from eBioscience. Cells were analyzed on an LSR-II (BD Pharmingen). All FACS data were analyzed with FlowJo v. 6.4.1 (TreeStar software).
0.5mg of BrdU was injected into the peritoneal cavity of Lyn−/− and wild-type mice every 12 hours. Mice were euthanized at the specified time period, and 2×106 cells from bone marrow or spleen were permeabilized and stained for DNA with BrdU incorporation as well as cell surface markers according to the manufacturer's instructions (BD Pharmingen, FITC BrdU Flow Kit).
We first tested whether there are significant differences in BrdU labeling kinetics between Lyn-deficient and WT mice, using GLM repeated measures, a method based on ANOVA. Each dependent variable is represented by as many variables as there are measurement times (e.g., T1_2, T1_3, … are the BrdU fractions of labeled T1 cells after 2 days, 3 days, etc). Mouse type is the categorical predictor variable.
The modeled subsets were divided into labeled and unlabeled B cell populations, and the dynamics of these subsets are given in equations 7–18 below. In these equations UBM1, UBM2, UT1, UT2, UT3 and UM denote the numbers of unlabeled BM1, BM2, T1, T2, T3 and mature Fo B cells, respectively, and LBM1, LBM2, LT1, LT2, UT3 and LM similarly denote the labeled subsets. In these equations, BM1=UBM1+LBM1, BM2=UBM2+LBM2, T1=UT1+LT1, T2=UT2+LT2, T3=UT3+LT3, and M=UM+LM. Cells in each unlabeled compartment move to the corresponding labeled subset upon dividing (as shown for the BM cell subsets in Figure 2B). When an unlabelled cell divides during the labeling period, it leave the unlabeled subpopulation and its two daughter cells both become labeled cells and join the labeled subpopulation, hence the number of unlabelled proliferating cells entering the labeled subset is multiplied by a factor of 2. We neglect BrdU toxicity (and thus assign the same death rates to labelled and unlabelled cells), as in previous work (unpublished) we found that when the BrdU experiment is very short – as in this case, only 8 days – the assumption holds.
The mathematical models were simulated and fitted to data using a C language program. The program receives as input the experimental data, and the ranges of parameter values within which the model should be run. The program divides each parameter range to very small intervals, thus providing a thorough coverage of the parameter space. This creates a set of millions of parameter combinations to be checked by the program. For each parameter value set, the program integrates the model equations as follows. The initial conditions are zero cells in all populations; labeling only starts after the populations have reached a steady state. After integration, the program first checks whether the total cell numbers in all populations are within the experimentally measured ranges. Runs in which this is not the case are discarded. For all other runs, the program records the value of the fit of the model to the data, and outputs the parameter values and the value of the fit (MLE, see below). This process was performed for each type of mice (Lyn KO and WT), and the acceptable parameter value ranges were determined using the AIC method (see below).
We used the maximum likelihood estimation (MLE) to determine the parameter values that maximize the probability (likelihood) of the data. The best fit parameter value set is defined as the set of parameter values that yields the Maximum Likelihood Estimates for all time points and for all experimental repetitions (mice). This is standard methodology and has been implemented as follows. If, for each time point i, (i=1,…, k), we define a random variable Xi normally distributed with variance σi, then the likelihood function is:
where xj (j=1,…,Ni) are the measurements obtained at the time point i, and ψ(ti) is the value of the solution of the differential equations at time ti.
Furthermore, we used "Akaike's Information Criterion" (AIC) to find sets of parameter values that have a >95% probability to give a fit as good as the best parameters set. The AIC score is defined as follows.
Where k is the number of parameters in the model and MLE is the maximized value of the likelihood function for the estimated model. In fact, since the number of data points was not exceedingly large, we used the corrected Akaike Information Criterion (AICc), which accounts not only for MLE and k, but also for sample size n:
Suppose AICc(A) is the score of the parameter value set with the MLEA, and AICc(B) is the score of the best set of parameters with MLEB. In this case, the difference between the AICc scores is given by equation 22, and has a negative value, since MLEA / MLEB <1.
The probability that we have chosen the correct set of parameters is then computed from equation 22 according to the following expression.
If we reject parameters sets which, in comparison with the set with the maximum MLE, have |ΔAIC|≤6, then according to equation 23 we are left only with parameter sets which have a probability of at least 0.95 to give as good a fit as the best set. Thus, this method defines an equivalent of the 95% confidence intervals for the parameter values, and we used it to find the differences between the parameter ranges of WT and Lyn−/− mice. These ranges were considered significantly different only where the parameter confidence intervals did not overlap.
The work was supported in parts by the following grants: Israel Science Foundation grant number 270/09, and a Human Frontiers Science Program Research Grant (to R.M.); NIH grant R01AI20038 (to ALD); NIH grant K08AI52249 and the Rosalind Russell Medical Research Center for Arthritis, San Francisco (to AJG).
Conflict of interests
The authors declare no financial or commercial conflict of interest.