Effect of Promoter Architecture on the Cell-to-Cell Variability in Gene Expression

According to recent experimental evidence, promoter architecture, defined by the number, strength and regulatory role of the operators that control transcription, plays a major role in determining the level of cell-to-cell variability in gene expression. These quantitative experiments call for a corresponding modeling effort that addresses the question of how changes in promoter architecture affect variability in gene expression in a systematic rather than case-by-case fashion. In this article we make such a systematic investigation, based on a microscopic model of gene regulation that incorporates stochastic effects. In particular, we show how operator strength and operator multiplicity affect this variability. We examine different modes of transcription factor binding to complex promoters (cooperative, independent, simultaneous) and how each of these affects the level of variability in transcriptional output from cell-to-cell. We propose that direct comparison between in vivo single-cell experiments and theoretical predictions for the moments of the probability distribution of mRNA number per cell can be used to test kinetic models of gene regulation. The emphasis of the discussion is on prokaryotic gene regulation, but our analysis can be extended to eukaryotic cells as well.


Introduction
A fundamental property of all living organisms is their ability to gather information about their environment and adjust their internal physiological state in response to environmental conditions. This property, shared by all organisms, includes the ability of single-cells to respond to changes in their environment by regulating their patterns of gene expression. By regulating the genes they express, cells are able to survive, for example, changes in the extracellular pH or osmotic pressure, switch the mode of sugar utilization when the sugar content in their medium changes, or respond to shortages in key metabolites by adapting their metabolic pathways. Perhaps more interesting is the organization of patterns of gene expression in space and time resulting in the differentiation of cells into different types, which is one of the defining features of multicellular organisms. Much of this regulation occurs at the level of transcription initiation, and is mediated by simple physical interactions between transcription factor proteins and DNA, leading to genes being turned on or off. Understanding how genes are turned on or off (as well as the more nuanced expression patterns in which the level of expression takes intermediate levels) at a mechanistic level has been one of the great challenges of molecular biology and has attracted intense attention over the past 50 years.
The current view of transcription and transcriptional regulation has been strongly influenced by recent experiments with single-cell and single-molecule resolution [1][2][3][4][5][6][7][8][9][10][11]. These experiments have confirmed the long-suspected idea that gene expression is stochastic [12,13], meaning that different steps on the path from gene to protein occur at random. This stochasticity also causes variability in the number of messenger RNAs (mRNA) and proteins produced from cell-to-cell in a colony of isogenic cells [11,[14][15][16][17]. The question of how transcriptional regulatory networks function reliably in spite of the noisy character of the inputs and outputs has attracted much experimental and theoretical interest [18,19]. A different, but also very relevant, question is whether cells actually exploit this stochasticity to fulfill any physiologically important task. This issue has been investigated in many different cell types and it has been found that stochasticity in gene expression plays a pivotal role in processes as diverse as cell fate determination in the retina of Drosophila melanogaster [20], entrance to the competent state of B. subtilis [7], resistance of yeast colonies to antibiotic challenge [17], maintenance of HIV latency [21], promoting host infection by pathogens [22] or the induction of the lactose operon in E. coli [23]. Other examples have been found, and reviewed elsewhere [24,25]. The overall conclusion of all of these studies is that stochasticity in gene expression can have important physiological consequences in natural and synthetic systems and that the overall architecture of the gene regulatory network can greatly affect the level of stochasticity.
A number of theoretical and experimental studies have revealed multiple ways in which the architecture of the gene regulatory network affects cell-to-cell variability in gene expression. Examples of mechanisms for the control of stochasticity have been proposed and tested, including the regulation of translational efficiency [8], the presence of negative feedback loops [26,27,28], or the propagation of fluctuations from upstream regulatory components [29]. Another important source of stochasticity in gene expression is fluctuations in promoter activity, caused by stochastic association and dissociation of transcription factors, chromatin remodeling events, and formation of stable pre-initiation complexes [5,15,16,23,30]. In particular, it has been reported that perturbations to the architecture of yeast and bacterial promoters, such as varying the strength of transcription factor binding sites [17], the number and location of such binding sites [11,31], the presence of auxiliary operators that mediate DNA looping [23], or the competition of activators and repressors for binding to the same stretch of DNA associated with the promoter [32], may strongly affect the level of variability.
Our goal is to examine all of these different promoter architectures from a unifying perspective provided by stochastic models of transcription leading to mRNA production. The logic here is the same as in earlier work where we examined a host of different promoter architectures using thermodynamic models of transcriptional regulation [33,34]. We now generalize those systematic efforts to examine the same architectures, but now from the point of view of stochastic models. These models allow us to assess the unique signature provided by a particular regulatory architecture in terms of the cell-to-cell variability it produces.
First, we investigate in general theoretical terms how the architecture of a promoter affects the level of cell-to-cell variability. The architecture of a promoter is defined by the collection of transcription factor binding sites (also known as operators), their number, position within the promoter, their strength, as well as what kind of transcription factors bind them (repressors, activators or both), and how those transcription factors bind to the operators (independently, cooperatively, simultaneously). We apply the master-equation model of stochastic gene expression [35,36,37,38] to increasingly complex promoter architectures [30], and compute the moments of the mRNA and protein distributions expected for these promoters. Our results provide an expectation for how different architectural elements affect cell-to-cell variability in gene expression.
The second point of this paper is to make use of stochastic kinetic models of gene regulation to put forth in vivo tests of the molecular mechanisms of gene regulation by transcription factors that have been proposed as a result of in vitro biochemical experiments. The idea of using spontaneous fluctuations in gene expression to infer properties of gene regulatory circuits is an area of growing interest, given its non-invasive nature and its potential to reveal regulatory mechanisms in vivo. Different theoretical methods have recently been proposed, which could be employed to distinguish between different modes (e.g. AND/OR) of combinatorial gene regulation, and to rule out candidate regulatory circuits [27,39,40] based solely on properties of noise in gene expression, such as the autocorrelation function of the fluctuations [27] or the three-point steady state correlations between multiple inputs and outputs [39,40].
Here, we make experimentally testable predictions about the level of cell-to-cell variability in gene expression expected for different bacterial promoters, based on the physical kinetic models of gene regulation that are believed to describe these promoters in vivo. In particular, we focus on how varying the different parameters (i.e., mutating operators to make them stronger or weaker, varying the intracellular concentration of transcription factors, etc.) should affect the level of variability. This way, cell-tocell variability in gene expression is used as a tool for testing kinetic models of transcription factor mediated regulation of gene expression in vivo.
The remainder of the paper is organized as follows: First we describe the theoretical formalism we use to determine analytic expressions for the moments of the probability distribution for both mRNA and protein abundances per cell. Next, we examine how the architecture of the promoter affects cell-to-cell variability in gene expression. We focus on simple and cooperative repression, simple and cooperative activation, and transcriptional regulation by distal operators mediated by DNA looping. We investigate how noise in gene expression caused by promoter activation differs from repression, how operator multiplicity affects noise in gene expression, the effect of cooperative binding of transcription factors, as well as DNA looping. For each one of these architectures we present a prediction of cell-to-cell variability in gene expression for a bacterial promoter that has been well characterized experimentally in terms of their mean expression values. These predictions suggest a new round of experiments to test the current mechanistic models of gene regulation at these promoters.

Methods
In order to investigate how promoter architecture affects cell-tocell variability in gene expression, we use a model based on classical chemical kinetics (illustrated in Figure 1A), in which a promoter containing multiple operators may exist in as many biochemical states as allowed by the combinatorial binding of transcription factors to its operators. The promoter transitions stochastically between the different states as transcription factors bind and fall off. Synthesis of mRNA is assumed to occur stochastically at a constant rate that is different for each promoter state. Further, transcripts are assumed to be degraded at a constant rate per molecule. This kind of model is the kinetic counterpart of the so-called ''thermodynamic model'' of transcriptional regulation [41], and it is the standard framework for interpreting the kinetics of gene

Author Summary
Stochastic chemical kinetics provides a framework for modeling gene regulation at the single-cell level. Using this framework, we systematically investigate the effect of promoter architecture, that is, the number, quality and position of transcription factor binding sites, on cell-to-cell variability in transcription levels. We compare architectures resulting in transcriptional activation with those resulting in transcriptional repression. We start from simple activation and repression motifs with a single operator sequence, and explore the parameter regime for which the cell-to-cell variability is maximal. Using the same formalism, we then turn to more complicated architectures with more than one operator. We examine the effect of independent and cooperative binding, as well as the role of DNA mechanics for those architectures where DNA looping is relevant. We examine the interplay between operator strength and operator number, and we make specific predictions for single-cell mRNA-counting experiments with well characterized promoters. This theoretical approach makes it possible to find the statistical response of a population of cells to perturbations in the architecture of the promoter; it can be used to quantitatively test physical models of gene regulation in vivo, and as the basis of a more systematic approach to designing new promoter architectures.
regulation in biochemical experiments, both in vivo [2,23] and in vitro [42,43]. This class of kinetic models can easily accommodate stochastic effects, and it leads to a master equation from which the probability distribution of mRNA and protein copy number per cell can be computed. It is often referred to as the standard model of stochastic gene expression [38,44,45]. The degree of cell-to-cell variability in gene expression can be quantified by the stationary variance, defined as the ratio of the standard deviation and the mean of the probability distribution of mRNA or protein copy number per cell [35], or else by the Fano factor, the ratio between the variance and the mean. These two are the two most common metrics of noise in gene expression, and the relation between them will be discussed later.
In order to compute the noise strength from this class of models, we follow the same approach as in a previous article [30], which extends a master equation derived elsewhere [36,37,46] to promoters with arbitrary combinatorial complexity. The complexity refers to the existence of a number of discrete promoter states corresponding to different arrangements of transcription factors on the promoter DNA. Promoter dynamics are described by trajectories involving stochastic transitions between promoter states which are induced by the binding and unbinding of transcription factors. A detailed derivation of the equations which describe promoter dynamics can be found in the Text S1, but the essentials are described below.
There are only two stochastic variables in the model: the number of mRNA transcripts per cell, which is represented by the unitless state variable m, and the state of the promoter, which is defined by the pattern of transcription factors bound to their operator sites. The promoter state is described by a discrete and finite stochastic variable (s) (for an example, see Figure 1A). The example in Figure 1A illustrates the simplest model of transcriptional activation by a transcription factor. When the activator is not bound (state 1), mRNA is synthesized at rate r 1 . When the activator is bound to the promoter (state 2), mRNA is synthesized at the higher rate r 2 . The promoter switches stochastically from Simple two-state bacterial promoter undergoing stochastic activation by a transcriptional activator binding to a single operator site. The rates of activator association and dissociation are given by k on A and k off A , respectively and the rates of mRNA production for the basal and active states are r 1 and r 2 respectively. The mRNA degradation rate is assumed to be constant for each molecule, and is given by the parameter c. (B) List of all possible stochastic transitions affecting either the copy number of mRNA (m) or the state of the promoter (s) and their respective statistical weight. State 1 has the operator free. State 2 is the activator bound state. The weights represent the probability that each change of state will occur during a time increment Dt. The master equation is constructed based on these rules. doi:10.1371/journal.pcbi.1001100.g001 state 1 to state 2 with rate k on A , and from state 2 to state 1 with rate k off A . Each mRNA molecule is degraded with rate c. The time evolution for the joint probability of having the promoter in states 1 or 2, with m mRNAs in the cell (which we write as p(1,m) and p(2,m), respectively), is given by a master equation, which we can build by listing all possible reactions that lead to a change in cellular state, either by changing m or by changing s (Figure 1b) we can rewrite the system of equations (1)  This has several advantages, but the most important one is that the matrix approach reduces the task of obtaining analytical expressions for the moments of the steady state mRNA distribution for an arbitrarily complex promoter to solving two simple linear matrix equations (more details are given in the Text S1).
The matrices appearing in equation (4) all have simple and intuitive interpretations. The matrixK K describes the stochastic transitions between promoter states: The off-diagonal elements of the matrixK K ij are the rates of making transitions from promoter state (j)to promoter state (i).The diagonal elements of the matrix K K jj are negative, and they represent the net probability flux out of state (j):K K jj~{ P i=jK K ij . The matrixR R is a diagonal matrix whose elementR R jj gives the rate of transcription initiation when the promoter is in state (j). Finally, the matrixÎ I is the identity matrix. An example of matricesK K andR R is presented pictorially in Figure 1 in Text S1. It is straightforward to see that even though equation (4) has been derived for a two-state promoter, it also applies to any other promoter architecture. What will change for different architectures are the dimensions of the matrices and vectors (these are given by the number of promoter states) as well as the values of the rate constants that make up the matrix elements of the various matrices.
An important limit of the master equation, which is often attained experimentally, is the steady state limit, where the probability distribution for mRNA number per cell does not change with time. Although the time dependence of the moments of the mRNA distribution can be easily computed from our model, for the sake of simplicity and because most experimental studies have been performed on cells in steady state, we focus on this limit. As shown in Text S1, analytic expressions for the first two moments of the steady state mRNA probability distribution are found by multiplying both sides of equation (4) by m and m 2 respectively, and then summing m from 0 to infinity. After some algebra (elaborated in an earlier paper and in Text S1), we find that the first two moments can be written as: The vectorr r contains the ordered list of rates of transcription initiation for each promoter state. For the two-state promoter shown in Figure 1,r r~r 1 ,r 2 ð Þ. The vectorm m (0) contains the steady state probabilities for finding the promoter in each one of the possible promoter states, whilem m (1) is the steady-state mean mRNA number in each promoter state. The vectorm m (0) is the solution to the matrix equation while the vectorm m (1) is obtained from Construct the matricesK K andR R, and the vectorr r, (Figure 1 in Text S1).
The normalized variance of the mRNA distribution in steady state is then computed from the equation: Equation (9) reveals that, regardless of the specific details characterizing promoter architecture, the intrinsic noise is always the sum of two components, and it can be written as The first component is due to spontaneous stochastic production and degradation of single mRNA molecules, it is always equal to the Poissonian expectation of 1=SmT, and is independent of the architecture of the promoter. For an unregulated promoter that is always active and does not switch between multiple states (or does so very fast compared to the rates of transcription and mRNA degradation), the mRNA distribution is well described by a Poisson distribution [45,47], and the normalized variance is equal to 1=SmT. The second component (''promoter noise'') results from promoter state fluctuations, and captures the effect of the promoter's architecture on the cell-to-cell variability in mRNA: In order to quantify the effect of the promoter architecture in the level of cell-to-cell variability in mRNA expression, we define the deviation in the normalized variance caused by gene regulation relative to the baseline Poisson noise for the same mean (see Figure 2): Therefore, the deviation in the normalized variance caused by gene regulation is equal to the ratio between the variance and the mean. This parameter is also known as the Fano factor. Thus, for any given promoter architecture, the Fano factor quantitatively characterizes how large the mRNA noise is relative to that of a Poisson distribution of the same mean (i.e. how much the noise for the regulated promoter elevates with respect to the Poisson noise). This is the parameter that we will use throughout the paper as the metric of cell-to-cell variability in gene expression.

Promoter noise and variability of mRNA and protein numbers
For proteins, the picture is only slightly more complicated. As shown in the Text S1, in the limit where the lifetime of mRNA is much shorter than that of the protein it encodes for (a limit that is often fulfilled [30]), the noise strength of the probability distribution of proteins per cell takes the following form (where we define n as a state variable that represents the copy number of proteins per cell): where d stands for the protein degradation rate, and the constant b is equal to the protein burst size (the average number of proteins produced by one mRNA molecule). The mean protein per cell is given by and the vectorñ n (1) is the solution of the algebraic equation: The reader is referred to the Text S1 for a detailed derivation and interpretation of these equations. In the previous section we have shown that the noise for proteins and mRNA take very similar analytical forms. Indeed, if we definer r n~br r andR R n~bR R, as the vector and matrix containing the average rates of protein synthesis for each promoter state, it is straightforward to see that equations (8) and (15) are mathematically equivalent, with the only difference being that in equation (15) the matrixR R n represents the rates of protein synthesis, so all the rates of transcription are multiplied by the translation burst size b. Therefore, the vectors m m (1) andñ n (1) are only going to differ in the prefactor b multiplying all the different transcription rates. We conclude that the promoter contribution to the noise takes the exact same analytical form both for proteins and for mRNA, with the only other quantitative difference being the different rates of degradation for proteins and mRNA. Therefore, promoter architecture has the same qualitative effect on cell-to-cell variability in mRNA and protein numbers. All the conclusions about the effect of promoter architecture on cellto-cell variability in mRNA expression are also valid for proteins, even though quantitative differences do generally exist. For the sake of simplicity we focus on mRNA noise for the remainder of the paper.

Parameters and assumptions
In order to evaluate the equations in our model, we use parameters that are consistent with experimental measurements of rates and equilibrium constants in vivo and in vitro, which we summarize in Table 1. Although these values correspond to specific examples of E. coli promoters, like the P lac or the P RM promoter, we extend their reach by using them as ''typical'' parameters characteristic of bacterial promoters, with the idea being that we are trying to demonstrate the classes of effects that can be expected, rather than dissecting in detail any particular promoter. The rate of association for transcription factors to operators in vivo is assumed to be the same as the recently measured value for the Lac repressor, which is close to the diffusion limited rate [48]. In order to test whether the particular selection of parameters in Table 1 is biasing our results, we have also done several controls (See Figures 2-4 in Text S1) in which the kinetic parameters were randomly sampled. We found that the conclusions reached for the set of parameters in Table 1 are valid for other parameter sets as well.
Operator strength reflects how tightly operators bind their transcription factors, and it is quantitatively characterized by the equilibrium dissociation constant K O{TF . The dissociation constant has units of concentration and is equal to the concentration of free transcription factor at which the probability for the operator to be occupied is 1/2. K O{TF is related to the association and dissociation rates by K O{TF~koff k 0 on , where k off is the rate (i.e., the probability per unit time) at which a transcription factor dissociates from the promoter, and k 0 on is a second order rate constant, which represents the association rate per unit of concentration of transcription factors, i.e., k on~k 0 on N TF ½ . Note that in the last formula k on , which has units of s 21 , is written as the product of two quantities: ½N TF , which is the concentration (in units of (mol/liter)) of transcription factors inside the cell, and k 0 on , a second order rate constant that has units of (mol/ liter) 21 s 21 . For simplicity, we assume that the binding reaction is diffusion limited, namely, k 0 on is already close to its maximum possible value, so the only parameter that can differ from operator to operator is the dissociation rate: strong operators have slow dissociation rates, and weak operators have large dissociation rates.
Throughout this paper, we also make the assumption that the mean expression level is controlled by varying the intracellular concentration of transcription factors, a scenario that is very  [49,50,51]. We also assume that changing the intracellular concentration of transcription factors only affects the association rate of transcription factors to the operators, but the dissociation rate and the rates of transcription at each promoter state are not affected. In other words, k off is a constant parameter for each operator, and it is not changed when we change the mean by titrating the intracellular repressor level. All of these general assumptions need to be revisited when studying a specific gene-regulatory system. Here our focus is on illustrating the general principles associated with different promoter architectures typical of those found in prokaryotes.

Simulations
To generate mRNA time traces, we applied the Gillespie algorithm [52] to the master equation described in the text. A single time step of the simulation is performed as follows: one of the set of possible trajectories is chosen according to its relative weight, and the state of the system is updated appropriately. At the same time, the time elapsed since the last step is chosen from an exponential distribution, whose rate parameter equals the sum of rate parameters of all possible trajectories. This process is repeated iteratively to generate trajectories that exactly reflect dynamics of the underlying master equation. For the figures, simulation lengths were set long enough for the system to reach steady state and for a few promoter state transitions to occur.
To generate the probability distributions, it is convenient to reformulate the entire system of mRNA master equations in terms of a single matrix equation. To do this, we first define a vectorP where p(i,m)is the joint probability of having mmRNAs while in the ith promoter state. Then the master equation for time evolution of this probability vector is : : mRNA near the population average. In contrast, the stronger operator with a slower repressor dissociation rate, leads to larger promoter noise and strongly non-Poissonian mRNA statistics. (B) Kinetic mechanism of repression for an architecture involving a single repressor binding site. The repressor turns off the gene when it binds to the promoter (with rate k on R ), and transcription occurs at a constant rate r when the repressor falls off (with rate k off R ). (C) Normalized variance as a function of the fold-change in mean mRNA copy number. The parameters used are drawn from Table 1. The value of k off R~0 :0023s {1 from Table 1 corresponds to the in vitro dissociation constant of the Lac repressor from the Oid operator (black). The results for an offrate 10-times higher are also plotted (red). As a reference for the size of the fluctuations, we show the normalized variance for a Poisson promoter. (D) Fano factor for two promoters bearing the same off-rates as in (B). Inset. Prediction for the Fano factor for the D O3 D O2 P lacUV5 promoter, a variant of the P lacUV5 promoter for which the two auxiliary operators have been deleted. The fold-change in mRNA noise is plotted as a function of the fold-change in mean mRNA copy number for mutants of the promoter that replace O1 for Oid, O2 or O3. The parameters are taken from Table 1 and [33]. Lifetimes of the operator-repressor complex are 7 min for Oid, 2.4 min for O1, 11s for O2 and 0.47 s for O3. (E) Fold-change in protein noise as a function of the foldchange in mean expression. As expected, the effect of operator strength is the same as observed for mRNA noise. doi:10.1371/journal.pcbi.1001100.g002 Table 1. Kinetic parameters used to make the quantitative estimates in the text and plots in the figures.

Kinetic Rate Symbol Value Reference
Unregulated promoter transcription rate r 0.33s 21 [99] Repressor and activator associations rates Repressor and activator dissociation rates mRNA decay rate c 0:011 s {1 [10] Ratio between transcription rates due to activation f~r 1 =r 2 11 [50] Cooperativity in repression V repression 0.013 [50] Cooperativity in activation V activation 0:1 [33] Looping J-factor ½J 660 nM [33] Protein translation burst size b 31.2 proteins/mRNA [5] Protein decay rate ½J 0.00083s 21 [99] These parameters are all measured for model systems such as the P lac promoter or the P RM in E. coli, and are here considered representative for promoter-transcription factor interactions. doi:10.1371/journal.pcbi.1001100.t001 where each element of the matrix is itself an N by N matrix as described in the text. Then finding the steady-state distributioñ P P SS is equivalent to finding the eigenvector of the above matrix associated with eigenvalue 0. To perform this calculation numerically, one must first choose an upper bound on mRNA copy number in order to work with finite matrices. In this work, we chose an upper bound six standard deviations above mean mRNA copy number as an initial guess, and then modified this bound if necessary. Computations were performed using the SciPy (Scientific Python) software package.

Promoters with a single repressor binding site
We first investigate a promoter architecture consisting of a single repressor binding site, and examine how operator strength affects intrinsic variability in gene expression. Although this particular mode of gene regulation has been well studied theoretically before [1,16,36,37,45], it is a useful starting point for illustrating the utility of this class of models. Within this class of models, when the repressor is bound to the operator, it interferes with transcription initiation and transcription does not occur. When the repressor dissociates and the operator is free, RNAP can bind and initiate transcription at a constant rate r. The probability per unit time that a bound repressor dissociates is k off R , and the probability per unit time that a free repressor binds the empty operator is k on R~k 0 on N R ½ , where k 0 on is the second-order association constant and N R ½ is the intracellular repressor concentration. The rate of mRNA degradation per molecule is c. This mechanism is illustrated in Figure 2B.
We compute the mean and the Fano factor for this architecture following the algorithm described in the Mathematical Methods section. The kinetic rate and transcription rate matricesK K andR R are shown in Table S1 in Text S1. For this simple architecture, the mean of the mRNA probability distribution and the normalized variance take simple analytical forms: Using the relationship between k on R and the intracellular concentration of repressor, we can write the mean as: Here we have defined the equilibrium dissociation constant between the repressor and the operator as: k 0 on . It is interesting to note that equation (20) could have been derived repression for a dual-repression architecture. The parameters k off R and k on R are the rates of repressor dissociation and association to the operators, and V is a parameter reflecting the effect of cooperative binding on the dissociation rate. For independent binding, V~1 and for cooperative binding V~0:013 (see Table 1). (B) Fold-change in the mRNA noise caused by gene regulation for independent (red) and cooperative (black) repression as a function of the mean mRNA copy number. Inset: Prediction for a variant of the l P R promoter where the upstream operators O L1 , O L2 and O L3 are deleted. The promoter mRNA noise is plotted as a function of the mean mRNA number for both wild-type cI repressor (blue line) and a repressor mutant (Y210H) that abolishes cooperativity (red line). Parameters taken from [43,97]. The lifetime of the O R1 -cI complex is 4 min. Lifetime of O R2 -cI complex is 9.5s. (C) mRNA distribution for the same parameters used in (B). doi:10.1371/journal.pcbi.1001100.g003 The Fano factor is plotted as a function of the fold-change in the mean expression level, in the absence (blue) and presence (black) of the auxiliary operator, and assuming that dissociation of the operator from Om is the same in the looped and the unlooped state (c = 1). The presence of the auxiliary operator, which enables repression by DNA looping, increases the cell-to-cell variability. The regions over which the state with two repressors bound, the state with one repressor bound, or the looped DNA state are dominant are indicated by the shading in the background. The noise is larger at intermediate repression levels, where only one repressor is found bound to the promoter region, simultaneously occupying the auxiliary and main operators through DNA looping. The rate of DNA loop formation is k loop~6 60nM ð Þ k 0 R [33]. We also show the effect of DNA looping in the case where the kinetics of dissociation from the looped state are 100 times faster than the kinetics of dissociation from the unlooped state: c~k unloop . k off R~1 00(red). In this limit, the presence of the auxiliary operator leads to less gene expression noise. (C) Prediction for a library of P lacUV5 promoter variants, harboring an O2 deletion, and with the position of O3 moved upstream by multiples of 11 bp while keeping its identity (red), or replaced by the operator by Oid (black). Parameters are taken from the analysis in [33] of the data in [98]. We assume a concentration of 50 Lac repressor tetramers per cell. The association rate of the tetrameric repressor to the operators is taken from Table 1. The lifetimes of the operator-repressor complex are given in the caption to Figure 2. The dependence of the rate of DNA looping on the inter-operator distance is taken from [33], and equal to: using the thermodynamic model approach [33,34,41,53]. In particular we see that this expression is equal to the product of the maximal activity in the absence of repressor SmT max~r =c, and the so-called fold-change in gene expression: 1zk on [34]. The fold-change is defined as the ratio of the level of expression in the presence of the transcription factor of interest, and the level of expression in the absence of the transcription factor.
The Fano factor for the mRNA distribution can be computed from equation (12) and we obtain: which is also shown as the first entry of Table S2 in Text S1. In many experiments [4,15,31,50], the concentration of repressor inside the cell N R ½ (and therefore the association rate k on R~k 0 on N R ½ ) can be varied by either expressing the repressor from an inducible promoter, or by adding an inducer that binds directly to the repressor rendering it incapable of binding specifically to the operators in the promoter region. When such an operation is performed, the only parameter that is varied is typically k on R , and all other kinetic rates are constant. The Fano factor can thus be rewritten as a function of the mean mRNA, and we find: Therefore, for any given value of the mean, the Fano factor depends only on two parameters: the maximal mRNA or protein expression per cell, and a parameter that reflects the strength of binding between the repressor and the operator: k off R .
c. Equations (20) and (21) reveal that changes in the mean due to repressor titration affect the noise as well as the mean. Since neither the repressor dissociation rate k off R nor the mRNA degradation rates are affected by the concentration of repressors, c is a constant parameter that will determine how large the cell-to-cell variability is: The Fano factor is maximal for promoters with very strong operators, (k off R ,,c, and it goes to 1 (i.e., the distribution tends to a Poisson distribution) when the operator is very weak and the rate of dissociation extremely fast (k off R ..c). In the latter limit of fast promoter kinetics, the fast fluctuations in promoter occupancy are filtered by the long lifetime of mRNA. Effectively, mRNA degradation acts as a low-pass frequency filter [54,55], and fast fluctuations in promoter occupancy are not propagated into mRNA fluctuations. Therefore, promoters with strong operators are expected to be noisier than promoters with weak operators [56]. From this discussion it should also be clear that the mRNA degradation rate critically affects cell-to-cell variability. Any processes that tend to accelerate degradation will tend to increase noise, and mRNA stabilization (i.e., protection of the transcript by RNA binding proteins) leads to reduction of variability. However, the focus of this article is on promoter architecture and transcriptional regulation. Therefore, we do not consider regulation of transcription by mRNA degradation, and assume that all the promoters transcribe the same mRNA as is often the case in experimental studies.
The effect of operator strength on the output of transcription and translation is illustrated in Figure 2A, where we show results from a stochastic simulation of the model depicted in Figure 2B, for the case of a weak and a strong operator. The simulation yields trajectories in time for the promoter state, the mRNA, and protein number, as well as the steady state distribution of mRNA number. Concentrations of repressor in the simulations were chosen so that the mean expression level was equal for the two different promoter architectures. As expected from the general arguments presented above, we clearly see that the level of variability is smaller for the weak operator than for the strong operator, due to faster promoter switching leading to smaller mRNA fluctuations and a more Poisson-like mRNA distribution (Figure 2A, weak promoter). Slow dissociation from a strong operator, on the other hand, causes slow promoter state fluctuations and a highly non-Poissonian mRNA distribution, with few cells near the mean expression level (see Figure 2A, strong promoter).
In order to show that the effect of operator strength on the cellto-cell variability is general and does not depend on the particular set of parameters chosen in the simulation, in figures 2C and 2D, we show the normalized variance and the Fano factor as a function of the fold-change in the mean mRNA concentration for a strong operator whose dissociation rate is k off R~0 :0027s {1 (a value that is representative of well characterized repressoroperator interactions such as the Lac repressor-lacOid, or the cI2-lO R1 ), and for a single weak operator whose dissociation rate k off R is 10 times larger.
The Fano factor has a characteristic shape whereby it takes values approaching 1 at low and high transcription levels with a peak at intermediate values. The reason for this shape is that for very low transcription levels the promoter is nearly always inactive, firing only very rarely. In this limit successive transcription events become uncorrelated and the time in between them is exponentially distributed, leading to a distribution of mRNA per cell that approaches a Poisson distribution characterized by a Fano factor equal to 1. In contrast, for very high transcription levels the promoter is nearly always active, switching off very rarely and staying in the off state for short times. In this limit, transcription events are again uncorrelated and exponentially distributed, leading once again to a Poisson distribution of mRNA number. It is only for intermediate values of the mean that the promoter is switching between a transcriptionally active and an inactive state. This causes transcription to occur in bursts, and the mRNA distribution to deviate from Poisson, leading to a Fano factor that is larger than 1.
In Figure 2E we plot the fold-change in protein noise due to gene regulation for the simple repression architecture. As expected, we find that the effect of operator strength in protein noise is qualitatively identical to what we found for mRNA. Since the same can be said of all the rest of architectures studied, we will limit the discussion to mRNA noise for the rest of the paper, with the understanding that for the class of models considered here, all the conclusions about the effect of promoter architecture in cell-tocell variability that are valid for mRNA, are true for intrinsic protein noise as well.
In Figure 2, and throughout this paper, we plot the Fano factor as a function of transcription level, which is characterized by the fold-change in gene expression. The fold-change in gene expression is defined as the mean mRNA number in the presence of the transcription factor, normalized by the mean mRNA in the absence of the transcription factor. For architectures based on repression, the fold-change in gene expression is always less than 1, since the repressor reduces the level of transcription. For example, a fold-change in gene expression of 0.1 means that in the presence of repressor, the transcription level is 10% of the value it would have if the repressor concentration dropped to 0. For the case of activators, the fold-change is always greater than 1, since activators raise the level of transcription.
An example of the single repressor-binding site architecture is a simplified version of the P lacUV5 promoter, which consists of a single operator overlapping with the promoter. Based on a simple kinetic model of repression, in which the Lac repressor competes with RNAP for binding at the promoter, we can write down theK K andR R matrices and compute the cell-to-cell variability in mRNA copy number. The matrices are presented in Table S1 in Text S1. Based on our previous analysis, we know that stronger operators are expected to cause larger noise and higher values of the Fano factor than weaker operators. Therefore, we expect that if we replace the wild-type O1 operator by the 10 times weaker O2 operator, or by the ,500 times weaker operator O3, the foldchange in noise should go down. Using our best estimates and available measurements for the kinetic parameters involved, we find that noise is indeed much larger for O1 than for O2, and it is negligible for O3. This prediction is presented as an inset in Figure 2C.

Promoters with two repressor-binding operators
Dual repression occurs when promoters contain two or more repressor binding sites. Here, we consider three different scenarios for architectures with two operators: 1) repressors bind independently to the two operators, 2) repressors bind cooperatively to the two operators and 3) one single repressor may be bound to the two operators simultaneously thereby looping the intervening DNA. At the molecular level, cooperative repression is achieved by two weak operators that form long-lived repressor-bound complexes when both operators are simultaneously occupied. Transcription factors may stabilize each other either through direct proteinprotein interactions [53], or through indirect mechanisms mediated by alteration of DNA conformation [57].
Cooperative and independent repression. The kinetic mechanisms of gene repression for both the cooperative and independent repressor architectures are reproduced in Figure 3A. For simplicity, we assume that both sites are of equal strength, so the rates of association and dissociation to both sites are equal. Cooperative binding is reflected in the fact that the rate of dissociation from the state where the two operators are occupied is slower (by a factor Vvv1) than the dissociation from a single operator. This parameter is related to the cooperativity factor v often found in thermodynamic models [54] by V~1=v. A typical value of , 10 {3 {10 {2 [50,53]. By way of contrast, independent binding is characterized by a value of V~1, which reflects the fact that the rate of dissociation from each operator is not affected by the presence of the other operator.
TheK K andR R matrices for these two architectures are defined in Table S1 in Text S1. Using these matrices, we can compute the mean gene expression and the Fano factor for these two architectures as a function of the concentrations of repressor. The resulting expression for the fold-change in noise is shown as entry number 3 of Table S2 in Text S1. As shown in Figure 3B, the noise for cooperative repression is substantially larger than for the independent repression architecture. The high levels of intrinsic noise associated with cooperative repression can be understood intuitively in terms of the kinetics of repressor-operator interactions. At low repressor concentration, the lifetime of the states where only one repressor is bound to either one of the two operators can be shorter than the time it takes for a second repressor to bind. This makes simultaneous binding of two repressors to the two operators a rare event. However, when it occurs, the two repressors stabilize each other, forming a very long-lived complex with the operator DNA. This mode of repression, with rare but long-lived repression events, is intrinsically very noisy, since the promoter switches slowly between active (unrepressed) and inactive (repressed) states, generating wide bimodal distributions of mRNA (see Figure 3C). On the other hand, independent binding to two operators causes more frequent transitions between repressed and unrepressed states, leading to lower levels of intrinsic noise and long-tailed mRNA distributions (see Figure 3C). In order to illustrate these conclusions, we have evaluated the model with a specific parameter set that is representative of this kind of bacterial promoters, and plotted the Fano factor as a function of the mean, under the assumption that we vary the mean by titrating the amount of repressor inside the cell. Furthermore, so as to demonstrate that our conclusions are not dependent on choice of parameters, we have randomly generated 10,000 different sets of kinetic parameters and compared the Fano factor for cooperative and independent binding. The result of this analysis is shown in Figure 2 in Text S1, where we demonstrate that cooperative binding always results in larger cell-to-cell variability than non-cooperative binding.
As an example of the two repressor-binding sites architecture, we consider a simplified version of the lytic phage-l P R promoter, which is controlled by the lysogenic repressor cI. The wild-type P R promoter consists of three proximal repressor binding sites, O R1 , O R2 and O R3 , with different affinities for the repressor (O R2 is ,25 times weaker than O R1 ) [58], and three distal operators O L1 , O L2 and O R3 . For simplicity, we consider a simpler version of P R , harboring a deletion of the three distal operators. In the absence of these operators, the O R3 operator plays only a very minor role in the repression of this promoter, and it can be ignored [50,59]. We are then left with only O R1 and O R2 . The cI repressor binds cooperatively to O R1 and O R2 , and that cooperativity is mediated by direct protein-protein interactions between cI bound at each operator [59]. Mutant forms of cI that are cooperativity deficient (i.e., not able to bind cooperatively to the promoter) have been designed [60]. In the inset in Figure 3B, we compare the normalized variance of the mRNA distribution, both for wild-type cI repressor, and for a cooperativity deficient mutant such as Y210H [60]. The cooperative repressor is predicted to have significantly larger promoter noise than the cooperativity deficient mutant.
Simultaneous binding of one repressor to two operators: DNA looping. Repression may also be enhanced by the presence of distant operators, which stabilize the repressed state by allowing certain repressors to simultaneously bind to both distant and proximal operators, forming a DNA loop [61,62]. The P lac promoter is a prominent example of this architecture. The kinetic mechanism of repression characterizing this promoter architecture is presented in Figure 4A. The repressor only prevents transcription when it is bound to the main operator Om, but not when it is only bound to the auxiliary operator Oa. DNA loop formation is characterized by a kinetic rate k loop~k 0 on J ½ where J ½ , the looping J-factor, can be thought of as the local concentration of repressor in the vicinity of one operator when the repressor is bound to the other operator [33,34]. The rate of dissociation of the operator-repressor complex in the looped conformation is given by k unloop~c k off R . The parameters J ½ and c have both been measured in vitro for the particular case of the Lac repressor [42,63], and also estimated from in vivo data [33,64]. Thê K KandR R matrices for this architecture are defined in Table S1 in Text S1. We use these matrices to compute the mean and the noise strength, according to equations (5)(6)(7)(8)(9)(10)(11)(12) resulting in the fifth entry of Table S2 in Text S1.
We first examine how the presence of the auxiliary operator affects the level of cell-to-cell variability in mRNA expression. In Figure 4B we compare the Fano factor in the absence of the auxiliary operator with the Fano factor in the presence of the auxiliary operator, which is assumed to be of the same strength as the main operator. We use parameters in Table 1, and we first assume that the dissociation rate of the operator-repressor complex in the looped state is the same as the dissociation rate in the unlooped state, so c~1 and k unloop~k off R . This assumption is supported by single-molecule experiments in which the two operators are on the same side of the DNA double-helix, separated by multiples of the helical period of DNA [42,63]. Under these conditions we find that the presence of an auxiliary operator results in a larger Fano factor, in spite of the fact that the auxiliary operator Oa does not stabilize the binding of the repressor to the main operator Om. Interestingly, we find that the Fano factor is maximal at intermediate concentrations of repressor for which only one repressor is bound to the promoter, making the simultaneous occupancy of the auxiliary and main operators mediated by DNA looping possible. In contrast, the Fano factor is identical to that of the simple repression case if the concentration of repressor is so large that it saturates both operators and looping never occurs. It had been previously hypothesized that DNA looping might be a means to reduce noise in gene expression, due to rapid re-association kinetics between Om and a repressor that is still bound to Oa, which may cause short and frequent bursts of transcription [64,65]. Here, by applying a simple stochastic model of gene regulation, we show that the presence of the auxiliary operator does not, by itself, decrease cell-to-cell variability. On the contrary, it is expected to increase it. The reason for this increase is that the rate of dissociation from the main operator is not made faster by DNA looping; instead the presence of the auxiliary operator causes the repressor to rapidly rebind the main operator, extending the effective period of time when the promoter is repressed.
Indeed, we find that only if the dissociation rate for a repressor in the looped state is faster than in the unlooped state, the presence of the auxiliary operator might reduce the cell-to-cell variability. To illustrate this limit, we have assumed a value of c = 100, so that k unloop~1 00 k off R , and find that the Fano factor goes down, below the expectation for the simple repression architecture. A modest increase in the dissociation rate in the looped conformation has been reported in recent single-molecule experiments for promoter architectures in which the two operators are out of phase (located on different faces of the DNA) [42]. In order to verify the general validity of these conclusions, we have randomly chosen 10,000 different sets of kinetic parameters and compared the Fano factor for an architecture with an auxiliary operator and an architecture without the auxiliary operator (simple repressor). In this analysis the operator strength, rate of transcription, rate of DNA loop formation and mean mRNA are randomly sampled over up to 4 orders of magnitude. The results are shown in Figure 4 in Text S1. In the limit where dissociation of the repressor from the operator is not affected by DNA looping c = 1, we find that the presence of the auxiliary operator leads to an increase in noise ( Figure 4A in Text S1). In contrast, we find that when this parameter c is allowed to be larger than 1, the presence of the auxiliary operator reduces cellto-cell variability in many instances ( Figure 4B in Text S1).
An example of this type of architecture is a simplified variant of the P lacUV5 promoter, which consists of one main operator and one auxiliary operator upstream from the promoter. The kinetic mechanism of repression is believed to be identical to the one depicted in Figure 4A [23,42,63,64]. We can use the stochastic model of gene regulation described in the theory section to make precise predictions that will test this kinetic model of gene regulation by DNA looping. We find that the kinetic model predicts that, if we move the center of the auxiliary operator further upstream from its wild-type location, in increments of distance given by the helical period of the DNA, such that both operators stay in phase, the fold-change in noise should behave as represented in Figure 4C. In order to model the effect of DNA looping, we assume that the dependence of the rate of DNA looping on the inter-operator distance D (in units of base-pairs) is given by [33], k loop~k where u~140:6, v~2:52, w~0:0014, z~19:9 [33], and we assume the same concentration of repressors (and therefore the same value for k on R ) for all of the different loop lengths. Note that in Figure 4C, the Fano factor is not plotted as a function of the mean, but as a function of the inter-operator distance D. That is, we keep the number of repressors constant, and instead we alter the distance between the two operators. In particular, as the operator distance is changed, both the mean and the variance will change, and therefore a direct comparison between Figures 4C and 4B cannot be made. If we had plotted the Fano factor as a function of the mean (as we do in Figure 4B) we would have seen that, for the same mean, the Fano factor for looping is always larger than for a simple repression motif, consistent with Figure 4B.

Simple activation
Transcriptional activators bind to specific sites at the promoter from which they increase the rate of transcription initiation by either direct contact with one or more RNAP subunits or indirectly by modifying the conformation of DNA around the promoter [57]. The simplest example of an activating promoter architecture consists of a single binding site for an activator in the vicinity of the RNAP binding site. When the activator is not bound, transcription occurs at a low basal rate. When the activator is bound, transcription occurs at a higher, activated rate. Stochastic association and dissociation of the activator causes fluctuations in transcription rate which in turn cause fluctuations in mRNA copy number.
This simple activation architecture is illustrated in Figure 1A. TheK KandR R matrices for this architecture are given in Table S1 in Text S1. Solving equations (5)(6)(7)(8) for this particular case, we find that the mean mRNA per cell for this simple mechanism takes the form: The mean mRNA can be changed by adjusting the intracellular concentration of the activator. The rate at which one of the activators binds to the promoter is proportional to the activator concentration: k on A~k 0 on N A ½ . Following the same argument as we used in the simple repression case, the equilibrium dissociation constant for the activator-promoter interaction is given by K OAk k 0 on . Finally, it is convenient to define the enhancement factor: the ratio between the rate of transcription in the active and the basal states f~r 2 =r 1 . The mean mRNA can be written in terms of these parameters as: The Fano factor can be computed using equations (5)(6)(7)(8)(9)(10)(11)(12) and it is shown as entry 2 of Table S2 in Text S1. We can rewrite the equation appearing in Table S2 in Text S1 by writing k on A as a function of the mean: With these equations in hand, we explore how operator strength affects noise in gene expression in the case of activation. Stronger operators bind to the activator more tightly than weak operators, leading to longer residence times of the promoter in the active state.
In Figure 5A we plot the Fano factor as a function of the foldchange in mean expression for a strong operator as well as a 10 times weaker operator. We have used the parameters in Table 1. Just as we saw for the simple repression architecture, it is also true for the simple activation architecture that stronger operators cause larger levels of noise for activators than weaker operators.
To get a sense of the differences between these two standard regulatory mechanisms, we compare simple repression with simple activation. In Figure 5B, we plot the Fano factor as a function of the mean for a repressor and an activator with identical dissociation rates. We assume that the promoter switches between a transcription rate r~0 in its inactive state (which happens when the repressor is bound in the simple repression case, or the activator is not bound in the simple activation case), and a rate equal to r~0:33 s {1 (see Table 1) in the active state (repressor not bound in the simple repression case, activator bound in the simple activation case). As shown in Figure 5B, at low expression levels the simple activation is considerably (.20 times) noisier than the simple repression promoter. At high expression levels both architectures yield very similar noise levels, with the simple repression architecture being slightly noisier. A low level of gene expression may be achieved either by low concentrations of an activator, or by high concentrations of a repressor. Low concentrations of an activator will lead to rare activation events. High concentrations of a repressor will lead to frequent but shortlasting windows of time for which the promoter is available for transcription. As a result, and as we illustrate in Figure 5C, the activation mechanism leads to bursty mRNA expression whereas the repressor leads to Poissonian mRNA production. This result suggests that in order to maintain a homogeneously low expression level, a repressive strategy in which a high concentration of repressor ensures low expression levels may be more adequate than a low activation strategy. We confirmed that this statement is true for other parameter sets in addition to the particular choice used above. We randomly sampled the rates of activator and repressor dissociation, as well as the rates of basal and maximum transcription. As shown in Figure 3 in Text S1, the statement that the simple activation architecture is noisier than the simple repression architecture at low expression (less than 10 mRNA/cell) levels is valid for a wide range of parameter values, with over 99% of the conditions sampled leading to this conclusion.
An example of simple activation is the wild-type P lac promoter, which is activated by CRP when complexed with cyclic AMP (cAMP). CRP is a ubiquitous transcription factor, and is involved in the regulation of dozens of promoters, which contain CRP binding sites of different strengths [66]. In the inset of Figure 5A we include CRP as an example of simple activation, and make predictions for how changing the wild-type CRP binding site in the P lac promoter by the CRP binding site of the P gal promoter (which is ,8 times weaker [67]), should affect the Fano factor. As expected from our analysis of this class of promoters, the noise goes down.

Dual activation: Independent and cooperative activation
Dual activation architectures have two operator binding sites. Simultaneous binding of two activators to the two operators may lead to a larger promoter activity in different ways. For instance, in some promoters each of the activators may independently contact the polymerase, recruiting it to the promoter. As a result, the probability to find RNAP bound at the promoter increases and so does the rate of transcription [33,68]. In other instances, there is no increase in enhancement factor when the two activators are bound. However, the first activator recruits the second one through protein-protein or protein-DNA interactions, stabilizing the active state and increasing the fraction of time that the promoter is active [59]. These two modes are not mutually exclusive, and some promoters exhibit a combination of both mechanisms [69].
We first investigate the effect of dual activation in the limit where binding of the two transcription factors is not cooperative. Assuming that activators bound at the two operators independently recruit the polymerase, we compare this architecture with the simple activation architecture. The mechanism of activation is depicted in Figure 6A, and matricesK K andR R are presented in Table S1 in Text S1. For simplicity, we assume that both operators have the same strength, and both have the same enhancement factor f~r 2 =r 1~r3 =r 1 . When the two activators are bound, the total enhancement factor is given by the product of the individual enhancement factors, which in this case is f |f~r 4 =r 1 [33]. All of the other relevant kinetic parameters are given in Table 1. The Fano factor is plotted in Figure 6B. We find that compared to the single operator architecture, the second operator increases the level of variability, even when binding to the operators is noncooperative.
We then ask whether this is also true when the binding of activators is cooperative. We assume a small cooperativity factor V~0:1. Just as we found for repressors, cooperative binding of activators generates larger cell-to-cell variability than independent binding, which in turn generates larger cell-to-cell variability than simple activation. This is illustrated in the stochastic simulation in Figure 6C. As expected the dual activation architectures are noisier than the simple activation, characterized by rare but long lived activation events that lead to large fluctuations in mRNA levels. In contrast, the simple activation architecture leads to more frequent but less intense activation events.
Together with the results from the dual repressor mechanism, these results indicate that multiplicity in operator number may introduce significant intrinsic noise in gene expression. Multiple repeats of operators commonly appear in eukaryotic promoters [1,70,71], but are often found in prokaryotic promoters as well [59,68,72]. It is interesting to note that this prediction of the model is in qualitative agreement with the findings by Raj et al. [2] who report an increase in cell-to-cell variability in mRNA when the number of activator binding sites was changed from one to seven.
An example of cooperative activation is the lysogenic phage-l P RM promoter [59]. This promoter contains three operators (O R1 , O R2 and O R3 ) for the cI protein, which acts as an activator. When O R2 is occupied, cI activates transcription. O R1 has no direct effect on the transcription rate, but it helps recruit cI to O R2 , since cI binds cooperatively to the two operators. Finally, O R3 binds cI very weakly, but when it is occupied, P RM becomes repressed. There are variants of this promoter [50] that harbor mutations in O R3 that make it unable to bind cI. In Figure 6D, we include one of these variants, r1-P RM [51] as an example of dual activation, and we present a theoretical prediction for the promoter noise as a function of the mean mRNA. We examine the role of cooperativity by comparing the wild-type cI, with a cooperativity deficient mutant. We find that the cooperative activator causes substantially larger cell-to-cell variability than the mutant, emphasizing our expectation that cooperativity may cause substantial noise in gene expression in bacterial promoters such as P RM .

Discussion
The DNA sequence of a promoter encodes the binding sites for transcriptional regulators. In turn, the collection of these regulatory sites, known as the architecture of the promoter, determines the mechanism of gene regulation. The mechanism of gene regulation, determines the transcriptional response of a promoter to a specific input, in the form of the concentration of one or more transcription factors or inducer molecules. In recent years we have witnessed an increasing call for quantitative models of gene regulation that can serve as a conceptual framework for reflecting on the explosion of recent quantitative data, testing hypotheses, and proposing new rounds of experiments [34,73,74]. Much of this data has come from bulk transcription experiments with large numbers of cells, in which the average transcriptional response from a population of cells (typically in the form of the level of expression of a reporter protein) was measured as a function of the concentration of a transcription factor or inducer molecule [50,75]. Thermodynamic models [34,41,53] of gene regulation are a general framework for modeling gene regulation and dealing with this kind of bulk transcriptional regulation experiments. This class of models has proven to be very successful at predicting gene expression patterns from the promoter architecture encoded in the DNA sequence [49,[73][74][75][76][77]. However, a new generation of experiments now provides information about gene expression at the level of single-cells, with single-molecule resolution [2,4,5,6,9,10,23,31,47,51]. These experiments provide much richer information than just how the mean expression changes as a function of an input signal: they tell us how that response is spread among the population of cells, distinguishing homogeneous responses, in which all cells express the same amount of proteins or mRNA for the same input, from heterogeneous responses in which some cells achieve very high expression levels while others maintain low expression. Thermodynamic models are unable to explain the single-cell statistics of gene expression, and therefore are an incomplete framework for modeling gene regulation at the single-cell level.
A class of stochastic kinetic models have been formulated that make it possible to calculate either the probability distribution of The Fano factor is plotted as a function of the fold-change gene expression (blue line). In red, we show the effect of reducing operator strength (i.e., reducing the lifetime of the operator-activator complex) by a factor of 10. Just as we observed with single repression, weak activator binding operators generate less promoter noise than strong activating operators. The parameters used are shown in Table 1 with the exception of r 1~0 :33 s {1 f , where f is the enhancement factor. Inset: Prediction for the activation of the P lac promoter. The fold-change in noise is plotted as a function of the fold-change in mean mRNA expression for both the wild-type P lac (CRP dissociation time = 8 min), represented by a blue line, and a P lac promoter variant where the lac CRP binding site has been replaced by the weaker gal CRP binding site (dissociation time = 1 min). The enhancement factor was set to f~50 [33]. These parameters are taken from [67] and [33]. The remaining parameters are taken from Table 1. (B) Fano factor as a function of SmRNAT= SmRNAT max for a repressor (black) and an activator (red) with the same transcription factor affinity. The transcription rate in the absence of activator is assumed to be zero. The transcription rate in the fully activated case is equal to the transcription rate of the repression construct in the absence of repressor and is r~0:33 s {1 as specified by Table 1. For low expression levels SmT=SmT max v0:5 simple activation is considerably noisier than simple repression. (C) The results of a stochastic simulation for the simple activation and simple repression architectures. We assume identical dissociation rates for the activator and repressor, and identical rates of transcription in their respective active states. As shown in (B), low concentrations of an activator result in few, but very productive transcription events, whereas high concentrations of a repressor lead to the frequent but short lived excursions into the active state. doi:10.1371/journal.pcbi.1001100.g005 mRNA or proteins per cell or its moments, for simple models of gene regulation involving one active and one inactive promoter state [36,37,45,78]. Recently, we have extended that formalism to account for any number of promoter states [30], allowing us to model any promoter architecture within the same mathematical framework. Armed with this model, we can now ask how promoter architecture affects not only the response function, but also how that response is distributed among different cells.
In this paper we have explored the feasibility of this stochastic analog of thermodynamic models as a general framework to understand gene regulation at the single-cell level. Using this approach we have examined a series of common promoter architectures of increasing complexity, and established how they affect the level of cell-to-cell variability of the number of mRNA molecules, and proteins, in steady state. We have found that, given the known kinetic rates of transcription factor association and dissociation from operators, the level of variability in gene expression for many well studied bacterial promoters is expected to be larger than the simple Poissonian expectation, particularly for mRNA and short-lived proteins. We have investigated how the level of variability generated by a simple promoter consisting of one single operator differs from more complex promoters containing more than one operator, and found that the presence of multiple operators increases the level of cell-to-cell variability even in the absence of cooperative binding. Cooperative binding makes the effect of operator multiplicity even larger. We also found that operator strength is one of the major determinants of cell-to-cell variability. Strong operators cause larger levels of cellto-cell variability than weak operators. We have also examined the case where one single repressor may bind simultaneously to two operators by looping the DNA in between. We have found that the stability of the DNA loop is the key parameter in determining whether DNA looping increases or decreases the level of variability, suggesting a potential role of DNA mechanics in regulating cell-to-cell variability.
We have examined the difference between activators and repressors, and found that repressors tend to generate less cell-tocell variability than activators at low expression levels, whereas at  Table 1 and r 1~0 :33 s {1 f , r 2~f |r 1 , r 3~f |r 1 , and r 4~f 2 |r 1 ; f is the enhancement factor. (C) A stochastic simulation shows the effect of independent and cooperative binding in creating a sustained state of high promoter activity, resulting in high levels of mRNA in the active state and large cell-to-cell variability. (D) Prediction for the r1-P RM promoter (a P RM promoter variant that does not exhibit O R3 mediated repression [51]). This promoter is activated by cI, which binds cooperatively to O R1 and O R2 . The prediction is shown for wild-type cI (V~0:013) and for a cooperativity deficient mutant (Y210H, V~1). Parameters are taken from [33,43,58,97]. The lifetime of O R1 -cI complex is 4 min. Lifetime of O R2 -cI complex is 9.5 s. doi:10.1371/journal.pcbi.1001100.g006 high expression levels repressors and activators generate similar levels of cell-to-cell variability. We conclude that induction of gene expression by increasing the concentration of an activator leads to a more heterogeneous response at low and moderate expression levels than induction of gene expression by degradation, sequestration or dilution of a repressor. In addition, we have used this model to make quantitative predictions for a few well characterized bacterial promoters, connecting the kinetic mechanism of gene regulation that we believe applies for these promoters in vivo with single-cell gene expression data. Direct comparison between the model and experimental data offers an opportunity to validate these kinetic mechanisms of gene regulation.

Intrinsic and extrinsic noise
There are two different classes of sources of cell-to-cell variability in gene expression. The first class has its origins in the intrinsically stochastic nature of the chemical reactions leading to the production and degradation of mRNAs and proteins, including the binding and unbinding of transcription factors, transcription initiation, mRNA degradation, translation, and protein degradation. The noise coming from these sources is known as intrinsic noise [79]. A different source of variability originates in cell-to-cell differences in cell size, metabolic state, copy number of transcription factors, RNA polymerases, ribosomes, nucleotides, etc. This second kind of noise is termed extrinsic noise [79]. The contributions from intrinsic and extrinsic sources can be separated experimentally, and the total noise can be written as the sum of intrinsic and extrinsic components [3]. In this paper we focus exclusively on intrinsic noise, and the emphasis is on bacterial promoters. This double focus requires us to discuss to what extent intrinsic noise is relevant in bacteria.
The experimental evidence gathered so far indicates that intrinsic noise is the dominant source of cell-to-cell variability in bacteria of the mRNA copy number. In a recent single-molecule study, transcription was monitored in real time for two different E. coli promoters, P RM and P lac/ara [4]. The authors measured the rates of mRNA synthesis and dilution, as well as the rates of promoter activation and inactivation in single cells. The intrinsic noise contribution was calculated from all of these rates. It was found to be responsible for the majority of the total cell-to-cell variability, accounting for over 75% of the total variance. Another recent experiment in B. subtilis [7] found that mRNA expressed from the ComK promoter is also dominated by intrinsic noise. Furthermore, this study indicated that intrinsic mRNA noise is responsible for activation of a phenotypic switch that drives a fraction of the cells to competence for the uptake of DNA [7]. A third recent report investigated the activation of the genetic switch in E. coli, which drives the entrance of a fraction of cells into a lactose metabolizing phenotype [23]. The authors of the study found evidence that stochastic binding and unbinding of the Lac repressor to the main operator was responsible for the observed cell-to-cell variability in gene expression and, consequently the choice of phenotype. Furthermore, the authors discovered that the deletion of an auxiliary operator that permits transcriptional repression by DNA looping, leads to a strong increase in the level of cell to cell variability in the expression of the lactose genes, indicating that promoter architecture plays a big role in determining the level of noise and variability in this system. Taken all together, these experiments suggest that intrinsic mRNA noise is dominant and may have important consequences for cell fate determination. In addition, at least in one case, promoter architecture has been shown to be of considerable importance.
At the protein level, the contribution of extrinsic and intrinsic noise to the total cell-to-cell variability has also been determined experimentally for a variety of promoters and different kinds of bacteria. The first reports examined intrinsic and extrinsic protein noise in E. coli and found that extrinsic noise was the dominant source of cell-to-cell variability in protein expressed from a variant of the P L promoter in a variety of different strains [3]. However, the intrinsic component was non-negligible and for some strains, dominant [3]. A second team of researchers examined a different set of E. coli promoters involved in the biosynthetic pathway of lysine [80]. The authors found that the intrinsic noise contribution was significant for some promoters (i.e. lysA), but not for others. In a third study the total protein noise was measured for a Lac repressor-controlled promoter in B. subtilis, and it was reported that the data could be well explained by a model consisting only of intrinsic noise [8]. The authors found that the rates of transcription and translation could be determined by directly comparing the total cell-to-cell variability to the predictions of a simple stochastic model that considered only intrinsic sources of noise. They also found that the model had predictive power, and that mutations that enhanced the rate of translation or transcription produced expected effects in the total noise.
In summary, all studies that have measured mRNA noise in bacteria so far report that intrinsic noise contributes substantially to the total cell-to-cell variability. This is further supported by observations that most of the mRNA variability comes from intrinsic sources in yeast [31] and mammalian cells [1]. The issue is less clear for protein noise. Some reports indicate that it is mostly extrinsic [3], but others suggest that intrinsic noise may also be important [8,23,80]. It seems likely that the relative importance of intrinsic and extrinsic noise depends on the context, and that for some promoters and genes extrinsic noise will be larger, whereas for others the intrinsic component may dominate. In any case, it is clear that both contributions are important, and both need to be understood.

Comparison with experimental results
The aim of this paper is to formulate a set of predictions that reflect the class of kinetic models of gene regulation in bacteria that one routinely finds in the literature [42,64,[81][82][83][84]. Our analysis indicate that if these models are correct, and if the kinetic and thermodynamic parameters that have been measured over the years are also reasonably close to their real values in live cells [85], the effect of promoter architecture in cell-to-cell variability in bacteria should be rather large and easily observable. In this sense, our intention is more to motivate new experiments than to explain or fit any currently available data. We only know of one published report in which the effect of perturbing the architecture of a bacterial promoter on the cell-to-cell variability in gene expression has been determined [23]. Given that there are several examples of promoters in bacteria for which a molecular kinetic mechanism of gene regulation has been formulated [42,64,[81][82][83][84]86], we hope that the computational analysis in this paper may serve as an encouragement for researchers to do for bacteria the same kind of experiments that have been already performed in eukaryotes [1,11,15,17,31]. Indeed, several different studies have examined the effect of promoter architectural elements in cell-to-cell variability in protein and mRNA in eukaryotic cells. Although our efforts in this paper have focused on bacterial promoters rather than eukaryotic promoters, it is worthwhile to discuss the findings of these studies and compare them (if only qualitatively) with the predictions made in this paper.
Two recent studies measured intrinsic mRNA noise in yeast [31] and mammalian cells [1]. Both papers concluded that stochastic promoter activation and inactivation was the leading source of intrinsic noise. While stochastic chromatin remodeling is suspected to be the origin of those activation events, neither one of these studies was conclusive about the precise molecular mechanism responsible for promoter activation. However, both studies found that promoter architecture had an important role and strongly affected the level of total mRNA noise. In both studies, the authors found that when the number of binding sites for a transcriptional activator was raised from one to seven, the normalized variance increased several-fold. This qualitative behavior is in agreement with our prediction that dual activation causes larger intrinsic mRNA noise than simple activation. It is possible that this agreement is coincidental, since the actual mechanism of gene regulation at these promoters could be much more complicated than the simple description of gene activation at a bacterial promoter adopted here.
Other studies [11,15,17] have measured the total protein noise from variants of the GAL1 promoter in yeast, and found that their data could be well explained by a model that considered only intrinsic noise sources. These studies also concluded that the main sources of intrinsic noise were stochastic activation and inactivation of the promoter due to chromatin remodeling. However, it was also found that the stable formation of pre-initiation complex at the TATA box and the stochastic binding and unbinding of transcriptional repressors contributed to the total noise [11,15,17]. The authors of these studies found that for point mutations in the TATA box of the GAL1 promoter in yeast, which made the box weaker, the level of cell to cell variability went down significantly. This is also in good agreement with our prediction that the stronger the binding site of a transcriptional activator, the larger the intrinsic noise should be. However, since this study measured the total noise strength, and did not isolate the intrinsic noise, the observed decrease in noise strength as a result of making the TATA box weaker may have other origins. These experiments were conducted under induction conditions that minimize repression by nucleosomes and activation by chromatin remodeling. A more recent report by the same lab [11] found that the copy number and location of a transcriptional repressor binding site greatly affects the total protein noise. The authors found that when they increased the number of repressor binding sites, the noise went up. This is also in qualitative agreement with our prediction that operator number positively correlates with intrinsic noise in the case of dual repression. However, the same caveat applies here as in the previous case studies, which is that only the total noise was measured. Although the authors of this study attributed all of the noise to intrinsic sources, it is still possible that extrinsic noise was responsible for the observed dependence of noise strength on operator number.
Finally, it is worth going back to bacteria, and discussing the only study that has yet examined the effect of a promoterarchitecture motif on cell-to-cell variability in gene expression. In this paper, the authors investigated the effect of DNA looping on the total cell-to-cell variability for the P lacUV5 promoter in E. coli [23]. Using a novel single-protein counting technique, Choi and co-workers measured protein distributions for promoters whose auxiliary operator had been deleted (leaving them with a simple repression architecture), and compared them to promoters with the auxiliary operator O3 present, which allows for DNA looping. They report a reduction in protein noise due to the presence of O3, which according to our analysis, may indicate that the dissociation of the repressor from the looped state is faster than the normal dissociation rate. The authors attributed this loopingdependent decrease in noise to intrinsic origins, related to the different kinetics of repressor binding and rebinding to the main operator in the presence of the auxiliary operator, and in its absence. However, their measurements also reflect the total noise, and not only the intrinsic part, so the explanation may lie elsewhere. These results emphasize the need for more experiments in which the intrinsic noise is isolated and measured directly.
More recently, several impressive experimental studies have measured the noise in mRNA in bacteria for a host of different promoters ( [87], and Ido Golding, private communication). In both of these cases, simplified low-dimensional models which do not consider the details of the promoter architecture have been exploited to provide a theoretical framework for thinking about the data. Our own studies indicate that the differences between a generic two-state model and specific models that attempt to capture the details of a given architecture are sometimes subtle and that the acid test of ideas like those presented in this paper can only come from experiments which systematically tune parameters, such as the repressor concentration, for a given transcriptional architecture.

Future directions
Some recent theoretical work has analyzed the effect of cooperative binding of activators in the context of particular examples of eukaryotic promoters [88,89]. The main focus of this study is bacterial promoters. The simplicity of the microscopic mechanisms of transcriptional regulation for bacterial promoters makes them a better starting point for a systematic study like the one we propose. However, many examples of eukaryotic promoters have been found whose architecture affects the cellto-cell variability [1,11,17,31,32]. Although the molecular mechanisms of gene regulation in these promoters are much more complex, with many intervening global and specific regulators [90], the stochastic model employed in this paper can be applied to any number of promoter states, and thus can be applied to these more complex promoters. Recent experimental work is starting to reveal the dynamics of nucleosomes and transcription factors with single-molecule sensitivity [91,92], allowing the formulation of quantitative kinetic and thermodynamic mechanistic models of transcriptional regulation at the molecular level [73,77]. The framework for analyzing gene expression at the single-cell level developed in this paper will be helpful to investigate the kinetic mechanisms of gene regulation in eukaryotic promoters, as the experimental studies switch from ensemble, to single-cell.

Shortcomings of the approach
Although the model of transcriptional regulation used in this paper is standard in the field, it is important to remark that it is a very simplified model of what really happens during transcription initiation. There are many ways in which this kind of model can fail to describe real situations. For instance, mRNA degradation requires the action of RNases. These may become saturated if the global transcriptional activity is very large, and degradation the becomes non-linear [55]. Transcription initiation and elongation are assumed to be jointly captured in a single constant rate of mRNA synthesis for each promoter state. This is an oversimplification also. When considered explicitly, and in certain parameter ranges, the kinetics of RNAP-promoter interaction may cause noticeable effects in the overall variability [46]. Similarly, as pointed out elsewhere [93,94,95], translational pausing, backtracking or road-blocking may also cause significant deviations in mRNA variability from the predictions of the model used in this paper. How serious these deviations are depends on the specifics of each promoter-gene system. The model explored in this paper also assumes that the cell is a well-mixed environment. Deviations from that approximation can significantly affect cell-to-cell variability [56,96]. Another simplification refers to cell growth and division, which are not treated explicitly by the model used in this paper: cell division and DNA replication cause doubling of gene and promoter copy number every cell cycle, as well as binomial partitioning of mRNAs between mother and daughter cells [3]. In eukaryotes, mRNA often needs to be further processed by the splicing apparatus before it becomes transcriptionally active. It also needs to be exported out of the nucleus, where it can be translated by ribosomes.
To study the effect of transcription factor dynamics on mRNA noise we assume that the unregulated promoter produces mRNA in a Poisson manner, at a constant rate. This assumption can turn out to be wrong if there is another process, independent of transcription factors, that independently turns the promoter on and off. In eukaryotes examples of such processes are nucleosome positioning and chromatin remodeling, while in prokaryotes analogous processes are not as established, but could include the action of non-specifically bound nucleoid proteins such as HU and HNS, or DNA supercoiling. Experiments that measure cell-to-cell distributions of mRNA copy number in the absence of transcription factors (say without Lac repressor for the lac operon case) can settle this question. In case the Fano factor for this distribution is not one (as expected for a Poisson distribution) this can signal a possible transcription factor-independent source of variability. The stochastic models studied here can be extended to account for this situation. For example, the promoter can be made to switch between an on and an off state, where the transcription factors are allowed to interact with promoter DNA only while it is in the on state. In this case the mRNA fluctuations produced by an unregulated promoter will not be Poissonian. One can still investigate the affect of transcription factors by measuring how they change the nature of mRNA fluctuations from this new baseline. Comparison of this extended model with single-cell transcription experiments would then have the exciting potential for uncovering novel modes of transctriptional regulation in prokaryotes.
For the purpose of isolating the effect of individual promoter architectural elements on cell-to-cell variability in gene expression, we have artificially changed the value of one of those parameters, while keeping the other parameters constant. For instance, we have investigated the effect of altering the strength of an operator on the total cell-to-cell variability. In order to do this, we ask how changes in the dissociation rate of the transcription factor alter the cell-to-cell variability, given that all other rates (say the rate of transcription, or mRNA degradation) remain constant. This assumption is not necessarily always correct, since very often the operator sequence overlaps the promoter, and therefore changes in the sequence that alter operator strength also affect the sequence from which RNAP initiates transcription, which can potentially affect the overall rates of transcription. As is usually the case, biology presents us with a great diversity of forms, shapes and functions, and promoters are no exception. One needs to examine each promoter independently on the basis of the assumptions made in this paper, as many of these assumptions may apply for some promoters, but not for others.
For the same reason of isolating the effect of promoter architecture and cis-transcriptional regulation on cell-to-cell variability in gene expression, when we compare different architectures we make the simplifying assumption that they are transcribing the same gene, and therefore that the mRNA transcript has the same degradation rate. Care must be taken to take this into account when promoters transcribing different genes are investigated, since the mRNA degradation rate has a large effect on the level of cell-to-cell variability.
We have also assumed that when transcription factors dissociate from the operator, they dissociate into an averaged out, wellmixed, mean-field concentration of transcription factors inside the cell. The possibility of transcription factors being recaptured by the same or another operator in the promoter right after they fall off the operator is not captured by the class of models considered here. Recent in vivo experiments suggest that this scenario may be important in yeast promoters containing arrays of operators [31].
In spite of all of the simplifications inherent in the class of models analyzed in this paper, we believe they are an adequate jumping off point for developing an intuition about how promoter architecture contributes to variability in gene expression. Our approach is to take a highly simplified model of stochastic gene expression, based on a kinetic model for the processes of the central dogma of molecular biology, and add promoter dynamics explicitly to see how different architectural features affect variability. This allows us to isolate the effect of promoter dynamics, and develop an intuitive understanding of how they affect the statistics of gene expression.
It must be emphasized, however, that the predictions made by the model may be wrong if any of the complications mentioned above are significant. This is not necessarily a bad outcome. If the comparison between experimental data and the predictions made by the theory for any particular system reveals inconsistencies, then the model will need to be refined and new experiments are required to identify which of the sources of variability that are not accounted for by the model are in play. In other words, experiments that test the quantitative predictions outlined stand a chance of gaining new insights about the physical mechanisms that underlie prokaryotic transcriptional regulation.

Supporting Information
Text S1 Mathematical derivations and supplementary information. A derivation of all equations in the text is presented, together with its corresponding tables and figures. Found at: doi:10.1371/journal.pcbi.1001100.s001 (2.04 MB DOC)