Characterizing non-exponential growth and bimodal cell size distributions in fission yeast: An analytical approach

Unlike many single-celled organisms, the growth of fission yeast cells within a cell cycle is not exponential. It is rather characterized by three distinct phases (elongation, septation, and reshaping), each with a different growth rate. Experiments also showed that the distribution of cell size in a lineage can be bimodal, unlike the unimodal distributions measured for the bacterium Escherichia coli. Here we construct a detailed stochastic model of cell size dynamics in fission yeast. The theory leads to analytic expressions for the cell size and the birth size distributions, and explains the origin of bimodality seen in experiments. In particular, our theory shows that the left peak in the bimodal distribution is associated with cells in the elongation phase, while the right peak is due to cells in the septation and reshaping phases. We show that the size control strategy, the variability in the added size during a cell cycle, and the fraction of time spent in each of the three cell growth phases have a strong bearing on the shape of the cell size distribution. Furthermore, we infer all the parameters of our model by matching the theoretical cell size and birth size distributions to those from experimental single-cell time-course data for seven different growth conditions. Our method provides a much more accurate means of determining the size control strategy (timer, adder or sizer) than the standard method based on the slope of the best linear fit between the birth and division sizes. We also show that the variability in added size and the strength of size control in fission yeast depend weakly on the temperature but strongly on the culture medium. More importantly, we find that stronger size homeostasis and larger added size variability are required for fission yeast to adapt to unfavorable environmental conditions.

Let V b and V d denote cell sizes at birth and at division in a particular generation, respectively, and let V s denote cell size in the septation phase, which is assumed to be a constant. In the main text, we have stated that (i) the generalized added size V α s − V α b in the elongation phase is Erlang distributed with shape parameter N 0 and mean M 0 = N 0 g 0 α/a and (ii) the generalized added size V α d − V α s in the reshaping phase is also Erlang distributed with shape parameter N 1 and mean M 1 = N 1 g 1 α/a. To see this, note that when V b is fixed, the cell size in the elongation phase is given by V (t) = V b e g0t . Since the transition rate from one stage to the next at time t is equal to aV (t) α , the distribution of the transition time T is given by This shows that Hence V α b (e g0αT − 1) is exponentially distributed with mean g 0 α/a. Note that V α b (e g0αT − 1) is the generalized added size in a particular cell cycle stage. As a result, the generalized added size in the elongation phase, V α s − V α b , is the independent sum of N 0 exponentially distributed random variables with mean g 0 α/a. This shows that V α s − V α b has an Erlang distribution with shape parameter N 0 and mean M 0 = N 0 g 0 α/a. The fact that V α d − V α s has an Erlang distribution with shape parameter N 1 and mean M 1 = N 1 g 1 α/a can be proved in the same way.

A.2 Connection between our model and the conventional model
In our model, we use the parameter α to characterize the strength of cell size control. However, in the conventional model given in Eq. (1) in the main text, the size control strength is characterized by the parameter β. This raises the following natural question: what is relationship between the two parameters α and β?
To answer this, we first focus on the cell size dynamics in the elongation phase. Recall that the generalized added size ∆ 0 = V α s − V α b in the elongation phase is independent of V b and has an Erlang distribution. When the fluctuations of V b and V s are small, we have the following Taylor expansion: where v b and v s are the typical values of V b and V s , respectively. Therefore, we obtain This shows that where 0 is a noise term independent of V b . Similarly, for the cell size dynamics in the reshaping phase, we have where v d is the typical value of V d and 1 is a noise term independent of V b . Combining the above two equations, we obtain Under symmetric division, we have v d = 2v b and thus where the right-hand side is a noise term independent of V b . Recall that in the conventional model, we have where 0 ≤ β ≤ 2 and γ ≥ 0 are two constants and is Gaussian white noise. This shows that V d − βV b is a noise term independent of V b . Comparing the above two equations, we conclude that which gives the relationship between α and β.
B Cell size distribution B.1 Model with deterministic partitioning: general case Here we compute the cell size distribution of lineage measurements for the model with deterministic partitioning. The microstate of the cell can be represented by an ordered pair (k, x), where k is the cell cycle stage and x is the cell size. Let p k (x) denote the probability density function of cell size when the cell is in stage k. Then the evolution of cell size dynamics is governed by the master equation For the first, second, and fourth equations, the first term on the right-hand side describe cell growth and the remaining two terms describe transitions between cell cycle stages. For the third equation, the two terms on the right-hand side describe cell cycle stage transitions. For the first equation, the middle term on the right-hand side describes partitioning of cell size at division. For convenience, we next focus on the αth power of cell size, y = x α . Letp k (y) denote the probability density function of the αth power of cell size when the cell is in stage k. It is easy to see that p k (x) andp k (y) are related byp In other words, we have αyp k (y) = xp k (x).
Based on this relation, the master equation can be rewritten as (2) To proceed, for each cell cycle stage k, we introduce the Laplace transform wherep(y) = N k=1p k (y). Then Eq.
(2) can be converted to the following differential equations: At the steady state, the above equations reduce to where A 0 = αg 0 /a and A 1 = αg 1 /a. Note the the above equations are recursive with respect to k, which means that F k−1 can be represented by F k for each 2 ≤ k ≤ N . Using the recursive relations repeatedly, we find that F k can be represented by F N as In particular, F 1 can be represented by F N as Inserting this equation into the first equality of Eq. (4), we obtain This is a functional equation satisfied by F N . Applying Eq. (6) repeatedly, we obtain Taking n → ∞ in the above equation yields On the other hand, summing over all the equalities in Eq. (5), we obtain where Combining Eqs. (8) and (9), we obtain Since F (∞) = 0, integrating the above equation yields Finally, using the fact that F (0) = 1, we obtain the explicit expression of the Laplace transform, which is given by where is a normalization constant. Taking the inverse Laplace transform of F (λ) gives the probability density functionp(y) of the αth power of cell size. Finally, the probability density function p(x) of the original cell size is given by Next we focus on how to compute the cell size distribution numerically. Taking the derivative with respect to λ on both sides of Eq. (11) yields Replacing λ by iλ in the above equation yields This shows that the Fourier transform of yp(y) is exactly G(iλ). Since the Fourier transform and inverse fourier transform are inverses of each other, we only need to take the inverse Fourier transform of G(iλ) to obtain yp(y). Finally, we use Eq. (12) to obtain the cell size distribution p(x).

B.2 Model with deterministic partitioning: limiting case
We next focus on the special case where the cell cycle variability is very small (N 1). In this limit, the function b(λ) reduces to and the function f (λ) reduces to where r 0 = N 0 /N is the proportion of the elongation phase and r 1 = N 1 /N is the proportion of the reshaping phase. This shows that Combining Eqs. (13) and (14), it is easy to check that where v b , v m , and v d are constants defined as Thus the Laplace transform F (λ) can be simplified as is the exponential integral and is a normalization constant. The normalization constant can be calculated explicitly as where T 1 , T 2 , and T 3 are constants defined as Taking the inverse Laplace transform finally gives the distribution of the αth power of cell sizẽ where I A (x) is the indicator function which takes the value of 1 when x ∈ A and takes the value of 0 otherwise, δ(x) is Dirac's delta function, and w 1 , w 2 , and w 3 are constants defined as Finally, the distribution of cell size is given by where w 1 , w 2 , and w 3 represent the proportions of subpopulations in the elongation, septation, and reshaping phases, respectively and v b , v m , and v d represent the typical birth, septation, and division sizes, respectively. The analytical distribution can be used to compute some other quantities of interest. For example, when N 1, the mean cell size is given by

B.3 Model with stochastic partitioning
Next we compute the cell size distribution of lineage measurements for the model with stochastic partitioning. In this case, the evolution of cell size dynamics is governed by the master equation For convenience, we next focus on the αth power of cell size, y = x α . Letp k (y) denote the probability density function of the αth power of cell size when the cell is in stage k. Then the master equation can be rewritten as Using the Laplace transform defined in Eq. (3), the above equations can be converted to the following differential equations: At the steady state, the above equations reduce to where A 0 = αg 0 /a and A 1 = αg 1 /a. In analogy to the derivation for model I, we have where b(λ) is the function given in Eq. (7). This is a functional integral equation satisfied by F N . To solve this, we expand both F N (λ) and b(λ) in power series as Using the functional form of b(λ) in Eq. (7), it is not hard to check that where C n,m = n!/m!(n − m)! is the combinatorial number and (x) m = x(x + 1) · · · (x + m − 1) is the Pochhammer symbol. Inserting Eq. (17) into Eq. (16), we obtain where Comparing the coefficients on both sides of Eq. (19), we obtain This can be rewritten as where a n can be determined using the following recursive relations: where b n is the sequence defined in Eq. (18) and c n is the sequence defined in Eq. (20). To summarize, we obtain where a n is the sequence defined in Eq. (21). On the other hand, in analogy to the derivation for model I, we have where f (λ) is the function given in Eq. (10). Combining Eqs. (22) and (23), we obtain ∞ n=0 a n λ n .
Since F (∞) = 0, integrating the above equation yields ∞ n=0 a n u n du.
Finally, using the fact that F (0) = 1, we obtain the explicit expression of the Laplace transform, which is given by ∞ n=0 a n u n du. where is a normalization constant. Taking the inverse Laplace transform of F (λ) gives the probability density functionp(y) of the αth power of cell size. Finally, the probability density function p(x) of the original cell size is given by For the special case of exponential growth of cell size, there is only the elongation phase and the septation and reshaping phases vanish. In this case, we have N 1 = 0 and N = N 0 ; the cell size distribution is still determined by Eq. (24) with the sequence b n and the function f (λ) being simplified as

C Birth size distribution
Here we compute the distribution of V b . To this end, let V b (k) and V d (k) denote the cell sizes at birth and at division in the kth generation, respectively. Under the assumption of deterministic partitioning, we have V b (k + 1) = pV d (k) and thus we obtain the recursive equation is the generalized added size in the kth generation. Using the recursive equation repeatedly, we obtain Recall that ∆ 0 , ∆ 1 , ∆ 2 , · · · are i.i.d. hypoexponentially distributed random variables with the Laplace transform of each term being given by It thus follows from (25) and the independence of Since the distribution of V b (k) converges to the steady-state distribution of the birth size as k → ∞, taking k → ∞ in Eq. (26) shows that the Laplace transform of V α b is given by Taking the inverse Laplace transform gives the probability density function of V α b , from which is the probability density function of V b can be obtained.

D Correlation between birth and division sizes D.1 Model with deterministic partitioning
Let V b and V d denote the cell sizes at birth and at division in a particular generation, respectively, and let V b and V d denote the the birth and division sizes in the next generation, respectively. We first focus on the correlation between the birth size V d and the division size V d for the model with deterministic partitioning (model I). Since the generalized added size as well as where Cov(X, Y ) denotes the covariance between random variables X and Y and Var(X) denotes the variance of X. This shows that where ρ(X, Y ) denotes the covariance between X and Y . Since the generalized added size ∆ is the independent sum of an Erlang distributed random variable with shape parameter N 0 and mean M 0 and another Erlang distributed random variable with shape parameter N 1 and mean M 1 , we have This shows that p α (V α b + ∆) and V α b have the same distribution. Thus we obtain It then follows from Eq. (35) that Similarly, it follows from Eq. (36) that Combining Eqs. (37) and (38) shows that Inserting this equation into Eq. (28) finally shows that We next focus on the correlation between two successive birth sizes and the correlation between two successive division sizes. Since V b = pV d , the correlation coefficient between V α b and V α b is exactly the same as that between V α b and V α d , i.e.
Finally, since V b = pV d , the correlation coefficient between V α d and V α d is exactly the same as that between V α b and V α d , i.e.

D.2 Model with stochastic partitioning
We next focus on the correlation between birth and division sizes for the model with stochastic partitioning (model II). In this case, Eqs. (28) and (29) still hold and thus the key is to compute the variance of V α b . Let R = V b /V d denote the partition ratio. Since V b = RV d , we have This shows that R α (V α b + ∆) and V α b have the same distribution. Thus we obtain where we have used the fact that the partition ratio R is independent of the birth size V b and generalized added size ∆. Since R has a beta distribution with mean p and sample size parameter ν, we have It then follows from Eq. (35) that where K 1 = ER α 1 − ER α = B(α + pν, qν) B(pν, qν) − B(α + pν, qν) .
As ν → ∞, model II reduces to model I. In this case, we have Using these two equations, it is easy to check that (2K 1 + 1)K 2 − K 2 1 = 0 and thus Therefore, the correlation coefficient given in Eq. (39) for model II reduces to the one given in Eq. (34) for model I.