Modeling stochastic gene expression: From Markov to non-Markov models

1.   Introduction

Mathematical models are a powerful tool for modeling and analysis of complex gene expression processes. An important task of systems biology is to design reasonable mathematical models that faithfully describe the dynamics of gene expression. Any theoretical model of gene expression established based on the first principles needs to make assumptions. This is because the complexity of gene expression would make it difficult to specify what transitions can be condensed into effective reaction steps and what concentrations can be absorbed into reaction rate constants. It seems always difficult to tell if the correct principles are addressed when modeling a complex process of gene expression. However, it is possible to understand simple mechanistic models of gene expression at the initial stage.

Most published gene models have independently focused on the fundamental genetic information flow described by the central dogma in biology (i.e., the information stored in DNA is first transcribed into mRNAs, which are then translated into functional proteins). These models include (but do not limit to) transcription models ^[1,2,3,4,5], transcription and translation models ^[6,7,8,9], gene regulatory models ^{[10,11,12,13,14]}, Markov models of gene expression ^[15], non-Markov (or queuing) models of gene expression ^[16,17]. The underlying fundamental principle may be extended, depending on the dissection of molecular details involved in gene expression. A common way of the extension is to implicitly include all other processes in effective reaction rates of gene models ^{[2,4,5,6,7,8,12,13,17]}. With the development of single-cell and single-molecule measurement technologies, more molecular details that can significantly impact gene-product levels have been revealed. Here we state key molecular events occurring in gene expression processes, which will be partially incorporated in gene models of this paper.

First, activation and inactivation of a gene can have many different molecular causes, such as alternative splicing, transcription initiation, recruitment of polymerases, dissociation of repressors, association of activators, histone modifications, and chromatin remodeling as well as recruitment of transcription factors (TFs). In particular, the activation of a gene depends on the chromatin template that accumulates over time until the promoter becomes active ^{[2,18,19,20,21,22]}, implying a multistep process. The details vary from gene to gene and from cell to cell, and can include random transitions between different promoter activity states. As a first approximation, the reaction kinetics of promoters in prokaryotic and eukaryotic cells may be described as a random telegraph process with constant rates or equivalently with exponential waiting-time distributions ^[23].

Second, transcription or translation is often assumed to follow a Poisson process, i.e., the number of mRNA or protein molecules produced per time unit is often assumed to follow a Poisson distribution. This assumption may or may not be biologically reasonable. Complex regulation may affect transcriptional processes: (1) The binding of TFs to DNA sites may change the promoter structure, thereby impacting the efficiency of transcription; (2) Replication may alter the chromatin structure, giving rise to abrupt changes in expression rates; (3) mRNAs’ entry into or exit from the cellular cytoplasm may impact reaction kinetics.

Third, feedback regulation is not exceptional but is ubiquitous in gene expression systems. Positive and negative feedback loops are common regulatory forms found in biological signaling systems. For example, in mammals, more than 3000 signaling proteins and more than 15 second messengers form hundreds and thousands of cell-specific signaling pathways ^[24]. Many of these signaling pathways may have numerous upstream regulators and downstream targets, which altogether constitute a huge web of connectivity within and between signaling pathways ^[24]. The presence of multiple feedback loops poses a challenge to understanding how receptor inputs control the cellular behavior.

Fourth, apart from heredity-relevant molecular cues that lead to the so-called non-genetic heterogeneity, other possible contributors include fluctuations in many enzymes and substrates involved, cell cycle effects, and random partitioning of copies at cell division as well as fluctuations in methylation, histone modification, or more generally, epigenetic regulation ^[25]. These and other complications, which can also affect mRNA or protein abundances, will not be explicitly considered in gene models of this article.

Given the above complexity, how gene expression is modeled is challenging. According to the existing literature, we can divide the proposed models of gene expression into two large classes: Markovian and non-Markovian, referring to Figure 1 that shows the discrepancy between a Markov and a non-Markov gene models. For Markov models, it is needed to make the memoryless hypothesis, i.e., the stochastic motion of the reactants is influenced only by the current state but not by previous states. This hypothesis implies that waiting times for reaction events occurring in gene expression obey exponential distributions ^[26,27], and the reaction kinetics of gene expression are thus Markovian ^{[26,27,28,29]}. Experimental and theoretical studies have shown that the characteristic parameters in waiting-time distributions can have a significant impact on gene expression levels. However, reaction events occurring in gene expression may happen in a non-Markov fashion. Possible reasons include: (1) the complex control process of transcription initiation, which would involve many kinds of regulators (e.g., repressors, mediators, and TFs) as well as chromatin remodeling and changes in supercoiling, can generate non-exponential time intervals between transcription windows ^{[18,30,31,32,33,34,35]}; (2) The synthesis of an mRNA in general involves multiple intermediate reaction steps that would not have been specified or cannot be specified due to experimental technologies, creating molecular memory between individual events ^[18,36]; and (3) medium heterogeneity of gene expression may lead to broad waiting-time distributions between reaction events ^[37], thus influencing the efficiency of reactants in exploring their environments ^{[38,39,40,41,42]}. We will mainly focus on the second case, i.e., molecular memory is created by multistep biochemical processes involved in gene expression.

In most situations, biochemical reactions involved in gene expression are essentially single-molecule events, leading to stochastic fluctuations in mRNA and further protein levels. Thus far, there have been two kinetic modes of stochastic transcription that have been experimentally observed in individual cells: one mode is the Poissonian way where mRNAs are synthesized in a probabilistic manner with a probability that is uniform over time ^[43,44], and the other mode is the bursty way where mRNAs are produced in a bursty fashion with burst size following a distribution ^{[45,46,47,48]}. Each of the two ways can result in temporal, stochastic fluctuations in mRNA and protein numbers. This cell-to-cell variability is often referred to as gene expression noise.

In order to help the reader understand how stochastic gene expression is modeled and analyzed, we begin with Markov models and then extend them to non-Markov cases. Our analysis focuses on statistical quantities of protein, such as the mean protein level (or expectation), the protein noise defined as the ratio of the variance over the squared mean, and the Fano factor defined as the ratio of the variance over the mean. We will use these statistical quantities to characterize the protein noise even though a random variable is best characterized by its distribution. We highlight the contributions of noisy sources to protein levels by deriving analytical formulae for the protein noise. Additionally, we discus commonalities between these models as well as the possible linking between theoretical predictions and experimental observations.

2.   Markov models of stochastic gene expression

As mentioned in the introduction, Markov models assume that waiting times for reaction events are exponentially distributed. Here we review results for Markov models of stochastic gene expression, which are established based on the random telegraph model ^[18,19,49]. The telegraph model has long been used as a basis for several studies focusing on analysis of probability distributions for gene products and their statistical quantities as well as on inference of gene expression parameters based on experimental observations ^[50,51,52]. It is well known that random telegraph model has a good property: in the limit of transcriptional bursting, the arrival of bursts is a Poisson process. We start by the common on-off model of gene expression and then extend it to more realistic cases including multiple off mechanisms and feedback regulation of various forms.

2.1.   The common on-off model

In prokaryotic and eukaryotic cells, most genes are expressed in a bursty manner, namely gene product molecules are produced in short periods of high transcriptional activity (on) followed by long periods of inactivity (off). In order to model this expression manner, on-off models have been proposed ^{[48,53,54,55]}, which assume that the gene promoter has two activity states: one off state where the gene is not expressed and one on state where the gene is expressed with burst size $B$ following a probability distribution: ${\rm{prob}}\left\{ {B = i} \right\}$ with $i = 0, 1, 2, \cdots$. Denote by $a$ and $b$ the mean switching rates from off to off states and vice versa, respectively. As done in previous works, we assume that proteins are produced instantaneously after mRNAs are synthesized, so that transcription and translation processes can be lumped into a single reaction step. In addition, we assume that proteins are produced with a constant rate denoted by $\varepsilon$ (representing the mean synthesis rate) and degrade with another constant rate denoted by $d$ (representing the mean degradation rate). Without loss of generality, we set $d = 1$ throughout this paper. For clarity, we list all the reactions as follows

where $X$ represents protein. According to experimental facts ^[53], we will always assume that $B$ follows a geometric distribution given by ${\rm{prob}}\left\{ {B = i} \right\} = {{{{\left\langle B \right\rangle }^i}} / {{{\left({1 + \left\langle B \right\rangle } \right)}^{i + 1}}}}$, $i = 0, 1, 2, \cdots$, where $\left\langle B \right\rangle$ represents the mean burst size. Then, we can show that the protein expectation, denoted by $\left\langle X \right\rangle$, is given by $\left\langle X \right\rangle = \left\langle B \right\rangle \mu {\rho _{on}}$, where ${\rho _{on}} = {{{\tau _{on}}} / {\left({{\tau _{on}} + {\tau _{off}}} \right)}}$ represents the probability that the promoter is at the on state. Here ${\tau _{on}} = {1 / b}$ and ${\tau _{off}} = {1 / a}$ represent the mean times that the gene dwells at on and off states respectively. The protein noise, denoted by $\eta _X^{\left({common} \right)}$, is then given by

where ${\eta _{burst}} = {{\left({1 + \left\langle B \right\rangle } \right)} / {\left\langle B \right\rangle }}$ represents the burst noise due to the bursty expression manner, and $\eta _{promoter}^{\left({common} \right)} = {{{{\left({{\tau _{off}}} \right)}^2}} / {\left({{\tau _{on}} + {\tau _{off}} + {\tau _{on}}{\tau _{off}}} \right)}}$ represents the promoter noise due to switching between on and off states. In Eq. (2), the first term on the right hand side represents spontaneous fluctuations (i.e., the internal noise) due to the birth and death of the proteins as well as due to bursty expression, whereas the second term represents the forced fluctuations (i.e., the promoter noise), which depends, only, on the two switching rates but is independent of the synthesis and degradation rates of the protein. Note that each of three statistical quantities $\left\langle B \right\rangle$, ${\tau _{on}}$ and ${\tau _{off}}$ is experimentally measurable.

We point out that Eq. (2) is a fundamental formula for the fast evaluation of the protein noise. One will see that if other factors or processes associated with gene expression such as a multistep transition from the off to on state, then term $\eta _{promoter}^{\left({common} \right)}$ needs to be modified but the first term on the right hand side of Eq. (2) is kept unchanged.

Note that both if $b \to 0$ (or in the limit of small $b$) and if bursting is not considered, then the common on-off model reduces to the gene model of constitutive expression, which is essentially a birth-death process. In this sense, Eq. (2) includes the result for the gene-product noise in the latter model.

2.2.   A gene model of promoter progress

The gene model analyzed in the previous section assumes that the stwiching from inactive to active states is a single-step process. However, biological evidence supports that the chromatin template determining transcription kinetics accumulates over time until the promoter becomes active, implying that the off-to-on transition occurs not in a single-step manner but in a multistep manner. Specifically, the gene activity would proceed sequentially through an on state and several off states (note: in these states, polymerases could be either absent from the promoter or present in a paused or inactive state), and then returns to the on state. As a result, all the promoter’s on and off states constitute a chromatin loop. The corresponding model is called the gene model of promoter progress ^[2,3].

In order to model the gene expression with promoter progress and derive analytical results, we still assume that proteins are produced instantaneously after mRNAs are produced and degrade in a linear manner. The promoter switching from the on state to the off state has been reported to occur essentially with a single rate-limiting step and can thus be modeled by a constant rate ^[36,56]. Therefore, we introduce the following set of reactions for conveniently writing down the chemical master equation (CME) that describes the time evolution of the probability distribution function for the protein

where $b$ is the mean rate of gene inactivation, ${a_M}$ is the rate of gene activation, ${a_k}$ is the transition rate from the $k$th off state to the $\left({k + 1} \right)$th off state with $k = 1, 2, \cdots, M - 1$, $\varepsilon$ and $d$ are the synthesis and degradation rates of mRNA, respectively. Note that if $M = 1$, the corresponding model reduces to the common on-off model analyzed above. In Eq. (3), $B$ represents burst size that is assumed to follow the geometric distribution as specified above.

Denote ${\tau _k} = {1 / {{a_k}}}$, which represents the time that the gene dwells on the $k{\rm{th}}$off state. Then, ${\tau _{off}} = \sum\nolimits_{k = 1}^M {{\tau _k}}$ represents the total time that the gene dwells on the inactive state (called the mean off time). Denote by ${\tau _{on}} = {1 / b}$, which represents the mean time that the gene dwells on the on state (called the mean on time). The protein expectation still takes the form of $\left\langle X \right\rangle = \left\langle B \right\rangle \mu {\rho _{on}}$ and the protein noise ($\eta _X^{\left({loop} \right)}$) still takes the form of Eq. (2), but the promoter noise ($\eta _{promoter}^{\left({loop} \right)}$) should be modified as $\eta _{promoter}^{\left({loop} \right)} = \frac{{\left({{\tau _{on}} + {\tau _{off}}} \right)\prod\nolimits_{k = 1}^M {\left({1 + {\tau _k}} \right)} }}{{\left({1 + {\tau _{on}}} \right)\prod\nolimits_{k = 1}^M {\left({1 + {\tau _k}} \right)} - 1}} - 1$. If the total mean off time ${\tau _{off}}$ is fixed, then we can show $\prod\nolimits_{k = 1}^M {\left({1 + {\tau _k}} \right)} \geqslant {\left({{{{\tau _{off}}} / M} + 1} \right)^M}$, where the equal sign holds if and only if ${\tau _k} = {{{\tau _{off}}} / M}$ for all $k$. In this case, the promoter noise reaches the lowest level. For a fixed mean protein level, the protein noise $\eta _X^{\left({loop} \right)}$ still takes the form Eq. (2) if $\eta _{promoter}^{\left({common} \right)}$ is replaced with $\eta _{promoter}^{\left({loop} \right)}$, i.e., the protein noise consists of the term representing spontaneous fluctuations and the term representing the forced fluctuations. Moreover, the protein noise in this case equals or is greater than that in the common on-off model. Similarly, we can analyze the cases of multiple off pathways even with crosstalk, and obtain a similar result for the protein noise.

The above analysis indicates that the number ($M$) of the inactive states can enlarge the protein noise. Moreover, the portein noise in the multi-off model is always higher than that in the common on-off model, implying that the protein noise is underestimated in previous studies. In addition, model (3) is essentially equivalent to model (1) in terms of the mean on and off times.

2.3.   A gene model with feedback of any form

Feedback is a ubiquitous mechanism that regulates gene expression in prokaryotic and eukaryotic cells. Here we further introduce feedbacks to the common on-off model analyzed above. Specifically, we assume that the produced proteins regulate, as TFs, two switching rates between on and off states, as well as the protein production rate. Thus, the mean switching rates introduced above, $a$ and $b$, should be replaced with two functions depending on the number ($n$) of the protein molecules, denoted by ${\Theta _1}\left(n \right)$ and ${\Theta _2}\left(n \right)$, whereas the mean protein production rate should be replaced with another function that also depends on the number of the protein molecules, denoted by ${\Theta _3}\left(n \right)$. The corresponding gene model includes many gene models (e.g., the common on-off model analyzed above) in the literature as its special case. To that end, the corresponding CME at steady state takes the form

where ${\rm{E}}$ is the step operator whereas ${\rm{I}}$ is the unit operator, and we define ${g_i} = {\rm{prob}}\left\{ {B = i} \right\}$ and use the assumption that the protein degradation rate is the unit. In Eq. (4), ${Q_0}\left(n \right)$ and ${Q_1}\left(n \right)$ represent the stationary probability distributions of protein when the gene is in off and on states respectively. Apparently, this CME is an extension of that for the above analyzed common on-off model since each ${Q_i}\left(n \right)$ may be any function of $n$.

In principle, the total stationary distribution, $Q\left(n \right) = {Q_0}\left(n \right) + {Q_1}\left(n \right)$, can be exactly obtained by solving Eq. (4) using an analytical method that we ever developed ^[57], but the form is very complex. Here we are interested in statistical quantities of the protein, and present approximate results. For example, the stationary mean protein level is approximately given by $\left\langle X \right\rangle \approx \frac{{\left\langle B \right\rangle {\Theta _3}\left(0 \right){\Theta _1}\left(0 \right)}}{{{\Theta _1}\left(0 \right) + {\Theta _2}\left(0 \right)}}$^[57], which is apparently an extension of the expression for the mean protein in the Markovian case. Similarly, the steady-state protein noise, denoted by $\eta _X^{\left({feedback} \right)}$, is approximately given by $\eta _{promoter}^{\left({feedback} \right)} = \frac{{{\Theta _2}\left(0 \right)}}{{{\Theta _1}\left(0 \right)}}\frac{1}{{{\Theta _1}\left(0 \right) + {\Theta _2}\left(0 \right) + 1}}$^[57], which is also an extension of the expression for the promoter noise in the Markovian case.

For clarity, let us consider a special case: ${\Theta _1}\left(n \right) = a + \rho n$, where parameter $\rho$ represents feedback strength, ${\Theta _2}\left(n \right) = b$ and ${\Theta _3}\left(n \right) = \varepsilon$. In this case, we can derive the following exact expressions ^[58,59]

for the mean protein level and

for the protein noise intensity, where $\tilde a = \frac{a}{{\rho + 1}}$, $\tilde b = \frac{b}{{\rho + 1}}$, and $\tilde \varepsilon = \tilde \rho \varepsilon$ with $\tilde \rho = \frac{\rho }{{\rho + 1}}$, and ${}_1{F_1}\left({\alpha, \beta; z} \right)$ is a confluent hypergeometric function ^[60]. Note that if $\rho \to 0$, Eqs. (5a) and (5b) can reproduce previous results, e.g., ${\left\langle X \right\rangle ^{\left({common} \right)}} = \frac{{\varepsilon a}}{{a + b}}$ and $\eta _X^{\left({common} \right)} = \frac{{1 + \left\langle B \right\rangle }}{{\left\langle X \right\rangle }} + \frac{b}{{a\left({a + b + 1} \right)}}$.

3.   Non-Markov models of stochastic gene expression

In the above section, we have reviewed and analyzed Markov models of stochastic gene expression, where waiting-time distributions for reaction events are assumed to be exponential. For many systems of gene expression, however, this assumption is too strong to describe the realistic case. In fact, the waiting times in general follow non-exponential distributions or reaction kinetics is in general non-Markovian, with reasons mentioned as in the introduction. Here, we introduce and analyze several gene models of bursting expression with molecular memory, which can be divided into two classes: one class described in terms of queuing theory ^[16,36,56] and the other class described in terms of continuous time random walk (CTRW) theory ^[61,62,63].

3.1.   A queuing model

In queuing theory, queuing models can be categorized into several classes, depending on assumptions on the numbers of customers and servers, and the ways of customer arrival and departure ^[64]. Here, we map gene expression processes into queuing models in the queuing theory ^[16,36,56]. In this mapping, proteins are the analogs of customers in a queue, and their production is analogous to arrival of customers whereas their degradation to customers who leave the queue after receiving service. We assume that waiting-time distribution for protein degradation is exponential. In addition, we assume that each protein degrades independently, implying that there are an infinite number of servers in the queuing model. If we let $GI$ represent that the waiting-time distribution for customer arrival is general, $M$ represent that the distribution for customer service time is exponential, and $\infty$ represent that there are infinite servers in the queuing system, then our queuing model of gene expression belongs to the $GI/M/\infty$ system in the queuing theory.

Denote by $\varepsilon$ (representing the mean production rate) the constant rate for protein synthesis, by $d$ (representing the mean degradation rate) the constant rate for protein degradation, and by $b$ (representing the mean on-to-off switching rate) the constant switching rate from active to inactive states of the promoter. Let $\xi \left(\tau \right)$ represent the probability density function for waiting time from inactive to active states. For this queuing model, we can derive exact expressions for the steady-state moments for the number of customers in the system. In this subsection, we do not explicitly consider bursty expression.

For the above queuing model of stochastic gene expression, it seems to us that the CME has not been established so far. In order to help the reader’s understanding, here we present the following steady-state CME without providing details of mathematical derivation (note: the CME can be derived using the total probability principle)

where $h\left(\tau \right) = {{\xi \left(\tau \right)} / {\Pi \left(\tau \right)}}$ is a hazard function ^[65,66,67], and $\Pi \left(\tau \right) = \int_\tau ^\infty {\xi \left(t \right)dt}$ represents the cumulant probability. The boundary condition is imposed by ${Q_0}\left({n, 0} \right) = b{Q_1}\left(n \right)$. Equation (6) can be solved using the characteristic line method combined with the binomial moment method that we ever developed ^[59].

3.1.1.   Statistical quantities for the protein noise

Denote by $\eta _X^{\left({queue} \right)}$ the protein noise. According to the noise definition and noting that the protein degradation rate has been set as unit, we can derive the following analytical expression for the protein noise from Eq. (6)

where ${\tau _{on}} = {1 / b}$ and ${\tau _{off}} = \tilde \xi \left(0 \right)$, and $\tilde \xi \left(s \right)$ is the Laplace transform of function $\xi \left(\tau \right)$, i.e., $\tilde \xi \left(s \right) = \int_0^\infty {{e^{ - s\tau }}\xi \left(\tau \right)d\tau }$. This result is accurate in contrast to previous works that gave approximate expressions for the protein noise ^[16,36,56]. In Eq. (7), the first term on the right hand side represents spontaneous fluctuations resulting from the production and degradation of the protein whereas the second term is only related to the promoter kinetics and can then be viewed as forced fluctuations originated from switching between the promoter states. Thus, the decomposition formula in Eq. (7) is an analog of Eq. (2), namely, the protein noise is decomposed into two parts: spontaneous fluctuations resulting from the birth and death of the protein and forced fluctuations originated from switching between on and off states of the promoter.

Next, we consider a special case of distribution $\xi \left(\tau \right)$ for the off-to-on switching time, i.e., $\xi \left(\tau \right)$ is a Gamma distribution with the form of $\xi \left(\tau \right) = \frac{{{a^k}}}{{\Gamma \left(k \right)}}{\tau ^{k - 1}}{e^{ - a\tau }}$, where $k$ is a shape parameter and will be called memory index. Note that $k = 1$ corresponds to the Markov case whereas $k > 1$ to the non-Markov case. If $k$ is a positive integer, then function $\xi \left(\tau \right)$ is an Erlang distribution. Also note that the mean of $\xi \left(\tau \right)$ is $\frac{k}{a}$ and the variance is $\frac{{k\left({k + 1} \right)}}{{{a^2}}}$.

For the above setting, the mean protein level $\left\langle X \right\rangle$ is exactly given by $\left\langle X \right\rangle = \mu {\rho _{on}}$, where ${\rho _{on}}$ is the probability that the promoter is in the on state and is exactly given by ${\rho _{on}} = \frac{{{a / k}}}{{{a / k} + b}}$. This indicates that molecular memory characterized by memory index $k$ can reduce the mean protein level. In addition, the promoter noise, denoted by $\eta _{promoter}^{\left({queue} \right)}$, is given by $\eta _{promoter}^{\left({queue} \right)} = \frac{{{\tau _{off}} - \tilde \xi \left(1 \right)}}{{{\tau _{on}} + \tilde \xi \left(1 \right)}}$. If memory index $k = 1$, then we can show that $\eta _{promoter}^{\left({queue} \right)}$ reduces to $\eta _{promoter}^{\left({common} \right)} = \frac{{{\tau _{off}}{\tau _{off}}}}{{{\tau _{on}} + {\tau _{off}} + {\tau _{on}}{\tau _{off}}}}$. If $k > 1$, we can prove that function $\kappa \left(I \right) \equiv \frac{{{\tau _{off}} - \tilde \xi \left(1 \right)}}{{{\tau _{on}} + \tilde \xi \left(1 \right)}}$ with ${\tau _{off}} = \frac{k}{a}$ and $\tilde \xi \left(1 \right) = 1 - {\left({\frac{a}{{1 + a}}} \right)^k}$ is a monotonically increasing function of index $k$. Thus, $\kappa \left(k \right) \geqslant \kappa \left(1 \right) = \frac{b}{{a\left({a + b + 1} \right)}} = \eta _{promoter}^{\left({common} \right)}$, indicating that the promoter noise in the queuing model is stronger than that in the common on-off model. This monotonicity implies that molecular memory always amplifies the protein noise if the protein mean is fixed.

3.1.2.   Bounds for the number of promoter states

The experimental observations ^[18,33,54] have indicated that the waiting time for promoter switching from inactive to active states in general follows an non-exponential distribution, but the specific form of this distribution is in general unknown since it depends, in general, on the number of promoter activity states, possibly changing from gene to gene. On the other hand, the current single-cell measurement technologies make it possible to count the number of transcripts in different cells. An interesting yet unsolved question is the estimation on the number of promoter internal states. In order to address this question, we will analyze the Fano factor of protein rather than its noise intensity as done above, focusing on deriving bounds on the number of promoter internal states based on experimentally accessible measurements.

For analysis convenience, we only consider a special case: the number of inactive states is fixed at a positive integer ($k$). In this case, the distribution for the off-to-on switching time, denoted by $\xi \left(\tau \right)$, can be modeled as a phase-type distribution of order $k$ ^[68]. Note that for a fixed mean switching time, the noise for the phase-type distribution of order $k$ is the minimum if the distribution is set as $\xi \left(\tau \right) = \frac{{{a^k}}}{{\Gamma \left(k \right)}}{\tau ^{k - 1}}{e^{ - a\tau }}$ (i.e., an Erlang distribution with shape parameter $k$ ^[68,69]), where parameter $a$ represents the rate for a single-step reaction whereas parameter $k$ represents the number of reaction steps. This setting will allow us to derive bounds on the Fano factor for protein with details being as follows.

Note that for the $\xi \left(\tau \right)$ set above, we can show that the exact expression for the Fano-factor, denoted by $FF\left(k \right)$, is given by $FF\left(k \right) = 1 + \left\langle B \right\rangle {\left({{\rho _0}} \right)^2}\left({1 + {1 / k}} \right)$, where ${\rho _0}$ represents the probability that the promoter is in an off state and is a fixed number if the protein expectation is fixed. This expression indicates that the Fano factor of protein is a monotonically decreasing function of $k$, and reaches the maximum at $k = 1$. Therefore, we have inequality $FF\left(\infty \right) \leqslant FF\left(k \right) \leqslant FF\left(1 \right)$. Notice that $2FF\left(\infty \right) = 1 + FF\left(1 \right)$ holds. Thus, both the lower and upper bounds can be specified in terms of $FF\left(1 \right)$, that is, ${{\left({1 + FF\left(1 \right)} \right)} / 2} \leqslant FF\left(k \right) \leqslant FF\left(1 \right)$. Furthermore, if we make the transform of ${\vartheta _{theor}} = {{\left[{FF\left(1 \right) - 1} \right]} / {\left[{\left\langle B \right\rangle \rho _0^2} \right]}}$, which can be rewritten as ${\vartheta _{theor}} = 2 + \left({{1 / k}} \right)$, then we have $1 \leqslant {\vartheta _{theor}} \leqslant 2$, which is independent of memory index $k$.

The significance of the above analysis is as follows. Consider an experimental analog of the $FF$ defined above, which has the experimental Fano factor ${\vartheta _{{\rm{exper}}}}$. In order to set upper and lower bounds on the number of promoter states, it is needed to determine the largest value of integer $k$ such that ${\vartheta _{{\rm{exper}}}} < {\vartheta _{theor}}$. Note that the Erlang distribution of order $k$ has the lowest variance among all phase-type distributions of order $k$. Also note that ${\vartheta _{theor}}$ is a monotonically decreasing function of $k$. The combination of both implies that the promoter must have at least $k + 1$ off states if ${\vartheta _{{\rm{exper}}}} < {\vartheta _{theor}}$ is satisfied. Thus, measurements of the Fano factor for protein can be used in setting bounds on the minimum number of promoter states.

3.2.   A CTRW model

As is well known, CTRWs can be further divided into two classes ^[62]: active CTRWs where waiting times need to be reset; passive CTRWs where waiting times are not reinitialized. The above analyzed queuing model of gene expression belongs to the latter class. From the viewpoint of stochastic simulation, however, it is more convenient to use active CTRWs to describe stochastic processes including gene expression. In particular, physicians frequently use the active CTRW framework to analyze and simulate stochastic processes ^[61,62,63]. Here we consider a stochastic gene expression model in the sense of active CTRWs.

Let ${\psi _1}\left({t; n} \right)$ and ${\psi _2}\left({t; n} \right)$ be intrinsic-event waiting-time distributions from off to on states and vice versa respectively, ${\psi _3}\left({t; n} \right)$ and ${\psi _4}\left({t; n} \right)$ be intrinsic-event waiting-time distributions for protein synthesis and degradation respectively, where $n$ represents the number of protein molecules. Note that self-regulation exists if ${\psi _1}\left({t; n} \right)$ depends on $n$ and there is no feedback otherwise. The similar case holds for ${\psi _2}\left({t; n} \right)$. If ${\psi _3}\left({t; n} \right)$ depends on $n$, this implies post-transcriptional or posttranslational regulation. Assume that proteins are generated in bursts with burst size ($B$) following the geometric distribution described above. To that end, we have finished the setting of our gene model. For convenience, we denote by ${R_i}$ ($1 \leqslant i \leqslant 4$) four reactions from off to on states and vice versa, synthesis and degradation of protein, respectively.

3.2.1.   Effective transition rates

Since ${\psi _i}\left({t; n} \right)$ ($i = 1, 2, 3, 4$) may be general distributions, which may lead to molecular memory, each of 4 reactions ${R_i}$ ($1 \leqslant i \leqslant 4$) has a memory function denoted by ${M_i}\left({t; n} \right)$ ^[62,73,74]. Interestingly, we can prove that for arbitrary waiting-time distributions ${\psi _k}\left({t; n} \right)$ for reaction ${R_k}$ ($1 \leqslant k \leqslant 4$), the limit $\mathop {\lim }\limits_{s \to 0} {\tilde M_i}\left({s; n} \right)$ ($1 \leqslant i \leqslant 4$) always exists ^[70]. If this limit is denoted by ${\Theta _i}\left(n \right)$, then we have ^[70]

where ${\bar \psi _k}\left({t; n} \right)$ is a modification of ${\psi _k}\left({t; n} \right)$ due to the consideration of on and off states. Specifically, ${\bar \psi _k}\left({t; n} \right) = {\psi _k}\left({t; n} \right)$ ($k = 1, 4$) and ${\bar \psi _k}\left({t; n} \right) = 0$ ($k = 2, 3$) when we calculate ${\Theta _1}\left(n \right)$ according to ${\Theta _1}\left(n \right) = \mathop {\lim }\limits_{s \to 0} {\tilde M_1}\left({s; n} \right)$; ${\bar \psi _1}\left({t; n} \right) = 0$, ${\bar \psi _k}\left({t; n} \right) = {\psi _k}\left({t; n} \right)$ ($k = 2, 3, 4$) when we calculate ${\Theta _2}\left(n \right)$ or ${\Theta _3}\left(n \right)$ according to ${\Theta _k}\left(n \right) = \mathop {\lim }\limits_{s \to 0} {\tilde M_k}\left({s; n} \right)$($i = 2, 3$); and ${\bar \psi _i}\left({t; n} \right) = {\psi _i}\left({t; n} \right)$ ($1 \leqslant i \leqslant 4$) when we calculate ${\Theta _2}\left(n \right)$ according to ${\Theta _4}\left(n \right) = \mathop {\lim }\limits_{s \to 0} {\tilde M_4}\left({s; n} \right)$.

Function ${\Theta _i}\left(n \right)$ is called the effective transition rate for the ${\rm{i}}th$ reaction ^[70,71], where $1 \leqslant i \leqslant 4$. Note that if ${\psi _i}\left({t; n} \right)$ is an exponential distribution of the form ${\psi _i}\left({t; n} \right) = {\lambda _i}\left(n \right){e^{ - {\lambda _i}\left(n \right)t}}$, where ${\lambda _i}\left(n \right)$ should be understood as the reaction propensity function for ${R_i}$, we can show ${\Theta _i}\left(n \right) = {\lambda _i}\left(n \right)$. This indicates that effective transition rates are extensions of reaction propensity function s. If all waiting-time distributions are exponential, which corresponds to the Markov reaction case, the model reduces to the common CME for the on-off model with feedback and bursting in the sense of Laplace transform. We point out that the introduction of effective transition rates will be a key for analyzing non-Markov behavior. In addition, we emphasize that our gene model includes almost the on-off models of gene expression in the existing literature as its special cases.

The promoter switching from the on state to the off state has been reported to occur essentially with a single rate-limiting step and can thus be modeled by a constant rate ^[39]. Therefore, ${\psi _2}\left({t; n} \right)$ can be assumed to be an exponential distribution, i.e., ${\psi _2}\left({t; n} \right) = b{e^{ - bt}}$, where $b$ represents the mean switching rate from on to off states. Assume that proteins degrade in a linear manner with a constant rate denoted by $d$ (implying that the waiting-time distribution for degradation takes the form ${\psi _4}\left({t; n} \right) = n{e^{ - nt}}$ due to the assumption of $d = 1$).

In order to show the explicit effect of molecular memory, we consider the following two special cases.

Case 1: ${\psi _1}\left({t; n} \right) = \left[{{{{a^{{k_1}}}} / {\Gamma \left({{k_1}} \right)}}} \right]{t^{{k_1} - 1}}{e^{ - at}}$, ${\psi _2}\left({t; n} \right) = b{e^{ - bt}}$ and ${\psi _3}\left({t; n} \right) = \varepsilon {e^{ - \varepsilon t}}$, where positive constants $a$, $b$ and $\varepsilon$ represent an average switching rates from off to on and vice versa, and an average transcriptional or translational rate, respectively. Note that ${k_1} = 1$ corresponds to the Markov reaction case whereas ${k_1} > 1$ to the non-Markov reaction case. Therefore, ${k_1}$ is also called memory index. According to Eq. (8), we can show ${\Theta _1}\left(n \right) = {{n{a^{{k_1}}}} / {\left[{{{\left({a + n} \right)}^{{k_1}}} - {a^{{k_1}}}} \right]}}$ with ${\Theta _1}\left(0 \right) = {{{\lambda _1}\left(0 \right)} / {{k_1}}}$, ${\Theta _2}\left(n \right) = b$, ${\Theta _3}\left(n \right) = \varepsilon$, and ${\Theta _4}\left(n \right) = n$. Therefore, the effect of molecular memory is equivalent to the introduction of a negative feedback if ${I_1} > 1$.

Case 2: ${\psi _1}\left({t; n} \right) = a{e^{ - at}}$, ${\psi _2}\left({t; n} \right) = b{\rm{ }}{e^{ - b{\rm{ }}t}}$, and ${\psi _3}\left({t; n} \right) = \left[{{{{\varepsilon ^{{k_3}}}} / {\Gamma \left({{k_3}} \right)}}} \right]{t^{{k_3} - 1}}{e^{ - \varepsilon t}}$, where ${k_3}$ is also called a memory index. Note that ${k_3} > 1$ corresponds to the non-Markov case and ${k_3} = 1$ corresponds to the Markov case. In this case, we can show ${\Theta _3}\left(n \right) = {\varepsilon ^{{k_3}}}\frac{{b + n}}{{{{\left({b + \varepsilon + n} \right)}^{{k_3}}} - {\varepsilon ^{{k_3}}}}}$, ${\Theta _1}\left(n \right) = a$, ${\Theta _2}\left(n \right) = b$, and ${\Theta _4}\left(n \right) = n$. Moreover, the effect of molecular memory is equivalent to the introduction of a negative posttranscriptional or posttranslational regulation if ${k_3} > 1$.

Similarly, we can analyze the combination of the above two cases. In a word, effective transition rates not only explicitly decode the effect of molecular memory but also can give us useful information on the underlying system. In particular, the use of effective transition rates can transform a non-Markov issue into a Markov one ^[70]. Therefore, we can use the results in section 2.3 to directly give results for the non-Markov gene model.

3.2.2.   Statistical quantities for the protein noise

Based on Eq. (4) above, we can first show that the stationary mean protein level $\left\langle X \right\rangle$ is approximately given by $\left\langle X \right\rangle \approx \left\langle B \right\rangle \frac{{{\Theta _1}\left(0 \right){\Theta _3}\left(0 \right)}}{{{\Theta _1}\left(0 \right) + {\Theta _2}\left(0 \right)}}$, which is apparently an extension of the expression for the mean protein in the Markov case. Then, we can show that the noise intensity for the protein at steady state, denoted by $\eta _X^{\left({CTRW} \right)}$, still takes the form of Eq. (2) if the expression of $\eta _{promoter}^{\left({on - off} \right)}$ is replaced with $\eta _{promoter}^{\left({CTRW} \right)} = \frac{{{\Theta _2}\left(0 \right)}}{{{\Theta _1}\left(0 \right)}}\frac{1}{{{\Theta _1}\left(0 \right) + {\Theta _2}\left(0 \right) + 1}}$, which is also an extension of the expression for the promoter noise in the Markov case.

In order to analytically show the effects of global delay and molecular memory, we consider the above two special cases.

For Case 1, the mean level and the noise intensity for the protein are approximated as $\left\langle X \right\rangle \approx \left\langle B \right\rangle \frac{{\varepsilon a}}{{a + b{k_1}}}$ and $\eta _X^{\left({CTRW} \right)} \approx \frac{{1 + \left\langle B \right\rangle }}{{\left\langle X \right\rangle }} + \frac{b}{a}\frac{{k_1^2}}{{a + \left({b + 1} \right){k_1}}}$, respectively. We observe that a larger value of memory index ${k_1}$ leads to a lower mean protein level but to the stronger protein noise (including the stronger promoter noise) if the mean protein level is fixed.

For Case 2, the mean level and the noise intensity for the protein can be approximated as $\left\langle X \right\rangle \approx \left\langle B \right\rangle \frac{a}{{a + b}}\frac{{a{\varepsilon ^{{k_3}}}}}{{{{\left({a + \varepsilon } \right)}^{{k_3}}} - {\varepsilon ^{{k_3}}}}}$ and $\eta _X^{\left({CTRW} \right)} \approx \frac{{1 + \left\langle B \right\rangle }}{{\left\langle X \right\rangle }} + \frac{b}{a}\frac{1}{{a + b + 1}}$, respectively. Apparently, $\left\langle X \right\rangle$ is a monotonically decreasing function of memory index ${k_3}$.

The above analysis indicates that molecular memory is an unneglectable factor affecting the gene expression level and noise.

4.   Conclusion and discussion

Although gene expression is complex biochemical process, the noise in mRNA or protein can be decomposed into two parts: spontaneous fluctuations resulting from the birth and death of the protein and the forced fluctuations originated from switching between the promoter states, as shown in Eq. (2). This formula is universal in form, suitable to various gene-expression models established based on the central dogma in Biology, which can be categorized into two classes: Markov and non-Markov models. But specific processes may fine-tune the level of the protein noise.

4.1.   On upstream and downstream dynamics of gene expression

The models analyzed in this paper were simplified, namely they only considered downstream dynamics and neglected upstream dynamics. However, chromatin regulators play a major role in establishing and maintaining gene expression states. For these regulators, silencing and reactivation occur in all-or-none events, enabling the regulators to modulate the fraction of cells silenced rather than the amount of gene expression. The regulators operate over different time scales and generate distinct types of epigenetic memory, through their individual transition rates. Relevant dynamics can be described by a three-state model involving stochastic transitions between active, reversibly silent, and irreversibly silent states ^[25].

All-or-none stochastic switching of gene expression states can occur at two levels: chromatin-mediated switching and transcriptional bursting. Active and silent chromatin states may transition to each other. However, even within an active chromatin state, promoters in general switch stochastically between transcriptionally active and inactive states (this phenomenon is called transcriptional bursting ^[72]). In ref. ^[25], the authors discussed the connection between these two modes of gene regulation. They also analyzed the dynamic regimes that lead to either graded or fractional all-or-none responses at the protein level. Interestingly, they found that chromatin-mediated switching generally produces bimodal protein distributions similar to those observed experimentally, but can also produce graded protein level distributions in some parameter regimes.

The two levels of gene regulation can be combined in a single model (Figure 2) ^[72]. In this model, a gene can switch between active and silent chromatin states. In the silent chromatin state, the promoter is always in an off state, so the mRNA level is zero. By contrast, the active chromatin state is permissive, allowing the promoter to switch between periods of active transcription, which produce bursts of multiple mRNA molecules, and periods of inactivity ^[4,73,74].

In general, transitions in chromatin state and transcriptional bursting are associated with different timescales. At the transcription level, switching between on and off promoter states has been attributed to short-lived interactions of TFs and core machinery with the promoter, occurring usually on the timescale of seconds to minutes ^[72]. By contrast, at the chromatin level, switching between active and silent chromatin states occur usually on the timescale of hours to days, an order of magnitude much slower than the timescales involved in transcriptional bursting.

The timescales of switching can determine the noise in protein levels at the population level. When switching timescales are fast in contrast to mRNA and protein half-lives, which is typically the case for transcription factor regulation, the level of protein expressed from the gene reaches a stationary unimodal distribution ^[75]. The mean value of this distribution depends, in a graded manner, on the occupancy of the TF at the promoter ^[76,77]. When switching timescales are slow relative to mRNA and protein half-lives, as is generally the case for chromatin-mediated regulation, one expects bimodal protein distributions. The relative ratio of the two peaks in this bimodal distribution depends not only on the occupancy of the chromatin remodeling at the promoter, but also on the time it has been recruited there.

4.2.   On division of noisy sources

Previous studies have used different classifications to clarify different sources of gene expression noise. Here, we focus on the distinction between intrinsic and extrinsic noise by considering four aspects: The statistical nature of fluctuations, correlations between different proteins in a single cell, how central a process is to gene expression, and experimental strategies for measuring noise:

(Ⅰ) Terms ‘intrinsic’ and ‘extrinsic’ have no specific meaning except that they mean ‘inside’ and ‘outside’ respectively. In fact, the classification of intrinsic and extrinsic noise depends on the definition of system versus environment and is therefore relative. For example, if proteins are taken as system, the first term in Eq. (2) is intrinsic and the promoter noise is extrinsic ^[53,78]. Some studies ^[16,31] define the two terms on the right hand side of Eq. (2) as the intrinsic noise to distinguish them from the extrinsic noise in the overall state of the cell. However, spontaneous fluctuations in the protein level are largely different from the noise that originates in enslavement by the promoter, but the latter noisy sources are not fundamentally different, e.g., ribosome-mediated noise also enslaves proteins. If we consider the statistical nature of fluctuations only, there is thus no reason to label the two componential noise terms in Eq. (2) as intrinsic.

(Ⅱ) Protein noise can be classified based on correlations between different types of proteins. Some noisy sources are shared by multiple genes in a single cell whereas others are exclusive to a particular gene or a small set of genes. The most specific componential noise comes from having low protein numbers, originating in the random births and deaths of individual molecules. additionally, there are many other specific sources of gene expression noise, which include: (a) spontaneous mRNA fluctuations are quite specific although some transcripts encode several different proteins; (b) in contrast, operator fluctuations are less specific and typically affect all genes in an operon; and (c) many DNA and RNA binding proteins act as repressors and activators, some of which are specific to a particular gene and others would regulate large classes of genes. This can form complicated correlation structures. For example, two genes may be regulated by the same repressor, but their mRNAs may be degraded by different RNases. A few central factors, e.g., ribosomes, core polymerases, tRNAs, and amino acids, are shared universally. However, different proteins are still affected differently by such fluctuations because they have different sequences and lifetimes.

From this perspective, the two terms in Eq. (2) or its other forms could all be intrinsic. However, it should be emphasized that such a classification can be made in many different ways. The mRNA term would be more intrinsic than the promoter term since the former is directly related to translation, but the latter is not necessarily more intrinsic than fluctuations in regulator concentrations.

(Ⅲ) From a biological viewpoint, sources of gene expression noise could be categorized according to how central the corresponding component is to gene expression. But classifying fluctuations in promoter activity as intrinsic and fluctuations in ribosomes as extrinsic would not separate the central parts of gene expression from more peripheral cell processes. The opposite classification seems more appropriate, e.g., ribosomes are inherent to gene expression whereas spontaneous changes in promoter activity can indirectly reflect regulation.

(Ⅳ) The classification of noisy sources can be tailor-made to the experimental methods available. A good strategy for separating noisy sources would be based on correlations between the expressions of two physically separate but identically regulated fluorescent reporter genes ^[74,79]. If the two proteins are kinetically independent, e.g., if the underlying mechanisms are linear, then the normalized covariance between them is equal to the sum of all the common noise terms ^[14]. Fluctuations in ribosomes, polymerases, RNases, etc., together end up in the ‘extrinsic’ covariance category while the two componential noise terms in Eq. (2) or its other forms are ‘intrinsic’ because each green fluorescent protein has its own operators and transcripts. The appeal of this approach is not only that some noisy sources are separated from others but also that the separation to some extent relates to specificity (see subsection 4.2.2 above). The only risk for this separation is that the terms are over-explained, comprehending biological or physical meaning where none is intended: The ‘intrinsic noise’ only partially relates to specificity, as demonstrated by experiments where most extrinsic noise comes from a repressor that is specific to that particular gene ^[74]. Realizing that the distinction is largely a side effect of the experimental set-up can in turn open the doors for other applications. If the two reporters were placed under the same operator, then operator fluctuations would become shared between the corresponding two genes, thus moving from the intrinsic to the extrinsic category. But if the reporters were encoded on the same transcript, the mRNA term would become extrinsic. In addition, if the reporters were regulated by different repressors, repressor fluctuations would move from the extrinsic to the intrinsic category. The dual reporter strategy does not separate fluctuations in target proteins based on a priori physical or biological principles. It is much more general and useful than that.

4.3.   On experimental measures of gene expression noise

There are many ways to quantify gene expression noise, e.g., autocorrelation analysis is an analytically tractable way since autocorrelations conveniently summarize both the magnitude and the frequency of fluctuations. However, most models of stochastic gene expression have focused on stationary averages and variances so far.

The above analytical results for the protein noise are formulated in terms of the ratio of the variance over squared average, which allows for the fast evaluation and a clear separation of different noise sources as long as the corresponding gene models are weakly nonlinear. Another common measure is the Fano factor, which is defined as the ratio of the variance over the expectation and has a good property: it equals one for Poisson distributions. The Fano factor can work well only for univariate discrete random processes where the variance is proportional to the average with a proportional constant that reflects the overall nature of the process. For a multivariate random process, however, the Poisson distribution holds no special position and using the Fano factor would be misleading. To illustrate this with a more extreme example than the bursts in Eq. (4), we assume that fluctuations in the protein abundance come, partly, from fluctuations in ribonuclease or protease concentrations. The ratio of the variance over the squared average would then contain an extra term that is more or less independent of transcription or translation rates. Multiplying by the average to obtain the Fano factor would thus force the measure to depend on anything that impacts the average.

The results given by Eq. (5a, 5b) or obtained by Eq. (8) suggest that no measures work well for all types of fluctuations: the componential noise due to non-exponential waiting-time distribution cannot be separated from the promoter noise; spontaneous fluctuations depend on the number of molecules whereas forced fluctuations do not. In spite of this, the ratio of the variance over the squared mean is still a suitable basis for various possible experimental interpretations. This is because (a) the average number of proteins per cell is too high to contribute substantial spontaneous fluctuations in most experimental studies so far; (b) by plotting the ratio of the variance over the squared average as a function of the inverse average, any univariate scaling behavior can easily be identified without introducing scaling problems for any extrinsic source of noise; and (c) the relevance of fluctuations depends, typically, on the size of the underlying system. The variance (i.e., a second-order moment) must be normalized by the squared average.

Acknowledgments

This work was supported by grants 11775314 and 11931019 from Natural Science Foundation of China, and 202007030004 from Key-Area Research and Development Program of Guangzhou.

Conflict of interest

All authors declare no conflicts of interest in this paper.