Stochastic Modeling of B Lymphocyte Terminal Differentiation and Its Suppression by Dioxin
© Zhang et al. 2010
Received: 15 October 2009
Accepted: 1 April 2010
Published: 1 April 2010
Skip to main content
© Zhang et al. 2010
Received: 15 October 2009
Accepted: 1 April 2010
Published: 1 April 2010
Upon antigen encounter, naïve B lymphocytes differentiate into antibody-secreting plasma cells. This humoral immune response is suppressed by the environmental contaminant 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD) and other dioxin-like compounds, which belong to the family of aryl hydrocarbon receptor (AhR) agonists.
To achieve a better understanding of the immunotoxicity of AhR agonists and their associated health risks, we have used computer simulations to study the behavior of the gene regulatory network underlying B cell terminal differentiation. The core of this network consists of two coupled double-negative feedback loops involving transcriptional repressors Bcl-6, Blimp-1, and Pax5. Bifurcation analysis indicates that the feedback network can constitute a bistable system with two mutually exclusive transcriptional profiles corresponding to naïve B cells and plasma cells. Although individual B cells switch to the plasma cell state in an all-or-none fashion when stimulated by the polyclonal activator lipopolysaccharide (LPS), stochastic fluctuations in gene expression make the switching event probabilistic, leading to heterogeneous differentiation response among individual B cells. Moreover, stochastic gene expression renders the dose-response behavior of a population of B cells substantially graded, a result that is consistent with experimental observations. The steepness of the dose response curve for the number of plasma cells formed vs. LPS dose, as evaluated by the apparent Hill coefficient, is found to be inversely correlated to the noise level in Blimp-1 gene expression. Simulations illustrate how, through AhR-mediated repression of the AP-1 protein, TCDD reduces the probability of LPS-stimulated B cell differentiation. Interestingly, stochastic simulations predict that TCDD may destabilize the plasma cell state, possibly leading to a reversal to the B cell phenotype.
Our results suggest that stochasticity in gene expression, which renders a graded response at the cell population level, may have been exploited by the immune system to launch humoral immune response of a magnitude appropriately tuned to the antigen dose. In addition to suppressing the initiation of the humoral immune response, dioxin-like compounds may also disrupt the maintenance of the acquired immunity.
In response to antigen stimulation, naïve B cells residing in lymphoid organs such as the spleen and lymph nodes differentiate terminally into antibody-secreting plasma cells [1, 2]. This adaptive humoral immune response can be adversely affected by exposure to some environmental chemicals [3–5]. The environmental contaminant 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD) and other dioxin-like compounds with similar structures suppress humoral immunity, mainly by interfering with B cell differentiation and subsequent antibody secretion [5–8]. These compounds, capable of producing a variety of additional toxic responses including cancer, liver damage, and developmental defects, pose a serious potential risk to human health .
Consistent with the idea of two mutually exclusive steady states for a bistable system, naïve B cells are characterized by a transcriptional profile of low Blimp-1 and high Bcl-6/Pax5 expression, whereas plasma cells feature the opposite profile: high Blimp-1 and low Bcl-6/Pax5 [1, 27]. Pax5, acting as a transcriptional repressor, is essential for maintaining the B cell identity . It actively suppresses transcription of the immunoglobulin molecular components including the heavy chain (IgH), κ light chain (Igκ), and J chain, as well as X box protein-1 (XBP-1) which promotes the formation of the cellular apparatus for immunoglobulin assembly and secretion . The transcriptional repression of these genes by Pax5 ensures that naïve B cells neither produce nor secrete antibodies. Upon antigen encounter, the transcriptional profile in naïve B cells is reversed. High expression of Blimp-1, the 'master regulator' of B cell terminal differentiation, activates a number of downstream genes responsible for the plasma cell phenotype . Most importantly, a simultaneous decrease in Pax5 gene expression releases IgH, Igκ, and J chain from repression, leading to enhanced production and secretion of immunoglobulin by plasma cells [28, 29]. Depending on the type of antigen or B cell stimulus and whether T cells are involved, distinct signal transduction pathways may be activated to drive B cells to differentiate into plasma cells [30–33]. In the case of the bacterial endotoxin lipopolysaccharide (LPS), a polyclonal B cell activator, the differentiation is driven, at least in part, via AP-1, a dimeric transcription factor that upregulates Blimp-1 gene transcription [18, 34].
The suppressive effect of TCDD and dioxin-like compounds on the humoral immune response is currently believed to occur via the following mode of action. By binding and activating the aryl hydrocarbon receptor (AhR), TCDD inhibits AP-1 gene transcription indirectly, leading to reduced AP-1 binding to the Blimp-1 promoter [35, 36]. Since AP-1 mediates the action of LPS, TCDD would thus impair the capability of LPS to drive the switching of the transcriptional gene circuit that underlies B-to-plasma cell differentiation, thereby suppressing the formation of antibody-secreting plasma cells . In addition, TCDD-activated AhR binds to the dioxin response element in the 3' α enhancer of the IgH gene, in a limited fashion inhibiting transcription of this immunoglobulin component .
Gene expression is inherently stochastic . This is particularly true for low-abundance proteins, such as transcription factors, for which the steady-state expression level may fluctuate significantly [39–44]. In many circumstances, this intrinsic noise in protein expression is undesirable for cellular homeostasis and as such needs to be minimized. On the other hand there are instances where noise in protein expression is exploited by cells to generate advantageous non-genetic variability in cellular phenotype [12, 45–47]. When an antigen or polyclonal B cell activator is encountered, the all-or-none differentiation of individual B cells occurs in a seemingly random fashion with new plasma cells appearing at various times [22, 32, 34]. This heterogeneity in the differentiation response is likely to be a key factor contributing to the dose response pattern for a population of B cells. It is possible that the intrinsically stochastic gene expression of Blimp-1, Bcl-6, and Pax5 is responsible for the observed heterogeneity in the differentiation response among otherwise identical B cells. In our previous work, we have suggested the emergence of bistability with a deterministic representation of the transcriptional network underlying B cell differentiation . In the current paper, we use a stochastic computational model of this transcriptional network (Figure 1), constrained by experimental data, to examine how stochastic fluctuations in the abundance of the key transcriptional repressors can make B cell differentiation a dose-dependent probabilistic event. The present study indicates that the stochastic effect plays a central role in modulating the shape, especially the steepness, of the dose response curves for both LPS and TCDD, a result with significant implications for quantitative risk assessment of dioxin immunotoxicity.
Virus-free, female B6C3F1 mice (6 weeks of age) were purchased from Charles River (Portage, MI). Mice were randomized, transferred to plastic cages containing bedding (five per cage), and quarantined for 1 week. Mice were given food (Purina certified laboratory chow) and water ad libitum and were not used until their body weight was 17-20 g. Mice were used in accordance with the guidelines set forth by the Michigan State University Institutional Animal Care and Use Committee (East Lansing, MI).
Spleens were isolated and made into single-cell suspensions aseptically. Primary B cells were isolated from the splenocytes by depleting all other cell types by magnetic separation using the B cell isolation kit from Miltenyi Biotec Inc. (Auburn, CA). B cell purity was always greater than 95% as assessed by flow cytometry. B cells were activated with different concentrations of LPS (Sigma-Aldrich Co., St. Louis, MO) at 3 × 106 cells/ml in RPMI 1640 supplemented with 10% BCS, 100 U/ml penicillin, 100 μg/ml streptomycin, and 50 μM 2-ME. Vehicle (0.02% DMSO) and appropriate TCDD concentrations (AccuStandard, New Haven, CT) were added just prior to LPS addition in a 48-well tissue culture plate.
ELISpot assay was performed as previously described with a few modifications . ELISpot filter plates (MultiScreen-HA) (Millipore, Billerica, MA) were precoated with goat anti-mouse IgM antibody (Sigma-Aldrich Co., St. Louis, MO) at 5 μg/ml in a sodium bicarbonate buffer (pH 9.6), washed and blocked with PBS containing 5% BSA. Purified B cells were washed twice and diluted in RPMI containing 10% BCS and loaded into the culture plate wells. Typically 1-4 × 103 cells in 100 μl cell culture media were incubated on precoated plates overnight at 37°C with 5% CO2. After incubation, plates were washed and biotin-conjugated goat anti-mouse IgM antibody (Sigma-Aldrich Co., St. Louis, MO) was added to the plates at 1 μg/ml and incubated for 2 h at room temperature. To develop spots, streptavidin-HRP was added (Sigma-Aldrich Co., St. Louis, MO) followed by the aminoethylcarbazole per staining kit instructions (Sigma-Aldrich Co., St. Louis, MO). Data was acquired and analyzed using the CTL ImmunoSpot system (Cellular Technology Ltd., Shaker Heights, OH).
IgM secretion was measured using an enzyme-linked immunesorbent assay protocol as previously described . Briefly, supernatants were collected from experimental cultures from which LPS-induced IgM antibody responses were evaluated. Anti-mouse IgM capture antibody (Sigma-Aldrich Co., St. Louis, MO) was added to wells of a 96-well microtiter plate (100 μl/well at 5 μg/ml) and incubated at 4°C overnight. The plate was washed three times with 0.05% Tween-20 in PBS and three times with dH2O. Blocking buffer (3% BSA in PBS) was added to each well and incubated at room temperature for at least 1.5 h. This was followed by the same washing steps described above. Standards (mouse IgM) or supernatant samples were then added to the blocked plate and then incubated at 37°C for 1.5 h. After the incubation, the plate was washed again, followed by addition of 100 μl HRP-conjugated goat anti-mouse IgM detection antibody. After 1.5 h at 37°C, the plate was again washed and 100 μl ABTS (1 mg/ml; Roche, Indianapolis, IN) was added. The detection of the HRP substrate reaction was conducted over a 1 h period using a Synergy HT automated Microplate reader with a 405-nm filter (Bio-Tek, Winooski, VT). The KC4 computer analysis program (Bio-Tek, Winooski, VT) calculated the concentration of IgM in each well based on a standard curve generated from the absorbance readings of known IgM concentrations.
The model structure is illustrated in Additional File 1: Figure S1. Simulation of gene expression control of the transcriptional repressors Bcl-6, Blimp-1, and Pax5 was based on the current understanding of eukaryotic gene regulation [43, 44, 50, 51]. At any given time, a gene could be in one of two discrete transcriptional states: inactive or active, corresponding to the compact and relaxed chromatin structure of the promoter, respectively. Once in the active state (GENE1, Additional File 1: Figure S1), the gene is transcribed at a relatively constant rate; in the inactive state (GENE0), no transcription occurs. Transitions between the inactive and active states (i.e., gene activation and deactivation) are controlled by transcriptional activators and repressors specifically targeting the promoter. In addition to promoter transition, subsequent steps including transcription, translation, and mRNA and protein degradation were modeled for each of the three transcriptional repressors. Since these key steps contribute, to varying degrees, to the stochastic fluctuation of gene expression [38–44], their explicit inclusion adds variability to the switching behavior of the bistable circuit. In turn, the stochastic behavior of the bistable switch likely forms the basis for the heterogeneous differentiation response observed for a population of B cells stimulated with an antigen or B cell activator [22, 32, 34].
The negative mutual regulation among Bcl-6, Blimp-1, and Pax5 is modeled based on the modes of repression reported in the literature. Blimp-1 represses gene expression of both Bcl-6 and Pax5 [15–17]. As a transcriptional repressor, Blimp-1 suppresses target genes by recruiting co-repressors such as histone deacetylase, histone methyltranferase, and those belonging to the Groucho family [52–54]. Functioning as chromatin-modifying enzymes, these co-repressors alter the local chromatin structure to a compact (transcriptionally inactive) state. Based on this active mode of repression, the transcriptional repression of Bcl-6 and Pax5 by Blimp-1 was implemented by having Blimp-1 promote the deactivation step of these two genes. Pax5 in turn represses Blimp-1 gene expression directly , also by an active mechanism mediated through recruitment of co-repressors from the Groucho family [55, 56]. Therefore, repression of Blimp-1 by Pax5 was implemented by having Pax5 promote the deactivation step of Blimp-1 gene. Although the transcriptional repressor Bcl-6 may directly bind to target genes, it exerts its effect on Blimp-1 indirectly through AP-1. By binding and sequestering AP-1, Bcl-6 blocks its transcriptional activity, thus functioning as a passive repressor . Since AP-1 positively regulates Blimp-1 gene expression, repression of Blimp-1 by Bcl-6 was implemented by having Bcl-6 impinge upon the activation step of Blimp-1 gene, thus curbing the maximal induction of Blimp-1 by AP-1. Through these specific transcriptional regulations, a gene circuit is established in the form of coupled double-negative feedback loops among Bcl-6, Blimp-1, and Pax5 (Additional File 1: Figure S1).
Activation of the transcriptional program underlying B cell terminal differentiation occurs in a stimulus-dependent manner, involving various cytokines and cell surface receptors . In the case of LPS, binding to its cognate Toll-like receptor 4 (TLR4) initiates a sequence of intermediate intracellular signaling events, culminating in the activation of the AP-1 protein complex [34, 35]. For simplicity, the activation of AP-1 by LPS-bound TLR4 was modeled as a single phosphorylation event. Like many other membrane-residing receptors, TLR4 is usually downregulated after occupied by LPS, possibly through a desensitization mechanism that includes receptor internalization [57, 58]. Downregulation of TLR4 was implemented here by ascribing a higher rate constant to the turnover of the LPS-bound receptor compared to that of the free receptor. Inclusion of TLR4 downregulation enables the model to recapitulate the attenuation of AP-1 signaling in time following its initial activation, as observed in B cells stimulated with LPS and other cytokines [34, 35].
The suppression of AP-1 activity by TCDD in LPS-activated B cells was assumed to occur in an AhR-dependent manner. Although the precise mechanism for this process is unknown, a decrease in the AP-1 subunit c-jun was observed in the AP-1 DNA-binding complex in the presence of TCDD . For simplicity, AP-1 was treated as a single entity in our model with its synthesis rate negatively regulated by TCDD. In the absence of detailed information at the molecular level on the signaling pathways by which TCDD-AhR affects AP-1, the specific regulation of AP-1 by TCDD-AhR was described in the equation as a continuous function of the concentration of TCDD-AhR. Immunoglobulin M (IgM), which consists of several subunits, was also modeled as a single entity subject to negative control by both Pax5 and AhR. A first-order secretion rate was assumed for intracellular IgM, with the rate constant large enough that intracellular degradation of IgM is negligible.
Reaction details including ordinary differential equations, parameter values, and initial steady-state conditions are provided in Additional File 1: Tables S1 - S3. Unless otherwise indicated, the unit of abundance used for state variables in the model is number of molecules per cell. The deterministic version of the model was implemented in PathwayLab (InNetics Inc., Linköping, Sweden), while the stochastic version was implemented in the BioNetS program , based on Gillespie's stochastic simulation algorithm . Both versions of the model were exported into MatLab (The Mathworks, Inc., Natick, MA, USA) for analysis and simulation. The model in SBML format is provided in Additional File 2. Bifurcation diagrams were generated using the XPP-AUT program . The effective Hill coefficients (nH) of simulation-generated dose response curves were estimated by using the equation nH = ln81/ln(X90/X10), where X90 and X10 are the doses that produce 90% and 10% of the maximal response, respectively.
By activating AP-1, LPS triggers the differentiation of B cells to plasma cells. Bifurcation analysis (Figure 2B) indicates that as the steady-state AP-1p (active form of AP-1) level increases (rightward arrows), Blimp-1 increases and Pax5 decreases initially by small amounts. However, once AP-1p reaches a threshold value of about 200 molecules, the B cell switches in a discrete fashion to a plasma cell state (vertical arrows). This abrupt transition underlies an all-or-none, switch-like response and signifies a true discontinuity that precludes the cell from settling at an intermediate state. The switch also exhibits hysteresis. Starting from the plasma cell state, as AP-1p is reduced (leftward arrows) to its basal level (~100), the system remains in the differentiated plasma cell state with high Blimp-1 and low Pax5/Bcl-6 levels. In the context of humoral immune response, this irreversibility of the switch ensures that the antibody-secreting plasma cell phenotype persists after the initial antigen stimulus recedes, thus serving as a maintenance mechanism for the acquired immunity.
In a deterministic model of a bistable system, the occurrence of switching from one state to the other and the time when it occurs are uniquely determined by the level of the input stimulus . If the biochemical circuits are nearly identical across a population of B cells, the individual responses of these cells would be expected to be similar. However, experimental observations indicate that the responses of individual cells to a given activator are rather heterogeneous [22, 32, 34], even in an isogenic B cell population . Some cells differentiate into antibody-secreting plasma cells at early times, some do so at late times, and many cells remain undifferentiated in the time window of observation. Thus the number of plasma cells that appear after B cell activation is time-dependent. This heterogeneous response pattern is likely to be important in shaping dose response behaviors, which are quantitative measurements necessary for evaluating the health risk of immune-suppressive chemicals. Therefore a more explicit description of the heterogeneous response of individual B cells is needed than is provided by a deterministic model.
Effects of noise in Blimp-1 protein expression on the steepness and position of the LPS dose response curves.
% plasma Cells
Many cellular processes in living organisms, including sex determination, lineage specification, cell differentiation, proliferation, and apoptosis, are irreversible, all-or-none (binary) phenomena. It is increasingly recognized that such binary cellular processes are controlled by bistable biochemical circuits, which ensure the discreteness and irreversibility of these processes [62, 70]. Most cellular bistable systems are based on molecular circuits of mutual activation and/or mutual inhibition between genes and proteins . These molecular circuits often involve a number of proteins forming multiple feedback loops, as observed in systems for cell cycle control, oocyte maturation in Xenopus laevis, cell lineage specification and differentiation, etc [11, 72–75]. The coupling of multiple feedback loops is believed to provide redundancy as well as robustness to the bistable property of the entire circuit [73, 76]. In the B cell transcriptional network, double negative feedback loops between Bcl-6, Blimp-1, and Pax5 are coupled at the Blimp-1 node. The resulting bistability produces two mutually exclusive transcriptional profiles: (1) high Bcl-6/Pax5 and low Blimp-1, representing the B cell state; and (2) low Bcl-6/Pax5 and high Blimp-1, representing the plasma cell state (Figure 2A). This mutual exclusivity of stable gene expression patterns ensures definitive separation between two discrete cellular phenotypes and an all-or-none response in individual B cells to antigen stimulation. Further, once a B cell differentiates into a plasma cell, the plasma cell phenotype is retained even if the stimulating antigen is removed. This irreversibility is underpinned by hysteresis, another important biological property conferred by the bistable system (Figure 2B).
The switching from the B cell to plasma cell state is stimulated by specific antigens or polyclonal activators. In the case of LPS, the switching is mediated by alterations in concentration of active AP-1 protein (AP-1p) and characterized by a threshold concentration of AP-1p (Figure 2B). This in turn suggests a threshold for LPS, the initiating stimulus, and a switch-like dose response in a deterministic setting. In reality, the bistable gene circuit in the B cell has to operate in an inherently stochastic environment within the nucleus, which is likely to produce significant cell-to-cell variability. Our experimental data shows that individual B cells differentiate into plasma cells in a rather heterogeneous manner, with the percentage of plasma cells formed appearing to be a graded function of the LPS dose, devoid of the sharp transitions expected for an all-or-none bistable switch (Figure 6A). By considering the stochasticity in gene expression, our simulations recapitulate these observations. Thus whether a particular B cell will respond or not within a given time window of LPS stimulation, as well as the timing of the response, are both chance events resulting from intrinsically stochastic gene expression.
Due to the low expression level of Blimp-1 in B cells , its abundance is likely to show considerable moment-to-moment and cell-to-cell variability. In our implementation of the B-cell differentiation network, Blimp-1 fluctuates in a nearly pulsatile manner (Figure 3), which is consistent with the phenomenon of pulsatile mRNA production and bursting of protein translation observed for genes with low expression [77–80]. In the phase space of the bistable system, mature B cells would reside on one side of an imaginary boundary separating the two basins of attraction for the B cell and plasma cell states. Even in the absence of LPS, an occasional pulse in Blimp-1 expression may be large enough to send the system across the boundary to the plasma cell state. This might explain the spontaneous appearance of a negligibly small fraction of plasma cells observed in a population of B cells in vitro in the absence of any antigen [22, 23, 34, 81]. Exposure to LPS, which activates the AP-1 protein, leads to increased Blimp-1 gene transcription. Due to the stochastic nature of Blimp-1 gene expression, the LPS-induced increase in Blimp-1 abundance takes the form of more frequent and larger Blimp-1 pulses (Additional File 1: Figure S7). These pulses increase the probability of the system crossing the fate-separating boundary and then being attracted to the plasma cell state eventually. Our analysis suggests that the cell-to-cell variability in Blimp-1 pulses at a given moment explains the heterogeneous switching behavior of individual B cells. Previous computational studies have shown that stochastic fluctuations could obscure the threshold of the switching behavior of a bistable system, making the dose response less switch-like [82, 83]. In our implementation of the B cell differentiation network, sufficiently large stochastic fluctuations in Blimp-1 expression result in substantially graded response for percentage plasma cell formation with an estimated Hill coefficient close to 1. This number, a measure of the steepness of the dose response curve, is inversely correlated with the level of noise in Blimp-1 protein expression (Figure 7 and Table 1).
In summary, the bistable gene circuit operating in a stochastic environment in B cells confers two essential properties to the humoral immune response. First, bistability ensures that distinct, mutually exclusive phenotypes are associated with the B cell and plasma cell, and that the switching from the B cell to plasma cell phenotype is irreversible so as to maintain the acquired immunity for a period of time following antigen encounter. Second, stochastic fluctuations in protein expression provide the necessary variability for the differentiation response to becomes probabilistic. This allows the number of plasma cells formed and the total amount of antibody produced to scale in a graded manner with the antigen dose. Numerous studies have shown that noise in protein expression can be exploited by cells to generate necessary non-genetic variability in cellular phenotype and fate [12, 45–47, 84–86]. For example, intrinsic noise in ComK protein expression in the soil bacterium Bacillus subtilis allows a small fraction of cells to reside in the competent state capable of DNA uptake rather than in a vegetative state. This strategy increases fitness of the species by increasing phenotypic diversity among a population of genetically identical bacteria [12, 87]. A recent study by Spencer et al. demonstrated that protein expression noise may be responsible for the variability in apoptotic response observed in clonal populations of mammalian cells . It is thus likely that noise in gene expression is utilized by B cells to launch a humoral immune response of appropriate magnitude (as measured by the extent of plasma cell formation and aggregated IgM secretion) in handling the amount of pathogen in the body. This proposed effect of variability at the molecular level on the overall dose response curves is similar to the hypothesis on chemical carcinogenesis at the human population level proposed by Lutz, which states that even if each human individual had a specific response threshold for a carcinogen, variability among individuals would result in a continuous dose response curve for a human population .
The molecular basis of the immunosuppressive effect of the environmental contaminant TCDD and similar compounds acting through AhR remains incompletely understood at the molecular level. A primary goal in developing cell-based quantitative network simulations of toxic response is to evaluate the shape of the dose response curve over a broad range of exposures. To this end, it is necessary to acquire mechanistic understanding of how dioxin-like compounds interact with the gene transcriptional network underlying the B cell terminal differentiation program. As we have tried to demonstrate in the current work, computational modeling approaches can be helpful in attaining this goal. Through activation of the AhR signaling pathway, TCDD, a potent dioxin compound, suppresses AP-1 signaling, a key mediator of B cell differentiation [35, 36]. This attenuation of antigen-stimulated AP-1 activation by TCDD (Additional File 1: Figure S6) makes it more difficult for a B cell to switch to the plasma cell state. In a stochastic gene expression context, this inhibitory effect of TCDD is manifested as a reduction in the probability of B cell differentiation. With fewer antibody-secreting plasma cells formed, the humoral immunity is compromised. Importantly, stochasticity in gene expression of transcription factors transforms the binary, all-or-none response occurring in individual B cells into a much more graded response to TCDD for a population of B cells (Figure 8). Our computational model, which incorporates as default a continuous description of TCDD affecting AP-1 and simulates in vitro experimental scenarios, suggests that the graded nature of the suppressive effect on B cell differentiation may extend well into the low-dose region, where no abrupt transitions were observed (Additional File 1: Figure S5).
Differentiation of mature B cells into plasma cells is a terminal, physiologically irreversible process. Some of the newly differentiated plasma cells migrate to the bone marrow, where, aided by survival signals from stromal cells, they can survive for several months [89–91]. Independent of memory B cells, these long-lived plasma cells contribute to acquired immunity against pathogens by continuing to secrete antibody molecules for an extended period of time. While the effect on B cell differentiation has been a focus of research concerning immunotoxicity of TCDD and related compounds, its possible effect on the fate of terminally differentiated plasma cells has not received much attention. In its present form, our B cell response network suggests that TCDD may alter the phase space landscape of the transcription program, bringing the plasma cell state closer to the boundary separating it from the B cell state (Figure 11A). This would make it easier for stochastic fluctuations to switch plasma cells back to the B cell state. This behavior occurs in the stochastic implementation but not in the deterministic mode with an identical parameter set. Stochastic simulations suggest that TCDD destabilizes plasma cells in a dose-dependent manner, possibly allowing these cells to dedifferentiate back to a B cell phenotype (Figure 11B). While the predicted decline in the number of plasma cells due to phenotype reversal is a slow process, the impairment to the acquired humoral immunity could be tangible, given that long-lived plasma cells may survive in the bone marrow for up to several months [89–91].
Although cell differentiation is traditionally regarded as a physiologically irreversible process, gene manipulation in terminally differentiated cells has been shown to successfully reprogram these cells to progenitor or even stem cell states [92–95]. Among immune cells, the B cell lineage itself has considerable plasticity [96, 97]. B cells can be dedifferentiated by genetic manipulation to precursor cells, and thereafter diverted through intermediate states to fully differentiated T cells or macrophages [98, 99]. Transiently expressing ectopic Bcl-6 and associated co-repressor MTA3 in plasma cell lines reversed the phenotype to that of B cells . This finding is consistent with our model of the underlying bistable gene circuit, where TCDD first represses AP-1 and Blimp-1, leading to upregulation of Bcl-6 and Pax5, and consequent destabilization of the plasma cell state. While this prediction requires further study, it could have important implications for the immunotoxicity of TCDD. If confirmed, it would imply that besides suppressing initiation of the humoral immune response by inhibiting plasma cell formation, TCDD and similar compounds may disrupt the maintenance of the humoral response through destabilization of long-lived plasma cells.
Although the stochastic model presented here reproduces some experimental observations, the entire transcriptional network of B cell differentiation is potentially far more complex. Besides Bcl-6, Blimp-1, and Pax5, additional transcription factors such as Bach2, IRF-4, Mad1, and MTA3 are also involved in regulation of B-to-plasma cell differentiation [26, 33, 101–103]. Additional feedback loops may also exist, for example between Blimp-1 and IRF-4 . The molecular pathways mediating the toxicity of TCDD are also likely to be more extensive than the impairment of AP-1 signaling. For instance, the role of the AhR repressor protein in TCDD signaling may need to be considered . An important purpose of risk assessment for TCDD and other dioxin-like compounds is to establish specific "safe" values of exposure in the low-dose region, below which adverse biological effects are negligibly small. However, low-dose effects are generally more difficult to measure precisely and economically, as the variability associated with measured endpoints is usually large, obscuring the change in curvature of the low-dose region. Estimation of low-dose effects is therefore frequently made by extrapolation from high-dose data or computer simulations. In the present study, we show that a mechanistic, systems-level computational model of the B cell transcriptional network incorporating stochastic gene expression could generate dose response curves for a broad range of doses in agreement with experimental measurements. However, its value in accurately predicting low-dose effects for both LPS and TCDD is currently limited. A more confident prediction would likely require greater structural and kinetic details of pertinent signaling pathways leading from exposure to LPS and TCDD up to the bistable circuit, and perhaps a more complete description of the bistable circuit itself. Finally, for in vivo toxicity prediction, the cellular model has to be linked to dosimetry models that calculate the concentrations of dioxin-like compounds in immune tissues for environmental exposures.
In conclusion, we demonstrate that a bistable gene regulatory network comprising three genes - Bcl-6, Blimp-1, and Pax5 - exhibits behavior consistent with the terminal differentiation of B lymphocytes to antibody-secreting plasma cells. In a stochastic gene expression environment, the response to both antigen and TCDD at the cell population level appears to be substantially graded even though activation of individual cell is a binary, all-or-none phenomenon.
We would like to thank the Superfund Research Program (P42ES04911) from the National Institute of Environmental Health Sciences for supporting this work.