Modelling and Analysis of miRNA Regulation

—MiRNA regulation is involved in many important biological processes such as cell proliferation, apoptosis and metabolism. Computational predictions of miRNA targets suggest that up 30 % of human proteins coding genes may be regulated by miRNAs. This makes miRNAs one of the most abundant classes of regulatory genes in humans. In the present paper we develop a time delay model of a feedback system regulated via miRNA. The model resulted in three DDEs with three discrete time delays. Since this system is a classical case study, covering several essential features of miRNA and genetic regulatory mechanisms, general conclusions about design principles and role of time delays in the stability of gene circuits can be suggested. The basic view that time delays are a key factor in the dynamic behaviour of the system was conﬁrmed by the analytical calculations and numerical simulations.


I. INTRODUCTION
Cells are the structural and functional units of all living organisms as protein synthesis is an essential function of a cell. Protein synthesis is achieved in a two-step process. Firstly, when a protein is needed in a cell, a messenger RNA (mRNA) is created as a copy of the information from the appropriate gene (process called transcription). Then, the mRNA molecule is used as a template for the creation of the protein (process called translation). The processes of transcription and translation in eukaryotic cells are very complicated, since there are more levels of control of gene expression. Transcription and translation are processes which require a certain amount of time to complete. Thus, these two processes of gene expression induce time delays in the biochemical systems.
Smolen and coauthors in [28] describe a time delay associated to the translocation of proteins and mRNAs between 50 and 100 min. Rateitschak and Wolkenhauer in [24] describe a time delay for gene transcription between 10 and 40 min. Finally in [29] a time delay around 7 min is defined and estimated for nucleocytoplasmic shuttling.
Recently, an additional level of regulation in protein expression has been introduced following the identification of short non-protein-coding RNAs [22]. MicroRNAs (miRNAs) are small regulatory RNAs with length of 18-25 nucleotides (typically ≈ 22nt size in human), which function is to regulate the activity and stability of specific messenger RNA targets. A large body of data suggests that miRNAs are supposed to account for 1-5% of animal genes [3,12] (bind through imperfect base paining to the 3 ' UTR of their target mRNAs and interfere with translational output [16]), making them one of the most abundant gene product regulators.
MiRNAs are embedded in complex regulatory networks that involved gene activation, posttranslational regulation and protein-protein interactions. Therefore, this makes miRNAs as one of the most abundant classes of regulatory genes in animals. A key recurring function of miRNA in networks is to reinforce the gene expression of differentiated cellular states. For example, regulatory networks involving miRNAs have been described in the asymmetric differentiation of left-right neurons in C.elegans, eye and sensory organprecursor development in Drosophila, and granulocytic differentiation in human [8]. Thus, miRNA provide an ideal way to regulate rapid and localized protein synthesis. A recent study [31] classified gene regulatory networks involving miRNAs in two types of circuits: in type A circuits the transcription of the miRNAs and their targets are positively coregulated, while in type B circuits the transcription of the miRNAs and their targets is oppositely regulated by common upstream factors or processessee Fig.1. Type A circuits may be more prevalent in networks operating in homeostasis, to maintain in protein steady state and in general to reduce the noise in transcriptional/translational fluctuations. Type B circuits instead can serve as surveillance mechanisms to suppress 'leaky' transcription of target gene, and in general as positive feedback loops where a transient signal can be converted in a long-lasting cellular response such as irreversible cellular differentiation [8].
Time lags in continuous systems can produce complex dynamics and instabilities. In the recent years ordinary differential equations (ODEs) with retarded argument(s) have been widely used in modelling and analyzing regulatory systems involve many genes, factors and complex interac-tions (Gene Regulatory Networks). In these models, the gene and mRNA concentrations are timedependent variables, interactions are represented by functional and differential between variables, and retarded arguments are usually the time of transcription and translation, or time delay feedback loop. For the concentration x i (t) of the factor i at time t, the evolution in time is described by the system: where f i is a nonlinear real-valued function of the states (x j (t) , x j (t − τ )) of the vertices j ∈ N i , that interact with the factor i. An interesting and important problem appears due to the nonlinearity of f i , since an analytical solution of the system is usually not possible. In this sense, the recent molecular biological discoveries (like the miRNAs and their complex regulatory effects) clearly show the need to develop mathematical models that take into consideration the post-transcriptional regulation. Also, the models should take into account that the interactions can be complex, consisting of cyclic dependences, cooperative regulations (including Transcription Factors (TFs) that act simultaneously or that assemble in complex structures), DNA looping regulation, etc. These criteria seem difficult to satisfy with a standard modeling formalism, and answers may lie in hybrid approaches [11].
As an intuitive example of nonlinearity by characterizing the behavior of a system in terms of stimulus and response is the following: If we give the system a "kick" or "input signal" and observe a certain response to that kick, then we can ask what happens if we kick the system twice as hard. If the response is not twice as large (it might be larger or smaller), then we say the system's behavior is nonlinear [7]. What is the fuss over nonlinearity? The basic idea is that, for nonlinear systems, a small change in a parameter can lead to sudden and dramatic changes in both the qualitative and quantitative behavior of the system. For one value, the behavior might be stable, for another value only slightly different from the first, the behavior might be periodic or completely aperiodic. When time-delay, feedback-loops and nonlinearity appear combined in systems like miRNA networks, a suitable approach for investigating its dynamics is the systems bifurcation analysis. Bifurcation theory studies persistence and exchange of qualitative properties of dynamical systems under continuous perturbations. From the point of view of the dynamic systems theory, the Hopf bifurcation theorem [22] together with other elements of the bifurcation theory are basic analytical tools to investigate pathological conditions in biological systems.
A typical case is that of systems depending continuously on a single parameter. Under small variations of the parameter, the systems may lose stability, and re-stabilize near another equilibrium or a closed orbit, or a larger attractor. The type of transition can be continuous, gentle and smooth, with the new equilibrium being far from the "trivial" one [2,27]. The first case corresponds to a local bifurcation, the second one is a global bifurcation. Local bifurcations can be of two types: i) the system leaves its equilibrium state and reaches a new equilibrium state; ii) the system goes from an equilibrium state to an invariant subset generally composed of several equilibriums and curves connecting them, closed orbits, or tori, etc. For example, two component mechanisms with autocatalysis easily generate oscillations and bistability. They also exhibit a rich structure of bifurcations to more complicated behavior, for example pitchfork bifurcation, saddle-node bifurcation or Takens-Bogdanov bifurcation [26]. The most popular and elementary situation is the Andronov-Hopf bifurcation, characterized by the onset of a closed orbit, starting near the trivial equilibrium (from focus type), which is the phase portrait of a periodic solution with a period close to some fixed number [1,27]. In this work, we will restrict our investigation to the Andronov-Hopf bifurcation.
The aim of this article is to elucidate how the dynamics of the gene expression (associated with the time delays) is regulated by the miRNA. Thus, the plan of the paper is as follows: In Section 2 we formulate the mathematical model. It results in three delayed differential equations (DDEs). A qualitative analysis of the model equations is performed in Section 3. In Section 4 we explore numerically the model. Finally, in Section 5, we discuss and unify results from previous sections.

II. MODEL
There are two possible ways in which miRNAs can play a role in gene expression [33]. One is that they can accelerate the degradation of mRNAs, while the other is that they repress translation by forming silencing complexes (see Fig. 2). In the version of the model discussed here we suppose that the system acts as a feedback loop where the protein, y 2 , controls its own synthesis through the repression of mRNA, y 1 . Here we construct one equation, describing the rate of change in the concentration of miRNA, y 3 , taking into account only the first possible mechanism. We assume also (similarly to [33]) that the production of miRNA occurs at a constant rate, l, and we define k 6 as its degradation rate. According to the first mechanism (in which miRNAs are regarded as an enhancer of mRNA degradation) an extra degradation term, k 4 y 1 y 3 , to the rate equation for mRNA and miRNA (where k 4 is the degradation rate) is added. Based on these hypotheses, the system can be represented by the following mathematical model in time delayed differential equations where τ 1 is the time delay for translation, τ 2 is the time delay for transcription, τ 3 is the average time delay for degradation of miRNA, n 1 is often referred as a Hill coefficient or a cooperativity coefficient, k i (i = 1, ..., 6) and γ 1 , γ 2 are the kinetic rate constants. In the present work, according to [20,34], we consider one representative value for Hill coefficient, n 1 = 2.

III. QUALITATIVE ANALYSIS
In this section, we consider the system (2) when the Hill coefficient, n 1 , is equal to two and all constants of the model are real positive numbers. Furthermore, we investigate the bifurcation structure-particularly the Andronov-Hopf bifurcation-for system (2), using time delays τ 1 , τ 2 or τ 3 as bifurcation parameters. Here, we note that system (2) is a private case of the general system considered in [23], and for our analysis we use a specific theoretical approach obtained there.
The fixed points of the system, E = (y 1 , y 2 , y 3 ), represented by Eq. (2) can be analytically estimated and are defined by the following set of algebraic equations, including the rate constants of the model: According to Descarte's rule [10,14], the first equation in (3) has only one real positive root, which ensures that system has only one physiologically feasible fixed point.
Let us consider a small perturbation around the fixed point E of the system (2) defined as: In the case when n 1 = 2, the function k1 k2+k3y 2 2 (t−τ1) can be written as a MacLaurin series: where δ = k 2 + k 3 y 2 2 and η = 2y 2 x 2 (t − τ 1 ) + x 2 2 (t − τ 1 ). If we take only linear term from (5) and after substitution of (4) into differential equation (2) we have where The characteristic equation of (6) has the form that is where Remark 1. We note that χ = 0 is a root of (9) if and only if In the absence of delays ( Here we note that (9) is a third-degree exponential polynomial in χ with three discrete time delays. Because of the presence of three different discrete delays into (9), the analysis of the sign of the real parts of the eigenvalues is very complicated, and a direct approach cannot be considered. Thus, in our analysis we will use a method consisting of determining the stability of the steady state when firstly two delays are equal to zero, and secondly when one delay is equal to zero.
Generally, in eukariotes translation takes longer than transcription and one of the reasons is intron splicing. It has been suggested that the length of the introns is fundamental to cell timing [30]. Hence, we assume that the finite time delay τ 2 of transcription is longer than the time delay τ 1 of translation and finite time delay τ 3 of degradation of miRNA.
Setting τ 1 = τ 3 = 0 into (9), the characteristic equation becomes where It is well-known that that the stability of the equilibrium state E depends on the sign of the real parts of the roots of (12). We recall that a steady state is locally asymptotically stable if and only if all roots of (12) have negative real parts, and its stability can only be lost if these roots cross the vertical axis, that is if purely imaginary roots appear. Generally speaking, the transcendental equation (12) (for nonzero delay) cannot be solved analytically and has an indefinite number of roots. In essence, we have two main tools besides direct numerical integration; firstly, the linear stability analysis in the case of a small time delay, and secondly, the Hopf bifurcation theorem. Because from biological point of view it is known that time delay of transcription, τ 2 , in many cases is bigger than one [15,20,21,35] here we use Hopf bifurcation theorem. Thus, we let χ = m + in (m, n ∈ R), and rewrite (12) in terms of its real and imaginary parts as To find the first bifurcation point we look for purely imaginary roots χ = ±in, n ∈ R of (12), i.e. we set m = 0. Then the above two equations are reduced to −K 11 n 2 + K 31 = T 51 cosnτ 2 + T 4 nsinnτ 2 , −n 3 + K 21 n = T 4 ncosnτ 2 − T 51 sinnτ 2 , or another one It is clear that if the first bifurcation point is .., ∞. One can notice that if n is a solution of (15) (or (16)), then so is −n. Hence, in the following we only investigate for positive solutions n of (15), or (16) respectively. By squaring the two equations into system (15) and then adding them, it follows that: As E is locally asymptotically stable at τ 2 = 0, satisfies the Routh-Hurwitz conditions for stability for a cubic polynomial [9,19]. Equation (17) is a cubic in n 2 and the left-hand side is positive for very large values of n 2 and negative for n = 0 if and only if T 2 51 > K 2 31 , i.e. when Eq. (17) has at least one positive real root. Moreover, to apply the Hopf bifurcation theorem, according to [9], the following theorem in this situation applies: Theorem 1. Suppose that n b is the last positive simple root of (17). Then, in (τ b ) = in b is a simple root of (12) and m (τ 2 ) + in (τ 2 ) is differentiable with respect to τ 2 in a neighbourhood of τ 2 = τ b .
To establish an Andronov-Hopf bifurcation at τ 2 = τ b , we need to show that a pair of complex eigenvalues crosses the imaginary axis with nonzero speed, i.e. the following transversality condition d(Reχ) dt | τ =τb = 0 is satisfied. From (15) we know that τ bk corresponding to n b is Because for τ 2 = 0, equilibrium E is stable, by Butler's lemma [5], equilibrium remains stable for τ 2 < τ bk , where τ b = τ bk as k = 0. We have now to show that d(Reχ) dt | τ =τb = 0. Hence, if denote H (χ, τ 2 ) = χ 3 + K 11 χ 2 + K 21 χ + K 31 Evaluating the real part of this equation at τ 2 = τ b and setting χ = in b yield where Let θ = n 2 b , then, (17) reduces to Then, for g (θ) we have If n b is the least positive simple root of (17), then Hence, According to the Hopf bifurcation theorem [17], we define the following Theorem 2: Theorem 2. If n b is the least positive root of (17), then an Andronov-Hopf bifurcation occurs as τ 2 passes through τ b .
Corollary 2.1. When τ 2 < τ b , then the steady state E of system (2) is locally asymptotically stable.
We return to the study of (9) which with τ 2 , τ 3 > 0 has the form where τ = [τ 2 , τ 3 , τ 23 = τ 2 + τ 3 ] T denotes a point in the time delay space, i.e. τ ∈ Ω ⊂ R 3 + . Ω is the time delay space and R 3 + denotes the set of nonnegative real numbers. In order to assess the stability of E with respect to any delay τ , one should know where all χ roots of (27) lie on the complex plane. Eq. (27) has infinitely many roots on the complex plane due to the transcendental term −χτ . This makes the analytical stability assessment intractable.
where τ i ≥ 0 (i = 1, 2, ..., m) and p Obviously, in (n > 0) is a root of (27) if and only if n satisfies Separating the real and imaginary parts into (28), we obtain We square and add the equations (29), and after simplifying, we get that τ and n must be among the real solutions of We note that the right-hand side of (30) is always less than Hence if the inequality has no real solution on 0 < ω < n + , then (30) cannot be satisfied. Note that n + is the positive solution of first equation in (29), which we write as i.e.
Thus, for n + we have Rearranging terms, we write (31) as Hence, the following theorem can be formulated Remark 2. In the special case that τ 2 = τ 3 , the characteristic equation (27) becomes where T 35 = T 3 + T 5 and T 24 = T 2 + T 4 . Therefore, this case is a private one of Theorem 3.
Separating the real and imaginary parts into (36), we have where A 1 = −T 1 n 2 + T 3 cosnτ 3 + T 2 nsinnτ 3 , Adding up the squares of both equations into (37), we have where τ 5 − τ 4 = τ 3 . Clearly, the right-hand side of (38) is always less than Hence if the inequality has no real solution on 0 < λ < n + , then (39) cannot be satisfied. Similar to previous Section, we note that n + is the positive solution of first equation into (37), which is written as i.e. we obtain the same formula as (32) (and (33) respectively). Rearranging terms, we write (41) as Therefore, the following theorem can be formulated Remark 3. In the special case that τ 1 = τ 2 = τ 3 , the characteristic equation (9) becomes Hence, this case is a private one of Theorem 4.
In next Section, we illustrate numerically the existence of the behavior predicted for some values of the rate constants k i (i = 1, ..., 6) , γ 1 , γ 2 , l and the time delays τ 1 , τ 2 and τ 3 .

IV. NUMERICAL ANALYSIS
In the previous section, we proposed the analytical tools and used them for a qualitative analysis of the system, obtaining predictions about dynamics of the system, i.e. the stability and existence of periodic solutions via Andronov-Hopf bifurcation in time delay model (2). In this section, we perform a numerical analysis of model (2), based on the results previously obtained.
Some of the parameter values used in the numerical analysis were selected according to [4,20,24,29,33] in the form According to [6,32] proteins (or RNAs) degraded with a probability that depends on their structure because some of the degradation mechanisms involve multiple steps. Therefore, individual protein (or RNA) senesces through time. MiRNAs which are incorporated into the RISC complex, do not degrade with their targets but return to the cytosol a new round of target mRNA repression. It is plausible, however, that miRNA may be degraded after a few cycles of target mRNA binging [13]. Since we do not know the exact values of the average time delay for degradation of miRNA, we set τ 3 in large boundaries-from few minutes to few hours and its degradation rate k 6 ∈ [0.05, 0.3] min −1 . Our model include also one additional parameter for which no values are available and his estimations require further experimental studies. Thus, we assume here that l = 0.1 in minutes.
In order to compare the predictions with numerical results, the governing equations of the model (2), were solved numerically using MAT-LAB [18]. In Figure 3, the stable solutions for the concentration of mRNA ( y 1 ), the concentration of protein ( y 2 ) and the concentration of miRNA ( y 3 ) are shown for absence of time delay (i.e. τ 1 = τ 2 = τ 3 = 0) -see Fig. 3a, and for τ 1 = τ 3 = 0, τ 2 = 12 -see Fig. 3b. It is evident that after several physiological acceptance fluctuations, the solution of system (2) approaches a constant value (stable equilibrium state). In other words the system possesses a stable equilibrium state which corresponds to a normal miRNA regulation process. This conclusion is in accordance with the Theorem 2 (Corollary 2.1) proofed in previous Section 3.  It is seen that for larger values of τ 2 (see Fig. 4a) and τ 3 (see Fig.4b) than bifurcation one the stable limit cycle (self oscillations) occur and sustained oscillations take place. In other words, in these cases the conditions of Theorem 3 and Theorem 4 are not satisfied and the steady state of system (2) is unstable. From biological point of view, the occurrence of oscillation implies that if the average time delay for degradation of miRNA, τ 3 , can increase, then the effect of miRNA on gene expression is initially destabilizing. If the average time delay for degradation of miRNA increases further, then the effect of miRNA on gene expression can promote stability-see Fig.  4c, and again instability-see Fig. 4d. Thus gene expression follows a cyclic pattern (from stable to unstable behaviour and vice versa) as function of average time delay for degradation of miRNA. This cyclic regime is shown in Figure 5.  (2). Unstable regime (periodic oscillations) for τ 1 = τ 3 = 0, τ 2 = 20 (a); and τ 1 = 0, τ 2 = 13, τ 3 = 10 (b). Stable regime for τ 1 = 1, τ 2 = 13, τ 3 = 37 (c). Unstable regime (sustained oscillations) for τ 1 = 1, τ 2 = 13, τ 3 = 60 (d).
In Figure 5, the stable and unstable zones in the (τ 3 , τ 1 + τ 2 ) parameter space are shown. It is seen that these zones are with different size. It is also interesting to note that in unstable zones sustained oscillations with period one and quasiperiodic motion take place.

V. CONCLUSIONS
Gene expression in the human organism is posttranscriptionally regulated by miRNAs. MiRNAs are embedded in complex regulatory networks that involved gene activation, post-translational regulation and protein-protein interactions. Therefore, this makes miRNAs as one of the most abundant classes of regulatory genes in animals. In the present paper we develop a time delay model of a feedback system regulated via miRNA. Our hypothesis (according [33] is that miRNA can participate in the regulation of gene expression by accelerating the degradation of mRNA or by repressing the translation process. The model resulted in three DDEs with three discrete time delays. Since this system is a classical case study, covering several essential features of miRNA and genetic regulatory mechanisms, general conclusions about design principles and role of time delays in the stability of gene circuits can be suggested. The basic view that time delays τ 1 , τ 2 and τ 3 are a key factor in the dynamic behaviour of the system was confirmed by the analytical calculations and numerical simulations. In more details, under the assumption that an equilibrium exists, we have estimated the length of delays for which local asymptotic stability will be preserved. We have also derived criteria for which no change in stability will occur. If system (2) starts with a stable equilibrium, which for some delay(s) becomes unstable, it will likely destabilize by means of an Andronov-Hopf bifurcation leading to small amplitude periodic solutions. Our investigation of such a behaviour is devoted to the use of bifurcation analysis. Particularly, a Hopf bifurcation theorem was employed. From the viewpoint of the qualitative theory of DDEs, time delay(s) appears as a bifurcation parameter on whose values the altered (stable or unstable) behaviour of the model depends. For time delays longer than τ b , the gene expression system regulated by means of miRNA would present sustained oscillations with coupled periodic variations on the concentration of the mRNA, protein and the miRNA. In contrast, a time less than τ b would provoke damped oscillations around a stable steady state. We can say that in this case time delays have a destabilizing role. From a physiological point of view, the loss of stability might be related to emergence of new configurations in the regulatory gene circuit that could lead the system to a pathological state.
To conclude, mathematical modelling and analysis can enable to understand the mechanism underlying an observed biological process, and at the same time, provide a testable hypothesis for future studies. In this paper we investigate the main processes of the formation of a protein-transcription time (starts with splicing and polyadenylation of the initial transcript), translation time (the timespan from the emergence of mRNA) and average time delay for degradation of miRNA. Thus, our dynamical predictions can be tested in future experiments.