Transcriptional noise is known to play a crucial role in heterogeneity in bacteria and yeast. Mammalian macrophages are known to exhibit cell-to-cell variation in their responses to pathogens, but the source of this heterogeneity is not known. We have developed a detailed stochastic model of gene expression that takes into account scaling effects due to cell size and genome complexity. We report the results of applying this model to simulating gene expression variability in mammalian macrophages, demonstrating a possible molecular basis for heterogeneity in macrophage signalling responses. We note that the nature of predicted transcriptional noise in macrophages is different from that in yeast and bacteria. Some molecular interactions in yeast and bacteria are thought to have evolved to minimize the effects of the high-frequency noise observed in these species. Transcriptional noise in macrophages results in slow changes to gene expression levels and would not require the type of spike-filtering circuits observed in yeast and bacteria.
White blood cells, such as macrophages, recognize microbial pathogens via the innate immune receptors, of which the Toll-like receptors (TLRs) are critical members (Takeda et al. 2003). Each TLR recognizes specific features of pathogens and differentially activates inflammatory and other immune responses via a regulated series of events that tailor defences to deal with the particular threat. Inflammation can lead to multiple outcomes: resolution of the infection and complete restoration of tissue architecture, tissue destruction (scarring), ongoing (chronic) inflammation, initiation of new inflammatory responses (autoimmunity) and failure to control the infection. Moreover, past and concurrent signalling events can also influence these outcomes, depending on duration and intensity (Aderem et al. 1986; Sato et al. 2000; Gordon 2003).
TLR signalling is well-characterized in terms of inputs (microbial ligands) and outputs (macrophage coordination of immune responses), but not in terms of the signalling and transcriptional regulatory mechanisms that orchestrate the specific ligand-induced responses. These molecular mechanisms presumably account for stimulus- and time-dependent gene regulation via feedback loops and specific circuitry such as switches, multiplexers, amplifiers, oscillators, etc.
2. Specificity of the TLR-mediated macrophage responses to microbes
When a macrophage encounters microbes, TLR signalling alters the expression of hundreds of genes (Nau et al. 2002). The highly specific response of TLRs to different molecular identifiers of pathogens (e.g. lipopolysaccharide, dsRNA, bacterial flagellin, CpG DNA) is well documented. Signalling is induced by the recruitment of proximal adapter molecules to the TLR cytoplasmic domains (figure 1). Several TLR adapters have been identified, and it is hypothesized that selective adapter recruitment underlies the specificity of signal transduction (Horng et al. 2002; Yamamoto et al. 2002, 2003; Hoebe et al. 2003). However, the majority of the approximately 30 known signal transduction proteins downstream of the adapters are used by all TLRs. Thus, the manner in which different cellular responses are generated by a common set of signalling elements remains mysterious. The macrophage's recent history also regulates the TLR pathway. Prior exposure to microbial products, such as lipopolysaccharide, dramatically alters responses, and depending on dose and timing, leads to the phenomena of priming (increased sensitivity subsequent to a small initial stimulation; Aderem et al. 1986) and tolerance (adaptation to long-term stimulation; Sato et al. 2000).
3. Macrophage diversity
An additional complexity in modelling the TLR pathway is that macrophage populations are known to be heterogeneous in terms of cellular state and response (Hume 2000; Ravasi et al. 2002; Hoebe et al. 2003) as illustrated in the example cells in figure 2. Such heterogeneity could be due to a number of factors including: (i) inherited differences in DNA structure or cellular protein profile; (ii) a unique past history that drives individual cells to be in different stable states that are maintained through regulatory networks controlling mRNA and protein levels; and (iii) intrinsic noise in gene expression patterns due to the stochastic nature of transcription and translation, evident particularly at low copy numbers.
It is well established that measurement of cell-averaged responses in heterogeneous populations of cells can be highly misleading (McAdams & Arkin 1997; Vilar et al. 2003). For example, a recent study of the p53 pathway (Lahav et al. 2004), showed that what appear to be damped oscillations in cell population assays are actually different numbers of equal-sized pulses in individual cells. It is, therefore, important to develop single-cell models that provide experimentally testable, mechanistic and quantitative explanations of cellular heterogeneity (Rosenfeld et al. 2005).
4. Stochastic noise in gene expression
In bacteria, seminal studies (McAdams & Arkin 1997; Ozbudak et al. 2002; Swain et al. 2002) have demonstrated that inherent stochastic noise in gene expression (due to small numbers of molecules, thermal noise, etc; referred to as intrinsic variability), as well as variability in cellular transcription factor activity levels (e.g. inherited factor concentrations; referred to as extrinsic variability), can result in cellular heterogeneity. A recent study (Raser & O'Shea 2004) measured comparable levels of gene expression variability in yeast. The study also showed that the kinetics of individual steps in gene expression, such as transcription factor complex formation, RNA polymerase recruitment, and translational efficiency, can vary the amount of intrinsic noise in gene expression several fold.
As Raser and O'Shea point out, intrinsic gene expression noise has been invoked as a source of phenotypic variation in a number of very different settings: the lambda phage lysis–lysogeny switch (Arkin et al. 1998), phase variation in bacteria (Hallet 2001), receptor expression on mammalian olfactory neurons (Serizawa et al. 2003) and tumour formation (Magee et al. 2003). Does it also underlie macrophage heterogeneity?
As we (Orrell & Bolouri 2004; de Atauri et al. 2005) and others (Becskei & Serrano 2000; Morishita & Aihara 2004) have shown, intrinsic noise can be filtered out by some intra-cellular biochemical networks (e.g. negative feedback, dimerization). Extrinsic variability leading to cells being in different states can arise from stochastic and deterministic origins. Stochastic differences in cellular content could arise from non-deterministic processes such as the numbers of molecules, molecular complexes, or organelles inherited during cell division. Deterministic differences could occur through cell–cell interactions, or through processes such as DNA re-arrangement, or chromatin remodelling during cell division. Another extrinsic source of variability in macrophages could be that each macrophage senses a slightly different amount of ligand (for example due to expression of different numbers of receptors), and multiple ultra-sensitive thresholds in the signal transduction pathway act as a decision tree, distributing macrophage responses according to the amounts of ligand-sensed.
Intuitively, one may expect much less intrinsic noise in larger mammalian cells. This intuition is principally based on a scaling argument. Figure 3 compares the volumes of four cell types. Yeast cells are approximately 10 times larger (by volume) than bacteria. Macrophages are approximately 1000 times larger (by volume) than Escherichia coli. Literature reports of mRNA and protein copy numbers in each cell type suggest that these copy numbers scale with the cell volume (Neidhardt 1996; Ghaemmaghami et al. 2003; Lehner & Cresswell 2004). The coefficient of variation of the Poisson distribution decreases with increasing average copy number, so we might expect macrophages to be virtually free of intrinsic variability. In addition, if the half-life of a protein is very long (as it can be in mammalian cells), then gene expression (within a single-cell context) can stop after a short burst of synthesis. From that point onwards, there will be no noise due to stochastic gene expression. Moreover, the number of copies of extremely slowly degrading proteins tends to be very high, so one might argue that the average behaviour is essentially noise-free. Finally, we note that the number of copies of most genes is limited to one in E. coli and haploid yeast, and two in diploid yeast and in macrophages; this is especially true for transcription factor genes (Nelson et al. 2004).
As we will demonstrate, the intuition that macrophage cells must have very low cell-to-cell variability due to the above size-scaling arguments turns out to be incorrect. There are two reasons for this. Firstly, as shown in figure 3, genome size grows even faster than volume as we go from E. coli to yeast to macrophages. This sharply increases the number of non-specific transcription factor–DNA interactions. As a result, we may expect transcription factor occupancy on cis-regulatory DNA to be approximately at the same level in all four organisms considered. Secondly, transcription factors, because they need to be regulated dynamically, tend to have shorter half-lives and fewer copies per cell than most structural genes, and are thus a source of extrinsic noise for the genes they regulate.
5. Model of gene expression at the single-gene scale
We have constructed stochastic models of gene expression in macrophages, yeast and bacteria. Our models differ from previous models (McAdams & Arkin 1997; Hasty et al. 2000; Ozbudak et al. 2002; Swain et al. 2002; Raser & O'Shea 2004) in that, instead of Langevin equations and/or formulations based purely on first order kinetics, we stochastically simulate multiple steps in transcription and translation and perform statistically large numbers of single-cell in silico experiments. Furthermore, we take into account the effect of non-specific binding of transcription factors on the average fractional activation of a gene (Bolouri & Davidson 2003). A summary of the steps in our model of gene expression is given in figure 4. Due to the time delays for transcription and translation, our stochastic model is non-Markovian (Gibson & Bruck 2000). The complete description of our model is presented, in the Dizzy model definition language (Ramsey et al. 2005), in the electronic supplementary material. Our bacteria, yeast and macrophage models differ from one another by more than cell and genome size; as far as possible, we employ species-specific parameters as described in table 1.
To distinguish between intrinsic and extrinsic noise, our model contains two similar reporter genes (ostensibly coding for cyan and yellow fluorescent protein, respectively) with identical promoters and identical fractional activation of gene expression. The rates of initiation of transcription and translation, and the rates of degradation for mRNA and proteins were assumed to be the same for the two reporters. The example promoter used in the studies reported herein has cis-elements for two transcription factors that bind cooperatively. Both factors must be bound for transcription to occur, in the model used in this study. However, our method is general and can capture regulation by multiple transcription factors. We used a Monte Carlo technique to simulate the stochastic dynamics of this system, under steady-state conditions, at different fractional levels of gene activation. With parameters appropriate to bacteria and yeast, our theoretical model of expression noise captures the reported characteristics of transcriptional noise well. Figure 5 shows a side-by-side comparison of our simulation results with experimental observations for E. coli and yeast. These results demonstrate that our model of transcription for yeast and bacteria captures the relationship between intrinsic and extrinsic noise, and the dependency of the intrinsic noise magnitude on the cell-averaged gene expression level.
We adapted our eukaryotic model of transcription and translation to the macrophage system by altering the parameters in the model to values appropriate to a mammalian macrophage. Table 1 lists the parameters which we extracted from the literature. Many of these parameters are known to adequate accuracy for our needs (e.g. the size of the genome, codon-lengths of specific proteins). Others appear to have little effect on our predictions (e.g. the numbers of polymerase and ribosome molecules reported are well in excess of the numbers of transcribing genes and transcripts, respectively). By far the most critical parameters of the models are (i) the rate of initiation of transcription and translation; and (ii) the half-lives of mRNAs and proteins. Half-lives are known to vary considerably across the transcriptome and proteome (Pratt et al. 2002; Wang et al. 2002). However, for the studies reported here the important issue is ratio of half-lives of similar proteins in protozoans and mammals. We note that—due to dilution caused by rapid cellular growth (cellular volume approximately doubling prior to every cell division)—the half-lives of mRNA and proteins in budding yeast and bacteria during exponential growth is necessarily less than the cell division time. Since mammalian cells such as macrophages grow and divide an order of magnitude more slowly, their protein half lives can be correspondingly longer, which leads to qualitatively different transcriptional noise characteristics.
Using our model system of two reporter genes, we investigated whether a macrophage gene (at an equivalent fractional level of activation) has less intrinsic noise than a yeast or bacterial gene. With a fractional occupancy of approximately 10–15% in all three models (bacteria, yeast and macrophage), we found that the coefficient of variation (the standard deviation divided by the mean) for the intrinsic noise contribution to protein abundance was not significantly different. The results for all three model organisms are shown in figure 6; the spread of data transverse to the axis of equal-gene expression represents variation due to intrinsic gene expression noise.
Next, we studied the stochastic dynamics of transcriptional activation. We found that in bacteria and yeast, at low gene activation, protein abundances can exhibit large transient spikes. This is due to the comparatively short lifetime for protein and mRNA in those systems, and the predominance of intrinsic noise in the mRNA concentration at low expression levels. Figure 7b shows an extreme example in the yeast model for a gene at basal (very low) activation. The protein abundance is seen to exhibit significant transient spikes that are highly correlated (with a fixed delay) with the previous production of a completed mRNA transcript. Note the simulation was performed over an artificially long period in order to capture a few examples of such activity spikes.
In the macrophage, we found that the protein abundance showed generally lower-frequency (slowly varying) noise. Figure 8 compares representative simulated protein levels for a macrophage gene with parameters typical of transcription factors. Transcription was activated at time zero. We simulated the dynamics of the model over a period of approximately 48 h, typical of periods over which macrophages are experimentally assayed. Note that the variability in each trace is considerable, and that it occurs on a very slow time-scale. The lack of high-frequency noise, as compared to simulated single-cell protein abundances in prokaryotes, is due to the relatively slow mRNA and protein degradation in mammalian cells. The upper two curves (red and blue) represent the gene at just 5% more transcriptional activity (higher transcription factor occupancy) than the bottom two curves (black and green). Note that, as indicated by the areas highlighted in red and blue, the two slightly more active genes have protein levels more than 10 times in excess of the 5% less activated genes for periods exceeding a day. The long-term average expression level of the less active gene pair is only about 30% lower than that of the more active gene pair. Thus, the intrinsic expression noise for each gene, when added to small, potentially stochastic variations in cellular content (here 5%) can result in very clear heterogeneity in cellular states as defined by gene expression levels.
To better understand the role of transcriptional noise in macrophages, we systematically compared the steady-state noise profile of a single gene for the three model systems (bacteria, yeast and macrophage) at steady-state. A large ensemble of 7500 stochastic simulations was used, to ensure accuracy in computing the standard deviation and average protein abundance. Figures 9 and 10 summarize our findings for the cross-species comparison of single-gene expression. In figure 9, we display the average mRNA level, protein level, coefficient of variation of the protein level (standard deviation divided by the mean) and the Fano factor of the protein level (the ratio of the variance to the mean). While the coefficient of variation in the protein level is slightly less in macrophages than in yeast and bacteria, the magnitude of protein abundance noise, given by the Fano factor is twofold higher in macrophages. A Poisson process will result in a Fano factor of 1; the large Fano factor for macrophages indicates a high degree of variability. For protein abundance, a related measure of noise is the ‘burst size’ of protein production, which is the number of proteins produced over the lifetime of an mRNA transcript (Ozbudak et al. 2002). In our model, the protein burst size is determined parametrically, and is 99.25, 1004.6 and 5437.4, respectively, for bacteria, yeast and macrophage (irrespective of activation level). Considering the data in figure 2e, the twofold higher noise scale in macrophages could lead to considerable long-term heterogeneity, despite the large average abundance of a typical structural macrophage protein.
Using the formula of (Ozbudak et al. 2002), we estimated the intrinsic noise contribution to the variance in the protein abundance as the square root of the ratio of the ‘burst parameter’ to the average protein abundance. We compared this intrinsic noise estimate to the total average variation observed in the simulation, for the three species and conditions of 10, 50 and 90% gene activation. The results are summarized in figure 10 and confirm our observation that intrinsic noise coefficients in yeast, bacteria and macrophage are comparable at 10% activation. They also show that extrinsic noise makes an important contribution at 10% activation, but is insignificant at 90% activation.
6. Implications of transcriptional noise for organizational principles in macrophages
We showed above that macrophages can exhibit significant gene expression noise, and that they have fundamentally different noise characteristics from yeast and bacteria. In this section, we speculate that this difference can lead to different evolutionary pressures and hence different organizational principles in macrophage genetic regulatory networks. We illustrate our argument with a specific example: the feed-forward loop (FFL) network motif, which has been found to be prevalent in bacterial and yeast regulatory networks (Milo et al. 2002; Shen-Orr et al. 2002; Luscombe et al. 2004). In the context of a genetic regulatory network, the FFL can consist of three genes arranged in a three-gene cascade, with the first gene additionally acting as an input to the third gene.
Figure 11a shows a gene network diagram with a FFL. Depending on whether each link in the motif upregulates or downregulates its target, and whether the third gene's inputs act as a logical AND or OR, the FFL can act as either a sign-sensitive delay or a sign-sensitive accelerator, and the output can either be inverted or non-inverted (Mangan & Alon 2003). Sign sensitivity means that the delay effect depends on the sign of the input transition, i.e. whether the input changes from the ‘high’ to the ‘low’ state, or vice versa.
In E. coli and yeast, the most abundant FFL motif has ‘coherent’ upregulation for all three links (Mangan & Alon 2003). In this configuration, the FFL acts as a delay when the input to the first gene goes from a low (inactive) state to a high (active) state. When the third gene's cis-regulatory inputs combine in a logical AND, this configuration is a sign-sensitive delay, without output inversion. It can act as a filter to suppress spikes when the first gene's input is normally in the low state. Simulation results demonstrating this effect, for a simulated gene cascade implanting a FFL in yeast, are shown in figure 11b. The simulation results illustrate the noise suppression that a FFL can exert in the case of high-frequency burst-like noise on top of a low steady state input. We note, however, that unlike a low-pass filter, the FFL as described above would not delay a rapid transition of the input from the high state to the low state; in this sense, the FFL is functionally quite distinct from a feedback-type noise filter. Given that E. coli and yeast tend to have high intrinsic transcriptional noise, it is possible that the prevalence of this particular FFL type may be explained by the need to suppress propagation of spurious noise spikes through a regulatory cascade, so that downstream genes will not be activated unless the upstream controller genes switch to the high state for a minimum duration.
Recent studies have indicated that the FFL motif is also prevalent in mammalian genetic regulatory networks. A recent chromatin–immunoprecipitation study (Odom et al. 2004) mapped out the regulatory network controlled by the HNF family of transcription factors, in human pancreas and liver cells. The authors demonstrated that HNF6 is the controller of a FFL regulating the gene PCK1. Another study (Hinz et al. 2000) reported evidence that PGE2 acts through a FFL to regulate the expression of its own synthesizing enzyme (COX-2) in RAW 264.7 cells. The authors speculate that the purpose of this feed-forward interaction may be to modulate COX-2 expression in macrophages within inflamed tissues. These studies lead us to consider the general implications of a FFL in a genetic regulatory network within a macrophage. The propensity of gene expression within a macrophage to exhibit a slowly varying protein abundance (as shown in figure 7 of this paper), as opposed to the high-frequency burst-type noise observed in prokaryotes, would seem to obviate the need for a FFL implementing a sign-sensitive delay. The delay effect mediated by the middle gene (gene Y in figure 11a) would need to be of the order of several hours to act as a filter for the type of expression noise we see in macrophages. Thus, it seems unlikely that the FFL motif acts a transcriptional noise filter in macrophages. Perhaps the FFL within a macrophage regulatory network helps accelerate the de-activation of the downstream gene (Z), after the upstream inputs have switched to a low state. In this way, the FFL motif would ameliorate the effect of long protein half lives in mammalian cells.
7. Material and methods
The human IκBα cDNA was amplified from peripheral-blood leucocytes and cloned into the pEF6/V5-His vector (Invitrogen), in frame with EGFP (Clontech), and expressed under the control of the EF1-a promoter. This construct was electroporated into RAW 264.7 macrophage cells and a high-expressing clone (#28) was selected by flow cytometry and expanded. Cells were cultured in RPMI 1640/10% FCS including 25 mM HEPES and stimulated with 1 μg ml−1 pI.pC dsRNA (Sigma). Time course images were taken on a Leica TCS–SP2 laser-scanning confocal microscope. Fluorescence was excited at 488 nm and emission measured from 500 to 550 nm at 40× magnification with the detection pinhole at 600 μm diameter to give a 14 μm depth-of-field.
Whole-cell fluorescence was measured as follows. A mask was created using the entire series of images to define the extent of cellular motion over the course of the experiment. All contiguous pixels which were more than 10% above the background in any image in the series were included in the mask region. After background subtraction, the total flux within this boundary was taken as a measure of the total amount of IκBα-GFP within each cell.
All simulations were carried out with the Dizzy kinetic simulation modelling environment (Ramsey et al. 2005). An ensemble size of 7500 was used for the stochastic simulations. The stochastic simulations were carried out using the Gibson–Bruck simulator option, on an Intel Pentium 4 workstation running the Sun Java virtual machine v.1.4.2 under Fedora Core 1 Linux with kernel v.2.4.22. Steady-state was obtained using a fifth order Dormand–Prince ODE solver with fourth order error estimation. The relative and absolute error tolerances for the ODE solver were 0.0001.
8. Dizzy stochastic simulation software
Dizzy (Ramsey et al. 2005) is a stochastic simulation software package and model description language for systems of interacting biochemical species. Dizzy's simulation engine includes Gillespie's original exact stochastic algorithm (Gillespie 1976), the Gibson–Bruck optimization of this algorithm (Gibson & Bruck 2000), and the ‘Tau-leap’ approximate accelerated algorithm (Gillespie 2001). Dizzy is a stable software package which we have used extensively over the past year to model and simulate multi-gene networks stochastically (de Atauri et al. 2005). Dizzy is written in the Java programming language, and will run on any computing platform that has a Java 1.4 runtime environment. It is available under a free software and open-source license (the Lesser General Public License) and can be found on our group's web page (magnet.systemsbiology.net/dizzy).
An especially pertinent feature of Dizzy is that it allows user-defined library elements: models can be constructed in a modular, hierarchical fashion. We used this feature of Dizzy to construct an organism-specific, parameterized model of all the key steps in gene expression. Changing the parameters of the model allows us to model genes with different characteristics such as: number of regulatory factors, their concentrations over time, and the nature of their kinetic interactions; the rate of RNA polymerase recruitment and transcription initiation by the transcription factor complex; the rate of transcription along DNA; the length of transcribed DNA; the length of the mRNA, mRNA half-life; the rate of ribosome recruitment and initiation of translation; the rate of translocation along the mRNA; and protein half life. A second feature of Dizzy is the ability to estimate, using ODE simulation and the symbolic Jacobian matrix, the steady-state stochastic fluctuations for all species in the model (Orrell et al. 2005). This feature has been successfully applied to estimate the steady-state fluctuations in a nonlinear biochemical model (with 55 reaction channels and 26 chemical species) describing the galactose uptake pathway in yeast.
This work was supported in part by grant #10830302 from the National Institute of Allergy and Infectious Disease.
- © 2006 The Royal Society