Original article Biomath 3 (2014), 1407231, 1–15 B f Volume ░, Number ░, 20░░ BIOMATH ISSN 1314-684X Editor–in–Chief: Roumen Anguelov B f BIOMATH h t t p : / / w w w . b i o m a t h f o r u m . o r g / b i o m a t h / i n d e x . p h p / b i o m a t h / Biomath Forum Modelling and Analysis of miRNA Regulation Svetoslav G. Nikolov Institute of Mechanics, Bulgarian Academy of Sciences, Bulgaria Department of Systems Biology and Bioinformatics, University of Rostock University of Transport, Sofia, Bulgaria e-mail: S.Nikolov@imbm.bas.bg Received: 16 May 2014, accepted: 23 July 2014, published: 15 September 2014 Abstract—MiRNA regulation is involved in many important biological processes such as cell prolif- eration, 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 dis- crete 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 confirmed by the analytical calculations and numerical simulations. Keywords-miRNA regulation; time delay; Andronov-Hopf bifurcation; numerical analysis 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 ap- propriate 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 follow- ing the identification of short non-protein-coding RNAs [22]. MicroRNAs (miRNAs) are small reg- ulatory RNAs with length of 18–25 nucleotides (typically ≈ 22nt size in human), which function Citation: Svetoslav G. Nikolov, Modelling and Analysis of miRNA Regulation, Biomath 3 (2014), 1407231, http://dx.doi.org/10.11145/j.biomath.2014.07.231 Page 1 of 15 http://www.biomathforum.org/biomath/index.php/biomath http://dx.doi.org/10.11145/j.biomath.2014.07.231 S. Nikolov, Modelling and Analysis of miRNA Regulation 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 regula- tory networks that involved gene activation, post- translational regulation and protein-protein inter- actions. 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 dif- ferentiated cellular states. For example, regulatory networks involving miRNAs have been described in the asymmetric differentiation of left-right neu- rons in C.elegans, eye and sensory organprecursor development in Drosophila, and granulocytic dif- ferentiation in human [8]. Thus, miRNA provide an ideal way to regulate rapid and localized protein synthesis. A recent study [31] classified gene reg- ulatory networks involving miRNAs in two types of circuits: in type A circuits the transcription of the miRNAs and their targets are positively co- regulated, while in type B circuits the transcription of the miRNAs and their targets is oppositely reg- ulated by common upstream factors or processes- see 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 in- volve many genes, factors and complex interac- tions (Gene Regulatory Networks). In these mod- els, the gene and mRNA concentrations are time- dependent variables, interactions are represented by functional and differential between variables, and retarded arguments are usually the time of transcription and translation, or time delay feed- back loop. For the concentration xi (t) of the factor i at time t, the evolution in time is described by the system: dxi (t) dt = axi (t) + bxi (t− τ) = fi (xj (t) , xj (t− τ) , j ∈ Ni) (i = 1, 2, ..., n) , (1) where fi is a nonlinear real-valued function of the states (xj (t) , xj (t− τ)) of the vertices j ∈ Ni, that interact with the factor i. An interesting and important problem appears due to the nonlinearity of fi, 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 de- velop mathematical models that take into consid- eration the post-transcriptional regulation. Also, the models should take into account that the in- teractions 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 diffi- cult to satisfy with a standard modeling formalism, and answers may lie in hybrid approaches [11]. As an intuitive example of nonlinearity by char- acterizing 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 Biomath 3 (2014), 1407231, http://dx.doi.org/10.11145/j.biomath.2014.07.231 Page 2 of 15 http://dx.doi.org/10.11145/j.biomath.2014.07.231 S. Nikolov, Modelling and Analysis of miRNA Regulation miR-X Pol IIEnhancer promoter miRNA Type A miRNA circuit Type B miRNA circuit Gene Z Enhancer Promoter Z mRNA No Z protein Steady state levels of Z protein miRNA degradation/inhibition Pol II Gene Y Pol IIEnhancer Promoter Transcription off Y mRNA Promoter No Y protein Fig. 1. Integrating miRNAs into regulatory networks. Type A circuits might be involved in maintaining protein steady-state levels and homeostasis. Type B circuits can serve as a surveillance mechanism to suppress leaky transcription of target genes. 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 non- linearity appear combined in systems like miRNA networks, a suitable approach for investigating its dynamics is the systems bifurcation analysis. Bi- furcation theory studies persistence and exchange of qualitative properties of dynamical systems un- der continuous perturbations. From the point of view of the dynamic systems theory, the Hopf bifurcation theorem [22] together with other el- ements of the bifurcation theory are basic analyt- ical 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 “triv- ial” 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 sub- set 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 bista- bility. They also exhibit a rich structure of bifur- cations to more complicated behavior, for exam- ple 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 Biomath 3 (2014), 1407231, http://dx.doi.org/10.11145/j.biomath.2014.07.231 Page 3 of 15 http://dx.doi.org/10.11145/j.biomath.2014.07.231 S. Nikolov, Modelling and Analysis of miRNA Regulation 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, y2, controls its own synthesis through the repression of mRNA, y1. Here we construct one equation, describing the rate of change in the concentration of miRNA, y3, 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 k6 as its degradation rate. According to the first mechanism (in which miRNAs are regarded as an enhancer of mRNA degradation) an extra degrada- tion term, k4y1y3, to the rate equation for mRNA and miRNA (where k4 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 dy1 dt = k1 k2 + k3y n1 2 (t− τ1) −γ1y1 −k4y1y3, dy2 dt = k5y1 (t− τ2) −γ2y2, dy3 dt = l−k6y3 (t− τ3) −k4y1y3, (2) mRNA protein Delayed negative feedback miRNA mRNA degradation RISC τ1 τ2 gene τ3 Fig. 2. Scheme of transcription ( τ1) and transla- tion ( τ2) time delay negative feedback control of protein synthesis 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, n1 is often referred as a Hill coefficient or a cooperativity coefficient, ki (i = 1, ..., 6) and γ1, γ2 are the kinetic rate constants. In the present work, ac- cording to [20, 34], we consider one representative value for Hill coefficient, n1 = 2. III. QUALITATIVE ANALYSIS In this section, we consider the system (2) when the Hill coefficient, n1, is equal to two and all constants of the model are real positive numbers. Furthermore, we investigate the bifur- cation 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 = (y1, y2, y3), represented by Eq. (2) can be an- alytically estimated and are defined by the follow- ing set of algebraic equations, including the rate constants of the model: Biomath 3 (2014), 1407231, http://dx.doi.org/10.11145/j.biomath.2014.07.231 Page 4 of 15 http://dx.doi.org/10.11145/j.biomath.2014.07.231 S. Nikolov, Modelling and Analysis of miRNA Regulation y42 + k5 (k4l + k6γ1) k4γ1γ2 y32 + k2 k3 y22 + k5 k3k4γ1γ2 [k2 (k4l + k6γ1) −k1k4] y2 − k1k 2 5k6 k3k4γ1γ 2 2 = 0,y1 = γ2 k5 y2, y3 = k5l k5k6 + γ2k4y2 . (3) 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 physio- logically feasible fixed point. Let us consider a small perturbation around the fixed point E of the system (2) defined as: yk = yk + xk (k = 1, 2, 3) (4) In the case when n1 = 2, the function k1 k2+k3y 2 2 (t−τ1) can be written as a MacLaurin series: k1 k2 + k3y 2 2 (t− τ1) = k1 δ + k3η = k1 δ ( k3 δ η + 1 ) = k1 δ ( 1 − k3 δ η + ( k3 δ )2 η2 − ( k3 δ )3 η3 + ... ) (5) where δ = k2 + k3y22 and η = 2y2x2 (t− τ1) + x22 (t− τ1). If we take only linear term from (5) and after substitution of (4) into differential equation (2) we have dx1 dt = −c1x1 − c2`−χτ1x2 − c3x3 − c4`−2χτ1x22 −k4x1x3, dx2 dt = k5` −χτ2x1 −γ2x2, dx3 dt = −c5x1 − ( k6` −χτ3 + c3 ) x3 −k4x1x3, (6) where c1 = γ1 + k4y3, c2 = 2k1k3y2 δ2 , c3 = k4y1, c4 = k1k3 δ2 , c5 = k4y3. (7) The characteristic equation of (6) has the form det ( −c1 −χ −c2`−χτ1 −c3 k5` −χτ2 −γ2 −χ 0 −c5 0 −(c3 + k6`−χτ3 ) −χ ) = 0 (8) that is χ3 + K1χ 2 + K2χ + K3 = ` −χτ3 (T1χ 2 + T2χ +T3) + ` −χ(τ1+τ2) (T4χ + T5) + T6` −χ(τ1+τ2+τ3) (9) where K1 = γ2 + c1 + c3, K2 = γ2c1 + c3 (γ2 + c1 − c5) ; K3 = γ2c3 (c1 − c5) , T1 = −k6, T2 = −k6(γ2 + c1), T3 = −γ2k6c1, T4 = −k5c2, T5 = −k5c2c3, T6 = −k5k6c2. (10) Remark 1. We note that χ = 0 is a root of (9) if and only if K3 = T3 = T5 = T6 = 0. In the absence of delays ( τ1 = τ2 = τ3 = 0), E is locally asymptotically stable if p = γ2 + k6 + c1 + c3 > 0, q = γ2 (k6 + c1 + c3) + k5c2 + c3 (c1 − c5) > 0, r = γ2c3 (c1 − c5) + γ2c1k6 + k5c2 (k6 + c3) > 0, R = pq −r > 0. (11) Biomath 3 (2014), 1407231, http://dx.doi.org/10.11145/j.biomath.2014.07.231 Page 5 of 15 http://dx.doi.org/10.11145/j.biomath.2014.07.231 S. Nikolov, Modelling and Analysis of miRNA Regulation Here we note that (9) is a third-degree ex- ponential 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. The case τ1 = τ3 = 0, τ2 > 0. 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 χ3 + K11χ 2 + K21χ + K31 = ` −χτ2 (T4χ + T51) , (12) where K11 = K1 −T1, K21 = K2 −T2, K31 = K3 −T3, T51 = T5 + T6. (13) 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 ∣∣∣∣∣∣∣∣∣∣ m3 − 3mn2 + K11 ( m2 −n2 ) + K21m + K31 = `−mτ2 [(T4m + T51) cosnτ2 + T4nsinnτ2] , −n3 + 3m2n + 2K11mn + K21n = `−mτ2 [T4ncosnτ2 − (T4m + T51) sinnτ2] . (14) 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 ∣∣∣∣ −K11n2 + K31 = T51cosnτ2 + T4nsinnτ2,−n3 + K21n = T4ncosnτ2 −T51sinnτ2, (15) or another one cos nτ2 = T51(K31−K11n2)+T4n2(K21−n2) T251+T 2 4 n 2 , sin nτ2 = n[T4(K31−K11n2)−T51(K21−n2)] T251+T 2 4 n 2 . (16) It is clear that if the first bifurcation point is( n0b, τ 0 b ) , then the other bifurcation points (nb, τb) satisfy nbτb = n0bτ 0 b + 2νπ, ν = 1, 2, ..., ∞. 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: Biomath 3 (2014), 1407231, http://dx.doi.org/10.11145/j.biomath.2014.07.231 Page 6 of 15 http://dx.doi.org/10.11145/j.biomath.2014.07.231 S. Nikolov, Modelling and Analysis of miRNA Regulation n6 + ( K211 − 2K21 ) n4 + (K221 − 2K11K31 −T 24 )n 2 + K231 −T 2 51 = 0. (17) 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 n2 and the left-hand side is positive for very large values of n2 and negative for n = 0 if and only if T 251 > 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 nb is the last positive simple root of (17). Then, in (τb) = inb 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 non- zero speed, i.e. the following transversality condi- tion d(Reχ)dt |τ=τb 6= 0 is satisfied. From (15) we know that τbk corresponding to nb is τbk = 1 nb arccos [ (−T4n4b + (T4K21 −T51K11)n 2 b +T51K31)/(T 2 51 + T 2 4 n 2 b) ] + 2kπ nb , k = 0, 1, 2, ... (18) 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 6= 0. Hence, if denote H (χ, τ2) = χ 3 + K11χ 2 + K21χ + K31 − `−χτ2 (T4χ + T51) , (19) then dχ dτ2 = − ∂H ∂τ2 /∂H ∂χ = (−χ`−χτ2 (T4χ + T51)) /(3χ2 + 2K11χ + K21 + τ2` −χτ2 (T4χ + T51) − `−χτ2T4). (20) Evaluating the real part of this equation at τ2 = τb and setting χ = inb yield dm dτ2 |τ2=τb = d (Reχ) dt |τ2=τb = (n2b[3n 4 b + 2 ( K211 − 2K21 ) n2b + K 2 21 − 2K11K31 −T 24 ]) / (L2 + I2) (21) where L = −3n2b + K21 + τ2 ( −K11n2b + K31 ) −T4cosnbτ2, I = 2K11nb + τ2 ( −n3b + K21nb ) + T4sinnbτ2. (22) Let θ = n2b, then, (17) reduces to g (θ) = θ3 + ( K211 − 2K21 ) θ2 + ( K221 − 2K11K31 −T 2 4 ) θ + K231 −T 2 51. (23) Then, for g′ (θ) we have g′(θ)|τ2=τb = dg dθ |τ2=τb = 3θ2 + 2 ( K211 − 2K21 ) θ + K221 − 2K11K31 −T 2 4 . (24) If nb is the least positive simple root of (17), then Biomath 3 (2014), 1407231, http://dx.doi.org/10.11145/j.biomath.2014.07.231 Page 7 of 15 http://dx.doi.org/10.11145/j.biomath.2014.07.231 S. Nikolov, Modelling and Analysis of miRNA Regulation dg dτ2 |θ=n2b > 0. (25) Hence, dm dτ2 |τ2=τb = d (Reχ) dτ2 |τ2=τb = n2bg ′ ( n2b ) L2 + I2 > 0. (26) According to the Hopf bifurcation theorem [17], we define the following Theorem 2: Theorem 2. If nb 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. The case τ1 = 0; τ2, τ3 > 0. We return to the study of (9) which with τ2, τ3 > 0 has the form χ3 + K1χ 2 + K2χ + K3 = `−χτ3 ( T1χ 2 + T2χ + T3 ) + `−χτ2 (T4χ + T5) + T6` −χτ23, (27) where τ = [τ2, τ3, τ23 = τ2 + τ3] T denotes a point in the time delay space, i.e. τ ∈ Ω ⊂ R3+. Ω is the time delay space and R3+ 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. Previously, we obtain that in the absence of delays, E is locally asymptotically stable if the conditions (11) are valid. By Remark 1, this im- plies that χ = 0 is not root of (27). Further, we introduce the following simple result (which is proved in [25]) using Rouche’s theorem Lemma 1. Consider the exponential polynomial P (χ, `−χτ1, ..., `−χτm ) = χn + p (0) 1 χ n−1 + ... + p (0) n−1χ + p (0) n + [p (1) 1 χ n−1 + ... + p (1) n−1χ + p(1)n ]` −χτ1 + ... + [p (m) 1 χ n−1 + ... + p (m) n−1χ + p(m)n ]` −χτm, where τi ≥ 0 (i = 1, 2, ..., m) and p (i) j (i = 0, 1, ..., m; j = 1, 2, ..., n) are constants. As (τ1, τ2, ..., τm) vary, the sum of the order of the zeros of P (χ, `−χτ1, ..., `−χτm )on the open right half plane can change only if a zero appears on or crosses the imaginary axis. Obviously, in (n > 0) is a root of (27) if and only if n satisfies −n3i−K1n2 + K2ni + K3 = (cosnτ3 − isinnτ3) ( −T1n2 + T2ni + T3 ) + (cosnτ2 − isinnτ2) (T4ni + T5) + T6 (cosnτ23 − isinnτ23) . (28) Separating the real and imaginary parts into (28), we obtain ∣∣∣∣∣∣∣∣∣∣∣∣∣∣ −K1n2 + K3 = T4nsinnτ2 + T5cosnτ2 + ( −T1n2 + T3 ) cosnτ3 + T2nsinnτ3 +T6cosnτ23, −n3 + K2n = T4ncosnτ2 −T5sinnτ2 − ( −T1n2 + T3 ) sinnτ3 + T2ncosnτ3 −T6sinnτ23 . (29) We square and add the equations (29), and after simplifying, we get that τ and n must be Biomath 3 (2014), 1407231, http://dx.doi.org/10.11145/j.biomath.2014.07.231 Page 8 of 15 http://dx.doi.org/10.11145/j.biomath.2014.07.231 S. Nikolov, Modelling and Analysis of miRNA Regulation among the real solutions of n6 + ( K21 −T 2 1 − 2K2 ) n4 + (K22 −T 2 2 −T 2 4 − 2K1K3 + 2T1T3)n2 + T 26 + K 2 3 −T 2 3 −T 2 5 = 2{[T2T4n2 + T5(T1n2 + T3)]cosn(τ2 − τ3) + n[T4(−T1n2 + T3) −T2T5]sinn(τ2 − τ3) −T6[(K1n2 −K3)cosnτ23 + (−n3 + K2n)sinnτ23]} (30) We note that the right-hand side of (30) is always less than 2[|T2T4|n2 + |T5|| − T1n2 + T3| + (|T4|| − T4n2 + T3| − |T2T5|)n − |T6| ( |K1n2 −K3| + |−n3 + K2n| ) ]. Hence if the inequality ω6 + ( K21 −T 2 1 − 2K2 ) ω4 + ( K22 −T 2 2 −T 2 4 − 2K1K3 + 2T1T3 ) ω2 + T 26 + K 2 3 −T 2 3 −T 2 5 > 2[|T2T4|ω 2 + |T5||−T1ω2 + T3| + |T4||−T4ω2 + T3|ω −|T6| ( |K1ω2 −K3|− |−ω3 + K2ω| ) ] (31) 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 K1n 2 = ψ (n) = [K + (T1cosnτ3) n 2 − (T4sinnτ2 + T2sinnτ3) n−T5cosnτ2 −T3cosnτ3 −T6cosnτ23] ≤ K3 + |T1|n2 − (|T2| + |T4|)n−|T3|− |T5|− |T6|, i.e. (K1 −|T1|) n2 + (|T2| + |T4|) − (K3 −|T3|− |T5|− |T6|) = 0. (32) Thus, for n+ we have n+ = 1 2a ( −b + √ b2 + 4ac ) , (33) where a = K1 −|T1| 6= 0, b = |T2| + |T4|, c = K3 −|T3|− |T5|− |T6|. It is clear that n ≤ n+. Rearranging terms, we write (31) as ( ω|−ω2 + K2|− |T6| )2 + ( |K1ω2 −K3| + |T6| )2 + ( −T4ω2 + T3 )2 + (|T2|ω + |T5|) + T 24 ω 2 >( |−T1ω2 + T3| + |T5| )2 + (|−T4ω2 + T3| + T4ω) 2 + (T2 + T4) 2 ω2 + T 22 ω 2 + T 25 . (34) Hence, the following theorem can be formulated Theorem 3. Let K3 −|T3|− |T5|− |T6| 6= 0 and (34) hold. Then there is no change in stability of E. Remark 2. In the special case that τ2 = τ3, the characteristic equation (27) becomes χ3 + K1χ 2 + K2χ + K3 = = `−χτ2 ( T1χ 2 + T24χ + T35 ) + T6` −2χτ2, (35) where T35 = T3 + T5 and T24 = T2 + T4. Therefore, this case is a private one of Theorem 3. Corollary 3.1. If conditions of Theorem 3 are not valid and τbif2́ defined as in Theorem 2, then according to Lemma 1 for any τ2 ∈ [0, τb), there exists a τbif23 (τ2) > 0 ( τ bif 1 (τ2) > 0 respectively) such that the steady state E of system (2) is unstable when τ3 ∈ [0, τbif3 (τ2)) ( τ1 ∈ [0, τbif1 (τ2)) respectively), and an Andronov- Hopf bifurcation takes place. The general case τ1, τ2, τ3 > 0. Similar to previous Section, we set that χ = in (n > 0) is a root of (9) if and only if n satisfies Biomath 3 (2014), 1407231, http://dx.doi.org/10.11145/j.biomath.2014.07.231 Page 9 of 15 http://dx.doi.org/10.11145/j.biomath.2014.07.231 S. Nikolov, Modelling and Analysis of miRNA Regulation n3i−K1n2 + K2in + K3 = (cosnτ3 − isinnτ3) ( −T1n2 + T2in + T3 ) + (cosnτ4 − isinnτ4) (T4in + T5) + (cosnτ5 − isinnτ5) T6, (36) where in this case τ = [τ3, τ4 = τ1 + τ2, τ5 = τ1 + τ2 + τ3] T denotes a point in the time delay space, i.e. τ ∈ Ω ⊂ R3+. Here, Ω is the time delay space and R3+ denotes the set of nonnegative real numbers. Separating the real and imaginary parts into (36), we have ∣∣∣∣ −K1n2 + K3 −A1 = B1 + T6cosnτ5,−n3 + K2n + A2 = B2 −T6sinnτ5 (37) where A1 = ( −T1n2 + T3 ) cosnτ3 + T2nsinnτ3, A2 = ( −T1n2 + T3 ) sinnτ3 −T2ncosnτ3, B1 = T4nsinnτ4 + T5cosnτ4, B2 = T4ncosnτ4 −T5sinnτ4. (38) Adding up the squares of both equations into (37), we have n6 + ( K21 + T 2 1 − 2K2 ) n4 + T 23 −T 2 5 −T 2 6 + ( K22 + T 2 2 −T 2 4 − 2K1K3 − 2T1T3 ) n2 = 2{T6 (−T4nsinnτ3 + T5cosnτ3) − ( −n3 + K2n ) [ ( −T1n2 + T3 ) sinnτ3 −T2ncosnτ3] − ( K1n 2 −K3 ) [ ( −T1n2 + T3 ) cosnτ3 + T2nsinnτ3]}, (39) where τ5−τ4 = τ3. Clearly, the right-hand side of (38) is always less than 2[T6|−T4n + T5|− (−n3 −K1n2 + K2n + K3)|−T1n2 + T2n + T3|]. (40) Hence if the inequality λ6 + ( K21 + T 2 1 − 2K2 ) λ4 + ( K22 + T 2 2 −T 2 4 − 2K1K3 − 2T1T3 ) λ2 + T 23 −T 2 5 −T 2 6 > 2[T6|−T4λ + T5| − ( −λ3 −K1λ2 + K2λ + K3 ) |−T1λ2 + T2λ + T3|], (41) 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 K1n 2 = Ψ (n) = K3 − [ ( −T1n2 + T3 ) cosnτ3 + T2nsinnτ3 + T4nsinτ4 + T5cosnτ4 +T6cosnτ5] ≤ K3 + |T1|n2 − (|T2| + |T4|) n −|T3|− |T5|− |T6|, (42) i.e. we obtain the same formula as (32) (and (33) respectively). Rearranging terms, we write (41) as ( |−T1λ2 + T2λ + T3|−λ3 −K1λ2 + K2λ + K3 )2 + ( |T2λ + T3| + T1λ2 )2 + ( λ2 + K1λ )2 + (T2 −T3λ)2 + ( T 25 + 2K 2 2 ) λ2 + K23 > (|−T4λ + T5| + T6)2 + (T5λ + T4)2 + (K2λ + K3) 2 + (T2λ + T3) 2 + ( λ3 −K3 )2 + ( K1λ 2 −K2λ )2 + T 21 λ 4 + T 23 λ 2. (43) Therefore, the following theorem can be formu- lated Biomath 3 (2014), 1407231, http://dx.doi.org/10.11145/j.biomath.2014.07.231 Page 10 of 15 http://dx.doi.org/10.11145/j.biomath.2014.07.231 S. Nikolov, Modelling and Analysis of miRNA Regulation Theorem 4. Let K3 −|T3|− |T5|− |T6| 6= 0 and (43) hold. Then there is no change in stability of E. Remark 3. In the special case that τ1 = τ2 = τ3, the characteristic equation (9) becomes χ3 + K1χ 2 + K2χ + K3 = ` −χτ3 (T1χ 2 + T2χ + T3) + ` −2χτ3 (T4χ + T5) + T6` −3χτ3. (44) Hence, this case is a private one of Theorem 4. Corollary 4.1. If conditions of Theorem 4 are not valid and τbif3 is defined as in Theorem 2, then according to Lemma 1 for any τ3 ∈ [0, τb), there exists a τbif1 ((τ3)) > 0 ( τ bif 2 (τ3) > 0; τbif4 (τ3) > 0; τ bif 5 (τ3) > 0 respectively) such that the steady state E of system (2) is unstable when τ1 ∈ [0, τbif1 (τ3)) ( τ2 ∈ [0, τ bif 2 (τ3)); τ4 ∈ [0, τbif4 (τ3)); τ5 ∈ [0, τ bif 5 (τ3)) respectively), and an Andronov-Hopf bifurcation takes place. In next Section, we illustrate numerically the existence of the behavior predicted for some values of the rate constants ki (i = 1, ..., 6) , γ1, γ2, l and the time delays τ1, τ2 and τ3. IV. NUMERICAL ANALYSIS In the previous section, we proposed the analyti- cal 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 per- form a numerical analysis of model (2), based on the results previously obtained. Some of the parameter values used in the nu- merical analysis were selected according to [4, 20, 24, 29, 33] in the form k1 = 0.3 [ min−1 ] , k2 = k3 = 1 [ min−1 ] , k4 = 0.3 [ min−1 ] , k5 = 0.5 [ min−1 ] , γ1 = 0.1 [ min−1 ] , γ2 = 0.2 [ min−1 ] , τ1 ∈ [1, 8] , τ2 ∈ [12, 35] . (45) According to [6, 32] proteins (or RNAs) de- graded with a probability that depends on their structure because some of the degradation mecha- nisms 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 de- graded 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 k6 ∈ [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 nu- merical 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 ( y1), the concentration of protein ( y2) and the concentration of miRNA ( y3) 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. Biomath 3 (2014), 1407231, http://dx.doi.org/10.11145/j.biomath.2014.07.231 Page 11 of 15 http://dx.doi.org/10.11145/j.biomath.2014.07.231 S. Nikolov, Modelling and Analysis of miRNA Regulation (a) (b) Fig. 3. Stable regime of system (2) for a) τ1 = τ2 = τ3 = 0 and b) for τ1 = τ3 = 0, τ2 = 12. Model parameters: k1 = k4 = 0.3; k2 = k3 = 1; γ1 = 0.1; γ2 = 0.2; k5 = 0.5; k6 = 0.05; n1 = 2; l = 0.1. Figure 4 depicts the cases when τ2 and τ3 are varied. 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. (a) (b) Biomath 3 (2014), 1407231, http://dx.doi.org/10.11145/j.biomath.2014.07.231 Page 12 of 15 http://dx.doi.org/10.11145/j.biomath.2014.07.231 S. Nikolov, Modelling and Analysis of miRNA Regulation (c) (d) Fig. 4. Dynamic behaviour of system (2). Unsta- ble 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 quasi- periodic motion take place. 21 ττ + 10 3τ25 49 78 99 135 147 unstable 300 Fig. 5. Stable and unstable zones of system (2) at τ1 = 1, τ2 = 13 and τ3 ∈ [10, 300]. V. CONCLUSIONS Gene expression in the human organism is post- transcriptionally regulated by miRNAs. MiRNAs are embedded in complex regulatory networks that involved gene activation, post-translational regu- lation 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 par- ticipate in the regulation of gene expression by ac- celerating the degradation of mRNA or by repress- ing 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 Biomath 3 (2014), 1407231, http://dx.doi.org/10.11145/j.biomath.2014.07.231 Page 13 of 15 http://dx.doi.org/10.11145/j.biomath.2014.07.231 S. Nikolov, Modelling and Analysis of miRNA Regulation in stability will occur. If system (2) starts with a stable equilibrium, which for some delay(s) be- comes 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 bifur- cation 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 anal- ysis can enable to understand the mechanism un- derlying an observed biological process, and at the same time, provide a testable hypothesis for future studies. In this paper we investigate the main pro- cesses of the formation of a protein- transcription time (starts with splicing and polyadenylation of the initial transcript), translation time (the time- span from the emergence of mRNA) and average time delay for degradation of miRNA. Thus, our dynamical predictions can be tested in future ex- periments. REFERENCES [1] A. Andronov, A. Witt and S. Chaikin Theory of Oscilla- tions, Addison-Wesley, 1966. [2] O. Arino, M. Hbib and E. AitDads, Delay Differential Equations and Applications, Springer, 2006. [3] M. Bueno, I. Peres de Castro, M. Malumbres, Control of cell proliferation pathways by microRNAs, Cell Cycle 7(20) (2008) 3143–3148. [4] H. de Jong, Modeling and simulation of gene regulatory systems: a literature review, J. of Comp. Biol. 9(1) (2002) 67–103. [5] H. Freedman, and V. Rao, The trade-off between mutual interference and time lags in predator-prey systems. Bull. Math. Biol. 45(6) (1983) 991-1004. [6] A. Hershko and A. Ciechanover, The ubiquitin system. Annu. Rev. Biochem. 67 (1998) 425–479. [7] R. Hilborn, Chaos and Nonlinear Dynamics: an Introduc- tion for Sciencists and Engineers, Second edtion, Oxford University Press, 2011. [8] C. Kanellopoulou and S. Monticelli, A role for microR- NAs in the development of the immune system and in the pathogenesis of cancer, Seminars in Cancer Biol., 18 (2008) 79–88. [9] Q. Khan and D. Greenhalgh, Hopf bifurcation in epi- demic models with a time delay in vaccination. IMA J. of Math. Applied in Medicine and Biology 16 (1999) 113–142. [10] Q. Khan, Hopf bifurcation in multiparty political sys- tems with time delay in switching. Appl. Math. Lett. 13 (2000) 45–52. [11] F.M. Khan, U. Schmitz, S. Nikolov and et al., Hybrid modeling of the crosstalk between signaling and tran- scriptional networks using ordinary differential equations and multi-valued logic, BBA (Biochimica et Biophysica Acta)- Proteins and Proteomics, 1844(1) (2014) 289–298. [12] R. Khanin and V. Vinciotti, Computational modeling of post-transcriptional gene regulation by microRNAs, J. Comput Biol. 15(3) (2008) 305–316. [13] R. Khanin and D. Higham, 2010. Mathematical and computational modelling of post-transcriptional gene reg- ulation by miRNAs. In: MiRNA profiling in Cancer. A Bioinformatics Perspective, Pan Stanford Publishing Pte. Ltd. Chapter 10 (2010) 197–216. [14] G. Korn and T. Korn, Mathematical Handbook for Scientists and Engineers. McGraw-Hill, 1968. [15] H. Lodish, A. Berk, P. Matsudaira and et al., Molecular Cell Biology, Freeman and Company, 2004. [16] J. Lytle, Th. Yario and J. Steitz, Target mRNAs are repressed as efficiently by microRNA-binding sites in the 5’UTR as in the 3’UTR, PNAS 104(23) (2007) 9667– 9672. [17] J. Marsden and M. McCracen, The Hopf Bifurcation and its Applications, Springer-Verlag, 1976 [18] Matlab, The MathWorks, Inc. 2010. http://www.mathworks.com/, [19] R. May, Stability and Complexity in Model Ecosystems. Princeton University Press, 1973. [20] N. Monk, Oscillatory expression of Hes1, p53, and NF- kappaB driven by transcriptional delay, Curr. Biol. 13 (2003) 1409–1413. [21] D. Nicholl, An Introduction to Genetic Engineering. Cambridge Univ. Press, 2008. http://dx.doi.org/10.1017/CBO9780511800986 [22] S. Nikolov, J. Vera, U. Schmitz and O. Wolkenhauer, A model-based strategy to investigate the role of microRNA regulation in cancer signalling networks. Theory in Bio- sciences 130(1) (2011) 55–69. Biomath 3 (2014), 1407231, http://dx.doi.org/10.11145/j.biomath.2014.07.231 Page 14 of 15 http://www.mathworks.com/ http://dx.doi.org/10.1017/CBO9780511800986 http://dx.doi.org/10.11145/j.biomath.2014.07.231 S. Nikolov, Modelling and Analysis of miRNA Regulation [23] S. Nikolov, Stability and Andronov-Hopf bifurcation of a system with three time delays, Journal of Mathematics, 2013 (2013) art. ID 347071(11 pages). [24] K. Rateitschak and O. Wolkenhauer, Intracellular delay limits cyclic changes in gene expression, Math. Biosci. 205 (2007) 163–179. [25] Sh. Ruan and J. Wei, On the zeros of transcendental functions with applications to stability of delay differen- tial equations with two delays, Dynamics of Continuous, Discrete and Impulsive Systems, Series A: Mathematical Analysis 10 (2003) 863–874. [26] A. Sensse and M. Eiswirth, Feedback loops for chaos in activator-inhibitor systems, The Journal of Chemical Physics 122 (2005) 044516. [27] L. Shilnikov, A. Shilnikov, D. Turaev and L. Chua, Methods of Qualitative Theory in Nonlinear Dynamics. Part II, World Scientific, 2001. [28] P. Smolen, D. Baxter and J. Byrne, Effects of macro- molecular transport and stochastic fluctuations on dynam- ics of genetic regulatory systems, Am. J. Physiol. 277(4 Pt 1) (1999) C777–C790. [29] I. Swameye, T. Mueller, J. Timmer, and et al., Identifi- cation of nucleocytoplasmic cycling as a remote sensor in cellular signaling by data-based modelling, PNAS 100 (2003) 1028–1033. [30] I.A. Swinburne, P.A. Silver, Intron delays and transcrip- tional timing during development. Dev Cell. 14(3) (2008) 324–330. [31] J. Tsang, J. Zhu and A. van Ondenaarden, MicroRNA- mediated feedback and feedforward loops are recurrent networks motifs in mammals, Mol. Cell., 26(5) (2007) 753–767. http://dx.doi.org/10.1016/j.molcel.2007.05.018 [32] X. Wang, Y. Li, X. Xu and Y. Wang, Toward a system level understanding of microRNA pathway via mathemat- ical modelling, BioSystems 100 (2010) 31-38. [33] Zh. Xie, H. Yang, W. Liu, and M. Hwang, The role of microRNA in the delayed negative feedback regulation of gene expression, Bioch. and Bioph. Research Comm. 358 (2007) 722–726. [34] S. Zeiser, H. Liebscher, H. Tiedemann, and et al., Number of active transcription factors binding sites is essential for the Hes7 oscillator. Theor. Biol. Med. Model 23 (2006) 11. [35] S. Zeiser, O. Rivera, C. Kuttler and et al., Oscillations of Hes7 caused by negative autoregulation and ubiquiti nation, Comput. Biol. Chem. 32 (2008) 48–52. Biomath 3 (2014), 1407231, http://dx.doi.org/10.11145/j.biomath.2014.07.231 Page 15 of 15 http://dx.doi.org/10.1016/j.molcel.2007.05.018 http://dx.doi.org/10.11145/j.biomath.2014.07.231 Introduction Model Qualitative Analysis Numerical Analysis Conclusions References