Quantitative Understanding of Probabilistic Behavior of Living Cells Operated by Vibrant Intracellular Networks

Yu Rim Lim, Ji-Hyun Kim, Seong Jun Park, Gil-Suk Yang, Sanggeun Song, Suk-Kyu Chang, Nam Ki Lee, and Jaeyoung Sung Department of Chemistry and Institute of Innovative Functional Imaging, Chung-Ang University, Seoul 156-756, Korea Department of Chemistry, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA Department of Physics and School of Interdisciplinary Bioscience and Bioengineering, POSTECH, Pohang 790-784, Korea (Received 9 November 2014; revised manuscript received 18 March 2015; published 10 August 2015)


I. INTRODUCTION
To be viable, living cells must control the chemical fluctuation in intracellular reaction processes within a certain range [1][2][3]. Meanwhile, ubiquitous chemical fluctuations originating from various sources cause phenotypic variations in the morphological or embryonic development, the cellular decision making, and the cell fate [4][5][6][7]. Over a dozen years, impressive progress has been made in fluorescence imaging techniques that visualize the chemical fluctuation among living cells [8][9][10][11][12][13][14][15][16][17]. Achieving a quantitative understanding of the intracellular chemical fluctuation and its consequence for probabilistic biological behaviors of living cells is one of the challenging goals in modern biophysical sciences.
In quantitative investigation of intracellular chemical fluctuation and its consequence, Poisson network models have been employed and analyzed by the master equation devised by Pauli, which has been the most popular theoretical model and has thus been exploited for many years [18][19][20][21][22][23][24][25][26][27]. The key assumption in the conventional network model is that the rate coefficient of an elementary reaction in intracellular networks is constant, or rarely a function of time, that is the same across reaction vessels [ Fig. 1(a)]. However, according to modern single-molecule studies, the rate coefficient of an enzyme fluctuates with the stochastic conformational dynamics even under highly homogeneous environments [28,29]. For biopolymer reactions occurring in living cells, the rate coefficients are stochastic variables that differ from cell to cell and fluctuate over time because of their coupling to heterogeneous and dynamic environments [30,31], which adds additional complexities [ Fig. 1(b)].
The coupling of an intracellular network to cell environments has profound effects on the chemical fluctuation produced by the network [33][34][35][36][37]. For example, gene expression variability among isogenic cells is greatly affected by the coupling of the gene network to environmental states, such as the configuration of chromosomes, the gene regulation state, the nutrition state of cells, phase in the cell cycle, and populations of gene expression machinery proteins, which are often beyond direct measurement. The dual reporter method was proposed for a convenient separate estimation of the "intrinsic noise" and "extrinsic noise" [30], which has been widely used, and there were other propositions for the most appropriate decomposition scheme or definition of the intrinsic and extrinsic noises [38][39][40][41][42]. However, to our knowledge, an accurate quantitative description of a stochastic intracellular network interacting with the cell environment is still missing. This is because it is difficult to construct a quantitative model that correctly describes the interactions of an intracellular network with complex and mostly hidden cell environments [ Fig. 1(c)]. Because of the lack of information about the coupling of the intracellular network to the cell environment with enormous degrees of freedom, it is a state-of-the-art task to find the correct model of the environmental coupling of intracellular networks in terms of a few discrete chemical states and Poisson transition processes between them, the only balls and sticks in the conventional model. For this reason, a successful quantitative analysis of stochastic gene expression with the use of the conventional network models has been rare [24][25][26][27]. The chemical fluctuation-dissipation theorem was introduced to provide a general description of the chemical noise produced by the conventional gene network models [21]. A modification of the conventional gene network model and the master equation was made to describe the non-Poisson transcription or translation process [16,[43][44][45]. A few pioneering computational studies were carried out for specific models of the dynamically fluctuating rate coefficient in the gene expression network, which provided insight about the effects of the dynamic cell environment on stochastic gene expression [31,36].
However, yet to be developed is a general model or theory that makes it possible to take into account the effects of hidden cell environment on the stochastic outcome of intracellular networks in a complete manner [39,41], and an accurate quantitative analysis of the gene expression noise still remains a challenging task [39].
Here, we introduce a new type of theoretical model and stochastic kinetics for an efficient quantitative analysis of intracellular networks interacting with hidden cell environment. Hereafter, a reaction process whose rate coefficient is a stochastic variable coupled to a hidden environment will be designated as a vibrant reaction process and represented by a wavy arrow in distinction from the usual Poisson reaction process that is represented by a plain arrow [ Fig. 1(c)]. As a vibrant intracellular reaction is coupled to hidden cell environment, it is difficult to construct an explicit, quantitative model for the coupling of the vibrant reaction process to the dynamic hidden cell state. We could show that the stochastic outcome of the vibrant reaction process cannot be accurately described merely by introducing a time-dependent rate coefficient into the conventional master equation (see Figs. 1 and 2 in the Supplemental Material [32]). The rate coefficient of a vibrant reaction process is a stochastic variable coupled to a large number of hidden environmental variables, which cannot be represented by a time-dependent rate coefficient or any deterministic function of time. However, the present analysis reveals that various environmental effects on the product number fluctuation of intracellular reaction networks can be collectively and completely characterized by the timecorrelation function of the product creation rate fluctuation, more specifically, its Laplace transform with the Laplace variable replaced by the product decay rate. Therefore, an FIG. 1. Concept of vibrant reaction process coupled to complex and hidden cell environment. Schematic representations of (a) an elementary reaction process under a homogeneous environment and (b) a vibrant reaction process with the rate coefficient being a stochastic variable that differs from cell to cell and fluctuates over time ( Fig. 1 of Supplemental Material [32]). (c) Coupling of a reaction system to a complex and hidden cell environment makes the rate-coefficient stochastic variable dependent on the cell environmental state, the correct, explicit model of which cannot be easily constructed for lack of available information. A reaction process whose rate coefficient is coupled to hidden environment will be designated by vibrant reaction process and represented by the wavy arrow. accurate quantitative analysis of the product number fluctuation of a vibrant reaction process requires only a correct modeling of the dependence of the time-correlation function of the rate fluctuation on the control parameter in our experiment, which is quite feasible, as we demonstrate in the latter part of this work, and it does not require an explicit modeling of the environmental dynamics and its coupling to the system network.

II. RESULTS AND DISCUSSION
A. Response of the product number fluctuation to environmental fluctuations As the first example of a vibrant reaction process, we consider an enzyme reaction with product creation rate RðΓÞ dependent on environmental state Γ. Here, Γ collectively represents all the environmental variables, including the number of reactant molecules and enzymes, which completely specifies a state of cell environment. The product creation rate R may be dependent on part or all of the components of the environmental state vector Γ. Because the population of reactants and enzymes and other cell state variables differ from cell to cell and fluctuate over time, the time evolution of the environment-coupled reaction rate is a stochastic dynamic process. We assume that RðΓÞ can be factored into RðΓÞ ¼ θðN X ; N E ÞkðΓ 0 Þ, where θðN X ; N E Þ denotes the rate factor dependent on the numbers N X and N E of the reactant molecules X and enzymes E, and kðΓ 0 Þ denotes the rate coefficient dependent on environmental state variables Γ 0 other than N X and N E . For simplicity, we assume that θ and k are independent from each other and that product molecules decay with a constant rate γ.
For the reaction network, the joint probability density ψ z ðΓ; tÞ that a cell contains z product molecules and the hidden cell environment is at state Γ at time t obeys the following generalized master equation [48]: Here, E AE1 z is the mathematical operator defined by E AE1 z gðzÞ ¼ gðz AE 1Þ for any function g of z. LðΓÞ denotes the mathematical operator describing the time evolution of ψ z ðΓ; tÞ due to dynamics of the hidden environment. The time evolution of environmental state Γ can be an arbitrary dynamic process. For example, if the environmental state evolves according to the Langevin equation, LðΓÞ is the Fokker-Planck operator. If the environmental state evolves according to Newton's equation of motion, LðΓÞ is the Liouville operator. LðΓÞ also inclusively describes changes in ψ z ðΓ; tÞ due to hidden reaction processes coupled to the system network. The explicit mathematical forms of LðΓÞ and RðΓÞ are unavailable because of the lack of information in most cases. However, one can still obtain a general result for the product number fluctuation arising from the influence from hidden environment dynamics and its coupling to the system network, as we show below. ψ z ðΓ; tÞ satisfies the normalization condition R dΓ P ∞ z¼0 ψ z ðΓ; tÞ ¼ 1. The nth moment hz n i of the steady-state product number distribution is defined as lim t→∞ R dΓ P ∞ z¼0 z n ψ z ðΓ; tÞ. From Eq. (1), we obtain the exact results for the mean hzi and the relative variance or noise η 2 z ð¼ hδz 2 i=hzi 2 Þ of the number z of product molecules in the steady state (see Appendix A for the derivation procedure): Throughout, hxi represents the average of quantity x over an ensemble of cells in the steady state; η 2 x designates the relative variance, or the noise, in quantity x. The result has an intuitively appealing structure; the product number fluctuation is proportional to the fluctuation in the two factors of the reaction rate. The proportionality coefficient or the susceptibility is dependent on dynamics of the rate fluctuation and the product decay rate, which is found to be χ zs ¼ γ R ∞ 0 dte −γt ϕ s ðtÞ½≡γφ s ðγÞ ðs ∈ fk; θgÞ, with ϕ s ðtÞ being the normalized time-correlation function defined by ϕ s ðtÞ ≡ hδsðtÞδsð0Þi=hδs 2 i and δs ¼ s − hsi. Similarly, the susceptibility of the product number fluctuation to the bilinear rate factor fluctuation η 2 k η 2 θ is found to be χ zðk;θÞ ¼ γ R ∞ 0 dte −γt ϕ k ðtÞϕ θ ðtÞ. The magnitude of the susceptibility χ zs or χ zðk;θÞ increases with the relaxation time scale of the rate factor fluctuation. The susceptibility has its maximum value, unity, when the rate factor is statically distributed, i.e., when ϕ s ðtÞ ¼ 1 ðs ∈ fk; θgÞ, but it vanishes when the relaxation time scale of the rate fluctuation is much smaller than the product lifetime, γ −1 (see Appendix B). Our result shows that the effect of environmental dynamics and its coupling to the reaction rate is manifested in the product number fluctuation through the time-correlation function of the environment-coupled rate fluctuation. That is to say, the product noise is dependent on the two-time correlation function of the product creation rate, but it is not sensitive to other microscopic details of the environmental dynamics and its coupling to the reaction rate, so there can be numerous environment-coupled network models that yield the same product noise.
We emphasize that Eq. (2) exactly holds regardless of the environmental dynamics including the population dynamics of reactants and enzymes, and mathematical forms of the environmental-state-dependent rate coefficient kðΓÞ, and the reactant population-dependent rate factor θðN X ; N E Þ. The correctness of Eq. (2) could be confirmed by stochastic simulation for a couple of exactly solvable models of environment-coupled reaction processes (Fig. 2). As shown in Fig. 2, the prediction of Eq. (2) for Fano factor F z ð¼ η 2 z hziÞ is in excellent agreement with accurate stochastic simulation results for two different models of the environment-coupled reaction process. In the simulation of the vibrant reaction network models with the rate coefficients being stochastic processes, one cannot use Gillespie's stochastic simulation method because the conventional simulation method was developed to simulate Poisson reaction network models in which the rate coefficients of the elementary reactions do not fluctuate at all [49]. For an accurate stochastic simulation of the vibrant reaction model whose rate coefficient is dependent on a stochastic variable, we use the numerical algorithm developed in Ref. [50]. Details of the algorithm we use here are presented in Note 1 of the Supplemental Material [32].
The results in Figs. 2(b) and 2(c) show that F z − 1, which measures the deviation of the product number distribution from the Poisson distribution, approaches zero when the value of γ is large enough. In the large γ limit, all of χ zk , χ zθ , and χ zðk;θÞ in the right-hand side (rhs) of Eq. (2b) approach unity because of the Tauberian theorem, i.e., lim γ→∞ γfðγÞ ¼ lim t→0 fðtÞ, so that the non-Poisson product noise arising from the rate fluctuation saturates to the maximum value, η 2 k þ η 2 θ þ η 2 k η 2 θ . In contrast, the first term on the rhs of Eq. (2b), hzi −1 , which is the Poisson noise independent of the rate fluctuation, linearly increases with γ for any value of γ because of Eq. (2a) (see Fig. 3 in the FIG. 2. Effects of environment-coupled fluctuations in the rate coefficient on the fluctuation in the product number. (a) Schematic representation of the vibrant reaction and ensuing decay of the product molecule with constant rate γ. The reaction rate R is factored into the reactant population-dependent rate factor θ and the rate coefficient k, which are dependent on environmental variable r. η 2 z ð≡hδz 2 i=hzi 2 Þ is given by Eq. (2). (b),(c) Comparison between the prediction of Eq. (2) for Fano factor F z ð¼ η 2 z hziÞ and stochastic simulation results on two different models of environment-coupled reaction. We focus on the situation where the fluctuation in the reaction rate R is dominated by the fluctuation in the rate coefficient k by setting θ equal to unity. F z − 1 is proportional to F k ð¼ hδk 2 i=hkiÞ with the proportionality coefficient given byφ k ðγÞ ¼ R ∞ 0 dte −γt ϕ k ðtÞ with ϕ k ðtÞ being the normalized time-correlation function of rate-coefficient fluctuation [Eq. (C1)]. Examples of stochastic time trace of k, which is a function of stochastic variable rðtÞ (left-hand panel), ϕ k ðtÞ (middle panel), and F z ðγÞ − 1 (right-hand panel), are displayed for the following models. (b) kðrÞ ¼ k 0 þ κr 2 with stochastic dynamics of r being the stationary Gaussian Markov process known as Ornstein-Uhlenbeck process [46]. λ stands for the relaxation rate of the fluctuation in r defined by hrðtÞrð0Þi ¼ hr 2 ie −λt=2 . F z is a monotonically decaying function of γ. (c) kðrÞ ¼ κ exp½−βr with stochastic dynamics of r being the stationary Gaussian non-Markov process describing the position of the weakly damped harmonic oscillator with mass m under the frictional force 2mλdrðtÞ=dt [47]. For this model, hrðtÞrð0Þi is given by hr 2 ie −λt ½cosðωtÞ þ ðλ=ωÞ sinðωtÞ with ω ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where ω 0 is the natural frequency of the oscillator. When ϕ k ðtÞ is an oscillatory function, F z − 1 can be a nonmonotonic function of γ. See Appendix C for the details of the models. Supplemental Material [32]). Therefore, when the value of γ is large enough, the Poisson product noise dominates over the non-Poisson product noise; i.e., the ratio of the latter to the former, which is the same as F z − 1, is approximately given by hkihθiðη 2 k þ η 2 θ þ η 2 k η 2 θ Þ=γ and vanishes in the large γ limit.
When the fluctuation of the rate coefficient is absent and the rate factor θ is simply the reactant number, Eq. (2b) correctly reduces to the chemical fluctuation-dissipation theorem in Ref. [21].
B. Gene expression system: Synthetic Escherichia coli with controllable RNA polymerase level To demonstrate the applicability of our new stochastic kinetics to quantitative analysis of the intracellular chemical fluctuation, we analyze gene expression statistics in our synthetic Escherichia coli system [51]. In the E. coli system, T7 RNA polymerase (RNAP) fused with yellow fluorescent protein (YFP) is expressed, the expression level of which could be monitored by fluorescence microscopy and controlled by the concentration of inducer anhydrotetracycline (aTc). In the chromosome of the E. coli system, the genes encoding two fluorescent proteins, cyan fluorescent protein (CFP) and mCherry, are incorporated under the control of identical copies of a T7 RNAP promoter. This experimental system was previously used to show that the cell-to-cell variation in the upstream RNAP level actually propagates to the extrinsic noise component of the downstream protein level fluctuation but not to the intrinsic noise component, when both components are estimated by the dual reporter method [51].
Here, we quantitatively analyze various gene expression statistics of the system with the use of our model and stochastic kinetics in order to establish the mathematical relationship of various gene expression statistics to the environment-coupled rate fluctuations in the three major elementary steps composing gene expression: the RNAPpromoter association, the transcriptional elongation [synthesis of messenger RNA (mRNA) from DNA by RNAP], and the translation (synthesis of protein from mRNA by ribosome). For this purpose, we investigate the single gene expression variability and the correlation between the expression levels of CFP and mCherry genes, controlling the mean T7 RNAP level among the cells by changing the inducer concentration. In doing so, we keep the relative variance in the T7 RNAP level among the sampled cells maintained at around the typical value, ∼0.1, of the abundant protein noise in E. coli, using a careful cell sampling method (Fig. 4 in the Supplemental Material [32]). A comparison between our predictions and experimental results is made also for the dependence of gene expression statistics on the noise in the RNAP levels. It is observed that gene expression statistics mostly depend on the mean and variance of the T7 RNAP level among the sampled cells, and they are quite robust with respect to changes in the detailed shape of the T7 RNAP level distribution or other details of our cell sampling method (Figs. 5 and 6 in the Supplemental Material [32]).

C. Vibrant gene expression network model
For a quantitative analysis of the gene expression statistics among the synthetic E. coli system with controllable RNAP level, we construct the vibrant gene expression network model shown in Fig. 3(a), in which the rate parameters associated with the transcription and translation processes are stochastic variables coupled to hidden cell environment. As the transcription begins with the reversible association of RNAP to the promoter, which is followed by the transcriptional elongation, we model the transcription rate as R TX ¼ θk TX , where θ and k TX denote, respectively, the RNAP-occupation fraction of the promoter and the rate coefficient of the transcriptional elongation. As in the Michaelis-Menten enzyme kinetics, θ is related to the RNAP level N Rp by θ ¼ KN Rp =ð1 þ KN Rp Þ (Fig. 7 in the Supplemental Material [32]). However, in contrast to the conventional enzyme kinetics, the affinity K of the promoter for RNAP is not a constant but a stochastic variable coupled to hidden environmental state, which may include the promoter regulation state, the levels of transcription initiation factors, and the configuration of DNA. KN Rp ð≡ζÞ is the dimensionless variable measuring the promoter-RNAP interaction strength, which we will refer to often. We explicitly model θ because it is dependent on the RNAP level, the control variable in our experiment. The translation rate is modeled as R TL ¼ mk TL , where m and k TL denote the intracellular mRNA level and the translation rate coefficient, respectively. The rate coefficients of the transcriptional elongation and translation processes are dependent on hidden environmental variables including the cellular levels of transcription elongation factors and ribosomes. However, in the present approach, the dependence of the rate coefficients on those hidden environmental variables is not modeled explicitly because of a lack of information. Nevertheless, one can still develop an accurate mathematical description of the vibrant gene expression network model, neither making a conjecture about the dependence of K, k TX , or k TL on hidden environment nor assuming a particular environmental network, or Markov chain of the rate parameters for the description of the environment-coupled rate fluctuations. For simplicity, we do not consider any feedback regulation or any correlation between k TX and θ or k TL and m; we assume the decay rate of mRNA (protein) is constant γ m (γ p ). p;θ ð¼ β pθ η 2 θ Þ and η 2 p;k ð¼ β pk TX η 2 k TX þ χ pk TL η 2 k TL Þ, propagated from environment-coupled fluctuations in θ and from that in k TX and k TL . The RNAP noise among the sampled cells is controlled to be 0.1 AE 0.002. The error bar or the standard deviation of each experimental result is obtained over 10 different samples, each of which contains about 1500 cells. The size of the error bar is smaller than or comparable to the symbol size (see Fig. 6 of Supplemental Material [32] we obtain accurate analytic results for the mean and variance of the mRNA and protein levels, m and p, in the steady state (Appendix D): The analytic results have an intuitively appealing mathematical structure that could be conjectured from Eq. (2). Noting that Eq. (2) is derived for the case with the product creation and decay rates being R ¼ θkðΓÞ and γ, respectively, and that the creation rate and decay rate of mRNA (protein) are given by R TX ¼ θk TX ðΓÞ ½R TL ¼ mk TL ðΓÞ and γ mðpÞ , respectively, one can guess that the mean and variance of the mRNA (protein) level are given by Eqs. (2a) and (2b) with rate coefficient k, reactant populationdependent rate factor θ, product decay rate γ, and product level z being replaced by k TX ðk TL Þ, θðmÞ, γ mðpÞ , and mðpÞ, respectively. We find the conjecture is actually correct (Appendix D). The meaning and definition of each symbol in Eqs. (4) and (5) are the same as those of the corresponding symbol in Eq. (2).
In the limit where the fluctuations in θ, k TX , and k TL are negligible, Eq. (5b) reduces to η 2 p ¼ hpi −1 þ hmi −1 γ p = ðγ p þ γ m Þ, with γ p =ðγ p þ γ m Þ being χ pm in this limit [see Eq. (S6-2) in the Supplemental Material [32]]. The latter result was previously obtained from the conventional gene network model without rate-coefficient fluctuations [20,21] and was first identified as "intrinsic noise" in Ref. [20]. The intrinsic noise may be defined as the part of gene expression noise that can be described by a gene network model without any environment-coupled fluctuation in the rate parameters, while the extrinsic noise may be defined as the remaining part of the gene expression noise that arises from a coupling of the gene network model to external environments, i.e., the gene expression noise arising from the fluctuations in the rate parameters of the gene network. However, there has been controversy over the most appropriate definition of intrinsic noise or extrinsic noise [21,30,38,39], which is discussed in Appendix E in detail. An account of the relationship of the present theory to previous theories in this field is given in Note 5 of the Supplemental Material [32].
By substituting Eq. (4b) into Eq. (5b), we finally obtain the protein noise as the sum of the intrinsic noise, the extrinsic noise dependent on the RNAP-promoter interaction, and the extrinsic noise independent of the RNAPpromoter interaction: pm stands for χ pm in the absence of transcription rate fluctuation. β pθ or β pk TX denotes the susceptibility of the protein noise to the fluctuation in the transcription rate factor θ or k TX . Analytic expressions of α, β pθ , and β pk TX are presented in Appendix D. Equation (6) [16]. In addition, we could make an experimental estimation between the intrinsic noise and the extrinsic noise of the single gene expression noise for each of CFP and mCherry in our synthetic E. coli system (see Fig. 8 in the Supplemental Material [32] and Appendix F).
The first term on the rhs of Eq. (6) is designated by ðη int p Þ 2 and it may correspond to the previously identified intrinsic noise, which persists in the absence of environmental fluctuations. However, ðη int p Þ 2 is slightly dependent on the environmental fluctuations as discussed in Appendix D. The first term on the rhs of Eq. (6) is inversely proportional to the mean RNAP occupancy hθi of the promoter because the mean protein level hpi is given by with hpi max ½≡hk TX ihk TL i=ðγ m γ p Þ. Equation (7) can be obtained by substituting the expression of hmi in Eq. (4a) into Eq. (5a). Figures 4(a) and 4(b) show the best global fit of Eq. (7) to the experimental data for the dependence of the mean protein level on the mean RNAP level and RNAP noise. As the mean RNAP occupancy hθðζÞi of the promoter is slightly dependent on fluctuation in the RNAP-promoter interaction strength ζ, so are hpi and ðη int p Þ 2 . When η 2 ζ ¼ 0, hpi obeys the Michaelis-Menten formula; i.e., hpi=hpi max 1 ζ=ð1þζÞ ½≡θðζÞ, withζ ¼KN Rp . However, when η 2 ζ ≠ 0, hpi=hpi max or hθi is slightly smaller than θðζÞ [Eq. (S2-5) in the Supplemental Material [32]], which is consistent with our experimental results [ Fig. 4(c) and Supplemental Material Fig. 9 [32]].
The sum of the second and third terms on the rhs of Eq. (6) may be identified as the extrinsic noise ðη ext p Þ 2 that originates from the environmental fluctuations. The second term, η 2 p;θ , particularly designates the extrinsic protein noise originating from environment-coupled fluctuation in RNAP occupancy of the promoter. When the fluctuation in the RNAP-promoter interaction is not too large, the noise in the RNAP occupancy of the promoter can be related to the mean and noise of the promoter-RNAP interaction strength by η 2 θ ≅ η 2 ζ ½1 − θðζÞ 2 [Eq. (S4-6) of the Supplemental Material [32]], so that the protein noise originating from the fluctuation in RNAP occupancy of the promoter is given by which is in excellent agreement with our experimental results [ Fig. 4(d)]. The result tells us that the extrinsic protein noise originating from the RNAP-promoter association step is proportional to the fluctuation in RNAPpromoter interaction strength ζ and the square of the probability that the promoter is not occupied by RNAP. η 2 p;θ has the maximum value η 2Ã p;θ ¼ β pθ η 2 ζ in the small RNAP-promoter interaction limit where θðζÞ ¼ 0 orζ ¼ 0.
The third term, η 2 p;k , collectively designates the protein noise propagated from environment-coupled fluctuations in the rate coefficients of transcriptional elongation and translation processes, so this term is independent of the RNAP-promoter interaction. In our system, the RNAPpromoter association step rather than the ensuing gene expression steps makes the major contribution to the extrinsic protein noise for a wide range of the RNAP level ( Fig. 10 of the Supplemental Material [32]); η 2 p;θ is greater than η 2 p;k as long as θðζÞ < 1 − ffiffiffi q p with q ≡ η 2 p;k =η 2Ã p;θ . The value of q is estimated to be 0.098 for CFP and 0.12 for mCherry (Table I). Parameter q serves as a measure of the relative magnitude between the protein noise originating from the RNAP-promoter association step and that (c) Dependence of the mean protein level scaled by its maximum or the mean RNAP occupancy of promoter, hpi=hpi max ½¼ hθðζÞi, on θðζÞ. Inset: Dependence of hθðζÞi onζ (see Fig. 9 of Supplemental Material [32]). hθi ¼ θðζÞ for η 2 ζ ¼ 0 (dashed line). η 2 ζ ¼ 0.44 and 0.59 for CFP and mCherry (Table I). (d) Dependence of η 2 p;θ scaled by its maximum η 2 * p;θ on θðζÞ. For both CFP and mCherry, given by Eq. (F1) and η 2 p;k given in Table I originating from the transcriptional elongation and translation steps. The small value of parameter q indicates that the extrinsic noise originating from the RNAP-promoter association step makes the dominant contribution to the extrinsic noise in our synthetic E. coli system with T7 RNAP transcription machinery.
E. Relative contribution of two extrinsic noise components η 2 p;θ and η 2 p;k is manifested in the dependence of Fano factor on the mean RNAP occupancy of promoter In the absence of environment-coupled fluctuations in the RNAP-promoter bound fraction and the rate coefficients of the transcriptional elongation and translation processes, η 2 p;θ and η 2 p;k in Eq. (6) vanish so that the protein level Fano factor is simply given by F p ð≡hδp 2 i=hpiÞ ¼ η 2 p hpi ¼ 1 þ α, a constant independent of the mean RNAP occupancy hθi of the promoter. However, it has been experimentally observed that the protein level Fano factor is far from being constant in hθi and its dependence on hθi varies depending on the gene expression system [24]. An important question here is what information about extrinsic noise can be extracted from the dependence of the protein level Fano factor F p ðhθiÞ on the mean RNAP occupancy. We find that the shape of F p ðhθiÞ contains information about the relative ratio of the two extrinsic noise components, η 2 p;θ and η 2 p;k , on the rhs of Eq. (6) or the value of parameter qð≡η 2 p;k =η 2Ã p;θ Þ. From Eqs. (6) and (7), we obtain F p as the following third order polynomial in hθi: given that θðζÞ ≅ hθi, which is the case in our experimental system On the other hand, when q > 1=3 or when 3η 2 p;k > η 2Ã p;θ , the Fano factor monotonically increases with hθi [Figs. 5(a) and 5(c); see also Fig. 11 in the Supplemental Material [32]].
According to the theoretical prediction, the Fano factor of CFP or mCherry expression in our E. coli system with T7 RNAP transcription machinery should have the single maximum at about hθi * ≅ 1=3 because the value of parameter q is about 0.1, far smaller than 1=3 in the system. Indeed, this is found to be the case [Fig. 5(d)]. The nonmonotonic dependence of the protein level Fano factor on the mean RNAP occupancy of the promoter was previously observed also for yEGFP gene expression under various transcription control mechanisms in Saccharomyces cerevisiae with similar size and shape [24] Table I in the Supplemental Material [32]]. The strong nonmonotonic shape of the Fano factor data suggests that the RNAPpromoter association step rather than the ensuing transcription or translation step makes the dominant contribution to the yEGFP gene expression variability in the S. cerevisiae system investigated in Ref. [24].
It is interesting that Eqs. (6) and (7) also provide a unified quantitative explanation of the global trend in the mean protein level dependence of the protein noise previously reported for the entire genome of wild-type E. coli and a comprehensive set of S. cerevisiae genes [ Fig. 13(a) of the Supplemental Material [32]] [14,16], with only three adjustable parameters: α, η 2Ã p;θ , and η 2 p;k (Table II of the Supplemental Material [32]). The extracted value of TABLE I. Optimized values with the standard error of the adjustable parameters in Eqs. (6), (7), and (10) for our synthetic E. coli system with T7 RNAP transcription machinery. Each of the equations is compared with 54 experimental data points to extract the values of the three parameters in each row. Details of the data analysis are presented in Appendix I. The unit of K M in Table I  parameter q ð¼η 2 p;k =η 2Ã p;θ Þ is found to be much greater than 1=3 for these systems (Table III in the Supplemental Material [32]), for which the Fano factor should monotonically increase with the mean RNAP occupancy of the promoter or the mean gene expression level, as mentioned above. We find that this is indeed the case as shown in Fig. 13(b) of the Supplemental Material [32], where we display the experimental data for the protein level Fano factor for these systems and the result of Eq. (8) with the predetermined parameter values in Table II of the Supplemental Material. The linear dependence of the Fano factor on the mean protein level suggests that the extrinsic protein noise originating from the RNAP-promoter association step is far smaller than that from the ensuing transcriptional elongation and translation steps.
According to the previous work investigating the dependence of gene expression variability on the promoter strength in wild-type E. coli, the Fano factor of the protein level expressed in wild-type E. coli almost linearly increases with the mean RNAP occupancy of the promoter [10], in contrast to the Fano factor of CFP or mCherry expressed in our genetically modified E. coli with T7 RNAP transcription machinery. We confirm that, when a mCherry gene is expressed with E. coli RNAP under various E. coli RNAP promoters, the Fano factor almost linearly increases with the mean RNAP occupancy of the promoter or the mean mCherry level (Fig. 14 of the Supplemental Material [32]). Note that the Fano factor given in Eq. (8) becomes almost linear in the mean RNAP occupancy hθi of promoter when η 2 p;k is much greater than η 2Ã p;θ . Consistently, the value of parameter q ð¼η 2 p;k =η 2Ã p;θ Þ extracted from the Fano factor of mCherry expressed in wild-type E. coli is estimated to be of the order of 10 4 , which indicates that the protein noise originating from the RNAP-promoter association step is quite small for the mCherry expression by wild-type E. coli gene expression machinery.
It was recently reported that the noise of abundant proteins in wild-type E. coli is dominantly contributed from the transcription process consisting of the RNAPpromoter association and the ensuing transcriptional elongation [51]. Therefore, the result of our quantitative analysis in the previous paragraph indicates that the transcriptional elongation step rather than the RNAP-promoter association step makes the dominant contribution to the  6) and (7) with the predetermined parameters values extracted from the analysis in Fig. 3 (Table I). Equations (6) and (7) also provide a unified quantitative explanation of all the Fano factor data shown in (e) with a small number of adjustable parameters (Table I of Supplemental Material [32]). The extracted value of q is far smaller than 1=3 for data in (d) and (e), where the Fano factor has the maximum at about hθi The small value of parameter q means that the extrinsic noise originating from RNAP-promoter association step makes the dominant contribution to the total extrinsic protein noise. In the opposite case, q becomes very large and the Fano factor becomes a linear function of hθi, which is the case for the abundant proteins in the wild-type E. coli in (f) [see also Figs. 13(b) and 14 in Supplemental Material [32]]. See text for the analytic expression of F p ðhθiÞ.
cell-to-cell variation in the abundant protein levels in wildtype E. coli.
From our analysis of the global trends in the expression statistics for a comprehensive set of S. cerevisiae genes, shown in Fig. 13 in the Supplemental Material, we also find that the fluctuation in the size and shape of S. cerevisiae cells contributes to the protein noise originating from the post-transcriptional or translation process much more than it does to the protein noise originating from the RNAPpromoter association process. More details are presented in Appendix M. The relevance of the cell size and shape fluctuation on the cellular control over the gene expression levels varies depending on the cell type. In E. coli, the protein noise is found to be not so sensitive to the fluctuation in the cell size and shape [16].
F. Intrinsic (extrinsic) noise estimated by the dual reporter method is found to be close to the average of single gene intrinsic (extrinsic) noises of two reporters We make a separate estimation between the intrinsic noise ðη int p Þ 2 and extrinsic noise ðη ext p Þ 2 for each of CFP and mCherry by analyzing the dependence of the protein noise on the mean RNAP occupancy of the promoter in Sec. II D. The definition of the intrinsic noise and extrinsic noise composing the single gene expression noise is given below Eq. (6). In comparison, the dual reporter method [30], which has been widely used for separate estimation between the intrinsic and extrinsic noises, relies on different definitions of the intrinsic and extrinsic noises. In the dual reporter method, the extrinsic noise is defined as where p A and p B represent the expression levels of two reporter genes. Recently, a legitimate question was raised regarding the consistency between the intrinsic (extrinsic) noise of the single gene expression and the intrinsic (extrinsic) noise estimated by the dual reporter method [39] (Appendix E).
To investigate this issue, we obtain the analytic results for the intrinsic and extrinsic noises defined in the dual reporter method by considering the dual gene network model shown in Fig. 6(a) (Appendix J), and find the following two inequalities: where ðη int p Þ 2 and ðη ext p Þ 2 designate the intrinsic noise part and the extrinsic noise part of the single gene expression noise defined below Eq. (7). The latter inequality is consistent with the simulation result reported in Ref. [39]. The intrinsic (extrinsic) noise defined in the dual reporter method can be the same as the average single gene intrinsic (extrinsic) noise when the environment-coupled fluctuation of rate parameters in the expression network of gene A has the perfect positive correlation with that of the corresponding rate parameters in the expression network of gene B (Appendix J). To our surprise, the dual reporter method provides a good estimation of the intrinsic and extrinsic parts of the single gene expression noise in our system (Fig. 15 of the Supplemental Material [32]). This result indicates that, in our system, the primary source of the extrinsic noise is the fluctuations in the global environmental factors that influence the expression of every gene in the same fashion, which produce a positive correlation between the expression levels of two independent genes.
G. Environment-induced correlation between the expression levels of reporter genes has two components with different origins and dependence on the promoter strength Fluctuations in the global environmental variables such as the cellular levels of RNAP, transcription factors, and ribosomes produce a positive correlation between the expression levels of two independent genes without any cross-regulation mechanism [17,30]. To quantitatively investigate the environment-induced correlation between the expression levels of two independent genes, we consider the dual gene expression model shown in Fig. 6(a) and compare the prediction of the theoretical model with our experiment results for the dependence of the environmentinduced correlation between CFP and mCherry on the RNAP level.
The stochastic kinetic equation setup for the dual gene expression network in Fig. 6(a) is given by where ψ m A ;m B ;p A ;p B ðΓ; tÞ denotes the joint probability density that a cell contains m A copies of mRNA and p A copies of protein created from gene A and m B copies of mRNA and p B copies of protein created from gene B and that the hidden environment coupled to the gene expression network is at state Γ at time t. L R X ðm X ; p X ; ΓÞ is the mathematical operator describing the time evolution of the joint probability density due to the creation and annihilation processes of mRNA and protein generated from gene X∈fA; Bg, which is given by R TXðTLÞ (blue surface). (f)-(h) Dependence of C p , ϕ p , and hδp A δp B i on θðζ A Þ. In our system, θðζ B Þ is dependent on θðζ A Þ, as shown in the inset of (f). They are not the same because the value ofK A is different from that ofK B (Table I). The red (blue) line represents the contribution of the red (blue) surface in (c)-(e). In (f), the theory curve represents the best fit of Eq. (10) to experimental data. In (g) and (h), the theory curves represent the prediction of Eqs. (6), (7), and (10) with predetermined parameter values in Table I. (i) C p linearly increases with the probability that both promoters of CFP and mCherry are unoccupied by RNAP. (j) The joint distribution of CFP and mCherry levels represented by heat map (experiment) and the contour plot (joint gamma distribution). See Fig. 16 of Supplemental Material [32] and Appendix K for more details. Near the most probable protein level, ðp * A ; p * B Þ, the joint distribution becomes the two-dimensional Gaussian (Note 8 in and Fig. 17 in Supplemental Material [32]) with standard deviation σ þ (σ − ) along the positive (negative) diagonal direction.
The value of the RNAP noise η 2 N Rp is the same as in Figs. 3(d) and 3(e). analogous to that used in the single gene network model in Fig. 3(a). As in Eq. (3), LðΓÞ denotes the mathematical operator describing the time evolution of ψ m A ;m B ;p A ;p B ðΓ; tÞ due to dynamics of a hidden environment coupled to the gene expression network, whose explicit mathematical form is unknown. Applying P ∞ m B ;p B ¼0 ðÁ Á ÁÞ or P ∞ m A ;p A ¼0 ðÁ Á ÁÞ on both sides of Eq. (9), one can recover Eq. (3) for ψ m A ;p A ðΓ; tÞ or ψ m B ;p B ðΓ; tÞ. From Eq. (9), we obtain the analytic result of the mean scaled correlation C p ½≡hδp A δp B i=ðhp A ihp B iÞ between the expression levels of two independent genes (Appendix H). The result reads as with C p;θ ¼ β C pθ C θ and C p;k ¼ β C pk TX C k TX þ χ C pk TL C k TL . Here, C x denotes the mean scaled correlation between quantity x A associated with gene A and quantity x B associated with gene B; i.e., Note that the correlation between the protein levels has two different components: C p;θ originating from the correlation C θ between the RNAP occupancies of the promoters and C p;k originating from the correlation C k TX or C k TL between the rate coefficients of transcriptional elongation or translation processes of the two reporter genes. The proportionality coefficient β C ps or χ C ps ðs ∈ fθ; k TX ; k TL gÞ denotes the susceptibility of the product correlation to the source correlation C s , whose analytic expression is presented in Appendix H.
Equation (10) provides a quantitative explanation of our experimental data for the dependence of C p on the mean and noise of RNAP level variability [ Fig. 6(b)]. The values of the three parameters extracted from the quantitative analyses are collected in the third row of Table I. The first term on the rhs of Eq. (10) originates from the environment-induced correlation C θ between the RNAP occupancy of the promoter of the one reporter gene and that of the promoter of the other reporter gene. Because C θ is approximately given by C θ ≅ C ζ ½1 − θðζ A Þ½1 − θðζ B Þ withζ AðBÞ ¼K AðBÞNRp [see Eq. (S7-6) in Note 6 of the Supplemental Material [32]], C p;θ ð¼β C pθ C θ Þ can be written as with C Ã p;θ ¼ β C pθ C ζ . C p;θ given in Eq. (11) are shown as the red surface in Fig. 6(c) and the red lines in Figs. 6(f) and 6(i). Equation (11) tells us that C p;θ is directly proportional to ½1 − θðζ A Þ½1 − θðζ B Þ. The latter quantity is approximately equal to the probability that both promoters of the reporter genes are simultaneously unoccupied by RNAP, because θðζ X Þ ≅ hθðζ X Þi [ Fig. 4(c)]. It makes sense because, unless both promoters are unoccupied by RNAP at the same time, intracellular RNAP level fluctuations could not affect the expression of both genes simultaneously to produce any correlation in the gene expression. It can be shown that η 2 p;θ should be the same as C p;θ for the ideal, symmetric dual reporter system with perfectly correlated fluctuations in the environmentcoupled parameters (Appendix J). Therefore, the quadratic dependence of η 2 p;θ on 1 − θðζÞ shown in Fig. 4(d) is consistent with Eq. (11).
In comparison, the second term, C p;k , on the rhs of Eq. (10) originates from environment-coupled correlated fluctuations in the transcriptional elongation and translation processes, so it is independent of the RNAP unbound probability of the promoter [blue surface in Fig. 6(c) and the blue lines in Figs. 6(f) and 6(i)]. Consequently, C p comprising C p;θ and C p;k should be linear in the probability that both promoters of the reporter genes are simultaneously unoccupied by RNAP, as demonstrated in Fig. 6(f).
Very recently, it was reported that C p is a linearly increasing function of RNAP noise [51] [ Fig. 6(b)]. The experimental can be easily understood by noting that the first term on the rhs of Eq. (10) or C p;θ given in Eq. (11) is linearly proportional to the RNAP noise, because C Ã p;θ ∝ C ζ ¼ ð1 þ C K Þη 2 N Rp þ C K (Note 7 in the Supplemental Material [32]), and the second term, C p;k , on the rhs of Eq. (10) is independent of the RNAP noise because it originates from the environment-coupled fluctuations in the transcriptional elongation and translation processes.
With Eqs. (6), (7), and (10) at hand, one can also calculate the standardized correlation ϕ p and the covariance hδp A δp B i, which are other measures of the environmentinduced correlation between p A and p B , by ϕ p ¼ C p =ðη p A η p B Þ and hδp A δp B i ¼ C p hp A ihp B i [Figs. 6(d) and 6(e)]. Using the latter equations and the values of the parameters in Table I (Table I), θðζ A Þ is also different from θðζ B Þ, as shown in the inset of Fig. 6(f).
Our analyses of the data in Figs. 6(f)-6(h) show that the environment-induced correlation between the expression levels of the dual reporter genes can be decomposed into two components, one originating from environmental fluctuations in the RNAP-promoter association step, which is represented by the red lines, and the other originating from the environmental fluctuations in the ensuing gene expression steps, represented by the blue lines. The relative contributions of the two components to the environmentinduced correlation between the expression levels of the dual reporter genes are various depending on the RNAP-promoter interaction strength and the type of the experimental measure of the environment-induced correlation between the gene expression levels. However, to any measure of the environment-induced correlation between the expression levels of the dual reporter genes, the contribution of environmental fluctuations in the RNAP-promoter association step can be greater than that of environmental fluctuations in the ensuing gene expression steps only when the probability that both promoters of the reporter genes are unbound to RNAP is greater than some critical value, i.e., only when ½1−θðζ A Þ½1−θðζ B Þ>q C ð≡C p;k =C Ã p;θ Þ. The value of parameter q C is found to be 0.12 in our system [ Fig. 6(i) and Table I].
It is interesting that the standardized correlation ϕ p is quite robust with respect to changes in the RNAP-promoter interaction strengthζ X or in θðζ X Þ in contrast with the mean scaled correlation C p or the absolute covariance hδp A δp B i between the expression levels of dual reporter genes. The result of the present analysis shown in Fig. 6(g) indicates that ϕ p contributed from the environmental fluctuations in the RNAP-promoter interaction step (red line) mostly decreases with the RNAP-promoter interaction strength, whereas ϕ p contributed from the correlated environmental fluctuation in the ensuing gene expression steps (blue line) increases with the RNAP-promoter interaction strength, so that the sum of the two contributions becomes nearly independent of the RNAP-promoter interaction strength.
The entire shape of the joint distribution is well described by the joint gamma distribution [ Fig. 6(j) herein and Fig. 16 of the Supplemental Material [32]; see also Appendix K]. Near its maximum, the joint distribution of the standardized protein levels are approximated by the two-dimensional Gaussian (Note 8 of the Supplemental Material [32]), which elongates along the positive diagonal direction as long as (Fig. 17 of the Supplemental Material [32]); the ratio of the standard deviation σ þ along the positive diagonal direction to σ − along the negative diagonal direction is given by σ . We find that σ þ =σ − is always greater than the ratio of the extrinsic noise to the intrinsic noise as shown in Fig. 6(k), although the two quantities were considered to be the same in quite a bit of the literature [30,33,[53][54][55]. In an ideally symmetric dual reporter system in which hp A i and η 2 p A are exactly the same as hp B i and η 2 p B , σ þ =σ − is given by which is always greater than η ext dual =η int dual [see Eq. (L4)]. In an asymmetric dual reporter system, σ þ =σ − is greater than . That is to say, σ þ =σ − is greater than η ext dual =η int dual in any case. The stochastic kinetics introduced herein can be easily extended to investigate the more complex gene networks or other types of intracellular networks interacting with a hidden environment, which is expected to contribute to advances in diverse fields of modern science where the quantitative understanding of the chemical fluctuation and its consequence is or will be of great importance. Some predictions and suggestions for new experiments are presented in Appendix M.

APPENDIX B: SUSCEPTIBILITY OF THE PRODUCT NOISE TO THE RATE FLUCTUATION INCREASES WITH THE TIME SCALE OF RATE FLUCTUATION
In Eq. (2b), the magnitude of the susceptibility, χ zs or χ zðk;θÞ , increases with the fluctuation time scale of the rate factor fluctuation. We consider the simple example where ϕ s ðtÞ ¼ expð−t=τ s Þ with s ∈ fk; θg, where τ s is the relaxation time scale of s. We obtain χ zs ¼ γτ s =ð1 þ γτ s Þ and χ zðk;θÞ ¼ γτ ðk;θÞ =ð1 þ γτ ðk;θÞ Þ, where τ ðk;θÞ denotes the reduced lifetime defined by τ ðk;θÞ ¼ τ k τ θ =ðτ k þ τ θ Þ. Since τ ðk;θÞ −1 ¼ τ k −1 þ τ θ −1 , τ ðk;θÞ increases with τ k or τ θ . The result shows that the susceptibilities χ zk , χ zθ , and χ zðk;θÞ increase with the fluctuation time scales, τ k and τ θ , of the rate factors k and θ. In the statically disordered limit where ϕ k ðtÞ ¼ ϕ θ ðtÞ ¼ 1, χ zs and χ zðk;θÞ assume the maximum value, unity. However, due to the dynamic fluctuation of the rate factors, the susceptibilities are less than unity and dependent on the detailed shape of the time-correlation function ϕ s ðtÞ of the rate factor s ∈ fk; θg and the inverse lifetime γ of the product molecule.
It is general characteristics of dynamic heterogeneity that effects of heterogeneity decrease with the speed of dynamics. For example, in spectroscopy, the absorption line shape of a chromophore broadened by heterogeneous environments surrounding the chromophore gets narrower as the speed of environmental dynamics increases [56,57]. As another example, we mention Zwanzig's enzyme model with a fluctuating bottleneck whose catalytic rate is determined by the area of the bottleneck. The effects of the heterogeneity in the reaction rate on the shape of enzymatic turnover time distribution decreases with the speed of the rate fluctuation [58].

APPENDIX C: DETAILS OF THE VIBRANT REACTION MODELS CONSIDERED IN FIG. 2
In the absence of the rate fluctuation, the product noise becomes the Poisson noise; i.e., η 2 z ¼ 1=hzi. However, the reaction rate fluctuation makes the product noise deviate from the Poisson noise. The Fano factor F z ð≡hδz 2 i=hzi ¼ η 2 z hziÞ, which measures deviation of the product noise from the Poisson noise, depends on the dynamic fluctuation of the product creation rate through the single quantity,φ R ðγÞF R ; i.e., When ϕ R ðtÞ is a monotonically decaying function of t, F z ðγÞ is also a monotonically decaying function of γ [ Fig. 2(b)]. However, when ϕ R ðtÞ is an oscillating function of t [59], F z ðγÞ can be a nonmonotonic function of γ [ Fig. 2(c)].
Note that the response of the product number fluctuation to the environment-coupled rate-coefficient fluctuation is determined only by the time-correlation function of the rate-coefficient fluctuation and the product decay rate irrespective of microscopic details of environmental dynamics or its coupling to the rate coefficient [Eq. (C1)]. The correctness of the theoretical result is tested against stochastic simulation results for a couple of exactly solvable models.
In the model considered in Fig. 2(b), the rate coefficient is given by kðrÞ ¼ k 0 þ κr 2 , where r is the stationary Gaussian Markov process, the mean and the timecorrelation function of which are given by hrðtÞi ¼ 0 and ϕ r ðtÞ ¼ hrðtÞrð0Þi=hr 2 i ¼ e −λ r t , respectively. hr 2 i denotes the variance of r. The stochastic process rðtÞ is known as the Ornstein-Uhlenbeck (OU) process [46]. The time-evolution operator LðrÞ describing the stochastic dynamics of the OU process is well known as For example, the time evolution of the probability distribution ψðr; tÞ describing the position of a strongly damped harmonic oscillator at time t is given by ∂ t ψðr; tÞ ¼ LðrÞψðr; tÞ [60]. The reaction model considered in Fig. 2(b) is similar to Zwanzig's enzyme, whose catalytic rate is controlled by the product escape process out of the enzyme's active site though a bottleneck channel with a fluctuating radius [58]. For both models, the time-correlation function ϕ k ðtÞ ½¼hδkðtÞδkð0Þi=hδk 2 i of the rate-coefficient fluctuation is given by ϕ k ðtÞ ¼ e −λt , with γ þ λ. By substituting the latter equation into Eq. (C1), we get which is a monotonically decreasing function of γ þ λ.
In the model considered in Fig. 2(c), the rate coefficient is given by kðrÞ ¼ κ exp½−βr. κ and β are constant parameters. r is a stationary Gaussian non-Markov process with the mean and the time-correlation function given by hrðtÞi ¼ 0 and ϕ r ðtÞ ¼ e −λt ½cosðωtÞ þ ðλ=ωÞ sinðωtÞ with ω ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ω 2 0 − λ 2 p , respectively. The best known example of the latter stochastic process may be the position rðtÞ of a weakly damped harmonic oscillator with mass m and natural frequency ω 0 under a weak frictional force, 2mλdrðtÞ=dt [61]. This non-Markov process can be represented by two-dimensional coupled Markov processes [62]. In this model, our environmental state vector has two components, Γ ¼ fr; vg, and the corresponding timeevolution operator is given by The rate coefficient of an electron transfer reaction between two molecular units of which separation is fluctuating around a mean value may be represented by the model. For this model, ϕ k ðtÞ is given by Substituting Eq. (C5) into Eq. (C1), we calculate the Fano factor F z ðγÞ numerically. For this model, it is not easy to obtain the exact expression of F z ðγÞ; however, when parameter β 2 hr 2 i is small, F z ðγÞ can be approximated by F z ðγÞ in Eq. (C6) is a nonmonotonic function of γ that has a peak centred around γÃ ¼ ω 0 − 2λ when ω 0 > 2λ.
When ω 0 ≤ 2λ, F z ðγÞ becomes a monotonically decaying function of γ. In Note 10 in the Supplemental Material [32], we present details of the stochastic simulation method that we use to confirm the correctness of our theoretical results for the two exactly solvable models discussed here.

APPENDIX D: DERIVATION PROCEDURE OF EQS. (4) AND (5)
For the gene expression network model in Fig. 3(a), the joint probability density ψ m;p ðΓ; tÞ that a cell contains m mRNAs and p proteins and the hidden environment coupled to the gene expression network is at state at time t satisfies Eq. (3) [48]. ψ m;p ðΓ; tÞ satisfies the normalization condition R dΓ P ∞ m;p¼0 ψ m;p ðΓ; tÞ ¼ 1. The nth moment hq n i of the steady-state product number distribution is defined as lim t→∞ R dΓ P ∞ m;p¼0 q n ψ m;p ðΓ; tÞ, with q ∈ fm; pg. Applying P ∞ p¼0 ðÁ Á ÁÞ to both sides of Eq. (3), we obtain the following equation: Equation (D1) is formally equivalent to Eq. (1) so that the derivation procedure leading to Eq. (4) is also equivalent to Appendix A. Subsequently, by denoting fΓ; mg as Γ, we can rewrite Eq. (3) as whereLðΓÞ designates the operator governing dynamics of Γ, i.e., Equation (D3) also has exactly the same mathematical structure as Eq. (1) or Eq. (D1). Therefore, we obtain the analytic results of hpi and hpðp − 1Þi as Eq. (5) simply by replacing ðm; R TX ; k TX ; θÞ in Eq. (4) by ðp; R TL ; k TL ; mÞ. The microscopic expression of susceptibility χ zq is given by where ϕ q ðtÞ denotes the normalized time-correlation function of the fluctuation in quantity q. In addition, we define ϕ ðq 0 ;q 00 Þ ðtÞ as ϕ ðq 0 ;q 00 Þ ðtÞ ≡ ϕ q 0 ðtÞϕ q 00 ðtÞ. Before going further, we note that the susceptibilities involving ϕ m ðtÞ, χ pm , and χ pðm;k TL Þ require special attention, because they depend on the RNAP level, the control variable in our experiment. We obtain the analytic expression of ϕ m ðtÞ in the Laplace domain, from Eq. (A1). In Eq. (D5), η 2 R TXφ R TX ðsÞ or η 2 R TX ϕ R TX ðtÞ is related to η 2 θ or RNAP level fluctuation [see Eqs. (S4-6) and (S5-3) in the Supplemental Material [32]]. With Eq. (D5) at hand, we obtain the analytic expression for χ pm η 2 m and χ pðm;k TL Þ η 2 m in Eq. (5b) as follows: where f q ðtÞ and g q ðtÞ denote, respectively, The time-domain expression of ϕ m is given by ϕ m ðtÞ ¼ e −tγ m þ ðη 2 R TX =η 2 m Þf R TX ðtÞ, because Laplace transforms of e −γ m t and f q ðtÞ are given by ðs þ γ m Þ −1 and f q ðsÞ ¼ γ 2 m ½φ q ðγ m Þ −φ q ðsÞ=ðs 2 − γ 2 m Þ, respectively, both of which appear in Eq. (D5). Using Eqs. (D6) and (D7), χ pm η 2 m þ χ pðm;k TL Þ η 2 m η 2 k TL on the rhs of Eq. (5b) can be rewritten as pk TX η 2 χ ð1Þ pq ¼ γ p ½f q ðγ p Þ þ η 2 k TLĝ q ðγ p Þ½q ∈ fθ; k TX ; ðθ; k TX Þg: ðD12Þ By substituting Eq. (D10) into Eq. (5b), we reach Eq. (6), in which α ¼ χ pm and χ pk TL are defined by Eqs. (D11) and (D4), respectively; in addition, β pθ and β pk TX are given by pðθ;k TX Þ Þη 2 k TX ; ðD13Þ pk TX : ðD14Þ The term in the first set of parentheses on the rhs of Eq. (D13) and the first term on the rhs of Eq. (D14) show how transcriptional rate fluctuation propagates into the protein noise; there exist two different pathways: mRNA copy number fluctuation-mediated and direct pathways. Equation (13) also tells us that the propagation of the noise in the RNAP-promoter association step to the protein noise is affected by the rate fluctuation in the transcriptional elongation step.

APPENDIX E: DEFINITIONS OF THE INTRINSIC AND EXTRINSIC NOISES
Over the past decade, gene expression noise has often been decomposed into two parts: intrinsic noise and extrinsic noise. However, several different definitions of the terminology have been used without clear distinctions [17,20,21,30,31,36]. The intrinsic noise, coined by Thattai and Oudenaarden, can be defined as the part of gene expression noise that can be described by a gene network model without any environment-coupled fluctuation in the rate parameters, while the extrinsic noise can be defined as the remaining part of the gene expression noise that arises from a coupling of the gene network model to external environments, i.e., the gene expression noise arising from the fluctuations in the rate parameters of the gene network. According to the latter definition, the first term on the rhs of Eq. (6) in the main text can provide an estimation of the intrinsic noise, which can also be obtained from the conventional master equation approach, and the remaining two terms can be identified as the extrinsic noise for our gene expression model.
In some literature, the protein noise originating from the mRNA level fluctuation is identified as the extrinsic noise even if it is produced by the conventional gene expression network without any fluctuation in the rate parameters; the intrinsic noise is used to designate the Poisson protein noise that is persistent in the absence of the mRNA level fluctuation.
In the dual reporter method, which has been widely used for separate estimation of the intrinsic and extrinsic noises, the extrinsic noise is defined as the mean scaled correlation C p ½≡ðη ext dual Þ 2 between two reporter protein levels so that the intrinsic noise can be defined as [30,31]. The sum of the intrinsic and extrinsic noises defined in the dual reporter method yields the mean total protein noise; i.e., In literature referring to the dual reporter method, another definition of the intrinsic noise, has been frequently used [30]. However, the latter definition has a drawback; its sum with extrinsic noise C p does not yield the mean total noise: In real experiments, hp A i cannot be exactly the same as hp B i. Therefore, given that the total noise should be the sum of the intrinsic and extrinsic noises, i is a more appropriate definition of the intrinsic noise in the dual reporter method.

APPENDIX F: ESTIMATION OF THE INTRINSIC
AND EXTRINSIC NOISES, ðη int p Þ 2 AND ðη ext p Þ 2 By analyzing the dependence of the protein noise on the RNAP level or the mean protein level, we extract the values of adjustable parameters, ð1 þ αÞ=hpi max , β pθ , and η 2 p;k , for each of CFP and mCherry in our E. coli system (see Appendix I and Table 1). With the value of ð1 þ αÞ=hpi max at hand, we estimate the intrinsic noise by where hIi denotes the mean fluorescence intensity that is proportional to the mean protein level, hIi ∝ hpi. The value of hIi max ð∝hpi max Þ can be extracted from the dependence of the mean fluorescence intensity on the mean RNAP level. From Eq. (S2-4) of the Supplemental Material [32], we obtain where ζ is the dimensionless RNAP-promoter interaction strength defined by ζ ¼ KN Rp . The noise η 2 ζ in ζ is given by

in the Supplemental
Material [32]). By analyzing the dependence of hIi onN Rp and η 2 N Rp with the use of Eq. (F2), we extract the values of hIi max ,K, and η 2 K (Table I). The extrinsic noise ðη ext p Þ 2 is estimated simply by subtracting the intrinsic noise estimated above from the experimentally measured value of the total protein noise η 2 p ; i.e., APPENDIX G: FANO FACTOR OF THE PROTEIN LEVEL From Eqs. (6) and (7), we obtain the analytic results for the Fano factor F p ð≡hδp 2 i=hpiÞ of the protein level as a function of the RNAP-occupation probability hθi of the promoter: where θðζÞ is recursively related to hθi by Equation (G2) is obtained from Eq. (S2-3) of the Supplemental Material [32] and the definition of θðζÞ, θðζÞ ≡ζ=ð1 þζÞ. When η 2 ζ ≪ 1, θðζÞ is approximately the same as hθi (Fig. 9 of the Supplemental Material [32]). When η 2 ζ is not negligible, one can calculate θðζÞ from Eq. (G2) for given values of hθi and η 2 ζ , either by using successive iteration with initial trial, θðζÞ ¼ hθi, or by finding the positive real root of the cubic equation, The full analytic expression of θðζÞ in terms of η 2 ζ and hθi is omitted here. In the case where η 2 ζ is so small that θðζÞ can be approximated by hθi, the Fano factor given in Eq. (G1) becomes a simple third-order polynomial of ð1 − hθiÞ: where F p ð1Þ and q are given by F p ð1Þ ¼ 1 þ α þ hpi max η 2 p;k and q ¼ η 2 p;k =η 2Ã p;θ , respectively. Note that the shape of the regularized Fano factor is determined by q only in the small η 2 ζ regime. When q < 1=3, the Fano factor given in Eq. (G3) has the maximum value at hθi ¼ 3 −1 ð2 − ffiffiffiffiffiffiffiffiffiffiffiffiffi 1 − 3q p Þ and the minimum value at hθi ¼ 3 −1 ð2 þ ffiffiffiffiffiffiffiffiffiffiffiffiffi 1 − 3q p Þ. On the other hand, when q > 1=3, the Fano factor is a monotonically increasing function of hθi (Fig. 11 of the Supplemental Material [32]).
In experiment, the Fano factor F I of the fluorescence intensity emitted from the fluorescent protein is measured, which is related to the Fano factor F p of the protein level by F I ¼ hδI 2 i=hIi ¼ chδp 2 i=hpi, with c being the proportionality constant defined by c ¼ hIi=hpi. For this reason, the absolute value of F I is not the same as that of F p . However, Eq. (G3) can still be used to analyze F I ðhθiÞ because Therefore, the relative shape of F I ðhθiÞ or the regularized Fano factorFðhθiÞ provides us with the information about the value of q ¼ η 2 p;k =η 2Ã p;θ , the ratio of the protein noise originating from the transcriptional elongation and ensuing translation to the maximum of the protein noise originating from the stochastic RNAP-promoter interaction [see Eq. (6) in the main text].
In the case where η 2 ζ is not negligible compared with unity, we use the more accurate equation for the analysis of F I ðhθiÞ: where θðζÞ is given by Eq. (G2).

APPENDIX H: DERIVATION PROCEDURE OF EQ. (10)
By applying P ∞ p A ;p B ¼0 ðÁ Á ÁÞ on both sides of Eq. (9), we obtain The method of solution for Eq. (H2) is similar to the method of solution for Eq. (A2) presented in Appendix A. By following a similar line of derivation, one can obtain the steady-state expression of hδm A δm B i as follows: which seems analogous to Eq. (A12). Using Eqs. (H3) and (H4), we obtain the expression for hδm A δm B i=hm A ihm B i. By identifying fΓ; m A ; m B g asΓ 0 , one can obtain the evolution equation of ψ p A ;p B ðΓ 0 ; tÞ from Eq. (9), which has exactly the same mathematical structure as Eq. (H1). Therefore, by following a similar line of derivation to above, one can obtain the expression for The resulting expressions for C q ð≡hδq A δq B i=hq A ihq B iÞ with q ∈ fm; θ; k TL ; k TX g are given by Here, χ C zq with ðz; qÞ ∈ fðm; k TX Þ; ðp; k TL Þ; ðp; mÞ; ðm; θÞg denotes the susceptibility of correlated noise C m in the mRNA level ðz ¼ mÞ or the susceptibility of correlated noise C p in the protein level ðz ¼ pÞ to the upstream correlated noise C q between quantities q A and q B .
The microscopic expression of susceptibility χ C zq is given by whereγ z is defined byγ z ¼ γ ðAÞ z γ ðBÞ z =ðγ ðAÞ z þ γ ðBÞ z Þ. For q ∈ fθ; k TX ; m; k TL g, ϕ XY q ðtÞ denotes the normalized timecorrelation function defined by ϕ XY q ðtÞ≡hδq ðXÞ ðtÞδq ðYÞ ð0Þi= hδq ðXÞ δq ðYÞ i. For q ¼ ðθ; k TX Þ or q ¼ ðm; k TL Þ, ϕ XY q ðtÞ is defined by ϕ XY q ðtÞ ¼ ϕ XY θ ðtÞϕ XY k TX ðtÞ or ϕ XY q ðtÞ ¼ ϕ XY m ðtÞϕ XY k TL ðtÞ, respectively. As in the case of the single gene expression considered in Appendix D, the susceptibilities χ C pm and χ C pðm;k TL Þ are dependent on ϕ XY m ðtÞ, which is, in turn, dependent on the RNAP level, the control parameter in our experiment. The Laplace transform of ϕ XY m ðtÞ is given aŝ Note that ϕ XY m ðtÞ can be simply expressed in terms of f XY q ðtÞ as ϕ XY m ðtÞ ¼ e −tγ ðXÞ m þ ðC R TX =C m Þf XY R TX ðtÞ. Using Eqs. (H11) and (H12), we can rewrite χ C pm C m þ χ C pðm;k TL Þ C m C k TL on the rhs of Eq. (H6) as where χ Cð0Þ pm and χ By substituting Eq. (H13) into Eq. (H6), we obtain Eq. (10), in which β C pθ and β C pk TX should read as pðθ;k TX Þ ÞC k TX ; ðH16Þ Note that Eqs. (H16) and (H17) have similar structures to Eqs. (D13) and (D14) (see also Notes 6 and 7 in the Supplemental Material [32]).

APPENDIX I: DATA ANALYSIS
The raw experimental data used in the present work were obtained in Ref. [51] for the synthetic E. coli in which T7 RNAP fused with yellow fluorescent protein (YFP) was expressed from a tetracycline (tet) promoter in a low-copy number plasmid (pNL001). The expression level of T7 RNAP could be controlled by varying the concentration of inducer, anhydrotetracycline (aTc), inhibiting the tet-regulatory action of the tet repressor (tetR). In the chromosome of the E. coli system, the genes encoding two reporter proteins, cyan fluorescent protein (CFP) and a kind of red fluorescent protein (RFP), mCherry, were incorporated under the control of identical copies of T7 RNAP promoter Φ10.
Each cell data point contains the YFP, CFP, and mCherry fluorescence levels in each cell.
Cell data are collected from a clonal population of cells in the steady state at each of the six different inducer (aTc) concentrations. At each inducer concentration, a large collection of cell data are prepared by merging the cell data obtained from three independent experiments (Fig. 18 in the Supplemental Material [32]).
(1) The raw cell data obtained are obtained at six different inducer concentrations. At each inducer concentration, the YFP distribution in the raw cell data could well be approximated by a gamma distribution [51]. We first calculate the mean YFP fluorescence level of the total cell data at each inducer concentration. (2) From the raw cell data prepared at each inducer concentration, we sample cell data by using a gamma distribution of YFP fluorescence levels with nine different values of the relative variance or noise between 0.03 and 0.15. Each sample contains 1500 different cell data (see Fig. 5 in the Supplemental Material [32]). The nine subsets sampled from the raw cell data have the same mean YFP fluorescence level as that of the raw cell data from which they are sampled. In total, we prepare 54 different cell data subsets, each of which has different mean and noise in the YFP levels. We choose a gamma distribution for the YFP fluorescence level according to the previous report that the distribution of the expression level of E. coli genes could be well be approximated by a gamma distribution [16]. (3) For each of 54 different sets of the mean and noise values of YFP fluorescence levels, we repeat the sampling process to prepare ten different subsets of cell data with the same mean and noise in the YFP level. For each of the ten different subsets thus prepared, we calculate the mean hpi and noise η 2 p of the CFP and mCherry levels, the mean scaled correlation C p between the two fluorescence levels. Then we calculate the average value of each statistical measure over the ten different subsets. When the sample size is as large as 1500, the relative fluctuation in the values of each statistical measure is so small that it is barely noticeable (Fig. 6 of the Supplemental Material [32]). In this manner, for each statistical measure, we prepare 54 different data points, each of which represents the cell data subset with a particular set of the mean and noise in the YFP level or the T7 RNAP level. The 54 data points are used to determine values of parameters contained in Eqs. (6), (7), and (10). The extracted values of the parameters are collected in Table I. Other measures such as Fano factor of protein level ðF p ¼ η 2 p hpiÞ, normalized correlation ½ϕ p ¼C p =ðη p A η p B Þ or covariance ðhδp A δp B i ¼ C p hp A ihp B iÞ between p A and p B , and two different definitions of intrinsic and extrinsic noise for the dual reporter system ½ðη intðextÞ p Þ 2 ; ðη intðextÞ dual Þ 2 are not separately fitted to theory but predicted by using the values of parameters given in Table I. (4) The mean and noise (hpi and η 2 p ) of the CFP and mCherry levels, and the mean scaled correlation ðC p Þ between the CFP and mCherry levels obtained for each of 54 subsets with different mean and noise of T7 RNAP levels are quantitatively analyzed by Eqs. (6), (7), and (10). One can rewrite these equations in terms of our control variables, the mean and noise (N Rp and η 2 N Rp ) of T7 RNAP levels, as In the quantitative analysis of the protein noise, the least-squares fit of Eq. (I1) to the CFP noise data and that of Eq. (I1) to mCherry noise data are simultaneously made under the constraint that ðη int p Þ 2 ≤ ðη int dual Þ 2 [see Appendix J, Eq. (J10)]. By analyzing the mean protein level with the use of Eq. (7) or (I2) for cell subsets with various mean and noise,N Rp and η 2 N Rp , of the T7 RNAP (YFP) level, we extract the values of the meanK and noise η 2 K of the RNAPpromoter affinity K, because the two independent variables,ζ and η 2 ζ , in Eq. (7) are related toK and η 2 K byζ ¼KN Rp and η 2 ζ ¼ η 2 K þ η 2 N Rp þ η 2 K η 2 N Rp , respectively. The estimation of the value of η 2 K in the analysis of the mean protein level can be refined by analyzing the protein noise with the use of Eq. (6) or (I1) because the protein noise is more sensitive to η 2 K than the mean protein level is. In accordance with the refinement in the estimated value of η 2 K by the protein noise analysis, the values of the other parameters are refined in a self-consistent manner. All of the estimated parameter values are settled down at the third round of the iterations and do not show any more refinement in the next round. The analysis of the mean scaled correlation C p between CFP and mCherry levels on the basis of Eq. (10) is performed separately from the analyses of the mean and noise of the protein levels. The best-fit values of the adjustable parameters are collected in Table I. In many figures in the present paper, we show dependence of the experimental data on hθðζÞi or θðζÞ, where hθðζÞi and θðζÞ are calculated by hIi exp =hIi max andN exp Rp =ðN exp Rp þ K M Þ, respectively, with the values of hIi max and K M ð¼ 1=KÞ given in Table I. For RNAP, the value of the intensity-to-number conversion factor is estimated to be about 0.13 by visually counting the number of YFP fluorescence spots at the most dilute case (zero inducer concentration). Therefore, unlike hIi assigned for reporter proteins, the separate notation N Rp is assigned for RNAP as the number of RNAP copies per cell.

APPENDIX J: RELATIONS BETWEEN TWO DIFFERENT DEFINITIONS OF THE INTRINSIC AND EXTRINSIC NOISES
Recently, a legitimate question has been raised about consistency between the different definitions of the intrinsic and extrinsic noises [39]. To shed some light on this issue, we make a direct comparison between the different decomposition schemes of the gene expression noise (see Fig. 8 in the Supplemental Material [32]). The extrinsic noise C p , defined in the dual reporter method, would be exactly the same as ðη ext p Þ 2 for an ideal, symmetric dual reporter system with perfectly correlated fluctuations in the environment-coupled parameter. For this system, we have γ ðAÞ pðmÞ ¼ γ ðBÞ pðmÞ ; ðJ1Þ and χ C zk , χ C pm , and χ C mθ given in Eq. (H7) are exactly the same as χ zk , χ pm , and χ mθ given in Eq. (D4), respectively. In this case, β C pθ ð≡χ C pm χ C mθ Þ, C θ , C p;k ð ≡ χ C pk TL C k TL þ χ C pm χ C mk TX C k TX Þ on the right-hand side in Eq. (10) are the same as β pθ ð≡χ pm χ mθ Þ, η 2 θ , and η 2 p;k ð≡χ pk TL η 2 k TL þ β pk TX η 2 k TX Þ, so that Eq. (10) reduces to For the ideal, symmetric dual reporter system with perfectly correlated fluctuations in the environment-coupled parameters, the intrinsic noise defined in the dual reporter method becomes Substituting Eqs. (10) and (J6) on the right-hand side of Eq. (J7), we obtain In general, ðη int dual Þ 2 and ðη ext dual Þ 2 defined in the dual reporter method would assume different values from ½ðη int p A Þ 2 þ ðη int p B Þ 2 =2 and ½ðη ext p A Þ 2 þ ðη ext p B Þ 2 =2 because the expression processes of the two genes, A and B, are not exactly symmetric (see Table I) and also because the fluctuations of the environment-coupled parameters, δq A ðtÞ and δq B ðtÞ, would not be perfectly correlated. Because C q ¼ ϕ q η 2 q ≤ η 2 q even in the perfectly symmetric case, ðη ext dual Þ 2 or C p given in Eqs. (H5) and (H6) is always smaller than ½ðη ext p A Þ 2 þ ðη ext p B Þ 2 =2 with ðη ext p i Þ 2 being the last two terms on the right-hand side of Eq. (6) for the expression noise of gene ið∈ fA; BgÞ, i.e., This means that the intrinsic noise defined in the dual reporter method is always greater than or equal to ½ðη int For our system, the intrinsic and extrinsic noises estimated by the dual reporter method are found to be quite close to ðη int p Þ 2 and ðη ext p Þ 2 , and the former have qualitatively the same dependence on the mean and noise of RNAP level as the latter (Fig. 15 of the Supplemental Material [32]).
It should be noted, at the same time, that the correlation among the expression levels of different genes can be produced not only by the global environmental fluctuations but also by regulatory interactions among the genes [63], which can be explained by employing suitable regulatory network models [31,39].

APPENDIX K: ANALYTIC EXPRESSION FOR THE JOINT PROBABILITY DISTRIBUTION OF THE EXPRESSION OF DUAL REPORTER GENES
The analytic expression of the joint gamma distribution shown in Fig. 6(j) and in Fig. 16 of the Supplemental Material [32], in which the unit of the protein level fluctuation is the standard deviation of the protein level, is given by where P Γ is given in Eq. (K8) and p Ã AðBÞ is given by Other parameters, b i and ϕ p , required to evaluate P Γ in Eq. (K8) can be calculated by Following is a brief account of the joint gamma distribution. The two-dimensional Gaussian distribution shows the key feature of the joint distribution of the two reporter proteins near the most probable population state (see Fig. 16 of the Supplemental Material [32]). However, the marginal distribution of each protein level is non-Gaussian, but better described by gamma distribution (see Fig. 19 of the Supplemental Material [32]) [16]. There are various ways to construct bivariate gamma distributions [64]. Among them, we choose the method proposed by Moran in 1969 to model rainmaking experiments, by which the joint Gaussian near the most probable population state can be naturally transformed into a bivariate gamma distribution [65]. The transformation from the joint Gaussian to the bivariate gamma distribution can be done as follows.
Near the most probable populations, p Ã A and p Ã B , any unimodal joint probability density of p A and p B could be approximated by the two-dimensional Gaussian (see Note 8 of the Supplemental Material [32]), where ξ AðBÞ ¼ ðp AðBÞ − p Ã AðBÞ Þ=σ AðBÞ , with σ AðBÞ being the standard deviation in the protein level p AðBÞ . ϕ p denotes the standardized cross-correlation between p A and p B defined by The one-dimensional probability density of ξ AðBÞ is given by the standard normal distribution, i.e., ð2πÞ −1=2 × expð−ξ 2 AðBÞ =2Þ ¼ R ∞ −∞ dξ BðAÞ P G ðξ A ; ξ B Þ. The cumulative distribution of ξ AðBÞ defines a new random variable u i , where erfcðzÞ ð¼ð2= ffiffiffi π p Þ R ∞ z dte −t 2 Þ denotes the complementary error function. Because the cumulative distribution of any continuous one-dimensional random variable is a uniform distribution in the interval (0,1), u i can also be defined as cumulative distributions of the gamma distributed stochastic variable p i , where ΓðzÞ ð¼ R ∞ 0 dte −t t z−1 Þ and Γða; zÞ ð¼ R ∞ z dte −t t a−1 Þ denote the gamma and incomplete gamma functions, respectively. a i and b i are the two parameters that characterize the gamma distribution, From Eqs. (K3) and (K4), ξ i can be expressed in terms of p i as where erfc −1 ðzÞ denotes the inverse complementary error function. Since we are looking for the joint gamma distribution P Γ ðp A ; p B Þ defined by P Γ ðp A ; p B Þdp A dp B ¼ P G ðξ A ; ξ B Þdξ A dξ B , the joint gamma distribution P Γ ðp A ; p B Þ can be written as From Eq. (K6) and the definitions given in Eqs. (K3) and (K4), we obtain ð2πÞ −1=2 expð−ξ 2 i =2Þdξ i ¼ P i ðp i Þdp i (i ∈ fA; Bg). Substituting the latter equation and Eq. (K2) into Eq. (K7), we obtain in which ξ i is the function of hp i i defined in Eq. (K6). Values of a i and b i specifying P i ðp i Þ in Eq. (K5) can be obtained from the mean and the variance of the protein level, which are given by hp i i ¼ a i b i and σ 2 i ¼ a i b 2 i (i ∈ fA; Bg), respectively, for gamma distributions. That is, We are ready to obtain the expression for the joint distribution , which is shown in Fig. 6(j) and in Fig. 16 of the Supplemental Material [32]. Pðξ A ; ξ B Þ can be calculated from P Γ ðp A ; p B Þ given in Eq. (K8) by with the following relationships: Fig. 6(j), we choose the values of σ A and σ B as the units of the protein level fluctuations, p A − p Ã A and p B − p Ã B , respectively. In the unit system, the theoretical curves for Pðξ A ; ξ B Þ displayed in Fig. 6(j) and in Fig. 16 of the Supplemental Material [32] can be calculated by the equally probable states ðξ A ;ξ B Þ satisfy ln½P G ðξ A ;ξ B Þ= P G ð0;0Þ¼−C 2 , with C being a constant. The equally probable states constitute a contour defined by The contour is elongated along the diagonal direction when ϕ p is positive, but it is elongated along the negative diagonal direction when ϕ p is negative (Fig. 17 of the Supplemental Material [32]). For the contour, C ffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ ϕ p p is the length of the projection onto the diagonal line ξ B ¼ ξ A , while the length of projection onto the negative diagonal line The ratio of the former to the latter is given by The correlated random variables, ξ A and ξ B , can be transformed into two normal components, ξ þ and ξ − , defined by ξ þ ≡ ðξ A þ ξ B Þ= ffiffi ffi 2 p and ξ − ≡ ðξ A − ξ B Þ= ffiffi ffi 2 p . IfPðξ þ ; ξ − Þ denotes the joint distribution of ξ þ and ξ − , one can show that Note that the transformation from ðξ A ; ξ B Þ to ðξ þ ; ξ − Þ corresponds to the counterclockwise rotation of the twodimensional Cartesian coordinate system by π=4 to obtain the new coordinate system. Equation (L2) tells us that the standard deviation σ AE of ξ AE is given by ffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 AE ϕ p p . σ þ is the standard deviation of the projection of the two-dimensional Gaussian given in Eq. (L1) onto the diagonal line, ξ B ¼ ξ A , whereas σ − is the standard deviation of the projection onto the negative diagonal line In the ideally symmetric dual reporter system in which hp A i and η 2 p A are exactly the same as hp B i and η 2 p B , σ þ =σ − is given by In an asymmetric dual reporter system, σ þ =σ − in Eq. (L3) can be written as ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Noting that ðx þ aÞ=ðx − aÞ monotonically decreases with x as long as x > a and thatη 2 p ¼ ðη 2 p A þ η 2 p B Þ=2 is always greater than or equal to η p A η p B , one obtains the following inequality: Because C p is defined as ðη ext dual Þ 2 andη 2 p is the same as ðη int dual Þ 2 þ ðη ext dual Þ 2 , the latter inequality reads as

APPENDIX M: PREDICTIONS AND SUGGESTIONS FOR NEW EXPERIMENTS
Here, we present some discussion and prediction about the effects of cell size fluctuations on the Fano factor of the gene expression level in S. cerevisiae and the nonexponential time-correlation function of mRNA levels. In addition, we present information about the application of the present approach to experimental analyses and make suggestions for new experiments.

Prediction 1: Effects of cell size fluctuations on the Fano factor of eukaryotic gene expression
Our theory predicts that the dependence F p ðhθiÞ of the protein level Fano factor on the mean RNAP occupancy by the promoter would undergo a transition from a nonmonotonic to a monotonic one, as the value of parameter q passes through 1=3 [see Figs. 5(a)-5(c) and also Fig. 11 of the Supplemental Material [32]]. We expect that the transition in F p ðhθiÞ could be observable in some eukaryotic gene expression as the fluctuation in size and granularity of the cells or the radius of cell gate increases.
From our analysis of the global trend in the mean expression level dependence of the expression statistics measured for a comprehensive set of S. cerevisiae genes, shown in Fig. 13 of the Supplemental Material [32], we find that the fluctuation in the size and shape of S. cerevisiae cells contributes to the protein noise originating from the post-transcriptional or translation process much more than it does to the protein noise originating from the RNAPpromoter association process.
In the present work, we introduce parameter qð¼η 2 p;k =η 2Ã p;θ Þ as a measure of the protein noise originating from the post-transcriptional and translation processes relative to the protein noise originating from RNAPpromoter interactions (see Sec. II E). For S. cerevisiae, the value of parameter q increases with the fluctuation in the size or granularity of cells, which can be experimentally controlled by adjusting the radius of the gate used to sample the cells with a certain size and granularity. The value of parameter q extracted from the mean expression level dependence of the expression noise, shown as green dots in Fig. 13 of the Supplemental Material [32], for a comprehensive set of S. cerevisiae genes is found to be 12.8. In comparison, the value of parameter q extracted from the data shown as red dots in Fig. 13 of the Supplemental Material [32] reduces to 0.86, which are the data obtained for the same set of S. cerevisiae genes but expressed in the gated S. cerevisiae cells with similar size and granularity. The yEGFP gene expression data shown in Fig. 5(e) are also obtained for the gated S. cerevisiae cells with similar size and shape, for which the extracted value of q is found to be quite small, 0.00 or 0.068. These results suggest that the reduction of fluctuation in the size and shape of S. cerevisiae cells diminishes the protein noise originating from the post-transcriptional or translation process much more than it does the protein noise originating from the RNAP-promoter association process. The relevance of cell size and shape fluctuation on the cellular control over the gene expression levels varies depending on the cell type. In E. coli, the protein noise was found to be not so sensitive to the fluctuation in the cell size and shape [16].
Since the parameter q, the value of which can be controlled by the gate radius for S. cerevisiae, is one of the crucial factors determining the dependence of protein level Fano factor F p ðhθiÞ on the mean RNAP occupancy hθi of the promoter (Fig. 5), the shape of F p ðhθiÞ should be dependent on the gate radius for S. cerevisiae. As discussed in Sec. II E, if q is less than 1=3, F p ðhθiÞ is a nonmonotonic function of hθi. Otherwise, F p ðhθiÞ is a monotonically increasing function of hθi. Therefore, we predict that the shape of F p ðhθiÞ undergoes a transition from a nonmonotonic to a monotonic one, as the gate radius increases from a fully gated one to an ungated one, so that the value of q passes through 1=3.
By making the comparison between values of the parameters extracted at various values of gate radius, one could investigate which parameters are sensitive to the variation of the cell size fluctuation. By analyzing the relative shape of F p ðhθiÞ on gate radius, one may be able to extract the value of η 2 ζ as well because the relative shape of F p ðhθiÞ depends on η 2 ζ as well as q [see Fig. 11(b) of the Supplemental Material [32]]. η 2 ζ is a quantitative measure of heterogeneity in the RNAP-promoter interaction strength. When the value of parameter q is far smaller than 1=3, hθi * at which the Fano factor reaches its maximum is quite sensitive to the value of η 2 ζ .

Prediction 2: Nonexponential time-correlation function of mRNA level
Our theory predicts that the time-correlation function ϕ m ðtÞ of mRNA level fluctuation deviates from the simple exponential function, expð−γ m tÞ, due to the effects of the transcription rate fluctuation [see Eq. (D5)], which has not been foreseen by any of the previously reported theories. According to Eq. (D5), even if the mRNA survival probability decays following the simple exponential function, expð−γ m tÞ [26], the time-correlation function ϕ m ðtÞ of mRNA level fluctuation deviates from the latter. The difference between the time-correlation function of mRNA level fluctuation and the mRNA survival probability results from the fluctuation in the transcription rate fluctuation, i.e., ϕ m ðtÞ ¼ expð−γ m tÞ þ ðη 2 R TX =η 2 m Þf R TX ðtÞ; ðM1Þ with f R TX ðtÞ defined by Eq. (D8). By analyzing the nonexponential dependence of the time-correlation function of the mRNA level, one could extract the quantitative information about the dynamics of transcription rate fluctuation, or the time-correlation function of the transcription rate fluctuation, which would provide us with the relaxation time scale of the environment-coupled transcription rate fluctuation.

Application of the present formulation to other biological networks or dynamically heterogeneous reaction networks
In the present work, we propose a new theoretical approach to quantitative analyses of the chemical fluctuations produced by intracellular reaction networks interacting with hidden cell environments and demonstrate applications of our approach to the single gene expression and the dual gene expression networks shown in Figs. 3(a) and 6(a). Our theoretical approach is applicable to other intracellular networks or reaction networks coupled to hidden, dynamically heterogeneous environment. In the application, our approach requires an explicit modeling only for the dependence of the system reaction network on the experimental control parameter, but it does not require an explicit modeling of the hidden cell state dynamics and its coupling to the system network in taking into account the effects of the hidden cell environment on the chemical fluctuation produced by the system network. The latter fact is a great advantage of the present approach over the conventional one because it is difficult to construct an explicit and accurate model of hidden cell state dynamics and its coupling to the system network.
In the actual application of the present approach, one should first set up the vibrant reaction network model describing the experimental system and analyze the generalized master equation describing the vibrant reaction network, such as Eqs. (1), (3), and (9), to describe the stochastic dynamics of the system network of interest and hidden cell states. We will soon report a few more applications of the present approach to quantitative interpretation of the chemical fluctuation in living cells. For example, by a rather straightforward application of the present approach, one can provide an excellent quantitative explanation of the inducer concentration dependence of the steady-state mRNA level fluctuation in E. coli reported by Golding and co-workers [26].
As mentioned in the main text, we assume that the product decay process is a Poisson process. An extension of the present formulation is required to take into account the effects of non-Poisson product decay on the product number fluctuation [41]. An extension of the present formulation for a general model of non-Poisson product decay processes will soon appear elsewhere. In the present work, we focus on the case where the product creation rate is independent of the number of product molecules. Extensions of the present formulation to biological networks with various topologies and regulatory interactions are under way. Through the extensions of the present theory, we will be able to establish a set of general rules or analytic results that are directly applicable to quantitative analysis of the chemical fluctuations produced by various types of biological networks.

Suggestion of new experiments
Below, we suggest new experimental observables that could provide us with new quantitative information about the stochastic gene expression in living cells.

a. Time-correlation function of mRNA and protein levels
It would be nice if one could accurately measure the time-correlation function (TCF) of the mRNA level by systematically controlling either the transcription rate fluctuation or the mRNA decay rate. The result could be directly analyzed by Eq. (M1). As predicted in Prediction 2, the time-correlation function of the mRNA level would deviate from the simple exponential function, and the quantitative analyses of the discrepancy between the time-correlation function and the survival probability of the mRNA with use of Eq. (M1) could provide us with quantitative information about the stochastic dynamics of the transcription process, given that the decay process of the mRNA level is a Poisson process [26]. Experimentally, it may be easier to obtain the TCF of the protein level than to obtain the TCF of the mRNA level. The experimental result for the TCF of the protein level could also be quantitatively analyzed by an extension of the present theory, which will be published elsewhere.
For a large number N C of cells, the TCF of the gene expression level q can be estimated by where δqðtÞ stands for the deviation of the gene expression level q i ðtÞ of the ith cell at time t from the average hqðtÞi defined by hqðtÞi ¼ N −1 C P N C i¼1 q i ðtÞ. The variance σ 2 q ðtÞ in the gene expression level is the same as h½δqðtÞ 2 i. When the gene expression process is a stationary process, hδqðt þ t 0 Þδqðt 0 Þi is independent of t 0 , and hqðtÞi and σ 2 q ðtÞ are independent of t. The normalized TCF ϕ q ðtÞ of the gene expression level q can be calculated by ϕ q ðtÞ ¼ hδqðt þ t 0 Þδqðt 0 Þi=σ 2 q ðt 0 Þ. For an ergodic system, Eq. (M2) calculated for an ensemble of a system is the same as the following TCF calculated from a long enough time trace of q i ðtÞ for any single member of the ensemble: However, we expect that the living cell system is not an ergodic system so that the TCF of the gene expression level given in Eq. (M2) would be different from that given in Eq. (M3). It is because the lifetime of a cell is finite and the dynamics of a cell state is not fast enough to span the entire state phase space during the lifetime of the cell. Part of our future research efforts will also be directed to the nonergodicity of biological systems and its consequence to the functions of biological systems.

b. Transcription event counting statistics
It would be wonderful if one could monitor the number of transcription events of a single gene or dual genes with the use of a single-molecule imaging technique, because the time dependence of transcription event counting statistics could provide us with detailed microscopic information about the stochastic transcription dynamics in living cells. The control parameters of the experiment would include the strength of the promoter, the length of the genes to be transcribed, and the type of cells.
The transcription event counting statistics is intimately related to the steady-state mRNA level fluctuation. By analyzing the dependence of the steady-state mRNA level fluctuation on the inducer concentration, we could make a prediction about the transcription event counting statistics, which will be published soon elsewhere.

c. Dependence of the gene expression statistics on the ribosome level fluctuation
Our experimental method can be extended to investigate the effects of the ribosome level fluctuations on the expression statistics of a single reporter gene or dual reporter genes. In this case, the ribosome level becomes our control parameter, so we should explicitly model the dependence of the translation rate on the ribosome level; on the other hand, we do not have to model the dependence of the transcription rate on the RNAP level because the latter is no longer our control variable. It is easy to extend the present theory to analyze the propagation of the ribosome level fluctuation to the gene expression statistics.
For an accurate interpretation of the experiment, it is desirable to investigate the possible dependence of the mRNA survival probability on the ribosome level. Because the translation process is competing with the mRNA degradation process, the change in the ribosome level or the translation rate could induce a significant change in the decay rate of mRNA in cells.
d. Stochastic property of mRNA decay process and its effects on the cellular control over stochastic gene expression It is now well known that the mRNA decay process is important in the cellular control over stochastic gene expression. However, to our knowledge, quantitative investigation on the stochastic dynamics of the mRNA decay process and its effects on the gene expression variability are still missing.
There exist several pathways in the mRNA decay process, which is catalyzed by various enzymes [66]. To quantitatively investigate the effects of the stochastic mRNA decay process on the gene expression variability, one needs to investigate how the survival probability of mRNA and the cellular levels of mRNA or protein depend on the mean and variance of the key enzyme that catalyzes the rate-determining step of the mRNA decay process.
It is possible to extend the present formulation to the case where the mRNA decay process is a general non-Poisson stochastic process. For the general model of the mRNA decay process, we obtain the exact analytic expressions for the protein noise, which will be reported elsewhere.