Kinetic characteristics of transcriptional bursting in a complex gene model with cyclic promoter structure

1.   Introduction

Transcription is a core step of gene-expression processes. For most genes in eukaryotic cells, transcription often occurs in a bursty fashion characterized by mRNA synthesis in short periods of activity followed by typically longer periods of inactivity. Single-cell experimental observations have provided evidence for such transcription bursting ^[1,2,3,4]. This variability, which potentially gives rise to cell-to-cell variability in gene expression levels ^[5,6,7], can propagate to downstream proteins or target genes through a network. It is believed that the expression noise can not only contribute to alternative cellular fates ^[8] but also maintain cellular functioning ^[9,10,11]. Revealing the mechanism of transcriptional bursting using mathematical models is an important step toward understanding fundamental cellular processes.

Transcriptional bursting kinetics is often characterized by two parameters: burst size and burst timing (or burst frequency). The former is defined as the number of transcripts produced during an active period, whereas the latter involves waiting times between transcription initiation events. For rate-limiting transcription or for transcription with exponential waiting times, mRNA burst sizes follow geometric distributions, which were also observed in early experiments ^[1,3,10,12]. However, transcription is a complex process involving many regulatory proteins and molecular factors, such as transcription factors, RNA polymerase, and chromatin remodeling complexes ^{[13,14,15,16]}. Recent experiments have shown that these factors or processes can result in non-exponential waiting-time distribution for genes' active (ON) and inactive (OFF) durations ^{[17,18,19,20]}. In addition, some experimental studies have also suggested complex promoter structures ^[21,22,23], which often involve nucleosomes competing with or being removed by transcription factors ^[24], and epigenetic regulation via histone modifications ^[25,26,27]. For example, Dunham et al. ^[28] found that the transcriptional deactivation occurs with graded, stepwise decreases in transcription rate, implying that the gene promoter has variable transcription states. Sepúlveda et al. ^[29] demonstrated that the lysogeny maintenance promoter of phage lambda switches between multiple promoter states. Neuert et al. ^[30] showed that the transcription is a multistate process for the osmotic stress response pathway in Saccharomyces cerevisiae. In a word, transcription is finished not in a simple single-step manner but often in a complex multistep fashion ^{[31,32,33,34]}.

It has been verified that mathematical models are a powerful tool for estimating the impacts of molecular mechanisms on transcriptional bursting. In past decades, many theoretical models have been developed to characterize transcriptional bursting kinetics including burst size, burst frequency, and variability in the mRNA abundance ^[35,36]. Of these models, the most widely used model is the two-state model ^{[17,37,38,39,40]}, where the promoter is assumed to switch between transcriptionally active (ON) and inactive (OFF) states with constant switching rates. Queuing models have also been proposed to analyze the impact of transcription regulation on bursting kinetics ^{[31,41,42,43]}, where waiting-time distributions are general (either exponential or non-exponential). For example, Schwabe et al. used Erlang-distributed times to model ON-phase and transcription ignition, and found that burst sizes peaked distribution ^[44]. Kumar et al. derived steady-state moments of mRNA distribution and estimated burst parameters for non-geometric burst distribution ^[45]. Multistate gene models such as multi-OFF models ^[46,47], and multi-ON models ^[48,49], have been proposed to analyze transcription dynamics. Recently, Zoller et al. ^[50] combined mathematical modeling with experimental measurement to estimate burst parameters in a multi-OFF promoter cycle. In contrast, Rybakova et al. ^[51] used stochastic and ordinary differential equation models to estimate burst parameters in a multi-ON promoter cycle. In addition, Daigle et al. ^[52] inferred the number and configuration of promoter topologies from single-cell gene expression data. In spite of these efforts, how promoter kinetics affect transcription bursting kinetics remains poorly understood. Neither gene models with multi-ON states or with cyclic promoter structure, whose prototypes can be found in natural and synthetic systems ^[53], have been systematically studied, nor analytical results on burst size and timing have been reported.

In this paper, we introduce and analyze a model of stochastic transcription with cyclic promoter structure (i.e., ON and OFF states constitute a cycle), focusing on analysis of the effect of multi-ON mechanism on transcriptional bursting kinetics. Interesting results are obtained, e.g., burst size follows a mixed geometric distribution rather than a single geometric distribution that was assumed often in previous studies; each of ON and OFF times obeys a mixed exponential distribution; and the multi-ON mechanism can lead to bimodal burst-size distribution, antagonistic timing of ON and OFF. While our results reveal essential characteristics of transcriptional bursting kinetics, our analysis provides a framework for studying more complex transcription processes.

2.   Methods

2.1.   Model hypothesis

Here we introduce a biologically reasonable model for transcription bursting (Figure 1(A)), which considers the effects of many processes (or factors) such as recruitment of various polymerases and transcription factors, assembly of the pre-initiation complexes, chromatin remodeling, and histone modifications ^[13,14,15]. Then, we map this model into a theoretical model of multistate promoter cycle (Figure 1(B)) where the promoter proceeds sequentially through several reversible active (ON) states and inactive (OFF) states and returns to the active (ON) state, forming a cycle. We assume that transcription occurs in a bursting manner at each active state but is prohibited at the inactive states (Figure 1(C)). Specifically, we assume that the gene promoter has $N$ states in total, which consists of $L$ active states and $K$ inactive states. We call $L$ and $K$ as active and inactive indexes, respectively. Denoted by ${A_1}, \cdot \cdot \cdot, {A_L}$ ON states and by ${I_1}, \cdot \cdot \cdot, {I_K}$ OFF states. The corresponding biochemical reactions are listed in Table 1.

Note that both $L = 1$ and $N > 2$ correspond to the multi-OFF model (i.e., the promoter has more than one OFF state but only one ON state), and both $K = 1$ and $N > 2$ correspond to the multi-ON model (i.e., the promoter has more than one ON state but only one OFF state). Therefore, the model introduced here includes transcription models that were previously studied ^{[46,47,48,49]}, and can be used to model general transcription dynamics generated by gene promoters with multiple states. We point out that for the multi-OFF and multi-ON models, some analytical results on steady-state mRNA distribution and noise have been obtained ^[46,48], but for multistate promoter models, the analytical expressions on transcription bursting kinetics, such as distributions of burst size and timing, have not been derived so far.

2.2.   Analysis framework

2.2.1.   Burst-size distribution

Let $U$ represent the promoter's active state, which is a time-dependent random variable, i.e., $U\left(t \right) = u \in \left\{ {{A_1}, \cdots, {A_L}} \right\}$, $B\left(t \right)$ be burst size, i.e., $B\left(t \right) = n, {\text{ }}n = 0, 1, 2, \cdot \cdot \cdot$, and $T$ represent the exit time (i.e., burst termination time) from an active state. We assume that at time $t = 0$, the promoter begins to switch to active state ${A_1}$ from inactive state ${I_K}$, or to active state ${A_L}$ from inactive state ${I_1}$. Let ${S_{n, u}}\left(t \right)$ represents the survival probability that $n$ mRNAs are produced per burst in state $u$ at time $t$. Then, the transcriptional system's state at time $t$ can be described as ${S_{n, u}}\left(t \right) = \Pr {\text{ob}}\left\{ {T > t\left| {B(t) = n, U\left(t \right) = u} \right.} \right\}$, where we stipulate ${S_{n, u}}\left(t \right) \equiv 0\; \left({n < 0} \right)$. Note that the probability that burst termination time $T$ falls within an infinitesimal interval $\left({t, t + \Delta t} \right)$ can be expressed as

Therefore, the joint probability density function ${P_{{\text{burst}}}}\left({n, t} \right)$ for burst size $B\left(t \right)$ (a discrete random variable) and burst termination time $T$ (a continuous random variable) includes all incoming fluxes driving the system from active state ${A_1}$ to inactive state ${I_K}$, or from active state ${A_L}$ to inactive state ${I_1}$, that is, ${P_{{\text{burst}}}}\left({n, t} \right) = {\lambda _{1, N}}{S_{n, {A_1}}}\left(t \right) + {\lambda _{L, L + 1}}{S_{n, {A_L}}}\left(t \right).$ Then, the steady-state probability distribution for burst size $B$ is calculated according to

The key question is how ${S_{n, {A_1}}}\left(t \right)$ and ${S_{n, {A_L}}}\left(t \right)$ are obtained. However, if we rewrite ${S_{n, {A_l}}}\left(t \right)$ as ${S_{n, l}}\left(t \right)$ for convenience, the master equation for survival probability ${S_{n, l}}\left(t \right)$ can be described as

where ${\delta _{i, j}}$ is a Kronecker delta. Note that this equation can be solved by introducing probability-generating functions. Specifically, defining ${G_l}\left({z, t} \right) = \sum\limits_{n = 0}^\infty {{S_{n, l}}\left(t \right){z^n}}$ ($1 \leqslant l \leqslant L$) and denoting ${\bf{G}}\left({z, t} \right) = {\left({{G_1}\left({z, t} \right), \cdots, {G_L}\left({z, t} \right)} \right)^T}$, we can derive a linear partial differential equation group of the form: ${{\partial {\bf{G}}\left({z, t} \right)} \mathord{\left/ {\vphantom {{\partial {\bf{G}}\left({z, t} \right)} {\partial t}}} \right. } {\partial t}} = {\bf{H}}\left(z \right){\bf{G}}\left({z, t} \right)$ ^[54,55], where the coefficient matrix ${\bf{H}}(z)$ is a tridiagonal matrix with the diagonal elements taking the form ${h_l}\left(z \right) = \mu \left({z - 1} \right) - {\delta _{l, 1}}{\lambda _{l, N}} - \left({1 - {\delta _{l, 1}}} \right){\lambda _{l, l - 1}} - {\lambda _{l{\kern 1pt}, l + 1}}$ ($1 \leqslant l \leqslant L$) (seeing Appendix A for details).

2.2.2.   Dwell-time distributions and burst frequency

Note that the distribution of the time that the promoter dwells on active (or inactive) state, denoted by ${f^{{\text{ON}}}}\left(t \right)$ (or ${f^{{\text{OFF}}}}\left(t \right)$), is equal to the sum of the times dwelling on all ON (or OFF) states. In order to calculate ${f^{{\text{ON}}}}\left(t \right)$ and ${f^{{\text{OFF}}}}\left(t \right),$ we introduce two column vectors ${{\bf{S}}^{{\text{ON}}}}\left(t \right) = {\left({S_1^{{\text{ON}}}\left(t \right), \cdots, S_L^{{\text{ON}}}\left(t \right)} \right)^T}$ and ${{\bf{S}}^{{\text{OFF}}}}\left(t \right) = {\left({S_1^{{\text{OFF}}}\left(t \right), \cdots, S_K^{{\text{OFF}}}\left(t \right)} \right)^T}$, where $S_l^{{\text{ON}}}\left(t \right)$ and $S_k^{{\text{OFF}}}(t)\, \,$ are the survival probabilities that the promoter is still at ${l^t}^h$ ON state and ${k^t}^h$ OFF state, respectively.

Now, we establish the master equation for $S_l^{{\text{ON}}}\left(t \right)$ and $S_k^{{\text{OFF}}}(t)\, \,$. Denote by ${{\bf{T}}^{{\text{ON}}}}$ ($L \times L$ matrix) and ${{\bf{T}}^{{\text{OFF}}}}$ ($K \times K$ matrix) the transition matrices of promoter among the active states and among the inactive states, respectively. Then, the master equations of the survival probability for ON and OFF states can take the following unified form

where $\vartheta = {\text{ON}}$ or ${\text{OFF}}$.

Equation (4) is a linear equation group and is easily solved. In order to derive the analytical expression of ${{\bf{S}}^\vartheta }(t),$ we first introduce two matrices ${{\bf{\Lambda }}^{{\text{ON}}}} = {\text{diag}}\left({ - \lambda _1^{{\text{ON}}}, \cdots, - \lambda _L^{{\text{ON}}}} \right)$ and ${{\bf{\Lambda }}^{{\text{OFF}}}} = {\text{diag}}\left({ - \lambda _1^{{\text{OFF}}}, \cdots, - \lambda _K^{{\text{OFF}}}} \right)$, where $- \lambda _l^{{\text{ON}}}$ and $- \lambda _k^{{\text{OFF}}}$ are eigenvalues of matrices ${{\bf{T}}^{{\text{ON}}}}$ and ${{\bf{T}}^{{\text{OFF}}}},$ respectively. Then, we introduce two diagonal matrices ${{\bf{R}}^{{\text{ON}}}} = {\text{diag}}\left({\gamma _1^{{\text{ON}}}, \cdot \cdot \cdot, \gamma _L^{{\text{ON}}}} \right),$ where $\gamma _1^{{\text{ON}}} = 1$ and $\gamma _l^{{\text{ON}}} = \prod\limits_{i = 2}^l {\sqrt {{{{\lambda _{i, i - 1}}} \mathord{\left/ {\vphantom {{{\lambda _{i, i - 1}}} {{\lambda _{i - 1, i}}}}} \right. } {{\lambda _{i - 1, i}}}}} }$ ($2 \leqslant l \leqslant L$), and ${{\bf{R}}^{{\text{OFF}}}} = {\text{diag}}\left({\gamma _1^{{\text{OFF}}}, \cdot \cdot \cdot, \gamma _K^{{\text{OFF}}}} \right),$ where $\gamma _1^{{\text{OFF}}}{\text{ = }}1$ and $\gamma _k^{{\text{OFF}}} = \prod\limits_{i = 2}^k {\sqrt {{{{\lambda _{L + i, L + i - 1}}} \mathord{\left/ {\vphantom {{{\lambda _{L + i, L + i - 1}}} {{\lambda _{L + i - 1, L + i}}}}} \right. } {{\lambda _{L + i - 1, L + i}}}}} }$ ($2 \leqslant k \leqslant K$). In addition, we define ${{\bf{\tilde S}}^\vartheta }(t) = {{\bf{R}}^\vartheta }{{\bf{S}}^\vartheta }(t)$ and ${{\bf{\tilde T}}^\vartheta } = {{\bf{R}}^\vartheta }{{\bf{T}}^\vartheta }{\left({{{\bf{R}}^\vartheta }} \right)^{ - 1}}$, where $\vartheta = {\text{ON}}$ or ${\text{OFF}}$. Note that there is an orthogonal matrix ${{\bf{\tilde V}}^\vartheta }$ such that ${\left({{{{\bf{\tilde V}}}^\vartheta }} \right)^{\text{T}}}{{\bf{\tilde T}}^\vartheta }{{\bf{\tilde V}}^\vartheta } = {{\bf{\Lambda }}^\vartheta }$. Then, vector ${{\bf{S}}^\vartheta }(t)$ for the survival probabilities can be analytically expressed as

where ${{\bf{S}}^\vartheta }(0)$ is the vector of initial survival probabilities determined by the steady-state distribution of the promoter states. To that end, the dwell-time distribution at the active and inactive states can be calculated according to

where $\vartheta = {\text{ON}}$ or ${\text{OFF}}$.

After having obtained two dwell-time distributions at ON and OFF states, we can separately calculate the mean ON and OFF times, denoted by $\left\langle {{\tau ^{{\text{ON}}}}} \right\rangle$ and $\left\langle {{\tau ^{{\text{OFF}}}}} \right\rangle$. In addition, the ON dwell-time noise (denoted by $\eta _{{\text{ON}}}^2$) and OFF dwell-time noise (denoted by $\eta _{{\text{OFF}}}^2$) can also be calculated according to Eq (6). Furthermore, the mean burst frequency (BF) can be formally expressed as

2.2.3.   mRNA mean and noise

Note that the master equation for mRNA is different from that for burst size. Here we first establish the master equation for the former and then calculate the mean and noise of mRNA.

Let ${P_s}\left({m; t} \right)$ ($1 \leqslant s \leqslant N$) represent the probability that the number of mRNA molecules is $m$ at state ${A_1}, \cdot \cdot \cdot, {A_L}, {I_1}, \cdot \cdot \cdot, {I_K}$ at time $t$, respectively. Set $P\left({m; t} \right) = {\left[{{P_1}\left({m; t} \right), \cdot \cdot \cdot, {P_N}\left({m; t} \right)} \right]^{\text{T}}}$, and denote by $A = \left({{\lambda _{k, l}}} \right)$ the $N \times N$ transition matrix of a promoter. Note that ${\lambda _{k, l}} = 0$ means that there is no transition occurrence from state $k$ to state $l$. Let ${\bf{\Omega }} = {\text{diag}}\left({\mu, \mu, \cdot \cdot \cdot, \mu, 0, 0, \cdot \cdot \cdot, 0} \right)$ represent the transcription matrix. Then, the corresponding chemical master equation describing the mRNA dynamics takes the following form:

where ${\bf{E}}, {{\bf{E}}^{ - 1}}$ are shift operators and ${\bf{I}}$ is the identity operator.

Equation (8) can be solved using the binomial moment (BM) method that we previously developed ^[56]. The framework of this method is stated below. If the factorial BM for the probability ${P_s}\left({m; t} \right)$ is defined as $a_i^{\left(s \right)}\left(t \right) = \sum\limits_{m \geqslant i} {\left({\begin{array}{*{20}{c}} m \\ i \end{array}} \right){P_s}\left({m; t} \right)}$ ($1 \leqslant s \leqslant N$), then the total BM for the total probability defined as $P\left({m; t} \right) = \sum_{s = 1}^N {{P_s}\left({m; t} \right)\left(t \right)}$ is calculated by ${b_i}\left(t \right) = \sum_{s = 1}^N {a_{_i}^{(s)}\left(t \right)}$. Note that ${b_0}\left(t \right) = 1$ due to conservative condition: $\sum_{m = 0}^\infty {P\left({m; t} \right)} = 1$. Furthermore, if we introduce column vector ${a_i}\left(t \right) = {\left({a_i^{\left(1 \right)}\left(t \right), \cdot \cdot \cdot, a_i^{\left(N \right)}\left(t \right)} \right)^T}$, it follows from Eq (8) that

which is a linear ordinary differential equation group and is easily solved. Once all $a_i^{\left(s \right)}\left(t \right)$ and further the total BMs ${b_i}\left(t \right)$ are obtained, the mRNA probability distribution can be reconstructed according to

In addition, if the mRNA noise is quantified by the ratio of the variance over the squared mean (called noise intensity and denoted by $\eta _{_{{\text{mRNA}}}}^2$), then $\eta _{_{{\text{mRNA}}}}^2$ is calculated according to

3.   Analytical results

Here we present analytical results, obtained according to the above analysis framework, to help understand the mechanism of how the cyclic promoter structure affects transcriptional bursting kinetics.

3.1.   General results

First, we present analytical results for burst-size distribution. Let $P_s^{(0)}$ represent the steady-state probability that the promoter is in state $s$, where $s = 1, 2, \cdot \cdot \cdot, N,$ then these probabilities can be analytically derived (see Appendix B for details). Interestingly, we find that initial conditions ${S_{n, l}}\left(0 \right)$ are directly related to the steady-state distributions at the OFF states. Specifically, ${S_{0, 1}}\left(0 \right) = {{P_N^{(0)}{\lambda _{N, 1}}} \mathord{\left/ {\vphantom {{P_N^{(0)}{\lambda _{N, 1}}} {\left({P_N^{(0)}{\lambda _{N, 1}} + P_{L + 1}^{(0)}{\lambda _{L + 1, L}}} \right)}}} \right. } {\left({P_N^{(0)}{\lambda _{N, 1}} + P_{L + 1}^{(0)}{\lambda _{L + 1, L}}} \right)}},$ ${S_{0, L}}\left(0 \right) = {{P_{L + 1}^{(0)}{\lambda _{L + 1, L}}} \mathord{\left/ {\vphantom {{P_{L + 1}^{(0)}{\lambda _{L + 1, L}}} {\left({P_N^{(0)}{\lambda _{N, 1}} + P_{L + 1}^{(0)}{\lambda _{L + 1, L}}} \right)}}} \right. } {\left({P_N^{(0)}{\lambda _{N, 1}} + P_{L + 1}^{(0)}{\lambda _{L + 1, L}}} \right)}},$ and ${\kern 1pt} {\kern 1pt} {S_{n, l}}\left(0 \right) = 0$ for other subscripts. Owing to the special structure of matrix ${\bf{H}}(z)$, generating function ${G_l}\left({z, t} \right)$ can be analytically expressed as

where $- \lambda _l^h\left({1 \leqslant l \leqslant L} \right)$ are eigenvalues of matrix ${\bf{H}}\left(z \right)$ and thus are functions of $z$, and $\tilde V_{lj}^h\left({1 \leqslant l, j \leqslant L} \right)$ are the elements of orthogonal eigenvector ${\bf{\tilde V}}_{l}^h{\text{ = }}{\left({\tilde V_{1l}^h, \tilde V_{2l}^h, \cdots, \tilde V_{Ll}^h} \right)^T}$ associated with the eigenvalues $- \lambda _l^h$. Furthermore, generating function ${G_{{\text{burst}}}}\left(z \right)$ for burst-size distribution ${P_{{\text{burst}}}}\left(n \right)$ can be expressed as

Thus, burst-size distribution ${P_{{\text{burst}}}}\left(n \right)$ can be given via ${P_{{\text{burst}}}}\left(n \right) = \left({{1 \mathord{\left/ {\vphantom {1 {n!}}} \right. } {n!}}} \right){\left. {{{{d^n}{G_{{\text{burst}}}}\left(z \right)} \mathord{\left/ {\vphantom {{{d^n}{G_{{\text{burst}}}}\left(z \right)} {d{z^n}}}} \right. } {d{z^n}}}} \right|_{z = 0}}$.

Second, we present analytical results for dwell-time distributions. Note that $S_l^{{\text{ON}}}\left(t \right)$ and $S_k^{{\text{OFF}}}\left(t \right)$ for survival probabilities at the ${l^{th}}$ active state and at the ${k^{th}}$ inactive state can be respectively expressed as

where $\tilde V_{lj}^{{\text{ON}}}\left({1 \leqslant l, j \leqslant L} \right)$ and $\tilde V_{kj}^{{\text{OFF}}}\left({1 \leqslant k, j \leqslant K} \right)$ are the elements of matrices ${{\bf{\tilde V}}^{{\text{ON}}}} = \left({\tilde V_{lj}^{{\text{ON}}}} \right)$ and ${{\bf{\tilde V}}^{{\text{OFF}}}} = \left({\tilde V_{kj}^{{\text{OFF}}}} \right)$, respectively.

Using Eqs (14a) and (14b), the dwell-time distributions are thus expressed as

where if $\vartheta = {\text{ON}}$, then

and if $\vartheta = {\text{OFF}}$, then

Equation (15) indicates that the dwell-time distribution at ON or OFF states is in general the sum of exponential distributions.

Third, we present analytical results for BMs. We find that the steady-state total BMs take the following form

where ${{\bf{u}}_N} = \left({1, 1, \cdot \cdot \cdot, 1} \right)$ is an $N -$ dimensional row vector, ${\left({i\delta {\bf{I}} - {\bf{A}}} \right)^*}$ and $\det \left({i\delta {\bf{I}} - {\bf{A}}} \right)$ are respectively the adjacency matrix and the determinant of matrix $\left({i\delta {\bf{I}} - {\bf{A}}} \right),$ and the initial vector ${a_0}$ is the steady-state distribution of promoter states (see Appendix B).

In particular, the first-order binominal moment, i.e., the mRNA mean is given by (see Appendix C)

and the second-order binomial moment takes the form

where $W_{kl}^{(1)}$ represents a determinant of an $\left({N - 1} \right) \times \left({N - 1} \right)$ matrix that is the left part of matrix $\left({\delta {\bf{I}} - {\bf{A}}} \right)$ by crossing out its ${k^{th}}$ row and ${l^{th}}$ column elements, and $- \alpha _{_1}^{(1)},$ $- \alpha _{_2}^{(1)},$ $\cdot \cdot \cdot, - \alpha _{_{N - 1}}^{(1)}$ are the nonzero eigenvalues of matrix ${{\bf{A}} \mathord{\left/ {\vphantom {{\bf{A}} \delta }} \right. } \delta }$. Then, we can show that the Fano factor for mRNA, defined as the ratio of the variance over the mean, is in general not equal to one, implying that the mRNA distribution is in general not Poissonian but may be super-Poissonian or sub-Poissonian.

We point out that all the above analytical results are formal and tell us the finite information on the effect of the cyclic promoter structure on transcriptional bursting kinetics. In order to clearly see how the cyclic promoter structure affects transcriptional bursting kinetics, we will consider special cases.

3.2.   Specific results

Here we consider two representative cases: (1) multi-OFF mechanism ^[46], i.e., $L = 1$ and $K = N - 1$; (2) mixture mechanism where all the forward (i.e., clockwise) transition rates of active states are identical (denoted by ${\lambda _f}$) and all the backward (i.e., counterclockwise) transitions of active states are also identical (denoted by ${\lambda _b}$); and all the forward transition rates of inactive states are the same (denoted by ${\rho _f}{\lambda _f}$ with ${\rho _f}$ being a common forward ratio) whereas all the backward transition rates of inactive states are also the same (denoted by ${\rho _b}{\lambda _b}$ with ${\rho _b}$ being a common backward ratio).

For the multi-OFF mechanism, according to general formula Eq (13), it is not difficult to show that the generating function for burst size is given by

Thus, the burst-size distribution takes the form

which is a geometric distribution. These results are in agreement with previous assumptions ^[44].

For the mixture mechanism, we denote $a = \frac{{{\lambda _f}}}{{{\lambda _b}}}$ and $b = \frac{{{\rho _f}{\lambda _f}}}{{{\rho _b}{\lambda _b}}}$. Using general formula Eq (15) and noting the fact: $\sum\limits_{l = 1}^A {{{\sin }^2}\frac{{l\pi }}{{A + 1}}} = \frac{{A + 1}}{2}\, \left({A = L, K} \right),$ we can show that

where $j = 1, 2, \cdots, L$, and

where $j = 1, 2, \cdots, K$. Correspondingly, generating function ${G_{{\text{burst}}}}\left(z \right)$ for burst size is given by (see Appendix D for derivation)

where

Thus, the burst-size distribution takes the form

where $n = 0, 1, 2, \cdots$, and

${\theta _j} = \frac{{1 + a - 2\sqrt a \cos \left({{{j\pi } \mathord{\left/ {\vphantom {{j\pi } {\left({L + 1} \right)}}} \right. } {\left({L + 1} \right)}}} \right)}}{{{\mu \mathord{\left/ {\vphantom {\mu {{\lambda _b}}}} \right. } {{\lambda _b}}} + 1 + a - 2\sqrt a \cos \left({{{j\pi } \mathord{\left/ {\vphantom {{j\pi } {\left({L + 1} \right)}}} \right. } {\left({L + 1} \right)}}} \right)}}$, ${w_j}{\text{ = }}\frac{{{{{q_j}} \mathord{\left/ {\vphantom {{{q_j}} {{\lambda _b}}}} \right. } {{\lambda _b}}}}}{{1 + a - 2\sqrt a \cos \left({{{j\pi } \mathord{\left/ {\vphantom {{j\pi } {\left({L + 1} \right)}}} \right. } {\left({L + 1} \right)}}} \right)}}$ with $\sum\limits_{j = 1}^L {{w_j}} = 1$.

Equation (23) indicates that burst-size distribution is a mixture of geometric distributions. At the same time, the mean burst size is given by

In addition, we can also calculate the burst-size noise, which is given by

where ${r_j} = {\lambda _b}\left({1 + a - 2\sqrt a \cos \left({{{j\pi } \mathord{\left/ {\vphantom {{j\pi } {\left({L + 1} \right)}}} \right. } {\left({L + 1} \right)}}} \right)} \right)$. The corresponding noise will contribute to the mRNA noise.

4.   Analytical results

In contrast to previous studies that focused on the effect of multi-OFF mechanisms on gene expression, here we numerically demonstrate influences of multi-ON mechanisms on transcriptional bursting kinetics.

4.1.   Multi-ON mechanism can lead to bimodal burst-size distribution

Bursting kinetics is commonly characterized by burst size (BS) and burst frequency (BF). Here, we analyze how the multi-ON mechanism affects burst size and its distribution. Previous studies of the two-state gene model assumed that burst size follows a geometrical distribution ^[1,12], but whether or not this assumption is reasonable is unclear. In particular, it is not clear whether the burst size still obeys a geometric distribution in the case of multi-ON mechanisms. By performing numerical calculations, we find that multi-ON mechanism can induce bimodal distribution of burst sizes, referring to Figure 2(A)–(C). Figure 2(A), (B) show two representative unimodal and bimodal burst-size distributions, respectively, and Figure 2(C) shows how the number of the most probable mRNA molecules obtained by a statistical method depends on active index $L$. In Figure 2(A), (B), histograms show stochastic simulations obtained using the Gillespie stochastic simulation algorithm (SSA) ^[57], and solid lines are theoretical predictions. Figure 2(C) shows that the bimodal mRNA distribution exists for a large range of $L$. Specifically, the burst-size distribution is unimodal at the origin with $L = 1$ (i.e., multi-OFF model), a geometric distribution as predicted by Eq (20). Then, the burst-size distribution becomes a single peak away from the origin if we increase $L$, e.g., $L = 2, 3$. Finally, the burst-size distribution becomes bimodal with one peak at the origin and another peak away from the origin for larger $L$ with $L \geqslant 4$. These results imply that burst-size distribution in the case of the multi-ON mechanism is non-geometric but a mixture of geometric distributions as predicted by Eq (23).

Next, we analyze the influences of the multi-ON mechanism on mean burst size and burst-size noise, referring to Figure 2(D)–(F), where three representative promoter loops are shown. In all loops, we observe that mean burst size is always a monotonically increasing function of the active index, implying that more active states can lead to a larger burst size. In addition, we observe that burst-size noise is a monotonically decreasing function of the active index in the loops where the forward ratio is equal to the backward ratio, and the forward ratio is smaller than the backward ratio. However, there is a non-monotonic relationship between the active index and burst-size noise if the forward ratio is larger than the backward ratio. These results indicate that burst-size noise in the case of the multi-ON mechanism can exhibit different characteristics.

Previous studies showed that genes with a multistep mechanism could have peaked burst-size distribution and reduce noise in burst size ^[44]. Here, we point out that the multi-ON mechanism can induce more complex behaviors of burst size, such as bimodal burst-size distribution and non-monotonic burst-size noise.

4.2.   Multi-ON mechanism can lead to antagonistic timing of ON and OFF

Here we focus on analyzing how the multi-ON mechanism influences burst timing between transcription initial events by performing numerical calculations, referring to Figure 3. The first and second rows of Figure 3 show two representative OFF dwell-time distributions (Figure 3(A), (B)) and ON dwell-time distributions (Figure 3(D), (E)), as well as how the times for peak appearance in OFF dwell-time distribution (Figure 3(C)) and ON dwell-time distribution (Figure 3(F)) depend on active index $L$, respectively. First, we observe that OFF dwell-time distribution changes from unimodality to bimodality, whereas ON dwell-time distribution varies from bimodality to unimodality with the increase of active index $L$ (Comparing Figure 3(A)–(C) with Figure 3(D)–(F)). For example, the OFF dwell-time distribution is bimodal with one peak at the origin and another peak away from the origin for a broad range of $L$ with $1 \leqslant L \leqslant 8$, and then this distribution becomes unimodal with a single peak away from the origin if $L = 9$ and further becomes unimodal with a peak at the origin if $L = 10$ (see Figure 3(C)). However, for the ON dwell time, the distribution is unimodal with a peak at the origin if $L = 1$, and is an exponential distribution corresponding to a multi-OFF model. The distribution changes from unimodality with one single peak away from the origin if $L = 2$ to bimodality with one peak at the origin and another peak away from the origin for a broader range of $L$ : $3 \leqslant L \leqslant 10$ (see Figure 3(F)).

The bottom row of Figure 3 shows the mean and noise of OFF dwell time (Figure 3(G)) and ON dwell-time (Figure 3(H)), respectively, as well as mean burst frequency for three representative loops (Figure 3(I)). We observe that with the rise of the active index, the mean OFF time decreases whereas mean ON time increases, but the OFF time noise increases and ON time noise decreases. In addition, we find that with increasing the active index, the mean burst frequency can exhibit different characteristics, such as no change with ${\rho _f} = {\rho _b} = 1,$ decrease with ${\rho _f} = 0.5, \; {\rho _b} = 2,$ and increase with ${\rho _f} = 2, \; {\rho _b} = 0.5$.

The above analysis implies that the multi-ON mechanism can give rise to antagonistic timing of ON and OFF, and induce different characteristics of burst frequency, hinting that the multi-ON mechanism has a considerable influence on transcriptional bursting timing.

4.3.   Propagation of bursting kinetics to mRNA variability

Single-molecule experiments have verified that mRNAs are produced often in a bursting manner ^[1,2]. It is of interest to understand how bursting kinetics affect the variations in mRNA levels. Previous studies focused on the influence of burst parameters such as mean burst size and burst frequency on gene expression levels and noise by analyzing queuing models ^{[32,33,44,45]}. However, how the burst size and timing affect gene expression dynamics in a multistate promoter remains elusive (see Figure 1(C)). Here, we explore the influences of burst-size distribution and dwell-time distribution on mRNA distribution, and the numerical results are shown in Figure 4(A)–(H).

From Figure 4(A), (B), we first observe that the burst-size distribution, OFF dwell-time distribution, and ON dwell-time distribution are bimodal, and the mRNA distribution is also bimodal. From Figure 4(C), (D), we find that the burst-size distribution, OFF dwell-time distribution, and ON dwell-time distribution are unimodal, but the mRNA distribution is bimodal. These results indicate that both unimodal and bimodal distributions of burst size and dwell time can lead to bimodal mRNA distribution. On the other hand, we find that the bimodal or unimodal distributions of burst size, OFF dwell time, and ON dwell time can lead to unimodal mRNA distribution, referring to Figure 4(E)–(H). A previous study showed that unimodal or bimodal gene expression has important biological implications, i.e., a cause of phenotypic diversity in a population of genetically identical cells ^[58]. These results indicate that the variation in the distributions of burst size and dwell-time can give rise to variability in the distribution of gene expression production and further phenotypic variation.

Next, we analyze the influence of burst-size noise, OFF dwell-time noise, and ON dwell-time noise on mRNA noise, with numerical results shown in Figure 5. We observe that the changes in burst-size noise, OFF dwell-time noise, and ON dwell-time noise have different characteristics such as increasing, decreasing, and non-monotonic with the increase of active index $L$, referring to Figure 5(A), (C), and (E). However, the mRNA noise is always a monotonically decreasing function of active index $L$ (see the last column of Figure 5). These results indicate that the integration of the burst-size noise and dwell-time noise may reduce gene expression noise. Our previous study also showed this by analyzing multi-OFF models ^[46]. Together, the results imply that a multistep process is a universal mechanism of reducing gene expression noise.

The above analysis shows that the multi-ON mechanism can induce complex transcription kinetics, e.g., unimodal or bimodal distribution of mRNA counts, different burst sizes and timing, and different characteristics of noise. The obtained results imply that the propagation of bursting kinetics to gene expression is complex, and would be important for understanding cellular processes and variation in cell populations.

5.   Discussion and conclusions

Transcription is an inherently stochastic process and often occurs in a bursting fashion. Because of involving recruitment of TFs and DNA polymerases, chromatin remodeling, and sequence of transitions between activity states of the promoter, the kinetics of transcriptional bursting is poorly understood. Previous studies showed that mRNA burst-size distribution is geometric with exponential waiting time for two-state promoter models ^[38,40], but how transcription burst size and timing in multistate promoter models remains unclear.

In the present work, we have analyzed a gene model that corporate the complexity of promoter structure, focusing on the effects of multi-ON mechanisms on transcriptional bursting kinetics. By establishing the survival probability master equations for mRNA bursting production and promoter dwelling at ON/OFF states, we have successfully derived analytical expressions for burst-size distribution and ON/OFF waiting time distributions, and have found that the burst size is a mixed geometric distribution and the ON/OFF dwell time is a mixed exponential distribution in the case of the mixture mechanism. These obtained distributions can be used to analyze transcriptional burst size and timing for arbitrarily complex promoter structure shown in Figure 1.

While previous studies on multi-OFF models ^[47] and multi-ON models ^[48] focused on analysis of mRNA distribution and transcriptional noise, in the present work we have explored the effects of the multi-ON mechanism on transcriptional bursting kinetics. We have shown that this mechanism can lead to bimodal burst-size distribution and different patterns of burst-size noise, which is different from previous results that the multistep mechanism can give rise to a peaked burst-size distribution and a reduction in the noise of burst size ^[44]. In addition, we found that the multi-ON mechanism can lead to antagonistic timing of ON and OFF. For example, the OFF waiting time distribution changes from bimodality to unimodality, whereas the ON waiting time distribution changes from unimodality to bimodality, as the active index increases. Finally, by exploring the distributions and noise for burst size, OFF/ON dwell time, and mRNA copy-number, we found that the propagation of bursting kinetics to mRNA variability is complex. In contrast to previous studies on two-state gene models, which showed that the lower expression noise could result from higher burst frequency, or smaller burst size ^[35,40,59], here we showed that the qualitative conclusions no longer hold in multistate promoter models. More specifically, we find that mRNA noise always reduces, but the burst-size noise, OFF/ON dwell-time noise, and burst frequency can exhibit complex behaviors in the case of the multi-ON mechanism.

Recently, many studies aimed at the inference of kinetic parameters from experimental data on single cells ^[50,52]. In a previous study, we showed that different ON/OFF waiting time distributions can exhibit identical steady-state mRNA distributions ^[48], implying that when one infers a gene system based on transcription dynamics obtained by measurements of fluctuations of mRNA levels, caution must be kept. However, our results obtained here would be used to infer the promoter structure and investigate the variability of the transcription output based on the experimental data. It should be pointed out that the intermediate processes of gene expression such as the partitioning of mRNA at cell division ^[60], the alternative splicing ^[61], and the feedback loop ^[62], would additionally influence the obtained results, but they are neglected in our study. How details of these factors affect transcriptional bursting kinetics is worth further investigation.

Acknowledgments

This work was supported by the National Nature Science Foundation of P. R. China (Nos. 11601094, 12171494, 11631005, 11775314, 11931019), Guangdong Key Research and Development Project (2019B0233002), and Guangdong Province Key Laboratory of Computational Science at the Sun Yat-sen University (2020B1212060032).

Conflict of interest

The authors have declared that no competing interest exists.

Appendix

A.   Expression for coefficient matrix ${\bf{H}}(z)$

According to the definition of generating function and Eq (3) in the main text, we can show the expression for the coefficient matrix

where ${{h}_{l}}\left(z \right) = \mu \left(z-1 \right)-{{\delta }_{l1}}{{\lambda }_{l, N}}-\left(1-{{\delta }_{l1}} \right){{\lambda }_{l, l-1}}-{{\lambda }_{l, l+1}}, \ l = 1, 2, \cdot \cdot \cdot, L.$

B.   Derivation of steady-state probability for promoter state

According to Table 1 in the main text, we can establish the following master equation for promoter kinetics

where $\bf{P}(t) = {{\left({{P}_{1}}\left(t \right), {{P}_{2}}\left(t \right), \cdots, {{P}_{N}}\left(t \right) \right)}^{\text{T}}}$ with ${{P}_{k}}\left(t \right)$ being the time-dependent probability of promoter at state k, and A is the transcription matrix of promoter state. Our interest is in deriving the steadystate probability ${{\bf{P}}^{\left(0 \right)}} = \left(P_{1}^{\left(0 \right)}, P_{2}^{\left(0 \right)}, \cdots, P_{N}^{\left(0 \right)} \right)$ satisfying the condition: $\bf{A}{{\bf{P}}^{\left(0 \right)}} = 0$.

First, note that the transition matrix ${\bf{A}}$ is an M-matrix (i.e., the sum of every column elements is equal to zero). So we can assume the eigenvalues of the matrix ${\bf{A}}$ as $0, - {\alpha _1}, - {\alpha _2} \cdot \cdot \cdot, - {\alpha _{N - 1}}$ with ${\alpha _s} \ne 0$, and the eigenvalues of the matrix ${A_k}$ (i.e., crossing out the ${k^{th}}$ row and ${k^{th}}$ column entry of matrix ${\bf{A}}$) as $- \beta _1^{\left(k \right)}, - \beta _2^{\left(k \right)}, \cdots, - \beta _{N - 1}^{\left(k \right)}$ with $\beta _s^{\left(k \right)} \ne 0$. This assumption implies that the characteristic polynomial of ${\bf{A}}$ and ${{\bf{A}}_k}$ takes the following form, respectively

Note that Laplace's formula for the M-matrix A gives

Since the null space of A is one-dimensional, we have ${{\bf{P}}^{\left(0 \right)}} = \rho {{\left[\det \left({{\bf{A}}_{1}} \right), \det \left({{\bf{A}}_{2}} \right), \cdot \cdot \cdot, \det \left({{\bf{A}}_{N}} \right) \right]}^{\text{T}}}$. From the conservative condition ${{\bf{u}}_{N}}{{\bf{P}}^{\left(0 \right)}} = 1$ with ${{\bf{u}}_{N}} = \left(1, 1, \cdot \cdot \cdot, 1 \right)$ being a row vector, we can obtain

From Eq (B2), we have ${{\left. \frac{d{{g}_{\bf{A}}}\left(x \right)}{dx} \right|}_{x = 0}} = \prod\limits_{s = 1}^{N-1}{{{\alpha }_{s}}}$. On the other hand, Jacobi's formula gives

where "tr" represents the trace of a matrix. The combination of both yields the equality $\sum\limits_{k = 1}^{N}{\prod\limits_{s = 1}^{N-1}{\beta _{s}^{\left(k \right)}}} = \prod\limits_{s = 1}^{N-1}{{{\alpha }_{s}}}$. Finally, we obtain the explicit expression of steady-state probability for promoter at every state

C.   Derivation of mRNA mean and noise

First, we compute the binomial moments (BMs). Let $W_{kl}^{(i)}$ represent a determinant of a $\left(N-1 \right)\times \left(N-1 \right)$ matrix, which is from matrix $\left(i\delta {\bf{I}}-{\bf{A}} \right)$ by crossing out its k^th row and l^th column elements. Denoting by $\bf{W}_{k}^{(i)} = {{\left(W_{k1}^{(i)}, \cdots, W_{kN}^{(i)} \right)}^{T}}$ a vector, then the adjacency matrix of $\left(i\delta {\bf{I}}-{\bf{A}} \right)$ can be written as

Therefore, it is easy to get the following expression

In addition, the determinant is

where $-\alpha _{_{1}}^{(i)}, \cdot \cdot \cdot, -\alpha _{_{N-1}}^{(i)}$ are the nonzero eigenvalues of M-matrix ${{\bf{\tilde{A}}}^{(i)}} = {\bf{A}}/{i\delta }\; .$ Substituting Eqs (C2) and (C3) into expression ${{a}_{i}} = \frac{{{\left(i\delta {\bf{I}}-{\bf{A}} \right)}^{*}}{\bf{\Omega }}}{\det \left(i\delta {\bf{I}}-{\bf{A}} \right)}{{{\bf{a}}}_{i-1}}$, which can be derived from Eq (9) in the main text, we can obtain the expressions of the factorial BMs

Finally, we can get the expression of the total BMs

Equations (C4) and (C5) imply that the vectors of BMs a_i and total BMs b_i can be analytically given based on a₀, i.e., on the steady-state probabilities of promoter states.

D.   Derivation the explicit expression of generating function for the mixture loop model

For the mixture loop model, matrix ${\bf{H}}\left(z \right)$ becomes

where $h\left(z \right) = \mu \left(z-1 \right)-{{\lambda }_{f}}-{{\lambda }_{b}}$. Note that matrix ${\bf{H}}\left(z \right)$ has eigenvalues $-\lambda _{1}^{h}, -\lambda _{2}^{h}, \cdots, -\lambda _{L}^{h}$ with

In addition, the eigenvectors associated with the eigenvalues $- \lambda _j^h$ are given by ${\bf{V}}_j^h = {\left({V_{1j}^h, \cdots, V_{Lj}^h} \right)^{\text{T}}}$ with $V_{lj}^h = \sqrt {\frac{2}{{L + 1}}} {\left({\frac{{{\lambda _f}}}{{{\lambda _b}}}} \right)^{{{\left({l - 1} \right)} \mathord{\left/ {\vphantom {{\left({l - 1} \right)} 2}} \right. } 2}}}\sin \frac{{lj\pi }}{{L + 1}}$. If denoting by ${{\bf{V}}^h} = \left({{\bf{V}}_1^h, \cdots, {\bf{V}}_L^h} \right),$ then we have ${\left({{{\bf{V}}^h}} \right)^{ - 1}}{\bf{H}}\left(z \right){{\bf{V}}^h} = {{\bf{\Lambda }}^h} = {\text{diag}}{\left({ - \lambda _1^h, \cdots, - \lambda _L^h} \right)^T}.$ Making the transform ${{\bf{\tilde V}}^h}{\text{ = }}{{\bf{R}}^{{\text{ON}}}}{{\bf{V}}^h},$ and we obtain the orthogonal matrix ${{\bf{\tilde V}}^h} = \left({\tilde V_{lj}^h} \right)$ with $\tilde V_{lj}^h = \sqrt {\frac{2}{{L + 1}}} \sin \frac{{lj\pi }}{{L + 1}}.$ Substituting this expression and Eq (D2) into Eq (13) in the main text, we obtain the explicit expression of generating function in the case of the mixture mechanism.