Quantifying gene expression variability arising from randomness in cell division times

Abstract

The level of a given mRNA or protein exhibits significant variations from cell-to-cell across a homogeneous population of living cells. Much work has focused on understanding the different sources of noise in the gene-expression process that drive this stochastic variability in gene-expression. Recent experiments tracking growth and division of individual cells reveal that cell division times have considerable inter-cellular heterogeneity. Here we investigate how randomness in the cell division times can create variability in population counts. We consider a model by which mRNA/protein levels in a given cell evolve according to a linear differential equation and cell divisions occur at times spaced by independent and identically distributed random intervals. Whenever the cell divides the levels of mRNA and protein are halved. For this model, we provide a method for computing any statistical moment (mean, variance, skewness, etcetera) of the mRNA and protein levels. The key to our approach is to establish that the time evolution of the mRNA and protein statistical moments is described by an upper triangular system of Volterra equations. Computation of the statistical moments for physiologically relevant parameter values shows that randomness in the cell division process can be a major factor in driving difference in protein levels across a population of cells.

Appendix

Proof

(of Proposition 2) Let $z_k:=\begin{bmatrix} m(t_k^-)&p(t_k^-) \end{bmatrix}^\intercal ,\,\forall k\in \mathbb {N}$, and note that

$$\begin{aligned} z_{k+1}=\mathcal {M}(z_k, \tau _k), \forall k\in \mathbb {N}, \end{aligned}$$

where

$$\begin{aligned} \mathcal {M}\left( \begin{bmatrix} a \\ b \end{bmatrix},h\right) =\begin{bmatrix} e^{-\gamma _m h}\frac{a}{2} + \frac{k_m(1-e^{-\gamma _m h})}{\gamma _m} \\ e^{-\gamma _p h}\frac{b}{2} + \kappa _p\int _0^h e^{-\gamma _p (h-s)}\left( e^{-\gamma _m s}\frac{a}{2} + \frac{k_m(1-e^{-\gamma _m s})}{\gamma _m}\right) ds \end{bmatrix} \end{aligned}$$

(41)

for positive $a$ and $b$ and $h\in (\underline{\tau },\bar{\tau })$. We can assume without loss of generality that $z_1 \in \mathcal {I}_1$ where $\mathcal {I}_1 =[0,\frac{k_m}{\gamma _m}]\times [0,\frac{k_p}{\gamma _p}]$. In fact, one can conclude from (1), (2) that if the initial mRNA and protein counts do not belong to $\mathcal {I}_1$ then $(m(t_k), p(t_k)) \in \mathcal {I}_1$ after a sufficiently large time $k\in \mathbb {N}$ with probability one. Moreover if $(m(0), p(0)) \in \mathcal {I}_1$ then $(m(t), p(t)) \in \mathcal {I}_1$ for every time $t \in \mathbb {R}_{\ge 0}$. In particular $z_1\in \mathcal {I}_1$. Now, if $z_1\in \mathcal {I}_1$ then $z_k \in \mathcal {I}_k$ where the sets $\mathcal {I}_k,\,k\in \mathbb {N}$ are defined recursively

$$\begin{aligned} \mathcal {I}_{k+1} = \{\mathcal {M}(y,h)| y\in \mathcal {I}_{k}, h\in (\underline{\tau },\bar{\tau }) \}, \ \forall k\in \mathbb {N}. \end{aligned}$$

From (41) we conclude that $\mathcal {I}_{k+1} \subseteq \mathcal {I}_{k},\,\forall k\in \mathbb {N}$, and we can also conclude that $\cap _{k=0}^\infty \mathcal {I}_{k}$ has non-empty interior. Choose then an open set $A$ in $\cap _{k=0}^\infty \mathcal {I}_{k}$. By construction the Markov chain $z_k$ is irreducible [see Meyn and Tweedie (2009), Ch. 4] with respect to the indicator function of set $A$, which is equivalent to saying that one can reach any open set in $A$ for any initial condition $z_1 \in \mathcal {I}_1$. One can also conclude that for any initial condition $z_1 \in \mathcal {I}_1$, the state $z_k,\,k\in \mathbb {N}$, visits the set $A$ in infinite number of times and that the Markov chain is aperiodic meaning that it does not take values in disjoint sets visited periodically. This implies that the chain is Harris recurrent [see Meyn and Tweedie (2009), Ch. 9] and admits a unique invariant measure $\pi _{\mathrm{MC}}$ (which can be made a probability measure by properly scaling since the chain takes values in a bounded set) and, from the aperiodic ergodic theorem [cf. Meyn and Tweedie (2009), p. 309], we have

$$\begin{aligned} \lim _{n\rightarrow \infty }\mathsf {P}^n(z_1,B) = \pi _{\text {MC}}(B), \end{aligned}$$

for any initial condition $z_1 \in \mathcal {I}_1$ and for any open set $B \subseteq \mathcal {I}_1$, where $\mathsf {P}^n(y,B) := \text {Prob}(z_{n+1}\in B|z_1=y)$ for $n\in \mathbb {N}$. Using the results in (Antunes et al. 2013), one can prove that the process $w(t): = (m(t),p(t),\xi (t))$ where $\xi (t):=t-t_{N(t)},\,N(t):= \max \{k \in \mathbb {N}_{0}:t_k\le t\}$ is the time since the last division, is a Piecewise deterministic process in the sense of (Davis 1993). Then $z_k$ is the so called imbedded Markov chain of this process and using the connection between stationary distribution of the embedded chain and of the piecewise deterministic process (Davis 1993) one can conclude that the piecewise deterministic process also has an invariant distribution $\pi _{PD}$. Then for any function $g(m,p) = m^{n_1}p^{n_2},\,n_1,n_2 \in \mathbb {N}$, we have

$$\begin{aligned} \lim _{t\rightarrow \infty } \mathbb {E}[g(m(t),p(t))] = \int g(m,p) \pi _{PD}(dw), \ \ w=(m,p,\xi ), \end{aligned}$$

which concludes the proof. $\square$