A Fokker-Planck approach to the study of robustness in gene expression

We study several Fokker-Planck equations arising from a stochastic chemical kinetic system modeling a gene regulatory network in biology. The densities solving the Fokker-Planck equations describe the joint distribution of the messenger RNA and micro RNA content in a cell. We provide theoretical and numerical evidences that the robustness of the gene expression is increased in the presence of micro RNA. At the mathematical level, increased robustness shows in a smaller coefficient of variation of the marginal density of the messenger RNA in the presence of micro RNA. These results follow from explicit formulas for solutions. Moreover, thanks to dimensional analyses and numerical simulations we provide qualitative insight into the role of each parameter in the model. As the increase of gene expression level comes from the underlying stochasticity in the models, we eventually discuss the choice of noise in our models and its influence on our results.


Introduction
This paper is concerned with a mathematical model for a gene regulatory network involved in the regulation of DNA transcription. DNA transcription is part of the mechanism by which a sequence of the nuclear DNA is translated into the corresponding protein. The transcription is initiated by the binding of a transcription factor, which is usually another protein, onto the gene's DNA-binding domain. Once bound, the transcription factor promotes the transcription of the nuclear DNA into a messenger RNA (further denoted by mRNA), which, once released, is converted into the corresponding protein by the ribosomes. This process is subject to a high level of noise due to the large variability of the conditions that prevail in the cell and the nucleus at the moment of the transcription. Yet, a rather stable amount of the final protein is needed for the good operation of the cell. The processes that regulate noise levels and maintain cell homeostasis have been scrutinized for a long time. Recently, micro RNA's (further referred to as µRNA's) have occupied the front of the scene. These are very short RNA's which do not code for proteins. Many different sorts of µRNA's are involved in various epigenetic processes. But one of their roles seems precisely the reduction of noise level in DNA transcription. In this scenario, the µRNA's are synthesized together with the mRNA's. Then, some of the synthesized µRNA's bind to the mRNA's and de-activate them. These µRNA-bound mRNA become unavailable for protein synthesis. It has been proposed that this paradoxical mechanism which seems to reduce the efficiency of DNA transcription may indeed have a role in noise regulation (see [9,17,10] and the review [19]). The goal of the present contribution is to propose a mathematical model of the µRNA-mRNA interaction and to use it to investigate the role of µRNA's as potential noise regulators.
Specifically, in this paper, we propose a stochastic chemical kinetic model for the mRNA and µRNA content in a cell. The production of mRNA's by the transcription factor and their inactivation through µRNA binding are taken into account. More precisely, our model is a simplified version of the circuit used in [25, Fig. 2A and 2A']. We consider a ligand involved in the production of both an mRNA and a µRNA, the µRNA having the possibility to bind to the mRNA and deactivate it. By contrast to [25], we disregard the way the ligand is produced and consider that the ligand is such that there is a constant production rate of both mRNA and µRNA. A second difference to [25] is that we disregard the transcription step of the mRNA into proteins. While [25] proposes to model the µRNA as acting on the transcription rate of the mRNA into proteins, we assume that the µRNA directly influences the number of mRNA available for transcription. Therefore, we directly relate the gene expression level to the number of µRNA-free mRNA also referred to as the number of unbound mRNA. In order to model the stochastic variability in the production of the RNAs, a multiplicative noise is added to the production rate at all time. From the resulting system of stochastic differential equations, we introduce the joint probability density for mRNA and µRNA which solves a deterministic Fokker-Planck equation. The mathematical object of interest is the stationary density solving the Fokker-Planck equation and more precisely the marginal density of the mRNA. The coefficient of variation (also called cell-to-cell variation) of this mRNA density, which is its standard deviation divided by its the expectation, is often considered as the relevant criteria for measuring the robustness of gene expression (see for instance [25]).
Our main goal in this contribution is to provide theoretical and numerical evidences that the robustness of the gene expression is increased in the presence of µRNA. At the theoretical level we derive a number of analytical formulas either for particular subsets of parameters of the model or under some time-scale separation hypotheses. From these formulas we can easily compute the cell to cell variation numerically and verify the increased robustness of gene expression when binding with microRNA happens in the model. For general sets of parameters, the solution cannot be computed analytically. However we can prove well-posedness of the model and solve the PDE with a specifically designed numerical scheme. From the approximate solution, we compute the coefficient of variation and verify the hypothesis of increased gene expression.
Another classical approach to the study of noise in gene regulatory networks is through the chemical master equation [27] which is solved numerically by means of Gillespie's algorithm [18], see e.g. [12,25]. Here, we use a stochastic chemical kinetic model through its associated Kolmogorov-Fokker-Planck equation. Chemical kinetics is a good approximation of the chemical master equation when the number of copies of each molecule is large. This is not the case in a cell where sometimes as few as a 100 copies of some molecules are available. Specifically, including a stochastic term in the chemical kinetic approach is a way to retain some of the randomness of the process while keeping the model complexity tractable. This ultimately leads to a Fokker-Planck model for the joint distribution of mRNAs and µRNAs. In [16], a similar chemical kinetic model is introduced with a different modelling of stochasticity. The effect of the noise is taken into account by adding some uncertainty in the (steady) source term and the initial data. The authors are interested in looking at how this uncertainty propagates to the mRNA content and in comparing this uncertainty between situations including µRNA production or not. The uncertainty is modeled by random variables with given probability density functions. Compared to [16], the Fokker-Planck approach has the advantage that the random perturbations do not only affect the initial condition and the source term, but are present at all times and varies through time. We believe that this is coherent with how stochasticity in a cell arises through time-varying ecological or biological factors.
While Fokker-Planck equations are widely used models in mathematical biology [26], their use for the study of gene regulatory network is, up to our knowledge, scarce (see e.g. [22]). Compared to other approaches, the Fokker-Planck model enables us to derive analytical formulas for solutions in certain cases. This is particularly handy for understanding the role of each parameter in the model, calibrating them from real-world data and perform fast numerical computations. Nevertheless, in the general case, the theoretical study and the numerical simulation of the model remains challenging because of the unboundedness of the drift and diffusion coefficients. We believe that we give below all the tools for handling these difficulties, and that our simple model provide a convincing mathematical interpretation of the increase of gene expression in the presence of µRNAs.
The paper is organized as follows. In Section 2, we introduce the system of SDEs and the corresponding Fokker-Planck models. In Section 3, we discuss the well-posedness of the Fokker-Planck equations and derive analytical formulas for solutions under some simplifying hypotheses. In Section 4, we use the analytical formulas for solutions to give mathematical and numerical proofs of the decrease of cell-to-cell variation in the presence of µRNA. In Section 5, we propose a numerical scheme for solving the main Fokker-Planck model and gather further evidences confirming the hypothesis of increased gene expression from the simulations. Finally, in Section 6 we discuss the particular choice of multiplicative noise (i.e. the diffusion coefficient in the Fokker-Planck equation) in our model. In the appendix, we derive weighted Poincaré inequalities for gamma and inverse-gamma distributions which are useful in the analysis of Section 3.

Presentation of the models
In this section, we introduce three steady Fokker-Planck models whose solutions describe the distribution of unbound mRNA and µRNA within a cell. The solutions to these equations can be interpreted as the probability density functions associated with the steady states of stochastic chemical kinetic systems describing the production and destruction of mRNA and µRNA. In Section 2.1 we introduce the main model for which the consumption of RNAs is either due to external factors in the cell (transcription, etc.) or to binding between the two types of mRNA and µRNA. Then, for comparison, in Section 2.2 we introduce the same model without binding between RNAs. Finally in Section 2.3, we derive an approximate version of the first model, by considering that reactions involving µRNAs are infinitely faster than those involving mRNAs, which amplifies the binding phenomenon and mathematically allows for the derivation of analytical formulas for solutions. The latter will be made explicit in Section 3.
2.1. Dynamics of mRNA and µRNA with binding. We denote by r t the number of unbound mRNA and µ t the number of unbound microRNA of a given cell at time t. The kinetics of unbound mRNA and µRNA is then given by the following stochastic differential equations with c r , c µ , k r , k µ , σ r , σ µ being some given positive constants and c being a given non-negative constant. Let us detail the meaning of each term in the modeling. The first term of each equation models the constant production of mRNA (resp. µRNA) by the ligand at a rate c r (resp. c µ ). The second term models the binding of the µRNA to the mRNA. Unbound mRNA and µRNA are consumed by this process at the same rate. The rate increases with both the number of mRNA and µRNA. In the third term, the parameters k r and k µ are the rates of consumption of the unbound mRNA or µRNA by various decay mechanisms. The last term in both equations represents stochastic fluctuations in the production and destruction mechanisms of each species. It relies on a white noise dB t / dt where B t = (B 1 t , B 2 t ) is a bidimensional standard Brownian motion. The intensity of the stochastic noise is quantified by the parameters √ 2σ r r t and 2σ µ µ t . Such a choice of multiplicative noise ensures that r t and µ t remain non-negative along the dynamics.
In this paper we are interested in the invariant measure of (1) rather than the time dynamics described by the above SDEs. From the modelling point of view, we are considering a large number of identical cells and we assume that mRNA and µRNA numbers evolve according to (1). Then we measure the distribution of both RNAs among the population, when it has reached a steady state f ≡ f (r, µ). According to Itô's formula, the steady state should satisfy the following steady Fokker-Planck equation where the Fokker-Planck operator is given by Since we do not model the protein production stage, we assume that the observed distribution of gene expression level is proportional to the marginal distribution of mRNA, i.e.
By integration of (2) in the µ variable, ρ satisfies the equation The quantity j µ (r) is the conditional expectation of the number of µRNA within the population in the presence of r molecules of mRNA and it is given by

Dynamics of free mRNA without binding.
In the case where there is no µRNA binding, namely when c = 0, the variables r and µ are independent. Thus, the densitites of the invariant measures satisfying (2) are of the form where λ(µ) is the density of the marginal distrubution of µRNA. From the modelling point of view, it corresponds to the case where there is no feed-forward loop from µRNA. Therefore, only the dynamics on mRNA, and thus ρ 0 , is of interest in our study. It satisfies the following steady Fokker-Planck equation obtained directly from (4), It can be solved explicitly as we will discuss in Section 3.2. (2) cannot be solved explicitly. However, one can make some additional assumptions in order to get an explicit invariant measure providing some insight into the influence of the binding mechanism with µRNA. This is the purpose of the model considered hereafter. Let us assume the µRNA-mRNA binding rate, the µRNA decay and the noise on µRNA are large. Since the sink term of the µRNA equation is large, it is also natural to assume that the µRNA content is small. Mathematically, we assume the following scaling

Dynamics with binding and fast µRNA. The Fokker-Planck equation
for some small constant ε > 0. Then (r t ,μ t ) satisfies , whose corresponding steady Fokker-Planck equation for the invariant measure then writes, dropping the tilde, In the limit case where ε → 0, one may expect that at least formally, the density f ε converges to a limit density f fast satisfying As r is only a parameter in the previous equation and since the first marginal of f ε still satisfies (4) for all ε, one should have (formally)

Well-posedness of the models and analytical formulas for solutions
In this section, we show that the three previous models are well-posed. For the Fokker-Planck equations (6) and (7), we explicitely compute the solutions. They involve inverse gamma distributions.

Gamma and inverse gamma distributions.
The expressions of the gamma and inverse gamma probability densities are respectively Observe that by the change of variable y = 1/x one has which justifies the terminology. Let us also recall that the first and second moments of the inverse gamma distribution are Interestingly enough, we can show (see Appendix A for details and additional results) that inverse gamma distributions with finite first moment (α > 1) satisfy a (weighted) Poincaré inequality. The proof of the following proposition is done in Appendix A among more general considerations.
Proposition 3.1. Let α > 1 and β > 0. Then, for any function v such that the integrals make sense, one has where for any probability density ν and any function u on (0, ∞), the notation u ν denotes uν.

Explicit mRNA distribution without binding.
In the case of free mR-NAs, a solution to (6) can be computed explicitely and takes the form of an inverse gamma distribution.

Lemma 3.2. The following inverse gamma distribution
is the only classical solution to (6).
Proof. First observe that Therefore a solution of (8) must be of the form for some constants C 1 , C 2 . The first term decays like 1/r at infinity, thus the only probability density ρ 0 of this form is obtained for C 1 = 0 and C 2 = 1.
The Poincaré inequality (12) tells us that the solution of Proposition 3.2 is also the only (variational) solution in the appropriate weighted Sobolev space. Indeed, we may introduce the natural Hilbert space associated with Equation (6), Then the following uniqueness result holds. Proof. If ρ 0 andρ 0 are two solutions of (6), a straightforward consequence of (12) is that ρ 0 −ρ 0 X Another consequence of the Poincaré inequality is that if we consider the time evolution associated with the equation (6) then solutions converge exponentially fast towards the steady state ρ 0 . This justifies our focus on the stationary equations. The transient regime is very short and equilibrium is reached quickly. We can quantify the rate of convergence in terms of the parameters.

Proposition 3.4. Let ξ solve the Fokker-Planck equation
Proof. Observe that ξ − ρ 0 solves the unsteady Fokker-Planck equation, so that by multiplying the equation by (ξ − ρ 0 ) g −1 1+ kr σr , cr σr and integrating in r one gets Then by using the Poincaré inequality (12) and a Gronwall type argument, one gets the result.
3.3. Explicit mRNA distribution in the presence of fast µRNA. Now we focus on the resolution of (7). The same arguments than those establishing Lemma 3.3 show that the only function M satisfying (7) is the following inverse gamma distribution Then an application of (10) yields It remains to find ρ fast which is a probability density solving the Fokker-Planck equation Arguing as in the proof of Lemma 3.2, one observe that integrability properties force ρ fast to actually solve which yields is a normalizing constant making ρ fast a probability density function.
3.4. Well-posedness of the main Fokker-Planck model. Now we are interested in the well-posedness of (2), for which we cannot derive explicit formulas anymore. Despite the convenient functional framework introduced in Section 3.2, classical arguments from elliptic partial differential equation theory do not seem to be adaptable to the case c > 0. The main obstruction comes from an incompatibility between the natural decay of functions in the space X α,β and the rapid growth of the term c r µ when |(r, µ)| → ∞.
However, thanks to probabilistic methods detailed in [21] and focused specifically on Fokker-Planck equations, we are able to prove well-posedness of the steady Fokker-Planck equation (2). The method is based on finding a Lyapunov function for the adjoint of the Fokker-Planck operator and relies on an integral identity proved by the same authors in [20].
First of all let us specify the notion of solution. A weak solution to (2) is an integrable function f such that where the adjoint operator is given by A reformulation and combination of [21, Theorem A and Proposition 2.1] provides the following result.
is a Lyapunov function with respect to L ( i.e. it is positive on Ω and it satisfies (19) and (20)). (19) is clearly satisfied. Also, U is minimal at

Proof. First observe that condition
where it takes the value 2 and thus it is positive on Ω. Finally a direct computation yields and (20) follows.
A combination of Proposition 3.5 and Lemma 3.6 provides the following result.

Noise reduction by binding : the case of fast µRNA
In this section we focus on the comparison between the explicit distributions (13) and (16). We are providing theoretical and numerical evidences that the coefficient of variation (which is a normalized standard deviation) of (16) is less than that of (13). This quantity called cell to cell variation in the biological literature [25] characterizes the robustness of the gene expression level (the lower the better). We start by performing a rescaling in order to extract the dimensionless parameters which characterize the distributions.

Dimensional analysis.
In order to identify the parameters of importance in the models, we rescale the variable r around a characteristic valuer chosen to be This choice is natural in the sense that it corresponds to the steady state of the mRNA dynamics without binding nor stochastic effects, that is dr t = (c r −k r r) dt. It is also the expectation of ρ 0 . By rescaling r into r/r both ρ 0 and ρ fast are rescaled into dimensionless densities where C ad 0 and C ad fast are normalizing constants depending on the parameters of the model and δ, p and γ are dimensionless parameters. The first parameter only depends on constants that are independent of the dynamics of µRNAs. The two other dimensionless parameters are Let us give some insight into the biological meaning of these parameters. The parameter γ measures the relative importance of the two mechanisms of destruction of µRNAs, namely the binding with mRNAs versus the natural destruction/consumption. A large γ means that the binding effect is strong and conversely. The parameter p compares the production rate of µRNAs with that of mRNAs. Large values of p mean that there are much more µRNAs than mRNAs produced per unit of time.

Cell to cell variation (CV).
For any integrable non-negative function ν, let us denote by m k (ν) = y k ν(y) dy its k-th moment. The coefficient of variation or cell to cell variation (CV) is defined by where Exp(·) and Var(·) denote the expectation and variance. One has the following result. (22) and (23). Then one has that

Proposition 4.1. Consider the dimensionless distributions defined in
where the variance and coefficient of variation are well-defined only for δ > 1.

Finally one has the bound
which holds for all δ > 2, γ > 0 and p ≥ 0. Observe that C δ ≥ CV(ρ δ 0 ) but asymptotically Proof. The formulas for the moments follow from (10) and (11) and the limits can be taken using dominated convergence. The bound (29) is a consequence of the Prékopa-Leindler inequality (see [13] and references therein) which states that if f, g, h : R d → [0, +∞) are three functions satisfying for some λ ∈ (0, 1) and for all x, y, We . The condition (30) is then equivalent to which is satisfied as the term between brackets is always greater than 1 and the function z → z[(z + z −1 )/2] δ−1 , z > 0 is bounded from below by (1 + C 2 δ ) −1/2 , where C δ is given in (29). Then with the change of variable x = 1/(γx) in the integrals of (31), one recovers (29).
The bound (29) does not confirm at this point that CV(ρ δ,γ,p fast ) ≤ CV(ρ δ 0 ), which is the theoretical result one would hope for. However observe that C δ is fairly close to CV(ρ δ 0 ) for large δ. In the next section we provide numerical evidences that it should be possible to improve (29) to C δ = CV(ρ δ 0 ). In order to evaluate numerically the cell to cell variation we need to compute m k (ρ δ,γ,p fast ), for k = 0, 1, 2. Observe that after a change of variable these quantities can be rewritten (up to an explicit multiplicative constant depending on parameters) with f k (s) = s 2−k (1+s/(γδ)) pγδ . For the numerical computation of these integrals, we use a Gauss-Laguerre quadrature which is natural and efficient as we are dealing with functions integrated against a gamma distribution. We refer to [24] and references therein for the definition of the coefficients ω N i and quadrature points x N i . The truncation order N is chosen such that the numerical error between the approximation at order N and N + 1 is inferior to the given precision 10 −8 when p ≤ 1. For p ≥ 1, the function f k may take large values and it is harder to get the same numerical precision. In the numerical results below the mean error for the chosen sets of parameters with large values of p is around 10 −4 and the maximal error is 10 −2 . This is good enough to comment on qualitative behavior.
We plot the relative cell to cell variation CV(ρ δ,γ,p fast )/CV(ρ δ 0 ) with respect to γ and p for two different values of δ. The results are displayed on Figure 1. Then, on Figure 2, we draw the explicit distributions ρ δ,γ,p fast for various sets of parameters and compare it with ρ δ 0 . The numerical simulations of Figure 1 suggest that the bound (29) is nonoptimal. Actually, we conjecture that for all γ > 0, p ≥ 0 and δ > 1 one has . Observe that from Proposition 4.1 we know that for p ∈ (0, 1) the inequality becomes an equality when γ → 0 and γ → ∞. For p > 1 it saturates when γ → 0 and we conjecture that the coefficient of variation tends to 0 when γ tends to infinity.
From a modeling point of view, these simulations confirm that for any choice of parameter, the presence of (fast) µRNA makes the cell to cell variation decrease compared to the case without µRNA. Moreover, the qualitative behavior with respect to the parameters makes sense. Indeed we observe that whenever enough µRNA is produced (p ≥ 1), the increase of the binding phenomenon (γ → ∞) makes the cell to cell variation decay drastically.

Noise reduction by binding for the main Fokker-Planck model: numerical evidences
In this section, we compute the gene expression level of the main model described by equation (2). In this case, as there is no explicit formula for the solution, we will compute an approximation of it using a discretization of the Fokker-Planck equation. In order to compute the solution in practice, we restrict the domain to the bounded domain Ω b = [r min , r max ] × [µ min , µ max ]. Because of the truncation, we add zero-flux boundary conditions in order to keep a conservative equation. It leads to the problem

Dimensional analysis and reformulation of the equation.
In order for the numerical scheme to be more robust with respect to the size of the parameters, we start by rewriting the equation in a dimensionless version. It will also allow for comparisons with numerical experiments of the previous sections. We where the characteristic numbers of mRNA and µRNA are, as before respectively defined byr After some computations one obtains that equation (33) can be rewritten with the corresponding no-flux boundary conditions and normalization. The parameters δ = k r /σ r , p = c µ /c r and γ = cr/k µ are those of the previous sections and the two new parameters are The parameter κ compares consumption of µRNA versus that of mRNA by mechanisms which are not the binding between the two RNAs. The parameter ν compares the amplitude of the noise in the dynamics of µRNA versus that of the mRNA.

Remark 5.1.
Observe that the approximation of fast µRNA leading to the model discussed in Section 2.3 and in Section 4 in its dimensionless form amounts to taking ν = κ = 1/ε and let ε tend to 0.
As the coefficients in the advection and diffusion parts of (34) grow rapidly in r, µ and degenerate when r = 0 and µ = 0, an efficient numerical resolution of (34) is not straightforward. Moreover a desirable feature of the scheme would be a preservation of the analytically known solution corresponding to γ = 0. Because of these considerations we will discretize a reformulated version of the equation in which the underlying inverse gamma distributions explicitly appear. It will allow for a better numerical approximation when r and µ are either close to 0 or large. The reformulation is the following (37) with the associated no-flux boundary conditions and where the functions h 1 and h 2 are given by

Presentation of the numerical scheme.
We use a discretization based this reformulation (37). It is inspired by [8] and is fairly close to the so-called Chang-Cooper scheme [15]. We use a finite-volume scheme. The rectangle Ω b is discretized with a structured regular mesh of size ∆r and ∆µ in each respective direction. The centers of the control volumes are the points (r i , µ j ) with r i = ∆r/2 + i∆r and µ j = ∆µ/2 + j∆µ for i ∈ [0, N r − 1] and j ∈ [0, N µ − 1]. We also introduce the intermediate points r i+1/2 with j ∈ [−1, N µ − 1] and µ j+1/2 with i ∈ [−1, N r − 1] defined with the same formula as before. The approximation of the solution on the cell (i, j) is denoted The scheme reads, for all i ∈ [0, N r − 1] and j ∈ [0, N µ − 1], where the fluxes are given by an upwind discretization of the convection terms and a centered discretization of the diffusion term, namely One can show that the scheme (40) possesses a unique solution which is non-negative by following, for instance, the arguments of [14, Proposition 3.1]. Moreover, by construction, the scheme is exact in the case γ = 0.
Remark 5.2 (Choice of r min , r max , µ min , µ max ). Clearly f decays faster at infinity than ρ 0 since the convection term coming the binding phenomenon brings mass closer to the origin. Therefore an appropriate choice for r max and µ max , coming from the decay of the involved inverse gamma distributions, should be (say) r −δ max ≤ 10 −8 and µ −δκ/ν max ≤ 10 −8 so that the error coming from the tails of the distributions in the computation of moments is negligible. Similarly, near the origin the distributions decay very quickly to 0 (as exp(−1/·)). Therefore µ min , r min can be taken not too small without influencing the precision in the computation of moments of the solution. In practice, we chose µ min = r min = 0.05. Observe that even if nothing prevents the choice µ min = r min = 0 on paper, one experiences in practice a bad conditioning of the matrix which has to be inverted for solving the scheme.

Remark 5.3 (Implementation).
Observe that the matrix which has to be inverted in order to solve the scheme is not a square matrix because of the mass constraint (which is necessary to ensure uniqueness of the solution). In practice, in order to solve the corresponding linear system M F = B where F = (f ij ) ij and B = (0, . . . , 0, 1) ∈ R NrNµ+1 and M ∈ R (NrNµ+1)×NrNµ we use the pseudo-inverse yielding F = (M t M ) −1 M t B. Finally the use of a sparse matrix routine greatly improves the computation time. Our implementation was made using Matlab.

Numerical results.
In our test cases we use the following parameters: r min = 0.05, r max = 5, µ min = 0.05, µ max = 5, δ = 8, N r = 70, N µ = 200, κ = 1, ν = 1. On Figure 3 we compare the distribution functions f (r, µ) obtained for various sets of parameters (p, γ). We also draw the corresponding marginal ρ(r) as well as ρ 0 and ρ fast . We observe that for small values of p, ρ fast is a good approximation of ρ. For larger values it tends to amplify the phenomenon of variance reduction.
In order to confirm that the main Fokker-Planck model reduces the coefficient of variation as soon as γ > 0 we draw on Figure 4 the coefficient of variation for each distribution ρ, ρ fast relatively to that of ρ 0 for several values of p. We observe that indeed, the coefficient of variation is reduced. As in the case of fast µRNA, the decay is more pronounced when the production of µRNA is higher than that of mRNA, namely when p > 1. Interestingly enough, one also notices that the approximation ρ fast increases the reduction of CV when p > 1 and diminishes it when p < 1. A transition at the special value p = 1 was already observed on Figure 1.

Comments on the choice of noise
In this section, we discuss the influence of the type of noise in the Fokker-Planck models. Let us go back to the system of stochastic differential equations considered at the beginning and generalize it as follows   with D some given function. In the models of the previous sections we chose D(x) = x 2 . On the one hand it is natural to impose that D(x) vanishes when x → 0 in order to preserve the non-negativity of r t and µ t . On the other hand it is clear that the growth at infinity influences the tail of the equilibrium distribution which solves the corresponding Fokker-Planck equation. With a quadratic D we obtained algebraically decaying distributions. Nevertheless one may wonder if the decay of cell to cell variation due to µRNA would still be observed if D is changed so that it involves distributions with faster decay at infinity. In order to answer this question, we choose a simple enough function D so that we can still derive analytical formulas for distributions of mRNA without binding and mRNA in the presence of "fast" µRNA. Let us assume that D(r) = r .
6.1. Explicit formulas for distribution of mRNAs. In terms of modeling we may argue as in Section 2 and Section 3 in order to introduce the stationary probability distribution of mRNA without bindingρ 0 which solves It may still be solved analytically and one finds a gamma distribution In the case of fast µRNA, we may once again follow the method of Section 2 and Section 3 and introduceρ fast solving where the conditional expectation of the number of µRNA within the population with r mRNA is given bỹ where the parameters p and γ are given by (26) and (25) respectively and still quantify the intensity of the binding and the respective production of µRNA versus mRNA. The new parameter η is given by In the context of a dimensional analysis, let us mention that it would be inaccurate to compare η and δ as the σ r (and σ µ ) do not represent the same quantity depending on the choice of D. For D(r) = r 2 it has the same dimension as k r so δ = k r /σ r is the right dimensionless parameter. Here it has the same dimension as c r , which justifies the introduction of η.
6.3. Numerical computation of the cell to cell variation. The expectation, variance and coefficient of variation ofρ 0 are explicitly given by As there is no explicit formula for the coefficient of variation ofρ η,γ,p fast we evaluate it numerically as in Section 4.3. The results are displayed on Figure 5. We observe that unlike the case of a quadratic diffusion coefficient the relative cell to cell variation, i.e. the cell to cell variation in the presence of µRNA relative to cell to cell variations of the free case, is not unconditionally less than 1. For a large enough production of µRNA, it eventually decays when the binding effect is very strong. However for smaller production of µRNA or when the binding is weak, the effect is the opposite as the relative cell to cell variation is greater than 1. This is not satisfactory from the modeling point of view.
In conclusion the choice of noise is important in this model. An unconditional cell to cell variation decay in the presence of µRNA is observed for quadratic noise only. While other choices of noise may still lead to similar qualitative results, the choice D(r) = r 2 allowed us to derive explicit formulas for the approximate density ρ fast which, as numerical simulations show, is fairly close to the marginal ρ corresponding to the solution of the main Fokker-Planck model.

Concluding remarks and perspectives
In this paper, we introduced a new model describing the joint probability density of the number of mRNA and µRNA in a cell. It is based on a Fokker-Planck equation arising from a system of chemical kinetic equations for the number of two RNAs. The purpose of this simple model was to provide a mathematical explanation of the increase of gene expression in a cell by the regulatory feed-forward loop due to the presence of µRNAs with their ability to bind to and deactivate mRNAs.
Thanks to the derivation of analytical formulas for solutions, and to numerical simulations, we showed that the increased robustness of gene expression is indeed explained by the model. Along the way we provided theoretical tools for the analysis of the Fokker Planck equation at play and robust numerical methods for simulations. As the main biological hypothesis for the usefulness of µRNA in the regulation of gene expression is based on their ability to reduce external noise, we also discussed the particular choice of stochasticity in the model.
There are several perspectives to this work. A first one would be the calibration of the parameters of the model from real-world data. This would allow to quantify more precisely the amount of cell-to-cell variation reduction due to µRNA, thanks to the thorough numerical investigation done in this contribution of the effects of the parameters of the model. Besides, another perspective would be an improvement of Inequality (29) into the conjecture (32). This would bring a definitive theoretical answer to the hypothesis of increased gene expression level in the simplified model of "fast" µRNAs. One may also look into establishing a similar inequality for the general model. Finally, as the gene regulatory network in a cell is considerably more complex than the simple, yet enlightening in our opinion, dynamics proposed in this paper. A natural improvement would be the consideration of more effects in the model, such as the production of the transcription factor, or the translation of mRNA into proteins, among many others.

Appendix A. Poincaré inequalities for gamma and inverse gamma distributions
In this section we give a elementary proof of the 1D version of the Brascamp and Lieb inequality (see [13,Theorem 4.1], [11]), which is an extension of the Gaussian Poincaré inequality in the case of log-concave measures. This allows us to derive a weighted Poincaré inequality for the gamma distribution and deduce, by a change of variable, a similar functional inequality for the inverse gamma distribution. Then, for any suitably integrable function u, one has where for a density ν the notation u ν denotes uν.
Proof. Without loss of generality, as one may replace u with u− u e −V , we assume that u e −V = 0. We also assume that u is of class C 1 and compactly supported in I and one can then extend a posteriori the class of admissible function by a standard density argument. Then, using (ii) one has One concludes by integrating the right-hand side by parts and observing that boundary terms vanish again by assumption (iii).

Remark A.2. The proof is an adaptation of the original proof of the (flat) Poincaré-Wirtinger inequality by Poincaré.
Observe that for I = R and V (x) = x 2 /2, one recovers the classical Gaussian Poincaré inequality. From the Brascamp-Lieb inequality, we now infer Poincaré inequalities for gamma and inverse gamma distributions.
Proposition A.4. Let α > 1 and β > 0. Then, for any functions u, v such that the integrals make sense, one has where for a probability density ν the notation u ν denotes uν.