1. Introduction

Transcription is one of the most fundamental cell biology processes. Gene transcription regulation is also one of the central topics in modern molecular biology. With the advent of MS2 tagging and single molecule RNA fluorescence in situ hybridisation technologies, detecting mRNA production in a single cell has come true. The measurement results reveal that mRNA copy number presented in an individual cell is low and mRNA level can vary significantly across cells in isogenic or clonal populations exposed to the same environment [5,10,12,13,15,26,32,40]. This variability may affect many important processes in cellular biology, such as response to apoptosis-inducing factors [3,45], bet-hedging in bacterial phenotypes [14,29], information processing [33], cellular decision-making [2,35,42,53], stem-cell fate decisions [8,9,10,35,36], the effectiveness of clinical treatment [30], and even cancer development [7]. Thus, understanding the origins of gene expression variabilities and establishing theoretical model to describe these variabilities have become very important in cell biology.

There is a large body of literature devoted to mean transcript level and transcription noise [5,6,12,13,14,19,28,31,32,33,46,47,48,49,54,55]. The models include Poisson expression statistics, two-state model of gene regulation [5,37,41], multi-state model of gene regulation [49], and so on. Additionally, reaction rates can also fluctuate because of stochastic variation in the global pool of housekeeping genes or because of fluctuations in environmental conditions that affect all genes [5,10,12,15]. For example, substrates, enzymes and regulatory molecules can also fluctuate and further randomize expression rates. All these lead to considerable intercellular variation in the mRNA levels in genetically identical cell populations.

Cell division is an integral part of life activity by which a single-celled fertilized egg develops into a mature organism, as well as the process by which blood cells, hair, skin, and some internal organs are renewed. Life growth and development are driven by the continuous cell divisions. It is a stochastic and complex process. The volume-based division and the time-based division are two common methods to describe cell division. For the former, when the volume of a cell is growing twice the original volume, the cell begins to divide [16,34,38]. But the time-based division is different from the volume-based cell division. A cell might divide before its volume reaching double the original volume or not divide although its volume reaches double the original volume. The time at which the cell divides is deterministic or random depending on cellular environments or cellular physiological states [17,18,38,50]. Cell division modes contain both symmetric and asymmetric ones, where asymmetric cell division contains binomial inheritance, random subtractive inheritance and random additive inheritance, etc. Symmetric cell division yields two daughter cells with equal cellular components. However, asymmetric cell division is different. Biological experiments have provided evidence for division times, but the mechanism how symmetric and asymmetric cell division, as well as cell cycle affect gene expression has not been fully elucidated. The cell cycle between adjacent divide times has global effects on mRNA synthesis. Recent experiments have suggested that it is also an important source of transcription noise [4,11,51]. At the end of cell cycle, stochastic partitioning of mRNA molecules can also create further fluctuations at cell division, it is also a vital source of transcription noise. For example, in Drosophila melanogaster, asymmetric cell division plays an important role in neural development. Elucidating the stochastic kinetics of cell division is crucial to comprehensively understanding gene expression.

Although there has been much work focused on cell divisions [4,11,16,17,18,24,34,38,43,50,51,52], there is little work of quantitative analysis exploring the effects of symmetric and asymmetric cell divisions as well as cell cycle on mean transcript level and transcription noise. In this paper, we consider the time-based division and merge cell division as well as a variety of inheritance regimes into a two-state model of gene transcription. We derive the dynamic mean transcript level and transcription noise with cell divisions. The analytical solutions agree exactly with stochastic simulations. The mathematical formalism we explored provides an effective method to get explicit moment formula for mRNA number and provides an insight into gene expression.

---

2. Model description

We construct the following theoretical model to describe the random transcription with cell division (Figure 1). The model is built on the most prevailing stochastic model for gene transcription, the so-called two-state model commonly used in the literature [15,23,26,44]. It has been postulated that the durations in the gene active (G) and inactive (G$'$) states are independent exponentially distributed, and mRNA copy numbers are controlled by a random birth and death process. Each copy of the gene switches stochastically with rates $k_{on}$ and $k_{off}$ between active and inactive states in a cell cycle. In the active state, transcription is initiated stochastically with rate $\mu_{on}$. Mature mRNA is degraded stochastically with rate $\delta$ (Figure 1. A). The model includes two additional features: the cell cycle and the partitioning of mRNA molecules effects. Since we are concerned with the effects of cell division on mean transcript level and transcription noise, we don't distinguish several phases in the cell cycle but merge them into one stage (Figure 1. B). It should be pointed out that this simplification does not influence the qualitative conclusions reached in this paper. The $i$th cell cycle is from ${\rm W}_{i-1}$ to ${\rm W}_{i}$, ${\rm T}_i = {\rm W}_{i}-{\rm W}_{i-1}$ is the $i$th cell cycle duration, where ${\rm W}_0 = 0$. ${\rm T}_1, {\rm T}_2, {\rm T}_3, \cdots$ are assumed to be the independent and identically distributed random variables and the division is a transient process, the duration for division is not counted. Cell goes from the $i$th cell cycle into the $(i+1)$th cell cycle. At the end of the cell cycle, each mRNA molecule has the same probability to be segregated to one daughter cell which is regarded as the mother cell in the next cell division (for convenience, we call it the aim daughter cell), or a number of mRNA molecules are lost or gained at each cell division, and this number is itself a random variable. We focus on exploring the analytic expressions of mean transcript level and transcription noise.

---

3. Analytical results

---

3.1. The mean and noise within a cell cycle

In this subsection, we derive the analytical expressions of the dynamic mean transcript level and the dynamic transcription noise corresponding to the reaction network (Figure 1. A) within a cell cycle. Let M$(t)$ be the copy number of the mRNA molecules at time $t$. Suppose that the initial copy number of the mRNA molecules is $m_0$ and the gene in a single cell is inactive initially. Let

$\label{eq:1} P(m,t) = {\rm Pr}\{{\rm M}(t) = m\}$

quantify the probability that there are $m$ mRNA molecules at time $t$. The mean transcript level $m(t)$ and the second moment $\mu(t)$ are defined as

$\begin{eqnarray*} m(t) = \sum\limits_{m = 0}^{\infty}mP(m,t),\ \mu(t) = \sum\limits_{m = 0}^{\infty}m^2P(m,t).\end{eqnarray*}$

The stochastic fluctuations of transcription in cell populations have been characterized by noise $\eta^2(t)$, the variance normalized by the square of $m(t)$ and noise strength, and the ratio between the variance and the mean [5,12,28,46,47]

$\label{eq:3} \eta^2(t) = \frac{\mu(t)-m^2(t)}{m^2(t)}.$

Clearly, the noise is completely determined by the mean $m(t)$ and the second moment $\mu(t)$.

For the convenience of the reader, we add Appendix A to give proofs of the next (1) and (2).

$\label{eq:4} m(t) = m_0h(t)+g(t),$

(1)

and

$\mu(t) = m_0^2h^2(t)+m_0r(t)+s(t),$

(2)

where

$\begin{split} h(t) = &e^{-\delta t},\\ g(t) = &\frac{\mu_{on}k_{on}}{(k_{on}+k_{off})(k_{on}+k_{off}-\delta)}e^{-(k_{on}+k_{off})t} -\frac{\mu_{on}k_{on}}{\delta(k_{on}+k_{off}-\delta)}e^{-\delta t}\\ &+\frac{\mu_{on}k_{on}}{\delta(k_{on}+k_{off})},\\ r(t) = &\left(\frac{2\mu_{on}k_{on}}{\delta(k_{on}+k_{off})}+1\right)\left(e^{-\delta t}-e^{-2\delta t}\right)\\ &-\frac{2 \mu_{on}k_{on}}{(k_{on}+k_{off})(\delta-k_{on}-k_{off})}\left(e^{-(k_{on}+k_{off}+\delta)t}-e^{-2\delta t}\right),\\ s(t) = &-\frac{\mu_{on}k_{on}(2\delta^2+(2\mu_{on}-k_{on}-k_{off})\delta-2\mu_{on}k_{off})}{\delta(k_{on}+k_{off})(k_{on}+k_{off}-\delta)(k_{on}+k_{off}-2\delta)} \left(e^{-(k_{on}+k_{off})t}-e^{-2\delta t}\right)\\ &-\frac{\mu_{on}k_{on}(\delta(k_{on}+k_{off})+2\mu_{on}k_{on})}{\delta^2(k_{on}+k_{off})(k_{on}+k_{off}-\delta)}\left(e^{-\delta t}-e^{-2\delta t}\right) \\ &+\frac{\mu_{on}k_{on}}{\delta(k_{on}+k_{off})}\left(1+\frac{\delta\mu_{on}+\mu_{on}k_{on}}{\delta(k_{on}+k_{off}+\delta)}\right)\left(1-e^{-2\delta t}\right)\\ &-\frac{2\mu_{on}^2k_{on}k_{off}}{\delta(k_{on}+k_{off})((k_{on}+k_{off})^2-\delta^2)}\left(e^{-(k_{on}+k_{off}+\delta) t}-e^{-2\delta t}\right). \end{split}$

(1) and (2) are the basis for the dynamic mean transcript level and the noise with cell divisions.

---

3.2. The mean and noise with cell divisions

In this section, we detail the approach by which we obtain the first two moments of ${\rm M}(t)$ with any number of cell cycles and stochastic partitioning mRNA molecules at cell divisions. For consistency reasons, we stipulate $\sum_{j = k}^{i}a_j = 0$ and $\prod_{j = k}^{i}a_j = 1$ for $i<k$ throughout this paper.

Now, we consider the stochastic model of transcription with a series of cell divisions. The $i$th cell cycle is from ${\rm W}_{i-1}$ to ${\rm W}_{i}$, ${\rm T}_i = {\rm W}_{i}-{\rm W}_{i-1}$ is the $i$th cell cycle duration, where ${\rm W}_0 = 0$. The most common distributions to model the cell cycle ${\rm T}_{i}$ are exponential distribution, Erlang distribution [1], log-normal distribution [18], uniform distribution, etc. Let $\tau_i$ be the value of the $i$th cell cycle ${\rm T}_{i}$, after which a division occurs. A variable $t_{i+1}$ measures the elapsed time since the $i$th cell division (Figure 2), then $t_{i+1}\in[0,\tau_{i+1})$. Time is accounted for throughout this hybrid model in the following manner. If the $i$th (the recent) cell division has occurred, the total elapsed time can be written as $\sum^{i}_{j = 1}\tau_j+t_{i+1}$. Hence for any $t\geq 0$ there exists an $(i+1)$-dimension vector $(\tau_1,\tau_2,\cdots,\tau_{i},t_{i+1})$ such that

$\begin{eqnarray*} t = \sum^{i}_{j = 1}\tau_j+t_{i+1},\end{eqnarray*}$

where $i$ stands for the number of the cell division before the time $t$, and $t_{i+1}$ measures the elapsed time since the $i$th (the recent) cell division.

In order to obtain mean mRNA level $m(t)$ and transcription noise $\eta^2(t)$ in each cell cycle stage, we need the following definitions and notations.

● $m_i^-$: the copy number of mRNA molecules just before the $i$th cell division;

● $m_i^+$: the copy number of mRNA molecules right after the $i$th cell division;

● $P_{\kappa}(m_{i}^+|m_{i}^-)$: the probability that there are $m_i^+$ mRNA molecules in the daughter cell after the $i$th cell division, given $m_i^-$ mRNA molecules in the mother cell before the $i$th cell division;

● $P_{j}(m,t_{i+1}|m_{i-j}^+)$: the probability that there are $m$ mRNA molecules at time $t$ after the $i$th cell division, given $m_{i-j}^+$ mRNA molecules right after the $(i-j)$th cell division;

● $\Phi_j = P_{\kappa}\left(m_{j}^+|m_{j}^-\right)$$P_{0}\left(m_{j}^-,\tau_j|m_{j-1}^+\right)$;

● $\theta_j({m_i^-}) = \sum_{{m_i^+} = 0}^{{m_i^-}} ({m_i^+})^j$$P_{\kappa}({m_i^+}|{m_i^-})$;

● $\vartheta_j({m_i^-}) = \sum_{{m_i^+} = {m_i^-}}^{\infty} ({m_i^+})^jP_{\kappa}({m_i^+}|{m_i^-})$;

● $\sum_{i,j}^{'} = \sum_{{m_i^-} = 0}^{\infty}\sum_{{m_i^+} = 0}^{{m_i^-}}\sum_{m_{i-1}^- = 0}^{\infty}$$\sum_{m_{i-1}^+ = 0}^{m_{i-1}^-}\cdots\sum_{m_{j}^- = 0}^{\infty}\sum_{m_{j}^+ = 0}^{m_{j}^-}$;

● $\sum_{i,j}^{''} = \sum_{{m_i^-} = 0}^{\infty}\sum_{{m_i^+} = {m_i^-}}^{\infty}\sum_{m_{i-1}^- = 0}^{\infty}$$\sum_{m_{i-1}^+ = {m_{i-1}^-}}^{\infty}\cdots\sum_{m_{j}^- = 0}^{\infty}\sum_{m_{j}^+ = {m_{j}^-}}^{\infty}$;

● $A_{i} = \sum_{i,1}^{'}m_0(\tau_{i+1}|{m_i^+})\prod^{i}_{j = 1}\Phi_j$, $A_{0} = m_0(\tau_{1})$;

● $C_{i} = \sum_{i,1}^{''}m_0(\tau_{i+1}|{m_i^+})\prod^{i}_{j = 1}\Phi_j$, $C_{0} = m_0(\tau_{1})$;

● $B_{i} = \sum_{i,1}^{'}\mu_0(\tau_{i+1}|{m_i^+})\prod^{i}_{j = 1}\Phi_j$, $B_{0} = \mu_0(\tau_{1})$;

● $D_{i} = \sum_{i,1}^{''}\mu_0(\tau_{i+1}|{m_i^+})\prod^{i}_{j = 1}\Phi_j$, $D_{0} = \mu_0(\tau_{1})$.

Define

$\label{eq:7} m_j(t_{i+1}|m_{i-j}^+) = \sum\limits_{m = 0}^{\infty}mP_{j}(m,t_{i+1}|m_{i-j}^+),\ i = 0,1,2,\cdots,\ j = 0,1,2,\cdots,i$

and

$\label{eq:8} \mu_j(t_{i+1}|m_{i-j}^+) = \sum\limits_{m = 0}^{\infty}m^2P_{j}(m,t_{i+1}|m_{i-j}^+),\ i = 0,1,2,\cdots,\ j = 0,1,2,\cdots,i.$

The aim in the following is to obtain the expressions of $m_i(t_{i+1}|m_0)$ and $\mu_i(t_{i+1}|m_0)$. For symmetric division and binomial inheritance, based on the total probability formula, $P_{i}(m,t_{i+1}|m_0)$ can be decomposed into the following form [25]

$P_{i}(m,t_{i+1}|m_0) = \sum\limits_{m_{i}^- = 0}^{\infty}\sum\limits_{{m_i^+} = 0}^{{m_i^-}}P_{i-1}({m_i^-},\tau_i|m_0^+)P_{\kappa}({m_i^+}|{m_i^-}) P_{0}(m,t_{i+1}|{m_i^+}).$

(3)

By induction, one can get

$\begin{eqnarray*} P_{i}(m,t_{i+1}|m_0) = \sum\limits_{i,1}^{'}\left(P_{0}(m,t_{i+1}|{m_i^+})\prod^{i}_{j = 1}\Phi_j\right).\end{eqnarray*}$

(4)

From (1) and (2), we have by noting $m_{0}^+ = m_{0}^- = m_0$

$\begin{eqnarray*} m_0(t_1|m_{0}) = m(t_1) = m_0h(t_1)+g(t_1)\end{eqnarray*}$

(5)

and

$\mu_0(t_1|{m_0}) = \mu(t_1) = m_0^2h^2(t_1)+m_0r(t_1)+s(t_1),$

(6)

respectively. Taking advantage of (4)-(6), we obtain

$m_i(t_{i+1}|m_0) = \sum\limits_{i,1}^{'}m_0(t_{i+1}|{m_i^+})\prod\limits_{j = 1}^{i} {\Phi_j} \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ = h(t_{i+1})\sum\limits_{i-1,1}^{'}\sum\limits_{{m_i^-} = 0}^{\infty}\theta_1({m_i^-}) P_{0}({m_i^-},\tau_i|m_{i-1}^+)\prod\limits_{j = 1}^{i-1} {\Phi_j} \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ +g(t_{i+1})\sum\limits_{i-1,1}^{'}\sum\limits_{{m_i^-} = 0}^{\infty}\theta_0({m_i^-}) P_{0}({m_i^-},\tau_i|m_{i-1}^+)\prod\limits_{j = 1}^{i-1} {\Phi_j}$

(7)

and

$\ \ \ \ \ \ \ \ \ \ \ \ \mu_i(t_{i+1}|m_0) = \sum\limits_{i,1}^{'}\mu_0(t_{i+1}|{m_i^+})\prod\limits_{j = 1}^{i} {\Phi_j} \\ = h^2(t_{i+1})\sum\limits_{i-1,1}^{'}\sum\limits_{{m_i^-} = 0}^{\infty} \theta_2({m_i^-})P_{0}({m_i^-},\tau_i|{m_{i-1}^+})\prod\limits_{j = 1}^{i-1} {\Phi_j}\\ \ \ \ \ +r(t_{i+1})\sum\limits_{i-1,1}^{'}\sum\limits_{{m_i^-} = 0}^{\infty} \theta_1({m_i^-})P_{0}({m_i^-},\tau_i|{m_{i-1}^+})\prod\limits_{j = 1}^{i-1} {\Phi_j}\\ \ \ \ \ +s(t_{i+1})\sum\limits_{i-1,1}^{'}\sum\limits_{{m_i^-} = 0}^{\infty} \theta_0({m_i^-})P_{0}({m_i^-},\tau_i|{m_{i-1}^+})\prod\limits_{j = 1}^{i-1} {\Phi_j}.$

(8)

---

3.3. The mean and noise with some specific divisions

In this subsection, we will derive the specific expressions of (7) and (8) for some specific cell divisions. Usually, there are two distinct types of cell divisions, i.e., symmetric division and asymmetric division, where asymmetric cell division includes binomial inheritance, random subtractive inheritance and random additive inheritance, etc.

---

3.3.1. Symmetric cell division

For symmetric cell division, the mRNA molecules in mother cell are equally divided into two daughter cells [20,21,22,27]. Note that ${m_i^+}$ is a nature number, $P_{\kappa}({m_i^+}|{m_i^-})$ can be expressed as

$\label{eq:15} P_{\kappa}({m_i^+}|{m_i^-}) = 1_{\{m_i^-/2\}}({m_i^+})+1_{\{({m_i^-}\pm1)/2\}}({m_i^+}),$

where $1_{\Omega}$ is the indicator function of $\Omega$.

For convenience, let ${m_i^+}$ take real numbers, we obtain

$\label{eq:16} m_i(t_{i+1}|m_0) = \frac{1}{2}h(t_{i+1})A_{i-1}+g(t_{i+1})$

and

$\label{eq:17} \mu_i(t_{i+1}|m_0) = \frac{h^2(t_{i+1})}{4}B_{i-1} +\frac{r(t_{i+1})}{2}A_{i-1}+s(t_{i+1}),$

where

$\begin{eqnarray*} A_{i-1} = 2m_0\prod ^{i}_{j = 1}\left( \frac{h(\tau_{j})}{2}\right) +\sum^{i}_{k = 1}\left(g(\tau_{k}) \prod ^{i}_{j = k+1}\left( \frac{h(\tau_{j})}{2}\right) \right)\end{eqnarray*}$

and

$\begin{eqnarray*} B_{i-1} = 4m_0^2\prod ^{i}_{j = 1} \left(\frac{h(\tau_{j})}{2}\right)^{2}+\sum^{i}_{k = 1}\left( \prod ^{i}_{j = k+1} \left(\frac{h(\tau_{j})}{2}\right)^{2} \left( \frac{r(\tau_{k})}{2}A_{k-2}+s(\tau_{k})\right)\right).\end{eqnarray*}$

---

3.3.2. Binomial inheritance

For binomial inheritance, every mRNA molecule in mother cell can be assigned into the daughter cell with the equal probability $p_i$. The number of mRNA molecules assigned into the daughter cell obeys binomial distribution. Hence

$\label{eq:20} P_{\kappa}({m_i^+}|{m_i^-}) = \binom{{m_i^-}}{{m_i^+}}p_i^{{m_i^+}}(1-p_i)^{{m_i^-}-{m_i^+}},\ i = 0,1,2,\cdots,$

where $p_i\in(0,1)$. In this case, we obtain

$\label{eq:21} m_i(t_{i+1}|m_0) = p_ih(t_{i+1})A_{i-1}+g(t_{i+1})$

and

$\label{eq:22} \begin{split} \mu_i(t_{i+1}|m_0) = &p_i^2h^2(t_{i+1})B_{i-1}+\left( (p_i-p_i^2)h^2(t_{i+1})+p_ir(t_{i+1})\right) A_{i-1}\\ &+s(t_{i+1}), \end{split}$

where

$\begin{eqnarray*} A_{i-1} = \frac{m_0}{p_i}\prod^{i}_{j = 1}\left(p_jh(\tau_{j})\right)+\sum\limits_{k = 1}^{i}\left(\frac{p_kg(\tau_{k})}{p_i}\prod^{i}_{j = k+1}p_jh(\tau_{j})\right)\end{eqnarray*}$

and

$\begin{equation*} \begin{split} B_{i-1} = &\frac{m_0^2}{p_i^2}\prod^{i}_{j = 1}\left(p_jh(\tau_{j})\right)^{2}\\ &+\sum\limits_{k = 1}^{i}\left(\left((p_{k-1}-p_{k-1}^2)h^2(\tau_{k})+p_{k-1}r(\tau_{k})\right)A_{k-2}\prod^{i}_{j = k+1}\left(p_jh(\tau_{j})\right)^{2}\right)\\ &+\sum\limits_{k = 1}^{i}\left(s(\tau_{k})\prod^{i}_{j = k+1}\left(p_jh(\tau_{j})\right)^{2}\right),\\\end{split} \end{equation*}$

where $p_0\equiv0.$

---

3.3.3. Random subtractive inheritance

We now consider the case where a number of mRNA molecules are lost at each cell division. This number is a random variable and obeys binomial distribution with population size $2\eta_i$ and probability $0.5$, then the average loss number is $\eta_i$. Hence

$\label{eq:25} P_{\kappa}({m_i^+}|{m_i^-}) = \binom{2\eta_i}{{m_i^-}-{m_i^+}}2^{-2\eta_i}.$

In this case, we obtain

$\label{eq:26} m_i(t_{i+1}|m_0) = h(t_{i+1})A_{i-1}+g(t_{i+1})-\eta_ih(t_{i+1})$

and

$\label{eq:27} \begin{split} \mu_i(t_{i+1}|m_0) = &h^2(t_{i+1})B_{i-1} +\left(r(t_{i+1})-2\eta_ih^2(t_{i+1})\right)A_{i-1}\\ &+h^2(t_{i+1})\left( \eta_i^2+\frac{\eta_i}{2}\right)-\eta_ir(t_{i+1})+s(t_{i+1}), \end{split}$

where

$\begin{eqnarray*} A_{i-1} = m_0 \prod^{i}_{j = 1}h(\tau_{j}) +\sum\limits_{k = 1}^{i}\left(\left( g(\tau_{k})-\eta_{k-1}h(\tau_{k})\right) \prod^{i}_{j = k+1}h(\tau_{j})\right)\end{eqnarray*}$

and

$\begin{eqnarray*} B_{i-1} = &m_0^2\prod ^{i}_{j = 1} h^{2}(\tau_{j})+\sum\limits_{k = 1}^{i}\left( \left(r(\tau_{k})-2\eta_{k-1}h^2(\tau_{k})\right)A_{k-2}\prod ^{i}_{j = k+1} \left(h(\tau_{j})\right)^{2}\right)\\ &+\sum\limits_{k = 1}^{i}\left(\left(h^2(\tau_{k})\left( \eta_{k-1}^2+\frac{\eta_{k-1}}{2}\right)\right)\prod ^{i}_{j = k+1} \left(h(\tau_{j})\right)^{2}\right)\\ &+\sum\limits_{k = 1}^{i}\left(\left(-\eta_{k-1}r(\tau_{k})+s(\tau_{k})\right)\prod ^{i}_{j = k+1} \left(h(\tau_{j})\right)^{2}\right). \end{eqnarray*}$

---

3.3.4. Random additive inheritance

We now consider the case where a number of mRNA molecules are gained at the $i$th cell division. This number is a random variable and obeys binomial distribution with population size $2\zeta _i$ and probability $0.5$, then the average loss number is $\zeta _i$. Hence

$\label{eq:30} P_{\kappa}({m_i^+}|{m_i^-}) = \binom{2\zeta _i}{{m_i^+}-{m_i^-}}2^{-2\zeta _i}.$

We firstly decompose $P_{i}(m,t_{i+1}|m_0)$ into the following form

$\begin{eqnarray*} P_{i}(m,t_{i+1}|m_0) = \sum\limits_{m_{i}^- = 0}^{\infty}\sum\limits_{{m_i^+} = {m_i^-}}^{\infty}P_{i-1}({m_i^-},\tau_i|m_0^+)P_{\kappa}({m_i^+}|{m_i^-}) P_{0}(m,t_{i+1}|{m_i^+}).\end{eqnarray*}$

It follows by induction that

$\begin{eqnarray*} P_{i}(m,t_{i+1}|m_0) = \sum\limits_{i,1}^{''}\left(P_{0}(m,t_{i+1}|{m_i^+})\prod^{i}_{j = 1}\Phi_j\right).\end{eqnarray*}$

Similarly to the random subtractive inheritance, we obtain

$\label{eq:33} m_i(t_{i+1}|m_0) = h(\tau_{i+1})C_{i-1}+g(\tau_{i+1})+\zeta _ih(\tau_{i+1})$

and

$\label{eq:34} \begin{split} \mu_i(t_{i+1}|m_0) = &h^2(t_{i+1})D_{i-1} +\left(r(t_{i+1})+2\zeta _ih^2(t_{i+1})\right)C_{i-1}\\ &+h^2(t_{i+1})\left( \zeta _i^2+\frac{\zeta _i}{2}\right)+\zeta _ir(t_{i+1})+s(t_{i+1}), \end{split}$

where

$\begin{eqnarray*} C_{i-1} = m_0 \prod^{i}_{j = 1}h(\tau_{j}) +\sum\limits_{k = 1}^{i}\left(\left( g(\tau_{k})+\zeta _{k-1}h(\tau_{k})\right) \prod^{i}_{j = k+1}h(\tau_{j})\right)\end{eqnarray*}$

and

$\begin{eqnarray*} D_{i-1} = &m_0^2\prod ^{i}_{j = 1} h^{2}(\tau_{j})+\sum\limits_{k = 1}^{i}\left( \left(r(\tau_{k})+2\zeta_{k-1}h^2(\tau_{k})\right)C_{k-2}\prod ^{i}_{j = k+1} \left(h(\tau_{j})\right)^{2}\right)\\ &+\sum\limits_{k = 1}^{i}\left(h^2(\tau_{k})\left( \zeta_{k-1}^2+\frac{\zeta_{k-1}}{2}\right)\prod ^{i}_{j = k+1} \left(h(\tau_{j})\right)^{2}\right)\\ &+\sum\limits_{k = 1}^{i}\left(\left(\zeta_{k-1}r(\tau_{k})+s(\tau_{k})\right)\prod ^{i}_{j = k+1} \left(h(\tau_{j})\right)^{2}\right). \end{eqnarray*}$

---

4. Numerical results

In this section, we will present our numerical results. For a given cell cycle distribution (exponential distribution, Erlang distribution [1], log-normal distribution [18], uniform distribution, etc.), we first randomly generate a series of division cycles $\tau_1,\tau_2,\cdots,\tau_{i}$ whose mean is denoted by $\tau$ (representing the mean cell-cycle length). $m_i(t_{i+1}|m_0)$ and $\mu_i(t_{i+1}|m_0)$ denote respectively the mean and noise of mRNA numbers with the generating fixed cell division cycles $\tau_1,\tau_2,\cdots,\tau_{i}$. By averaging the different $m_i(t_{i+1}|m_0)$ and $\mu_i(t_{i+1}|m_0)$ corresponding to the different series of division cycles $\tau_1,\tau_2,\cdots,\tau_{i}$, we obtain the time-evolutional analytic mean and noise of mRNA number, respectively. Similarly, with a series of division cycles $\tau_1,\tau_2,\cdots,\tau_{i}$, by using the Gillespie stochastic simulation algorithm, we generate 10000 run at every time point and calculate the stochastic trajectory. By averaging these trajectory data, we obtain the time-evolutional numerical mean mRNA number.

For the narrative convenience, we introduce the following abbreviations:

● S: symmetric cell division;

● BF: binomial inheritance with fixed probability $p = 0.5$;

● BR: binomial inheritance with random probability $p\in[0.4, 0.6]$;

● RA: random additive inheritance;

● AS: random additive and subtractive inheritance;

● RS: random subtractive inheritance.

For S, $m_i^+ = \frac{m_i^-}{2}$ if $m_i^-$ is even, otherwise $m_i^+$ takes randomly between $\frac{m_i^--1}{2}$ and $\frac{m_i^-+1}{2}$. For BF, $m_i^+$ obeys the binomial distribution with parameters $m_i^-$ and $p = 0.5$. For BR, each mRNA molecule is assigned to the daughter cell with the probability $p$ at the end of a given cell cycle, but $p$ may be different at the ends of the different cell cycles, which can be randomly selected from the interval [0.4, 0.6]. For AS, $m_i^+$ takes randomly between $m_i^-+n$ and $m_i^--n$ where $n$ obeys the binomial distribution with parameters $m_i^-$ and $p = 0.5$. For constant cycle, we take the average in a few cycles as the steady-state mRNA mean or steady-state mRNA noise, otherwise we average the values in the time interval [200,800] as the steady-state mRNA mean or steady-state mRNA noise. The parameters $k_{on} = k_{off} = 1$, $\mu_{on} = 10$, $\delta = 0.0529$ [52]. The initial mRNA molecule number $m_0 = 4.$ By (1) and (2) we know the steady-state mRNA mean and the steady-state mRNA noise with no cell division are $\overline{m} = 94.5180$ and $\overline{\eta} = 0.0363,$ respectively.

The analytic solutions agree exactly with Gillespie stochastic simulation using plausible parameters $k_{on} = k_{off} = 1$, $\mu_{on} = 10$, $\delta = 0.0529$ [52]. Figure 4 shows the analytic solutions and the numerical solutions of mean transcript level with BR and AS, where the cell cycle obeys log-normal distribution with the mean $\tau = 120$.

Table 1 shows the steady-state mean transcript level (the steady-state mRNA noise) with S, BF, BR and different cell cycle distributions, where the mean of cell cycle is $\tau = 120$.

Table 2 shows the steady-state mean transcript level (the steady-state mRNA noise) with BA, BS, AS and different cell cycle distributions, where the mean of cell cycle is $\tau = 120$.

From Table 1 and Table 2 we can see that the steady-state mRNA noise with S is less than that with BF, and the steady-state mean transcript levels with the two modes are same. The steady-state mean transcript levels (the steady-state mRNA noises) with BF and BR are almost same. The steady-state mean transcript levels with S, BF and BR fall by about 7.65%-8.5% compared with no cell division. The steady-state mean transcript level with RA rises by about 1.5%-1.9%. The steady-state mean transcript levels with AS and RS fall by 0.4%-0.43% and 2.18%-2.42% respectively, compared with no cell division. All division modes except RA raise the steady-state mRNA noises.

In the case of the same cell cycle distribution, Figure 5 show that the steady-state mRNA noise with S is less than that with BF, but the steady-state mean transcript level with S is always equal to that with BF; The steady-state mean transcript levels (the steady-state mRNA noises) with BF and BR are almost same; The steady-state mean transcript levels (the steady-state mRNA noises) with S, BF and BR increase (decrease) with the average cell cycle duration, while the steady-state mean transcript level (the steady-state mRNA noise) with RA decreases (increases) with the average cell cycle duration; The steady-state mean transcript level is almost constant although the steady-state mRNA noise with AS decreases with the average cell cycle duration; The steady-state mean transcript level with RS increases with the average cell cycle duration although the steady-state mRNA noise with RS is almost constant. By fitting we obtain the relationship between the steady-state mean transcript level (the steady-state mRNA noise) and the cell cycle duration is very close to the power function with the form $f(x) = ax^b+c$. From the fitting we can predict the steady-state mean transcript level (the steady-state mRNA noise) in different cycles.

---

5. Conclusion and discussion

Gene transcription is a complex stochastic process, which result in high variability in gene expression. The sources of transcription noise come about in two ways, intrinsic source and extrinsic source. Intrinsic noise results from the inherent stochasticity of biochemical processes. The extrinsic noise originates from the cellular environmental perturbations or from extracellular signals regulating intracellular processes. Although there is a large body of literature devoted to mean transcript level and transcription noise [5,6,12,13,14,19,28,31,32,33,46,47,48,49,54,55] even cell divisions [4,11,16,17,18,24,34,38,43,50,51,52], there is little work of quantitative analysis exploring the effects of symmetric and asymmetric cell divisions as well as cell cycle on mean transcript level and transcription noise. Here, we have constructed a theoretical model to describe the random transcription with cell division. We merge several phases in a cell cycle into one stage and derive the analytical formulas for mean transcript level and mRNA noise corresponding to our model. By analysis, we have shown that symmetric cell division and binomial inheritance not only decrease mean transcript level but also increase transcription noise. For symmetric cell division and binomial inheritance, the steady-state mean mRNA level increases with the average cell cycle duration. All results are confirmed by Gillespie stochastic simulation using plausible parameters. The mathematical formalism we explored provides a method for yielding explicit formula of each moment of mRNA. Our method can be extended to more complex cases of gene expression. Our next step is to take gene replication into account and investigate the effects of specific phase of cell cycle on mean transcript level and transcription noise.

---

Appendix A. Proofs of (1) and (2)

Let ${\rm U}(t)$ be the number of active gene in a single cell at time $t$. Then ${\rm U}(t)$ is a random process. ${\rm U}(t) = 0$ and ${\rm U}(t) = 1$ denote the gene is inactive and active, respectively. By the assumption, we know

$\label{eq:37} {\rm M}(0) = m_0,\ {\rm U}(0) = 0.$

Let

$\label{eq:38} P^{(j)}(m,t) = {\rm Pr}\{ {\rm M}(t) = m, {\rm U}(t) = j\}, \ m = 0,1,2,\cdots,\ j = 0,1$

denote the probability that there are $m$ mRNA molecules and $j$ active genes. The initial conditions are

$P^{(j)}(m,0) = \left\{ \begin{aligned} 1, &\ m = m_0,\ j = 0, \\ 0, &\ {\rm othere wise}. \end{aligned} \right.$

(9)

Then

$\label{eq:40} P(m,t) = {\rm Pr}\{{\rm M}(t) = m\} = P^{(0)}(m,t)+P^{(1)}(m,t).$

We can derive the following chemical master equations

$\left\{ \begin{aligned} \frac{\partial }{\partial t}\left(P^{(0)}(m,t)\right) = &-(k_{on}+\delta m)P^{(0)}(m,t)+k_{off}P^{(1)}(m,t)\\ &+\delta(m+1)P^{(0)}(m+1,t), \\ \frac{\partial}{\partial t}\left(P^{(1)}(m,t)\right) = &k_{on}P^{(0)}(m,t)-(k_{off}+\mu_{on}+m\delta )P^{(1)}(m,t)\\ &+\mu_{on}P^{(1)}(m-1,t)+\delta(m+1)P^{(1)}(m+1,t). \end{aligned} \right.$

(10)

To obtain the mean transcript level $m(t)$ and second moment $\mu(t)$, we define

$\begin{eqnarray*} m^{(0)}(t) = \sum\limits_{m = 0}^{\infty}mP^{(0)}(m,t),\ m^{(1)}(t) = \sum\limits_{m = 0}^{\infty}mP^{(1)}(m,t)\end{eqnarray*}$

and

$\begin{eqnarray*} \mu^{(0)}(t) = \sum\limits_{m = 0}^{\infty}m^2P^{(0)}(m,t),\ \mu^{(1)}(t) = \sum\limits_{m = 0}^{\infty}m^2P^{(1)}(m,t).\end{eqnarray*}$

Then we have

$\label{eq:44} m(t) = m^{(0)}(t)+m^{(1)}(t),\ \mu(t) = \mu^{(0)}(t)+\mu^{(1)}(t),$

which together with (9) yields

$m(0) = m_0,\ \mu(0) = m_0^2.$

(11)

Set

$\begin{eqnarray*} P^{(1)}(t) = \sum\limits_{m = 0}^{\infty}P^{(1)}(m, t).\end{eqnarray*}$

Then by [39,52], we have

$\label{eq:47} P^{(1)}(t) = \frac{k_{on}}{k_{on}+k_{off}}-\frac{k_{on}}{k_{on}+k_{off}}e^{-(k_{on}+k_{off})t}.$

(12)

Multiplying both sides of each equation in (10) by $m$, and summing $m$ from 0 to $\infty$, respectively, we obtain

$\left\{ \begin{aligned} {m^{(0)}}'(t)& = -(k_{on}+\delta)m^{(0)}(t)+k_{off}m^{(1)}(t),\\ {m^{(1)}}'(t)& = k_{on}m^{(0)}(t)-(k_{off}+\delta)m^{(1)}(t)+\mu_{on}P^{(1)}(t). \end{aligned} \right.$

(13)

It follows by summing its both sides that

$m'(t)+\delta m(t) = \mu_{on}P^{(1)}(t).$

(14)

Combining (14) with the initial condition (11) we get (1).

By eliminating $m^{(0)}(t)$ from (13) we can get the differential equation

$\label{eq:50} {m^{(1)}}'(t)+(k_{on}+k_{off}+\delta)m^{(1)}(t) = k_{on}m(t)+\mu_{on}P^{(1)}(t).$

(15)

Combining with (9) we obtain the initial condition of $m^{(1)}(t)$

$\label{eq:51} m^{(1)}(0) = 0.$

(16)

Solving equation (15) with the initial condition (16) we get

$\label{eq:52} \begin{split} m^{(1)}(t) = &\frac{m_0k_{on}}{k_{on}+k_{off}}\left( e^{-\delta t}-e^{-(k_{on}+k_{off}+\delta) t}\right) \\ &+\frac{\mu_{on}k_{on}(\delta-k_{off})}{\delta(k_{on}+k_{off})(k_{on}+k_{off}-\delta)}e^{-(k_{on}+k_{off})t}\\ &-\frac{k^2_{on}\mu_{on}}{\delta(k_{on}+k_{off})(k_{on}+k_{off}-\delta)}e^{-\delta t} \\ &-\frac{\mu_{on}k_{on}k_{off}}{\delta(k_{on}+k_{off})(k_{on}+k_{off}+\delta)}e^{-(k_{on}+k_{off}+\delta) t}\\ &+\frac{\mu_{on}k_{on}(k_{on}+\delta)}{\delta(k_{on}+k_{off})(k_{on}+k_{off}+\delta)}.\\ \end{split}$

Now multiplying both sides of each equation in (10) by $m^2$, and summing $m$ from 0 to $\infty$, respectively, we get the differential equations

$\left\{ \begin{aligned} {\mu^{(0)}}'(t) = &-(k_{on}+2\delta)\mu^{(0)}(t)+k_{off}\mu^{(1)}(t)+\delta m^{(0)}(t),\\ {\mu^{(1)}}'(t) = &k_{on}\mu^{(0)}(t)-(k_{off}+2\delta)\mu^{(1)}(t)+(2\mu_{on}+\delta)m^{(1)}(t)\\ &+\mu_{on}P^{(1)}(t). \end{aligned} \right.$

(17)

Adding the two equations in (17) yields

$\mu'(t)+2\delta \mu(t) = \delta m(t)+2\mu_{on} m^{(1)}(t)+\mu_{on}P^{(1)}(t).$

(18)

Combining (18) with (1), (12) and the condition (11), we can get (2).

---

Acknowledgments

We would like to thank the two anonymous reviewers for their valuable and precious comments and suggestions.

---