Exactly solvable model of gene expression in proliferating bacterial cell population with stochastic protein bursts and protein partitioning

TL;DR

This paper introduces an exactly solvable stochastic model of gene expression in dividing bacterial cells, accounting for stochastic protein bursts and partitioning, providing new analytical distributions and insights into noise sources.

Contribution

The model explicitly incorporates stochastic protein bursts and partitioning during cell division, offering exact solutions and improved understanding over classical decay-based models.

Findings

Analytical steady-state distributions similar to gamma distributions.

Explicit formulas for protein number variance and noise.

Differences from classical models in burst size and frequency estimates.

Abstract

Many of the existing stochastic models of gene expression contain the first-order decay reaction term that may describe active protein degradation or dilution. If the model variable is interpreted as the molecule number, and not concentration, the decay term may also approximate the loss of protein molecules due to cell division as a continuous degradation process. The seminal model of that kind leads to gamma distributions of protein levels, whose parameters are defined by the mean frequency of protein bursts and mean burst size. However, such models (whether interpreted in terms of molecule numbers or concentrations) do not correctly account for the noise due to protein partitioning between daughter cells. We present an exactly solvable stochastic model of gene expression in cells dividing at random times, where we assume description in terms of molecule numbers with a constant mean protein burst size. The model is based on a population balance equation supplemented with protein production in random bursts. If protein molecules are partitioned equally between daughter cells, we obtain at steady state the analytical expressions for probability distributions similar in shape to gamma distribution, yet with quite different values of mean burst size and mean burst frequency than would result from fitting of the classical continuous-decay model to these distributions. For random partitioning of protein molecules between daughter cells, we obtain the moment equations for the protein number distribution and thus the analytical formulas for the squared coefficient of variation.

Tables

Table 1. Table 1: Parameters of the distributions p n ​ ( x ; ω ​ ( n ) , b ) subscript 𝑝 𝑛 𝑥 𝜔 𝑛 𝑏 p_{n}(x;\omega{(n)},b) and p γ ​ ( x ; α , β ) subscript 𝑝 𝛾 𝑥 𝛼 𝛽 p_{\gamma}(x;\alpha,\beta) shown in Fig. 4 .

$n$	2	3	4
$\omega$	3	7	15
$\alpha$	1.90	2.52	2.88
$\beta /b$	3.07	5.28	9.74

Table 2. Table 2: Simulation parameter values for Fig. 2 in the main text. These values were the same in all simulations.

${\Delta }_{0\left(B\right)}$	1
$k_p$	250
$k_{dm}$	50
max_simulation_time	20000
time_cutoff	2000

Table 3. Table 3: Simulation parameter values for Fig. 2 in the main text. These values correspond to the specific curves in Fig. 2.

$n$	1	2	3	4
$k$	1	3	7	15

Equations

$$\frac{\partial p(x,t)}{\partial t}$$

$$-\gamma \frac{d}{dx}\left[xp(x)\right]$$

$$\gamma s\frac{dG(s)}{ds}=\Delta {\int }_0^1\eta (q)\left[G(qs)-G(s)\right]dq+k\hat{w}(s)G(s),$$

$$\left[(1-{\mathcal{M}}_r)\Delta +r\gamma \right]{\mu }_r=k\sum _{l=1}^r\left(\genfrac{}{}{0ex}{}{r}{l}\right){\mu }_{r-l}{m}_l,$$

$${\mu }_1=a{m}_1\equiv {\mu }_{1,\gamma },{\kappa }_2={\mu }_2-{\mu }_1^2=\frac{1}{2}a{m}_2\equiv {\kappa }_{2,\gamma },$$

$${\mu }_1=2\omega m_1\equiv {\mu }_{1,\Delta },{\kappa }_2=\frac{4{\omega }^2{\mathcal{M}}_2{m}_1^2+\omega m_2}{(1-{\mathcal{M}}_2)}\equiv {\kappa }_{2,\Delta },$$

$$c_v^2({\mu }_1)\equiv \frac{{\kappa }_{2,\Delta }}{{\mu }_{1,\Delta }^2}=\frac{m_2}{2m_1}\frac{1}{(1-{\mathcal{M}}_2)}\frac{1}{{\mu }_1}+\frac{{\mathcal{M}}_2}{(1-{\mathcal{M}}_2)}.$$

$$c_v^2({\mu }_1)=\frac{{\mathcal{E}}_1}{{\mu }_1}+{\mathcal{E}}_0,$$

$$\eta (q)={\eta }_d(q)\equiv \delta \left(q-1/2\right),$$

$$\nu (u)$$

$$\hat{\nu }(s)\equiv \mathcal{L}[\nu (u)]$$

$$c_v^2({\mu }_1)={\mathcal{E}}_0=\frac{{\mathcal{M}}_2}{1-{\mathcal{M}}_2}=const.$$

$$p_{\gamma }(x;\alpha ,\beta )\equiv \frac{x^{\alpha -1}\exp \left(-x/\beta \right)}{{\beta }^{\alpha }\Gamma (\alpha )}.$$

$$\alpha =\frac{{\mu }_1^2}{{\kappa }_2}=\frac{\omega \left(1-{\mathcal{M}}_2\right)}{\omega {\mathcal{M}}_2+1/2},\beta =\frac{{\kappa }_2}{{\mu }_1}=\frac{\left(2{\mathcal{M}}_2\omega +1\right)b}{1-{\mathcal{M}}_2}.$$

$$\beta (\alpha )$$

$$G(s)$$

$$\frac{(az;q{)}_{\infty }}{(z;q{)}_{\infty }}=\sum _{n=0}^{\infty }\frac{(a;q{)}_n}{(q;q{)}_n}{z}^n.$$

$${\kappa }_r=\frac{2^r{b}^r(r-1)![(1+\omega {)}^r-1]}{(2^r-1)}.$$

$${\mathcal{L}}^{-1}\{G(s)\}=p(x)=\sum _{i=0}^{\infty }\frac{C_i(\omega )}{b}\exp \left(\frac{-2^{i}x}{b(\omega +1)}\right),$$

$$C_i(\omega )=\frac{2^i\omega }{(1+\omega {)}^2}\frac{(2(1+\omega {)}^{-1};2{)}_i}{(2;2{)}_i}\frac{(\frac{1}{2}(1+\omega {)}^{-1};\frac{1}{2}{)}_{\infty }}{(\frac{1}{2};\frac{1}{2}{)}_{\infty }}.$$

$$p(x)\equiv p_n(x)=\frac{1}{b}\sum _{l=1}^n{A}_{n,l}\exp \left(-\frac{x}{2^{l}b}\right),$$

$$A_{n,l}=\frac{(-1{)}^{n-l}{2}^{\frac{l(l-3)}{2}}}{{\prod }_{i=1}^{l-1}(2^i-1){\prod }_{j=1}^{n-l}(2^j-1)}.$$

$${\mathcal{I}}_{12}={\int }_{x_1}^{x_2}\frac{k{h}_p^{\prime }(x)+h_d^{\prime }(x)+h_d(x)/b+R/b}{k{h}_p(x)+h_d(x)+R}dx,$$

$$\frac{\partial F(x,t)}{\partial t}$$

$$p(x,t)=\frac{F(x,t)}{\mathcal{N}(t)},$$

$$\mathcal{N}(t)={\int }_0^{\infty }F(x,t)dx$$

$$\frac{d\mathcal{N}(t)}{dt}=\mathcal{N}(t){\int }_0^{\infty }{h}_d(x)p(x,t)dx.$$

$$\frac{\partial F(x,t)}{\partial t}=\mathcal{N}(t)\frac{\partial p(x,t)}{\partial t}+p(x,t)\frac{d\mathcal{N}(t)}{dt}.$$

$$\mathcal{N}(t)={\mathcal{N}}_0{e}^{{\nu }_{m}t},$$

$${\nu }_m={\int }_0^{\infty }{h}_d(x)p(x)dx$$

$$\frac{\partial p(x,t)}{\partial t}$$

Full text

We present an exactly solvable stochastic model of gene expression in cells dividing at random times, where we assume description in terms of molecule numbers with a constant mean protein burst size. The model is based on a population balance equation supplemented with protein production in random bursts.

If protein molecules are partitioned equally between daughter cells, we obtain at steady state the analytical expressions for probability distributions similar in shape to gamma distribution, yet with quite different values of mean burst size and mean burst frequency than would result from fitting of the classical continuous-decay model to these distributions. For random partitioning of protein molecules between daughter cells, we obtain the moment equations for the protein number distribution and thus the analytical formulas for the squared coefficient of variation.

I Introduction

Genetically identical cells in a proliferating microbial population may differ not only in their cell-cycle age, volume or growth rate, but also in the copy numbers of molecules, in particular, the protein molecules. Even if the dependence of gene expression on growth rate Hintsche and Klumpp (2013); Klumpp and Hwa (2014); Keren et al. (2015); Shahrezaei and Marguerat (2015) and cell-cycle stage Walker et al. (2016) can be neglected, there still remain at least two factors causing fluctuations in the molecule numbers between individual cells: First, gene expression is an inherently stochastic process because of the small number of molecules involved in biochemical reactions. Second, molecules, including proteins, are partitioned between daughter cells at cell division in a more or less random manner Huh and Paulsson (2011a, b).

Yet, many stochastic models of gene expression neglect the protein partitioning Marathe et al. (2012). The classical models with a first-order decay reaction term Paulsson (2005); Cai et al. (2006); Friedman et al. (2006); Taniguchi et al. (2010) may be interpreted in two ways: (i) If the model describes the protein levels in terms of molecule numbers, then the degradation term describes not only a possible active protein degradation (rare in bacteria Maurizi (1992)) but it also mimicks the effect of cell division. However, this description neglects the fluctuations both due to the abrupt protein loss at cell division (which are present even at the half-by-half protein partitioning) and due to a possible random partitioning of the molecules between cells. (ii) If the same classical model describes the protein levels in terms of concentrations, then the continuous degradation term reflects, besides possible active degradation, the dilution due to cell growth. Under the assumption of a perfectly equal protein partitioning, the concentrations remain unchanged at cell division. However, if protein partitioning is random, then it should result in concentration fluctuations (at binomial partitioning, they extinct at high gene expression levels), and these fluctuations are neglected in the conventional model.

Based on the classical models of gene expression where proteins are assumed to be produced in random translational bursts Paulsson (2005); Cai et al. (2006); Friedman et al. (2006); Taniguchi et al. (2010) and on the population balance equations Mantzaris (2005, 2006); Friedlander and Brenner (2008), we propose a minimal analytically tractable model of gene expression in a population of dividing cells. This model directly accounts for the protein loss at cell division and random partitioning of the molecules between daughter cells.

The manuscript is organized as follows: Sec. II presents the theoretical model. In Sec. III.1 we calculate the moments of the steady-state protein number distribution for any symmetric random partitioning of proteins between daughter cells. In Sec. III.2 we show that a noise floor, i.e. the absolute lower bound for noise as a function of mean protein number, occuring for highly expressed proteins, arises in the model as a consequence of protein loss at cell division. This kind of a noise floor is a property of the description in terms of molecule numbers. In Sec. III.3 we compare the protein number distributions resulting from our model to the distributions from the seminal Friedman’s model Friedman et al. (2006). In Sec. III.3.1 we show that if both distributions have the same mean and variance then the underlying mean burst frequency and mean burst size may differ. In Sec. III.3.2 we obtain the analytical form of the protein number distribution for the half-by-half protein partitioning and we show that its shape is very similar to the shape of a gamma distribution from the Friedman’s model Friedman et al. (2006). In Sec. III.4 we analyze the behavior of the distribution tails for large protein numbers. In Sec. III.5 we compare the squared coefficient of variation of our model to the existing experimental data. Section IV contains discussion and conclusions.

II Model

We combine the analytical framework first introduced by Friedman et al. Friedman et al. (2006) with the population balance equation (PBE) Mantzaris (2005, 2006); Friedlander and Brenner (2008), a standard tool for modelling microbial populations.

In Ref. Friedman et al. (2006), gene expression was described as a compound Poisson process with random Poissonian events of protein production, where a random number of protein molecules was produced at once in each such event. The protein burst size $u$ was drawn from a distribution $\nu(u)$. Originally, $\nu(u)$ was assumed to be an exponential distribution Friedman et al. (2006), based on the experimental data for E. coli Yu et al. (2006); Cai et al. (2006), but here we will relax this assumption. Protein production was counterbalanced in the Friedman’s model Friedman et al. (2006) by a continuous, deterministic protein decay term. The model allowed self-regulation: The burst occurrence rate $k$ was modulated by a transfer function $h_{p}(x)$, dependent on the current protein level $x$. We note that the continuous variable $x$ was assumed to be protein concentration in Ref. Friedman et al. (2006). However, in this study, we interpret $x$ as the molecule number. This distinction is crucial because cell division exactly by half maintains a constant protein concentration, whereas the protein number drops by 1/2. If protein production bursts are constant in terms of molecule number, then their size in the units of concentration decreases in growing cells.

Therefore, we treat here the theoretical framework of Ref. Friedman et al. (2006) as a continuous approximation to discrete fixed cell size gene expression models, cf. e.g. Shahrezaei and Swain (2008), and not in its original interpretation, as a model describing evolution of protein concentration in growing and dividing cell. Within the model proposed in Ref. Friedman et al. (2006), the first order decay term describes both true protein degradation (if present) and protein dilution due to cell growth.

On the other hand, population balance equation (PBE) describes the time evolution of the cell number density in a proliferating microbial population Mantzaris (2005, 2006). At cell division, protein molecules are partitioned according to $x\to\{qx,(1-q)x\}$, where $x\in[0,\infty)$ is a continuous protein copy number, and not protein concentration, whereas the division ratio $0\leq q\leq 1$ is a random number, drawn from the division ratio probability density function (PDF) $\eta(q)=\eta(1-q)$. Within the PBE approach as used in Mantzaris (2005, 2006); Friedlander and Brenner (2008), protein production is described as a deterministic process. However, both protein production rate and cell division rate may depend on $x$ in a nontrvial way.

By merging both the above theoretical frameworks, we get the time-evolution equation for $p(x,t)$, a PDF of the protein number $x$ at time $t$ in the population of cells (Appendix A):

[TABLE]

In the above equation, $x-x^{\prime}$ is the burst size, $\nu(x-x^{\prime})$ is the burst size PDF, $w(x-x^{\prime})\equiv\nu(x-x^{\prime})-\delta(x-x^{\prime})$ and $\delta(x-x^{\prime})$ is Dirac delta distribution Friedman et al. (2006). The $x$-dependence of cell division and protein production rates are given by $h_{d}(x)$ and $h_{p}(x)$, respectively. Thus, $kh_{p}(x)$ is the burst frequency, whereas $g(x)$ describes either deterministic protein production or degradation. Eq. (1) is a generalization of both Eq. (1) of Ref. Friedman et al. (2006) and Eq. (1) of Ref. Friedlander and Brenner (2008).

If $x$ cannot be treated as proportional to the cell mass or volume, it is unclear how the division rate should depend on the copy number of a given protein, and various functional forms of $h_{d}(x)$ were used in the literature Mantzaris (2005, 2006); Friedlander and Brenner (2008). For that reason, and in order to obtain an analytically tractable model, we focus here on the simplest case of the constant cell division rate, $2h_{d}(x)\equiv\Delta=const$, which implies $\int_{0}^{\infty}h_{d}(\xi)p(\xi)d\xi=\Delta/2$.

III Results

III.1 Moments of the steady-state protein number distribution

We consider the steady-state solution of Eq. (1), $\partial p(x,t)/\partial t=0$. Assuming $g(x)=-\gamma x$, $\gamma\geq 0$ and unregulated protein production ($h_{p}(x)=1$), we have

[TABLE]

The Laplace transform ($\mathcal{L}\{\ldots\}$) of Eq. (2) yields

[TABLE]

where $G(s)=\mathcal{L}\{p(x)\}$, and $\hat{w}(s)=\mathcal{L}\{w(u)\}=\hat{\nu}(s)-1$. Moment equations follow from Eq. (3),

[TABLE]

where $\mu_{r}=\int_{0}^{\infty}x^{r}p(x)dx$ is the $r$-th moment of the protein number PDF, $m_{r}=\int_{0}^{\infty}u^{r}\nu(u)du$ is the $r$-th moment of the burst size PDF, and $\mathcal{M}_{r}=\int_{0}^{1}q^{r}\eta(q)dq$ is the $r$-th moment of the division ratio PDF (note that $\mathcal{M}_{1}=1/2$).

For $\gamma\neq 0$, $\Delta=0$ (the proteins are degraded, but there is no protein partitioning at cell division) Eq. (4) yields

[TABLE]

where $a=k/\gamma$ Jedrak and Ochab-Marcinek (2016). On the other hand, for $\gamma=0$ (the proteins are stable) and $\Delta\neq 0$, from Eq. (4) we get

[TABLE]

where $\omega=k/\Delta$. The variance of the protein number, $\kappa_{2,\Delta}=\text{var}[p(x)]$ (6) is a monotonically increasing function of $(\mathcal{M}_{2}-\mathcal{M}_{1}^{2})/\mathcal{M}_{1}^{2}=4\mathcal{M}_{2}-1$, i.e., of the noise due to protein partitioning.

For a given $\nu(u)$ and $\eta(q)$, if $a=2\omega$, i.e., if $\mu_{1,\gamma}$ (5) is equal to $\mu_{1,\Delta}$ (6), then from Eqs. (5) and (6) it follows that $\kappa_{2,\Delta}>\kappa_{2,\gamma}$. Therefore, within the present description in terms of molecule numbers, the first-order decay such as in the classical model by Friedman et al. Friedman et al. (2006) is not equivalent to any protein partitioning mechanism.

III.2 Squared coefficient of variation has a noise floor due to cell division

Protein degradation is important for eukaryotic organisms, but may be neglected in the case of bacteria; hence we put $\gamma=0$. From Eq. (6) we obtain the squared coefficient of variation

[TABLE]

The protein number noise, measured by the squared coefficient of variation (7) is thus of the form

[TABLE]

with $\mu_{1}$-independent $\mathcal{E}_{0}$ and $\mathcal{E}_{1}$, which leads to the characteristic boomerang-like shape of the log-log plot (Fig. 1A). In Fig. 1, the half-by-half protein partitioning is assumed,

[TABLE]

which does not depend on $x$ and which gives $\mathcal{M}_{2}=1/4$. However, for small $x$, the half-by-half division as given by (9) is much less realistic than other continuous approximations to discrete binomial distribution, such as the ’all or none’ partitioning ($\eta(q)=\eta_{b}(q)\equiv\frac{1}{2}\left[\delta(q)+\delta(1-q)\right]$) or even a uniform distribution of partition ratio ($\eta(q)=\eta_{u}(q)\equiv\Theta(q)\Theta\left(1-q\right)$), considered in Refs. Brenner and Shokef (2007); Friedlander and Brenner (2008).

$\mathcal{E}_{0}=\mathcal{M}_{2}/(1-\mathcal{M}_{2})$ provides the lower bound for extrinsic noise and does not depend on the details of protein production mechanism. However, the contribution of protein partitioning to the value of $\mathcal{E}_{0}$ predicted by the present model for a realistic level of protein partitioning noise is likely to be overestimated for the following reason: Eq. (2) is a Master equation with time-independent parameters, hence it describes the cell division events as a Poisson process Gardiner (2009). In consequence, the cell cycle length has an unrealistic, exponential distribution, contributing with too much noise (see Appendix C). Therefore, the age structure of the population is also unrealistic (Appendix D). This indicates a limitation of the simplest, single-variable version of the PBE approach as a candidate for a minimal model of gene expression as it assumes the "most random" distribution of cell cycle lengths (i.e., the one that naturally occurs from a Poisson process).

Thus, our model provides an "upper bound" for the noise contribution due to cell cycle length fluctuations (meaning the "most random" fluctuations that can be generated by a model based on Master equation Jafarpour et al. (2018)): In the case when $\eta(q)$ is given by Eq. (9) there is no protein partitioning noise ($\text{var}[\delta\left(q-\frac{1}{2}\right)]=0$) and $\mathcal{E}_{0}$ has the minimal value $\mathcal{E}_{0}=1/3\equiv\mathcal{E}_{0,d}$. Therefore, $\mathcal{E}_{0,d}$ is the contribution to the noise floor coming solely from fluctuations of the cell cycle length, whereas the remaining part, $\mathcal{E}_{0}-\mathcal{E}_{0,d}=(4\mathcal{M}_{2}-1)/[3(1-\mathcal{M}_{2})]$ is due to protein partitioning noise.

III.3 Our model vs. Friedman’s model Friedman et al. (2006): Similar distributions but different parameters

III.3.1 Mean burst frequency and mean burst size differ between the models despite the same mean and variance of protein number distribution

Let us consider a single-parameter burst size PDF, $\nu(u;b)=h\left(u/b\right)/b$, where $u\equiv x-x^{\prime}$ is a burst size, such that $m_{1}=b$ and $m_{2}=Cb^{2}$. $C=2$ for the exponential burst distribution:

[TABLE]

For $m_{2}/m_{1}=Cb=const$, the increase in mean protein number $\mu_{1}=2\omega b$ is caused solely by the increase in $\omega$ (i. e., by the increase in burst frequency $k$ or decrease in cell division rate $\Delta$). In such a case, we move along the boomerang-shaped trajectories for a fixed value of $b$, like those depicted in Fig. 1A. In general, $\mu_{1}$ may be varied both due to the change in $b$ and $\omega$.

If $k=0$ and $g(x)=-\gamma x$ term in Eq. (2) is replaced with $g(x)=\sigma>0$ (this is a special case of the model of Ref. Friedlander and Brenner (2008)), then instead of (7) we obtain

[TABLE]

In other words, the squared coefficient of variation is then described by some general equation in the form of (8) (but not (7)), where $\mathcal{E}_{0}\neq 0$, $\mathcal{E}_{1}=0$ if protein production in our model is deterministic and cell division is stochastic. The opposite situation, $\mathcal{E}_{1}\neq 0$, $\mathcal{E}_{0}=0$, is encountered for $\gamma\neq 0$, $\Delta=0$, in particular for the model of unregulated, bursty gene expression with continuous protein degradation Friedman et al. (2006), which predicts that protein numbers are gamma-distributed:

[TABLE]

Then, $c_{v}^{2}(\mu_{1})=1/\alpha=\beta/\mu_{1}$, i. e., $\mathcal{E}_{1}=\beta$ and $\mathcal{E}_{0}=0$. To make a comparison of $p_{\gamma}(x;\alpha,\beta)$ (12) to the protein number PDF of the present model, $p(x)$, we assume now that these two distributions have equal means and equal second moments, $\mu_{r,\gamma}=\mu_{r,\Delta}$, $r=1,2$. We also assume that the burst size PDF $\nu(u)$ is given by (10), but the mean burst size may be different for the two models, $m_{1,\gamma}=\beta\neq b=m_{1,\Delta}$. These assumptions yield the following mapping:

[TABLE]

From (13) we obtain

[TABLE]

For simplicity, we assume now that protein partitioning is deterministic, i.e., that each daughter cell obtains exactly half of the total number of protein molecules at cell division. In such a case, $\eta(q)=\eta_{d}(q)$ is given by Eq. (9), so that $\mathcal{M}_{2}=1/4$. In Fig. 1B we plot $\alpha$ vs. $\beta$ (14) to show that the mean burst frequency $\omega$ and mean burst size $b$ in our model take different values than the corresponding parameters $\alpha$ and $\beta$ of the gamma distribution (12). In particular, $\omega$ is not limited to several bursts per cell cycle as was $\alpha$, so the mean burst size $b$ does not need to take as high values as would $\beta$ need to take in order to obtain a high level of gene expression.

III.3.2 Analytical form of the protein number distribution for half-by-half protein partitioning. Apparent similarity to gamma distribution.

Below, we will show that equating the two first moments of our PDF $p(x)$ with those of the gamma distribution $p_{\gamma}(x)$ (12) yields similar overall shapes of the distributions. If so, any experimental gamma-shaped distributions with $c^{2}_{v}\geq\mathcal{E}_{0,d}=1/3$ may be as well fitted by the distributions given by our model (under the assumption that protein numbers, and not concentrations, were measured in experiment). With $\eta(q)=\delta(q-\frac{1}{2})$, we can rewrite Eq. (3) as $G(s)=R(s)G\left(s/2\right)$, where $R(s)\equiv\{1-\omega[\hat{\nu}(s)-1]\}^{-1}$, $R(0)=1$. Solving by iteration, we get $G(s)=\prod_{i=0}^{\infty}R\left(s/2^{i}\right)$. For $\hat{\nu}(s)$ (10), $G(s)$ reads

[TABLE]

where $(a;q)_{k}$ is a $q$-Pochhammer symbol. In (15) we used the $q$-binomial theorem Andrews (1986),

[TABLE]

The symbol $q$ appering in Eq. (16) should not be confused with protein partitioning ratio $q$ appearing e.g. in Eqs. (1), (2), (3) or (9). The letter $q$ is traditionally used in the branch of mathematics called $q$-theory or $q$-analogs ($q$-binomial theorem or $q$-Pochhammer symbol are examples of such $q$-analogs).

From Eq. (15) we obtain cumulants of $\mathcal{L}^{-1}\{G(s)\}$,

[TABLE]

We also have

[TABLE]

where

[TABLE]

For $\omega=\omega_{n}\equiv 2^{n}-1$, $n=1,2,3,\ldots$, $p(x)$ (18) can be written as a finite series

[TABLE]

where

[TABLE]

Equations (18) and (19) define a two-parameter family of PDFs, resembling the gamma PDF (Fig. 2). In particular, $p(x)$ (18) is right-skewed, unimodal for $\omega>1$, and monotonically decreasing for $\omega\leq 1$.

III.4 Distribution tail for large protein numbers

In both $p(x)$ (18) and $p_{\gamma}(x)$ (12), the mean burst size ($b$ and $\beta$, respectively) is the scaling parameter. The tail of $p(x)$ (18) is exponential: For large $x$ we have $p(x)\sim\exp[-x/(b(\omega+1))]$. However, in contrast to the gamma distribution, where the leading term is $\exp[-x/b]$, the exponent in our model depends not only on the mean burst size, but also on burst frequency.

The same asymptotic behaviour as for $p(x)$ (18) is present if instead of $\eta(q)=\delta(q-\frac{1}{2})$ any other $\eta(q)$ is used (except for some pathological cases, e.g. for $\eta_{b}(q)=\frac{1}{2}\left[\delta(q)+\delta(1-q)\right]$) (Appendix F). This is the special case of a yet more general result. Using Eq. (1), it can be shown in a similar manner as in Ref. Friedlander and Brenner (2008) that for $g(x)=0$ and $\nu(u)$ (10) the ratio $p(x_{2})/p(x_{1})$ at the tail of the PDF is well approximated by $\exp\left(-\mathcal{I}_{12}\right)$, where

[TABLE]

and $R=\int_{0}^{\infty}h_{d}(\xi)p(\xi)d\xi$ (Appendix F). However, even if the tails are universal, i.e., independent on the division-ratio PDF $\eta(q)$, the corresponding probability distributions are not. For example, for $\eta(q)=\Theta(q)\Theta\left(1-q\right)$ used in Refs. Friedlander and Brenner (2008) and Brenner and Shokef (2007), instead of (18) we obtain a statistical mixture of two broad ($c^{2}_{v}\geq 1/2$) gamma distributions (Appendix G). Therefore in the present case it is no longer true that the “steady-state population distribution (…) becomes insensitive to the division details” Brenner and Shokef (2007). For large $x$, $\eta(q)=\Theta(q)\Theta\left(1-q\right)$ is much less realistic than $\eta(q)=\eta_{d}(q)=\delta(q-\frac{1}{2})$ (9) used here, and it leads to a higher noise floor ($\mathcal{E}_{0}=1/2$).

III.5 Comparison to experimental data

In this section, we compare the squared coefficient of variation as a function of the mean protein level in our model to the existing experimental data, under the assumption of equal protein partitioning between daughter cells (the lowest possible noise due to cell division, Fig. 1A).

In Ref. Nordholt et al. (2017), the authors measured total cell fluorescence emitted by the green fluorescent protein (GFP), encoded for by a gene controlled by the hyper-spank promoter in B. subtilis. The promoter activity was modulated by the concentration of isopropyl-D-thiogalactoside (IPTG). An independent parameter varied in the experiment was the cell growth rate, depending on the carbon source in the medium. The data points shown in Fig. 1A were obtained in Ref. Nordholt et al. (2017) by variation of these two parameters (we use the symbols and colors as in that Reference; we manually extracted the data points from Fig. 4a therein using the xyscan software).

Assuming that the total cell fluorescence scales linearly with the number of reporter proteins, the $x$ variable in our model may be reinterpreted as the fluorescence level and the mean protein burst size $m_{1}\equiv b$ is then expressed in fluorescence units. Thus, the mean protein number is proportional to the mean total cell fluorescence, and the squared coefficient of variation does not depend on the units in which gene expression was measured (molecule number or fluorescence units). We note that the linear scaling between the GFP molecule number and total cell fluorescence is not necessarily an obvious assumption (see, e.g., Ref. Newman et al. (2006)). On the other hand, the authors of Ref. Nordholt et al. (2017) estimated that GFP maturation time did not affect strongly the fluorescence level, which may possibly exclude one source of non-linear scaling of the fluorescence detected vs. molecule number. The linear scaling was also assumed in Ref. Wolf et al. (2015).

The range of the mean total cell fluorescence measured in Ref. Nordholt et al. (2017) was too narrow to show a noise floor (Fig. 1A). However, the data provide a hint that if the noise floor exists, it would be lower than that predicted by our model. This suggests that the distribution of cell cycle lengths in our model is indeed too wide, as we pointed out in Section III.2.

For a rough comparison, we also plotted a fit of the squared coefficient of variation $c_{v}^{2}$ of our model (7) to the data. The slope of $c_{v}^{2}$ in our model seems to be in agreement with the experimental results. The fit yields an estimation of the mean burst size $m_{1}\equiv b=18.0\pm 2.1$ in the units of fluorescence, as used in Ref. Nordholt et al. (2017); however, this value should be treated with caution because some other model with a lower noise floor due to a narrower distribution of cell-cycle lengths may possibly yield a different value of $b$.

For completeness, we also plotted the data from Ref. Nordholt et al. (2017) in Fig. 1B (we omitted the colors and symbols used in Fig. 1A). We made the assumption that the total cell fluorescence in the experiment was gamma-distributed, as shown in the SI of Ref. Nordholt et al. (2017). We mapped the mean $\langle f\rangle$ and squared coefficient of variation $c_{v,f}^{2}=var(f)/\langle f\rangle^{2}$ of the total cell fluorescence onto the corresponding values of parameters of gamma distributions: $\alpha=(c_{v,f}^{2})^{-1}$, $\beta=\langle f\rangle\cdot c_{v,f}^{2}$.

A noise floor was observed in other experiments reported in literature Volfson et al. (2006); Keren et al. (2015); Taniguchi et al. (2010); Silander et al. (2012); Wolf et al. (2015); Newman et al. (2006); Bar-Even et al. (2006) but most of these data were not suitable for comparison to our model (see an overview in Appendix E) because they were normalized to cell volume or gated to make them independent on cell-cycle stage. Thus, molecule number fluctuations due to loss of proteins at cell division were absent in these data.

There are two studies on S. cerevisiae that reported the gene expression noise floor in ungated measurements. Its levels lower than in our model (where $\mathcal{E}_{0}=1/3$): In Ref. Volfson et al. (2006), Fig. 2 and S2 therein, the noise floor measured as the coefficient of variation $c_{v}$ took the values between 0.3 and 0.4, corresponding to the squared coefficient of variation $c_{v}^{2}$ between 0.09 and 0.16. In Ref. Keren et al. (2015), Fig. 2A and S9 therein, $c_{v}^{2}$ was between 0.1 and 0.2. One reason for that difference may be the too wide distribution of cell cycle lengths in our model. But it should also be noted that our model may be far too simplistic for eukaryotic cells, as it does not account for the discrete promoter activity states due to chromatin remodeling, the nuclear transport, etc.

IV Discussion and conclusions

We have proposed a model of gene expression in a proliferating cell population, which is a generalization of the models of Refs. Friedman et al. (2006) and Ref. Friedlander and Brenner (2008). In Ref. Friedman et al. (2006), the protein production was stochastic, whereas the decrease of protein concentration was due to a deterministic protein degradation process. Here, we treat the model of Ref. Friedman et al. (2006) as a continuous approximation (still in the units of protein number) to discrete gene expression models, such as in Ref. Shahrezaei and Swain (2008), that describe non-growing and non-dividing cells.

In Ref. Friedlander and Brenner (2008), the situation was just opposite: Protein partitioning was stochastic and protein production was deterministic. After combining the stochastic protein production (of Ref. Friedman et al. (2006)) with the stochastic protein partitioning (used in Ref. Friedlander and Brenner (2008)), we obtain the boomerang-like shape of the log-log plot of mean protein copy number vs. protein copy number noise.

Because the protein copy number and not protein concentration is the variable in our model, the protein partitioning at cell division along with an age structure of the cell population (cf. Eq. (41)) lead to the existence of the noise floor – an absolute lower bound for noise, present for highly expressed proteins.

Our results suggest that the values of mean burst size and burst frequency that may be obtained by fitting theoretical distributions to experimental data are strongly model-dependent. Naïve fitting of the gamma distribution Friedman et al. (2006) or its discrete counterpart, negative binomial distribution Shahrezaei and Swain (2008), to the data measured in terms of molecule numbers, would neglect protein loss at cell division because the underlying models neglected cell division, and thus such a fitting might overestimate the mean burst sizes and underestimate mean burst frequencies for higher gene expression levels.

Our model, directly accounting for the protein molecule numbers, may be important for interpretation of experimental results. In some studies, data correction for cell size was carried out: In Taniguchi et al. (2010), cell volume was measured by image recognition and the protein fluorescence was normalized by the volume. In other studies, the flow cytometry data were gated Silander et al. (2012); Wolf et al. (2015); Newman et al. (2006); Bar-Even et al. (2006). Gating filters out the data only for those cells that scatter light to a similar extent, so the observed cells are supposed to be of similar sizes. However, the gating procedure may be imperfect because setting the size range too narrow leaves too little data, whereas setting the gate too wide results in some variation in cell sizes. One can avoid these problems by measuring the total cell fluorescence Nordholt et al. (2017); Volfson et al. (2006); Keren et al. (2015) and then separating the noise contributions from varying cell sizes (due to cell-cycle progression and other sources of variability) using independent measurements Nordholt et al. (2017) or estimating these contributions theoretically. Our model brings us closer to theoretical determination of the magnitude of one of these contributions, namely, the protein loss due to cell division and partitioning between daughter cells.

We have compared our model to the existing literature data Nordholt et al. (2017), however, the possibilities of such a comparison are limited so far. The results for B. subtilis, reported in Ref. Nordholt et al. (2017), were measured in too narrow a range of gene expression levels, whereas the results for S. cerevisiae of Refs. Volfson et al. (2006); Keren et al. (2015) may not be comparable to our minimal model because of the complexity of gene expression mechanisms in eukaryotes. However, a cautious comparison of our model to the above data suggests that the main limitation of the model is its prediction of too large a contribution of protein loss at cell division to the noise floor. This is due to the exponential distribution of cell cycle lengths, since the cell divisions are modelled as Poissonian events (Appendix C). In a more realistic model, the cell cycle length distribution could be modeled, instead of the exponential distribution, as gamma distribution Powell (1956) or some other distribution peaked around the mean cell cycle length, with the limiting case of Dirac delta distribution describing a deterministic cell cycle. However, this is beyond the scope of the present paper. In preparation is our new paper Jędrak and Ochab-Marcinek (2018) that explores the noise levels in a more realistic model.

Another limitation of our theoretical approach may be the description of the messenger RNA (mRNA) as very short-lived molecules, only implicitly present in the model Friedman et al. (2006), and thus the neglect of their partitioning at cell division. Also, the model does not describe the factors that may significantly shape the gene expression noise levels in eukaryotes, e. g., the discrete on-off promoter switching or nuclear transport.

We may expect that negative gene autoregulation (not considered here) would suppress protein noise Friedman et al. (2006); Jia et al. (2017), and thus decrease the noise floor. Yet, positive autoregulation is expected to have the opposite effect Friedman et al. (2006); Jia et al. (2017). In a similar manner, nontrivial dependence of cell division rate on protein number ($h_{d}(x)\neq const$) would probably result in a lower noise floor. Still, the analytical solution to our model in presence of gene autoregulation or cell division regulation does not seem feasible. For that reason, we have studied here the non-regulated gene expression and cell division only.

Acknowledgments

AOM was supported by the National Science Centre SONATA BIS 6 grant no. 2016/22/E/ST2/00558.

Author contributions

JJ designed the study, performed the analytical calculations and wrote the manuscript. MK wrote the simulation code for cell division. AOM participated in the study design by the idea of comparison of the model with experimental data and by making the comparison; designed the simulation algorithm for cell division, supervised writing the simulation code and introduced minor modifications in the code; performed the simulations; performed the distribution fitting by maxima matching and observed heuristically the behavior of the probablity distributions described in Appendix J; participated in writing the manuscript.

Appendix A Derivation of Eq. (1)

Equation (1) follows from an apparently very similar, but more fundamental equation describing the number density of cells, $F(x,t)$, in a proliferating population,

[TABLE]

Namely, $F(x,t)dx$ is the number of those cells in a population, which at time $t$ contain exactly $x$ protein molecules. In order to derive Eq. (1) from Eq. (23), we define

[TABLE]

where

[TABLE]

is the total number of cells in the population Mantzaris (2006). Integrating both sides of Eq. (23) from $x=0$ to $x=\infty$, assuming that $F(\infty,t)=g(0)F(0,t)=0$ and making use of (24) and (25) we obtain

[TABLE]

Now, Eq. (1) follows from (23)-(26) and from obvious identity

[TABLE]

The above derivation is essentially that of Ref. Mantzaris (2006), the only difference is the presence of the terms responsible for bursty protein production in the present case.

In spite of their apparent similarity, there are important differences between Eq. (1) and Eq. (23). First, the terms related to protein partitioning in Eq. (23) conserve the total number of protein molecules (molecules are neither created, nor destroyed during cell division), whereas from Eq. (1) it follows that the cell division always decreases the mean protein number in the population. This is due to the $(-1)p(x,t)\int_{0}^{\infty}h_{d}(\xi)p(\xi,t)d\xi$ term in Eq. (1), not present in Eq. (23).

Second, $\mathcal{N}(t)$ grows indefinitely according to Eq. (26), therefore in contrast to Eq. (1), Eq. (23) does not have nontrivial stationary solutions.

Our task now is to determine, for the present model: (i) the population growth rate, $\nu_{m}$, (ii) probability distribution $f(\tau)$ of generation time (cell cycle length), $\tau$, for each newborn cell, and (iii) the cell age ($\tilde{a}$) distribution $\phi(\tilde{a})$ in the state of a balanced, exponential growth 111In this section we use notation of Ref. Powell (1956), but the cell age is denoted $\tilde{a}$ and not $a$ to avoid confusion with quantity $a=k/\gamma$ defined in the main text.. This is done in next three sections.

Appendix B Population growth rate $\nu_{m}$

The steady state solutions of Eq. (1), $p(x)$, correspond to the so called state of balanced growth, when the shape of $F(x,t)=\mathcal{N}(t)p(x)$ does not change but is only rescaled by $\mathcal{N}(t)$ Mantzaris (2006). In such a case, from Eq. (26) we obtain

[TABLE]

where $\nu_{m}$, given by

[TABLE]

is the population growth rate. $\mathcal{N}(t)$ as given by (28) is characteristic for (in fact, defines) the phase of exponential growth Powell (1956)222We assume that the number of cells in the population, $\mathcal{N}(t)$, is sufficiently large its fluctuations may be neglected, and that $\mathcal{N}(t)$ may be regarded as evolving according to deterministic equation (28)..

If $g(x)=-\gamma x$ and if neither the protein production rate nor the cell division rate depend on the number of protein molecules ($h_{p}(x)=1$, $h_{d}(x)\equiv\Delta/2$), Eq. (1) reads

[TABLE]

Central to this paper is the steady-state limit of Eq. (30),

[TABLE]

i.e., Eq. (2). Although time-independent, Eq. (31) describes the stationary protein distribution $p(x)$ in a growing population of dividing cells (i.e., the state of balanced growth), as discussed above. In particular, for $h_{d}(x)=\Delta/2\equiv\Delta_{0}$, from Eq. (29) it follows that

[TABLE]

i.e., the population growth rate $\nu_{m}$ appearing in Eq. (28) is equal to individual cell division rate $\Delta_{0}$, as should be expected.

It is instructive to show this result in an alternative way. Consider time evolution of the moments of $p(x,t)$, resulting from Eq. (30),

[TABLE]

In the absence of protein production ($k=0$) and degradation ($\gamma=0$), from Eq. (33) it follows that the time dependence of the first moment of $p(x,t)$, i.e., the mean protein number is given by 333In such a case the only possible steady state is a trivial one, i.e., $\lim_{t\to\infty}p(x,t)=p(x)=\delta(x)$.

[TABLE]

On the other hand, from the definition of the population average, the mean protein number $\mu_{1}(t)$ is equal to $X(t)/\mathcal{N}(t)$, where $X(t)\equiv\sum_{i=1}^{\mathcal{N}(t)}x_{i}(t)$ denotes the total number of protein molecules in a population, and $x_{i}(t)$ is a number of protein molecules in $i$-th cell at time $t$. If there is no protein production or degradation, $X(t)$ is constant and equal to its initial value, $X(t)=X(0)=\mathcal{N}_{0}\mu_{1}(0)$. In such a case, time evolution of the mean protein number $\mu_{1}(t)$ is caused solely by the increase in cell number,

[TABLE]

Comparing (34) and (35) we see that $\nu_{m}=\Delta_{0}$, i.e., the population growth rate $\nu_{m}$ is equal to the individual cell division rate $\Delta_{0}$, as it should.

So far, we have considered the whole proliferating cell population, for which $\Delta=2\Delta_{0}=2\nu_{m}$ (scenario A). However, Eq. (30), but not Eq. (1), may be also used to describe the time evolution of the protein number distribution $p(x,t)$ in a single cell lineage. In such a case, we discard one of the daughter cells at each cell division, and therefore $\Delta=\Delta_{0}$ (scenario B).

Appendix C Generation time distribution $f(\tau)$

In order to find the generation time distribution $f(\tau)$, consider Eq. (30). This is a special case of the differential Chapman-Kolmogorov equation Gardiner (2009) with $x$- and $t$-independent coefficients. We assume once again that there is no protein production ($k=0$) and that the protein is stable ($\gamma=0$). In such a case, (30) becomes the Master equation, describing ’jump process’ between different states of the system (there is no drift term), and these ’jumps’ are solely due to the cell division events. The probability that the system does not undergo such a ’jump’, and that it is still in the same state at $t=\tau$ as it was at $t=0$, is equal to $\exp(-\Delta\tau)$ Gardiner (2009). Therefore, the probability of a jump occurring in the infinitesimal interval $(t,t+dt)$ is given by van Kampen (2007)

[TABLE]

In the case of a single lineage (scenario B mentioned above), we have $\Delta=\Delta_{0}$, and $\pi(t)$ must be identified with the cell cycle length distribution, $\pi(\tau)=f(\tau)$. Therefore, we obtain

[TABLE]

If we deal with the whole proliferating population, then $\Delta=2\Delta_{0}$ (scenario A), and $\pi(t)$ (36) reads

[TABLE]

$\pi(\tau)$ (38) is identical with the quantity defined as

[TABLE]

and called the ’carrier distribution’ in Powell (1956).

Appendix D Cell age distribution $\phi(\tilde{a})$

To convince ourselves that (37) is valid, and to find the explicit form of the age distribution $\phi(\tilde{a})$, assume that $f(\tau)=\lambda\exp(-\lambda t)$, for $\lambda$ yet unspecified. From Eq. (9) of Ref. Powell (1956) it follows that $\phi(\tilde{a})$ is given by

[TABLE]

The normalization of $\phi(\tilde{a})$ (i.e., the condition $\int_{0}^{\infty}\phi(\tilde{a})d\tilde{a}=1$) yields $\nu_{m}=\lambda$, whereas from Eq. (32) or from (34) and (35) we obtain $\lambda=\Delta_{0}$, i.e., we recover Eq. (37). Hence, the age distribution $\phi(\tilde{a})$ is given by

[TABLE]

It should be emphasized that functional forms of $f(\tau)$ (37), $\mathcal{C}(\tau)$ (39) and $\phi(\tilde{a})$ (41) are rather unrealistic. This is a consequence of a Poissonian nature of cell division in the present model. In particular, both $f(\tau)$ and $\mathcal{C}(\tau)$ should be unimodal, and vanishing not only for $\tau=0$, but also for the values of $\tau$ sufficiently close to zero – there certainly must be a minimal length of generation time.

Because the functional forms of $f(\tau)$ (37) and $\mathcal{C}(\tau)$ (39) are identical, one can treat the results of numerical simulation of a single cell lineage as referring to the whole proliferating population, if only the division rate is rescaled. Namely, for a single call lineage (scenario B) we obtain identical protein number probability distribution as for the whole growing population (scenario A), provided that in the latter case the the true division rate $\Delta_{0}=\nu_{m}\equiv\Delta_{0(A)}$ is two times smaller then in the former, i.e., $2\Delta_{0(A)}=\Delta_{0(B)}$.

Note that the simulation curves shown in Fig. 2 in the main text were generated by simulation of a single cell lineage with the cell division rate $\Delta_{0(B)}$ (scenario B, see Appendix K). However, the theoretical curves $p_{n}(x;\omega(n),b)$ shown in Fig. 2 in the main text can be interpreted in two ways, depending on how we define the $\omega$ parameter: $\omega=k/\Delta_{0(A)}$ assumes Scenario A (whole population) and $\omega=k/\Delta_{0(B)}=k/(2\Delta_{0(A)})$ assumes Scenario B (single lineage). Therefore, the simulation results shown in Fig. 2 can also be reinterpreted as the results for Scenario A (whole population) where the cell division rate was twice smaller than the value set in the simulation algorithm: $\Delta_{0(A)}=\Delta_{0(B)}/2$.

The above discussion applies to the batch culture. Analogous formulas can be derived for the continuous cell culture Powell (1956), for which, in fact, conditions for state of balanced growth can be more easily reached and maintained.

Appendix E Noise floor: Overview of other experimental results in the literature

A noise floor was observed in a number of experiments reported in the literature, however, these data were not suitable for comparison to our model. Below we overview the existing studies, which we are aware of, and which report the noise floor in gene expression.

In Ref. Taniguchi et al. (2010), the noise floor in gene expression in E. coli manifested itself as a boomerang-shaped log-log plot of the squared coefficient of variation vs. mean gene expression. However, protein levels were measured in that Reference as concentrations and not absolute molecule numbers, which excluded the effect of protein loss at cell division inherent to our model. The description of the results in Ref. Taniguchi et al. (2010) may seem slightly misleading at first glance because the plots in that reference were shown in the units of protein numbers. In fact, the protein fluorescence was measured in each cell and then its level for that cell was normalized by the volume of that same cell to get the protein concentration. The protein concentrations were then again normalized by a mean cell volume to obtain the description in the units of molecule numbers. And therefore, the resulting plots in Ref. Taniguchi et al. (2010) show the protein numbers corresponding to the content of average-sized cells, but the underlying method of measurement intentionally removed the effects of protein loss at cell division from the data. Thus, the results of Ref. Taniguchi et al. (2010) should be interpreted using a model that describes protein levels in terms of concentrations. For that reason, the Friedman’s model Friedman et al. (2006) describing protein concentrations was used for data fitting in Ref. Taniguchi et al. (2010). Our model cannot be fitted to these data because it describes protein levels in terms of protein numbers and thus it explicitly accounts for protein loss at cell division.

In Refs. Silander et al. (2012) and Wolf et al. (2015), a noise floor was observed in gene expression in E. coli. These data were unsuitable for comparison to our model because they were gated to observe the gene expression levels in only those cells that were in similar cell-cycle stages. Thus, the effect of protein loss due to cell division would not be visible in these data. A noise floor was also observed in gated data in S. cerevisiae Newman et al. (2006); Bar-Even et al. (2006).

Appendix F Distribution tails

The protein number distribution $p(x)$ as given by Eqs. (18) and (19) in the main text, i.e.,

[TABLE]

where

[TABLE]

has exponential tail of the form $\exp\{-x/[b(\omega+1)]\}$ (the leading exponent in (42)). The same is true for $p_{u}(x)$ (56) (see the next section), i.e. the solution of Eq. (31) for $\eta(q)=\Theta(q)\Theta\left(1-q\right)=1$ given by Eq. (51). In fact, it can be shown that for most of functional forms of $\eta(q)$, the solution of Eq. (31) has exactly the same exponential tail (there are, however, exceptions: for example, $\eta(q)=\frac{1}{2}\left[\delta(q)+\delta(1-q)\right]$ yields the tail of the form $\exp\{-x/[b(2\omega+1)]\}$).

This is a special case of a yet more general result. We assume here that $g(x)=0$ and the burst size PDF is exponential, as given by Eq. (10),

[TABLE]

but we allow for arbitrary $x$-dependence of $h_{d}(x)$ and $h_{p}(x)$ in Eq. (1). Following Ref. Friedlander and Brenner (2008), we neglect the $q$-integral term in the steady-state limit of Eq. (1) (for large protein number $x$, the contribution of states with still larger $x$ may be neglected) and obtain

[TABLE]

Next, we combine the derivative of (45) with respect to $x$ with the original equation and obtain

[TABLE]

where

[TABLE]

c.f. Eq. (29). From (47) it follows that at the tail of the PDF the ratio $p(x_{2})/p(x_{1})$ is well approximated by

[TABLE]

where

[TABLE]

For $h_{d}(x)=\Delta/2$ and $h_{p}(x)=1$, from (48), (49) and (50) we obtain the tail of the form $\exp\{-x/[b(\omega+1)]\}$.

Note that the the above derivation is not universally valid. For example, for the ’all or none’ mode of protein partitioning, i.e., for $\eta_{b}(q)=\left[\delta(q)+\delta(1-q)\right]/2$, one cannot simply drop out the integral over $q$ in Eq. (1). However, apart from such rather pathological situations, it seems that the derivation of the large-$x$ asymptotic behaviour of $p(x)$ proposed in Ref. Friedlander and Brenner (2008) can be generalized to the present case of stochastic protein production, provided the burst size PDF is of the form (44) and there are is no protein decay.

Appendix G Distribution of protein numbers for the uniform protein partition ratio distribution.

In this section we solve Eq. (31) (Eq. (2) of the main text) for $\gamma=0$, exponential burst size PDF as given by Eq. (44), and for uniform partition ratio distribution

[TABLE]

for $0<q<1$, where $\Theta(q)$ denotes the Heaviside theta function. $\eta(q)$ (51) is not very realistic at large $x$, where any such distribution should be peaked around $q=1/2$, like, e.g., the symmetric beta distribution (a continuous counterpart of the binomial distribution), considered in Ref. Mantzaris (2006). $\eta(q)$ (51) becomes more acceptable for $x$ of order of few protein molecules. Yet, $\eta(q)$ (51) has been used in Refs. Friedlander and Brenner (2008); Brenner and Shokef (2007), probably due to its mathematical simplicity – it is one of the few examples of a partition ratio PDF, for which analytical solutions of PBE-like equations are known.

For $\gamma=0$ and $\eta(q)$ (51), Eq. (3) reads

[TABLE]

By differentiating (52) with respect to $s$ and combining such obtained equation with the original one, integrating by parts and using the identity

[TABLE]

in order to to cancel the terms containing $G(qs)$ (note that $G^{\prime}(qs)\equiv(dG(y)/dy)_{y=qs}$), we finally get

[TABLE]

ODE (54) is equivalent to the integral equation (52). For $\hat{\nu}(s)$ (44) we have $\hat{w}(s)=\hat{\nu}(s)-1=-sb/(sb+1)$ and Eq. (54) has the following solution:

[TABLE]

where $\epsilon=\omega/(\omega+1)$, $\xi=1/[b(\omega+1)]$. The inverse Laplace transform of $G(s)$ (55) reads

[TABLE]

Because $p_{u}(x)$ (56) is a statistical mixture of two broad gamma distributions (note that $0<\epsilon<1$), $p_{u}(x)$ itself is broad ($c^{2}_{v}\geq 1/2$).

For the case of deterministic protein production, instead of $p_{u}(x)$ (56) we obtain a gamma distribution Friedlander and Brenner (2008); Brenner and Shokef (2007), cf. next Section.

Appendix H Distribution of protein numbers for the deterministic protein production and half-by-half protein partitioning

We consider here the case of unregulated, deterministic protein production, i.e., in Eq. (1) we put $k=0$ and $g(x)=\sigma>0$. For $h_{d}(x)=\Delta/2$, in the steady-state limit instead of Eq. (31) (Eq. (2)) we obtain

[TABLE]

where $\lambda=\sigma/\Delta$. Laplace transform of Eq. (57) yields

[TABLE]

and we have assumed $p(0{+})\equiv\lim_{x\to 0{+}}p(x)=0$ (validity of this assumption is checked a posteriori). For the half-by-half protein partitioning, $\eta(q)=\delta(q-\frac{1}{2})$, Eq. (58) can be rewritten as

[TABLE]

Eq. (59) can be solved by iteration; we obtain

[TABLE]

where $(a;q)_{k}$ is a $q$-Pochhammer Symbol. In (60) we have used the following identity

[TABLE]

which is a special case of the $q$-binomial theorem Andrews (1986). Note that $G(s)$ (60) can be also rewritten as

[TABLE]

where $\left[r\right]_{q}!\equiv(1-q)^{-r}(q;q)_{r}$ denotes the $q$-factorial and $e_{q}(x)$ is the $q$-exponential function.

The inverse Laplace transform of (60) reads

[TABLE]

where the coefficients $D_{i}$ are given by

[TABLE]

For $i=0$, we obtain the leading-order term, i.e., the slowest-decaying exponent, $e^{-x/\lambda}$, which determines the behaviour of $p(x;\lambda)$ (63) at the $x\to\infty$ limit. Note that $p(x;\lambda)$ (63) can be regarded as a counterpart of the probability distribution given by (42) and (43) (Eqs. (14) and (15) in the main text) in the case of deterministic protein production.

From $G(s)$ (60) we obtain both the moments ($\mu_{r}$) and the cumulants ($\kappa_{r}$) of $\mathcal{L}^{-1}[G(s)]\equiv p(x;\lambda)$ (63), namely

[TABLE]

[TABLE]

In the present case, $\mu_{r}$ can be also easily obtained from moments equations.

For a uniform distribution of protein division ratio (51), the solution of Eq. (57) reads Friedlander and Brenner (2008)

[TABLE]

The gamma distribution $p_{\gamma}(x;\lambda)$ (67) is therefore a deterministic counterpart of $p_{u}(x)$ (56).

Both $p(x;\lambda)$ (63) and $p_{\gamma}(x;\lambda)$ (67) define a one-parameter family of probability distributions; we have $p(x,\lambda)=\lambda^{-1}f(x/\lambda)$, where $f(z)=p(z,1)$ is a function of the rescaled variable $x/\lambda$. Also, both $p(x;\lambda)$ and $p_{\gamma}(x;\lambda)$ are right-skewed, unimodal, and vanish in the $x\to 0{+}$ limit, as has been assumed for $p(x;\lambda)$, cf. Fig 3. However, in contrast to the gamma distribution, for $p(x;\lambda)$ we have $p^{(k)}(0;\lambda)=0$ for arbitrary $k$, hence $p(x;\lambda)$ is not analytic at $x=0$. In consequence, the Taylor expansion of $p(x;\lambda)$ at $x_{0}=0$ does not exists. This can be shown using properties of $G(s)$ (60). We have

[TABLE]

etc., due to the obvious relation

[TABLE]

valid for any $k<\infty$.

Appendix I Similarity of gamma distributions and the distributions resulting from our model

In this section we present an alternative method of fitting of gamma distributions to the distributions resulting from our model: We match the maxima of our distribution,

[TABLE]

where $A_{n,l}=(-1)^{n-l}2^{\frac{l(l-3)}{2}}/[\prod^{l-1}_{i=1}(2^{i}-1)\prod^{n-l}_{j=1}(2^{j}-1)]$ (c.f. Eqs. (20) and (21) of the main text) and of the gamma distribution $p_{\gamma}(x;\alpha,\beta)$, c.f. Eq. (12) of the main text, i.e.:

[TABLE]

The coordinates of the gamma distribution’s maximum for $\alpha>1$ are given by:

[TABLE]

The coordinates $x_{max,p_{n}}$ and $y_{max,p_{n}}$ of the maximum of $p_{n}(x;\omega{(n)},b)$ can be found analytically for small values of $n$ (in our case, the Maple software was able to find them explicitly for $n\leq 5$). Otherwise, they can be calculated numerically. In Fig. 4, we put $b=5$ and $n=2,3,4$ (corresponding to $\omega=3,7,15$) in $p_{n}(x;\omega{(n)},b)$. By numerically solving the equations

[TABLE]

we found the $\alpha$ and $\beta$ parameters of the gamma distributions whose maxima exactly match the maxima of $p_{n}(x;\omega{(n)},b)$.

Fig. 4 shows that this way of fitting yields the gamma distribution plots that are even more similar to $p_{n}(x;\omega{(n)},b)$ than those obtained by matching the first two moments of the distributions, as in the main text, Fig 2.

Appendix J A property of the distribution $p_{n}(x)$

We denote the Laplace transform of $p_{n}(x)$ (70) as $G_{n}(s)$. From Eq. (12) of the main text, it follows that

[TABLE]

One can show that

[TABLE]

Assuming

[TABLE]

it follows from Eq. (75) that

[TABLE]

We can solve Eq. (77) iteratively for $n\geq 1$ to obtain the formulas for consecutive $p_{n+1}(x)$, starting from $p_{1}(x)=(b/2)\exp(-bx/2)$ and using (76) as the boundary condition. From Eq. (77) it follows that $p_{n+1}(x)$ has a maximum in the point where its plot intersects with the plot of $p_{n}(x)$.

Appendix K Simulation

The simulation results shown in Fig. 2 in the main text were obtained using a custom extension of the StochPy package Maarleveld et al. (2013): The standard simulation using the direct Gillespie algorithm was supplemented with cell division. Histograms were calculated along a single trajectory, which mimics the observation of a single cell lineage. The initial part of the trajectory was cut off to obtain only the steady-state behavior. The reaction kinetics (see the file NonRegulatedGeneCellDivision.psc) was given by:

[TABLE]

where \ceO denotes the gene promoter, \ceM is mRNA, and \ceP is protein. $k$ is transcription rate and $k_{p}$ is translation rate. $k_{dm}$ denotes mRNA degradation rate and its value is chosen such that mRNA life time is much shorter than the mean cell cycle duration $1/\Delta_{0{(B)}}$, consistently with the assumption of the analytical model described in the main text. Additionally, cell division occurs at the rate $\Delta_{0{(B)}}$: The custom function multiple_division() draws a random cell cycle length $T$ from a specified distribution (here: exponential with the mean $1/\Delta_{0{(B)}}$). In a given cell cycle, the simulation runs until the cell age $T$ and at that time point the molecule numbers $M$ and $P$ are divided by 2. If the remainder of the division is 1, then the remaining molecule goes to the "observed" daughter cell with probability 1/2. The simulation for the next cell cycle is initialized with the resulting molecule numbers $M$ and $P$ for the daughter cell. Simulation parameter values are shown in Table 2 and 3. The simulation code, containing the custom functions that implement cell division can be found in the file NonRegulatedGeneCellDivision.ipynb (Jupyter notebook file, Supplemental Material) 444See Supplemental Material at [URL will be inserted by publisher] for the simulation code..