1. Introduction

According to genetic central dogma, gene expression includes two main processes of transcription of genetic information to mRNA, and translation of each mRNA to protein. It is basically a biochemical process, which involves recruitment of transcription factors and polymerases, transition between active state and inactive state of promoter and chromatin remodeling, etc [2,47,3,51,31,71,49,30,59]. The concentration of some specific factors, such as RNA polymerases, eRNA and transcriptional factors, is significantly different even in the isogenic cell population [33,28,48,61]. Variations in these specific factors and the inherent randomness of biochemical reactions can result in stochasticity in the number of protein level across identical cells [49,33,62]. This implies that the timing of a cellular event that protein level triggers at a critical value is stochastic in nature. So, many biologists have been paying close attention to the study of the critical thresholds of some special proteins and showed that phenotypic diversity and cell fate decision often depend on the number of the particular protein [1,10,15]. Phenotype of B. subtilis switching to competence depends on the number of ComK molecules up to a certain amount [49,57,35]. Some researchers showed that lysis time for bacteriophage $\lambda$ relies on the accumulation of holin protein in the cell membrane exceeding a threshold [10,17,53]. Naturally, stochastic switching of phenotypes or cellular functions is linked closely with stochasticity in the minimal time, at which the number of some certain proteins hits a critical threshold. In the context of gene expression, such minimal time is often called as the first passage time (FPT) [53,52,21]. FPT provides a frame-work for studying the time for cellular functions conversion caused by crossing a threshold.

So far, the study of dynamical behaviour of FPT has been receiving increasing attention [53,52]. In [53], authors revealed mechanisms that transcription and translation efficiencies independently modulate the mean and variation of FPT. The impacts of different models of transcriptional and translational bursts on the mean of FPT were discussed in [52]. However, these studies did not investigate the influence of external stimuli on FPT. In fact, cells are always in a fluctuating environment and regulated by different kinds of random factors (kinases, ligands, eRNA, etc.) [27,46,11,9,20,13,34]. Owing to the small number of such molecules and random births and deaths of molecules, the stochasticity of the input signal is unneglectable [20]. Increasing investigation has shown that external signals play an important role in cellular function, for instance, a new treatment of HIV was proposed by using input noise [8]. An oscillating signal modulation may increase the mean of protein and decrease the noise of protein compared with a constant signal modulation [64]. A random signal modulation may decrease the rate of gene state switching compared with a constant signal modulation [22]. Therefore, the fluctuating environment can indeed result in different conclusions compared with homogeneous environment. So, taking the effect of external stimulus into consideration, many of the previous conclusions would be modified and the mechanism how input signals impact on FPT remains elusive. In order to gain more insight into the regulation mechanism of input signal on FPT, we discuss the effects of random input signal on FPT.

Quantifying the effect of signal modulation on FPT is an important step towards understanding cellular functional variability. In order to make up for the lack that the previous studies ([10,57,17,53,52,45,66]) did not consider the influence of input signal on FPT moments, we are going to investigate the mechanism of FPT in the case of signal regulation by using gene expression models. Based on the fact that gene expression is almost in a geometric bursting manner both in mRNA synthesis [7,12,38,70,6] and in protein synthesis [70,36,4,44] from a single mRNA, further research [24] revealed that the total number of proteins produced in a single burst event follows conditional geometric distribution. We distinguish signal into two types: noiseless and noisy signals. The regulation ways of noisy signal on gene expression are further classified into two regulation ways: burst frequency modulation and burst size modulation. Here, burst frequency modulation means that input signal regulates the burst times per unit time, whereas burst size modulation means that input signal regulates the number of proteins in a bursty event [66]. The main results of this paper are as follows. Firstly, analytical calculations of FPT moments are derived in each case of noiseless signal modulation and noisy signal modulation. This is the first time to obtain theses analytical results in lately literatures. In addition, our numerical results show that random input signal tends to increases the mean and noise of FPT compared with constant input signal for a given protein threshold. Our numerical results also show that burst size modulation tends to produce a larger noise of FPT than that of burst frequency modulation. In conclusion, our results show that random signal may prolong the mean of FPT. This implies that the randomness of environment may prolong the latency of some diseases (such as HIV [5]). Given the prevalence of random signal, illuminating the effects of stochasticity of signal molecules on FPT can aid understanding of their regulatory roles in biological processes.

---

2. Gene expression models and analytical formulas for FPT moments

In order to clearly reveal the mechanism of how input signal modulates FPT, we distinguish input signal into two cases: noiseless and noisy signals. In the case of noiseless signal, a gene expression model is given in Figure 1(A): A gene produces mRNA in a bursty fashion with rate $\lambda$ and mean burst size $b_{1}$, and each mRNA degrades with rate $\delta_{m}$; Production of protein from a mRNA is in a bursty manner with rate $k_{p}$ and mean burst size $b_{2}.$ Based on the study [24], burst size of proteins from a transcription burst event, say $Y$, follows the following conditional geometric distribution,

$Q(0)=Pr\{Y=0\}=\frac{1}{1+b_{1}}+\frac{b_{1}}{1+b_{1}}\frac{1}{1+(b_{1}+1)b_{2}}, \ $

$Q(n)=Pr\{Y=n\}=\frac{b_{1}}{1+b_{1}}\frac{(b_{1}+1)^{n}b_{2}^{n}}{(1+(b_{1}+1)b_{2})^{n+1}}, n=1, 2, 3, \cdots$

The reference [41] indicated that the loss of highly stable proteins is mainly due to dilution through growth and cell division. This paper mainly discusses the dynamics of FPT before cell division. Therefore the degradation of protein is not considered in the following models.

According to references [53,56,58,25], the waiting time between two consecutive transcription burst events obeys exponential distribution. In addition, the lifetime of mRNA is far shorter than the cell cycle and mRNA degrades instantaneously after producing protein in a burst manner [53]. Thus, we only need to consider gene expression wherein the interval time between two consecutive protein burst events (the time of protein burst equal to transcription time) follows an exponential distribution with parameter $\lambda$, and the number of protein burst follows conditional geometric distribution with mean $b_{1}b_{2}$.

Randomness in the level or localization of regulation factors, such as Calcium [42], eRNA [33,40], Bicoid [16,14] and NF-kB [68,67], has been observed in diverse gene regulation network. So far, a number of studies have been focusing on the influences of input signal on regulated-gene product. But the impact of random signal on FPT has been poorly understood. Our main purpose is to quantify the impacts of input signal on FPT by employing a gene regulation model wherein random signal only regulates burst frequency (burst size) shown in ① (②) in Figure 1(B). Here, the number $z(t)$ of random signal molecules obeys the equilibrium distribution of a Markov birth-death process [22] characterized by birthing at a constant rate $\alpha$ per unit time and degrading with a constant rate $\beta$ per unit time.

In order to study the effects of random signal on FPT, we distinguish signals into two types: noiseless and noisy signals, to show theoretical analysis. For both cases, finding random characters of FPT, such as mean, variance and noise, become a common interest in understanding the stochastic properties of gene expression. We will mainly concentrate on the calculation of the analytical expression for the mean, variance and noise of FPT in each case of noiseless and noisy signal regulations.

The intrinsic randomness of biochemical reactions leads to that the protein count $P(t)$ at time $t$ is a stochastic process. In this paper, we set $P(0)=0.$ Notably, FPT can be defined as the minimal time that the protein count $P(t)$ up to a given threshold $m$ $(m\neq 0)$. Let $F_{m}$ represent FPT at which the number of protein reaches the critical threshold $m.$ Therefore, we can write $F_{m}$ as follows

$F_{m}=\inf\{t:P(t)\geq m\}.$

Further, let $P_n$ be the protein count after $n$ gene expression events, and $N_{m}$ be the random variable which represents the minimum number of protein burst event that takes for protein count to reach the given threshold $m.$ Then

$N_{m}=\inf\{n:P_{n}\geq m\}.$

Let $T_{i}$ be the waiting time from the $(i-1)^{th}$ to the $i^{th}$ protein burst event. Then $T_{i}, i=1, 2, 3, \cdots$ are independent identical distribution. Furthermore, according to the definitions of $N_{m}$ and $T_{i},$ we have

$F_{m}=\sum\limits_{i=1}^{N_{m}}T_{i}.$

Let $Pr\{D\}$ represent the probability of event $D$. Next, we focus on finding the mean, variance and noise of FPT. The basic idea is first to introduce the properties of conditional expectation, and then solve a set of recurrence equations. In addition, based on introducing the Kolmogorov backward equation, probability generating functions, we solve two differential equations. The overall procedure for finding such random characters is technical.

By using the property for conditional expectation, we can obtain

$\begin{split} E(F_{m}) =&E\left(\sum\limits_{i=1}^{N_{m}}T_{i}\right)=E\left(E\left( \sum\limits_{i=1}^{N_{m}}T_{i}|N_{m}\right)\right)\\ =&\sum^\infty_{n=1}E\left( \sum^{N_{m}}_{i=1}T_{i}|N_{m}=n\right) P_{r}(N_{m}=n) \\ =&E( T_{1})E(N_{m}) \end{split}$

(1)

and

$\begin{split} Var(F_{m}) =&E \left( \left(\sum^{N_{m}}_{i=1}T_{i}\right)^2\right) -E^2 \left(\sum^{N_{m}}_{i=1}T_{i}\right)\\ =&\sum^\infty_{n=1}\left[Var\left( \sum^n_{i=1}T_{i}\right)+E^{2}\left( \sum^n_{i=1}T_{i}\right) \right]P_{r}\{N_{m}=n\}-E^{2}( T_{i})E^{2}(N_{m})\\ =&Var(T_{1})E(N_{m})+E^2(T_{1})Var(N_{m}). \end{split}$

(2)

Let $\eta$ represent the noise of $F_{m}.$ Then

$\begin{split} \eta=\frac{Var(F_{m})}{E^{2}(F_{m})}=&\frac{Var(T_{1})}{E^{2}(T_{1})E(N_{m})}+\frac{Var(N_{m})}{E^{2}(N_{m})}. \end{split}$

(3)

In order to calculate the mean, variance and noise of FPT, we only need to give exact formulae for the mean and variance of both $T_{i}$ and $N_{m}$ by (1), (2) and (3).

Next, we focus on the mean and noise of FPT in each case of noiseless and noisy signal regulations.

---

2.1. Mean and noise of FPT with noiseless case

Genetically identical cell populations exposed to the same extracellular environment exhibit considerable variability in gene expression [54]. The same extracellular environment means that transcription rate $\lambda$ and the mean $b_{1}b_{2}$ of burst size in gene expression are constant. Hence variations in gene product result from random births and deaths of individual molecules. This part consider the random characters of FPT in the case of homogenous environment.

Since $T_{i}$ $(i=1, 2, 3, 4, \cdots)$ are independent and identically exponential distribution, we obtain

$E(T_{i})=\frac{1}{\lambda}, Var(T_{i})=\frac{1}{\lambda^{2}}.$

Next, we concentrate on the first and second moments of $N_{m}$. Let $X_{i}$ $(i=1, 2, \cdots)$ be the random variable denoting the number of proteins produced in the $i^{th}$ burst expression. Notably, $X_{i}$ are independent and identically distribution. Then we have

$Pr\{X_{i}=0\}=Q(0)=\frac{1}{1+b_{1}}+A,$

$Pr\{X_{i}=n\}=Q(n)=AB^{n}, n=1, 2, 3, \cdots$

where

$A=\frac{b_{1}}{1+b_{1}}\frac{1}{1+(b_{1}+1)b_{2}}~~\text{and}~~ B=\frac{(b_{1}+1)b_{2}}{1+(b_{1}+1)b_{2}}.$

By the definition of random variable $N_{m}$, the probability of stochastic event $\{N_{m}=n\}$ $(n=1, 2, 3, \cdots)$ can be denoted as

$\begin{array}{ll} Pr\{N_{m}=1\}=&Pr\{X_{1}\geq m\}=1-Pr\{X_{1} < m\}=1-\sum\limits_{k=0}^{m-1}Q(k) \end{array}$

and

$\begin{array}{lll} Pr\{N_{m}=n\}&=&Pr\{P_{n}\geq m, P_{n-1}\leq m-1\}\\ &=&\sum\limits_{k=0}^{m-1}Pr\{P_{n}\geq m, P_{n-1}\leq m-1, X_{1}= k\}\\ &=&\sum\limits_{k=0}^{m-1}Pr\{P_{n}\geq m, P_{n-1}\leq m-1|X_{1}= k\}Pr\{X_{1}= k\}\\ &=&\sum\limits_{k=0}^{m-1}Pr\{P_{n-1}\geq m-k, P_{n-2}\leq m-k-1\}Pr\{X_{1}= k\}\\ &=&\sum\limits_{k=0}^{m-1}Pr\{N_{m-k}=n-1\}Pr\{X_{1}= k\}, \end{array}$

for $n=2, 3, 4, \cdots.$ In addition, we can get $E(N_{m})$ by the following calculation

$\begin{split} E(N_{m})=&\sum\limits_{n=1}^{\infty}nPr\{N_{m}=n\}\\ =&Pr\{N_{m}=1\}+\sum\limits_{n=2}^{\infty} \sum\limits_{k=0}^{m-1}nPr\{N_{m-k}=n-1\}Pr\{X_{1}= k\}\\ =&1+\sum\limits_{k=0}^{m-1}E(N_{m-k})Q(k). \end{split}$

(4)

It is easy to show by induction for a given threshold $m,$

$\begin{split} E(N_{m})=& ma+ab_{2}+1. \end{split}$

(5)

Similarly to the aforementioned derivation, we can also get that

$E(N^{2}_{1})=2E^{2}(N_{1})-E(N_{1}), ~ Var(N_{1})=E(N_{1})(E(N_{1})-1).$

It implies that

$E(N^{2}_{1})=[2a(1+b_{2})+1][a(1+b_{2})+1], ~~ Var(N_{1})=a(1+b_{2})[a(1+b_{2})+1].$

For $k=2, 3, 4\cdots,$ we can derive the following recurrence formula of $E(N^{2}_{m}),$

$\begin{split} E(N^{2}_{m})=E(N_{1})\sum\limits_{k=1}^{m-1}E(N^{2}_{m-k})Q(k) +2E(N_{m})E(N_{1})-E(N_{1}). \end{split}$

(6)

Therefore, we obtain

$\begin{split} E(N^{2}_{m})=&E(N^{2}_{m-1})+2[ma^{2}+2a^{2}b_{2}+2a]-a. \end{split}$

(7)

By using mathematical induction, we can find the formula of $E(N^{2}_{m})$ for a given threshold $m,$

$\begin{split} Var(N_{m})=&(2a^{2}b_{2}+a^{2}+a)m+ab_{2}(ab_{2}+1). \end{split}$

(8)

Thus, we get

$\begin{split} E(F_{m})=&\frac{1}{\lambda}[ma+ab_{2}+1] \end{split}$

(9)

and

$\begin{split} Var(F_{m})=&\frac{1}{\lambda^{2}}[ma+ab_{2}+1] +\frac{1}{\lambda^{2}}[(2a^{2}b_{2}+a^{2}+a)m+ab_{2}(ab_{2}+1)]\\ =&\frac{1}{\lambda^{2}}[(2a^{2}b_{2}+a^{2}+2a)m+ab_{2}(ab_{2}+2)+1]. \end{split}$

(10)

Further, we can obtain the noise $\eta_{c}$ of $F_{m}$ as follows

$\begin{split} \eta_{c}=&\frac{(2a^{2}b_{2}+a^{2}+2a)m+ab_{2}(ab_{2}+2)+1}{a^{2}m^{2}+2a(1+ab_{2})m+ab_{2}(ab_{2}+2)+1}. \end{split}$

(11)

Interestingly, the noise of $F_{m}$ is independent of burst frequency. Note that (9), (10) and (11) are exact for any values of the system parameters, and hence can be directly used for numerical calculations.

---

2.2. Mean and noise of FPT with noisy case

In a single cell, the creation of mRNA and protein occurs in a bursty, intermittent manner. Burst frequency and burst size are two main indexes of burst dynamics [55]. The frequency and size of bursts affect the magnitude of noise [50] and the modality of probability distribution of protein [23], and even may play a critical role in the realization and switching of biological functions [55]. Meanwhile, recent advances in experimental technology have confirmed that the stochastic nature of cell-signaling molecules, such as tumor necrosis factor α(TNF) and adenosine triphosphate (ATP), influences burst frequency and/or burst size of gene in vitro and vivo cells [26]. The regulation ways of input signals on gene expression are generally classified into three different but common modes [54,66,55,32]: burst frequency regulation (without regard to burst size), burst size regulation (without regard to burst frequency) and simultaneous regulation on both burst frequency and burst size. There have been some theoretical studies on input signal modulations [22,54,45,66,43,39], for instance, for burst frequency modulation, random input signal can cause stochastic focusing [43], make regulated-gene product generate a bimodal steady state output [39] and increase the switch rate [22]. For burst size modulation, random input signal dramatically increases noise compared with burst frequency modulation [54]. Therefore, burst frequency regulation and burst size regulation on gene expression result in different effects. Since litter was known about the effects of such two regulating ways of random input signal on FPT, the study on impacts of input signal on FPT is of great significance. In this subsection, we only focus on different effects of burst frequency regulation and burst size regulation on FPT. Of course, burst frequency and burst size may be regulated by the same signal, such as, trichostatin A can regulate simultaneously burst frequency and burst size [60]. For the final case of simultaneous regulation in both burst frequency and burst size, we will investigate its regulation effects on FPT in another paper.

Since the number $z(t)$ of input signal molecules follows the equilibrium distribution of a Markov birth-death process characterized by birthing at rate $\alpha$ and degrading with rate $\beta$, $z(t)$ follows poisson distribution,

$P_{r}\{z(t)=k\}= e^{-\frac{\alpha}{\beta}}\frac{({\alpha}/{\beta})^{k}}{k!}, \ k=0, 1, 2, 3, \cdots.$

---

2.2.1. The case of burst frequency modulation

Cells are often exposed to changing environment. They sense such changing environment with cell-surface receptors and/or ion channels. This ultimately leads to the change of the concentration of regulation factors. Although some studies on regulation factors have been recently begun, such as long non-coding RNAs (lncRNAs) and microRNA, their potential functions and mechanisms on gene expression still is incompletely understood [48].

In this subsection, we consider the case when the random factor $z(t)$ only regulates burst frequency without being regarded to burst size. The modulation way reflects the fact that burst rate is not longer a constant $\lambda$ but given by $\lambda z(t)$ [54] and burst size is not affected by signal. That is to say, proteins are produced at rate $\lambda z(t)$ in a burst manner. Given that signal molecules do not affect burst size, the number of proteins produced in a burst event is independent of $z(t).$ Thus, we have

$\begin{split} E(N_{m})= ma+ab_{2}+1;~Var(N_{m})= (2a^{2}b_{2}+a^{2}+a)m+ab_{2}(ab_{2}+1). \end{split}$

(12)

Next, we calculate the first two moments of interval time $T_{i}$ $(i=1, 2, 3, \cdots)$. Note that input signal obeys the equilibrium distribution of a Markov birth-death process, this implies that $T_{i} (i=1, 2, 3, \cdots)$ are independent and identically distribution. Since the mean of $T_{i}$ cannot be expressed simply by $\lambda$ and $E(z(t))$, we should consider the probability density function of $T_{i}$. Now, suppose that the number $z(t)$ of signal molecules equals to $n$, we can write the Kolmogorov backward equation for $\rho_{i, n}(t)$ as follows,

$\begin{split} \rho_{i, n}'(t)=&-(\alpha+\beta n+\lambda n)\rho_{i, n}(t)+\alpha \rho_{i, n+1}(t)+\beta n \rho_{i, n-1}(t), \end{split}$

(13)

where $\rho_{i, n}(t)$ denotes the probability density function for $T_{i}$. We multiply both sides of (13) by $e^{-\frac{\alpha}{\beta}}\frac{({\alpha}/{\beta})^{n}}{n!}x^{n}t,$ integrate $t$ and sum over all nonnegative integer $n,$ then we can get

$\begin{split} (\beta x+\lambda x-\beta) G'_{1}(x)+\alpha (1-x )G_{1}(x)= e^{-\frac{\alpha}{\beta}(1-x)}, \end{split}$

(14)

where

$G_{1}(x)=e^{-\frac{\alpha}{\beta}}\sum\limits^{\infty}_{n=0} x^{n}\frac{1}{n!}\left(\frac{\alpha}{\beta}\right)^{n}\int_{0}^{\infty}t\rho_{i, n}(t)dt.$

From (14), we get

$\begin{split} G_{1}(x)=\frac{x_{0}}{\beta}(x-x_{0})^{-\frac{\alpha}{\beta}(1-x_{0})x_{0}} e^{\frac{\alpha}{\beta}x_{0}x}\int_{x_{0}}^{x}e^{-\frac{\alpha} {\beta}[1-t+x_{0}t]}(t-x_{0})^{\frac{\alpha}{\beta}(1-x_{0})x_{0}-1}dt \end{split}$

(15)

Note that $E(T_{i})=\sum\limits_{n=0}^{\infty}e^{-\frac{\alpha}{\beta}} \frac{({\alpha}/{\beta})^{n}}{n!}\int_{0}^{\infty}t\rho_{i, n}(t)dt=G_{1}(1)$, combing with (15) we have

$\begin{split} E(T_{i})=&\frac{x_{0}}{\beta}(1-x_{0})^{-\frac{\alpha}{\beta}(1-x_{0})x_{0}} e^{-\frac{\alpha}{\beta}(x_{0}-1)^{2}}\int_{x_{0}}^{1}e^{\frac{\alpha} {\beta}(1-x_{0})(x-x_{0})}(x-x_{0}) ^{\frac{\alpha}{\beta}(1-x_{0})x_{0}-1}dx, \end{split}$

(16)

where $x_{0}=\frac{\beta}{\beta+\lambda}.$

Therefore, by (1), (12) and (16) we have

$\begin{split} E(F_{m})=&G_{1}(1)(ma+ab_{2}+1). \end{split}$

(17)

To calculate the second moment of waiting time, we multiply both sides of (13) by $e^{-\frac{\alpha}{\beta}}\frac{(\frac{\alpha}{\beta})^{n}}{n!}x^{n}t^{2}$, integrate over all $t$ and sum over all integer $n$, then we get

$\begin{split} (\beta x +\lambda x-\beta) G'_{2}(x)+\alpha (1-x )G_{2}(x)=2G_{1}(1), \end{split}$

(18)

where $G_{2}(x)=e^{-\frac{\alpha}{\beta}}\sum\limits^{\infty}_{n=0}x^{n}\frac{1}{n!} (\frac{\alpha}{\beta})^{n}\int_{0}^{\infty}t^{2}\rho_{i, n} dt.$ From (18), we get

$\begin{split} G_{2}(x)=&\frac{2x_{0}G_{1}(1)}{\beta}(x-x_{0})^{-\frac{\alpha}{\beta}(1-x_{0})x_{0}} e^{\frac{\alpha}{\beta}x_{0}x}\int_{x_{0}}^{x}e^{-\frac{\alpha} {\beta}x_{0}t}(t-x_{0})^{\frac{\alpha}{\beta}(1-x_{0})x_{0}-1}dt \end{split}$

(19)

Note that $E(T^{2}_{i})=\sum\limits_{n=0}^{\infty}e^{-\frac{\alpha}{\beta}} \frac{({\alpha}/{\beta})^{n}}{n!}\int_{0}^{\infty}t^{2}\rho_{i, n}(t)dt=G_{2}(1),$ combing with (19) we have

$\begin{split} E(T^{2}_{i})=&2G_{1}(1)\int_{x_{0}}^{1}\frac{1}{((\beta+\lambda)x-\beta)} e^{\frac{\alpha}{\beta+\lambda}(1-x)}\left(\frac{x-x_{0}} {1-x_{0}}\right)^{\frac{\alpha \lambda}{(\beta+\lambda)^{2}}}dx\\ =&G_{1}(1)\int_{x_{0}}^{1}\frac{2x_{0}}{\beta(x-x_{0})} e^{\frac{\alpha}{\beta}x_{0}(1-x)}\left(\frac{x-x_{0}} {1-x_{0}}\right)^{\frac{\alpha x_{0}(1-x_{0})}{\beta}}dx. \end{split}$

(20)

On the basis of (2), we obtain

$\begin{split} Var(F_{m})=&(G_{2}(1)-G^{2}_{1}(1))E( N_{m}) +G_{1}^{2}(1)Var(N_{m}). \end{split}$

(21)

All the analytical results in this part are exact but some of them are not intuitive because of these integrals.

---

2.2.2. The case of bursting size modulation

In this subsection, we consider the case when the random factor $z(t)$ only regulates the burst size. The modulation way reflects the fact that the mean of transcription burst size is $b_{1}z(t)$ and burst rate is a constant $\lambda$. That is, proteins are produced at rate $\lambda$ in a burst manner and the number of proteins synthesized at time $t$ when $z(t)=k$ follows the following probability distribution:

$\begin{align*} Q_{k}(0)=\frac{1}{1+b_{1}k}+A_{k}, ~~ Q_{k}(n)=A_{k}B^{n}_{k}, ~n=1, 2, 3, \cdots \end{align*}$

where $A_{k}=\frac{kb_1}{1+kb_1}\frac{1}{1+(kb_1+1)b_{2}}, ~~ B_{k}=\frac{(kb_1+1)b_{2}}{1+(kb_1+1)b_{2}}, ~~ k=1, 2, 3, \cdots.$

Given that signal molecules do not affect the burst rate, the interval time between two consecutive burst events is independent of $z(t).$ So we have

$\begin{split} E(T_{i})=&\frac{1}{\lambda}, ~~Var(T_{i})=\frac{1}{\lambda^{2}}. \end{split}$

(22)

In order to obtain the mean, variance and noise of $F_{m}$, we will derive recurrence formulae for calculating $E(N_{m})$ and $E(N^{2}_{m})$ by means of conditional expectation. Now, suppose that the $i^{th}$ protein burst event occurs at time $t_{i}$, we have

$\begin{split} E(N_{m})=&E[E(N_{m}|z(t_{1}))] =\sum\limits_{k=0}^{\infty} e^{-\mu}\frac{\mu^{k}}{k!}E[N_{m}|z(t_{1})=k], \end{split}$

(23)

where $\mu=\frac{\alpha}{\beta}$. For $z(t_{1})=0,$ we have $P_{r}\{N_{m}=0|z(t_{1})=0\}=0$ and $P_{r}\{N_{m}=1|z(t_{1})=0\}=0.$ According to the properties of conditional expectation, we have

$\begin{split} &E[N_{m}|z(t_{1})=0]=\sum\limits_{n=2}^{\infty}n P_{r}\{N_{m}=n|z(t_{1})=0\}\\ &=\sum\limits_{n=2}^{\infty}nP_{r}\left\{\sum\limits_{i=1}^{n}X_{i}\geq m, \sum\limits_{i=1}^{n-1}X_{i}\leq m-1|z(t_{1})=0\right\}\\ &=\sum\limits_{n=2}^{\infty}n P_{r}\left\{\sum\limits_{i=1}^{n}X_{i}\geq m, \sum\limits_{i=1}^{n-1}X_{i}\leq m-1|X_{1}=0, z(t_{1})=0\right\}P_{r}\{X_{1}=0|z(t_{1})=0\}\\ &+\sum\limits_{n=2}^{\infty}n P_{r}\left\{\sum\limits_{i=1}^{n}X_{i}\geq m, \sum\limits_{i=1}^{n-1}X_{i}\leq m-1|X_{1}\geq 1, z(t_{1})=0\right\}P_{r}\{X_{1}\geq 1|z(t_{1})=0\}\\ &=\sum\limits_{n=2}^{\infty}n P_{r}\left\{\sum\limits_{i=2}^{n}X_{i}\geq m, \sum\limits_{i=2}^{n-1}X_{i}\leq m-1\right\}\\ &=1+E(N_{m}). \end{split}$

(24)

Combing (23) and (24), we obtain the following formula

$\begin{split} E(N_{1})=&\sum\limits_{k=0}^{\infty}e^{-\mu}\frac{\mu ^{k}}{k!}E(N_{1}|z(t_{1})=k)\\ =&e^{-\mu}(1+E(N_{1}))+e^{-\mu}\sum\limits_{k=1}^{\infty}\frac{\mu ^{k}}{k!} (1-Q_{k}(0))\\& +e^{-\mu}\sum\limits_{k=1}^{\infty}\frac{\mu ^{k}}{k!}\sum\limits_{n=2}^{\infty}n[P_{r}\{N_{1}=n-1\}Q_{k}(0)]\\ =&e^{-\mu}(1+E(N_{1}))+e^{-\mu}\sum\limits_{k=1}^{\infty}\frac{\mu ^{k}}{k!} (1-Q_{k}(0))+e^{-\mu}\sum\limits_{k=1}^{\infty}\frac{\mu ^{k}}{k!} (1+E(N_{1}))Q_{k}(0)\\ =&1+\sum\limits_{k=0}^{\infty}a_{k}(\mu)E(N_{1}), \end{split}$

(25)

where $a_{k}(\mu)=e^{-\mu}\frac{\mu^{k}}{k!}(\frac{1} {1+kb_1}+A_{k}).$ Furthermore, we can obtain $E(N_{1})$ as follows

$\begin{split} E(N_{1})=\frac{1}{g(\mu)}, \end{split}$

(26)

where $g(\mu) =e^{-\mu}\sum\limits_{k=1}^{\infty}\frac{\mu^{k}}{k!}\frac{kb_1b_{2}} {1+(1+kb_1)b_{2}}.$ From (23), we obtain the next formula

$\begin{split} E(N_{m})=&e^{-\mu}\sum\limits_{k=0}^{\infty}\frac{\mu ^{k}}{k!}E(N_{m}|z(t_{1})=k)\\ =&e^{-\mu}[1+E(N_{m})] +e^{-\mu}\sum\limits_{k=1}^{\infty}\frac{\mu ^{k}}{k!}[1-\sum\limits_{j=0}^{m-1}Q_{k}(j)]\\ &\ \ \ +e^{-\mu}\sum\limits_{k=1}^{\infty}\frac{\mu ^{k}}{k!}\sum\limits_{n=2}^{\infty}n\sum\limits_{j=0}^{m-1} P_{r}\{N_{m-1}=n-1\}Q_{k}(j) \end{split}$

$\begin{split}=&e^{-\mu}[1+E(N_{m})]+e^{-\mu}\sum\limits_{k=1}^{\infty}\frac{\mu ^{k}}{k!} [1-\sum\limits_{j=0}^{m-1}Q_{k}(j)]\\ &\ \ \ +e^{-\mu} \sum\limits_{k=1}^{\infty}\frac{\mu ^{k}}{k!} \sum\limits_{j=0}^{m-1}[1+E(N_{m-j})]Q_{k}(j)\\ =&1+\sum\limits_{k=0}^{\infty}a_{k}(\mu)E(N_{m}) +\sum\limits_{k=1}^{\infty}\sum\limits_{j=1} ^{m-1}c_{k, j}(\mu)E(N_{m-j}), \end{split}$

(27)

where $c_{k, j}(\mu)=e^{-\mu}\frac{\mu^{k}}{k!}A_{k}B^{j}_{k}.$ On the basis of (26) and (27), we can obtain the recurrence formula of $E(N_{m})$ as follows

$\begin{split} E(N_{m})=&E(N_{1})+E(N_{1})\sum\limits_{k=1}^{\infty} \sum\limits_{j=1}^{m-1}c_{k, j}(\mu)E(N_{m-j}). \end{split}$

(28)

Similar to the previous analyses, we obtain the exact formula of $E(N^{2}_{1})$

$\begin{split} E(N^{2}_{1})=&\frac{2-g(\mu)}{g^{2}(\mu)} \end{split}$

(29)

and the recurrence formula of $E(N^{2}_{m})$

$\begin{split} E(N^{2}_{m})=&-E(N_{1})+2E(N_{1})E(N_{m}) +\sum\limits_{k=1}^{\infty}\sum\limits_{j=1}^{m-1}c_{k, j}(\mu)E(N^{2}_{m-j})E(N_{1}). \end{split}$

(30)

Now, in terms of (1), (2), (28) and (30), we can obtain the analytical expressions for the mean and variance of $F_{m}$. The recurrence formulae of $E(F_{m})$ and $Var(F_{m})$ are omitted because of their complex forms.

---

3. Main results

The above analytical results about the mean and noise of FPT, in principle, lay a solid foundation for understanding of how the two regulations ways of input signal affect FPT moments. Now we perform numerical calculations to give intuitive results for these impacts (shown in Figures 2 and 3).

---

3.1. Comparison between random and constant input signals

Recently, some studies have confirmed that biological fate selections are driven by the levels of protein. For example, differentiation in B. subtilise [57], lysis in bacteriophage $\lambda$ [10,53], viral latency of human immunodeciency virus (HIV)-1 [5,37], latency of herpes simplex virus (HSV, subfamily-α) in infected neurons [63,29], latency in Kaposi's sarcoma-associated herpesvirus (KSHV) [19,65], viral latency in cytomegalovirus infection in the lung [18] and so on, these all depend on the accumulation of proteins reaching a given threshold. Cell-to-cell variations in protein level mainly result from the randomness of external environment and biochemical reactions involved in gene expression [45]. However, the previous studies [10,53,52] rarely considered the impacts of random signal on FPT.

Now, we will interpret how random signal regulation affects the mean and noise of FPT. More precisely, compared with constant input signal, random input signal tends to increase the mean and noise of FPT in the case of burst frequency modulation. For burst size modulation, the conclusion holds only when the mean of transcription times is larger than a threshold. These results provide theoretical guidance for studies of cell fate decision caused by protein level upping to the critical threshold, such as lysis time of bacteriophage $\lambda$ and latency of HIV and so on. It reveals that both lysis time of bacteriophage $\lambda$ and latency of HIV in homogeneous environment may become smaller and stabler than in random environment.

Next, we will perform numerical calculations to further reveal quantitative effects of random signal on FPT shown in Figure 2(A, B). For burst frequency modulation, we observe the following three conclusions. The first is that the mean of FPT under random input signal modulation is larger than that under constant input signal modulation. This implies that random input signal may suppress expression of protein via prolonging waiting time. The second is that the mean of FPT is a monotonically decreasing function of input signal intensity, without regard to modes of input signals. The difference in the mean of FPT between noiseless modulation and noisy modulation is smaller with the increase of input signal strength. The conclusion is qualitatively invariant, independent of its related parameters. The final one is that the noise of FPT under random input signal modulation becomes larger than that under constant input signal modulation. The difference in the noise of FPT between them tends to become larger with the increase of small input signal strengths and become smaller after input signal strength exceeding a threshold.

For burst size modulation, we have the following four conclusions shown in Figure 2(C, D). The first is that the mean of FPT under random input signal modulation is larger than that under constant input signal modulation. The second is that the mean of FPT is a monotonically decreasing function of input signal intensity without being regarded to modes of input signals. The difference in the mean of $F_{m}$ between them becomes smaller with the increase of input signal intensity. The third is that the noise of FPT under random input signal modulation tends to be larger than that under constant input signal modulation. The final one is that the noise of FPT is a monotonically increasing function of input signal strength in each case of noiseless and noisy signal modulation. The difference in the noise of FPT between them becomes very smaller when the strength is beyond a threshold.

Summarizing the above analyses, we can conclude that random input signals play an important role in increasing the mean and noise of FPT no matter how input signal regulates gene expression. Therefore, constant input signal regulation can better modulate FPT, compared with the corresponding random signal regulation.

---

3.2. Comparison between random burst size modulation and random burst frequency modulation

The effects of burst size regulated by random signal on gene product are somewhat different from that burst frequency regulated by random signal. For example, Singh et al. [54] found that burst size regulation can enlarge both intrinsic and extrinsic noises but burst frequency regulation only increases extrinsic noise. However, their effects of the two regulations ways on FPT in gene regulation is still not clear. So, we compare the effects on FPT caused by burst frequency modulation and burst size modulation shown in Figure 3.

We observe the following two results by performing numerical calculations. On one hand, we observe that the mean of $F_{m}$ in the case of burst frequency modulation is larger than that in the case of burst size modulation and the difference between them becomes smaller with the increase of input signal intensity. On the other hand, the noise of $F_{m}$ in the case of burst frequency modulation is smaller than that in the case of burst frequency regulation and the difference between them becomes larger with the increase of input signal intensity. In conclusion, burst frequency modulation can better modulate FPT than burst size modulation, remarkably in large strength cases. These conclusions may reveal that both lysis time of $\lambda$ phage and latency of HIV in the case of burst frequency modulation may be stabler than in the case of burst size modulation.

---

4. Discussion

The inherent randomness of biochemical reactions can lead to cell-to-cell variability in the timing of proteins crossing a given threshold even in homogenous environment, let alone cells are often exposed to the changing environment. Hence, cells are always affected by various random signal molecules, but how signals quantitatively and qualitatively influence the timing of cellular key events, such as lysis, B.subtilis differentiation and HIV latency, is till not clear. To investigate the impact of signal on the timing of cellular key events, we considered the first-passage time of proteins up to a given threshold. Here, we have systematically analyzed a gene expression model where input signal only regulates burst frequency (or burst size) of gene expression. By analysis, we obtained the following main results:

a) Analytical calculations for FPT moments, either noiseless or noisy signal regulation, are derived.

b) Compared with constant input signal, random input signal tends to increase the mean and noise of FPT.

c) Compared with burst size modulation, burst frequency modulation tends to increase the mean of FPT and decrease the noise of FPT.

The analytical and numerical methods used in this paper allow us to explore how stochastic fluctuations of input signals affect FPT and can be extended to similar multistate gene regulation models. Realization and changing of some biological functions tend to depend on the fact that protein count reaches a given threshold, such as bacterial cell division in biological systems [1,15]. In the future, we will study how to combine our theoretical research with specific biological problems, such as cell division and latency of viruses.

---

Acknowledgments

We are greatly indebted to Dr. Zhanjiang Yuan, Linchao Hu, Qi Wang and Kunwen Wen for their kind suggestions in the course of our research.

---