Model formulationFor the two-stage gene expression model (1), the probability \(p_{m,n}(t)\) of observing m mRNA and n protein molecules at time t satisfies the CME$$\begin{aligned} \begin{aligned} \frac{\textrm{d}}{\textrm{d} t} p_{m,n} =&\lambda ^m (p_{m-1,n} – p_{m,n}) + \gamma ^m ((m+1) p_{m+1,n} – m p_{m,n} ) \\&+ \lambda ^p m (p_{m,n-1} – p_{m,n}) + \gamma ^p ((n+1) p_{m,n+1} – n p_{m,n}), \end{aligned} \end{aligned}$$
(11)
subject to initial condition$$\begin{aligned} p_{m,n}(0) = \delta _{m,m_0} \delta _{n,n_0}, \end{aligned}$$where \(\delta _{i,j}\) represents the Kronecker delta symbol, which is one if \(i=j\) and zero otherwise; \(m_0\) and \(n_0\) are the initial mRNA and protein amounts, respectively.Our aim is to obtain a PDE rather than working with the CME (11). To this end, we introduce the probability generating function defined by$$\begin{aligned} G(x,y,t) = \sum _{m} \sum _{n} x^m y^n p_{m,n}(t). \end{aligned}$$
(12)
Multiplying the CME (11) by the factor \(x^m y^n\) and summing over all m and n, and using (12), we arrive at the generating function which satisfies the linear first-order PDE$$\begin{aligned} \frac{\partial G}{\partial t} = (\gamma ^m (1-x) + \lambda ^p x (y-1) ) \frac{\partial G}{\partial x} + \gamma ^p (1-y) \frac{\partial G}{\partial y} + \lambda ^m (x-1) G. \end{aligned}$$
(13)
Equation (13) has been used in [24] to derive mRNA and protein moments; it has been solved at steady state in [18]. Here we shall derive and study a generalisation of (13).Without loss of generality, for the generalised model (2), the probability \(P(\textbf{m}, n, t)\) of observing \(m_1\) mRNA copies in state 1, \(m_2\) mRNA copies in state 2, and so on, at given time t satisfies the following CME,$$\begin{aligned} \begin{aligned} \frac{\text {d} P(\textbf{m},n, t)}{\textrm{dt}} =&\sum _{i=1}^{K} \left( \lambda _i^m (\mathbb {E}^{-1}_i-1)P + \gamma _i^m (\mathbb {E}_i-1) m_i P + \sum _{j=1}^{K} q_{ij} (\mathbb {E}_i\mathbb {E}^{-1}_j-1)\right. \\&\times m_i P + \left. \lambda _i^p (\mathbb {E}^{-1}_{K+1}-1) m_i P \right) + \gamma ^p (\mathbb {E}_{K+1}-1) n P, \end{aligned} \end{aligned}$$
(14)
where \(\textbf{m} = \begin{bmatrix} m_1&m_2&m_3&\ldots&m_K \end{bmatrix}\) is a vector of species copy numbers. Note that the step operator [36] \(\mathbb {E}_i\) in (14) is in the variable \(m_i\), whereas \(\mathbb {E}_{K+1}\) is in the variable n; \(\mathbb {E}_i \mathbb {E}^{-1}_j-1 = 0\) for \(i=j\).The multivariate probability generating function is given by$$\begin{aligned} G(\textbf{x}, y,t) = \sum _{m_1}\cdots \sum _{m_K}\sum _{n}P(\textbf{m},n,t) x_1^{m_1} x_2^{m_2}\cdots x_K^{m_K} y^n, \end{aligned}$$
(15)
where \(\textbf{x} = \begin{bmatrix} x_1&x_2&x_3&\ldots&x_K \end{bmatrix}\). Multiplying (14) by \(x_1^{m_1} x_2^{m_2}\ldots x_K^{m_K} y^n\) and summing over all \(m_1, m_2, \ldots , m_K, n\), and employing (15), we arrive at the PDE$$\begin{aligned} \begin{aligned} \frac{\partial G(\textbf{x},y,t)}{\partial t} =&\sum _{i=1}^{K} \left( \lambda _i^m (x_i-1) G + \gamma _i^m (1 – x_i) \frac{\partial G}{\partial {x_i}} + \sum _{j=1}^{K} q_{ij} (x_j – x_i) \frac{\partial G}{\partial {x_i}} \right. \\&+ \left. \lambda _i^p (y-1) x_i \frac{\partial G}{\partial {x_i}} \phantom {\sum _{i=1}^{K}} \right) + \gamma ^p (1-y) \frac{\partial G}{\partial {y}}. \end{aligned} \end{aligned}$$
(16)
Note that the step operators \(\mathbb {E}_i^{\pm 1}\) in (14) coincide with the variables \(x_i^{\mp 1}\) while the copy number of species \(m_i\) correspond to the terms \(x_i \partial _{x_i}\) in (16) for the generating function. In the next section, we will seek a solution to the PDE (16).SolutionIn this section, we shall provide a step-by-step breakdown of our solution method for solving the PDE (16). We are interested in the steady state; therefore, we set the time derivative in (16) to zero and rearrange the resulting equation to obtain$$\sum\limits_{i = 1}^K {\left( {\lambda _i^m({x_i} – 1)G + \left( {\gamma _i^m(1 – {x_i}) + \sum\limits_{j = 1}^K {{q_{ij}}} ({x_j} – {x_i}) + \lambda _i^p(y – 1){x_i}} \right)\frac{{\partial G}}{{\partial {x_i}}}} \right)} + {\gamma ^p}(1 – y)\frac{{\partial G}}{{\partial y}} = 0$$
(17)
for the time-independent generating function \(G(\textbf{x},y)\) of the stationary distribution. The probability normalisation condition translates to \(G(1,\ldots ,1) = 1\). Changing the variables according to$$\begin{aligned} x_i = 1 + u_i, \quad y = 1 + v, \quad G = \exp (\varphi ) \end{aligned}$$
(18)
allows us to transform (17) into$$\begin{aligned} \sum _{i=1}^{K} \left( \lambda _i^m u_i + \left( \lambda _i^p v(1+u_i) – \gamma _i^m u_i + \sum _{j=1}^{K} q_{ij} (u_j – u_i) \right) \frac{\partial \varphi }{\partial u_i} \right) = \gamma ^p v \frac{\partial \varphi }{\partial v}, \end{aligned}$$
(19)
which is subject to the normalisation condition$$\begin{aligned} \varphi (\textbf{0}) = 0. \end{aligned}$$
(20)
Below, we focus on seeking a solution to (19)–(20) using a suitable ansatz.Let us first consider that the solution is of the form$$\begin{aligned} \varphi (u_1, u_2, u_3, \ldots ,u_K, v) = \varphi _0(v) + u_1 \varphi _1(v) + \ldots + u_K \varphi _K(v) . \end{aligned}$$
(21)
With this in mind, we obtain from (21) that$$\begin{aligned} \frac{\partial \varphi }{\partial u_i} = \varphi _i(v), \quad \frac{\partial \varphi }{\partial v} = \varphi _0^{\prime }(v) + u_1 \varphi _1^{\prime }(v) + \ldots + u_K \varphi _K^{\prime }(v). \end{aligned}$$
(22)
Inserting the partial derivatives (22) into (19), we get$$\begin{aligned} \begin{aligned} \sum _{i=1}^{K} \left( \lambda _i^m u_i + \left( \lambda _i^p v (1+u_i)- \gamma _i^m u_i + \sum _{j=1}^{K} q_{ij} (u_j – u_i) \right) \varphi _i – \gamma ^p v u_i \varphi _i^\prime \right) \\ = \gamma ^p v \varphi _0^\prime . \end{aligned} \end{aligned}$$
(23)
Equation (23) can be rewritten as$$\left( {{\gamma ^p}\varphi _0^\prime – \sum\limits_{i = 1}^K {\lambda _i^p} {\varphi _i}} \right)\nu + \sum\limits_{i = 1}^K {\left( {{\gamma ^p}v\varphi _i^\prime + \left( {\gamma _i^m – \lambda _i^p\nu + \sum\limits_{j = 1}^K {{q_{ij}}} } \right){\varphi _i} – \sum\limits_{j = 1}^K {{q_{ji}}} {\varphi _j} – \lambda _i^m} \right)} {u_i} = 0.$$
(24)
In order that (24) hold, we must necessarily have$$\begin{aligned}&\sum _{i=1}^{K} \lambda _i^p \varphi _i – \gamma ^p \varphi _0^\prime = 0, \end{aligned}$$
(25)
$$\begin{aligned}&\gamma ^p v \varphi _i^\prime + \left( \gamma _i^m – \lambda _i^p v + \sum _{j=1}^{K} q_{ij} \right) \varphi _i – \sum _{j=1}^{K} q_{ji} \varphi _j = \lambda _i^m. \end{aligned}$$
(26)
Thus far, we have converted the system of PDEs (17) into the system of ODEs (25)–(26). Next, we provide a detailed explanation of solving this system using the power series method.Let us assume that the functions \(\varphi _0\) and \(\varphi _i\) are of the power series form, i.e.,$$\begin{aligned} \varphi _0(v) = \sum _{n=0}^{\infty } a_n v^n, \quad \varphi _i(v) = \sum _{n=0}^{\infty } b_n^{(i)} v^n \end{aligned}$$
(27)
for \(i \in \{1,\ldots ,K \}\). Differentiating (27) term by term we get$$\begin{aligned} \varphi ^{\prime }_0(v) = \sum _{n=1}^{\infty } n a_n v^{n-1}, \quad \varphi ^{\prime }_i(v) = \sum _{n=1}^{\infty } n b_n^{(i)} v^{n-1}. \end{aligned}$$
(28)
Inserting (27) and (28) into (26), and collecting same powers of v, we obtain the following system of recurrence relations$$\begin{aligned} \left( \gamma _i^m + \sum _{j=1}^{K}q_{ij} + n \gamma ^p \right) b_n^{(i)} – \sum _{j=1}^{K} q_{ji} b_n^{(j)} = \lambda _i^p b_{n-1}^{(i)} \end{aligned}$$
(29)
for the coefficients \(b_n^{(i)}\), where \(i=1,\ldots ,K\). For the sake of simplicity, equations (29) can be rewritten in matrix form as$$\begin{aligned} (\mathbf {A-Q^\top } + n \gamma ^p \textbf{I}) X_n = \textbf{B} X_{n-1}, \quad n \ge 1, \end{aligned}$$
(30)
where I is the identity matrix and the vector \(X_n\) is defined as$$\begin{aligned} X_n = \begin{bmatrix} b_n^{(1)},&b_n^{(2)},&b_n^{(3)},&\ldots ,&b_n^{(K)} \end{bmatrix}^\top . \end{aligned}$$In (30), A is a \(K \times K\) matrix defined by$$\begin{aligned} \textbf{A}_{ij} := {\left\{ \begin{array}{ll} \gamma _i^m & \text { for } i=j,\\ 0 & \text { for } i \ne j, \end{array}\right. } \end{aligned}$$
(31)
\(\textbf{Q}\) is a \(K \times K\) matrix defined by$$\begin{aligned} \textbf{Q}_{ij} := {\left\{ \begin{array}{ll} – \sum \limits _{k \ne i} q_{ik} & \text { for } i=j,\\ q_{ij} & \text { for } i \ne j, \end{array}\right. } \end{aligned}$$
(32)
and \(\textbf{B}\) is a \(K \times K\) matrix defined by$$\begin{aligned} \textbf{B}_{ij} := {\left\{ \begin{array}{ll} \lambda _i^p & \text { for } i=j,\\ 0 & \text { for } i \ne j. \end{array}\right. } \end{aligned}$$
(33)
In order to solve the recurrence relations (30) initial conditions are needed. These can be obtained from (26) by setting \(v=0\) for each \(i \in \{1,2,\ldots ,K \}\). The resulting system of linear equations is given in matrix form as$$\begin{aligned} \mathbf {(A-Q^\top )} X_0 = C, \end{aligned}$$
(34)
where C is a column vector defined as \(C = \begin{bmatrix} \lambda _1^m&\lambda _2^m&\ldots&\lambda _K^m \end{bmatrix}^\top\).Solving the system of algebraic equations (30) under the initial conditions (34) yields the terms of \(b_{n}^{(i)}\); the sequence \(a_n\) can be obtained by substituting (27) and (28) into (25) and collecting same powers of v. By doing so, we get$$\begin{aligned} a_{n} = \frac{1}{n \gamma ^p} \sum _{i=1}^{K} \lambda _i^p b_{n-1}^{(i)}, \quad n \ge 1. \end{aligned}$$
(35)
Note that the normalisation condition (20) implies that \(a_0 = \varphi _0(0) = \varphi (\textbf{0}) = 0\). Having found the sequences \(a_n\) and \(b_n^{(i)}\), we combine (21) and (27) to obtain$$\begin{aligned} \varphi (u,v) = \sum _{n=1}^{\infty } a_n u^n + \sum _{i=1}^{K} v_i \sum _{n=0}^{\infty } b_n^{(i)} u^n. \end{aligned}$$
(36)
We return to the original variables in (36) via (18) to obtain the generating function of the stationary distribution of mRNA and protein amounts, which is given by$$\begin{aligned} G(\textbf{x},y) = \exp \left( \sum _{n=1}^{\infty } a_n (y-1)^n + \sum _{i=1}^{K} (x_i – 1) \sum _{n=0}^{\infty } b_n^{(i)} (y-1)^n \right) . \end{aligned}$$
(37)
Equation (37) provides the sought-after steady-state solution to the PDE (16) and will be used in the following section.Marginal distributions and momentsIn this section, we use the analytical formula for the generating function (37) to obtain marginal mRNA distributions. We determine the moments of the protein distribution by way of the factorial cumulants, which allow us to recover the protein distribution. Additionally, we derive the protein Fano factor (variance-to-mean ratio) and express it in terms of the first two factorial moments.Marginal mRNA distributions In the generating function (37), if we take \(y=1\), then we obtain the marginal mRNA distributions as$$\begin{aligned} G^m(\textbf{x}) = G(\textbf{x}, 1) = \exp \left( \sum _{i=1}^{K} b_0^{(i)} (x_i – 1) \right) = \prod _{i=1}^{K} \exp \left( b_0^{(i)} (x_i – 1) \right) , \end{aligned}$$
(38)
from which we conclude that the steady state mRNA distributions are independent Poissons with means$$\begin{aligned} \langle m_i \rangle = b_0^{(i)}. \end{aligned}$$Marginal protein distribution Likewise, by inserting \(x_i = 1\) (\(i=1,\ldots ,K\)) into (37), we can recover the generating function of the marginal protein distribution$$\begin{aligned} G(y) = G(\textbf{1}, y)= \exp \left( \sum _{n=1}^{\infty } a_n (y-1)^n \right) , \end{aligned}$$
(39)
where \(\textbf{1}\) is a K-dimensional row vector of ones.Next, we determine the moments of the protein distributions. The factorial (combinatorial) moments \(h_n\) are obtained by expanding the generating function into a power series around \(y=1\):$$\begin{aligned} G(y) = \sum _{n=0}^{\infty } h_n (y-1)^n. \end{aligned}$$We aim to calculate the factorial moments \(h_n\) by way of the factorial cumulants \(a_n\). To that end, we first differentiate (39) to obtain$$\begin{aligned} \textrm{D}G(y) = G(y) \textrm{D} \ln G(y), \end{aligned}$$
(40)
where \(\textrm{D}\) denotes the differential operator \(\textrm{d}/\textrm{d}y\). Then, taking the \((n-1)\)th derivative of (40), we get$$\begin{aligned} \textrm{D}^n G(y) = \sum _{i=0}^{n-1} \left( {\begin{array}{c}n-1\\ i\end{array}}\right) \textrm{D}^i G(y) \textrm{D}^{n-i} \ln G(y), \end{aligned}$$which can be recast as$$\begin{aligned} \frac{\textrm{D}^n G(y)}{n!} = \sum _{i=0}^{n-1} \left( 1 – \frac{i}{n} \right) \frac{\textrm{D}^i G(y)}{i!} \frac{\textrm{D}^{n-i} (\ln G(y))}{(n-i)!}. \end{aligned}$$
(41)
Evaluating (41) at \(y=1\) gives the factorial moments of the protein distribution$$\begin{aligned} h_n = \sum _{i=0}^{n-1} \left( 1 – \frac{i}{n} \right) a_{n-i} h_i, \quad \text { for } n \ge 1, \end{aligned}$$
(42)
where \(h_0 = 1\). The terms of \(h_n\) can be recursively obtained by inserting (35) into (42). Subsequently, by employing the recurrence method proposed in [37], we recover the protein distribution$$\begin{aligned} p(n) = \sum _{j=1}^{\infty } \frac{(j+1)_n}{n!} h_{n+j}(-1)^j, \end{aligned}$$where \((x)_n\), n being a nonnegative integer, denotes the rising factorial or namely Pochhammer symbol.Moments Clearly, the mRNA distributions in (38) are Poissonian. Therefore, mRNA Fano factor is equal to 1. The protein mean and Fano factor can be derived from the factorial moments (42). The first two factorial moments are given by$$\begin{aligned} \langle n \rangle = h_1 = a_1 \quad \text {and} \quad \langle n(n-1) \rangle = 2h_2 = 2 a_2 + a_1^2, \end{aligned}$$
(43)
respectively. The Fano factor,$$\begin{aligned} \textrm{F} = \frac{\langle n^2 \rangle }{\langle n \rangle } – \langle n \rangle = \frac{\langle n(n-1) \rangle }{\langle n \rangle } + 1 – \langle n \rangle = \frac{2a_2}{a_1} + 1, \end{aligned}$$
(44)
is thus expressed in terms of the first two factorial cumulants \(a_1\) and \(a_2\).The mRNA inactivation loop modelIn this section, we present a particular example of the generalised model (2), which we refer to as the inactivation loop model, whose reaction scheme is given by (3). Specifically, we provide an explicit representation of the stationary solution using the cumulants. Furthermore, we calculate the steady-state protein Fano factor and express it as a function of the model parameters. Let us note that a possible biological scenario that can implement this model is by a regulatory RNA that temporarily blocks mRNA function [38].The inactivation loop model (3) can be readily obtained from the generalised model (2) by taking \(K=2\), which accounts for only two mRNA states denoting the active mRNA state \(m_1\) and the inactive mRNA state \(m_2\). In what follows, we assume that a newly produced mRNA is active, i.e. that the transcription rate satisfies$$\begin{aligned} \lambda _i^m = \lambda ^m \delta _{i,1}, \quad \text {for } i=1,2. \end{aligned}$$
(45)
Additionally, we assume that proteins are translated only from an active mRNA, so that we have$$\begin{aligned} \lambda _i^p = \lambda ^p \delta _{i,1}, \quad \text {for } i=1,2, \end{aligned}$$
(46)
for the translation rate. Here, \(\delta _{i,j}\) denotes the Kronecker delta symbol. Cumulants We aim to recover expressions for the inactivation loop model from the generalised model. The system of algebraic equations for this model follows from (34), taking the form of$$\begin{aligned} (\gamma _1^m + q_{12}) b_0^{(1)} – q_{21} b_0^{(2)}&= \lambda ^m, \\ (\gamma _2^m + q_{21}) b_0^{(2)} – q_{12} b_0^{(1)}&= 0, \end{aligned}$$from which we recover$$\begin{aligned} b_0^{(1)} = \frac{\lambda ^m (\gamma _2^m + q_{21})}{(\gamma _1^m + q_{12})(\gamma _2^m + q_{21}) – q_{12}q_{21}}. \end{aligned}$$
(47)
Combining (47) with (39) we find$$\begin{aligned} \langle m_1 \rangle = \frac{\lambda ^m}{\gamma _{\textrm{eff}}^m}, \end{aligned}$$where$$\begin{aligned} \gamma _{\textrm{eff}}^m = \gamma _1^m + \frac{q_{12} \gamma _2^m}{\gamma _2^m+q_{21}} \end{aligned}$$
(48)
is the effective rate of mRNA decay. The recurrence relations (30) read$$\begin{aligned}&(\gamma _1^m + q_{12} + n \gamma ^p) b_n^{(1)} – \lambda ^p b_{n-1}^{(1)} – q_{21} b_n^{(2)} = 0, \end{aligned}$$
(49)
$$\begin{aligned}&(\gamma _2^m + q_{21} + n \gamma ^p) b_n^{(2)} – q_{12} b_n^{(1)} = 0, \end{aligned}$$
(50)
for \(n \ge 1\). Solving the algebraic system (49)–(50) in \(b_n^{(1)}\) yields$$\begin{aligned} b_n^{(1)} = \frac{\lambda ^p (\gamma _2^m + q_{21} + n \gamma ^p)}{{\gamma ^p}^2 n^2 + \gamma ^p(\gamma _2^m + \gamma _1^m +q_{21}+q_{12}) n + \gamma _2^m \gamma _1^m +\gamma _1^m q_{21} + \gamma _2^m q_{12}} b_{n-1}^{(1)}, \end{aligned}$$
(51)
which is a recursive expression whose first term (i.e. zeroth) is given by (47).Explicit representation The recursive formula (51) can further be simplified by factorising its denominator as$$\begin{aligned} b_n^{(1)} = \lambda ^p \frac{\gamma _2^m+ q_{21} + n \gamma ^p}{{\gamma ^p}^2 (n+r_1)(n+r_2)} b_{n-1}^{(1)} \quad \text {for } n \ge 1, \end{aligned}$$
(52)
where$$\begin{aligned} r_{1,2} = \frac{\gamma _1^m + q_{12} + \gamma _2^m + q_{21}\pm \sqrt{(\gamma _2^m + q_{21} – \gamma _1^m – q_{12})^2 + 4 q_{21} q_{12}}}{2\gamma ^p} \end{aligned}$$
(53)
are the opposite numbers to the roots of the quadratic in the denominator of (51). The sequence (52) can be rewritten as$$\begin{aligned} b_n^{(1)} = \frac{\lambda ^m \left( 1+ \tau \right) _n }{\gamma _{\textrm{eff}}^m (1+r_1)_n (1+r_2)_n } \left( \frac{\lambda ^p}{\gamma ^p}\right) ^n , \quad {n\ge 1}, \end{aligned}$$
(54)
where$$\begin{aligned} \tau = \frac{\gamma _2^m + q_{21}}{\gamma ^p}, \end{aligned}$$
(55)
and \((x)_n\) represents the rising factorial. Thus, \(a_n\) can be obtained from (35) as$$\begin{aligned} a_{n} = \frac{\lambda ^p}{n \gamma ^p} b_{n-1}^{(1)}, \quad n \ge 1. \end{aligned}$$
(56)
Inserting (54) into (56) gives$$\begin{aligned} a_n = \frac{\lambda ^m r_1 r_2}{\gamma _{\textrm{eff}}^m \tau } \frac{\left( \tau \right) _n}{n (r_1)_n (r_2)_n}\left( \frac{\lambda ^p}{\gamma ^p}\right) ^n, \quad n \ge 1, \end{aligned}$$
(57)
and, similarly, inserting (54) into (50) gives$$\begin{aligned} b_n^{(2)} = \frac{q_{12} \lambda ^m \left( \tau \right) _n}{\gamma _{\textrm{eff}}^m(\gamma _2^m + q_{21}) (1+r_1)_n (1+r_2)_n} \left( \frac{\lambda ^p}{\gamma ^p}\right) ^n, \quad n \ge 1. \end{aligned}$$
(58)
Substituting (54), (57), and (58) into (37), we obtain an explicit representation of the stationary solution$$\begin{gathered} G\left( {x_{1} ,x_{2} ,z} \right) = \exp \left( {\frac{{\lambda ^{m} \lambda ^{p} }}{{\gamma _{{eff}}^{m} \gamma ^{p} }}\int\limits_{1}^{z} {{}_{2}F_{2} \left( {\begin{array}{*{20}c} {1,1 + \tau } \\ {1 + r_{1} ,1 + r_{2} } \\ \end{array} ;\frac{{\lambda ^{p} }}{{\gamma ^{p} }}(s – 1)} \right)ds} } \right. \hfill \\ \quad \quad \quad \quad \quad \quad \quad + \frac{{\lambda ^{m} \left( {x_{1} – 1} \right)}}{{\gamma _{{eff}}^{m} }}{}_{2}F_{2} \left( {\begin{array}{*{20}c} {1,1 + \tau } \\ {1 + r_{1} ,1 + r_{2} } \\ \end{array} ;\frac{{\lambda ^{p} }}{{\gamma ^{p} }}(z – 1)} \right) \hfill \\ \quad \quad \quad \quad \quad \quad \quad \left. { + \frac{{q_{{12}} \lambda ^{m} \left( {x_{2} – 1} \right)}}{{\gamma _{{eff}}^{m} \left( {\gamma _{2}^{m} + q_{{21}} } \right)}}{}_{2}F_{2} \left( {\begin{array}{*{20}c} {1,\tau } \\ {1 + r_{1} ,1 + r_{2} } \\ \end{array} ;\frac{{\lambda ^{p} }}{{\gamma ^{p} }}(z – 1)} \right)} \right), \hfill \\ \end{gathered}$$where$$_p F_q\left( \begin{array}{c} a_1, \ldots , a_p \\ b_1, \ldots , b_q, \ldots \end{array}; \tilde{z}\right) =\sum _{n=0}^{\infty } \frac{\left( a_1\right) _n \ldots \left( a_p\right) _n}{\left( b_1\right) _n \ldots \left( b_q\right) _n} \frac{\tilde{z}^n}{n!}$$is the generalised hypergeometric function [33]. Furthermore, combining (47) and (56) yields an equivalent expression$$\langle n \rangle = \frac{\lambda ^p \lambda ^m (\gamma _2^m + q_{21})}{\gamma ^p((\gamma _1^m + q_{12})(\gamma _2^m + q_{21}) – q_{12}q_{21})} = \frac{\lambda ^p \lambda ^m}{\gamma ^p \gamma _{\textrm{eff}}^m}$$for the protein mean given by (43) in terms of the model parameters. Likewise, substituting (56) and (51) into (44) and simplifying gives$$\begin{aligned} \textrm{F} = 1 + \frac{b_n^{(1)}}{b_n^{(0)}} = 1 + \frac{\lambda ^p}{\gamma ^p+\gamma _1^m+ \frac{q_{12}(\gamma ^p+\gamma _2^m)}{\gamma ^p+\gamma _2^m+q_{21}}} \end{aligned}$$
(59)
for the steady-state protein Fano factor as function of the model parameters.Multiphasic mRNA lifetimeIn this section, we consider that mRNA molecules posses \(K > 2\) stages of their lifetime, where the transition rates correspond to the ageing of an mRNA molecule. The chemical reaction system for this multiphasic model was given in (4). We note that kinetic proof reading cascades can be an interesting application of our multiphasic model [39].By (4), there are K stages of an mRNA’s molecule lifetime, each of which lasts \(1 / K \gamma _{\textrm{eff}}^m\) on average. The total mRNA lifetime is then \(1 / \gamma _{\textrm{eff}}^m\); \(\gamma _{\textrm{eff}}^m\) is thereby interpreted as the effective mRNA decay rate. The multiphasic mRNA decay in K steps leads to an Erlang-distributed lifetime with mean \(1 / \gamma _{\textrm{eff}}^m\) and variance \(1 / (\gamma _{\textrm{eff}}^m)^2\), whereas the lifetime distribution is exponential in the standard model (1).The multiphasic model (4) can be obtained by making the following choices in the general model statement (2):$$\lambda _{i}^{m} = \left\{ {\begin{array}{*{20}l} {\lambda ^{m} } \hfill & {{\text{ for }}i = 1,} \hfill \\ 0 \hfill & {{\text{ for }}i \ne 1,} \hfill \\ \end{array} } \right.$$and$$\begin{aligned} \gamma _i^m = {\left\{ \begin{array}{ll} K \gamma _{\textrm{eff}}^m & \text { if } i=K,\\ 0 & \text { otherwise}. \end{array}\right. } \end{aligned}$$The transition matrix \(\textbf{Q}\) (32) for the multiphasic model takes the form of$$\begin{aligned} \textbf{Q} = K \gamma _{\textrm{eff}}^m \begin{pmatrix} -1 & 1 & & & \\ & -1 & 1 & & \\ & & \ddots & \ddots & \\ & & & -1 & 1\\& & & & 0 \end{pmatrix}, \end{aligned}$$
(60)
and the matrix \(\textbf{A}\) (31) is given by$$\begin{aligned} \textbf{A} = K \gamma _{\textrm{eff}}^m \begin{pmatrix} 0 & & & & \\ & 0 & & & \\ & & \ddots & & \\ & & & 0 & \\ & & & & 1 \end{pmatrix}. \end{aligned}$$
(61)
Inserting (61) and (60) into (34), we obtain the system of recurrence equations$$\begin{aligned} K \gamma _{\textrm{eff}}^m \begin{pmatrix} 1 & & & & \\ -1 & 1 & & & \\ & -1 & \ddots & & \\ & & \ddots & 1 & \\ & & & -1 & 1 \end{pmatrix} \begin{pmatrix} b_0^{(1)}\\ b_0^{(2)}\\ \vdots \\ b_0^{(i)}\\ \vdots \\ b_0^{(K)} \end{pmatrix} = \begin{pmatrix} \lambda ^m\\ 0\\ \vdots \\ 0\\ \vdots \\ 0 \end{pmatrix}, \end{aligned}$$
(62)
from which, upon taking the i-th row of (62) and solving the recursive equations$$\begin{aligned} – K \gamma _{\textrm{eff}}^m b_0^{(i-1)} + K \gamma _{\textrm{eff}}^m b_0^{(i)} = 0, \quad \text {for } 2 \le i \le K, \end{aligned}$$where \(b_0^{(1)} = \lambda ^m / K \gamma _{\textrm{eff}}^m\), we recover$$\begin{aligned} b_0^{(i)} = \frac{\lambda ^m}{K \gamma _{\textrm{eff}}^m}. \end{aligned}$$
(63)
Formula (63) gives the mean of mRNA molecule in the i-th state of its lifetime. Note that the matrix \(\textbf{B}\) (33) takes the form of \(\textbf{B} = \lambda ^p \textbf{I}\), where \(\textbf{I}\) is the identity matrix.Having found the first moments (i.e. means) (63), we then determine the second moments. Taking \(n=1\) in (30), we have$$\begin{aligned} \begin{pmatrix} K \gamma _{\textrm{eff}}^m + \gamma ^p & & & & \\ -K \gamma _{\textrm{eff}}^m & K \gamma _{\textrm{eff}}^m + \gamma ^p & & & \\ & -K \gamma _{\textrm{eff}}^m & \ddots & & \\ & & \ddots & K \gamma _{\textrm{eff}}^m + \gamma ^p & \\ & & & -K \gamma _{\textrm{eff}}^m & K \gamma _{\textrm{eff}}^m + \gamma ^p \end{pmatrix} \begin{pmatrix} b_1^{(1)}\\ b_1^{(2)}\\ \vdots \\ b_1^{(i)}\\ \vdots \\ b_1^{(K)} \end{pmatrix} = \frac{\lambda ^p \lambda ^m}{K \gamma _{\textrm{eff}}^m} \begin{pmatrix} 1\\ 1\\ \vdots \\ 1\\ \vdots \\ 1 \end{pmatrix}, \end{aligned}$$
(64)
from which we obtain the first term of the sequence \(b_1^{(i)}\) as$$\begin{aligned} b_1^{(1)} := u = \frac{\lambda ^p \lambda ^m}{K \gamma _{\textrm{eff}}^m (K \gamma _{\textrm{eff}}^m + \gamma ^p)}. \end{aligned}$$
(65)
Equation (64) implies that$$\begin{aligned} – K \gamma _{\textrm{eff}}^m b_1^{(i-1)} + (K \gamma _{\textrm{eff}}^m + \gamma ^p) b_1^{(i)} = \frac{\lambda ^p \lambda ^m}{K \gamma _{\textrm{eff}}^m}, \quad \text {for } 2 \le i \le K, \end{aligned}$$which can equivalently be rewritten as$$\begin{aligned} b_1^{(i)} = u + v b_1^{(i-1)}, \quad 2 \le i \le K, \end{aligned}$$
(66)
where we set$$\begin{aligned} v = \frac{K \gamma _{\textrm{eff}}^m}{K \gamma _{\textrm{eff}}^m + \gamma ^p} \end{aligned}$$
(67)
for simplicity. Combining (66) and (65), we obtain$$\begin{aligned} b_1^{(i)} = \frac{u}{1-v} + v^{i-1} \left( u – \frac{u}{1-v}\right) , \quad 1 \le i \le K, \end{aligned}$$
(68)
from which all the elements of \(b_1^{(i)}\) (thereby the second moments) can be iteratively obtained. It is worth noting that one can derive higher moments using formula (30), but we limit our study to the first two moments.Next, we focus on calculating the first two terms of the sequence \(a_n\) (35). Setting \(n=1,2\) in (35) and inserting (63) and (68) into the resulting equations, respectively, we get$$\begin{aligned} a_1 = \frac{\lambda ^p \lambda ^m}{\gamma ^p \gamma _{\textrm{eff}}^m} \quad \text { and} \quad a_2 = \frac{\lambda ^p u\left( K+v \left( -1 – K + v^K \right) \right) }{2 \gamma ^p(1-v)^2}. \end{aligned}$$
(69)
Having found the first two terms of \(a_n\), we are now ready to calculate the Fano factor. Inserting (69) into (44), and substituting (65) and (67) into the resulting expression yields$$\begin{aligned} \mathrm {F_{m}} = 1 + \frac{\lambda ^p}{\gamma ^p} \left( 1 + \frac{\gamma _{\textrm{eff}}^m}{\gamma ^p} \left( -1 + \left( \frac{K \gamma _{\textrm{eff}}^m}{K \gamma _{\textrm{eff}}^m + \gamma ^p} \right) ^K \right) \right) , \end{aligned}$$where \(\mathrm {F_{m}}\) stands for the multiphasic Fano factor.
Translation regulation by RNA stem-loops can reduce gene expression noise | BMC Bioinformatics
