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 posttranscriptional 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 productcoding 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 selfregulation 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.
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 noncoding 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 SollnerWebb and Mougey51].
Our topic of interest is microRNA (miRNA) – short (approximately 21 nucleotides) noncoding RNA, which is engaged in transcriptional and posttranscriptional 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, GriffithsJones, 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 noncoding regions in between proteincoding 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 cotranscribed in a long transcript precursor (premRNA) 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].
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, BarJoseph, Gerber, Hannett, Harbison, Thompson, Simon and Zeitlinger32] and E. coli [Reference Khammash28, Reference ShenOrr, 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 OchabMarcinek25, 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 driftjump process, gradually increasing the complexity of the studied models. In section 2, we will start with the wellstudied model of NFB, where only a single species is involved; we will derive the stationary distribution of the protein and subsequently find its steadystate 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 steadystate 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 piecewisedeterministic 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
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 steadystate solution of the singlevalued ChapmanKolmogorov 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, OteroMuras and Vázquez38]
where as $J_{\text{in}}$ we denoted the probability influx, which is given by
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 nonincreasing function, which depends on concentration right before a burst. In particular, we define the response function $\lambda (x)$ as a step function:
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)$ :
The general solution of the above equation is [Reference Bokes and Singh7]
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,
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 realvalued exponential integral, that is, $E_1(a) = \int _a^\infty t^{1}e^{t} \mathrm{d}{t}$ .
The mean steadystate 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
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^{s1} 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:
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 twodimensional vector $(x, y)$ . Deterministic decay of concentrations between bursts is described by the massaction kinetics of reactions $R_2$ and $R_3$ :
The evolution of the probability density function (pdf) $p(x,y,t)$ is given by the equation:
where
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 integrodifferential equation (3.2) can be treated as a special twodimensional case of the ChapmanKolmogorov equation, which is reviewed in Appendix A.
According to the model requirements, bursts satisfy the following:

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

the burst size $b_z = zz'$ is also independent of the time; it is drawn from a probability density function $f(b_z)$ with the nonnegative 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_xb_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:
We proceed with simplification of the influx. Since an integral $\int _A f(x)\delta (xx_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
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 driftjump process (its sample paths are shown in Figure 3); an integrodifferential equation for $p(x,y,t)$ is the following:
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 abovementioned relation between $\theta (\cdot )$ and the integration limit $M$ .
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$ :
Note that the integral term in (3.5) is by definition a convolution about an axis [Reference Coon and Bernstein9]:
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
This allows us to apply a double Laplace transform to (3.5), which results in a PDE:
which is, in comparison with (3.5), more suitable for further analysis, as it is shown in the next sections.
4. Meanfield 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)$ :
In particular, applying this to (3.7), we obtain the system of moment equations:
As is typical for kinetics systems with bimolecular reactions, the moment equations are not closed [Reference Singh and Hespanha50]; that is, higherorder 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 meanfield model:
We take a look at steadystate mean concentrations of $x$ and $y$ ; this leads us to the solutions:
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 timedependent 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
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 largetime 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:
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)$ :
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 higherorder moments of $X$ :
Let us start with the differentiation of (5.1) with respect to $\phi$ , then equating $\phi$ to zero:
We define the function $f(\psi ) = \partial _\phi P(0, \psi )$ ; its expanded form is the following:
Afterwards, using substitution (5.5) and replacement (5.2) for $P(0, \psi )$ , we reduce the initial PDE (5.4) to the following ODE:
the general solution of which is
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 nonnegative 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):
Any nonzero $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
As was mentioned before, the definition of $f(\psi )$ makes it possible to use the Laplace image property (5.3) to obtain the steadystate mean value of X by setting $\psi =0$ in (5.9); this yields (Figure 5)
a graph of which is shown in Figure 6 (left panel). Clearly, in a highfrequency mode, the mean value in steady state only depends on the hazard rate of interaction between the species:
which is consistent with the meanfield 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
Although we do not provide a definite proof, we expect that an explicit formula can be derived for any $k$ th moment.
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:
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 steadystate 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
to ensure that both IFFL and NFB approach the same value (Figure 6, right panel).
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 nonzero 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 finetuning of factors is possible in a system, FFL provides a much narrower interval of possible steadystate 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 twodimensional 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 steadystate 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 holinantiholin [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. APVV180308 and the VEGA grant 1/0755/22.
Competing interests
Authors declare no competing interests.
A. Driftjump 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:
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
which we will refer to as a jump kernel. Then the trajectory of the process is [Reference Schuss47]:
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 ChapmanKolmogorov differential equation [Reference Gardiner16]:
where ${J_{\text{in}}}(\mathbf{x}, t)$ is given by
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{xx'}\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{xx'}\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 nonnegativity of trajectories $\mathbf{x}(t)$ :
where $\mathbb{R}^n_+$ is the nonnegative orthant of $\mathbb{R}^n$ . We also assume that jumps are nonnegative, 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:
where the integration domain $R(\mathbf{x})$ , where $\mathbf{x} = (x_1, \ldots, x_n)$ is $n$ dimensional rectangle given as Cartesian product:
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
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:
and the sample trajectories of this model are shown in Figure B1 (left panel).
Adding degradation of mRNA will change only the deterministic dynamics of mRNA level and will not affect miRNA:
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:
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)$ :
Since for any nonhomogeneous linear ODE
a general solution is given by
we obtain the homogeneous solution of (B.2), which has a similar form as one in (5.7):
The significant difference appears in the particular solution $f_p(\psi )$ , which is given by the integral:
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]:
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:
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
After using an argument transformation [Reference Abramowitz and Stegun1]:
the particular solution takes form:
Finally, we set $\psi = 0$ and use (5.3) to obtain the mean protein concentration:
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
Under the assumption about stable mRNA ( $\gamma = 0$ ), the hypergeometric function becomes a polynomial:
C. Timedependent solution of the meanfield system
Let us solve the system (4.3), provided that the initial conditions are given by
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:
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)$ :
Now we can find an integration factor $\mu (t)$ ; in this equation, it is given by
After multiplication of equation (C.1) by $\mu (t)$ and collapsing the lefthand side by the product rule, we find
that is,
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
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:
where the arbitrary constant is
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^{s1} 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
where the arbitrary constant remains the same.
Both of the solutions above have the same form at $K=0$ :
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
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:
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 preset 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$ :
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):
Subsequently, integrating (D.3) on each interval $\Delta t_i$ yields the terms $\bar{Y}_i$ of partial sums for $M_1(Y)$ :
Analogically, we define the second moments of $x(t)$ and $y(t)$ by
After repeating the same approach, we obtain the terms $\bar{X}^2_i$ and $\bar{Y}^2_i$ in form
Given that we know the piecewisedeterministic 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.