Implications of differential size-scaling of cell-cycle regulators on cell size homeostasis

Accurate timing of division and size homeostasis is crucial for cells. A potential mechanism for cells to decide the timing of division is the differential scaling of regulatory protein copy numbers with cell size. However, it remains unclear whether such a mechanism can lead to robust growth and division, and how the scaling behaviors of regulatory proteins influence the cell size distribution. Here we study a mathematical model combining gene expression and cell growth, in which the cell-cycle activators scale superlinearly with cell size while the inhibitors scale sublinearly. The cell divides once the ratio of their concentrations reaches a threshold value. We find that the cell can robustly grow and divide within a finite range of the threshold value with the cell size proportional to the ploidy. In a stochastic version of the model, the cell size at division is uncorrelated with that at birth. Also, the more differential the cell-size scaling of the cell-cycle regulators is, the narrower the cell-size distribution is. Intriguingly, our model with multiple regulators rationalizes the observation that after the deletion of a single regulator, the coefficient of variation of cell size remains roughly the same though the average cell size changes significantly. Our work reveals that the differential scaling of cell-cycle regulators provides a robust mechanism of cell size control.

Introduction Cells must coordinate growth and division to keep their sizes within a finite range, which is vital for various biological functions and cellular fitness. For example, a very large cell volume may remodel the proteome and trigger senescence [1][2][3]. On the other hand, a very small cell volume may also reduce cellular fitness due to the limited protein copy numbers and the resulting large fluctuation in protein concentrations [4][5][6]. Therefore, cells must evolve ways to measure their sizes to accurately decide the timing of division.
Nearly half a century ago, the seminal work by Fantes et al. [7] already systematically discussed the possible mechanisms for cell size control, though the molecular basis of cell-cycle regulation was largely unclear at that time. Fantes et al. [7] described two threshold types for cell-cycle progression. In one case, cell division is triggered once the number of one activator reaches a threshold. Despite its constant concentration, this activator may form a mitotic structure or titrate the genome to promote cell division. It has been proposed that FtsZ and DnaA in bacteria act through this mechanism [8]. In the other case, the regulator concentration matters. This hypothesis is supported by the observations of different eukaryotes across yeast [9][10][11][12][13][14], mammals [2,15], and plants [16,17] that the concentrations of many cell-cycle regulatory proteins change as the cell size increases.
Recently, Chen et al. [12] proposed that the timing of cell-cycle entry can be achieved through the differential scaling of multiple cell-cycle regulators. If the activators exhibit superlinear scaling with cell size while the inhibitors exhibit sublinear scaling, then at some sufficiently large cell size, the activators dominate over the inhibitors and trigger cell-cycle entry. RNA sequencing also showed that the cell-cycle activators (inhibitors) are enriched in the subset of genes that exhibit increasing (decreasing) mRNA concentrations as the cell size increases [12]. Importantly, this mechanism may shed light on a long-standing mystery, especially for budding yeast: cells with some cell-cycle regulators knocked out can still maintain a narrow cell size distribution [18]. While the mechanism of cell size control based on the differential scaling of cell-cycle regulators is plausible, it is unclear whether it leads to robust growth and division. More quantitatively, how does the differential scaling of regulators determine the cell size distribution? The answers to these questions can deepen our understanding of cell-cycle regulation and cell size control.
Regarding the mechanisms of differential scaling, one of us has recently proposed a gene expression model at the whole-cell level in which the promoters of all genes compete for the resource of RNA polymerases (RNAPs) [19]. The model's key assumption is that RNAP is ratelimiting for transcription, supported by observations including the direct correlation among RNAP subunit level, RNAP occupancy, and transcription rate [20][21][22], haploinsufficiency of RNAP subunit genes [23,24], and increased growth rate after overexpression of RNAP subunits [25]. The model predicts that the effective binding affinities of promoters to RNAPs determine the nonlinear scaling of protein numbers with cell size. A gene with a relatively strong promoter exhibits a decreasing protein concentration as the cell size increases, while a gene with a relatively weak promoter exhibits the opposite behavior. The model suggests that the nonlinear scaling between protein numbers and cell size is an inevitable consequence of heterogeneous promoter strengths in the genome. The relation between the promoter strength and nonlinear scaling of protein level raises the possibility that cell-cycle inhibitors may have strong promoters and therefore exhibit sublinear scaling with cell size, while activators may have weak promoters and exhibit superlinear scaling; thus, the ratio of their concentrations depends on cell size and allows a cell to determine its size and decide when to divide [12,19].
In the main text, we focus on the case where most proteins are nondegradable while the cell-cycle regulators are degradable, supported by the observations that most proteins are stable in yeast and the minority of short-lived proteins are mainly enriched in the cell-cycle category [26]. Our conclusions are qualitatively similar for nondegradable cell-cycle inhibitors (S1 Appendix and S7 Fig). We mainly consider a simple scenario with constant gene copy numbers, rationalized by imagining a cell that replicates its genes right before cell division. Therefore, the gene copy number is constant throughout the cell cycle. We also discuss a modified model in which the cell replicates its genes when the ratio of the activator to the inhibitor rises to a threshold and divides after a constant time (S1 Appendix and S6 Fig). Our conclusions are qualitatively similar between the simple scenario and the modified model. We mainly consider symmetric division for simplicity, but our model equally applies to asymmetric division (S1 Appendix and S11 Fig).
In the following, we first discuss the case of a single activator and a single inhibitor. When the ratio of the activator's concentration to the inhibitor's concentration rises to a threshold value θ, the cell divides ( Fig 1A). Based on the deterministic and stochastic versions of this simplified model, we derive the conditions of robust cell cycle and the expression for the cell size and its coefficient of variation (CV). We rationalize the proportional relation between cell size and ploidy, and the sizer behavior of cell size control. Later, we generalize our model to the case of multiple activators and inhibitors, and explain the experimental observations of mutants with regulators deleted.

Stable cell cycle exists between two critical division thresholds
We use the gene expression model at the whole-cell level which explicitly incorporates the competition between genes for the resource of RNAPs and between mRNAs for the resource of ribosomes [19,27]. In this model, the concentration of free RNAPs sets the transcription initiation rate through the Michaelis-Menten (MM) mechanism. The MM constant quantifies each gene's ability to recruit RNAPs. We first consider a model with one activator and one inhibitor and set the MM constants of the activator and inhibitor to be K n,act and K n,inh (K n,act > K n,inh ), while the MM constants of other genes are equal to K n for simplicity. The cell divides once the ratio of the activator's concentration to the inhibitor's concentration rises to the threshold value [12,19], Details of the gene expression model are included in Methods and also in [19]. One of our ultimate goals is to understand how cell size depends on the scaling behaviors of the cell-cycle regulatory proteins. We successfully derive the expression of the cell size at cell birth V b (see details in Methods). The essence is to find the RNAP number at birth n b . Because the total RNAP concentration c n is essentially constant, one immediately obtains whereK � yK n;act À K n;inh . Here a is the ratio between the cell volume and the nuclear volume, which is approximately constant, supported by experiments [28,29]. n c is the maximum number of RNAPs that the entire genome can hold (Eq 19). Intriguingly, we find that V b only depends onK , a linear combination of K n,act and K n,inh (Fig 1C and 1D). The comparison between the predictions and simulations shows good agreement (Fig 1D and 1E, see Methods for the simulation details). We note that the fraction of free RNAPs at cell division F n,d (see Eq 17 in Methods) must satisfy 0 < F n,d < 1, which is equivalent to requiring the numerator and denominator of Eq 2 to be positive. Therefore,K must satisfy from which we obtain the conditions of a stable cell cycle, K n;inh K n;act < y < c n þ K n;inh c n þ K n;act : ð4Þ The cell can grow and divide within a finite range of θ. The two critical threshold values y 1 ¼ K n;inh K n;act and y 2 ¼ c n þK n;inh c n þK n;act respectively set the lower bound and the upper bound. The predictions of a finite range of θ for a stable cell cycle are consistent with simulations (Fig 1B-1E).
We provide an intuitive explanation of the critical threshold values in the following. We note that the inhibitor's concentration keeps decreasing as the cell size increases since its advantage over other proteins becomes weaker and weaker with the increasing concentration of free RNAPs. The activator's concentration first increases as we predict. However, when the cell size is so big that genes get saturated by RNAPs, and mRNAs get saturated by ribosomes, the protein production rates will also reach their maximum values, which is called Phase 3 of gene expression in [27]. Mathematically, this corresponds to the situation in which the free RNAP concentration is above most genes' transcription MM constants (Eq 6) and the free ribosome concentration is above most mRNAs' translation MM constants (Eq 8). Thus, the protein numbers of the degradable activator and inhibitor both become constant, and their concentrations decrease inversely with the cell size in the large cell size limit. Therefore, for the two curves of c act and θc inh to cross, the threshold θ must be within a finite range (Fig 1F).
In the modified model incorporating gene replication (S1 Appendix and S6 Fig), we only need to slightly modify the expression of V b , and the two critical threshold values remain the same. For the case of nondegradable inhibitor, we also derive V b and provide an illustration for the critical threshold values (S1 Appendix and S7(A) and S7(B) Fig). It is also possible that one mechanism, such as activator-accumulation or inhibitor-dilution, dominates or that the cell integrates the information from each regulator through logic gates [30][31][32]. We also explore these cases and perform similar analysis (S1 Appendix and S8 Fig).
We define A as the amplitude of the oscillation of c act /c inh , i.e., the difference between θ and the minimum of c act /c inh in one cell cycle (Fig 1B, middle). We find that the amplitude A approaches 0 as θ approaches the two critical threshold values (Fig 1G and 1H). The critical point θ 1 corresponds to F n,d ! 0 (see details in Methods). In this limit, the fraction of free RNAPs at cell birth F n,b also approaches zero since the RNAP number decreases after cell division. The critical point θ 2 corresponds to F n,d ! 1 (Methods). In this limit, almost all RNAPs are free, which means that the genes are fully saturated by RNAPs at cell division, and therefore F n,b is still close to 1. In both cases, F n,b is close to F n,d . On the other hand, we note that the maximum and minimum of c act /c inh are respectively set by the maximum and minimum of F n (Eq 22), i.e., F n,d and F n,b . Therefore, the amplitude A approaches 0 as θ approaches the two critical threshold values. We derive an approximate expression of A (Methods), which agrees well with simulations ( Fig 1H). Later we show that the magnitude of A is crucial for a stable cell cycle when the threshold θ is subject to noise.

The cell size at birth is proportional to the ploidy
It is well known that the cell size is usually proportional to the ploidy [33][34][35]. Previous hypotheses about this observation were reviewed in [36,37]. Our model provides a new perspective on the linear relationship between cell size and ploidy. We remark that our model is invariant as long as the size-to-ploidy ratio is constant. This can be seen from Eq 7 where the doubling of cell size and gene copy numbers leaves the fraction of free RNAPs in total RNAPs F n invariant. Note that the above argument requires that the RNAP concentration c n is independent of ploidy, which is true in our model (Eq 14) and also supported by experiments [2]. Therefore, the ratio c act /c inh in a diploid cell is the same as that in a haploid cell as long as the size of the diploid cell is twice the size of the haploid cell. Thus, the cell volume at birth V b doubles in the diploid cell relative to the haploid cell (Fig 1E, inset). This can also be seen from Eq 2, where n c , the total number of RNAPs that can simultaneously bind to the genome, is proportional to ploidy (see Eq 19). In summary, the free RNAP fraction F n measures the size-toploidy ratio and passes on this information to c act and c inh , which then determine the cell size at division. Our model is based on the limiting role of RNAP, supported by experiments [20][21][22][23][24][25]. Our model predicts that F n increases with the size-to-ploidy ratio (see Eq 7), consistent with the observation that the number of RNAPs on the genome subscales with cell size [22].

The cell cycle model with stochastic division threshold reveals a sizer behavior
We extend the model by adding noise to the division threshold θ to mimic randomness in real biological systems [4][5][6]. We assume that the threshold of each generation is independent and identically distributed. Since the threshold value is unlikely to be zero or infinite in real biological systems, we assume θ to be within a finite interval ð � y À Dy; � y þ DyÞ where � y is the mean value. For a given parameter set, we simulate 5000 generations and trace one of the two daughter cells after division (Fig 2A-2C). Intriguingly, we find that the cell size at division V d is uncorrelated with the cell size at birth V b (Fig 2D), in concert with the sizer behavior found in experiments [38][39][40]. Therefore, ΔV � V d − V b and the doubling time T D are negatively correlated with V b (Fig 2E and 2F). The sizer behavior is robust against the parameter � y (Fig 3A). This can be seen from Eq 2, where the cell size at cell birth is independent of that of the previous cell cycle. Intuitively, the sizer behavior is due to the degradation of the regulatory proteins: their numbers mainly depend on the current cell size and are insensitive to the history. Sizer due to the fast degradation of cell-cycle regulators is also proposed in [41]. We also study the correlation between V d and V b for the case of nondegradable inhibitor and also find an approximately zero correlation between V d and V b (S1 Appendix and S7(C)-S7(E) Fig).
Eukaryotic cells such as the budding yeast Saccharomyces cerevisiae [38,42], the fission yeast Schizosaccharomyces pombe [43][44][45], and mammalian cells [40,46] are usually imperfect sizers. For budding yeast, this phenomenon could result from the unequal partitioning of Whi5 during asymmetric cell division [47], but it cannot explain the imperfect sizer observed in symmetrically dividing cells. To explain the imperfect sizer, we modify our model by assuming that the cell-cycle entry is probabilistic, with its probability increasing with the activatorto-inhibitor ratio, similar to some previous works [14,42,[48][49][50][51]. This modified model successfully reproduces the imperfect sizer (S1 Appendix and S9 Fig). In addition, we remark that different modes of cell size control in different cell-cycle phases also contribute to the imperfect sizer over the entire cell cycle [42,43,47].

Differential scaling behaviors of the activator and inhibitor narrow cell size distribution
We also derive the analytical expressions of the mean cell size at birth and the CV of cell size at birth (see details in Methods). The predictions agree well with the simulations under various � y, K n,act and K n,inh ( Since K n,act and K n,inh determine how nonlinear the activator and inhibitor are relative to the cell size, this observation suggests that the more differential the scaling behaviors of the cell-cycle regulators are, the smaller the CV of the cell size distribution is. In the presence of noise, the range of threshold that allows robust cell growth and division becomes smaller. We find that the mean threshold � y must satisfy (Methods) When Δθ = 0, Eq 5 becomes A(θ) > 0, which is precisely the condition of robust division for the deterministic model. We note that hV b i can not approach zero or infinity in the presence of noise (Fig 3B), in contrast to the deterministic model. This result suggests that too small or too large cell size is unstable against noise, and cell size must stay in an appropriate range. where � y ¼ 0:7, Δθ = 0.1. If y = 2 ð � y À Dy; � y þ DyÞ, we reset it as � y. https://doi.org/10.1371/journal.pcbi.1011336.g002

Multiple regulators render robust cell size distribution against gene deletion
Experimentally, it has been found that deleting one regulator's gene of budding yeast can change the mean cell size significantly, but the CV of cell size is approximately the same before and after deletion [12,52] (S3 Fig). This suggests that there should be multiple cell-cycle regulators working in concert. To explore the influence of gene deletion, we extend the model to multiple activators and inhibitors. Each regulator has one gene copy, and the numbers of activators and inhibitors are respectively g act and g inh . We assume that all the activators' promoters have the same effective binding affinity to RNAPs, i.e., K n,act , and all the inhibitors' promoters have the same effective binding affinity to RNAPs, i.e., K n,inh . The cell divides once the ratio between the concentration of all activators and the concentration of all inhibitors reaches the threshold θ, i.e., g act c act g inh c inh ¼ y, where c act (c inh ) is the concentration of a particular activator (inhibitor). We define the equivalent threshold asỹ � g inh g act y such that the division condition becomes c act c inh ¼ỹ. Therefore, all the results above for the model of one activator and one inhibitor hold as long as we substitute θ withỹ in the deterministic model and substitute � y and Δθ with g inh g act � y and g inh g act Dy in the stochastic model. After deleting one activator, the equivalent threshold changes to g inh g act À 1 y, larger than the wild type (WT) value g inh g act y. Because V b increases with θ (Fig 1E), In the following, we provide a more intuitive explanation for the observed constant CV after regulator deletion. We note a wide range of equivalent division thresholdỹ in which V b scales linearly withỹ (Fig 4A). In this range, because V b /ỹ, the CV of V b is equal to the CV ofỹ in the stochastic model, which is constant. Since deleting the cell-cycle regulator is equivalent to shifting the mean ofỹ while not altering the CV ofỹ, the CV of V b is robust against gene deletion (Fig 4B). Within our model, the simultaneous deletion of multiple activators and inhibitors may lead to a small shift of the equivalent threshold, so the CV in size can change mildly (S12(A) and S12(B) Fig). Meanwhile, our model also predicts that if the shift is drastic and goes beyond the linear regime (e.g., deleting multiple activators or multiple inhibitors), the mutant can have a larger CV in size or even become inviable (S12 Fig), in agreement with experiments [12,45,54]. Thus, our model provides a unified framework to understand these observations. We acknowledge that our model cannot explain the increase in CV for the fission yeast mutant wee1-50 cdc25Δ, which results from the quantized cell-cycle times [44,55,56].
We remark that given multiple regulators, the lack of one copy of an inhibitor's gene in a diploid cell does not necessarily halve the cell size. The redundancy of multiple inhibitors naturally explains the observation that the cell size of the diploid WHI5 heterozygous mutant is still bigger than the haploid WT cell [9].

Discussion
In this work, we study a model combining cell size control and the differential scaling of gene expression with cell size. As the cell grows, the ratio between the activator and inhibitor increases and eventually triggers cell division when the ratio is above a threshold value [12,19]. The activator has a relatively weak promoter and therefore scales superlinearly with cell size. Meanwhile, the inhibitor has a relatively strong promoter and therefore scales sublinearly. Our simplified model captures the essential features of cell size control, and the cell can periodically grow and divide, maintaining size homeostasis. We derive the analytical expression of the cell size and determine the conditions of stable cell cycle, revealing that global gene expression can impose a constraint on cell size control. We also provide a new perspective on the linear relationship between cell size and ploidy [33][34][35], showing that the free RNAP fraction F n measures the size-to-ploidy ratio and passes on this information to c act and c inh . All of our theoretical predictions are verified by simulations.
We also propose a stochastic version of our model in which the division threshold is random in each cell cycle. We find that the cell size at division is uncorrelated with the cell size at birth, known as a sizer [38][39][40]. The more differential the cell-size scaling of the cell-cycle regulators is, the narrower the cell-size distribution is, which unravels the role of differential scaling in size homeostasis.
We also extend our model to multiple activators and inhibitors. The cell size becomes larger after one activator is deleted, while it becomes smaller after one inhibitor is deleted, which agrees with experiments [12,52,53]. Interestingly, while the mean cell size can change significantly after regulator deletion, the CV of the cell size distribution can be approximately constant, consistent with experiments [12,52]. To our knowledge, no previous theoretical works explained this counter-intuitive phenomenon. Our model further predicts that though the CV of cell size is robust against deletion of one single regulator, it may increase after deletions of multiple activators or multiple inhibitors. It will be interesting to test our prediction by systematically measuring the CV of budding yeast strains with multiple activators or multiple inhibitors knocked out, which is still largely unexplored.
In the experiment that exchanged the promoters of CLN2 and WHI5, the mean and CV of the budding yeast cell size increased [12]. Our model can also reproduce the increase in the mean and CV of cell size after swapping the promoters of the activator and inhibitor (S10(E) Fig). Meanwhile, we note that the outcome of exchanging CLN2 and WHI5 promoters can be more complex due to factors such as regulations acting on the promoters. Whi5 inhibits the transcription of CLN2 by binding the SBF complex, and Cln2 can phosphorylate Whi5 and dissociate it from SBF, forming a positive feedback loop [57][58][59]. Therefore, exchanging the WHI5 promoter and CLN2 promoter may rewire the network and alter the regulatory relationship. Besides, systematic experiments of exchanging promoters are still lacking, and it is unclear whether the increase in cell size and its CV after promoter exchange is a universal phenomenon. Future experimental work on this issue will be very helpful. Additionally, experiments properly truncating the promoters of the activator and inhibitor might change the binding affinities of promoters without altering the regulatory network, which can help test the prediction that the nonlinear scaling behaviors of regulators influence the breadth of the cell size distribution.
Wang et al. have analyzed the RNA-seq data of [12] and showed that genes exhibiting sublinear scaling tend to have higher mRNA production rates. In particular, negative cell-cycle regulators are enriched in the sublinear regime [19]. It will be helpful to test further our assumption that cell-cycle activators tend to have lower RNAP occupancy due to their weak promoters, and vice versa for inhibitors, e.g., by combining RNA-seq and ChIP-seq data of budding yeast cells arrested in the G1 phase. Meanwhile, we acknowledge that this assumption may not apply to some cell-cycle regulators.
Conceptually, the ideas that limiting RNAP accounts for the increase of transcription rate with cell size [20,21,60,61] and that different recruitment abilities of promoters to RNAP generate nonlinear scaling of gene expression with cell size [19,47,62] have been brought forward before. Our work shares some similarities with previous works. Heldt et al. [47] studied a mathematical model of cell size control for budding yeast in which the sublinear scaling of Whi5 is introduced by a strong promoter while Cln3 exhibits linear scaling. We remark that their model was designed to understand the cell size control mechanism of budding yeast and, in particular, to explain the size of cells with abnormal WHI5 copy numbers and rule out alternative mechanisms. However, our model is not constructed for a specific organism and is more poised to obtain some general conclusions that may apply to multiple organisms. Moreover, our gene expression model is built at the whole-cell level and reveals the connection between global gene expression and cell size control. Furthermore, our model with multiple regulators explains the counter-intuitive phenomenon that the CV of cell size can change mildly after the deletion of key regulators [12,18]. While Heldt et al. [47] assumed a linear size-scaling for activators, we include the possibility of superlinear scaling for activators, given that it could also facilitate cell size homeostasis [12][13][14]. Other works also employed the idea of changing regulator concentrations [63,64].
Our results support the hypothesis that cells can measure their sizes through the differential scaling of cell-cycle regulatory proteins. Our model is not designed for specific cell-cycle regulators or organisms, so that some specific experimental phenomena could be beyond our model. Nevertheless, our model can be extended by including cell-cycle regulatory networks [54,[65][66][67].

Details of the gene expression model
We list the variables used in the main text in Table 1. A free RNAP can bind to a promoter and become an initiating RNAP, which can start transcribing or hop off the promoter without transcribing. Thus, the concentration of free RNAPs sets the transcription initiation rate through the Michaelis-Menten (MM) mechanism. The MM constant quantifies each gene's ability to recruit RNAPs. The time dependence of the mRNA copy number of one particular gene i follows [19], where c n,free is the free RNAP concentration in the nucleus and K n,i is the MM constant. A larger K n,i represents a weaker promoter, while a smaller K n,i represents a stronger promoter. g i is the copy number of gene i and τ m,i is the corresponding mRNA lifetime. Γ n,i is the rate for a promoter-bound RNAP to become a transcribing RNAP. Interaction between genes due to the limiting resource of RNAPs is through the free RNAP concentration. For one copy of each gene, the number of actively transcribing RNAPs is proportional to the probability that the promoter is bound by an RNAP, based on the assumption that the number of RNAPs starting transcription per unit time is equal to the number of RNAPs finishing transcription per unit time. This assumption is justified by the fact that it usually takes shorter than one minute for an RNAP to finish transcribing a gene [68] so that the system can quickly reach a steady state. We introduce the relative fraction of free RNAPs in total RNAPs F n such that c n,free = c n F n where c n is the concentration of total RNAPs in the nucleus. Using the conservation of the total RNAP number, we can find the free RNAP concentration self-consistently [19], The right side of Eq 7 represents RNAPs binding to the promoters or transcribing. Here, L i is the length of gene i in the unit of codon, and v n is the RNAP elongation speed. When the number of transcribing RNAPs on one copy of gene i, n i , does not change with time, we must have Γ n,i P b,i = v n n i /L i . Here the left side is the number of RNAPs starting transcription per unit time, where P b,i � F n c n /(F n c n + K n,i ) is the probability that the promoter is bound by an RNAP. The right side is the number of RNAPs finishing transcribing per unit time. Therefore, we obtain n i = P b,i Γ n,i L i /v n , which is included in the right side of Eq 7. A similar model as Eq 6 applies to the protein number, in which the free ribosome concentration sets the initiation rate of translation. In particular, we denote the protein number of RNAPs and ribosomes as n and r, and set their index as i = 1 and i = 2 respectively. Γ r,i is the translation initiation rate for a ribosome binding on the ribosome-binding site of mRNAs, c r is the concentration of total ribosomes in the cytoplasm, K r,i quantifies the binding strength of mRNAs to ribosomes. τ p,i is the protein lifetime, and F r is the fraction of free ribosomes in the total pool of ribosomes. F r satisfies the following equation due to the conservation of total ribosome numbers, The right side of Eq 9 represents ribosomes binding to the ribosome-binding sites of mRNAs or translating. v r is the ribosome elongation speed.
In this work, we assume that the cell size is simply proportional to the total protein mass, supported by experiments [69][70][71][72]. For simplicity, we assume that all genes except the activators and inhibitors share the same transcription MM constant K n . Heterogeneous K n,i among genes has no significant effects on our results (S5 Fig). All genes except those of RNAPs and ribosomes share the same transcription initiation rate Γ n (see Details of simulations in Methods). All the mRNAs share the same effective binding strength to ribosomes K r , the same translational initiation rate Γ r , and the same mRNA lifetime τ m . The cell-cycle regulators share the same protein lifetime τ p and other proteins are nondegradable.
In the simplified model with a single activator and a single inhibitor, the cell divides once the activator-to-inhibitor ratio c act /c inh rises to the threshold value θ. The cell can reach a steady state, growing and dividing periodically (Fig 1B).

The non-monotonic behavior of c act /c inh in one cell cycle
There is a transient increase in the inhibitor concentration and a transient decrease in the activator concentration right after cell division (Fig 1B). This is because of the reduced RNAP number at cell birth compared to cell division. The fewer total RNAP number leads to fewer free RNAPs given the same gene copy number, which means a decreasing concentration of free RNAPs (S1 Fig). The initial low concentration of free RNAPs lets the inhibitor increase faster than other proteins because of its strong promoter, leading to an increasing inhibitor concentration. However, the advantage of the inhibitor over other proteins becomes weaker as the free RNAP concentration increases since all promoters are fully saturated in the limit of very high free RNAP concentration. Therefore, the inhibitor's concentration decreases as the cell grows. The opposite situation applies to the activator. After the transient regime, c act /c inh increases, eventually reaches the threshold θ and triggers cell division.
A more quantitative explaination is as follows. Using Eq 6 about the activator and the inhibitor, we obtain The left-hand side of Eq 7 decreases monotonically with F n , while the right-hand side increases monotonically. Therefore, F n jumps to a smaller value after cell division since n is reduced by half after division. According to Eq 10, the drop in F n turns d dt m act m inh � � into a negative value, which means that m act /m inh starts to decrease. Consequently, the ratio of the production rate of the activator to that of the inhibitor declines, leading to the reduction of c act /c inh . We note that the transient decrease of c act /c inh at the beginning of cell cycle is abrupt in the limit of short lifetimes of mRNAs and regulatory proteins (S1 Appendix). Then, d dt m act m inh � � becomes positive again as F n increases, causing the increase in c act /c inh .

Derivation of the cell size at birth and the two critical threshold values
Assuming a short lifetime of m i [68] in Eq 6, we can approximate the mRNA copy number as For those proteins that are not cell-cycle regulators, substituting Eq 11 into Eq 8, we get d dt At steady state, d dt p i n À � ¼ 0. Therefore, The ratio of copy numbers between two proteins is determined by their gene copy numbers and transcription initiation rates. The total protein mass of the cell is ∑ i p i L i , in the number of amino acids. The cell size V(t) = ∑ i p i (t)L i /ρ, where ρ is the ratio between the total protein mass and cell size, which is a constant in our model, consistent with experiments [69][70][71][72]. Therefore, the concentration of total RNAPs in the nucleus is approximately constant and can be approximated as c n ¼ anðtÞ VðtÞ ¼ arG n;n g n P i G n;i g i L i ; ð14Þ Here we have neglected the contribution of the cell-cycle regulators to the total protein mass. Eq 14 shows that the RNAP concentration is independent of ploidy since the doubling of all genes' copy numbers leaves c n invariant. Assuming a short lifetime τ p , using Eq 11 and the quasi-steady state approximation about p act and p inh in Eq 8, we get p act ðtÞ ¼ G n g act t m G r t p F n ðtÞc n F n ðtÞc n þ K n;act In the simple case of g act = g inh = 1, substituting Eqs 15 and 16 into the division condition Eq 1 (note that c act /c inh = p act /p inh ), we obtain the fraction of free RNAPs at cell division F n;d ¼ yK n;act À K n;inh ð1 À yÞc n : We also simplify Eq 7 using the assumption that K n,i = K n for most genes, nð1 À F n Þ � n c F n c n F n c n þ K n ; ð18Þ where n c is the maximum number of RNAPs the entire genome can hold [19], Substituting Eq 17 into Eqs 11 and 18, we get the mRNA copy number and the RNAP number at cell division n d ¼ n c c n ð1 À yÞK ðð1 À yÞc n ÀK Þðð1 À yÞK n þK Þ : ð21Þ Using n b = n d /2 where n b is the RNAP number at cell birth and the RNAP concentration Eq 14, we obtain Eq 2 (see S1 Appendix for the discussion of asymmetric division). The condition 0 < F n,d < 1 in Eq 17 leads to the condition of stable cell cycle Eq 3 so that θ 1 corresponds to F n,d = 0 while θ 2 corresponds to F n,d = 1. We mathematically prove that the cell cycle is stable against infinitesimal perturbation (S1 Appendix).

Derivation of the amplitude A
Using Eqs 15 and 16, the ratio of the activator's concentration to the inhibitor's concentration at a given time t during the cell cycle can be approximated as c act ðtÞ c inh ðtÞ � F n ðtÞc n þ K n;inh F n ðtÞc n þ K n;act : ð22Þ The ratio reaches its minimum when F n is equal to its minimum value F n,b at cell birth. Therefore, A � y À F n;b c n þ K n;inh F n;b c n þ K n;act ; ð23Þ where F n,b satisfies n b ð1 À F n;b Þ ¼ n c F n;b c n F n ;b c n þ K n : ð24Þ As θ increases, n b increases (Eq 21). According to Eq 24, more RNAPs at cell birth generates more free RNAPs, leading to a larger F n,b . Therefore, F n;b c n þK n;inh F n;b c n þK n;act also increases with θ, consistent with the fact that A first increases and then decreases with θ.

Derivation of hV b i and the CV of V b
When Δθ is small, a function f of θ can be approximated as Therefore, we can approximate the average and variance of f(θ) as Ef ðyÞ � f ð � yÞ; Var f ðyÞ � ðf 0 ð � yÞÞ 2 s 2 : Here σ is the standard deviation of θ. Using Eqs 2 and 26, we obtain CV ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi where � K � � yK n;act À K n;inh .

The condition of robust cell division for the stochastic model
We denote θ k as the division threshold value of the k-th cell cycle. Given � y and Δθ, to ensure robust division, we require that the minimum value of the activator-to-inhibitor ratio in the kth cell cycle θ k−1 − A(θ k−1 ) is smaller than θ k , where A(θ k−1 ) is the amplitude of the k-th cell cycle, which is a function of θ k−1 . This condition must hold for any θ k−1 and θ k . Therefore, maxfy À AðyÞg < minfyg; y 2 ð � y À Dy; � y þ DyÞ: Because θ − A(θ) increases with F n,b (see Eq 23) and F n,b increases with θ, θ − A(θ) increases with θ. Therefore, Substituting Eq 30 into Eq 29, we obtain Eq 5.

Details of simulations
All simulations were done in MATLAB (version R2021b). We summarize some of the parameters used in simulations in Table A in S1 Appendix. To find reasonable values of Γ n,i , we first set n c = 10 4 as in [19]. Assuming the number of RNAPs is about 10% of ribosomes i.e. n = 0.1r, from Eq 13 we get G n;n ¼ 0:1G n;r g r g n : ð31Þ Assuming the length of genes except those of RNAPs and ribosomes is L, from Eq 19 we get, G n ¼ ðn c À P i g i Þv n À G n;n g n L n À G n;r g r L r L P i>2 g i : ð32Þ Substituting Eq 11 into Eq 9, we get Substituting Eqs 11 and 33 into Eq 8 for r, we find the growth rate as In this work, we set the attempted growth rate μ = (ln 2)/2 h -1 . Substituting Eqs 31, 32 to 34, and letting F r = 0, which is merely a numerical convenience, we obtain from which we find Γ n,n and Γ n . Our theoretical results do not rely on the initial values of simulations, but appropriate initial values make it easier to reach the steady state. We first set an attempted cell volume V 0 , and then we get n(0) = c n V 0 /a where c n is calculated using Eq 14. Finally, we determine p i (0) (i > 1) using Eq 13 and m i (0) using Eq 11. After symmetric division, mRNAs and proteins are equally distributed to the two daughter cells. We also have a requirement on the minimum cell-cycle duration to avoid cell division immediately after cell birth. Before we take the simulation results, we run the simulations for several generations to ensure that the cell reaches a steady state. the activator's threshold θ act given a fixed threshold θ inh = 350 μm -3 for the inhibitor. In this figure, K n,act = 12000 μm -3 , K n,inh = 4000 μm -3 . (TIF) S9 Fig. Imperfect sizer. (A) We assume that the probability for a cell with the activator-toinhibitor ratio θ in the infinitesimal interval ½y; y þ dyÞ to divide is kðyÞdy. The cell born with a smaller size is more likely to divide at a smaller size because it has a finite probability to divide before it reaches the birth size of the larger cell. (B-G) Simulations of V b and V d . We simulate 5000 generations and trace one of the two daughter cells after division. For each generation, we first let the cell grow for a minimal time (20 min), which has minor effect on the slope of V d vs. V b . We then use inverse transform sampling to generate the random variable θ, which is the activator-to-inhibitor ratio at cell division and has a probability density function p(θ) (see S1 Appendix for more details). In inverse transform sampling, we first take a random number x from a uniform distribution between 0 and 1. θ is then solved from F(θ) = x, where

Supporting information
Þdy 0 is the cumulative distribution function of θ. (B-D) C 1 = 1, C 2 = 2. The slope is similar to the imperfect sizer in [40,[43][44][45]. (E-G) C 1 = 200, C 2 = 1. The slope is similar to the near-adder in [46]. In (C-D,F-G), each bin includes an equal number of data points. Error bars: mean±SEM. In this figure, K n,act = 12000 μm -3 , K n,inh = 4000 μm -3 . (TIF) In (B-G), the inhibitor is distributed between the mother and daughter cells at division as Whi5 in budding yeast, with η = 0.45 [9]. α = 10. In this figure, K n,act = 12000 μm -3 , K n,inh = 1000 μm -3 .  4). (A) In the model of multiple regulators, the deletion of activators corresponds to an increase in the equivalent division threshold. If the threshold change is small, e.g., act1Δ, the cell size at birth V b is still proportional to the equivalent threshold (the yellow region) so that the CV of V b remains constant. If the threshold change is significant, e.g., act1Δ act2Δ, the equivalent threshold is beyond the critical division threshold (the white region), and the cell becomes inviable (the red hollow circle). However, additional deletion of inh1 rescues the cell. (B) The distributions of V b for WT, act1Δ and act1Δ act2Δ inh1Δ. (C) The deletion of inhibitors is equivalent to a decrease in the division threshold. If the threshold change is small, e.g., inh3Δ, V b is still proportional to the equivalent threshold so that the CV of V b remains constant. If the threshold change is significant, e.g., inh2Δ inh3Δ, V b becomes a nonlinear function of the equivalent threshold (the purple region), and the CV increases. (D) The distributions of V b for WT, inh3Δ, and inh2Δ inh3Δ. In (B, D), the mutants' relative changes in average cell size and CV compared with WT are shown at the top. In simulations, we introduce different weights χ act,i and χ inh,i for different regulators. The cell divides once ∑ i χ act,i c act,i /∑ χ inh,i c inh,i = θ so that the change of average cell size can be different after deleting different activators or inhibitors. The K n,act of 10 activators is an arithmetic sequence from 9000 to 15000 μm -3 . The K n,inh of 10 inhibitors is an arithmetic sequence from 750 to 1250 μm -3 . Meanwhile, we also introduce heterogeneous K n,i for other genes except the cell-cycle regulators, which follows a lognormal distribution with the mean K n = 6000 μm -3 and the CV equal to 0.5. Our conclusion regarding the CV of cell size after gene deletions remains valid in the presence of heterogeneous K n,i . χ act,1 = χ act,2 = 2, χ act,3 = . . . = χ act,6 = 1, χ act,7 = . . . = χ act,10 = 0.5. χ inh,1 = χ inh,2 = 2, χ inh,3 = . . . = χ inh,6 = 1, χ inh,7 = . . . = χ inh,10 = 0.5. � y ¼ 0:45, Δθ = 0.1. (TIF)