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.

Introduction

Despite the noisy nature of gene expression[1-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-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-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-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 time-dependent exponential growth rates λ(t) and λ(t) respectively, the time derivative of the logarithm of concentration then obeys d ln(c)/dt ~ λ(t) − λ(t). Even when the time-averaged growth rates of protein number and cell volume are equal, , 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 non-growing 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-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-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-20]. Moreover, even 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]).

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 k0 is a constant and the gene regulation is coarse-grained into the gene allocation fraction . g 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 (fr), the translation rate k, 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

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 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), , which we denote as the gene allocation fraction. Gene regulation is thus coarse-grained into the gene allocation fraction ϕ. 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-20,37].

In fact, explicit gene regulation can also be included in our model (Methods), with a time-dependent g. In such scenarios, g 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 k0, k are constants, characterizing the transcription (translation) rate of a single RNAP (ribosome). f 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 ϕ 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 ϕ’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 ϕ that will capture the essential features of the problem.

We assume the cell volume is approximately proportional to the total protein mass, i.e., , 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 asEquation (3) then leads to the time-dependence of average ribosome number, , reproducing the auto-catalytic nature of ribosome production and the growth ratedetermined by the relative abundance of active ribosomes in the proteome[30,32].

Similarly, the number of protein i grows as . 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 . 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 havewhere m(i) (p(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  = m/V and c = p/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 and , namely and . 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 gene i become , and . The results are independent of the cell volume and genome copy number agreeing with experimental data on various organisms[18-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 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-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ϕ/ϕ[47]. Since the protein number grows twofold during one cell cycle, the average protein number of gene i at cell birth is p(i) = Iϕ/ϕ and the corresponding average mRNA number at cell birth is m(i) = k0Iτϕϕ/ϕ. 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 b0 as the average number of proteins made per burst,  =  ≈ , independent of the cell volume. ϕ for N = 200 proteins are uniformly sampled in logarithmic space, with the sum over ϕ (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 = 104, n = 103, b0 = 0.8, I = 20, ϕ = 0.2, f = 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).

Fig. 2

Exponential growth of the cell volume, protein number, mRNA number; the homeostasis of protein and mRNA concentrations throughout the cell cycle. a Numerical simulated trajectories of cell volume, protein number, and mRNA number are shown (ϕ = 0.018). b The averaged values of protein and mRNA numbers of a highly expressed gene (ϕ = 0.04), are shown (circles) with 3 single trajectories in the background. The black lines are theoretical predictions of Eqs. (6a) and (6b). The average is over 130 cell cycles. The color band represents the standard deviation (same for (c)). c The averaged values of protein and mRNA concentrations of the same gene as in (b) are shown (circles). The black lines are theoretical predictions of Eqs. (7a) and (7b). Three trajectories are shown in the background. d Three trajectories of diverging concentrations in the scenario where the protein number and cell volume grow independently. See the numerical details in Methods. e The scatter plot of the protein numbers at cell division (P) v.s. the protein numbers at cell birth (P). The circles are binned data. The black line is a linear fit of the binned data with slope 1.03, consistent with the adder correlations

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-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 auto-correlation 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 ϕ 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 ϕ. The abrupt increase of ϕ corresponds to the replication of the specific gene i (Fig. 3a) ϕ → 2ϕ. 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 ϕ 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.

Fig. 3

Effects of finite duration of DNA replication. a The time trajectory of gene allocation fraction (triangles), mRNA concentration (squares) and protein concentration (circles) of a high copy number protein (μ ≈ 104, see (b)). The doubling time is T = 30 min, and we use the values of the C and D periods from ref.[57], namely, C = 35 min and D = 35 min. In this situation, the cell undergoes DNA replication throughout the cell cycle. Nevertheless, the noise in ϕ does not propagate to the noise in protein concentration significantly. The value of mRNA concentration is 5 times amplified for clarity. b An exponentially growing population is simulated (See Methods). The noise magnitude is quantified as the square of CV of protein concentrations. The mean protein number (μ) is the protein number per average cell volume. Gene dosage effects due to DNA replication do not generate a significant global extrinsic noise. Two different doubling times are considered

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 ϕ’s throughout the cell cycle), there is no global extrinsic noise (Supplementary Fig. 3). A global extrinsic noise may emerge after we introduce the time-dependent ϕ 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 yeast[32] support our assumption that ribosomes are limiting for translation. Experimental observations on plant and mammalian cells[18-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) 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) 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 (Methods),Here n is the upper bound of the number of RNAPs that can work on one gene and  ϕ 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 (). 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).

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 (ϕ). 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. n 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 fr). 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 is the upper bound of the number of ribosomes that can work on one mRNA simultaneously

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 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 , 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 (, 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 μm3, the protein density as 3 × 106proteins/μm3 and the total number of protein-coding 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 ~ 1500 and γ2 ~ 20,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 k0 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 ϕ 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 , and the concentration of mRNA and protein of gene i as  = m/V, c = p/V. We denote the RNAP and ribosome concentrations as c and c, respectively. According to Eqs. (1)–(3), the deterministic equations of the above variables then becomeUsing the condition that mRNA degradation time is much smaller than the doubling time (), we find the fixed points for the dynamics of f, c, and . These are,  and . Replacing f by ϕ and c by ϕ, 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, ξ(t), ξ(t) are white noise terms, with the auto-correlation function,  = . In Fig. 2d of the main text, we choose A = A = 1.

Simulations of growing cell model

We simulated Eqs. (1)–(3), fixing r, n, b0, ϕ, f, I, τ as well as the growth rate μ. Other parameters are inferred given the above values, e.g., ϕ = nϕ/r0, k = μ/(ϕf), k0 = kfr/(b0n). 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 ϕ 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 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 is constant, and the activation rate is set by the concentration of transcription factor through positive regulation, . Here, k is constant. When gene i is active, the corresponding gene allocation fraction follows , and when it becomes deactivated ϕ = 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 ϕ → ϕ/(1 − ϕ), with .

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 asHere n 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 and independent of the RNAP number. For one gene, the maximum load of RNAP that it can hold is gn, where n is the maximum number of RNAPs that a single copy of constitutively expressed gene (g = 1) can hold. n can be computed asWe also generalize the growth of protein number from Eq. (3) toHere r 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 is the maximum number of ribosomes one mRNA can hold. We can calculate r asGiven Eqs. (20) and (22), we obtain four possible phases: (i) n < n, r < r, (ii) n > n, r < r, (iii) n > n, r > r, and (iv) n < n, r > r. Given a fixed value of ϕ and ϕ, either (ii) or (iv) is possible. Realization of (ii) requires that and , thereforeIn 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 , we obtain the threshold PTD ratio for the transition from Phase 1 to Phase 2,In Phase 2, the average mRNA concentration becomeswhich is inversely proportional to the protein-to-DNA ratio.

From Eq. (23) and using , 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 (~103 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 ~ 30. A similar calculation for the upper bound of ribosome on a single mRNA leads to r ~ 10 since ribosome’s linear size is about 3 times larger than RNAP[62].

We take ϕ ≈ 0.2 according to the ref.[30], and estimate the gene allocation fraction of RNAP to be ϕ ~ 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 μm2/(0.2 μm2/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 k0 ≈ 0.05 s−1. Finally, we compute γ1 and γ2 using the above estimated parameters, and obtain γ1 ~ 1500, γ2 ~ 20,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 (ii) RNAPs bound to a promoter, n (iii) RNAPs nonspecifically bound to DNA, n (iv) free RNAPs, n. We assume a Michaelis-Menten relation for the number of promoter-bound RNAPs and nonspecifically bound RNAPs[40,65],Here c is the concentration of free RNAPs and K, K are the Michaelis constants. G is the total number of genes and g is the number of nonspecific binding sites per gene. Note that c = cF, with c the concentration of total RNAPs and F 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 is the transition rate from promoter-bound RNAPs to elongating RNAPs, τ0 is the time for transcribing a gene, and Λ = kτ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 + n + n + n, we can find a self-consistent equation for F,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 and nonspecifically bound RNAPs, F.

Since the left side of Eq. (31) monotonically decreases from cγ to 0 as F increases from 0 to 1 and the right side of Eq. (31) monotonically increases as F increases, we find that as the protein-to-DNA ratio increases, F increases. Previous studies have shown that F ≈ 0.1 for bacteria[40,41] and support the assumption that c is smaller or comparable to K[40,65] (note that nonspecific binding is characterized by K larger than K). In the following, we first assume a small F and , and show this leads to a behavior qualitatively equivalent to Phase 1 of the main text, albeit with a renormalization of the transcription constant k0. 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 can be written asHere g/G = ϕ is the gene allocation fraction of gene i, and F is the fraction of elongating RNAPs in the total RNAP pool. In the limit of a small F and , F becomes a constant:Therefore, the transcription rate corresponds to Phase 1 of our model, with .

In the limit that all RNAPs are actively transcribing and no nonspecifically bound RNAPs, and g = 0, we go back to the situation that all RNAPS are working with . 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 k0 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 = 1000, k0 = 1 min−1, G = 2000 (before gene duplication), K = 0.02〈c〉, K = 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〉 = ϕ = 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 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/G = n (where n is the upper bound of the number of RNAPs that can work on one gene simultaneously), and obtainWe 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 is small and c 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 (assuming that still holds). Because the transcription rate for a specific gene iswhen Ffree is comparable to 1, the transcription rate will be saturated as well and proportional to the gene copy number g. From Eq. (31), we find that  ≈ , therefore the threshold PTD ratio for F to be comparable to 1 is

Second, we can relax the assumption that (assuming still holds). When c ≈ K, the transcription rate can be saturated as well. Using  ≈ , 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), (36) and (37) shows that when n/ϕ < Λ/K and n < Λ, genes get saturated due to steric hindrance.

Reference[68] shows that for wild type E. coli in fast growth conditions the mRNA levels in the cell do not change when the DNA amount is lowered. Within our model this is consistent with Phase 1, and inconsistent with Phase 2, thus suggesting that and as discussed earlier.