Homeostasis of protein and mRNA concentrations in growing cells

Many experiments show that the numbers of mRNA and protein are proportional to the cell volume in growing cells. However, models of stochastic gene expression often assume constant transcription rate per gene and constant translation rate per mRNA, which are incompatible with these experiments. Here, we construct a minimal gene expression model to fill this gap. Assuming ribosomes and RNA polymerases are limiting in gene expression, we show that the numbers of proteins and mRNAs both grow exponentially during the cell cycle and that the concentrations of all mRNAs and proteins achieve cellular homeostasis; the competition between genes for the RNA polymerases makes the transcription rate independent of the genome number. Furthermore, by extending the model to situations in which DNA (mRNA) can be saturated by RNA polymerases (ribosomes) and becomes limiting, we predict a transition from exponential to linear growth of cell volume as the protein-to-DNA ratio increases.

D espite the noisy nature of gene expression [1][2][3][4][5][6] , various aspects of single cell dynamics, such as volume growth, are effectively deterministic. Recent single-cell measurements show that the growth of cell volume is often exponential. These include bacteria 7-10 , archaea 11 , budding yeast 10,12-15 and mammalian cells 10,16 . Moreover, the mRNA and protein numbers are often proportional to the cell volume throughout the cell cycle: the homeostasis of mRNA concentration and protein concentration is maintained in an exponentially growing cell volume with variable genome copy number [17][18][19][20][21][22] . The exponential growths of mRNA and protein number indicate dynamical transcription and translation rates proportional to the cell volume, rather than the genome copy number. However, current gene expression models often assume constant transcription rate per gene and constant translation rate per mRNA (constant rate model) 1,5,[23][24][25] . Assuming a finite degradation rate of mRNAs and non-degradable proteins, these models lead to a constant mRNA number proportional to the gene copy number and linear growth of protein number [26][27][28] , incompatible with the proportionality of mRNA and protein number to the exponentially growing cell volume.
Since the cell volume, protein copy number and mRNA copy number grow exponentially throughout the cell cycle, one may expect a sufficient condition to achieve a constant concentration is to let them grow with the same exponential growth rate. However, mathematical analysis suggests this is insufficient. Let us consider the logarithm of protein concentration c, which can be written as ln(c) = ln(p) − ln(V). Here p is the protein number and V is the cell volume. If one assumes the protein number and the cell volume grow exponentially but independently, with timedependent exponential growth rates λ p (t) and λ v (t) respectively, the time derivative of the logarithm of concentration then obeys d ln(c)/dt~λ p (t) − λ v (t). Even when the time-averaged growth rates of protein number and cell volume are equal, hλ p ðtÞi ¼ hλ v ðtÞi, any fluctuations in the difference between them will accumulate and lead to a random walk behavior of the logarithm of concentration. The homeostasis of protein and mRNA concentrations implies that there must be a regulatory mechanism in place to prevent the accumulation of noise over time.
The main goal of this work is to identify such a mechanism by developing a coarse-grained model taking into account cell volume growth explicitly. Specially, we only consider continuously proliferating cells and do not take account of nongrowing cells, e.g., bacterial cells in stationary phase 29 . The ubiquity of homeostasis suggests that the global machinery of gene expression, RNA polymerases (RNAPs) and ribosomes, should play a central role within the model. Based on the assumption that the number of ribosomes is the limiting factor in translation, we find that the exponential growth of cell volume and protein number originates from the auto-catalytic nature of ribosomes [30][31][32][33] . The fact that ribosomes make all proteins ensures that the protein concentrations do not diverge. Based on the assumption that the number of RNAP is the limiting factor in transcription, we find that the mRNA number also grows exponentially and the mRNA concentration is independent of the genome copy number because of the competition between genes for this global resource [18][19][20] . We also study the effects of genome replication. Due to the heterogeneous timing of gene replication, the transcription rate of one gene has a cell cycle dependence. Within our model, it doubles immediately after the gene is replicated and decreases gradually as other genes are replicated. Nevertheless, we find that this leads to a small effect on protein levels. Finally, we extend our model to more general situations in which an excess of RNAP (ribosome) leads to the saturation of DNA (mRNA). We propose a phase diagram of gene expression and cellular growth controlled by the protein-to-DNA ratio. We predict a transition from exponential growth to linear growth of cell volume as the protein-to-DNA ratio passes a threshold.

Results
Model of stochastic gene expression. In constant rate models, the transcription rate per gene and the translation rate per mRNA are constant 1,5,24 (Fig. 1a). Constant rate models predict a constant mRNA number proportional to the gene copy number and independent of the cell volume. However, experimental observations on plant and mammalian cells have revealed a proportionality between mRNA number and cell volume for cells with a constant genome copy number [18][19][20] Fig. 1 The growing cell model of stochastic gene expression in comparison with constant rate models. a In the constant rate model, the transcription rate is proportional to the gene copy number, and the translation rate is proportional to the mRNA number. These assumptions imply that the gene number and mRNA number are the limiting factors in gene expression. b In Phase 1 of the growing cell model, we introduce as limiting factors RNA polymerases (RNAPs) and ribosomes. Genes with different colors are transcribed with different rates. Here k 0 is a constant and the gene regulation is coarse-grained into the gene allocation fraction ϕ i ¼ g i = P j g j . g i is the effective copy number of gene i (also accounting for the promoter strength). n is the total number of RNAPs. Translation rates of mRNA depend on the number of active ribosomes (f a r), the translation rate k t , and the fraction of mRNA i in the total pool of mRNA. In a later section (A unified phase diagram of gene expression and cellular growth), we will relax our assumptions and consider situations in which the limiting factors of gene expression become the gene number and the mRNA number comparing the cells before and after the genome replication (S phase), the proportionality coefficient between mRNA and cell volume does not exhibit any obvious change. In contrast, a constant transcription rate per gene would predict a doubled transcription rate after the replication of the whole genome, leading to a higher mRNA concentration. In one class of constant rate models 26,27,34 , a deterministic exponential growth of cell volume is explicitly considered. The resulting perturbation on the concentrations due to genome replication is suppressed in the long lifetime limit, but still significant for short lifetime molecules, e.g., mRNA (see Fig. 1 in ref. 27 ).
Considering translation, various experiments have shown that the number of ribosomes is the limiting factor rather than the number of mRNAs. The most direct evidence is the growth law: the growth rate of cells is proportional to the fraction of ribosomal proteins in the total proteome (with a constant factor depending on the growth condition) 35 both for bacterial cells 30,31,36 and budding yeast cells 32 . This means a constant fraction of ribosomes are actively translating mRNAs. These results suggest that in general cells are below the saturation limit in which there are too many ribosomes that the mRNAs can bind. We will therefore assume the biological situation in which mRNAs in the cell compete for the limiting resource of actively translating ribosomes, therefore the translation rate of one type of mRNA is proportional to the number of active ribosomes times its fraction in the total pool of mRNAs.
Considering transcription, experiments have shown that mutants of fission yeasts altered in cell size regulated global transcription to maintain similar transcription rates per cell volume regardless of the cellular DNA content. The changes in total transcription correlated with coordinated changes in gene occupancy by RNA polymerases 37 . These results suggest that the number of RNAPs may be the limiting factor in transcription rather than the gene number, and similar evidence has been shown for bacterial cells 38 and mammalian cells 39 . However, in the same experiments on fission yeast 37 , it has also been found that in cell-cycle-arrested mutants, total transcription rates stopped increasing as the cell volume exceeded a certain value, which suggested DNA became limiting for transcription at low DNA concentration. This result suggests that an excess of RNAPs may lead the gene number to become the limiting factor in certain conditions. In this section, we will focus on the scenario that both RNAP and ribosome are limiting in gene expression, which we denote as Phase 1. In this phase, we will show that the mRNA number and the protein number are proportional to the cell volume and grow exponentially. In a later section (A unified phase diagram of gene expression and cellular growth), we will consider a more general model in which the limiting nature of RNAPs and ribosomes may break down and the dynamics of mRNA and protein number is different.
To address the limiting nature of RNAP, we define an effective gene copy number g i for each gene to account for its copy number and the binding strength of its promoter, which determines its ability to compete for RNAPs. The transcription rate for one specific gene i is proportional to the fraction of RNAPs that are working on its gene(s), ϕ i ¼ g i = P j g j , which we denote as the gene allocation fraction. Gene regulation is thus coarse-grained into the gene allocation fraction ϕ i . The transcription rate is independent of the genome copy number since a change in the genome number leaves the allocation fraction of one gene invariant, a conclusion which is consistent with a number of experimental results on various organisms [18][19][20]37 .
In fact, explicit gene regulation can also be included in our model (Methods), with a time-dependent g i . In such scenarios, g i may be a function of protein concentrations (for instance, the action of transcription factors modifies the transcription rate). Such models will lead to more complex dynamics of mRNA and protein concentrations. However, since we are interested in the global behavior of gene expression and cell volume growth, we do not focus on these complex regulations in this manuscript. Our conclusions regarding the exponential growth of mRNA and protein number for constitutively expressed genes and the exponential growth of cell volume on the global level are not affected by the dynamics of gene expression of particular genes.
In the following, m, p, r, n represents the numbers of mRNAs, proteins, ribosomes and RNA polymerases, respectively. Proteins (p) also include RNAPs (n) and ribosomes (r) 30 . We consider the degradation of mRNA with degradation time τ for all genes. The protein number decreases only through cell divisions (though adding a finite degradation rate for proteins does not affect our results). The stochastic dynamics of gene expression within Phase 1 of our model are summarized in the following sets of equations and Fig. 1b, Here k 0 , k t are constants, characterizing the transcription (translation) rate of a single RNAP (ribosome). f a is the fraction of active ribosomes, which we assume to be constant in a given nutrient environment 30,32 . We note that nonspecifically bound RNAPs have been reported in bacteria 40,41 . We will discuss their effect later. For simplicity, we first assume the values of ϕ i do not change in time. This can be formally thought of as corresponding to an instantaneous replication of the genome. In reality, a finite duration of DNA replication and the varying time of replication initiation for different genes lead to ϕ i 's that change during the DNA replication. We later analyze a more complete version of the model which includes these gene dosage effects, but we first consider the simplified scenario of constant ϕ i that will capture the essential features of the problem. We assume the cell volume is approximately proportional to the total protein mass, i.e., V / M ¼ P j p j , which is a good approximation for bacteria 42,43 and mammalian cells 17 . To simplify the following formulas, we consider each protein has the same mass and set the cell density as 1.
Due to the fast degradation of mRNA compared with the cell cycle duration 44,45 , the mRNA number can be well approximated as being in steady state. We can express the ensemble-averaged number of mRNA of gene i as Equation (3) then leads to the time-dependence of average ribosome number, d r h i=dt ¼ k t f a ϕ r r h i, reproducing the autocatalytic nature of ribosome production and the growth rate determined by the relative abundance of active ribosomes in the proteome 30,32 . Similarly, the number of protein i grows as d p i h i=dt ¼ k t f a ϕ i r h i. As the cell grows and divides, the dynamics becomes insensitive to the initial conditions, so the protein number will grow exponentially as well 21 . The ratio between the averages of two protein numbers in the steady state is set by the ratio of their production rate, therefore hp i i=hp j i ¼ ϕ i =ϕ j . The average number of mRNA traces the number of RNA polymerases according to Eq. (4), and therefore also grows exponentially. Throughout the cell cycle we have where m b (i) (p b (i)) is the number of mRNA (protein) of gene i at cell birth.
We denote the concentrations of mRNA and protein of gene i as c m i = m i /V and c i = p i /V respectively. According to Eqs. (1)-(3), the deterministic equations of the above variables become (see details in Methods) A fixed point exists for the dynamics of c i and c m i , namely This fixed point is stable due to the global nature of RNAPs and ribosomes: any noises arising from the copy number of RNAPs (ribosomes) equally affect all mRNAs (proteins), and therefore leave the relative fraction of one type of mRNA (protein) in the total pool of mRNAs (proteins) invariant. The average concentrations of mRNA and protein of The results are independent of the cell volume and genome copy number agreeing with experimental data on various organisms [18][19][20]22 .
We take cell division explicitly into account and, for concreteness, use the "adder" model for cell division by considering an initiator protein I. The initiator protein accumulates from cell birth, triggers the cell division once it reaches the division threshold I c and is then destroyed (or "reset", e.g., after initiation of DNA replication in bacteria, the ATP-bound DnaA is dephosphorylated to the ADP-bound form) [46][47][48] . During a division event, we assume proteins and mRNAs are divided between the two daughter cells following a binomial distribution 49 . The initiator protein sets the scale of absolute protein number, and the average number of proteins produced in one cell cycle is equal to Δ(i) = I c ϕ i /ϕ I 47 . Since the protein number grows twofold during one cell cycle, the average protein number of gene i at cell birth is p b (i) = I c ϕ i /ϕ I and the corresponding average mRNA number at cell birth is m b (i) = k 0 I c τϕ i ϕ n /ϕ I . We remark that the exact molecular mechanism of cell division does not affect our results.
We corroborate the above analytical calculations with numerical simulations. These will also capture the stochastic fluctuations in gene expression levels, which are not included in the previous analysis. Due to the short lifetime of mRNAs, the production of proteins can be approximated by instantaneous bursts 24 . We introduce the burst size parameter b 0 as the average number of proteins made per burst, b 0 = k t f a hrðtÞi=h P j m j i τ ≈ k t f a ϕ r = k 0 ϕ n À Á , independent of the cell volume. ϕ i for N = 200 proteins are uniformly sampled in logarithmic space, with the sum over ϕ i (including ribosome and RNAP) constraint to be precisely one. We choose the parameters to be biologically relevant for bacteria: the doubling time T = ln(2)/μ = 150 min, r b = 10 4 , n b = 10 3 , b 0 = 0.8, I c = 20, ϕ r = 0.2, f a = 0.7 and τ = 3.5 min, see other numerical details in Methods. Our conclusions are independent of the specific choice of parameters.
In Fig. 2a, we show the typical trajectories from our simulations of cell volume, protein number and mRNA number for the same gene over multiple generations. To verify the exponential growth of protein and mRNA, we average the protein and mRNA numbers given a fixed relative phase in the cell cycle progression, which is normalized by the generation time and changes from 0 to 1. The averaged values of protein and mRNA numbers (circles) are well predicted by exponential growth, Eqs. (6a) and (6b) (black lines) without any fitting parameters, as shown in Fig. 2b with 3 single trajectories in the background. We also simulate a regulated gene with a time-dependent gene copy number and obtain qualitatively similar results (Methods, Supplementary Fig. 1).
The corresponding trajectories of protein and mRNA concentrations are shown in Fig. 2c, with bounded fluctuations around the predicted averaged values (black lines). In contrast, if the protein number and cell volume grow exponentially but independently, the ratio between them will diverge as the effects of noise accumulate, exhibiting a random walk behavior (Fig. 2d). Considering the cell cycle dependence of mRNA number and the homeostasis of protein concentration throughout the cell cycle, the experimental observation in Escherichia coli showing negligible correlations between mRNA number and protein concentration 50 is consistent with our model, and not contradictory to the strong correlation of mRNA concentration and protein concentration 51 .
Within our model, we may also study the protein number dynamics: how does the protein number at cell division correlate with that at cell birth? We find that the correlations follow an "adder" (i.e. the number of new proteins added is uncorrelated with the number at birth), as shown in Fig. 2e. While this has been quantified in various organisms with respect to cell volume 8,9,11,[52][53][54] , checking correlations between protein content at cell birth and division has received significantly less attention 55,56 . Related to this, we study the auto-correlation function of protein concentration in time. We find that the autocorrelation function is approximately exponential, with a correlation time bounded from below by the doubling time ( Supplementary Fig. 2). Both of these results provide experimentally testable predictions.
Effects of finite duration of gene replication. So far, we considered a constant ϕ i throughout the cell cycle assuming an instantaneous replication of the genome. In this section, we relax this condition and study the effects of finite DNA replication time. We consider the bacterial model of DNA replication, specifically, E. coli, for which the mechanism of DNA replication is well characterized 57 . The duration of DNA replication is constant, and defined as the C period. The corresponding cell division follows after an approximately constant duration known as the D period. Details of the DNA replication model are in the Methods. In Fig. 3a, we show the time trajectories of the gene allocation fraction, mRNA concentration and protein concentration of one gene for a doubling time of T = 30 min with C + D = 70 min. The DNA replication introduces a cell cycle dependent modulation of ϕ i . The abrupt increase of ϕ i corresponds to the replication of the specific gene i (Fig. 3a) However, as other genes are replicated, the relative fraction of gene i in the total genome decreases. This modulation propagates to the mRNA concentration which essentially tracks the dynamics of ϕ i due to its short lifetime. The modulation of mRNA concentration affects the protein concentration as well, yet with a much smaller amplitude. These results can be tested experimentally by monitoring the DNA replication process and mRNA concentration simultaneously. ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-06714-z Noise in gene expression can be classified as intrinsic and extrinsic noise 58 . While intrinsic noise is due to the stochastic nature of the chemical reactions involved in gene expression, extrinsic noise is believed to be due to the fluctuations of external conditions and common to a subset of proteins. Experiments have revealed a global extrinsic noise that affects all protein concentrations in the genome 50,59,60 . Because all genes are subjected to the finite duration of DNA replication, it is tempting to attribute the finite duration of DNA replication as one of the main sources of global extrinsic noise 34 . Within our model in the previous section (constant ϕ i 's throughout the cell cycle), there is no global extrinsic noise ( Supplementary Fig. 3). A global extrinsic noise may emerge after we introduce the timedependent ϕ i due to DNA replication. However, we find that the coefficient of variation (CV, the ratio between standard deviation and mean) of the most highly expressed proteins is only about 0.02 within the growing cell model (Fig. 3b), much smaller than that found in experiments 50,59 . We note that a small extrinsic noise due to gene replication is also observed in constant rate models 26,27 . Moreover, recent experiments and modeling have suggested that a significant part of the extrinsic noise of mRNA expression level can be attributed to the fluctuations of RNAP copy number 28 . Within our model, RNAP level fluctuations will lead to extrinsic noise in mRNA concentrations.
A unified phase diagram of gene expression and cellular growth. Experimental observations on E. coli 30 and budding  [18][19][20] and fission yeast 37 are also consistent with our assumption that RNA polymerase is limiting for transcription. However, as we discussed in the first section, in the same experiments on fission yeast 37 DNA became limiting for transcription at low DNA concentration. Therefore, we cannot exclude the possibility that in some cases because RNAPs are too abundant, DNA becomes the limiting resource for transcription rather than the number of RNAPs. Similarly, when ribosomes are too abundant relative to the transcript number, the limiting factor for translation becomes the transcript number rather than ribosome number.
In this section, we generalize our model by assuming that each gene has an upper bound on the number of RNAPs (n s ) than can simultaneously work on it. A possible extreme case is that the gene is fully loaded with RNAPs, on which RNAPs are only constrained by steric hindrance. The same assumption is made for mRNA with an upper bound of ribosomes (r s ) that can work on it simultaneously. We remark that the exact mechanism of DNA and mRNA saturation is beyond our coarse-grained model. If the number of RNAP (ribosome) is above the upper bound, the transcription (translation) rate is limited by the gene (mRNA) number, in a similar fashion to the constant rate models.
We define the protein-to-DNA ratio (PTD ratio) as the sum of protein numbers divided by the sum of effective gene numbers, As the PTD ratio becomes larger, e.g., due to a sufficiently large cell volume with a fixed number of gene, the number of RNAPs (ribosomes) will exceed the maximum load the total genes (mRNAs) can hold. We have discussed thoroughly Phase 1 (neither DNA nor mRNA is saturated) earlier and we summarize our predictions on the transition from Phase 1 to other phases in the following. Phase 2: In phase 2, the limiting factor in transcription becomes the gene copy number and the transcription rate is proportional to the gene copy number (Fig. 4b). The threshold PTD ratio for the transition from Phase 1 to Phase 2 is  Fig. 4 Phases of gene expression and cell volume growth. a Theoretical phase diagram of gene expression and cellular growth within our model. The x axis is the protein-to-DNA ratio (γ). When γ < γ 1 , neither DNA nor mRNA is saturated. The mRNA number, the protein number and the cell volume all grow exponentially with the growth rate set by the fraction of ribosomal gene in the total genome (ϕ r ). When γ 1 < γ < γ 2 , DNA is saturated but mRNA is not. The protein number and the cell volume still grow exponentially while the mRNA number is a constant proportional to the gene number. When γ > γ 2 , both DNA and mRNA are saturated. The protein number and cell volume grow linearly, and the cell volume growth rate is set by the genome copy number. b The gene expression dynamics in phase 2. In this phase, DNA is saturated by RNAPs, therefore, the transcription rate is proportional to the effective gene copy number, g i . n s is the upper bound of the number of RNAPs that can work on one gene simultaneously. The translation rate is the same as in phase 1. To simplify the formula, we assume all ribosomes are active (to include the effect of an inactive fraction, r should be replaced by f a r). c The gene expression dynamics in phase 3, in which both DNA and mRNA are saturated. The translation rate is proportional to the mRNA number. r s is the upper bound of the number of ribosomes that can work on one mRNA simultaneously Here n s is the upper bound of the number of RNAPs that can work on one gene and ϕ n is the gene allocation fraction of RNAP. Because mRNA is not saturated, the protein number and the cell volume grow exponentially with the same growth rate as Phase 1, Eq. (5), and the homeostasis of protein concentration is still valid. However, because the production rate of mRNA is now proportional to the gene copy number, the mRNA concentration is not constant anymore as the cell volume grows (Methods). In Phase 2, even though the transcription rate doubles after the genome is replicated, the translation rate is proportional to the relative fraction of mRNA in the total pool of mRNAs and therefore still independent of the genome copy number. The average protein concentration is equal to the gene allocation fraction ( c i h i ¼ ϕ i ). Recent proposed theoretical models of gene expression are consistent with this phase 61 . In terms of transcription, our model in Phase 2 is equivalent to constant rate models and we have confirmed that for both bacteria and mammalian cells, the typical lifetime of mRNA is short enough compared with the doubling time to distinguish Phase 1 and Phase 2 ( Supplementary Fig. 4).
Phase 3: As the cell volume becomes larger, mRNA may get saturated as well. The limiting factor for translation is now the mRNA copy number (Fig. 4c). The threshold PTD ratio for the transition from Phase 2 to Phase 3 is (Methods) Here r s is the upper bound of the number of ribosomes that can work on one mRNA. In this phase, the translation rate is proportional to the mRNA number and the protein number grows linearly as _ p i ¼ k t k 0 g i τn s r s , with a linear growth rate proportional to the gene number. Therefore, within the assumption that the cell volume is determined by the total protein number, the cell volume grows linearly as well with the linear growth rate proportional to the total gene number, and therefore proportional to the genome copy number. As in Phase 2, the mRNA concentration decreases as the cell volume grows, however, the protein concentration is still constant with the average protein concentration equal to the gene allocation fraction ( c i h i ¼ ϕ i , Methods). In Phase 3, even though the cell volume grows linearly, the population still grows exponentially with a population growth rate. However, there is no general relation between the ribosomal fraction in the proteome and the population growth rate, in contrast to the growth law in Phase 1 and 2. We summarize the predicted phase diagram of cellular growth in Fig. 4a.
To gain some sense regarding the parameters associated with our proposed phase diagram, we estimate the PTD ratio of E. coli. Considering the typical cell volume of E. coli as 1 μm 3 , the protein density as 3 × 10 6 proteins/μm 3 and the total number of proteincoding genes in E. coli as 4000 62 , we estimate the protein-to-DNA ratio for E. coli as γ~1000. Estimates of the two threshold values of PTD ratios (see Methods) suggest that γ 1~1 500 and γ 22 0,000. These estimates suggest that wild-type E. coli cells are found in Phase 1, but close to Phase 2. We remark that the actual threshold values of PTD ratio for the transitions between different growth phases may be affected by other factors, e.g., the heterogeneous size of genes, but we propose that the general scenario of the transition from Phase 1 to Phase 3 as the protein-to-DNA ratio increases should be generally applicable. As the PTD ratio increases, we predict a transition from exponential growth to linear growth for protein number and cell volume (Supplementary Fig. 5). We propose future experiments to study the potential transition from exponential to linear growth of cell volume, for example using filamentous E. coli where cell division and gene replication are inhibited. Similar experiments can also be done for larger cells, e.g., mammalian cells, in which the transition from exponential growth to linear growth of cell volume may be easier to achieve. Preliminary results from experiments measuring the growth of cell mass of mammalian cells 63 and yeast cells 64 indeed show a crossover from exponential growth to linear growth when the cell mass is above a threshold value, consistent with our prediction.
It has been shown in bacteria that there are excess RNAPs nonspecifically bound to DNA 40,41 . In the Methods, we consider a modified model taking into account the partitioning of RNAPs to free RNAPs, elongating RNAPs, promoter-bound RNAPs and nonspecifically bound RNAPs. The transcription rate is determined by the concentration of free RNAPs through Michaelis-Menten kinetics 40,65 . We find that our conclusions remain intact with an approximately constant fraction of actively transcribing RNAPs in the total RNAPs for Phase 1 (Supplementary Fig. 6). The effect of nonspecifically bound RNAPs is therefore to renormalize the transcription constant k 0 in Phase 1 (Eq. (1)) by a constant factor. The transition from Phase 1 to Phase 2 is qualitatively unaffected (Supplementary Fig. 7) and the threshold PTD ratio γ 1 (Eq. (9)) from Phase 1 to Phase 2 is changed by a constant factor (Methods). We note that alternative mechanisms of gene saturation can occur upon introducing the different classes of RNAPs, through the saturation of free RNAPs and the Michaelis-Menten kinetics (Methods).

Discussion
In this work, we propose a coarse-grained model of stochastic gene expression incorporating cell volume growth and cell division. In the first part, we consider the biological scenario that RNAPs are limiting for transcription and ribosomes are limiting for translation. In other words, neither DNA nor mRNA is saturated. We find that the limiting nature of ribosomes in the translation process leads to the exponential growth of protein numbers. The limiting nature of RNA polymerase and its exponential growth lead to the exponential growth of mRNA numbers. Homeostasis of protein concentrations originates from the fact that ribosomes make all proteins. Homeostasis of mRNA concentration comes from the resulting bounded concentration of RNAPs. Our model is consistent with the constancy of mRNA and protein concentration as the genome copy number varies since the transcription rate depends on the relative fraction of genes in the genome rather than its absolute number 22 .
During DNA replication, we find that the gene allocation fraction ϕ i for one specific gene doubles after the gene is replicated but decreases afterwards since other genes are replicated as well and compete for RNAPs. This prediction can be tested by monitoring the mRNA concentration and the copy number of one gene throughout the cell cycle. Furthermore, we extend our model to more general cases in which DNA and mRNA can be saturated by an excess of RNAP and ribosome. We find three possible phases of cellular growth as the protein-to-DNA ratio γ increases. A transition from exponential growth to linear growth of protein number and cell volume is predicted. In the future, it will be interesting to study the interplay between the global interactions which are the focus of this work and local interactions between genes. Our model provides an alternative model to constant rate models to study genetic networks, which would be advantageous when cell cycle effects are important.

Methods
Derivation of protein and mRNA concentrations. We define the fraction of mRNA i in the total mRNA pool as f i ¼ m i = P j m j , and the concentration of mRNA and protein of gene i as c m i = m i /V, c i = p i /V. We denote the RNAP and ribosome concentrations as c n and c r , respectively. According to Eqs. (1)-(3), the deterministic equations of the above variables then become Using the condition that mRNA degradation time is much smaller than the doubling time (μτ ( 1), we find the fixed points for the dynamics of f i , c i , and c m i . These are, f i ¼ c i ¼ ϕ i and c m i ¼ k 0 ϕ i c n τ. Replacing f i by ϕ i and c n by ϕ n , we obtain the approximate version of the above equations, Eqs. (7a) and (7b).
Simulations of independent growth model. In the growth model corresponding to Fig. 2d, we assume the protein number and cell volume grow exponentially and independently, Here, ξ p (t), ξ V (t) are white noise terms, with the auto-correlation function, ξ p;V ð0Þξ p;V ðtÞ D E = A p;V δðtÞ. In Fig. 2d of the main text, we choose A p = A V = 1.
Simulations of growing cell model. We simulated Eqs. (1)-(3), fixing r b , n b , b 0 , ϕ r , f a , I c , τ as well as the growth rate μ. Other parameters are inferred given the above values, e.g., ϕ n = n b ϕ r /r 0 , k t = μ/(ϕ r f a ), k 0 = k t f a r b /(b 0 n b ). We fix the time step δt so that the probability for one event to happen during a time step is smaller than 0.1. We track one of the daughter cells after cell division.
Gene dosage effects. In reality, the gene allocation fraction ϕ i changes during the cell cycle due to the finite duration of DNA replication. In this section we introduce the modified version of the gene expression model incorporating DNA replication. Although our model is general, we focus on DNA replication in bacteria for concreteness, specifically E. coli where this process is very well characterized. We expect our conclusions to be generally valid. Furthermore, we refine our model for cell division, assuming that the initiator protein triggers the initiation of DNA replication rather than cell division, with the threshold I c proportional to the number of origins of replication 57,66 (the number of which doubles at each initiation). We assume that the cell division takes place a fixed time C + D after initiation of the DNA replication, where C, D are respectively the time for DNA replication and the time between the completion of DNA replication and cell division. The number of origins reduce by half at each cell division. Other details are the same as in the main text. Each gene doubles its copy number during the C period, and we choose this gene replication time to be randomly and uniformly distributed across all genes. When a gene i replicates, where the second equation accounts for the normalization of the gene allocation fraction. We choose the experimentally reported C and D and cell doubling time from ref. 57 . In Fig. 3a, we simulate the model by tracking one daughter cells. In Fig. 3b, we track all the cells in an exponentially growing population, which starts from 100 cells to 5000 cells.
Simulations of gene activation. We generalize the constitutive expressed genes considered in the main text to include a single regulated gene by considering a random telegraph process of the effective gene copy number 1 , Here the gene deactivation rate k À g is constant, and the activation rate is set by the concentration of transcription factor through positive regulation, k þ g ¼ k g0 c TF . Here, k g0 is constant. When gene i is active, the corresponding gene allocation fraction follows ϕ i ¼ g i0 = P j g j , and when it becomes deactivated ϕ i = 0. Note that here we only consider one regulated gene i, but the changing gene allocation of gene i also affects other genes' allocation fraction. We simulate the model in Phase 1, and the deactivation of gene i increases other genes' allocation fraction as ϕ j → ϕ j /(1 − ϕ i ), with ϕ i ¼ g i0 = P j g j . Simulated trajectories of gene allocation fraction, mRNA number, protein number and cell volume are shown in Supplementary Fig. 1.
General model of gene expression. We consider the generalized equation of mRNA number, Eq. (1) in the deterministic limit as Here n c is the threshold number of RNAPs above which DNA starts to be saturated, in which case the transcription rate becomes proportional to the effective gene copy number g i and independent of the RNAP number. For one gene, the maximum load of RNAP that it can hold is g i n s , where n s is the maximum number of RNAPs that a single copy of constitutively expressed gene (g i = 1) can hold. n c can be computed as We also generalize the growth of protein number from Eq. (3) to Here r c is the maximum number of ribosomes above which mRNA starts to be saturated. We drop the fraction of actively working ribosomes since it is often a constant depending on the growth condition 30 . r s is the maximum number of ribosomes one mRNA can hold. We can calculate r c as Given Eqs. (20) and (22), we obtain four possible phases: (i) n < n c , r < r c , (ii) n > n c , r < r c , (iii) n > n c , r > r c , and (iv) n < n c , r > r c . Given a fixed value of ϕ r and ϕ n , either (ii) or (iv) is possible. Realization of (ii) requires that n> P i g i n s and r<k 0 τr s P i g i n s , therefore ϕ n ϕ r > 1 k 0 τr s : ð24Þ In cases where Eq. (24) breaks down, a finite fraction of ribosomes are not utilized. This requires a large fraction of genes in the genome making ribosomes that cannot translate because mRNAs are saturated. Since ribosomes are typically more expensive to make than other proteins 30,33 , we assume the biological scenario, Eq. (24) will be satisfied. From Eq. (21) and using n= P i p i ¼ ϕ n , we obtain the threshold PTD ratio for the transition from Phase 1 to Phase 2, In Phase 2, the average mRNA concentration becomes which is inversely proportional to the protein-to-DNA ratio. From Eq. (23) and using r= P i p i ¼ ϕ r , we obtain the transition PTD ratio from Phase 2 to Phase 3 as, In Phase 3, the mRNA concentration is the same as Phase 2. Because the cell volume is the sum of all proteins, the protein concentration is the same as Phase 2 and Phase 1, Estimation of the threshold protein-to-DNA ratios for E. coli. We approximate the upper bound of RNAP number working on a single gene as roughly equal to the number of RNAPs on a typical gene (~10 3 base pairs) when half of the gene is occupied. The linear size of RNAP is about 5 nm, and the length of one base pair is about 0.3 nm, leading to the estimate n s~3 0. A similar calculation for the upper bound of ribosome on a single mRNA leads to r s~1 0 since ribosome's linear size is about 3 times larger than RNAP 62 . We take ϕ r ≈ 0.2 according to the ref. 30 , and estimate the gene allocation fraction of RNAP to be ϕ n~0 .02 since the number of RNAPs in E. coli is roughly 10% of the number of ribosomes 62 . We estimate the life time of mRNA as 5 min 62 .
We estimate the transcription rate of one RNAP by considering two potential limiting steps in transcription and taking the slower one. First, assuming the initiation of transcription is diffusion limited, we could estimate the time scale for one RNAP to bind the transcription site as Δt~1 μm 2 /(0.2 μm 2 /s)~5 s using the measured diffusion constant of RNAP 41,67 . Second, we could also estimate the elongation time as the typical length of gene divided by the elongation rate of RNAP, Δt~1000 nt/50(nt/s)~20 s 62 . Taking the slower time scale from the above two calculations, we estimate k 0 ≈ 0.05 s −1 . Finally, we compute γ 1 and γ 2 using the above estimated parameters, and obtain γ 1~1 500, γ 2~2 0,000.
Effect of nonspecifically bound RNAPs. Previous studies on bacteria have shown that there are excess RNAPs bound nonspecifically to the genome and modeled their kinetics 40,41 . In this section, we consider a modified model to take into account nonspecifically bound RNAPs. For our purpose, we consider a simplified model with four classes of RNAPs, namely, (i) elongating RNAPs, n e (ii) RNAPs bound to a promoter, n p (iii) RNAPs nonspecifically bound to DNA, n ns (iv) free RNAPs, n free . We assume a Michaelis-Menten relation for the number of promoterbound RNAPs and nonspecifically bound RNAPs 40,65 , n ns ¼ Gg ns c free Here c free is the concentration of free RNAPs and K s , K ns are the Michaelis constants. G is the total number of genes and g ns is the number of nonspecific binding sites per gene. Note that c free = c n F free , with c n the concentration of total RNAPs and F free the fraction of free RNAPs in the total RNAP pool. For simplicity, we assume one promoter for each gene. The number of elongating RNAPs is then given by: where k cat is the transition rate from promoter-bound RNAPs to elongating RNAPs, τ 0 is the time for transcribing a gene, and Λ = k cat τ 0 is the ratio between the number of elongating RNAPs and promoter-bound RNAPs 40 . We consider parameter regimes motivated by typical biological scenarios, in particular (1) the number of sites for nonspecific binding is much larger than the number of promoters, (2) nonspecific binding of RNAPs to DNA is much weaker than the specific binding of RNAPs to promoters, (3) the number of promoter-bound RNAPs is small compared with elongating RNAPs. Using n = n e + n p + n ns + n free , we can find a self-consistent equation for F free , Here γ is the protein-to-DNA (PTD) ratio, i.e., the ratio between the total number of proteins (equivalent to cell volume V within our model) and the total gene number, γ = V/G. We can use Eq. (31) to compute the fraction of free RNAPs given a PTD ratio and use Eqs. (28)- (30) to compute the fraction of elongating RNAPs, F e and nonspecifically bound RNAPs, F ns .
Since the left side of Eq. (31) monotonically decreases from c n γ to 0 as F free increases from 0 to 1 and the right side of Eq. (31) monotonically increases as F free increases, we find that as the protein-to-DNA ratio increases, F free increases. Previous studies have shown that F free ≈ 0.1 for bacteria 40,41 and support the assumption that c free is smaller or comparable to K s 40,65 (note that nonspecific binding is characterized by K ns larger than K s ). In the following, we first assume a small F free and c free ( K s , and show this leads to a behavior qualitatively equivalent to Phase 1 of the main text, albeit with a renormalization of the transcription constant k 0 . We also later discuss the situations when these conditions break down, and show that they lead to behavior consistent with Phase 2. The transcription rate for one specific gene with an effective gene copy number g i can be written as Here g i /G = ϕ i is the gene allocation fraction of gene i, and F e is the fraction of elongating RNAPs in the total RNAP pool. In the limit of a small F free and c free ( K s , F e becomes a constant: F e % n e n e þ n p þ n ns ¼ Λ 1 þ Λ þ g ns K s =K ns : Therefore, the transcription rate corresponds to Phase 1 of our model, TR i 1 k 0 ϕ i n withk 0 ¼ k 0 F e . In the limit that all RNAPs are actively transcribing and no nonspecifically bound RNAPs, Λ ) 1 and g ns = 0, we go back to the situation that all RNAPS are working withk 0 ¼ 1=τ 0 ¼ k 0 . Therefore we conclude that the introduction of nonspecifically bound RNAPs does not affect our model qualitatively, and its effect is to renormalize the transcription constant k 0 in Eq. (1) by a constant factor, We simulate a single lineage of growing cells using the full model (with partitioning of RNAPs and gene replication). We set the parameters as: Λ = 50, g ns = 1000, k 0 = 1 min −1 , G = 2000 (before gene duplication), K s = 0.02〈c〉, K ns = 0.8 〈c〉, where 〈c〉 is the total protein concentration within the cell which we set to 1 throughout the paper. Our conclusions are independent of the specific values of parameters. The gene allocation fractions are the same as the main text and the average RNAP concentration during the cell cycle 〈c n 〉 = ϕ n = 0.02. The fractions of elongating RNAPs and nonspecifically bound RNAPs are approximately constant with a small cell cycle modulation (the coefficient of variation is of the order of 0.01), consistent with the above results since F free is small ( Supplementary  Fig. 6a). We also find a linear scaling between mRNA number and cell volume, consistent with Phase 1 of our model (Supplementary Fig. 6b).
We next consider the transition from Phase 1 (RNAP limiting) to Phase 2 (gene limiting). Assuming the saturation of genes is due to the steric hindrance of elongating RNAPs with a minimum distance between two successive RNAPs, we can find the threshold PTD ratio from Phase 1 to Phase 2. Since the fraction of elongating RNAPs is constant, we can compute the threshold PTD ratio using nF e /G = n s (where n s is the upper bound of the number of RNAPs that can work on one gene simultaneously), and obtaiñ We find that the introduction of nonspecifically bound RNAPs does not affect our main results qualitatively and it changes the threshold PTD ratio from Phase 1 to Phase 2, Eq. (9), by a constant factor, which is the inverse of the fraction of elongating RNAPs. In Supplementary Fig. 7, we show the transition from Phase 1 to Phase 2 in the proposed experiment in which cell division is blocked, in the presence of nonspecific binding. At each point in time we solve the self-consistent equation of the partitioning of RNAPs, Eq. (31). The plot shows both the deterministic dynamics of average mRNA number (black lines) as well as the results of the stochastic simulations (red/blue solid lines). We compare the dynamics for two cells with different fixed genome sizes. Initially, the two cells are in Phase 1 and therefore have identical mRNA number proportional to cell volume. When they exceed their respective threshold PTD ratio, the mRNA number begins to saturate and becomes limited by the gene number, therefore the cell with a twice larger genome size has twice more mRNAs (Phase 2).
So far we assumed that F free is small and c free is smaller than the Michaelis constant. In the following we relax these assumptions. We find that the introduction of the different RNAP classes introduces alternative mechanisms for gene saturation, that lead to behavior consistent with Phase 2.
We first relax the assumption that F free ( 1 (assuming that c free ( K s still holds). Because the transcription rate for a specific gene is when F free is comparable to 1, the transcription rate will be saturated as well and proportional to the gene copy number g i . From Eq. (31), we find that F free ≈ γ= 1þΛ K s þ g ns K ns þ γ , therefore the threshold PTD ratio for F free to be comparable to 1 is Second, we can relax the assumption that c free ( K s (assuming F free ( 1 still holds). When c free ≈ K s , the transcription rate can be saturated as well. Using F free ≈ γ= 1þΛ K s þ g ns K ns , we can obtain the corresponding threshold PTD ratio as We remark that the particular mechanism which drives the cell to Phase 2 (gene limiting) is the one with the smallest threshold. Comparison between Eqs. (34),