Maintenance of steady-state mRNA levels by a microRNA-based feed forward loop in the presence of stochastic gene expression noise

Abstract

All vital functions of living cells rely on the production of various functional molecules through gene expression. The production periods are burst-like and stochastic due to the discrete nature of biochemical reactions. In certain contexts, the concentrations of RNA or protein require regulation to maintain a fine internal balance within the cell. Here we consider a motif of two types of RNA molecules – mRNA and an antagonistic microRNA – which are encoded by a shared coding sequence and form a feed forward loop (FFL). This control mechanism is shown to be perfectly adapting in the deterministic context. We demonstrate that the adaptation (of the mean value) becomes imperfect if production occurs in random bursts. The FFL nevertheless outperforms the benchmark feedback loop in terms of counterbalancing variations in the signal. Methodologically, we adapt a hybrid stochastic model, which has widely been used to model a single regulatory molecule, to the current case of a motif involving two species; the use of the Laplace transform thereby circumvents the problem of moment closure that arises owing to the mRNA–microRNA interaction. We expect that the approach can be applicable to other systems with nonlinear kinetics.

1. Introduction

All events in a cell life cycle are primarily governed by changes in the expression (or production) of a wide variety of functional products, such as RNA. They arise in a process called gene expression, which is decoding the information in the gene and managing the synthesis of a gene product. We focus on its first stage – transcription, during which mRNA is produced based on one of the two strands of DNA by an enzyme called RNA polymerase, which binds to the DNA exactly at the beginning of the coding sequence. The promoter determines this starting site; its accessibility to the RNA polymerase is controlled by special proteins called transcriptional factors. While the promoter is accessible, the RNA polymerase may produce not only a single but a portion of mRNA molecules, which is called a transcriptional burst [10]. Then follows a post-transcriptional stage, where mRNA gains final modifications (becomes mature); only later on it can finally leave the nucleus and be further processed. Thus, mRNA is delivered via mobile export receptors to a ribosome, where translation is performed [30].

Significant fluctuations in RNA quantities are mostly inherited from the randomness of transcription and accompanying biochemical reactions [18, 40]. The main source of stochasticity is the low copy number of product-coding regions and peculiarity of their processing [4, 5, 39]. Moreover, transcriptional bursts are stochastic in nature [31, 37]. These facts motivate the study of production stability; in particular, due to presented randomness, it is imperative to find control mechanisms, which are not sensitive to fluctuations in biochemical parameters.

One way to regulate the process is to influence directly when and how much mRNA is produced. A common example of such a mechanism is negative feedback (NFB), which is based on the self-regulation principle; that is, the protein is capable to reduce the rate of its own expression by preventing transcription [3, 20]. Such proteins are inhibiting transcriptional factors to their own gene promoter (Figure 1a); thereby, the production explicitly depends on the current concentration of the protein in the cell.

Figure 1. (a) Scheme of the protein synthesis with implemented negative autoregulation: bold arrows correspond to reactions presented to the right, and the inhibiting arrow corresponds to the regulatory function $\lambda (x)$ (the figure was created with BioRender.com); (b) formalisation of the model in terms of chemical kinetics: B is a jump process associated with bursts, and the empty set symbol indicates that reactants or products are out of our interest.

Another way was found after qualitative improvement of laboratory technologies, which gave specific information about the aforementioned biochemical reactions [52, 53]. In particular, it was discovered that some RNAs, called non-coding RNA, are released into the cell immediately after transcription, avoiding further stages of processing. Mainly, they have support functions such as transferring amino acids (tRNA), forming ribosomes (rRNA), supporting splicing (snRNA) and maintaining spermatogenesis (piRNA) [21, 36, 51].

Our topic of interest is microRNA (miRNA) – short (approximately 21 nucleotides) non-coding RNA, which is engaged in transcriptional and post-transcriptional regulation of gene expression [13, 27]. Each miRNA is complementary, that is, has perfect pairing, to a specific type of an untranslated region (UTR) of an mRNA molecule [12, 43]. MiRNA recognises the target UTR and binds to it, causing degradation or suppressing further translation of the mRNA molecule (Figure 2a). In essence, microRNAs can be classified into two fundamental types based on the genomic location of their coding sequence. Intergenic miRNAs are located in the non-coding regions in between protein-coding sequences and have their own promoter; thus, they are usually transcribed independently [33, 59]. Intragenic miRNAs, which are located in the intron (intronic miRNA) or the exon (exonic miRNA) of host genes, are co-transcribed in a long transcript precursor (pre-mRNA) and then spliced into separate molecules [34]. In this work, we focus on the intronic miRNAs, the specific role of which is to directly regulate the production of the host gene [8, 22].

Figure 2. (a) Scheme of the protein synthesis in which miRNA is included; bold arrows correspond to reactions presented to the right (the figure was created with BioRender.com); (b) the studied model in terms of chemical kinetics; the empty set symbol indicates that either reactants or product of reaction is out of our interest.

A feed forward loop (FFL) is a frequently arising control motif in gene regulatory networks [42]. Its basic principle is to anticipate external disturbance in a controlled quantity and immediately provide a counterbalancing response. A subtype of FFL, called incoherent FFL (IFFL), appears when a controlled quantity Z is amplified by X and dampened by Y indirectly through X. The miRNA–mRNA interaction is also an IFFL since an upstream transcriptional activator amplifies the level of mRNA directly and at the same time dampens it indirectly through amplification of its miRNA antagonist [54]. In particular, it is present in transcriptional networks of S. cerevisiae [32] and E. coli [28, 49]. This control motif was modelled as a deterministic system and proved to be perfectly adaptating [2, 55]; that is, IFFL retains expression at a required level regardless of the strength of an upstream signal.

Usually, studies concerning the distribution of a protein/RNA treat the amount of molecules as a discrete quantity [6, 23, 48], yet it can be treated as a continuous variable [11, 14]. We will use a hybrid model, which allows to study the IFFL as a process with random production jumps and deterministic decay [24, 25, 35]. It is thereby supposed that the lengths of periods, when the gene is not available for transcription, are distributed exponentially. The aim of our work is to derive analytical equation for raw moments of mRNA concentration distribution; then based on it, we will explore the impact of the IFFL on the mRNA level in a steady state.

We model the bursty production of mRNA in an IFFL as a Markov drift-jump process, gradually increasing the complexity of the studied models. In section 2, we will start with the well-studied model of NFB, where only a single species is involved; we will derive the stationary distribution of the protein and subsequently find its steady-state mean concentration. Then we will extend the model in section 3 so that it describes the joint probability function of two species involved in the FFL (i.e. mRNA and miRNA). We will prosecute a solution approach using the Laplace transform, which broadens the number methods of applicable to the model. Sections 4–5 will be dedicated to the steady-state mean concentrations of mRNA and miRNA both under the low noise assumption and in the bursty stochastic case. Finally, this will allow us to compare the performance of the NFB and FFL models in the context of adapting to external disturbances. Our work is focused on the accuracy of the analytical results, which will be validated using computer simulations; the generality of the approach makes it possible to use the methods of this work in further research.

2. Negative feedback model

In this section, we consider a model with negative autoregulation, which includes only the protein X itself. The mathematical model for the single protein dynamics, as well as for mRNA–miRNA interaction in the next section, is studied in terms of a piecewise-deterministic Markov process [15] in a continuous time and state space. Continuous space state is achieved by studying protein concentration–– the number of protein molecules per unit of volume – as opposed to the number of protein molecules [58]. In Figure 1b, the process is expressed as a sequence of chemical reactions, where X is produced in random bursts with frequency governed by the response function $\lambda (x)$ in the reaction $R_1$ and naturally degrades in $R_2$ .

We suppose that the burst size is drawn from exponential distribution with mean $\beta$ . The value one of the degradation rate constant corresponds to measuring time in units of the protein lifetime, which can be done without any loss of generality. According to the reaction channel $R_2$ , the deterministic motion is given by

\begin{equation*}\dot {x} = -x,\end{equation*}

that is, the protein concentration trajectory $x(t)$ is proportional to $e^{-t}$ inside any time interval which does not contain a burst event.

The (stationary) distribution of protein concentration is obtained as a steady-state solution of the single-valued Chapman-Kolmogorov equation associated to the process described above (the general form (A.2) of which is considered in Appendix A). Here the master equation reads [38]

(2.1) \begin{equation} \frac{\partial p(x,t)}{\partial t} = \frac{\partial }{\partial x}(xp(x,t)) - \lambda (x)p(x,t) +{J_{\text{in}}}(x,t), \end{equation}

where as $J_{\text{in}}$ we denoted the probability influx, which is given by

\begin{equation*} {J_{\text {in}}}(x,t) = \int _{0}^{x} \omega (b_x) \, \lambda (x-b_x) p(x-b_x,t) \mathrm {d}{b_x}. \end{equation*}

In the probability flux ${J_{\text{in}}}(x,t)$ , describing an instant growth of the protein level in bursts, $\omega (b_x)$ is a jump kernel, which defines the probability that concentration will jump up by $b_x$ to $x$ in an infinitesimal period of time, and the jump rate $\lambda (\cdot )$ is the non-increasing function, which depends on concentration right before a burst. In particular, we define the response function $\lambda (x)$ as a step function:

(2.2) \begin{equation} \lambda (x) = \begin{cases} \lambda, & x \lt K,\\ \frac{\lambda }{2}, & x=K, \\ 0, & x \gt K, \end{cases} \end{equation}

where the parameter $K$ limits the possibility of production according to the following: bursts cannot occur as long as the concentration level is higher than $K$ ; once $x$ falls beneath $K$ , bursts occur with a constant frequency $\lambda$ . This type of response function is a limiting case of the Hill function [56].

At steady state, the master equation (2.1) reduces to a Volterra integral equation. Since bursts follow the exponential distribution with mean $\beta$ , the probability of a burst exceeding the size of $b_x$ is given by the jump kernel $\omega (b_x) = e^{-b_x/\beta }$ . Then the integral kernel is given by the product of the response function $\lambda (x)$ and the jump kernel $w(b_x)$ :

\begin{equation*}xp(x) = \int _{0}^{x} \lambda (x-b_x) e^{-b_x/\beta } p(x-b_x) \mathrm {d}{b_x}.\end{equation*}

The general solution of the above equation is [7]

\begin{equation*} p(x) = C x^{-1} e^{-x/\beta } \exp \left (\int \frac {\lambda (x)}{x}\mathrm {d}{x} \right ). \end{equation*}

After we substitute the specific choice of the response function defined in (2.2), the solution becomes a function with discontinuity at $x=K$ , that is,

(2.3) \begin{equation} p(x) = \begin{cases} C K^{-\lambda } e^{-x/\beta } x^{\lambda -1}, & x \lt K,\\ C e^{-x/\beta } x^{-1}, & x \geq K, \end{cases} \end{equation}

where the normalisation constant is $C = \left (\gamma (\lambda, K/\beta )(K/\beta )^{-\lambda } + E_1(K/b) \right )^{-1},$ the function $E_1$ is a real-valued exponential integral, that is, $E_1(a) = \int _a^\infty t^{-1}e^{-t} \mathrm{d}{t}$ .

The mean steady-state concentration of X is obtained by multiplying the probability density function $p(x)$ by the state $x$ and integrating over positive values [57], leading to

(2.4) \begin{equation} {\mathbb{E}}(X) = \frac{\lambda \beta \gamma (\lambda,K/\beta )}{\gamma (\lambda, K/\beta ) + (K/\beta )^\lambda E_1(K/\beta )}, \end{equation}

where $\gamma (\cdot, \cdot )$ is an incomplete gamma function, which is a special function given by an improper integral $\gamma (s, x)= \int _0^x t^{s-1} e^{-t}\mathrm{d}{t}$ .

Equation (2.4) can be used to obtain the mean concentration in the limiting case of the infinitely high frequency:

(2.5) \begin{equation} \lim _{\lambda \to \infty }{\mathbb{E}}(X) = \frac{\beta e^{-K/\beta }}{E_1(K/\beta )}. \end{equation}

In the infinite frequency mode, the concentration of protein never falls beneath the threshold value $K$ . The value (2.5) thus corresponds to the mean level of protein in homeostasis with the strongest NFB response to concentration changes. We use it in section 6 as a benchmark for the evaluation of IFFL under similar conditions.

3. Feed forward model

We proceed with a chemical system with two species X and Y (mRNA and miRNA), which are subject to reactions given in Figure 2b. The bursty production $R_1$ is modelled by stochastic jumps, while inhibition and degradation reactions $R_2$ and $R_3$ are modelled deterministically. We assume that X has a relatively long lifetime, so we neglect its spontaneous degradation. Therefore, a molecule of X can be eliminated only via the interaction $R_2$ . While a molecule of Y survives in $R_2$ , it degrades independently due to natural decay via the reaction channel $R_3$ . As it follows from $R_3$ , all reaction rates are normalised to the lifetime of Y; also the hazard rate of the interaction between X and Y is constant and equal to $\delta$ . The model, in which mRNA is unstable and degrades at a constant rate $\gamma$ , is further studied in Appendix B.

Since two species X and Y are involved, a state of their concentrations is given by a two-dimensional vector $(x, y)$ . Deterministic decay of concentrations between bursts is described by the mass-action kinetics of reactions $R_2$ and $R_3$ :

(3.1) \begin{equation} \dot{x} = -\delta xy, \quad \dot{y} = -y. \end{equation}

The evolution of the probability density function (pdf) $p(x,y,t)$ is given by the equation:

(3.2) \begin{align} \frac{\partial p(x,y,t)}{\partial t} = \frac{\partial }{\partial x} \left ( \delta x y p(x,y,t) \right ) + \frac{\partial }{\partial y} \left (y p(x,y,t) \right ) - \lambda (x,y,t) p(x,y,t) +{J_{\text{in}}}(x,y,t), \end{align}

where

\begin{equation*} {J_{\text {in}}}(x,y,t) = \lambda (x,y,t) \int _0^x \int _0^y \omega (b_x, b_y) p(x-b_x, y-b_y, t) \mathrm {d}{b_x}\mathrm {d}{b_y}, \end{equation*}

gives the probability influx, $\lambda (x,y,t)$ is the burst rate, and $\omega (b_x, b_y)$ is the burst size distribution, which will be made explicit below. Note that we do not incorporate feedback in our model; the regulation is only of feed forward type. The integro-differential equation (3.2) can be treated as a special two-dimensional case of the Chapman-Kolmogorov equation, which is reviewed in Appendix A.

According to the model requirements, bursts satisfy the following:

• the burst rate is time-homogeneous and does not depend on the state: $\lambda (x,y,t) \equiv \lambda$ .

• the burst size $b_z = z-z'$ is also independent of the time; it is drawn from a probability density function $f(b_z)$ with the non-negative support.

• both species are increased simultaneously and by the same amount so that the jump kernel takes the form:

  \begin{equation*} \omega (b_x, b_y) = f(b_x)f(b_y)\delta (F(b_x)-F(b_y)), \end{equation*}
  where $F(b_z)$ is the cumulative distribution function corresponding to $f(b_z)$ and $\delta (\cdot )$ is the Dirac delta function. Furthermore, using a property of the Dirac delta function, having an argument that is a differentiable function [29] allows us to simplify $\delta (F(b_x) - F(b_y))$ :
  \begin{equation*} \omega (b_x, b_y) = f(b_x)\delta (b_x-b_y). \end{equation*}
  In our problem, $b_z$ is drawn from an exponential distribution with mean $\beta$ , which leads to $f(b_z) = e^{-b_z/\beta }/\beta$ .

Application of the conditions above to (3.2) leads to the influx of the form:

(3.3) \begin{equation} {J_{\text{in}}}(x,y,t) =\, \frac{\lambda }{\beta } \int _0^{x} \int _0^{y} e^{-b_y/\beta } p(x- b_x, y - b_y, t) \delta (b_x-b_y) \,db_y \, db_x. \end{equation}

We proceed with simplification of the influx. Since an integral $\int _A f(x)\delta (x-x_0)dx = f(x_0)$ , only if $x_0 \in A$ ; in the inner integral of (3.3) with respect to $b_y$ , we set the root of the delta function $x_0 = b_x$ and obtain

(3.4) \begin{equation} {J_{\text{in}}}(x,y,t) = \frac{\lambda }{\beta } \int _0^{x} \theta (y-b_x) e^{-b_x/\beta }p(x- b_x, y- b_x, t) \, db_x. \end{equation}

The appearance of the Heaviside step function $\theta (\cdot )$ can be explained as follows: if a burst size of one RNA is greater than the final value of another one, then the condition of equally sized bursts cannot be satisfied. Thereby permissible values of $b_x$ are in the interval $(0, y]$ , and at the same time, the integration interval is $(0, x]$ . After we eliminate $\theta (\cdot )$ , the upper limit becomes minimum of $x$ and $y$ .

Combining production jumps with decay drift, we obtain the dynamics of the Markovian drift-jump process (its sample paths are shown in Figure 3); an integro-differential equation for $p(x,y,t)$ is the following:

(3.5) \begin{equation} \begin{split} \frac{\partial }{\partial t}p(x,y,t) = & \frac{\partial }{\partial x}(\delta x y \ p(x,y,t)) + \frac{\partial }{\partial y}(y p(x,y,t)) - \lambda p(x,y,t) \\ & + \frac{\lambda }{\beta } \int _0^M e^{-b_x/\beta } p(x- b_x, y - b_x, t) \,db_x. \end{split} \end{equation}

where $M = \min \{x,y\}$ gives the upper bound. Note that the equation will remain the same if we change the order of integration in (3.3) and then let the root of the delta function $x_0$ be equal to $b_x$ . This is due to the evenness of the Dirac delta function and the above-mentioned relation between $\theta (\cdot )$ and the integration limit $M$ .

Figure 3. Sample trajectories of X and Y concentrations. Between stochastic bursts, X and Y decay deterministically as per (3.1). Despite the fact that X degrades faster than Y, once level of Y is close to zero, X also stops decreasing; this indicates that the only way for X to degrade is interaction with Y. In the simulation, the parameters are the following: $\lambda = 0.25$ , $\delta = 1$ , $\beta = 1$ .

Proceeding to solve (3.5): we apply a double Laplace transform to $p(x, y, t)$ with respect to variables $x$ and $y$ and obtain its image as a function of Laplace variables $\phi$ and $\psi$ :

(3.6) \begin{equation} P(\phi, \psi, t) = \int _0^\infty \int _0^\infty e^{-\phi x -\psi y}p(x,y,t) \,dxdy. \end{equation}

Note that the integral term in (3.5) is by definition a convolution about an axis [9]:

\begin{equation*}f\underset {a \ b}{*} g = \int _0^M f(x-a\nu, y-b\nu ) g(\nu ) \,d\nu \text {, where } M = \min \left \{\frac xa, \frac yb\right \},\end{equation*}

where parameters $a$ and $b$ are equal to $\beta$ , $f(\cdot )$ is joint pdf $p(\cdot )$ , and $g(\cdot )$ is an exponential function. Hence, its Laplace transform is

\begin{equation*}\mathcal {L}\left [f\underset {a \ b}{*} g \right ](\phi, \psi ) = \mathcal {L}[f](\phi, \psi ) \mathcal {L}[g](a\phi +b\psi ).\end{equation*}

This allows us to apply a double Laplace transform to (3.5), which results in a PDE:

(3.7) \begin{equation} \frac{\partial P}{\partial t} = \delta \phi \frac{\partial ^2 P}{\partial \phi \partial \psi } - \psi \frac{\partial P}{\partial \psi } - P\frac{\lambda \beta (\phi + \psi )}{\beta (\phi + \psi ) + 1}, \end{equation}

which is, in comparison with (3.5), more suitable for further analysis, as it is shown in the next sections.

4. Mean-field model

From definition (3.6), it follows that we can derive the mixed moments of random processes $X = X(t)$ and $Y = Y(t)$ using the image of $p(x, y, t)$ . It is done by differentiating (3.6) and setting $(\phi,\psi ) = (0,0)$ :

(4.1) \begin{equation} {\mathbb{E}}(X^k Y^m) = \left ({-}1\right )^{k+m} \frac{\partial ^{k+m} P}{\partial \phi ^k \partial \psi ^m}(0,0). \end{equation}

In particular, applying this to (3.7), we obtain the system of moment equations:

(4.2) \begin{equation} \begin{split} \dot{{\mathbb{E}}}(X) &= \lambda \beta - \delta{\mathbb{E}}(XY)\\ \dot{{\mathbb{E}}}(Y) & = \lambda \beta -{\mathbb{E}}(Y). \end{split} \end{equation}

As is typical for kinetics systems with bi-molecular reactions, the moment equations are not closed [50]; that is, higher-order moments appear in the equations for the means. This difficulty can be eliminated provided that the noise in the concentrations of $x$ and $y$ is low; then we can remove the expectation operators from (4.2), obtaining the mean-field model:

(4.3) \begin{equation} \begin{split} \dot{x} &= \lambda \beta - \delta xy \\ \dot{y} &= \lambda \beta - y. \end{split} \end{equation}

We take a look at steady-state mean concentrations of $x$ and $y$ ; this leads us to the solutions:

(4.4) \begin{equation} x = \frac{1}{\delta }, \quad y = \lambda \beta . \end{equation}

Interpreting transcription as input, which is characterised by production rate $\lambda \beta$ , and the concentration X as the ultimate output of the system, the first equality in (4.4) indicates that the model is perfectly adaptating: the steady state of X does not depend on the rate of transcription. Thus, the model in the small noise regime successfully balances out any disturbance of the input. This statement is reinforced by deriving a time-dependent solution (see Appendix C), which is drawn as a red curve in Figure 4. In what follows, we investigate how the inclusion of a moderate noise affects this perfect adaptation property.

5. Stationary moments

Figure 4. Response of FFL to changes in the input signal. A positive shift of the production rate leads to a significant short-term increase in the concentration of the species X. Otherwise, a negative shift leads to a short-term decrease in X. Note that the lower the production rate is, the longer time is needed to bring X to the steady level. FFL, feed forward loop.

The smallness of noise is guaranteed only in the asymptotic regime of small but frequent bursts, that is, provided that $\beta \rightarrow 0$ , $\lambda \rightarrow \infty$ , while assuming that the production rate $\lambda \beta = O(1)$ . However, these assumptions are restrictive, and we are interested in general large-time behaviour of $p(x,y,t)$ . We assume $\frac{\partial }{\partial t} p(x,y,t) = 0$ ; then the probability density function $p(x, y)$ as well as its image $P(\phi, \psi )$ do not depend on time:

(5.1) \begin{equation} \delta \phi \frac{\partial ^2 P}{\partial \phi \partial \psi } - \psi \frac{\partial P}{\partial \psi } - P \frac{\lambda \beta (\phi + \psi )}{\beta (\phi + \psi ) + 1} = 0. \end{equation}

Notice that according to (3.6), the image $P(\phi,\psi )$ at $\phi = 0$ is equal to the single Laplace transform of the marginal distribution of the species Y, which is $p_Y(y) = \int _0^\infty p(x,y)\mathrm{d}{x}$ . From this follows that we can derive $p_Y(y)$ directly by equating $\phi = 0$ in (5.1). The solution of the derived equation provides an image of $p_Y(y)$ :

(5.2) \begin{equation} P(0, \psi ) = (\beta \psi +1)^{-\lambda }. \end{equation}

The obtained function is the Laplace image of the probability density function of unregulated gene expression [57]. This is not surprising since within our model, the synthesis of Y is not affected by any regulatory mechanisms.

Adhering to this logic, the marginal distribution of $X$ can in principle be obtained from $P(\phi, 0)$ , yet it does not seem possible to obtain an analytical solution as setting $\psi = 0$ in (5.1) does not yield a closed equation. Nevertheless, it is possible to derive explicit expressions for ${\mathbb{E}}(X)$ , cf. equation (4.1), as well as higher-order moments of $X$ :

(5.3) \begin{equation} {\mathbb{E}}(X^k) = \left ({-}1\right )^{k} \frac{\partial ^{k} P}{\partial \phi ^k}(0,0). \end{equation}

Let us start with the differentiation of (5.1) with respect to $\phi$ , then equating $\phi$ to zero:

(5.4) \begin{equation} (\delta - \psi ) \frac{\partial }{\partial \psi } \Big ( \frac{\partial P}{\partial \phi }\, (0,\psi ) \Big ) - \frac{\lambda \beta \psi }{\beta \psi + 1} \frac{\partial P}{\partial \phi }\,(0,\psi ) - \frac{\lambda \beta }{(\beta \psi + 1)^2} P(0, \psi ) = 0. \end{equation}

We define the function $f(\psi ) = \partial _\phi P(0, \psi )$ ; its expanded form is the following:

(5.5) \begin{equation} f(\psi ) = - \int _0^\infty \int _0^\infty x p(x,y) e^{-\psi y} \mathrm{d}{y}\mathrm{d}{x}. \end{equation}

Afterwards, using substitution (5.5) and replacement (5.2) for $P(0, \psi )$ , we reduce the initial PDE (5.4) to the following ODE:

(5.6) \begin{equation} (\delta - \psi )f'(\psi ) - f(\psi )\frac{\lambda \beta \psi }{\beta \psi + 1} - \frac{\lambda \beta }{(\beta \psi + 1)^{\lambda + 2}} = 0, \end{equation}

the general solution of which is

(5.7) \begin{equation} \begin{split} f(\psi ) = & \, C \left ((\beta \psi + 1)^{-q} (\delta - \psi )^{-\delta \beta q} \right ) \\ & - \frac{\beta \psi + \lambda \beta \delta + 1}{\delta (\delta \beta (\lambda +1) +1)} \frac{1}{(\beta \psi +1)^{\lambda +1}}, \end{split} \end{equation}

where $q = \frac{\lambda }{1+ \beta \delta }$ .

However, we do not have initial conditions for equation (5.6) to determine the integration constant $C$ . Instead, we can see from definition (5.5) the solution $f(\psi )$ as $\psi = 0$ is the negative value of the first moment of $X$ , that is, $f(0) = -{\mathbb{E}}(X)$ ; therefore, there are indirect conditions on (5.7) due to the studied model. First, concentration of the protein must be non-negative value; next, the model is placed in a cell, whose capacity is bounded, so ${\mathbb{E}}(X)$ must be a finite value. Using the limit comparison theorem for improper integrals [45], the formulated restrictions can also be applied to the solution (5.7):

(5.8) \begin{equation} -\infty \lt f(\psi ) \leq 0. \end{equation}

Any non-zero $C$ leads to a singularity at a point $\psi = \delta$ and a contradiction with the conditions (5.8); then the only satisfying value is $C = 0$ , and the solution of (5.6) is

(5.9) \begin{equation} f(\psi ) = - \frac{\beta \psi + \lambda \beta \delta + 1}{\delta (\delta \beta (\lambda +1) +1)} \frac{1}{(\beta \psi +1)^{\lambda +1}}. \end{equation}

As was mentioned before, the definition of $f(\psi )$ makes it possible to use the Laplace image property (5.3) to obtain the steady-state mean value of X by setting $\psi =0$ in (5.9); this yields (Figure 5)

(5.10) \begin{equation} {\mathbb{E}}(X) = \frac{1}{\delta } - \frac{\beta }{\delta \beta (\lambda +1) + 1}, \end{equation}

a graph of which is shown in Figure 6 (left panel). Clearly, in a high-frequency mode, the mean value in steady state only depends on the hazard rate of interaction between the species:

(5.11) \begin{equation} \lim _{\lambda \rightarrow \infty }{\mathbb{E}}(X) = \frac 1\delta, \end{equation}

which is consistent with the mean-field result (4.4).

Let us take the second derivative of (5.1) as $\phi = 0$ and repeat the same approach. The second moment of $x$ is

(5.12) \begin{equation} \begin{split} {\mathbb{E}}(X^2) = &\, \frac{(\xi - \delta \beta )^2}{\delta ^2 \xi (\xi + \delta \beta )} + \frac{2\beta (\xi - \delta \beta )}{\delta (\lambda +2) (2\xi -1) (\xi + \delta \beta )} \\ & + \frac{\lambda \beta }{\delta (\lambda +2)}, \ \ \ \text{ where } \xi = \beta \delta (\lambda +1) + 1. \end{split} \end{equation}

Although we do not provide a definite proof, we expect that an explicit formula can be derived for any $k$ -th moment.

Figure 5. (Left) Sample trajectories of the mRNA concentrations given that initial conditions are drawn from the uniform distribution on the interval $[0, 2{\mathbb{E}}(X)]$ . In the long-term observation, the initial conditions are not influencing the dynamics of the process. (Right) Convergence of the numerically computed $M_1$ and $M_2$ (dashed lines) to the corresponding analytical values of ${\mathbb{E}}(X)$ and ${\mathbb{E}}(X^2)$ (solid lines), respectively. The parameters of the simulation are the following: $\delta = 1$ , $\lambda = 2$ , $\beta = 2$ .

The verification of the obtained results (5.10) and (5.12) was performed using a stochastic simulation approach. We constructed an algorithm, in which burst events are simulated using the inversion sampling method and between two consecutive bursts; the trajectories of species X and Y are determined by the solutions of (3.1). Let us denote the generated trajectory of X on a time interval $[0, T]$ as the function $x_T(t)$ . By the mean value theorem for integrals, we calculate the $i$ -th sample moment $M_i(X)$ as follows:

(5.13) \begin{equation} M_i (X) = \frac{1}{T} \int _0^T x_T^i (t) \mathrm{d}{t}. \end{equation}

The simulations show that the first moment $M_1(X)$ (sample mean) and the second moment $M_2(X)$ consistently converge to the obtained values of ${\mathbb{E}}(X)$ and ${\mathbb{E}}(X^2)$ (Figure 5, right panel). Deviation of the computational results are due to the influence of the initial conditions; thus, it is crucial that the simulation duration is sufficiently long and the process is stabilised (Figure 5, left panel). Further details concerning simulations of a piecewise continuous trajectory $x(t)$ and the computations of $M_i(X)$ are discussed in Appendix D.

6. Imperfect adaptation

In the previous sections, we constructed a stochastic model of gene expression that implements the FFL motif. This allowed us to obtain analytical values of mean concentration in a steady state. According to the obtained results, the IFFL exhibits a nearly perfect adaptation in the low noise regime, that is, if the mean burst size $\beta$ is relatively small. This version of FFL is well suited to support homeostasis in cell processes. Nevertheless, it is not obvious how robust it is in the presence of high noise.

One way to evaluate is by comparison with another control mechanism. We have chosen gene expression, where the production of gene product X is regulated by applying the NFB loop to its production. For this, we have obtained an explicit formula for the steady-state mean concentration of X (2.4).

Despite the different mechanics of regulation and parameter sets, comparing equations (2.5) and (5.11) allows us to set

\begin{equation*}\delta = \beta ^{-1} e^{K/\beta } E_{1}({-}K/\beta )\end{equation*}

to ensure that both IFFL and NFB approach the same value (Figure 6, right panel).

Figure 6. (Left) Influence of production rate on the steady-state mean concentration of $X$ as $\delta = 1$ ; (Right) influence of IFFL (bold lines) and negative feedback (dashed lines) on ${\mathbb{E}}(x)$ ; both models were tuned so that in the high-frequency mode ${\mathbb{E}}(x)$ approaches to the same limit (dotted lines). We let dissociation constant $K=2$ .

As also shown in the right panel of Figure 6, in the presence of high noise, the FFL becomes insensitive to changes in production rate much earlier than the NFB. Although the system’s parameters have a greater impact as the mean burst size grows, the FFL still provides a non-zero concentration even for vanishingly small burst frequencies, unlike the NFB. This can be explained by noting that the main source of decreasing mRNA concentration X in FFL is miRNA Y (which is unstable), but not natural degradation.

Based on these observations, we can conclude that if fine-tuning of factors is possible in a system, FFL provides a much narrower interval of possible steady-state concentrations; it can be an efficient way of maintaining homeostasic expression of a gene.

7. Conclusion

We described the mRNA–miRNA interaction by a Markov process with stochastic jumps and deterministic drift which incorporates the main biochemical features of the regulatory motifs. In contrast with previous uses of the framework, the current model is two-dimensional and requires, in particular, the use of the double Laplace transform. Another distinction is the appearance of the convolution about an axis, which is the nonlocal term of the master equation. This appears because of a perfect correlation between mRNA and miRNA bursts; that is, both species are produced simultaneously and in equal proportions, which notably affected the solution approach.

We have studied the steady-state moments of the species in FFL, with the main focus on their analytical expressions. Since the model has nonlinear kinetics, the moment equations are not closed. Nevertheless, the Laplace transform still provided a convenient means to derive the explicit moments. Thus, using the current framework a variety of interactions between species forming an IFFL [41], for example, pair holin-antiholin [11], can be also studied.

It is expected that the current approach can be extended to general burst size distributions or further regulation of the underlying processes. Furthermore, we believe that the simple model can help in understanding expanded descriptions that reflect realistic cell scenarios, especially ones involving cell growth and cell division [26, 46].