The B‐lymphocyte lineage is a leading system for analyzing gene regulatory networks (GRNs) that orchestrate distinct cell fate transitions. Upon antigen recognition, B cells can diversify their immunoglobulin (Ig) repertoire via somatic hypermutation (SHM) and/or class switch DNA recombination (CSR) before differentiating into antibody‐secreting plasma cells. We construct a mathematical model for a GRN underlying this developmental dynamic. The intensity of signaling through the Ig receptor is shown to control the bimodal expression of a pivotal transcription factor, IRF‐4, which dictates B cell fate outcomes. Computational modeling coupled with experimental analysis supports a model of ‘kinetic control’, in which B cell developmental trajectories pass through an obligate transient state of variable duration that promotes diversification of the antibody repertoire by SHM/CSR in direct response to antigens. More generally, this network motif could be used to translate a morphogen gradient into developmental inductive events of varying time, thereby enabling the specification of distinct cell fates.
The generation of a diverse set of pathogen‐specific antibodies, with differing affinities and effector functions, by B lymphocytes is essential for efficient protection from many microorganisms. Antibody gene diversification in B cells is mediated by two molecular processes termed class‐switch recombination and somatic hypermutation (CSR/SHM) (F1A). The former enables the generation of antibodies with the same antigen‐binding specificity, but different effector domains, whereas the latter results in a repertoire of antibodies with a range of affinities for a given antigen containing the same effector domain. CSR/SHM occurs in antigen‐activated B cells before their terminal differentiation into plasma cells. The transcription factor IRF‐4 is required for CSR/SHM as well as plasma‐cell differentiation, with its highest levels of expression being necessary for the latter. IRF‐4 acts in the context of a network of regulators that include Blimp‐1, Pax5, Bach2 and Bcl‐6 (F1B). Despite extensive characterization of these individual factors, how the network responds to sensing of antigen by the B cell antigen receptor (BCR, antibody molecule expressed on cell surface) to regulate the extent of antibody gene diversification and plasma‐cell differentiation remains to be addressed.
To address this issue, we assemble a computational model. The model reveals two contrasting scenarios that can underlie B cell fate dynamics. In one case, the initial rate of IRF‐4 production controls a binary cell fate choice that involves either going to the CSR/SHM state or to the plasma‐cell state; the time spent in the CSR state is relatively insensitive to the initial rate of IRF‐4 production (herein called ‘basic bistability’). In the other case, IRF‐4 drives all cells through a transient CSR/SHM state, but the initial rate of IRF‐4 production sets its duration (‘kinetic control’). Both scenarios predict that increasing the initial rate of IRF‐4 production favors the generation of plasma cells at the expense of CSR/SHM, but they differ fundamentally with respect to the underlying gene expression patterns.
To distinguish between these two scenarios experimentally, we utilize two different genetic models. The first involves the B1‐8i transgenic mouse whose B cells express a rearranged V187.2 VDJ Ig heavy chain gene segment that is specific for the hapten nitrophenol (NP). The second is a newly developed mouse model that allows exogenous control of IRF‐4 expression in naive primary B cells using a tet‐inducible allele. Using these models, we show that (i) BCR signal strength sets the initial rate of IRF‐4 accumulation and (ii) the concentration of IRF‐4 is sensed by an incoherent gene regulatory network architecture to regulate the extent of CSR/SHM prior to plasma‐cell differentiation. Our results are consistent with the ‘kinetic‐control model’ in which the levels of BCR‐induced IRF‐4 expression control the duration of an obligate CSR/SHM state that enables B cell diversification before terminal differentiation into plasma cells. Evidence for the transient CSR/SHM state is corroborated by both patterns of gene expression and the presence of AID‐dependent mutations in individual non‐switched plasmablasts.
Our results provide a molecular framework for understanding how B cells balance the competing demands for Ig CSR and SHM with the secretion of antibodies during humoral immune responses. The key feature of the network architecture that allows IRF‐4 to coordinate the two competing states of gene expression in a temporal manner is that it simultaneously but asymmetrically activates both sides of a bistable mutual repression circuit. Because the two effects of the primary regulator antagonize each other, we describe the circuit as being based on an ‘incoherent’ regulatory motif. Other incoherent regulatory motifs in varied biological systems are also associated with the acquisition of transient cell states, and we consider how the kinetic‐control mechanism proposed by us could more generally serve to translate developmental cues into elaborate morphogenetic patterns.
The intensity of signaling through the B cell receptor controls the level of expression of IRF‐4, a transcription factor required for B cell differentiation. The rate of IRF‐4 production dictates the extent of antibody gene diversification that B cells undergo upon antigen encounter before differentiating into antibody‐secreting plasma cells.
Computational modeling and experimental analyses substantiate a model, whereby IRF‐4 regulates B cell fate trajectories via a ‘kinetic‐control’ mechanism.
Kinetic control is a process by which B cells pass through an obligate state of variable duration that sets the degree of cellular diversification prior to their terminal differentiation.
An incoherent regulatory network architecture, within which IRF‐4 is embedded, is the basis for realization of kinetic control.
Over the last decade, considerable progress has been made in the combined experimental and theoretical analyses of regulatory networks that are utilized by cells to process information about their environments and regulate their function (Alon, 2007). Much of this work has focused on gene regulatory networks (GRNs) in unicellular organisms (Aurell et al, 2002; Ozbudak et al, 2004; Santillan and Mackey, 2004; Acar et al, 2005; Suel et al, 2006). These studies have demonstrated that certain network topologies (motifs) are conserved and appear to have evolved to perform similar functions in a variety of contexts (Milo et al, 2002). The cells of multicellular organisms proceed through many developmental transitions and can acquire a multitude of differentiation states, each associated with a unique set of functions. The network topologies and dynamics underlying these more complex developmental programs are being intensively explored (Davidson, 2006). Understanding metazoan GRNs from an integrated theoretical and experimental perspective is of great interest given the elaboration of cellular states that exhibit highly complex spatial and temporal dynamics.
The mammalian immune system represents a powerful system for exploring GRNs that regulate cellular differentiation (Singh et al, 2005). These networks culminate in the acquisition of effector functions for host defense against pathogens. Immunoglobulin (Ig) molecules (antigen receptors) are the hallmark of B lymphocytes and have a central role in their effector functions. Whereas membrane Ig enables antigen uptake and its presentation to T cells, secreted Ig (antibody) functions in neutralization and phagocytosis of pathogens. Upon antigen encounter, B cells can undergo an exceptional developmental dynamic involving diversification of their antigen receptor (membrane Ig) by targeted somatic hypermutation (SHM) of Ig genes and/or expression of distinct Ig isotypes (IgM versus IgG, E, A) via class switch DNA recombination (CSR) (Rajewsky, 1996) (Figure 1A). SHM‐mediated diversification of the initial Ig repertoire enables affinity maturation by selection of mutations during the course of the humoral response. The molecular processes of SHM and CSR occur prior to terminal differentiation of B cells into Ig‐secreting plasma cells. It should be noted that whereas CSR can occur during B cell responses that do not require interaction with T cells, SHM occurs in the context of T cell‐dependent B cell responses within specialized structures termed germinal centers (MacLennan et al, 2003). Despite extensive characterization of B cell responses in the context of distinct immunization and infection regimes, the molecular circuits that regulate the extent and nature of Ig receptor diversification by SHM and/or CSR upon antigen encounter remain poorly understood.
CSR and SHM are dependent on the transient expression of AID (Muramatsu et al, 2000), an activation‐induced cytidine deaminase that targets the C residue, predominantly at ‘WRC’ sequences. Coordinated AID catalyzed deamination at donor and acceptor switch regions within the IgH locus followed by double‐strand DNA cleavage and non‐homologous end joining results in the functional juxtaposition of a given variable gene segment with a new IgH constant region. In SHM, targeted deamination of the Ig variable gene region by AID followed by error‐prone repair leads to the generation of new variants, some of which may confer higher‐affinity binding to antigen.
The state that enables CSR and SHM requires the transcription factors Pax5, Bcl‐6, Bach2 and IRF‐4 (Dent et al, 1997; Fukuda et al, 1997; Muto et al, 2010). Pax5 and IRF‐4 regulate the induction of the AID gene (Klein et al, 2006; Sciammas et al, 2006; Tran et al, 2010). Pax5 additionally positively regulates the expression of the Bach2 and Bcl‐6 genes, the products of which function to repress the plasma‐cell‐specific program (Shaffer et al, 2000; Tunyaplin et al, 2004; Nera et al, 2006; Ochiai et al, 2006; Schebesta et al, 2007). Whereas Bach2 is required within extrafollicular as well as germinal center B cells, Bcl‐6 functions specifically in the latter context. Unlike Pax5, Bcl‐6 and Bach2, the transcription factor IRF‐4 is also required for the plasma‐cell state (Klein et al, 2006; Sciammas et al, 2006). This state is additionally dependent on a distinct set of regulators that include Blimp‐1 and Xbp1 (Shapiro‐Shelef et al, 2003; Sciammas and Davis, 2004; Shaffer et al, 2004). In addition to promoting plasma‐cell differentiation, Blimp‐1 represses key components of the CSR/SHM program (Pax5 and Bcl‐6), thereby creating a mutual antagonism (Lin et al, 2002; Shaffer et al, 2002; Kallies et al, 2004; Sciammas and Davis, 2004; Cimmino et al, 2008). Although these transcription factors have been individually characterized and also assembled into a GRN (Sciammas et al, 2006) (Figure 1B), it remains to be determined how their collective interactions give rise to distinct developmental trajectories and cell fates.
As noted above, IRF‐4 promotes two mutually antagonistic programs of B cell gene expression that underlie CSR/SHM or plasma‐cell differentiation (Klein et al, 2006; Sciammas et al, 2006). IRF‐4 not only positively regulates the AID gene, but also the gene encoding Blimp‐1 (Sciammas et al, 2006). We proposed that differing concentrations of IRF‐4 regulate the transition between the competing programs of gene expression with low levels favoring CSR/SHM and higher levels inducing plasma‐cell differentiation (Sciammas et al, 2006). This remains to be experimentally tested. While our earlier results established the importance of IRF‐4 for both fates, it remained to be elucidated how IRF‐4 enables the GRN to modulate its response to signaling through the antigen or other receptors, so as to control the extent and nature of Ig diversification before differentiation of B cells into Ig‐secreting plasma cells. We now use computational modeling to analyze the network architecture and then test its predictions experimentally.
We demonstrate that B cell antigen receptor (BCR) signaling initially induces IRF‐4 expression in an antigen‐dependent manner, and this pattern resolves into a bimodal distribution of cells with low and high IRF‐4 levels. These differing levels of IRF‐4 promote B cell differentiation into Ig‐switched and Ig‐secreting cells, respectively. Computational modeling is used to contrast two scenarios that could underlie B cell fate dynamics. In one case, the initial IRF‐4 production rate controls a binary choice between going to the CSR/SHM state or to the plasma‐cell state; the time spent in the CSR state is relatively insensitive to the initial IRF‐4 production rate (herein called ‘basic bistability’). In the other case, IRF‐4 drives all cells through a transient CSR/SHM state, but the initial IRF‐4 production rate sets its duration (‘kinetic control’). To distinguish between these two scenarios experimentally, we control IRF‐4 expression levels via differential BCR signaling or by using a novel inducible transgenic Irf4 allele. Gene expression patterns corresponding with the two competing developmental states as well as cellular fate outcomes are quantitatively analyzed in response to varying IRF‐4 concentrations. Importantly, we quantify the fraction of B cells that pass through the CSR/SHM state by determining the frequency of AID‐dependent mutations experienced by individual cells. Collectively, our data are consistent with the kinetic‐control model in which the initial IRF‐4 production rate dictates the duration of an obligate CSR/SHM state that activated B cells transit through before differentiating into Ig‐secreting plasma cells. These results imply that IRF‐4 serves as a sensor of antigen receptor signaling to control the duration of CSR/SHM and promote the exit of high‐affinity B cells from the germinal center. We propose that the novel network motif could be used in other developmental contexts for translating a graded inductive signal into discrete temporally regulated programs of gene expression, thereby specifying cell fates.
Architecture of GRN that regulates effector B cell fate choice
The GRN that underlies the generation of effector B cells (Figure 1B) centers on the mutual repression between transcription factors of the plasma‐cell program (Blimp‐1 and IRF‐4hi) (Shaffer et al, 2002; Kallies et al, 2004; Sciammas and Davis, 2004; Sciammas et al, 2006; Teng et al, 2007) and those needed for CSR/SHM (Pax5, Bcl‐6, Bach2 and IRF‐4lo) (Shaffer et al, 2000; Tunyaplin et al, 2004; Nera et al, 2006; Ochiai et al, 2006; Schebesta et al, 2007). Mutual repression is a common strategy underlying exclusive realization of competing binary outcomes (Laslo et al, 2006; Alon, 2007). However, there are two crucial differences in the architecture of the GRN here in comparison with ones that were previously studied: (i) a single‐regulator (IRF‐4) activates genes on both sides of the genetic switch and (ii) in addition to the mutual repression, there is a positive feedback loop based on mutual activation (Blimp‐1 and IRF‐4; Figure 1B) (Kallies et al, 2004; Sciammas and Davis, 2004).
In this regard, it is important to note that the key positive regulatory connection between IRF‐4 and Blimp‐1 (Sciammas et al, 2006) was not supported by another study (Klein et al, 2006). To corroborate our earlier findings and validate the GRN, we crossed a Blimp‐1:GFP knock‐in reporter allele (Kallies et al, 2004) with the Irf4−/− alleles. Splenic B cells from these mice were stimulated with LPS or CD40L in the presence of IL‐4. These conditions efficiently induced Blimp‐1 expression in a large fraction of the cells and their differentiation into Ig‐secreting plasma cells. Whereas the Blimp‐1:GFP allele was strongly induced in Irf4+/− B cells, its activation was considerably impaired in Irf4–/– B cells (Supplementary Figure S1). Thus, in combination with our previous analyses and subsequent reports in the literature (Shaffer et al, 2008), these data demonstrate that IRF‐4 positively regulates Blimp‐1 gene expression during plasma‐cell differentiation.
Mathematical modeling of B cell fate trajectories
The GRN interactions do not in themselves explain the observed IRF‐4 expression dynamics and the control of B cell fate. To explore these issues, we translated the network into a set of rate laws for the expression of the molecular species in Figure 1B by standard means (see Materials and methods; Supplementary information). Activation is treated as catalysis of gene expression and repression as (non‐competitive) inhibition, and corresponding Michaelis–Menten‐like forms are used to represent the production of each species. We then add degradation terms and noise, and integrate the resulting equations to simulate the gene regulatory dynamics of activated B cells (here, by ‘activated B cells’ we mean cells that have been stimulated by antigen recognition through the BCR, but have not yet transitioned to either the CSR/SHM or plasma‐cell state). In an activated B cell, AID and Blimp‐1 are not expressed to measurable extents, and their initial concentrations and expression rates are thus set to zero in the model (marked ‘initiation’ in Figure 1C and D); Bach2 and Pax5 are expressed, and we set the initial rate of Pax5 production to be such that its concentration rises upon activation (Wakatsuki et al, 1994) (data not shown). As Irf4 is an immediate early gene downstream of antigen receptor signaling (Matsuyama et al, 1995), and its level of expression correlates with cell fate (Sciammas et al, 2006), we examined the dynamics of the model as a function of the initial IRF‐4 production rate (the initial concentration of IRF‐4 is set to zero).
In general, the kinetic parameters for the various elementary steps contributing to the gene regulatory dynamics have not been measured. To limit the free parameters to a number that can be exhaustively explored, we set the maximal rates of activated expression (kx in Equation (1)), the rates of protein degradation (b) and the binding affinities of the regulators (Ka=Kr) to be the same for as many species as allowed capturing the known biological features. We then verified that our conclusions were robust to non‐uniform choices of the parameters by randomizing them. The model revealed two possibilities regarding the manner in which increasing IRF‐4 production could regulate cell fate. In the first, IRF‐4 production controls an initial decision between the plasma‐cell and CSR/SHM states (Figure 1C). Greater production of IRF‐4 favors the plasma‐cell state. Importantly, exit from the CSR/SHM state in this scenario is relatively insensitive to the initial IRF‐4 production rate (Figure 1E); the plasma‐cell state is irreversible owing to the positive feedback loop between Blimp‐1 and IRF‐4. Thus, in this basic bistability scenario, IRF‐4 production modulates the cell fate decision, while having little effect on its kinetics. In the second possibility, the kinetic‐control scenario, all cells proceed through the CSR/SHM state with variable duration before differentiating into plasma cells (Figure 1D). The main effect of raising IRF‐4 production is to reduce the average time that cells dwell in the CSR/SHM state (Figure 1F). In this scenario, increases in the initial IRF‐4 production rate diminish the likelihood that B cells have sufficient time to undergo CSR/SHM. The mathematical basis for the two scenarios is discussed further in Box 1.
Box 1 Mathematical basis for the basic bistability and kinetic‐control scenarios
Here, we qualitatively illustrate how different behaviors emerge from the GRN in Figure 1B depending on the choice of parameters. To emphasize the generic nature of the behaviors, we use a simplified GRN in which an inductive signal or pivotal regulator (X) activates the expression of two counter‐acting components (Y and Z; see Warmflash and Dinner (2009) for a more quantitative discussion of this motif). As indicated in the figure, different rates of X production result in different developmental trajectories as follows: (i) At low rates, the system is bistable: there are two stable steady states (red and green circles) and the boundary (or tipping point) separating their basins of attraction is determined by an unstable steady state (blue). The initial conditions indicated at the bottom left of the YZ plane are in the basin of attraction of the green steady state, and trajectories flow to the green steady state in which Y represses Z. Initial conditions further to the right (higher Z, not shown) could lead directly to the red steady state. (ii) The green stable steady state and the blue unstable steady‐state coalesce as X production increases; trajectories starting at the bottom left still flow to the green state, but small perturbations cause them to continue to the red state. (iii) Trajectories follow the same flow pattern, but now there is a route out of the region where the green and blue steady states came together (the ‘ghost’ of the stable state (Strogatz, 1994), dashed circle); because the direction of the flow previously reversed at the blue steady state, the dynamics are slow in this region. All flows eventually lead to the red steady state, in which Z represses Y. (iv) Qualitatively similar to (iii), but now the previous steady states have less impact on cell fate dynamics (X production has changed more) and flow to the red steady state is faster; Z more strongly represses Y.Different choices of parameters cause the threshold point in (ii), known as a reverse saddle‐node bifurcation, to fall at different rates of X production relative to the range that is accessed experimentally. In the basic bistability scenario, the bifurcation point is high compared with the range of rates of X production, such that the behavior in (i) is always operative. This, together with choices of initial conditions that are poised at the boundary between the basins of attraction of the two steady states, allows early partitioning of the developmental trajectories. In the kinetic‐control scenario, the bifurcation point is in or just below the range of rates of X production; as a result, the green stable state is lost, but trajectories starting in the bottom left are still attracted to the corresponding pattern of gene expression and spend a finite time in its vicinity as in (iii). A robust feature is that the time required to pass through the region of the transient state (which here permits AID expression) decreases significantly as the controlling parameter (here, the initial IRF‐4 production rate) is raised further (iv). It is of course also theoretically possible that the bifurcation is well below the experimental range of rates of X production, but we exclude this third scenario from further consideration as it implies in the present study that cells never have an opportunity for CSR, which is inconsistent with the biology. Further analysis (Materials and methods, Supplementary information, Supplementary Figures S7–9, and data not shown) shows that, because the basic bistability and kinetic‐control scenarios differ primarily with respect to their early dynamics, the parameters that most strongly impact whether one or the other is realized are those that strengthen the regulation of Bcl‐6/Bach2 by Pax5 (s and nC in Equation (1d), as well as the initial levels of Bcl‐6/Bach2 and Pax5).
Gene regulatory and cell fate dynamics for a prototypic incoherent activation structure. Colors correspond to those used throughout the main text to indicate the CSR/SHM and plasma‐cell states.
In general, measurements on populations of cells can reflect not only differentiation but also proliferation. To account for both factors in making experimentally testable predictions, we developed a multiscale simulation. We consider a population of proliferating cells, in each of which the concentrations of the molecules shown in Figure 1B evolve independently according to the model described above. We assume that antigen is in equilibrium with the BCR to set the extent of ligand binding, and, in turn, the mean initial production rate of IRF‐4 (normally distributed). Based on the concentrations of AID and Blimp‐1, cells probabilistically undergo CSR/SHM and/or differentiation into plasma cells (see Supplementary information). Simulations based on either the basic bistability or kinetic‐control scenarios both predict that increasing the initial IRF‐4 production rate favors the generation of plasma cells at the expense of CSR/SHM (Figure 2A and B).
To discriminate between the basic bistability and kinetic‐control scenarios, we used the multiscale model to predict the mRNA dynamics for the genes in the network as a function of time and the initial IRF‐4 production rate (Supplementary Figure S6). In the basic bistability scenario, higher IRF‐4 production increases the likelihood that a cell transitions to the plasma‐cell state without visiting the CSR state. In contrast, for the kinetic‐control scenario, IRF‐4 production governs the time a cell spends in the CSR/SHM state. Thus, the two scenarios can be distinguished by examining the peak mRNA expression levels of the antagonistic regulators (Bcl‐6/Bach2 and Blimp‐1) as a function of the initial IRF‐4 production rate (Figure 2C and D). In the basic bistability scenario, the curve for Bcl‐6/Bach2 expression trends down and that for Blimp‐1 trends up (Figure 2C) because an increase in the production rate of IRF‐4 favors direct differentiation into plasma cells without visiting the CSR/SHM state. In contrast, for the kinetic‐control scenario, the expression curve for Bcl‐6/Bach2 is flat (Figure 2D) because changes in IRF‐4 production impact the duration rather than the peak levels of gene expression associated with the CSR/SHM state (see Supplementary Figure S8 for additional analysis of the sensitivity of the Bcl‐6/Bach2 gene expression curve with the specific choice of parameters). In summary, although both scenarios predict that increasing the initial IRF‐4 production rate favors the generation of plasma cells at the expense of CSR/SHM, they differ fundamentally with respect to their transcriptional dynamics.
Strength of BCR signaling determines the dynamics of IRF‐4 expression
With a view toward testing the mathematical modeling and distinguishing between the above scenarios, we utilized B cells from the B1‐8i transgenic mouse, which expresses a rearranged V187.2 VDJ gene segment knocked in to the Ig heavy chain locus (Sonoda et al, 1997). This rearranged IgH allele increases the precursor frequency of B cells that respond to nitrophenol (NP)‐conjugated antigens. Splenic B1‐8i B cells were activated in vitro with NP(12)‐Ficoll antigen; IRF‐4 expression levels and cell division (CFSE dilution) were monitored over time using flow cytometry (Figure 3A). IRF‐4 was induced as early as 12 h following antigen stimulation (data not shown) and exhibited a dose‐dependent pattern of expression in undivided cells (see day 1 distribution). Cell division was accompanied by the generation of two stable IRF‐4‐expressing populations, referred to as IRF‐4lo and IRF‐4hi (see day 3 distribution).
We next measured IRF‐4 levels in antigen‐stimulated B1‐8i B cells as a function of antigen dose and avidity [NP(12)‐versus NP(109)‐Ficoll] (Figure 3B). Notably, at day 1 of stimulation, increased antigen concentration shifted the IRF‐4 distribution toward the higher expression state (Figure 3B; Supplementary Figure S2). At later times (e.g. day 3), the distribution of cells in the IRF‐4lo and IRF‐4hi states was similarly altered as a function of antigen concentration. Since we had previously correlated IRF‐4hi‐expressing cells with expression of the plasma‐cell marker Syndecan‐1 (Sciammas et al, 2006), these data suggest that increased BCR signaling and, in turn, increased initial IRF‐4 expression correlates with enhanced plasma‐cell generation (see below for a more detailed analysis). A comparison of the dose response of NP(12)‐Ficoll with NP(109)‐Ficoll revealed that the higher avidity antigen shifted the IRF‐4 distribution toward the higher expression state at lower concentrations (data not shown; Figure 3B). In summary, these data show that BCR signal strength sets the initial pattern of IRF‐4 expression, which in turn dictates the frequencies of IRF‐4lo and IRF‐4hi cells.
BCR signal strength determines cell fate outcome
Having established the relationship between BCR signaling and IRF‐4 expression levels, we sought to validate the prediction common to both scenarios that increasing the initial IRF‐4 production rate favors the generation of plasma cells at the expense of CSR. To this end, we analyzed the outcome of B1‐8i B cell differentiation as a function of antigen dose or avidity. Purified splenic B1‐8i cells were stimulated with antigen, CD40 ligand and the cytokines interleukin‐2, ‐4 and ‐5 and their differentiation was analyzed by comparing the frequency of plasmablast (Syndecan‐1‐expressing) or switched (surface IgG1‐expressing) cells by flow cytometry (Figure 4A) and ELISpot (Supplementary Figure S3). Strikingly, higher antigen concentrations led to an increase in the frequency of Syndecan‐1‐expressing plasma cells, and that increase occurred at the expense of CSR to IgG1. As predicted, stimulation of B1‐8i cells with higher avidity antigen (NP(109)‐Ficoll) resulted in increased generation of plasmablasts and a corresponding decrease in the frequency of cells that had undergone CSR (Figure 4A). Thus, as predicted by the modeling (Figure 2A and B), the dynamic patterns of IRF‐4 expression, which vary with the intensity of BCR signaling, correlate with the modulation of B cell fate choice.
Care must be taken in interpreting cell frequency data in the context of a population of cells undergoing division; therefore, we enumerated absolute numbers of plasmablasts as well as cells that had undergone CSR to IgG1 as a function of NP(12)‐Ficoll dose on day 3. Consistent with the frequency data, there is a dramatic increase in the absolute numbers of plasma cells and a decrease in CSR (Supplementary Figure S3). We note, however, that the increase in plasma‐cell numbers is steeper than the decline in cells that have undergone CSR to IgG1, suggesting that IRF‐4 may exert effects on cell growth in addition to differentiation. To address this possibility, we quantitatively analyzed cell division rates (Supplementary information; Supplementary Table S2). Although the rate of cell division of IRF‐4hi cells was two‐ to three‐fold higher than that of their IRF‐4lo counterparts, simulations based on these rates yielded results similar to those in Figure 2A and B (Supplementary Figure S10). Taken together, these analyses support the idea that the GRN in Figure 1B functions to translate BCR signal strength into cell fate choice.
Gene expression dynamics support the kinetic‐control scenario
While the two scenarios suggested by the modeling lead to similar predictions for cellular fates, they differ in their predictions for the expression patterns of the components of the network (Figure 2). Thus, we analyzed the RNA expression patterns of the network components in B1‐8i cells by quantitative RT–PCR as a function of time and antigen concentration (Figure 4B) and avidity (data not shown). Consistent with the outcome of differentiation, increasing BCR signal strength results in a strong upregulation of the plasma‐cell program component Blimp‐1. With regard to the model predictions, we see that AID is induced to similar levels at all BCR signal intensities (see day 2 expression levels), but then decreases at day 3 in a manner strongly dependent on BCR signal intensity. These dynamics reflect uniform peak levels of CSR regulators (Bach2 and Pax5 at day 1) over a wide range of NP concentrations (Figure 4B and C) and are consistent with the kinetic‐control scenario, but not the basic bistability scenario (compare Figure 4C with Figure 2C and D). In other words, the CSR/SHM pattern of gene expression appears to reflect an obligate developmental state, and the role of IRF‐4 is to control its duration.
We additionally fit the core gene regulatory model (Equation (1)) to the mRNA data, as detailed in the Supplementary information. Given the minimal nature of the model, the expression patterns are reproduced with reasonable accuracy (Supplementary Figure S11). Although the data are not sufficient to constrain the parameters to unique values, in all cases, the resulting bifurcation points, which determine the model behavior as described in Box 1, are within the range of signals assumed, as in the kinetic control, but not the basic bistability scenario. This quantitative result, which is based on the complete mRNA data rather than only the peak levels, corroborates the qualitative predictions in Figure 2.
Finally, we compared the stochastic multiscale simulation behavior with the flow cytometry data in Figure 3B (Supplementary Figure S12). In both scenarios, the variations in signal are translated to shifts in the IRF‐4lo and IRF‐4hi populations by the positive feedback resulting from mutual repression between Bcl‐6/Bach2 and Blimp‐1 and mutual activation between IRF‐4 and Blimp‐1 (Figure 1B). However, there are significant numbers of cells with IRF‐4 levels between the IRF‐4lo and IRF‐4hi populations in Figure 3B, and this pattern is again better captured by the kinetic‐control scenario (Supplementary Figure S12). This scenario exhibits a greater number of such cells because they correspond to cells that pass through the CSR/SHM state rather than directly differentiating to the plasma‐cell state. We thus interpret the data in Figure 3B to lend support to the idea that a substantial fraction of cells pass through the CSR/SHM state, and we sought a means to examine this proposition directly.
Mutational analysis of IgH switch regions demonstrates that IgM‐secreting plasma cells transiently visit a CSR/SHM state
A complementary means of evaluating the model scenarios would be to quantify the fraction of cells that pass through the CSR/SHM state. We focus our analysis on AID‐dependent mutations associated with CSR as under the in vitro conditions, the frequency of SHM is low. Clearly, all IgG1+ cells pass through the CSR/SHM state; the question is how many IgM+ plasmablasts cells do as well. In a bistability scenario in which CSR/SHM is effectively coupled to AID induction, this fraction would be negligible, whereas in the kinetic‐control scenario with probabilistic opportunity for CSR/SHM, it would be appreciable and only limited by the approach used to detect the state. We thus sought to quantify the fraction of IgM+ plasmablasts cells that pass through the CSR/SHM state. The Sμ region is known to accumulate AID‐dependent mutations following B cell activation (Petersen et al, 2001; Nagaoka et al, 2002; Xue et al, 2006). We purified wild‐type B cells that were stimulated with CD40L and IL‐4, and sorted them based on IgM and Syndecan‐1 or IgG1 expression to screen for terminally differentiated cells that induced AID and had or had not undergone functional CSR (Supplementary Figure S4). Genomic DNA was purified from the sorted cells and used to generate PCR products of the Sμ region for subsequent cloning and sequencing.
Analysis of the sequences from individual IgM+, Syndecan‐1‐expressing plasmablasts revealed the presence of mutations in the Sμ region, the majority of which were C to T transitions involving the AID hotspot motif WRC (Table I; Supplementary Figure S4). Importantly, this spectrum of mutations was dependent on AID, as similarly sorted cells from Aicda−/− cells exhibited no such mutations. The background mutations likely reflect Taq polymerase error. A total of 7% of IgM+ plasmablasts have AID‐dependent mutations; given that only 21% of IgG1+ cells have AID‐dependent mutations, we estimate the fraction of IgM+ cells passing through the CSR/SHM state to be 33% (7% normalized by 21%). This is a lower bound since we expect that sorting for cells that have not switched experimentally selects for a population that expressed AID for a shorter time. These single‐cell data substantiate the kinetic‐control scenario in the sense that a large fraction of IgM+ plasma cells (cells that never completed CSR) followed a developmental trajectory through a state of transient AID activity accompanied by a probabilistic opportunity for CSR/SHM prior to terminal differentiation.
Induced levels of IRF‐4 dictate B cell fate dynamics
To directly test a causal role for Irf4 expression levels in controlling B cell fate, we developed an experimental system that would afford control of Irf4 expression in naive primary B cells. We constructed a tet‐inducible system (Beard et al, 2006) as follows; a cassette containing a tetracycline‐regulatable Irf4 transgene was recombined at a neutral, defined position downstream of the CollagenA1 locus using Flp‐mediated homologous recombination in ES cells (Supplementary Figure S5). These cells were engineered to express the tetracycline‐dependent transactivator (M2rtTA) from the Rosa26 locus. Suitable ES cell clones (Supplementary Figure S5) were used to transmit the tet‐regulatable Irf4 allele through the germline. When crossed to the Irf4+/− background, both endogenous and transgene expression of Irf4 combine to attain higher than wild‐type expression levels. In contrast, when crossed to the Irf4−/− background, Irf4 gene expression can be uncoupled from BCR as well as CD40 signaling, and a lower range of expression of IRF‐4 can be induced and maintained because the transgene is the sole source. In this way, the inducible‐Irf4 allele enables a wide dynamic range of IRF‐4 production.
We first tested the dynamics of alternate cell fate choice as a consequence of inducing IRF‐4 expression in cells containing an endogenous wild‐type Irf4 allele. We purified B cells from the Irf4‐inducible mouse (on the Irf4+/− background) and stimulated them with CD40L and IL‐4 in the presence of varying concentrations of a high‐affinity analog of tetracycline, doxycycline (DOX). Strikingly, increasing Irf4 expression levels induced by DOX resulted in markedly enhanced plasma‐cell generation that occurred at the expense of CSR (Figure 5A; compare with Supplementary Figure S12). The lower panel shows the amount of IRF‐4 expression achieved by DOX at day 3 of differentiation. We next measured the effects of exogenously controlled Irf4 expression on the dynamics of network components (Figure 5B and C). At early time points, Pax5, Bach2 and AID are all induced to comparable levels regardless of DOX concentration, consistent with the kinetic‐control scenario. Downregulation of the CSR factors occurs only as Blimp‐1 levels rise at later time points.
We also tested the effects of Irf4 levels on B cell fate dynamics in the absence of endogenous Irf4 expression. DOX‐inducible IRF‐4 was capable of complementing the Irf4 genetic deficiency (Figure 6A; compare with Supplementary Figure S13). In the absence of endogenous Irf4, we note that the extent of CSR initially increases and then decreases with dose. In the model, the increase results from the fact that IRF‐4 is needed for cells to initiate CSR, while the decrease results from the kinetic‐control mechanism. We note that the magnitude of the decrease in CSR frequency is less in Irf4−/− cells than Irf4+/− cells, likely due to subtle differences in Irf4 expression levels that were achieved in the two experimental conditions. The expression dynamics mirror those obtained in the presence of endogenous Irf4, although the fold change in Blimp‐1 at later times is much more dramatic because IRF‐4 is required for its expression. Importantly, the peak expression of CSR/SHM‐associated genes is the same across DOX conditions (Figure 6B and C). Collectively, these data confirm that varying IRF‐4 levels is sufficient to control the duration of the CSR/SHM state and in turn the cell fate dynamics of activated B cells, and the associated gene expression patterns provide additional support for the kinetic‐control scenario.
Our results strongly suggest that the levels of BCR‐induced IRF‐4 expression control the duration of an obligate CSR/SHM state that enables B cell diversification before their terminal differentiation into plasma cells. We have shown that (i) BCR signal strength sets the initial IRF‐4 concentration and (ii) the GRN architecture enables sensing IRF‐4 concentration to regulate the extent of CSR prior to plasma‐cell differentiation. Evidence for the transient CSR/SHM state is corroborated by both patterns of gene expression and the presence of AID‐dependent mutations in non‐switched plasmablasts. We note that our mutation data are consistent with recent in vivo observations that IgM‐expressing memory cells have experienced AID activity (Dogan et al, 2009). Given that Bach2 expression was the most sensitive discriminator between the basic bistability and the kinetic‐control scenarios, it is interesting to note that Muto et al (2010) recently showed that Bach2 promotes CSR indirectly by repressing Blimp‐1, and they proposed that this regulatory interaction functions as a time delay to enable CSR/SHM before plasma‐cell differentiation.
Given that AID is required for CSR/SHM and is regulated by IRF‐4, we propose that BCR signal strength sets the duration of SHM in vivo by dictating the induced concentration of IRF‐4. Thus, a low‐affinity BCR would induce a lower concentration of IRF‐4 upon antigen binding and sustain AID expression for a longer duration, thereby increasing the likelihood of acquiring mutations in variable regions of Ig loci that confer higher affinity. The higher‐affinity BCR would then induce higher levels of IRF‐4 expression and a qualitative change in the dynamics of the GRN that would terminate AID expression and the germinal center program and cause the activated B cell to undergo plasma‐cell differentiation (Figure 1A). This proposal does not preclude the possibility that there is an additional checkpoint, whereby follicular T helper cells ‘license’ high‐affinity B cells to exit the germinal center (Allen et al, 2007); in this scenario, a T cell‐derived signal could be integrated with one emanating from the BCR to regulate overall IRF‐4 concentration.
Our findings suggest a molecular explanation for the changing nature of B cell responses in BCR transgenic mice immunized with cognate antigens of differing affinity or avidity (O'Connor et al, 2006; Paus et al, 2006; Phan et al, 2006; Fink et al, 2007). In these experiments, increasing antigen affinity promotes the robust generation of IgM‐secreting plasma cells. The germinal center response is either diminished or unaffected. We propose that these cell fate dynamics result from strongly responding B cells that rapidly accumulate high levels of IRF‐4 and differentiate into IgM‐secreting plasma cells. Weakly responding B cells that induce lower levels of IRF‐4 may preferentially undergo a germinal center response. Consistent with our proposed role for IRF‐4, it was recently shown that during a germinal center response, high‐affinity transgenic B cells recognizing hen egg lysozyme are actively selected to become plasma cells (Phan et al, 2006). While earlier studies (Shih et al, 2002a, 2002b) reached a different conclusion to that of Brink and co‐workers (O'Connor et al, 2006; Paus et al, 2006; Phan et al, 2006; Fink et al, 2007), the former analysis was performed over a narrower range of BCR affinity (Tarlinton, 2008). Our results provide a molecular framework for understanding B cell fate dynamics, which balance the competing demands for Ig CSR/SHM with the secretion of antibodies during humoral immune responses.
In fact, extension of our cellular multiscale stochastic simulations to model in vivo plasma‐cell and germinal center responses (Supplementary information; Supplementary Figure S14) shows that the network dynamics described here are sufficient to account for the observations reported in immunized BCR transgenic mice (O'Connor et al, 2006; Paus et al, 2006) (Supplementary Figure S15). We note that while models of B cell fate dynamics (Kleinstein and Singh, 2001; Hawkins et al, 2007; Figge et al, 2008) were developed previously, to the best of our knowledge, they were wholly cellular in nature and did not link dynamic behavior to an underlying GRN, as we do here. By evaluating the cell autonomous capacity for differentiation, we have derived a ‘core’ module that operates, likely in cooperation with additional factors in physiological contexts (such as differential co‐receptor signaling and cell competition), to control effector B cell fate dynamics.
The ability of IRF‐4 to coordinate the two competing patterns of gene expression is a direct consequence of the network architecture (Box 1). When a determinant represses one side of a mutual repression while activating the other, it is termed ‘coherent’ (i.e. the two actions reinforce each other) and its effect is to irreversibly express one set of genes at the expense of another (Alon, 2007). Because the key feature that allows IRF‐4 to coordinate the two competing states of gene expression in a temporal manner is that it simultaneously but asymmetrically activates both sides of a bistable mutual repression (i.e. the two actions antagonize each other), we term it ‘incoherent’ (see Goentoro et al (2009) for an analogous use in reference to a different regulatory motif). Here, we showed how an incoherent architecture can give rise to both basic bistability and kinetic‐control behaviors depending on the parameters, and we identified means to distinguish between these competing scenarios. A related motif was proposed to enable an activin morphogen gradient to induce alternate states of gene expression for spatial patterning of the Xenopus embryo in a manner analogous to the basic bistability scenario discussed (Saka and Smith, 2007). An incoherent structure that gives rise to dynamics more akin to the ghost attraction is the GRN governing the transient acquisition of competence for DNA uptake by Bacillus subtilis in response to nutrient limitation (Suel et al, 2006, 2007). This response is controlled by a factor that positively regulates expression of its own gene and represses a gene encoding an inhibitor of its degradation. The incoherent structure of the network enables transient expression of the genes associated with competence before returning the cells to the original vegetative state. Recently, an incoherent feed forward loop was also suggested to give rise to a transient pulse of expression of a glucose transporter in yeast (Kuttykrishnan et al, 2010). To the best of our knowledge, our study is the first to provide experimental evidence for a network structure that enables a cell to sense an environmental or developmental signal and transiently realize a pattern of gene expression that promotes cellular diversification prior to realization of a terminally differentiated state.
In the context of antigen‐specific B cells, transient AID catalyzed CSR/SHM gives rise to a diverse repertoire of plasma cells. We speculate that a similar developmental strategy could enable cellular diversification for biological advantage in other contexts. Indeed, the competence GRN discussed above exhibits substantial variation in the period allowing DNA uptake (Suel et al, 2007). Engineered counterparts of this circuit show less variation, suggesting that the natural one is optimized for enhancing cellular diversification resulting from DNA uptake. The discovery that diversity generating retroelements are widespread in bacterial genomes (Medhekar and Miller, 2007) suggests that cellular diversification could be a general feature of adaptive responses, and it will be interesting to determine whether the extent of mutation in these systems is also regulated by transient patterns of gene expression that are of variable duration.
To the best of our knowledge, the mechanism of kinetic control is novel, and it could have important implications for a fundamental process in developmental biology, namely the sensing of a morphogen gradient. In this process, a spatially graded inductive signal is used to instruct a uniform field of progenitors to adopt distinct cell fates based on concentration of the morphogen at a given position along the gradient. The specification of two distinct cell fates based on low versus high concentrations of the morphogen can be readily understood in molecular terms based on differential affinity or accessibility of target genes that are activated by a morphogen‐induced transcription factor. However, this mechanism would require an ever increasing number of distinct transcription factor occupancy states to translate the continuously varying concentration of the morphogen into a large set of cell fates. In contrast, the GRN architecture described here enables the conversion of the spatial information manifest in the morphogen gradient into transient developmental states of continuously varying duration. If such developmental states are in turn based on successive layers of time‐dependent gene expression programs, then one can generate a large set of cell fates simply by varying the timing of exit from the transient inductive state and overt differentiation.
Materials and methods
Core computational model
Because the qualitative but not quantitative features of the regulatory network are known, we use a common generic form for the rate laws (Supplementary information), as we did to model development of myeloid lineages previously (Laslo et al, 2006). Denoting the concentrations of IRF‐4, Blimp‐1, Pax5, Bcl‐6/Bach2 and AID by I, B, P, C and A, respectively, Figure 1B translates to the equations
where kx (for x = I, B, P, C and A) is the maximum rate of protein production stimulated by an activator, Ka and Kr are the half saturation constants for activation and repression, b is the rate of protein degradation, nx and nr are Hill coefficients for activation and repression and ex is the rate of synthesizing x in the absence of activator. The ex terms are included to avoid extending the network to events that are upstream of the species considered. In Equation (1), ex=0 for all species other than Pax5 and IRF‐4 (see further discussion below). As noted in the Supplementary information, a slightly modified form for the rate of AID expression is required to incorporate synergistic activation. The variable s is a scale factor that allows for changing the Michaelis constants for the CSR regulators, while keeping their overall level of expression relatively constant (see below). We add to the equation for dx/dt a stochastic contribution that scales as x1/2 as is typically expected for intrinsic noise.
Both the basic bistability and kinetic‐control scenarios described above can be simulated using this set of equations with different parameters. The kinetic‐control scenario requires only that the initial conditions representing an activated B cell are closer to the CSR state than the plasma‐cell state, allowing each cell to initially adopt a CSR pattern of expression. Because the activated B cell and CSR states share molecular components (Pax5 and Bcl‐6/Bach2 are expressed in both), this is true for the majority of parameter combinations, and we choose generic ones here. On the other hand, for IRF‐4 to control the initial partitioning between cell fates, as in the basic bistability scenario, the regulation of the CSR genes must be strongly cooperative, and the initial conditions must not be strongly biased toward either state. To obtain these features, we set the Michaelis constants and Hill coefficients for the activation of the CSR genes AID and Bcl‐6/Bach2 sufficiently high that they are relatively insensitive to the amount of Pax5 expressed in activated B cells. Specifically, we increase s, nC and nA in Equation (1d) and (1e). See Table II for parameter values employed.
To understand the consequences of the dynamics discussed above for the populations of switched and antibody‐secreting cells, we performed multiscale stochastic simulations. In these simulations, we begin with a population of 10 cells, each of which is in a state corresponding to a naive B cell. The evolution of molecular populations in each cell was generated using the model described above. In addition to these dynamics, we also simulated cellular‐level events: division, death, class switching and differentiation into a plasma cell. These events and the parameters associated with them are detailed in Supplementary Table S1. Each such event was treated as a first order process and was coupled to the fluctuating molecular populations as follows: (1) A cell was assigned a probability per cell division to undergo CSR that is proportional to its AID concentration. (2) A cell with a Blimp‐1 concentration greater than a threshold concentration was assumed to have differentiated into a plasma cell. Additional details concerning the specific means by which we modeled signaling as well as sensitivity analyses are given in the Supplementary information; Supplementary Figures S16–18.
The use and genotyping of the Irf4−/− mice has been previously described (Sciammas et al, 2006). The B1‐8i knock‐in mice were a generous gift of K Rajewsky (Sonoda et al, 1997). The Blimp‐1+/GFP mice were provided by S Nutt (WEHI) (Kallies et al, 2004). Aicda−/− mice were provided by H Huang (University of Chicago) (Muramatsu et al, 2000). Mice were housed in specific pathogen‐free conditions and were used and maintained in accordance of the Institutional Animal Care and Use Committee guidelines.
The Irf4‐inducible mice were generated based on the strategy described by Beard et al (2006). Briefly, the Irf4 cDNA was blunt‐end cloned into the EcoRI site present in the pBS31 flp‐in vector. Next, this vector was cotransfected with an FLPe expression vector into KH2 ES cells. KH2 ES cells had been previously engineered to incorporate frt sites and a defective hygromycin resistance gene in the CollagenA1 locus that would be complemented by successful targeting. In all, 9 of 10 hygromycin resistant clones were properly targeted as judged by Southern blotting using a 3′‐probe specific to the ColA1 locus. Clone O2 was chosen for further analysis and used to generate chimeric mice by blastocyst injection. Chimeric mice were then backcrossed to the C57Bl/6 background and to the Irf4‐null allele.
Cells and culture conditions
Splenic B cells were isolated by sequential hypotonic lysis and magnetic bead separation. Briefly, a single‐cell suspension of spleen cells was treated with Tris‐buffered ammonium chloride to lyse red blood cells and then a negative selection, bead‐based immunopurification strategy was employed that depletes T, NK and macrophage cells using a cocktail of biotinylated antibodies specific for CD4, CD8, CD3, CD49b, Gr1, CD11b, Ter119 (BD Biosciences) and adsorption to Streptavidin microbeads (Miltenyi Biotec Inc.). The unbound cells (typically >90% B220 positive) were washed and placed into culture at cell densities ranging between 0.25 and 0.5 × 106 cells/ml, depending on the experiment. For CFSE labeling, the isolated B cells were washed with PBS and loaded with 1.5 μM CFSE (Molecular Probes) prior to culturing. Cells were stimulated with recombinant human IL‐2 (Roche) 100 U/ml, recombinant murine IL‐4 (R&D Systems) 5 ng/ml, recombinant mouse IL‐5 (R&D Systems) 1.5 ng/ml, as indicated. CD40L membranes were prepared from insect cells and used at 0.1 μg/ml. NP(12)‐ and NP(109)‐Ficoll were obtained from BioSearch Technologies Inc. used at concentrations ranging from 0.01 to 0.00001 ng/ml.
Cells were washed in staining buffer, PBS containing BSA (0.5%, w/v) and sodium azide (0.05%, w/v) and blocked with anti‐FcR (clone 2.4G2) hybridoma supernatant. Cells were stained with biotin‐ or fluorochrome‐coupled antibodies in staining buffer using standard procedures. Anti‐mouse antibodies specific for Syndecan‐1, IgG1 and isotype controls were purchased from BD Pharmingen. Dead cells were excluded from the analysis using DAPI (Molecular Probes). Flow cytometric data were processed using FlowJo software (Tree Star, Inc.). Intracellular detection of IRF‐4 was performed fixing cells with 1% (w/v) paraformaldehyde for 10 min and then permeabilized in staining buffer containing 0.03% (w/v) saponin (Sigma). Staining buffer containing 0.3% Saponin was used to incubate the cells with a polyclonal goat anti‐IRF‐4 or normal goat IgG (Santa Cruz Biotechnology) in the presence of 5% (v/v) donkey serum. Labeling was detected using a Cy5‐coupled donkey anti‐goat (Jackson Immunoresearch).
Mutational analysis of Sμ switch regions
Cells prepared as described above were sorted on day 4 of differentiation based on IgM, IgG1 and Syndecan‐1 into two subsets (i) IgM+ plasmablasts: IgM+, IgG1−, Syndecan‐1+ and (ii) IgG1‐switched cells: IgM−, IgG1+, Syndecan‐1−. Genomic DNA was purified and used as a template for PCR using the M2 primers and Pfu turbo polymerase as described (Xue et al, 2006). The PCR product was purified using gel electrophoresis cloned into the PCR4 Blunt‐TOPO vector (Invitrogen) for DNA sequencing.
We are grateful to members of the Singh and Dinner laboratories for meaningful discussions and to H Shen and U Storb for advice on analyzing mutation frequency. We thank K Rajewsky (HMS), S Nutt (WEHI) and H Huang (University of Chicago) for the B1‐8i heavy chain knock in, the Blimp‐1+/GFP and Aicda−/− mice, respectively. C Beard and R Jaenisch (Whitehead Institute, MIT) provided plasmids, cells and important advice for the generation of the Irf4‐inducible mouse. The University of Chicago Transgenic/ES cell technology mouse core performed the blastocyst injections and enabled the generation of the Irf4‐inducible mouse. Immunology Applications Core and DNA sequencing facilities at The University of Chicago provided flow cytometry and DNA sequencing services, respectively. RS acknowledges support from the Leukemia Research Foundation. AW, YL and ARD acknowledge support from the National Science Foundation, the Chicago Biomedical Consortium (CBC) and the Chicago NIH Systems Biology Center (P50 GM081892). HS acknowledges support from the CBC, NIH (P50 GM081892) and HHMI.
Author contributions: RS, YL, AW, ARD and HS designed the research; RS and YS performed the biological experimentation; YL and AW performed the computation. RS, YL, AW, ARD and HS wrote the manuscript.
Conflict of Interest
The authors declare that they have no conflict of interest.
Supplementary Information, Supplementary Figures S1–18, Supplementary Tables S1–3
Code with Comments
This is an open‐access article distributed under the terms of the Creative Commons Attribution License, which permits distribution, and reproduction in any medium, provided the original author and source are credited. This license does not permit commercial exploitation without specific permission.
- Copyright © 2011 EMBO and Macmillan Publishers Limited