Protein Distributions from a Stochastic Model of the lac Operon of E. coli with DNA Looping: Analytical solution and comparison with experiments

Although noisy gene expression is widely accepted, its mechanisms are subjects of debate, stimulated largely by single-molecule experiments. This work is concerned with one such study, in which Choi et al., 2008, obtained real-time data and distributions of Lac permease in E. coli. They observed small and large protein bursts in strains with and without auxiliary operators. They also estimated the size and frequency of these bursts, but these were based on a stochastic model of a constitutive promoter. Here, we formulate and solve a stochastic model accounting for the existence of auxiliary operators and DNA loops. We find that DNA loop formation is so fast that small bursts are averaged out, making it impossible to extract their size and frequency from the data. In contrast, we can extract not only the size and frequency of the large bursts, but also the fraction of proteins derived from them. Finally, the proteins follow not the negative binomial distribution, but a mixture of two distributions, which reflect the existence of proteins derived from small and large bursts.


Introduction
Data from many independent experiments show that the abundance of any given protein varies among individual cells of isogenic populations growing under identical conditions [1][2][3]. Early experiments with fluorescent reporters showed that such non-uniformity in protein abundance was due to the inherent stochasticity of gene expression (intrinsic noise) and various forms of cell-to-cell variation (extrinsic noise) [4,5]. The subsequent development of single-molecule techniques has led to deeper insights into the molecular mechanisms generating the noise [6,7]. By measuring the number of mRNAs in single cells, Golding et al. showed that transcription was too bursty to be modeled as a Poisson process [8]. Cai et al. [9] and Yu et al. [10] developed two different methods for measuring the number of proteins in single cells. The real-time data of both studies showed that protein synthesis was bursty, and the burst size was exponentially distributed. Under this condition, the steady state protein distribution follows the Gamma distribution, p n~n a{1 e {n=b =b a C(a), where a and b denote the mean burst frequency and burst size [11]. Cai et al. and Yu et al. showed that the Gamma distribution could fit their steady state data, and the values of the mean burst frequency and size derived from the steady state data agreed well with those obtained from real-time measurements.
Armed with these results, Choi et al. [12] attacked a longstanding problem. When non-induced cells of E. coli are exposed to small concentrations of the gratuitous inducer TMG, the lac operon is induced by stochastic switching of individual cells from the non-induced to the induced state [13]. Choi et al. sought the molecular mechanism of this stochastic switching. To this end, they first quantified the minimum number of LacY molecules required to switch a cell to the induced state, and found this threshold to be 375 molecules. They then suggested a molecular mechanism capable of yielding this threshold by appealing to the known mechanisms of repression and transcription of the lac operon. Repression is mediated by the stable DNA loops formed when the Lac repressor is simultaneously bound to the main and auxiliary operators (Fig. 1). Transcription can take place either due to partial dissociations, which occur when a repressor trapped in a DNA loop dissociates from the main operator, but not the auxiliary operator; or complete dissociations, which occur when the repressor dissociates completely from the DNA. Choi et al. hypothesized that since a partially dissociated repressor remains attached to the DNA, it rapidly rebinds to the main operator, thus limiting the number of transcription events. Although the evidence suggests that no more than one mRNA is made during a partial dissociation, it is conceivable that multiple transcripts are made during a partial dissociation despite its short lifetime, thus leading to a small transcriptional burst. In contrast, a completely dissociated repressor takes a relatively long time to find an operator, which results in a large transcriptional burst. These large transcriptional bursts can provide enough proteins to cross the threshold for stochastic switching.
Choi et al. tested the foregoing hypotheses as follows. The statistics of small transcriptional bursts were obtained with strain SX701, a lacY { strain that exhibits mostly small bursts. To capture the statistics of large bursts, they deleted the auxiliary operators of their lacY { cells, thus creating strain SX703 which yields only large bursts. The statistics of the small and large bursts were quantified by measuring the steady-state protein distributions for both strains at various inducer concentrations. They then concluded, based on the model of Friedman et al. [11], that if m,s 2 denote the mean and variance of a protein distribution obtained with strain SX701, then the Fano factor, F :s 2 =m, and the reciprocal of the noise, g {2 : m=s ð Þ 2 , represent the size and frequency of the small bursts. Likewise, if m m, s s 2 denote the mean and variance for SX703, then F F : s s 2 = m m, g g {2 : m m= s s ð Þ 2 represent the size and frequency of the large bursts. Analysis of the data for SX703 with this method showed that g g {2 did not change with inducer levels, but F F increased dramatically (Fig. 2a), thus confirming their hypothesis that large bursts can generate enough proteins to trigger stochastic switching. Surprisingly, analysis of the data for SX701 also yielded similar trends (Fig. 2b), but this was attributed to the distortions created by the few cells exhibiting large bursts. Indeed, if the data were filtered by removing the contribution of large bursts, g {2 and F did not change much with the inducer concentration ( Fig. 2c), leading the authors to conclude that the small burst frequency and size were independent of the inducer level.
Choi et al. also explained these results by appealing to the known states of the lac operon (Fig. 1). However, the mathematical model of Friedman et al., which forms the basis of their data analysis, does not account for these complexities -it only considers a constitutive (unregulated) promoter. Consequently, there is no strong support for the assumption that the proteins follow the Gamma distribution; F , F F represent the size of small and large bursts; and g {2 , g g {2 represent the frequency of small and large bursts. The goal of this study is to verify the validity of these assumptions by formulating a stochastic model accounting for the known states of the operon, and deriving analytical expressions for the steady state protein distribution, Fano factor, and noise.
There are stochastic models accounting for the details shown in Fig. 1 [14][15][16], but these studies do not give analytical expressions for the steady state protein distribution. The literature also contains several stochastic models of gene regulation for which analytical solutions were obtained [11,[17][18][19][20][21][22][23][24], but they do not account for the presence of multiple auxiliary operators and DNA looping. Our model fills this gap in the theoretical literature, and its analysis yields deeper insights into the experimental data. Specifically, we show that the size and frequency of small bursts cannot be extracted from the data for strain SX701 because they are averaged out. However, we can extract not only the size and frequency of the large bursts, but also their contribution to total protein synthesis, provided the data is not filtered (Fig. 2d). This result also yields tests for the consistency of the model by providing relationships between the size and frequency of large bursts in strains SX701 and SX703. Finally, we show that neither one of the two strains follow the negative binomial (or Gamma) distribution.
The paper is organized as follows. In the Analysis section, we describe the model, derive the master equation, and explain the key approximations used to obtain the steady state protein distribution. In the Results section, we perform simulations to check the validity of the analytical expression for the protein distribution, and we derive the expressions for mean and the variance of the distribution. We also show that the mean, variance, and hence, the Fano factor and the reciprocal of the noise, can be expressed in terms of the size and frequency of the transcriptional and translational bursts. In the Discussion section, the latter are compared with the assumptions of Choi et al. We also show that negative binomial distributions are obtained only if the size of the large transcriptional bursts is relatively small.

Analysis
The model scheme, shown in Figure 1, is based on the following facts enunciated by Oehler et al. [25,26]. The lac operon of E. coli contains three operators, namely the main operator O 1 , and the two auxiliary operators O 2 ,O 3 , lying downstream and upstream of O 1 . The lac operon rarely entertains more than one Lac repressor, and this single repressor R can bind to any one of the operators, thus forming the operon states, O 1 : R, O 2 : R, and O 3 : R. Since the tetrameric repressor is a "dimer of dimers,'' it has a free dimer even after it is bound to one of the operators. This free dimer can bind to one of the remaining two free operators, thus forming a DNA loop. In principle, three looped states are feasible, namely,  O 3 : R. The first two states permit full transcriptional activity. The last state can be neglected since it permits only 3-5% of the full transcriptional activity.
The model kinetics are based on the following assumptions. All cells have the same number of repressors, N, which is tantamount to neglecting extrinsic noise [4]. Since association of a cytosolic repressor to an operator is diffusion-limited, we assume that a cytosolic repressor has the same propensity, k a N, for association with each of the operators. In contrast, the propensity for dissociation of operator-bound repressor does depend on the identity of the operator, and we denote the propensity for dissociation of O i -bound repressor by k Oi . Next, we consider the kinetics of looping. The looped state O 1 : R : O 2 can be formed from either O 1 : R or O 2 : R, but both pathways have the same propensity because they are driven by the same local concentration effect [26]. Thus, we denote the propensity for formation of : R by the same symbol, k O1O3 . Finally, we let v 0 ,d 0 denote the propensities for mRNA synthesis and degradation, and v 1 ,d 1 denote the propensities for protein synthesis and dilution.

Equations
We take a master equation approach to describe the system, our state variables being the number of mRNAs, m, the number of proteins, n, and the six states of the operon shown in Figure 1. We let p s m,n denote the probability of m mRNAs and n proteins when the operon is in state s. Here, s~f when the operon is free, and s~i or s~ij when the operon is repressor-bound, where i,j are integers identifying the operator(s) to which the repressor is bound (e.g., s~1 denotes the state O 1 : R and s~12 denotes the state O 1 : R : O 2 ). Then the master equations for the kinetic scheme in Figure 1 are Our goal is to derive the steady state protein distribution corresponding to these equations. Table 1 shows the parameter values in the absence of the inducer. The parameters d 0 and d 1 reflect the experimental values measured by Yu et al. [10]. The parameter v 1 was chosen such that the the mean burst size, b:v 1 =d 0 , agreed with the measured value b~4, reported by Yu et al. The parameter v 0 was estimated by assuming that the mean burst frequency of fully induced cells, a:v 0 =d 1 , is 600. The rationale for this assumption is as follows. An uninduced cell contains, on average, 0.5 molecules of the tetrameric LacZ [9], and hence, is expected to contain 2 molecules of the monomeric LacY. Since the number of LacY and LacZ molecules increases *1200-fold in fully induced cells [25], there are 2400 LacY molecules in such cells, i.e., ab~2400, which implies that a~600. All other parameter values were estimated using the method of Vilar & Leibler [15]. They estimated all the equilibrium constants using the repression data of Oehler et al. [26]. Then, given an experimental estimate of any one parameter, they could find all other parameter values. They took that one parameter to be the dissociation rate constant, k O1 , and assigned to it the value obtained from in vitro data [27]. Based on this procedure, the association rate, k a N, was found to be 0.73 s {1 . However, recent in vivo measurement show that the association rate for a dimeric repressor is 0.014 s {1 [28]. If the dimeric and tetrameric repressor associate at the same rate, and each cell contains 10 repressors [29], the estimated value of k a N from these measurements is 0.14 s {1 . We assumed k a N~0:07 s {1 , and chose k O1 , k O2 , k O3 , k O1O2 , k O1O3 to ensure consistency with the repression data. As we show later, these parameter values yield good fits of the experimental data.

Parameter values
Since we are also concerned with protein distributions in the presence of the inducer, it is necessary to identify the parameters that change under these conditions. We assume that v 0 , d 0 , v 1 , and d 1 are independent of the inducer level. The propensities for looping, k O1O2 ,k O1O3 , are also unlikely to change in the presence of small inducer concentrations because a partially dissociated repressor has too little time to interact with the inducer: In the presence of 10 mM IPTG (considered equivalent to 100 mM TMG), the pseudo-first-order rate constant for repressor-inducer binding is 0.1 s {1 [30], which is negligible compared to the looping rate constant of 4 s {1 . Thus, the only parameters that can change with the inducer concentration are the association rate, k a N, and the dissociation rates, k Oi . Based on the analysis of their experimental protein distributions, Choi et al. concluded that the dissociation rates are independent of the inducer concentration, while the association rate decreases with the inducer concentration. We shall also assume that this is the case. This assumption holds only if the concentration of TMG is significantly below 1 mM [31,32]

Model reduction
The determination of the steady state protein distribution corresponding to eqs. (1)-(6) is facilitated by the fact that loop formation and mRNA degradation are relatively fast.
Rapid loop formation. Table 1 shows that in the absence of the inducer, k O1O2 ,k O1O3 are much greater than all other propensities, and as explained above, this persists even in the presence of low inducer concentrations. It follows that the repressor-bound states rapidly equilibrate on the fast time scale which represents the probability of m mRNAs and n proteins when the operon is repressor-bound. We then apply the quasisteady state approximation to the fast variables, p 2 m,n , p 3 m,n , p 12 m,n , p 13 m,n , and find that the probabilities of the equilibrated bound states are given by the relations which express the physical fact that after the bound states reach quasi-equilibrium, they obey the principle of detailed balance and are almost always in one of the looped states (Table 2). Moreover, the slow variables follow the equations where l: Equations (13)- (14) describe the evolution of the reduced model containing only two operon states -the free and the equilibrated bound states -between which are transitions with propensities, k 0 , k 1 , which are slow compared to the propensities for looping ( Table 2). This is highlighted in Figure 1 by enclosing the free and bound states in dashed boxes, and drawing dashed arrows with labels, k 0 and k 1 , to denote the transitions between them. The reduced model is similar to Shahrezaei & Swain's three-stage model for a regulated promoter [22], but there is an important difference. Both operon states are transcriptionally active: The transcription rates in the free and bound states are v 0 and v 0 l, respectively, where l is the probability of the O 2 : R state. Even though l%1 (Table 2), we cannot neglect the transcription from the bound state, since it captures the effect of the small transcriptional bursts, which can account, as we show later, for almost 80% of the mRNAs synthesized per cell cycle. Table 2 shows that in the absence of the inducer, k 0 %k 1 , so that the free state occurs infrequently and lasts for very short periods of time, i.e., p b m,n &1. We shall show later that this persists in the presence of the low inducer concentrations (ƒ 200 mM TMG) used by Choi et al. Hence, under the experimental conditions of interest, the conditional probabilities in (8)- (12) are essentially equal to the absolute probabilities.
Rapid mRNA degradation. The second approximation appeals to the fact that mRNA degradation is rapid compared to protein dilution, i.e., d 0 &d 1 . To apply this approximation, we follow Shahrezaei & Swain [22]. Thus, we begin by rescaling time with respect to the time scale for protein degradation. Letting  July 2014 | Volume 9 | Issue 7 | e102580 lac f b (z,z',t)~P m,n z' m z n p b m,n , to obtain the partial differential equations where u~z'{1 and v~z{1. Since c&1, we have the quasisteady state approximation, bv(1zu){u&0. The steady state protein distribution is therefore given by the equations Since we are interested in the generating function, which reduce to the second-order differential equation We solve this equation with the initial condition, f (0)~1, and revert to z as the independent variable, to obtain the following generating function for the steady state protein distribution where 2 F 1 denotes the Gaussian hypergeometric function and a,b: As expected, if l~0, (27) reduces to the generating function of the negative hypergeometric distribution [22]. In general, however, (27) is the generating function for a mixture of the negative binomial and negative hypergeometric distributions, which reflects, as we show below, the existence of two sub-populations of proteins, namely those derived from small and large transcriptional bursts.

Analytical expressions for the statistics of the protein distributions
Strain with auxiliary operators. The generating function (27) yields the following expressions for the mean, m, and variance, s 2 , of the protein distribution m~a r b, a r :a lz k 0 k 0 zk 1 , ð29Þ Since b represents the mean number of proteins synthesized per mRNA, (29) implies that a r is the mean frequency of regulated transcription. The two terms of a r also have simple physical interpretations: Since l and k 0 = k 0 zk 1 ð Þare the probabilities of the O 2 : R and free states, al and ak 0 = k 0 zk 1 ð Þrepresent the mean number of mRNAs produced per cell cycle due to small and large transcriptional bursts.
Expanding f (z) about z~0 yields the steady state protein distribution  Figure 3 shows that the protein distributions obtained from this expression agree well with those obtained by simulating the full model with the Optimized Direct Method implementation of Gillespie's Stochastic Simulation Algorithm [33] provided in the simulation package StochKit2 [34]. The protein distribution in the absence of the inducer, shown in Fig. 3a, was obtained with the parameter values in Table 1. The distributions in the presence of the inducer were obtained by decreasing the association rate, k a N, 10-fold (Fig. 3b) and 20-fold (Fig. 3c). Evidently, (31) is a good approximation to the exact solutions in all three cases. We conclude that our approximate solution is accurate down to a 20fold reduction of the association rate. Table 2 shows that in the absence of the inducer, l,k 0 =k 1 %1. These relations remain valid at the relatively low inducer levels studied by Choi et al. (ƒ200 mM TMG). Indeed, under these conditions, the operon is expressed to no more than 1% of the fully induced level [12], i.e., a r a~l z k 0 k 0 zk 1 ƒ0:01[l, k 0 k 0 zk 1 ƒ0:01[l, and (29) It is worth noting that due to rapid loop formation, small transcriptional bursts are very bursty (pulsatile). Moreover, under the weakly inducing conditions used in the experiments (ƒ 200 mM TMG), k a N is relatively large, and hence, the large transcriptional bursts are also quite bursty. It follows that under these conditions, (33)- (34) should be expressible in terms of the size and frequency of the small and large transcriptional bursts. We shall show below that this is indeed the case.
Strain without auxiliary operators. In the absence of auxiliary operators, the operon fluctuates between the free and the O 1 -bound state, and only the former allows transcription. This is identical to Shahrezaei & Swain's 3-stage model of a regulated promoter [22], and corresponds to the special case, l~0, k 0~kO1 , k 1~ka N of our model. It follows that the generating function for the steady state protein distribution is the Gaussian hypergeometric function and a a, b b: Moreover, the protein distribution is given by the expression and the mean and variance are m m~ a a r b, a a r~a k k 0 k k 0 z k k 1 , At TMG concentrations of ƒ 100 mM, which are equivalent to an IPTG concentration of ƒ 10 mM, the operon is expressed to no more than 5% of the fully induced level [35]. It follows that under the experimental conditions of interest k k 0 k k 0 z k k 1 ƒ0:05[ k k 0 k k 1 ƒ 0:05 0:95 and m m, s s 2 can be approximated by the expressions m m~ a a r b, a a r &a Expressing the statistics in terms of the burst size and frequency Strain with auxiliary operators. To express m~a r b in terms of the size and frequency of the transcriptional bursts, we begin by recalling that a r consists of two terms, al and ak 0 =k, which represent the mean frequency of transcription due to partial and complete dissociations of the repressor, respectively. Since partial dissociations occur when a repressor trapped in the O 1 O 2loop dissociates from O 1 , we define the number of the partial dissociations per cell cycle as where we have appealed to the detailed balance between the operon states O 2 : R and O 1 : R : O 2 . We also define the number of mRNAs synthesized per partial dissociation as since the time for rebinding of a partially dissociated repressor to O 1 is on the order of k {1 . It follows from these definitions that i.e., we have successfully expressed the first term of a r in terms of frequency and mRNA burst size due to partial dissociations. We now proceed to express the second term of a r in terms of the frequency and mRNA burst size due to complete dissociations. Since complete dissociations occur whenever the operon becomes repressor-free, it is natural to define the number of complete dissociations per cell cycle as We also define the number of mRNAs synthesized per complete dissociation as because the time for rebinding of a completely dissociated repressor to an operator is on the order of k {1 1 . Evidently and we conclude that Hence, (33)- (34) can be rewritten as which imply that F: where f c : is the fraction of proteins derived from complete dissociations. It follows from (53) that the total burstiness, F , is entirely due to translational and large transcriptional bursts. Moreover, the burstiness of large transcriptional bursts depends on their intrinsic burstiness, b c b, suitably weighted by f c , the fraction of proteins derived from such bursts. Importantly, f c is completely determined by k O 2 =k a N, the equilibrium constant for dissociation of the repressor from O 2 . In the absence of the inducer, this equilibrium constant is 0.25 [25,26], and hence, f c~0 :2, i.e., 20% of the proteins are derived from large transcriptional bursts. As the inducer concentration increases, f c increases because k a N decreases.
Strain without auxiliary operators. In this case, if we define the number of complete dissociations per cell cycle as a a c : and the number of mRNAs synthesized per complete dissociation as the mean frequency of regulated transcription can be rewritten as It follows that (42)-(43) can be rewritten as m m~ a a r b& a a c b b c b, which imply that F F : We are now ready to address questions concerning the physical meaning of the parameters of the distribution and their variation with inducer concentration [12].

Interpretation of the protein distribution data
Strain with auxiliary operators. Interpretation of F and g {2 derived from filtered data. Choi et al. assumed that F and g {2 derived from the filtered data (Fig. 2c) represent the size and frequency of small transcriptional bursts. In terms of our model, these assumptions have the form However, (53)-(54) imply that this F and g {2 , obtained by eliminating the contribution of the large transcriptional bursts, have a different physical meaning. Indeed, (53) implies that the Fano factor obtained from the filtered data has the form, F~1zb, which represents the size of the translational, rather than small transcriptional, bursts. Similarly, (54) implies that the reciprocal of the noise derived from the filtered data has the form, g {2~a p b p b=(1zb), which is proportional to a p b p , the average number of mRNAs derived from small bursts, rather than the frequency of the small bursts. Since g {2 &1 (Fig. 2c) and b&4, our interpretation of the filtered data implies that a p b p &1:25, which is close to the estimate obtained from the model (Table 3).
Evidently, there is a discrepancy between the assumptions of Choi et al. and the implications of our model. To understand its origin, observe that their assumptions are equivalent to the relations Table 3. Burst frequency and size in uninduced cells with and without auxiliary operators.  Table 1.
i.e., they assumed, in effect, that both the mean and the variance are dominated by contributions from small transcriptional bursts. In contrast, (51)-(52) show that small bursts contribute to the mean, but not to the variance. This difference arises because we assumed that looping is so fast that the rapid fluctuations due to partial dissociations are averaged out on the slow time scale of the other processes. This averaging process preserves the contribution of small transcriptional bursts to the mean, but eliminates their contribution to the variance.
The assumption F &b p b appears to be implausible. Indeed, (53) implies that translational bursts contribute the term b to the Fano factor. For the small bursts to make a significant, let alone dominant, contribution to the Fano factor, it is clear that b p *1, i.e., on average, approximately one mRNA must be synthesized per partial dissociation. However, looping is so fast compared to transcription that b p :v 0 =k O1O2 &0:03 in the absence of the inducer (Table 3). Moreover, b p is unlikely to change even in the presence of the inducer since v 0 and k O1O2 are constant over the range of inducer concentrations used in the experiments. We conclude that the bursts due to partial dissociations are so small that they cannot be the dominant source of burstiness.
Interpretation of F and g {2 derived from raw data. Choi et al. rejected the raw data shown in Fig. 2b since the occurrence of large bursts in a few cells distorted the statistics of the small bursts. We show below that these data are a valuable source of information about the statistics of large bursts. Specifically, (53)-(54) predict the observed variation of F and g {2 derived from the raw data, and thus provide a method for estimating not only the size and frequency of the large transcriptional bursts, but also the fraction of proteins derived from them. This method is particularly useful because, as we show below, there are simple relationships between the size and frequency of the large bursts in strains SX701 and SX703, but they are not identical.
The analysis of the raw data shows that the total burstiness, F , increases with inducer concentration (Fig. 2b). Eq. (53) implies that this is due to the growing burstiness of the large transcriptional bursts: Since both b c and f c increase with inducer level, so does f c b c b. This increase occurs so rapidly that at 100 mM TMG, large trancriptional bursts become the dominant source of burstiness, i.e, F &f c b c b. Indeed, assuming b&4, (53) implies F&f c b c b whenever F &5. Inspection of Fig. 2b shows that at 100 mM TMG, F&25, and hence, F&f c b c b. We shall show below that at such inducer levels, f c &1 and F&b c b.
In contrast to the total burstiness, F , the reciprocal of the total noise, g {2~m =F , decreases with inducer concentration until it reaches a constant value (Fig. 2b). The model suggests that this is because both m and F increase with inducer level, but F increases faster than m: Indeed, both b c and f c increase with inducer level, and Eq. (54) shows that m is proportional to the ratio b c =f c , whereas F increases with the product f c b c . The decreasing trend of g {2 continues until the inducer levels become so high that large bursts account for all the proteins (f c &1) and burstiness (F &b c b). Under these conditions g {2 approaches a c , the frequency of large bursts, which is independent of inducer concentration. Comparison with the data in Fig. 2b then implies that a c &0:2.
Given a c &0:2 and b&4, (53)-(54) provide a method for estimating the variation of b c b and f c with inducer levels from the raw data for SX701. To see this, it is convenient to rewrite (53)-(54) in the form Since the variation of F and g {2 with the inducer concentration is known (Fig. 2b), we can solve the above equations to obtain b c b and f c as a function of the inducer concentration. These calculated profiles, shown in Fig. 2d, agree with the claims above: Both b c b and f c increase with the inducer level, and the latter approaches 1 at 100 mM TMG.
Strain without auxiliary operators. Interpretation of F F and g g {2 . Choi et al. assumed that the F F and g g {2 shown in Fig. 2a represent the size and frequency of large transcriptional bursts, i.e., Our model implies that these relations are valid at all non-zero inducer concentrations used in the experiments. Indeed, since b&4, (61)-(62) imply that the above relations are valid whenever F F &5, which is satisfied ( F F * > 25) at all the non-zero inducer concentrations used in the experiments (Fig. 2a). In particular, comparison with the data in Fig. 2a implies that a a c &3.
Relationships between the statistics of large bursts in the strains with and without auxiliary operators. The model predicts simple relationships between the size and frequency of the large transcriptional bursts in strains SX701 and SX703, which provide tests for checking the consistency of the model. Indeed, it follows from (48) and (57) that b c = b b c~1 =3, a relationship that is also mirrored by the data (compare full and dashed lines in Fig. 2d). Similarly, (47) and (56) imply that a c a a c~3 a ratio estimated to be 1/80 based on the values in Table 1, which is of the same order of magnitude as the value 1/15, obtained from the experimentally determined values of a c &0:2 and a a c &3.

Condition for the negative binomial distribution
Choi et al. assumed that the protein distributions of both strains follow the Gamma distribution, the continuous analog of the negative binomial distribution. We have shown above that neither one of the strains follows the negative binomial distribution. Here, this was not true, (77) and (80) imply that the burstiness would be independent of inducer concentration, which contradicts the data. The negative binomial distribution is therefore unlikely to provide good fits to the raw data for both strains, but will fit the filtered data well, since the contribution of large bursts has been eliminated from it. The fits in Choi et al. are consistent with this conclusion. The Gamma distribution fits the filtered data for strain SX701 rather well. However, this is less so for the protein distributions obtained with strain SX703, which exhibits only large bursts. Figure 4 shows that better fits are obtained with the negative hypergeometric distribution (38).

Conclusions
We formulated and solved a stochastic model of lac expression accounting for auxiliary operators and DNA looping. Based on a comparison of our expressions for the Fano factor, noise, and protein distribution of strains SX701 (with auxiliary operators) and SX703 (without auxiliary operators) with those proposed by Choi et al., we arrive at the following conclusions: 1. The physical interpretations of the Fano factor F F and reciprocal noise g g {2 for strain SX703 are identical to those proposed by Choi et al., namely F F and g g {2 represent the size and frequency of (large) transcriptional bursts. 2. The physical interpretations of the Fano factor F and reciprocal noise g {2 derived from the filtered data for SX701 differ from those given by Choi et al., namely F and g {2 represent the size and frequency of small transcriptional bursts. Instead, we find that F represents the size of translational bursts, and g {2 is proportional to the mean number of mRNAs derived from small transcriptional bursts. Our interpretation is different because we assume that looping is so fast that fluctuations due to small transcriptional bursts are averaged out -small bursts therefore contribute to the mean, but not the burstiness, of the protein distribution. This has two consequences: (a) The information lost due to the averaging implies that the small burst size and frequency cannot be separately extracted from the data. At best, we can only determine the product of the small burst size and frequency, which represents the mean number of mRNAs derived from small bursts. (b) The burstiness is entirely due to translational and large transcriptional bursts. In particular, the burst size derived from the filtered data for strain SX701, from which the contribution of the large bursts has been deliberately eliminated, yields the size of translational, rather than small transcriptional, bursts.
3. Choi et al. did not consider the raw data for SX701 because large bursts, although rare, contributed significantly to protein synthesis. This is consistent with our model: Even in uninduced cells, 20% of the proteins are derived from large bursts. We find that the raw data contains valuable information about the statistics of large bursts. By analyzing this data with our model, we isolate not only the size and frequency of large bursts, but also the fraction of proteins derived from them. The large burst size obtained in this manner is consistent with another prediction of the model, namely, it is one-third of the (large) burst size in strain SX703. The model also predicts that the fraction of proteins derived from large bursts is completely determined by a measurable quantity, namely the dissociation constant for binding of the repressor to the auxiliary operator O 2 . 4. The protein distributions for both strains are not negative binomial: SX703 follows a negative hypergeometric distribution, and SX701 follows a mixture of the negative binomial and negative hypergeometric distributions that reflects the existence of two sub-populations of proteins, namely, those derived from small and large bursts. Negative binomial distributions are attained only if large bursts are insignificant, a condition that holds only if the data are filtered by eliminating the contribution of such bursts.
These results imply that interpretation of the steady state protein distributions depends crucially on the details of the regulatory mechanisms.