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 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 [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]
\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:
\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,
\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
\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:
\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$
:
\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:
\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:
where
\begin{equation*} \omega (b_x, b_y) = f(b_x)f(b_y)\delta (F(b_x)-F(b_y)), \end{equation*}
$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))$
:In our problem,
\begin{equation*} \omega (b_x, b_y) = f(b_x)\delta (b_x-b_y). \end{equation*}
$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:
\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
\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:
\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$
:
\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:
\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)$
:
\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:
\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:
\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:
\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:
\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)$
:
\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$
:
\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:
\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:
\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:
\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
\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):
\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
\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)
\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:
\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
\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:
\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:
\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]:
\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
\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:
\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
\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:
\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)$
:
\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
\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]:
\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:
\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*}
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)$
:
\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,
\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:
\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
\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.
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
\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:
\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}
\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):
\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)$
:
\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
\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.
Open access
