Hostname: page-component-77c89778f8-rkxrd Total loading time: 0 Render date: 2024-07-25T04:51:07.019Z Has data issue: false hasContentIssue false

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

Published online by Cambridge University Press:  29 May 2024

Iryna Zabaikina*
Affiliation:
Department of Applied Mathematics and Statistics, Comenius University, Bratislava, 84248, Slovakia
Pavol Bokes
Affiliation:
Department of Applied Mathematics and Statistics, Comenius University, Bratislava, 84248, Slovakia
*
Corresponding author: I. Zabaikina; Email: iryna.zabaikina@gmail.com
Rights & Permissions [Opens in a new window]

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.

Type
Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2024. Published by Cambridge University Press

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 [Reference Deng, Ramsköld, Reinius and Sandberg10]. 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 [Reference Köhler and Hurt30].

Significant fluctuations in RNA quantities are mostly inherited from the randomness of transcription and accompanying biochemical reactions [Reference Golding, Paulsson, Zawilski and Cox18, Reference Plesa, Erban and Othmer40]. The main source of stochasticity is the low copy number of product-coding regions and peculiarity of their processing [Reference Blake, Kærn, Cantor and Collins4, Reference Bokes5, Reference Paulsson39]. Moreover, transcriptional bursts are stochastic in nature [Reference Kumar, Singh and Kulkarni31, Reference Nicolas, Phillips and Naef37]. 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 [Reference Becskei and Serrano3, Reference Grönlund, Lötstedt and Elf20]. 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 [Reference Spiller, Wood, Rand and White52, Reference Suter, Molina, Gatfield, Schneider, Schibler and Naef53]. 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) [Reference He, Kokkinaki, Pant, Gallicano and Dym21, Reference Matera, Terns and Terns36, Reference Sollner-Webb and Mougey51].

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 [Reference Ferro, Bena, Grigolon and Bosia13, Reference John, Enright, Aravin, Tuschl, Sander and Marks27]. Each miRNA is complementary, that is, has perfect pairing, to a specific type of an untranslated region (UTR) of an mRNA molecule [Reference Fabian, Sonenberg and Filipowicz12, Reference Rodriguez, Griffiths-Jones, Ashurst and Bradley43]. 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 [Reference Lee, Jeon, Lee, Kim and Kim33, Reference Zhdanov59]. 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 [Reference Lin, Miller and Ying34]. In this work, we focus on the intronic miRNAs, the specific role of which is to directly regulate the production of the host gene [Reference Bosia, Osella, Baroudi, Corà and Caselle8, Reference Hinske22].

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 [Reference Reeves42]. 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 [Reference Tsang, Zhu and van Oudenaarden54]. In particular, it is present in transcriptional networks of S. cerevisiae [Reference Lee, Rinaldi, Robert, Odom, Bar-Joseph, Gerber, Hannett, Harbison, Thompson, Simon and Zeitlinger32] and E. coli [Reference Khammash28, Reference Shen-Orr, Milo, Mangan and Alon49]. This control motif was modelled as a deterministic system and proved to be perfectly adaptating [Reference Barkai and Leibler2, Reference Tyson, Chen and Novak55]; 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 [Reference Bokes, Hojcka and Singh6, Reference Hoessly, Wiuf and Xia23, Reference Shahrezaei and Swain48], yet it can be treated as a continuous variable [11, Reference Frei and Khammash14]. We will use a hybrid model, which allows to study the IFFL as a process with random production jumps and deterministic decay [Reference Holehouse, Cao and Grima24, Reference Jȩdrak and Ochab-Marcinek25, Reference Lin and Doering35]. 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 45 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 [Reference Friedman, Cai and Xie15] 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 [Reference Zabaikina, Bokes, Singh, Pang and Niehren58]. 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 [Reference Pájaro, Alonso, Otero-Muras and Vázquez38]

(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 [Reference Weiss56].

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 [Reference Bokes and Singh7]

\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 [Reference Zabaikina, Bokes and Singh57], 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 [Reference Khuri29] 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 [Reference Coon and Bernstein9]:

\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 [Reference Singh and Hespanha50]; 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 [Reference Zabaikina, Bokes and Singh57]. 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 [Reference Schinazi45], 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 [Reference Re, Corá, Taverna and Caselle41], 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 [Reference Jia, Singh and Grima26, Reference Schoenberg and Maquat46].

Acknowledgements

This work has been supported by the Slovak Research and Development Agency under the contract No. APVV-18-0308 and the VEGA grant 1/0755/22.

Competing interests

Authors declare no competing interests.

A. Drift-jump framework

In the main text, we study two distinct types of gene expression regulation using a hybrid model. Here we consider a Markov process with drift and jumps in general.

On the one hand, we assume that a deterministic motion $\mathbf{x}=\mathbf{x}(t)$ is given by a dynamical system:

(A.1) \begin{equation} \dot{\mathbf{x}}(t) = \mathbf{A}(\mathbf{x}(t), t), \ \mathbf{A}\;:\;\mathbb{R}^n \times [0, +\infty ) \to \mathbb{R}^n. \end{equation}

On the other hand, random jumps are governed by a random process $\mathbf{B}(t)$ , where a conditional probability to jump from $\mathbf{x'}$ to $\mathbf{x} = \mathbf{x'} + \mathbf{b}$ in an infinitesimal interval $\Delta t$ is

\begin{equation*} Pr\{\mathbf {B}(t) \in [\mathbf {b};\; \mathbf {b}+d\mathbf {b}] \,|\, \mathbf {x}(t) = \mathbf {x'}\} = \omega (\mathbf {b}|\mathbf {x'}, t)d\mathbf {x}, \end{equation*}

which we will refer to as a jump kernel. Then the trajectory of the process is [Reference Schuss47]:

\begin{equation*} \mathbf {x}(t+\Delta t) = \begin {cases} \mathbf {x}(t) + \mathbf {A}(\mathbf {x}, t)\Delta t, \text { w.p. } 1-\lambda (\mathbf {x}, t)\Delta t \\ \mathbf {x}(t) + \mathbf {A}(\mathbf {x}, t)\Delta t + \mathbf {B}(t), \text { w.p. } \lambda (\mathbf {x}, t)\Delta t, \end {cases}\end{equation*}

where $\lambda (\mathbf{x}, t)$ is a jump rate. Thus, particles can either proceed drifting along solutions of $\mathbf{A}$ or perform the random jump. Then a joint pdf of the vector $\mathbf{x}(t)$ is governed by the Chapman-Kolmogorov differential equation [Reference Gardiner16]:

(A.2) \begin{equation} \begin{split} \frac{\partial }{\partial t} p(\mathbf{x}, t) = -\sum _{i =1}^n \frac{\partial }{\partial x_i} \left ( A_i (\mathbf{x},t) p(\mathbf{x},t) \right ) - \lambda (\mathbf{x}, t)p(\mathbf{x}, t) +{J_{\text{in}}}(\mathbf{x}, t), \end{split} \end{equation}

where ${J_{\text{in}}}(\mathbf{x}, t)$ is given by

(A.3) \begin{equation} {J_{\text{in}}}(\mathbf{x}, t) = \int _{\mathbb{R}^n} \lambda (\mathbf{x'}, t)\omega (\mathbf{x}-\mathbf{x'}|\mathbf{x'}, t)p(\mathbf{x'}, t)\,\mathrm{d}{\mathbf{x'}}. \end{equation}

The integral term gives the probability that the process jumps from state a $\mathbf{x'}$ to the state $\mathbf{x}$ in an infinitesimal interval of time; in the notation $J_{\text{in}} (\mathbf{x}, t)$ , the subscript indicates that the term relates to the influx of probability.

It is usual that the coefficients of (A.2) do not explicitly depend on time; that is, we have $A_i (\mathbf{x},t) = A_i (\mathbf{x})$ for the drift, $\omega (\mathbf{x-x'}|\mathbf{x'}, t) = \omega (\mathbf{x}-\mathbf{x'}|\mathbf{x'})$ for the burst kernel and $\lambda (\mathbf{x}, t) = \lambda (\mathbf{x})$ for the burst frequency. The jump kernel may also not be conditioned by the outgoing state $\mathbf{x'}$ , in which case we have $\omega (\mathbf{x-x'}|\mathbf{x'}, t) = \omega (\mathbf{x}-\mathbf{x'})$ .

In the main section, the formulated equation is used in an applied problem, in which restriction is imposed on permissible values of the vector $\mathbf{x}(t)$ . Therefore, let us further assume that the dynamical system (A.1) maintains non-negativity of trajectories $\mathbf{x}(t)$ :

\begin{equation*}\forall t \geq 0, \forall \mathbf {x} = (x_1, \ldots, x_n) \in \mathbb {R}^n_+, \forall i \in {1, \ldots n}\; :\; x_i = 0 \Rightarrow A_i(\mathbf {x}, t) \geq 0,\end{equation*}

where $\mathbb{R}^n_+$ is the non-negative orthant of $\mathbb{R}^n$ . We also assume that jumps are non-negative, that is, $\omega (\mathbf{b}|\mathbf{x'}, t) = 0, \ \forall x \not \in \mathbb{R}^n_+$ ; then the influx in (A.2) becomes a definite integral:

(A.4) \begin{equation} J_{\text{in}} (\mathbf{x}, t) = \int _{R(\mathbf{x})} \lambda (\mathbf{x'}, t)\omega (\mathbf{x}-\mathbf{x'}|\mathbf{x'}, t)p(\mathbf{x'}, t)\,\mathrm{d}{\mathbf{x'}}, \end{equation}

where the integration domain $R(\mathbf{x})$ , where $\mathbf{x} = (x_1, \ldots, x_n)$ is $n$ -dimensional rectangle given as Cartesian product:

\begin{equation*}R(\mathbf {x}) = \prod _{i=1}^n [0, x_i] \subset \mathbb {R}^n.\end{equation*}

In this situation, equation (A.2) holds for $x \in \mathbb{R}_n^+$ , whereas $p(\mathbf{x}, t) = 0$ for $x \not \in \mathbb{R}_n^+$ .

It is sometimes useful to use the jump size $\mathbf{b}$ instead of outgoing state $\mathbf{x'} = \mathbf{x} - \mathbf{b}$ as the integration variable; performing this substitution in (A.4) gives an alternative expression of the influx term

(A.5) \begin{equation} J_{\text{in}} (\mathbf{x}, t) = \int _{R(\mathbf{x})} \lambda (\mathbf{x} - \mathbf{b}, t)\omega (\mathbf{b}|\mathbf{x} - \mathbf{b}, t) p(\mathbf{\mathbf{x} - \mathbf{b}}, t)\,\mathrm{d}{\mathbf{b}}. \end{equation}

We made use of it in a following way: first, we set the dimension of the process $n=1$ and obtained the master equation (2.1) for the protein distribution in the NFB model. Shortly afterwards, we set $n=2$ to derive the master equation (3.2) for the feed forward model.

B. Model expanded to the natural degradation of mRNA

In the main text, we assume that mRNA is stable and thus does not degrade. In this section, we study a more general case, where mRNA is unstable and degrades at a constant rate.

The reaction set corresponding to the new model (compared to the one in Figure 2b) is as follows:

\begin{equation*}\begin{aligned} R_1\ : &\ \emptyset \stackrel{\lambda }{\longrightarrow } B \times (X + Y), & R_3\ :\ Y \stackrel{1}{\longrightarrow } \emptyset, \\ R_2\ : &\ X + Y \overset{\delta }{\longrightarrow } Y, & R_4\ :\ X \overset{\gamma }{\longrightarrow } \emptyset, \end{aligned}\end{equation*}

and the sample trajectories of this model are shown in Figure B1 (left panel).

Figure B1. (Left) Sample trajectories of naturally degrading mRNA as $ \gamma = \delta/2$ (violet line) compared to the stable mRNA as $\gamma = 0$ (red line); miRNA concentration is not affected by and its dynamics remains the same (blue line). (Right) Violet line is (B.5), and red dashed line correspond to (5.10). Green dots are values obtained by simulation, which is constructed and tuned according to Appendix D. The parameters are the following: $\delta = 0.4, \ \lambda = 0.1, \ \beta = 0.5$ .

Adding degradation of mRNA will change only the deterministic dynamics of mRNA level and will not affect miRNA:

\begin{equation*} \dot {x} = -\delta xy -\gamma x, \quad \dot {y} = -y, \end{equation*}

where $\gamma$ is the natural degradation rate of mRNA. The evolution of the probability density function $p(x,y,t)$ is given by the following equation:

(B.1) \begin{align} \frac{\partial{p(x,y,t)}}{\partial t} = \frac{\partial }{\partial x} \left ( (\delta x y + \gamma x) 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 $J_{\text{in}}$ is identical to (3.4) since stochastic dynamics are not influenced by the $\gamma$ . The subsequent approach remains the same to one in sections 3 and 5. It allows us to convert (B.1) into ODE with respect to the function $f(\psi )$ (defined in (5.5)) with the property $f(0) = -{\mathbb{E}}(X)$ :

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

Since for any non-homogeneous linear ODE

\begin{equation*}f'(\psi ) + f(\psi ) f_0(\psi ) = g(\psi )\end{equation*}

a general solution is given by

(B.3) \begin{equation} f(\psi ) = e^{-F(\psi )} \left ( \int _0^\psi e^{F(x)} g(x) \mathrm{d}{x} + C \right ), \quad F(\psi ) = \int f_0(\psi ) \mathrm{d}{\psi }, \end{equation}

we obtain the homogeneous solution of (B.2), which has a similar form as one in (5.7):

\begin{equation*}f_H (\psi ) = C \left (\beta \psi + 1\right )^{-q} \left (\psi - \delta \right )^{ -\gamma -\delta \beta q}. \end{equation*}

The significant difference appears in the particular solution $f_p(\psi )$ , which is given by the integral:

\begin{align*} \begin{split} f_p(\psi ) & = (\beta \psi +1 )^{-P-\lambda -2} (\psi - \delta )^{-Q-1} \int _0^\psi (\beta x +1 )^{P} (x- \delta )^{Q} \mathrm{d}{x} \\ & = (\beta \psi +1 )^{-P-\lambda -2} (\psi - \delta )^{-Q-1} I(\psi ), \end{split} \end{align*}

where $P = -\frac{\lambda \beta \delta }{ \delta \beta +1} - 2$ and $Q = \frac{\lambda \beta \delta }{ \delta \beta +1} + \gamma - 1$ .

The hypergeometric function $_2F_1$ has the following representation as an integral with variable upper bound [Reference Gradshteyn and Ryzhik19]:

(B.4) \begin{equation} _2F_1(\nu, \mu, \mu +1; -bu) = \frac{\mu }{u^\mu } \int _0^u \frac{x^{\mu -1}}{(bx+1)^\nu } \mathrm{d}{x}, \quad \textrm{Re}\{\mu \} \gt 0, \ |\arg (bu+1)| \lt \pi . \end{equation}

We use substitution $y = x - \delta$ to bring the integral $I(\psi )$ to the required form. It leads to $\mu = Q+1$ and $bu+1 = \frac{\beta \psi +1 }{\delta \beta +1}$ , which satisfy conditions in (B.4). Then we evaluate the integral:

\begin{equation*}I(\psi ) = (1+\delta \beta )^P \frac { (\psi - \delta )^{Q+1}}{Q+1} \, _2F_1 \left ( -P, Q+1, Q+2;\; \frac {\beta (\delta - \psi )}{\delta \beta +1} \right ) + C_1,\end{equation*}

where appearance of a constant $C_1$ is caused by the substitution $y = x - \delta$ , which changed the interval of integration to $({-}\delta, \psi - \delta )$ . Then $I(\psi )$ was split into two integrals with middle boundary equal to zero: the first one was rewritten as the hypergeometric function as per (B.4); the second one is independent of $\psi$ and denoted by $C_1$ , and afterwards $C_1$ is absorbed by the arbitrary constant $C$ in (B.3). We obtain

\begin{align*} \begin{split} f_p(\psi ) = &- \frac{\lambda \beta (\beta \psi +1)^{-(\lambda +1)}}{\gamma +\beta \delta (\lambda +\gamma )}{\left (\frac{\beta \psi +1 }{\delta \beta +1}\right )}^{\xi +1} \\ & \times \, _2F_1 \left ( \xi + \gamma, \xi + 2, \xi + \gamma +1 ;\; \frac{\beta (\delta - \psi ) }{\beta \delta + 1 }\right ), \qquad \xi = \frac{\lambda \beta \delta }{1+\beta \delta }. \end{split} \end{align*}

After using an argument transformation [Reference Abramowitz and Stegun1]:

\begin{equation*} (1-u)^{a+b-c} \, _2F_1(a,b,c;\; u) = \, _2F_1(c-a,c-b,c;\; u), \end{equation*}

the particular solution takes form:

\begin{equation*} f_p(\psi ) = - \frac {\lambda \beta (\beta \psi +1)^{-(\lambda +1)}}{\gamma +\beta \delta (\lambda +\gamma )} \, {}_2 F_1 \left ( \gamma -1, 1, \xi + \gamma +1 ;\; \frac {\beta (\delta - \psi ) }{\beta \delta + 1 }\right ). \end{equation*}

Finally, we set $\psi = 0$ and use (5.3) to obtain the mean protein concentration:

(B.5) \begin{equation} {\mathbb{E}}(X) = \frac{\lambda \beta }{\gamma +\beta \delta (\lambda +\gamma )} \ _2F_1 \left ( 1, \gamma -1, \xi + \gamma +1 ;\; \frac{\xi }{\lambda }\right ). \end{equation}

In particular, if $\gamma = 1$ , that is, if miRNA and mRNA have the same degradation rate, then the hypergeometric function is equal to one and mean concentration of mRNA is

\begin{equation*}{\mathbb {E}}(X) = \frac {\lambda \beta }{1+\beta \delta (\lambda +1)}.\end{equation*}

Under the assumption about stable mRNA ( $\gamma = 0$ ), the hypergeometric function becomes a polynomial:

\begin{equation*}_2F_1 \left ( 1, -1, \xi +1 ;\; \frac {\xi }{\lambda } \right ) = \frac {\lambda \beta \delta +1}{\lambda \beta \delta }, \end{equation*}

and (B.5) reduces to (5.10).

C. Time-dependent solution of the mean-field system

Let us solve the system (4.3), provided that the initial conditions are given by

\begin{equation*} x(0) = x_0, \quad y(0) = y_0, \end{equation*}

where $x_0$ and $y_0$ are positive constants. The solution for $y(t)$ is derived directly from the second equation in (4.3) and the initial condition above:

\begin{equation*}y(t) = \lambda \beta \left ( \left (\frac {y_0}{\lambda \beta } - 1 \right )e^{-t} + 1 \right ).\end{equation*}

We substitute it into the first equation in (4.3), then define new parameters as $\rho = \lambda \beta \delta$ and $K=\frac{y_0}{\lambda \beta } - 1$ . We obtain a following ODE with respect to $x(t)$ :

(C.1) \begin{equation} \dot{x} + \rho x (Ke^{-t} +1) = \frac{\rho }{\delta }. \end{equation}

Now we can find an integration factor $\mu (t)$ ; in this equation, it is given by

\begin{equation*}\mu (t) = e^{\int \rho (Ke^{-t} +1) \mathrm {d}{t}} = e^{\rho t - \rho K e^{-t}}.\end{equation*}

After multiplication of equation (C.1) by $\mu (t)$ and collapsing the left-hand side by the product rule, we find

\begin{equation*}(x \mu )' = \frac {\rho }{\delta }\mu,\end{equation*}

that is,

(C.2) \begin{equation} x(t) = \frac{\rho }{\delta } \frac{1}{\mu (t)} \left ( \int _0^t \mu (\tau )\mathrm{d}{\tau } + C\right ), \end{equation}

where $C$ is an arbitrary constant. Proceeding with the integral term in (C.2), we integrate it by parts and perform a substitution $u(\tau ) = \rho K e^{-\tau }$ , obtaining

\begin{equation*}\int _0^t \mu (\tau )\mathrm {d}{\tau } = \frac {1}{\rho } \left ( e^{\rho (t-Ke^{-t})} - e^{-\rho K}\right ) + \frac {(\rho K)^\rho }{\rho } \int _{u(0)}^{u(t)} u^{-\rho } e^{-u} \mathrm {d}{u}. \end{equation*}

In case that $K\gt 0$ (i.e. the initial level of Y is greater then the production rate $\lambda \beta$ ), the concentration of X evolves in time as follows:

(C.3) \begin{equation} \begin{split} x(t) = & \ \frac{1}{\delta }\left (1-e^{-\rho (t+K-Ke^{-t})}\right ) \\ & + \frac{(\rho K)^\rho }{\delta } \frac{\gamma (1-\rho, \rho K e^{-t}) - \gamma (1-\rho, \rho K)}{e^{\rho (t-Ke^{-t})}} + C, \end{split} \end{equation}

where the arbitrary constant is

\begin{equation*}C = x_0 e^{-\rho (t+K-Ke^{-t})}.\end{equation*}

By $\gamma (\cdot, \cdot )$ is denoted the incomplete gamma function, which is a special function given by the integral $\gamma (s, x) = \int _0^x t^{s-1} e^{-t} \mathrm{d}{t}$ . The solution (C.3) is valid for $\rho \in (0, 1)$ due to the requirement for positivity of the parameter $s$ . However, by performing $\lceil \rho \rceil$ times integration by parts in (C.2) and subsequent following the approach above, the solution $x(t)$ can be expressed in terms of the special functions for any $\rho \in (k, k+1], \ k\in \mathbb{N}$ .

In case that $K \lt 0$ , we perform integration by a substitution $u(\tau )=-\rho Ke^{-\tau }$ , and the solution is

(C.4) \begin{equation} x(t) = \frac{1}{\delta }(1 - e^{-\rho (t+K-Ke^{-t})}) - \frac{({-}\rho K)^\rho }{\delta e^{\rho (t-Ke^{-t})}} \int _{u(0)}^{u(t)} u^{-\rho } e^{-u}\mathrm{d}{u} + C, \end{equation}

where the arbitrary constant remains the same.

Both of the solutions above have the same form at $K=0$ :

\begin{equation*}x(t) = \frac {1}{\delta }(1 - e^{-\rho t}) + x_0 e^{-\rho t},\end{equation*}

which can be also seen directly from (C.1).

D. Simulation

This section is aimed at verification of theoretical results obtained in section 5. This is done by constructing Algorithm 1, which simulates the process and calculates the first two moments of mRNA and miRNA concentrations.

Algorithm 1

Referring to section 3, a trajectory of X concentration in time interval $[0, T]$ has a finite number of discontinuity points, which are due to the random producing process $B(t)$ . Between bursts X decays deterministically as per dynamical system (3.1).

Let us consider the process step by step. We define $\Delta t_i$ as a waiting time after the $i$ -th burst $b_i$ until the next one; thereby the time of the next burst occurrence is $t_{i+1} = t_i+\Delta t_i$ . The quantity $\Delta t_i$ is drawn from the exponential distribution with mean $1/\lambda$ ; moreover, the size of the burst itself is also drawn from the exponential distribution with mean $\beta$ . Then, according to the theorem of average value, the mean value of the piecewise continuous trajectory $x(t)$ is given by

(D.1) \begin{equation} M_1(X) = \frac{1}{T}\sum _{i=0}^{N} \bar{X}_i, \quad \bar{X}_i = \int _{t_i}^{t_{i+1}} x(t;\; b_i) \mathrm{d}{t}, \end{equation}

where $N$ is number of observed bursts and $T$ is the total time of observation and equal to $\sum _{i=0}^N \Delta t_i$ . Also, despite that above we mainly appeal to X, the same is applied to Y.

The functions of deterministic decay $x(t)$ and $y(t)$ are obtained by solving the dynamical system (3.1). First, we obtain the solution of $y(t)$ from the second equation (since it is separated and does not involve $x(t)$ ); next, we substitute it into the equation of $x(t)$ . Due to discontinuities caused by bursts, we denote by $x_i(t)$ and $y_i(t)$ the separate solutions on the each interval $[t_i, t_{i+1})$ , which are the following:

(D.2) \begin{align} x_{i}(t) & = \hat{x}_i \, \text{exp}\left \{-\delta \hat{y}_{i}(1-e^{-(t-t_i)})\right \}, &\hat{x_i}&=x_{i-1}(t_i)+b_i, \end{align}
(D.3) \begin{align} y_{i}(t) & = \hat{y}_{i} \, e^{-(t-t_i)}, &\hat{y}_{i}&= y_{i-1}(t_i)+b_i, \end{align}

where $\hat{x}_{i}$ and $\hat{y}_{i}$ are initial conditions for the solution on each interval $\Delta t_i$ ; clearly, $\hat{x}_{0}$ and $\hat{y}_{0}$ are the pre-set concentrations at the beginning of observation at time $t=0$ .

First, let us find $\bar{X}_i$ terms by integrating (D.2) on $\Delta t_i$ :

\begin{equation*} \bar {X}_i = \hat {x}_ie^{-\delta \hat {y}_i} \int _{t_i}^{t_{i+1}} e^{ ae^{-t}}\mathrm {d}{t}, \quad a = \delta \hat {y}_{i} e^{t_i}.\end{equation*}

After performing a substitution $u(t) = ae^{-t}$ , the integral takes form of special function $\text{Ei}(x)$ [Reference Geller and Ng17]. Its evaluating on the interval $(u(t_i); u(t_{i+1}))$ gives the terms $\bar{X}_i$ in (D.1):

(D.4) \begin{equation} \bar{X}_{i} = \hat{x}_i \, \text{exp}\{ -\delta \hat{y}_i\} \big (\text{Ei}(\delta \hat{y}_i) - \text{Ei}(\delta \hat{y}_i e^{-\Delta t_i})\big ). \end{equation}

Subsequently, integrating (D.3) on each interval $\Delta t_i$ yields the terms $\bar{Y}_i$ of partial sums for $M_1(Y)$ :

(D.5) \begin{equation} \begin{split} \bar{Y}_{i} & = \hat{y}_i \, (1 - e^{-\Delta t_i}). \end{split} \end{equation}

Analogically, we define the second moments of $x(t)$ and $y(t)$ by

\begin{equation*} M_2(X) = \frac {1}{T}\sum _{i=0}^{N} \bar {X}^2_i, \quad \bar {X}^2_i = \int _{t_i}^{t_{i+1}} (x(t; b_i))^2 \mathrm {d}{t}. \end{equation*}

After repeating the same approach, we obtain the terms $\bar{X}^2_i$ and $\bar{Y}^2_i$ in form

(D.6) \begin{equation} \begin{split} \bar{X}^2_{i} & = \hat{x}^2_i \, \text{exp} \{ -2\delta \hat{y}_i\} \big (\text{Ei}(2\delta \hat{y}_i) - \text{Ei}(2\delta \hat{y}_i e^{-\Delta t_i})\big ), \\ \bar{Y}^2_{i} & = \frac{\hat{y}^2_i}{2} \, (1 - e^{-2\Delta t_i}). \end{split} \end{equation}

Given that we know the piecewise-deterministic behaviour of the trajectories, the last missing ingredient is to find an appropriate method for generation random values of $\Delta t_i$ and $b_i$ . Since they both are drawn from the exponential distribution, we choose the inverse transform method, which is simple, precise and computationally efficient [Reference Ross44].

Finally, the first and the second moments of X, obtained by executing the simulation algorithm, converge to values of corresponding analytical expressions (5.10) and (5.12). As shown in Figure 5, a satisfactory number of bursts to observe is $N \approx 10^5$ , after which the computational time significantly increases without a notable gain in the accuracy of result.

References

Abramowitz, M. & Stegun, I. A. (1968) Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. US Government printing office. Google Scholar
Barkai, N. & Leibler, S. (1997) Robustness in simple biochemical networks. Nature 387(6636), 913917.CrossRefGoogle ScholarPubMed
Becskei, A. & Serrano, L. (2000) Engineering stability in gene networks by autoregulation. Nature 405(6786), 590593.CrossRefGoogle ScholarPubMed
Blake, W. J., Kærn, M., Cantor, C. R., & Collins, J. J. (2003) Noise in eukaryotic gene expression. Nature 422(6932), 633637.Google ScholarPubMed
Bokes, P. (2022) Stationary and time-dependent molecular distributions in slow-fast feedback circuits. SIAM J. Appl. Dyn. Syst. 21(2), 903931.CrossRefGoogle Scholar
Bokes, P., Hojcka, M. & Singh, A. (2018) Buffering gene expression noise by microRNA based feedforward regulation. In International Conference on Computational Methods in Systems Biology, Springer, pp. 129145.CrossRefGoogle Scholar
Bokes, P. & Singh, A. (2017) Gene expression noise is affected differentially by feedback in burst frequency and burst size. J. Math. Biol. 74(6), 14831509.CrossRefGoogle ScholarPubMed
Bosia, C., Osella, M., Baroudi, M. E., Corà, D., & Caselle, M. (2012) Gene autoregulation via intronic microRNAs and its functions. BMC Syst. Biol. 6(1), 116.CrossRefGoogle ScholarPubMed
Coon, G. A. & Bernstein, D. L. (1953) Some properties of the double Laplace transformation. Trans. Am. Math. Soc. 74(1), 135176.CrossRefGoogle Scholar
Deng, Q., Ramsköld, D., Reinius, B., & Sandberg, R. (2014) Single-Cell RNA-seq reveals dynamic, random monoallelic gene expression in mammalian cells. Science 343(6167), 193196.CrossRefGoogle ScholarPubMed
Dey, S., Kannoly, S., Bokes, P., Dennehy, J. J., & Singh, A. (2021) Feedforward genetic circuits regulate the precision of event timing. In 2021 European Control Conference (ECC), IEEE, pp. 2127–2132.CrossRefGoogle Scholar
Fabian, M. R., Sonenberg, N. & Filipowicz, W. (2010) Regulation of mRNA translation and stability by microRNAs. Annu. Rev. Biochem. 79(1), 351379.CrossRefGoogle ScholarPubMed
Ferro, E., Bena, C. E., Grigolon, S., & Bosia, C. (2020) microRNA-mediated noise processing in cells: A fight or a game? Comput. Struct. Biotech. J. 18, 642649.CrossRefGoogle Scholar
Frei, T. & Khammash, M. (2021) Adaptive circuits in synthetic biology. Curr. Opin. Syst. Biol. 28, 100399.Google Scholar
Friedman, N., Cai, L. & Xie, X. S. (2006) Linking stochastic dynamics to population distribution: An analytical framework of gene expression. Phys. Rev. Lett. 97(16), 168302.CrossRefGoogle ScholarPubMed
Gardiner, , C. W.(1985) Handbook of Stochastic Methods, Springer, Berlin, Vol. 3.Google Scholar
Geller, M. & Ng, E. W. (1969) A table of integrals of the exponential integral. J. Res. Natl. Bur. Stand. 73(3), 191210.CrossRefGoogle Scholar
Golding, I., Paulsson, J., Zawilski, S. M., & Cox, , E. C. (2005) Real-Time kinetics of gene activity in individual bacteria. Cell 123(6), 10251036.CrossRefGoogle ScholarPubMed
Gradshteyn, I. S. & Ryzhik, I. M. (2014) Table of Integrals, Series, and Products, Academic Press.Google Scholar
Grönlund, A., Lötstedt, P. & Elf, J. (2013) Transcription factor binding kinetics constrain noise suppression via negative feedback. Nat. Commun. 4, 15.CrossRefGoogle ScholarPubMed
He, Z., Kokkinaki, M., Pant, D., Gallicano, G. I., & Dym, M. (2009) Small RNA molecules in the regulation of spermatogenesis. Reproduction 137(6), 901911.CrossRefGoogle ScholarPubMed
Hinske, L. C. G., Galante, P. A., Kuo, W. P., & Ohno-Machado, L. (2010) A potential role for intragenic miRNAs on their hosts’ interactome. BMC Genomics 11(1), 1–13.CrossRefGoogle ScholarPubMed
Hoessly, L., Wiuf, C. & Xia, P. (2021) On the sum of chemical reactions. Eur. J. Appl. Math. 34(2), 325325.Google Scholar
Holehouse, J., Cao, Z. & Grima, R. (2020) Stochastic modeling of autoregulatory genetic feedback loops: A review and comparative study. Biophys. J. 118(7), 15171525.CrossRefGoogle ScholarPubMed
Jȩdrak, J. & Ochab-Marcinek, A. (2016) Influence of gene copy number on self- regulated gene expression. J. Theor. Biol. 408, 222236.CrossRefGoogle ScholarPubMed
Jia, C., Singh, A. & Grima, R. (2022) Concentration fluctuations in growing and dividing cells: Insights into the emergence of concentration homeostasis. PLoS Comput. Biol. 18(10), e1010574.CrossRefGoogle ScholarPubMed
John, B., Enright, A. J., Aravin, A., Tuschl, T., Sander, C., & Marks, , D. S. (2004) Human microRNA targets. PLoS Biol. 2(11), e363.Google Scholar
Khammash, M. H. (2021) Perfect adaptation in biology. Cell Syst. 12(6), 509521.CrossRefGoogle ScholarPubMed
Khuri, A. I. (2004) Applications of Dirac’s delta function in statistics. Int. J. Math. Edu. Sci. Technol. 35(2), 185195.CrossRefGoogle Scholar
Köhler, A. & Hurt, E. (2007) Exporting RNA from the nucleus to the cytoplasm. Nat. Rev. Mol. Cell Biol. 8(10), 761773.CrossRefGoogle ScholarPubMed
Kumar, N., Singh, A. & Kulkarni, R. V. (2015) Transcriptional bursting in gene expression: Analytical results for general stochastic models. PLoS Comput. Biol. 11(10), e1004292.CrossRefGoogle ScholarPubMed
Lee, T. I., Rinaldi, N.J., Robert, F., Odom, D.T., Bar-Joseph, Z., Gerber, G.K., Hannett, N.M., Harbison, C.T., Thompson, C.M., Simon, I. & Zeitlinger, , J., (2002a) Transcriptional regulatory networks in Saccharomyces cerevisiae. Science 298(5594), 799804.CrossRefGoogle ScholarPubMed
Lee, Y., Jeon, K., Lee, J. T., Kim, S., & Kim, V. N. (2002b) MicroRNA maturation: Stepwise processing and subcellular localization. EMBO J. 21(17), 4663–4670.CrossRefGoogle ScholarPubMed
Lin, S. L., Miller, J. D., & Ying, S. Y. (2006) Intronic microrna (mirna). BioMed Res. Int. 2006, 1–13.Google Scholar
Lin, Y. & Doering, C. R. (2016) Gene expression dynamics with stochastic bursts: Construction and exact results for a coarse-grained model. Phys. Rev. 93(2).Google ScholarPubMed
Matera, A. G., Terns, R. M. & Terns, M. P. (2007) Non-coding RNAs: Lessons from the small nuclear and small nucleolar RNAs. Nat. Rev. Mol. Cell Biol. 8(3), 209220.CrossRefGoogle ScholarPubMed
Nicolas, D., Phillips, N. E. & Naef, F. (2017) What shapes eukaryotic transcriptional bursting? Mol. Biosyst. 13(7), 12801290.CrossRefGoogle ScholarPubMed
Pájaro, M., Alonso, A. A., Otero-Muras, I., & Vázquez, C. (2017) Stochastic modeling and numerical simulation of gene regulatory networks with protein bursting. J. Theor. Biol. 421, 5170.CrossRefGoogle ScholarPubMed
Paulsson, J. (2004) Summing up the noise in gene networks. Nature 427(6973), 415418.CrossRefGoogle ScholarPubMed
Plesa, T., Erban, R. & Othmer, H. G. (2019) Noise-induced mixing and multi-modality in reaction networks. Eur. J. Appl. Math. 30(5), 887911.CrossRefGoogle Scholar
Re, A., Corá, D., Taverna, D., & Caselle, M. (2009) Genome-Wide survey of microRNA–transcription factor feed-forward regulatory circuits in human. Mol. Biosyst. 5(8), 854867.CrossRefGoogle ScholarPubMed
Reeves, G. T. (2019) The engineering principles of combining a transcriptional in-coherent feedforward loop with negative feedback. J. Biol. Eng. 13, 111.CrossRefGoogle Scholar
Rodriguez, A., Griffiths-Jones, S., Ashurst, J. L., & Bradley, A. (2004) Identification of mammalian microRNA host genes and transcription units. Biotechfor 14(10a), 19021910.Google ScholarPubMed
Ross, S. M. (2013) Simulation, Academic Press, Chap. 5.Google ScholarPubMed
Schinazi, R. B. (2011). From Calculus to Analysis, Springer Science & Business Media Switzerland.Google Scholar
Schoenberg, D. R. & Maquat, L. E. (2012) Regulation of cytoplasmic mRNA decay. Nat. Rev. Genet. 13(4), 246259.CrossRefGoogle ScholarPubMed
Schuss, Z. (2009) Theory and Applications of Stochastic Processes: An Analytical Approach, New York, Springer Science & Business Media.Google Scholar
Shahrezaei, V. & Swain, P. S. (2008) Analytical distributions for stochastic gene expression. Proc. Natl. Acad. Sci. 105, 17256–17261.Google Scholar
Shen-Orr, S. S., Milo, R., Mangan, S., & Alon, , U.(2002) Network motifs in the transcriptional regulation network of Escherichia coli. Nat. Genet. 31(1), 6468.CrossRefGoogle Scholar
Singh, A. & Hespanha, J. P. (2010) Approximate moment dynamics for chemically reacting systems. IEEE Trans. Autom. Control 56(2), 414418.CrossRefGoogle Scholar
Sollner-Webb, B. & Mougey, E. B. (1991) News from the nucleolus: RRNA gene expression. Trends Biochem. Sci. 16, 5862.CrossRefGoogle ScholarPubMed
Spiller, D. G., Wood, C. D., Rand, D. A., & White, M. R. (2010) Measurement of single-cell dynamics. Nature 465(7299), 736745.CrossRefGoogle ScholarPubMed
Suter, D. M., Molina, N., Gatfield, D., Schneider, K., Schibler, U., & Naef, F. (2011) Mammalian genes are transcribed with widely different bursting kinetics. Science 332(6028), 472474.CrossRefGoogle ScholarPubMed
Tsang, J., Zhu, J. & van Oudenaarden, A. (2007) MicroRNA-mediated feedback and feedforward loops are recurrent network motifs in mammals. Mol. Cell 26(5), 753767.CrossRefGoogle ScholarPubMed
Tyson, J. J., Chen, K. C. & Novak, B. (2003) Sniffers, buzzers, toggles and blinkers: Dynamics of regulatory and signaling pathways in the cell. Curr. Opin. Cell Biol. 15(2), 221231.CrossRefGoogle ScholarPubMed
Weiss, J. N. (1997) The Hill equation revisited: Uses and misuses. FASEB J. 11(11), 835841.CrossRefGoogle ScholarPubMed
Zabaikina, I., Bokes, P. & Singh, A. (2020) Optimal bang–bang feedback for bursty gene expression. In 2020 European Control Conference (ECC), IEEE, pp. 277282.CrossRefGoogle Scholar
Zabaikina, I., Bokes, P. & Singh, A. (2023) Joint distribution of protein concentration and cell volume coupled by feedback in dilution. In Pang, J. & Niehren, J. (eds.), Computational Methods in Systems Biology, Springer Nature Switzerland.Google Scholar
Zhdanov, V. P. (2011) Kinetic models of gene expression including non-coding RNAs. Phys. Rep. 500(1), 142.CrossRefGoogle Scholar
Figure 0

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.

Figure 1

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.

Figure 2

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$.

Figure 3

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.

Figure 4

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$.

Figure 5

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$.

Figure 6

Figure B1. (Left) Sample trajectories of naturally degrading mRNA as $ \gamma = \delta/2$ (violet line) compared to the stable mRNA as $\gamma = 0$ (red line); miRNA concentration is not affected by and its dynamics remains the same (blue line). (Right) Violet line is (B.5), and red dashed line correspond to (5.10). Green dots are values obtained by simulation, which is constructed and tuned according to Appendix D. The parameters are the following: $\delta = 0.4, \ \lambda = 0.1, \ \beta = 0.5$.

Figure 7

Algorithm 1