# Stability and Andronov-Hopf Bifurcation of a System with Three Time Delays

Authors: Svetoslav Nikolov
Journal: Journal of Mathematics
 (2013)

---

## Abstract

A general system of three autonomous ordinary differential equations with three discrete time delays is considered. With respect to the delays, we investigate the local stability of equilibria by analyzing the corresponding characteristic equation. Using the Hopf bifurcation theorem, we predict the occurrence of a limit cycle bifurcation for the time delay parameters. Thus, some new mathematical results are obtained. Finally, the above mentioned criteria are applied to a system modelling miRNA regulation.


---

## Body

## 1. Introduction

The future status of many systems arising from engineering, physics, mechanics, biochemistry, or systems in biology is determined not only by their current behavior but is also by their history. Such phenomena are called delay or genetic effects. Time lags in continuous systems can produce complex dynamics and instabilities. In the recent years, many mathematical models that have appeared in the literature involve either a single discrete delay to investigate the role of phosphorylation, negative or positive feedback regulation of transcription factors [1–4], and gene expression multistability [5]. Several papers [6–13] among others consider systems with two or three discrete delays. When there is more than one delay in the equations, the local theory for stability is not fully complete.A dynamical system which is not finite dimensional is called an infinite dimensional dynamical system. A class of infinite dimensional system can be determined by functional differential equations of the retarded type. To introduce such equations, we letCr denote the set C[[-r,0],Rn] with the norm defined by ∥ϕ∥=max{|ϕ(t):-r≤t≤0|} where ∥·∥ denotes a norm (e.g., the Euclidean norm). Given a function x(·) defined on -r≤t≤0, let xt be the function determined by xt(s)=x(t+s) for -r≤s≤0. A retarded functional differential equation (with delay r) is an equation of the form
(1)x.(t)=F(t,xt),
where F∈C[Ω,Rn] and Ω is an open set in (R×Cr). The notations of the dynamical system are determined by varying (t0,xt0=ψ) over some appropriate subset of Ω. As far as we know, the first statement similar to the Hopf theorem (concerning the bifurcation of periodic solutions from a singular point of an ordinary differential equation) [14] for retarded functional differential equations was given in [15].In this paper, we consider stability and bifurcation behavior of a three-dimensional system of autonomous ordinary differential delay equations of the form(2)x′(t)=f(x(t),x(t-ζ),y(t),y(t-ξ),z(t),z(t-η)),y′(t)=g(x(t),x(t-ζ),y(t),y(t-ξ),z(t),z(t-η)),z′(t)=h(x(t),x(t-ζ),y(t),y(t-ξ),z(t),z(t-η)),
where(′=d/dt);f,g,h:R×R×R×R×R×R→R, are such that solutions to initial value problems exist and are continuable.


## 2. Derivation of the Characteristic Equation

We assume thatf,  g, and h have continuous first partial derivatives with respect to their arguments and that there exists unique x->0, y->0, and z->0, such that f(x-,x-,y-,y-,z-,z-)=g(x-,x-,y-,y-,z-,z-)=h(x-,x-,y-,y-,z-,z-)=0. We utilize the following notation: fi,i=1,…,6 that represents the partial derivative of f with respect to its ith argument evaluated E-(x-,x-,y-,y-,z-,z-). A similar meaning is given to gj,j=1,…,6, and hk,k=1,…,6. Hence, for example, f2=(∂f(x-,x-,y-,y-,z-,z-)/∂x(t-ζ)) and h2=(∂h(x-,x-,y-,y-,z-,z-)/∂x(t-ζ)). The variational system with respect to E- is
(3)(u(t)v(t)ρ(t))′=(f1f3f5g1g3g5h1h3h5)(u(t)v(t)ρ(t))+(f2f4f6g2g4g6h2h4h6)(u(t-ζ)v(t-ξ)ρ(t-η)).
We seek a solution of (3) of the form u(t)=c1ℓχt, v(t)=c2ℓχt and ρ(t)=c3ℓχt. For nontrivial solutions, this leads to the following characteristic equation in χ:
(4)χ3+a1χ2+a2χ+a3+ℓ-χζ(a4χ2+a7χ+a13)+ℓ-χξ(a5χ2+a8χ+a14)+ℓ-χη(a6χ2+a9χ+a15)+ℓ-χ(ζ+ξ)(a10χ+a16)+ℓ-χ(ζ+η)(a11χ+a17)+ℓ-χ(ξ+η)(a12χ+a18)+a19ℓ-χ(ζ+ξ+η)=0,
where
(5)a1=-(f1+g3+h5),a2=f1(g3+h5)-f3g1-f5h1+g3h5-g5h3,a3=f1(g5h3-g3h5)+f3(g1h5-g5h1)+f5(g3h1-g1h3),a4=-f2,a5=-g4,a6=-h6,a7=f2(g3+h5)-f3g2-f5h2,a8=g4(f1+h5)-f4g1-g5h4,a9=h6(f1+g3)-f6h1-g6h3,a10=f2g4-f4g2,a11=f2h6-f6h2,a12=g4h6-g6h4,a13=f2(g5h3-g3h5)+f3(g2h5-g5h2)+f5(g3h2-g2h3),a14=f1(g5h4-g4h5)+f4(g1h5-g5h1)+f5(g4h1-g1h4),a15=f1(g6h3-g3h6)+f3(g1h6-g6h1)+f6(g3h1-g1h3),a16=f2(g5h4-g4h5)+f4(g2h5-g5h2)+f5(g4h2-g2h4),a17=f2(g6h3-g3h6)+f3(g2h6-g6h2)+f6(g3h2-g2h3),a18=f1(g6h4-g4h6)+f4(g1h6-g6h1)+f6(g4h1-g1h4),a19=f2(g6h4-g4h6)+f4(g2h6-g6h2)+f6(g4h2-g2h4).Remark 1.
We note thatχ=0 is a root of (4) if and only if a3=a13=a14 = a15=a16 = a17=a18 = a19=0.It is well known that the stability of the equilibrium stateE- depends on the sign of the real parts of the roots of (4). If χ=m+in(m,n∈R) satisfies (4), then m and n are real solutions; that is, we rewrite (4) in terms of its real and imaginary parts as
(6)m3-3mn2+a1(m2-n2)+a2m+a3+ℓ-mζ[a4(m2-n2)cosnζ+2a4mnsinnζ+a7(mcosnζ+nsinnζ)+a13cosnζ{m2}]+ℓ-mξ[a5(m2-n2)cosnξ+2a5mnsinnξ+a8(mcosnξ+nsinnξ)+a14cosnξ{m2}]+ℓ-mη[a6(m2-n2)cosnη+2a6mnsinnη+a9(mcosnη+nsinnη)+a15cosnη{m2}]+ℓ-m(ζ+ξ)[a10(mcosn(ζ+ξ)+nsinn(ζ+ξ))+ℓ-m(ζ+ξ)+a16cosn(ζ+ξ)]+ℓ-m(ζ+η)[a11(mcosn(ζ+η)+nsinn(ζ+η))+ℓ-m(ζ+η)+a17cosn(ζ+η)]+ℓ-m(ξ+η)[a12(mcosn(ξ+η)+nsinn(ξ+η))+ℓ-m(ξ+η)+a18cosn(ξ+η)]+a19ℓ-m(ζ+ξ+η)cosn(ζ+ξ+η)=0,3m2n-n3+2a1mn+a2n+ℓ-mζ[2a4mncosnζ-a4(m2-n2)sinnζ+a7(ncosnζ-msinnζ)-a13sinnζ{m2}]+ℓ-mξ[2a5mncosnξ-a5(m2-n2)sinnξ+a8(ncosnξ-msinnξ)-a14sinnξ{m2}]+ℓ-mη[2a6mncosnη-a6(m2-n2)sinnη+a9(ncosnη-msinnη)-a15sinnη{m2}]+ℓ-m(ζ+ξ)[a10(ncosn(ζ+ξ)-msinn(ζ+ξ))+ℓ-m(ζ+ξ)-a16sinn(ζ+ξ)]+ℓ-m(ζ+η)[a11(ncosn(ζ+η)-msinn(ζ+η))+ℓ-m(ζ+η)-a17sinn(ζ+η)]+ℓ-m(ξ+η)[a12(ncos(ξ+η)-msinn(ξ+η))+ℓ-m(ξ+η)-a18sinn(ξ+η)]-a19ℓ-m(ζ+ξ+η)sinn(ζ+ξ+η)=0.
In the absence of delays (ζ=ξ=η=0), E- is locally asymptotically stable if
(7)p=a1+a4+a5+a6>0,q=a2+a7+a8+a9+a10+a11+a12>0,r=a3+a13+a14+a15+a16+a17+a18+a19>0,R=pq-r>0.
Because of the presence of three different discrete delays in (2), 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 when secondly one delay is equal to zero. Previously, this approach is used for system with two delays [9, 12, 13, 16–18].


## 3. The Caseζ=ξ=0 and η>0

Settingζ=ξ=0 in (4), the system (6) becomes
(8)m3-3mn2+K1(m2-n2)+K2m+K3+ℓ-mη[{(T2+a6(m2-n2))}T1(mcosnη+nsinnη)+cosnη(T2+a6(m2-n2))+2a6mnsinnη]=0,3m2n-n3+2K4mn+K2n+ℓ-mη[{(T2+a6(m2-n2))}T1(ncosnη-msinnη)-sinnη(T2+a6(m2-n2))]=0,
where
(9)K1=a1+a4+a5,K2=a2+a7+a8+a10,K3=a3+a13+a14+a16,K4=a1+a4+a5+a6,T1=a9+a11+a12,T2=a15+a17+a18+a19.To find the first bifurcation point, we look for purely imaginary rootsχ=±in, n∈R of (4) (when ζ=ξ=0); that is, we set m=0. Then, the above two equations are reduced to
(10)-K1n2+K3=-nT1sinnη-(T2-a6n2)cosnη,-n3+K2n=-nT1cosnη+(T2-a6n2)sinnη
or another one
(11)cosnη=(n2-K2)n2T1+(T2-a6n2)(K1n2-K3)n2T12+(T2-a6n2)2,sinnη=(K1n2-K3)nT1-n(T2-a6n2)(n2-K2)n2T12+(T2-a6n2)2.
We note here that it is not possible for (nT1)2 and (T2-a6n2)2 to be both zero and n=0 can be a solution of (11) if K3=T2. If the first bifurcation point is (nb0,ηb0), then the other bifurcation points (nb,ηb) satisfy nbηb=nb0ηb0+2νπ, ν=1,2,…,∞.One can notice that ifn is a solution of (10) (or (11)), then so is -n. Hence, in the following we only investigate for positive solutions n of (10), or (11) respectively. By squaring the two equations into system (10) and then adding them, it follows that
(12)n6+(K12-2K2-a62)n4+(K22-2K1K3-T12+2a6T2)n2+K32-T22=0.
As E- is locally asymptotically stable at η=0, it satisfies the Routh-Hurwitz conditions for stability for a cubic polynomial [19, 20]. Equation (12) is a cubic in n2 and the left-hand side is positive for very large values of n2 and also at n=0. Suppose that conditions of Lemma 1(I) in [19] are satisfied; that is, (12) has at least one positive real simple root. Moreover, to apply the Hopf bifurcation theorem, according to [19], the following theorem in this situation applies.Theorem 2.
Suppose thatnb is the least positive simple root of (12). Then, in(ηb)=inb is a simple root of (4) (at ζ=ξ=0) and m(η)+in(η) is differentiable with respect to η in a neighbourhood of η=ηb.To establish an Andronov-Hopf bifurcation atη=ηb, we need to show that the following transversality condition dm/dη|η=ηb≠0 is satisfied.Hence, if we denote(13)H(χ,η)=χ3+K1χ2+K2χ+K3+ℓ-χη(a6χ2+T1χ+T2)
then
(14)dχdη=-∂H/∂η∂H/∂χ=χℓ-χη(a6χ2+T1χ+T2)×(3χ2+2K1χ+K2-ηℓ(a6χ2+T1χ+T2)+ℓ-χη(2a6χ+T1){3χ2+})-1.Evaluating the real part of this equation atη=ηb and setting χ=inb yield
(15)dmdη|η=ηb=Re(dχdη)|η=ηb=nb2[3nb2+2(K12-2K2-a62)nb2+K22-2K1K3+2a6T2-T12]×(L2+I2)-1,
where L=-3nb2+K2+ηb(-K1nb2+K3)+T1cosnbηb+2a6nbsinnbηb and I=2K1nb+ηb(-nb3+K2nb)+2a6nbcosnbηb-T1sinnbηb.Letθ=ηb2, then (12) reduces to
(16)g(θ)=θ3+(K12-2K2-a62)θ2+(K22-2K1K3+2a6T2-T12)θ+K32-T22.
Then for g′(θ), we have
(17)g′(θ)|η=ηb=dgdθ|η=ηb=3θ2+2(K12-2K2-a62)θ+K22-2K1K3+2aT2-T12.Ifηb is the least positive simple root of (12), then
(18)dgdη|θ=nb2>0.Hence,(19)dmdη|η=ηb=Re(dχdη)|η=ηb=nb2g′(nb2)L2+I2>0.According to the Hopf bifurcation theorem [21], we define the following theorem.Theorem 3.
Ifnb is the least positive root of (12), then an Andronov-Hopf bifurcation occurs as η passes through ηb.Corollary 4.
Whenη<ηb, then the steady state E- of system (2) is locally asymptotically stable.


## 4. The Caseζ=0;ξ,η>0

We return to the study of (4) which with ξ,η>0 has the form
(20)χ3+K41χ2+K5χ+K6+ℓ-χτ1(a5χ2+T3χ+T4)+ℓ-χτ2(a6χ2+T5χ+T6)+ℓ-χτ3(a12χ+T7)=0,
where K41=a1+a4, K5=a2+a7, K6=a3+a13, T3=a8+a10, T4=a14+a16, T5=a9+a11, T6=a15+a17, T7=a18+a19, and τ=[τ1=ξ,τ2=η,τ3=ξ+η]T denotes a point in the time delay space; that is, τ∈Ω⊂R+3. Ω is the time delay space and R+3 denotes the set of nonnegative real numbers. In order to assess the stability of (2) with respect to any delay τ, one should know where all χ roots of (20) lie on the complex plane. Equation (20) has infinitely many roots on the complex plane due to the transcendental term ℓ-χτ. This makes the analytical stability assessment intractable.In Section2, we obtain that in the absence of delays, E- is locally asymptotically stable if the conditions (7) are valid. By Remark 1, this implies that χ=0 is not root of (20). Further, we introduce the following simple result (which was proved by [18]) using Rouche’s theorem.Lemma 5.
Consider the exponential polynomialP(χ,ℓ-χτ1,…,ℓ-χτm)=χn+p1(0)χn-1+⋯+pn-1(0)χ+pn(0) + [p1(1)χn-1+⋯+pn-1(1)χ+pn(1)]ℓ-χτ1+⋯+[p1(m)χn-1+⋯+pn-1(m)χ+pn(m)]ℓ-χτm, where τi≥0(i=1,2,…,m) and pj(i)(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 there or crosses the imaginary axis.Obviously,in(n>0) is a root of (20) if and only if n satisfies
(21)-n3i-K41n2+K5ni+K6+(cosnτ1-isinnτ1)(-a5n2+T3ni+T4)+(cosnτ2-isinnτ2)(-a6n2+T5ni+T6)+(cosnτ3-isinnτ3)(a12ni+T7)=0.
Separating the real and imaginary parts into (21), we obtain
(22)-K41n2+K6+a12nsinnτ3+T7cosnτ3=(a5n2-T4)cosnτ1-T3nsinnτ1+(a6n2-T6)cosnτ2-T5nsinnτ2,-n3+K5n+a12ncosnτ3-T7sinnτ3=(T4-a5n2)sinnτ1-T3ncosnτ1+(T6-a6n2)sinnτ2-T5ncosnτ2.
We square and add (22), and after simplifying, we get that τ and n must be among the real solutions of
(23)n6+(K412-a52-a62-2K5)n4+[K52+a122-T32-T52+2(a5T4+a6T6-K4K6)]n2+K62+T72=2[(K41n2-K6)(T7cosnτ3+a12nsinnτ3)+n(n2-K5)(a12ncosnτ3-T7sinnτ3)+(a5n2-T4)(a6n2-T6)cosn(τ1-τ2)+T5n(a5n2-T4)sinn(τ1-τ2)+T3n(T6-a6n2)sinn(τ1-τ2)+T3T5n2cosn(τ1-τ2)].
We note that the right-hand side of (23) is always less than 2[|K41n2-K6||T7+a12n|+|n2-K5||a12n-T7|n+|a5n2-T4||a6n2-T6|+|T5||a5n2-T4|n+|T3||T6 − a6n2|n + |T3T5|n2]. Hence if the inequality
(24)ω6+(K412-a52-a62-2K5)ω4+[K52+a122-T32-T52+2(a5T4+a6T6-K41K6)]ω2+K62+T72>2[|T3T5|ω2+|K41ω2-K6||T7+a12ω|+|ω2-K5||a12ω-T7|ω+|a5ω2-T4||a6ω2-T6|+|T5||a5ω2-T4|ω+|T3||T6-a6ω2|ω]
has no real solution on 0<ω<n+, then (23) cannot be satisfied. Note that n+ is the positive solution of first equation in (6) (when m=0 and ζ=0), which we write as
(25)(a1+a4)n2=ψ(n)=[a3+a13-(a5cosnξ+a6cosnη)n2+(a8sinnξ+a9sinnη+a10sinnξ+a11sinnη+a12sinn(ξ+η))n+a14cosnξ+a15cosnη+a16cosnξ+a17cosnη+a18cosn(ξ+η)+a19cosn(ξ+η){n2}]≤a3+a13+|a14|+|a15|+|a16|+|a17|+|a18|+|a19|+(|a8|+|a9|+|a10|+|a11|+|a12|)n-(|a5|+|a6|)n2,i.e.,(a1+a4+|a5|+|a6|)n2-(|a8|+|a9|+|a10|+|a11|+|a12|)n-(a3+a13+⋯+|a19|)=0.
Thus, for n+, we have
(26)n+=12a(b+b2+4ac),
where a=a1+a4+|a5|+|a6|≠0, b=|a8|+⋯+|a12| and c=a3+a13+|a14|+⋯+|a19|. It is clear that n≤n+.Rearranging terms, we write (24) as
(27)(|K41ω2-K6|-|T7+a12ω|)2+(ω|ω2-K5|-|a12ω-T7|)2+(|a5|ω2-|T4|)2+(|a6|ω2-|T6|)2+(T32+T52)ω2>(|a5ω2-T4|+|a6ω2-T6|)2+(|a5ω2-T4|+|T5|ω)2+(|T3|ω+|T6-a6ω2|)2+(|T3|+|T5|)2ω2+a122ω2+T42+T62+T72.
Hence, the following theorem can be formulated.Theorem 6.
Leta3+a13+|a14|+⋯+|a19|≠0 and (27) hold. Then there is no change in stability of E-.Remark 7.
In the special case thatξ=η, the characteristic equation (20) becomes
(28)χ3+K41χ2+K5χ+K6+ℓ-χτ1(T56χ2+T35χ+T46)+ℓ-2χτ1(a12χ+T7)=0,
where T56=a5+a6, T35=T3+T5 and T46=T4+T6. Therefore, this case is a private one of Theorem 6.Corollary 8.
If conditions of Theorem6 are not valid and τ2bif is defined as in Theorem 3, then according to Lemma 5 for any τ2∈[0,τb), there exists a τ3bif(τ2)>0(τ1bif(τ2)>0 resp.) such that the steady state E- of system (2) is unstable when τ3∈[0,τ3bif(τ2))(τ1∈[0,τ1bif(τ2)) resp.), and an Andronov-Hopf bifurcation takes place.


## 5. The General Caseζ,ξ,η>0

Similar to Section4, we set that χ=in(n>0) is a root of (4) if and only if n satisfies
(29)-n3i-a1n2+a2ni+a3+(cosnτ1-isinnτ1)(-a4n2+a7ni+a13)+(cosnτ2-isinnτ2)(-a5n2+a8ni+a14)+(cosnτ3-isinnτ3)(-a6n2+a9ni+a15)+(cosnτ4-isinnτ4)(a6ni+a16)+(cosnτ5-isinnτ5)(a11ni+a17)+(cosnτ6-isinnτ6)(a12ni+a18)+a19(cosnτ7-isinnτ7)=0,
where τ=[τ1=ζ, τ2=ξ, τ3=η, τ4=ζ+ξ, τ5=ζ+η, τ6=ξ+η, τ7=ζ+ξ+η]T denotes a point in the time delay space; that is, τ∈Ω⊂R+7. Ω is the time delay space and R+7 denotes the set of nonnegative real numbers.Separating the real and imaginary parts into (29), we have
(30)-a1n2+a3-A1-B1-C1=-D1-E1-F1-a19cosnτ7,-n3+a2n+A2+B2+C2=D2+E2+F2+a19sinnτ7,
where
(31)A1=(a4n2-a13)cosnτ1-a7nsinnτ1,B1=(a5n2-a14)cosnτ2-a8nsinnτ2,C1=(a6n2-a15)cosnτ3-a9nsinnτ3,D1=a10nsinnτ4+a16cosnτ4,E1=a11nsinnτ5+a17cosnτ5,F1=a12nsinnτ6+a18cosnτ6,A2=(a4n2-a13)sinnτ1+a7ncosnτ1,B2=(a5n2-a14)sinnτ2+a8ncosnτ2,C2=(a6n2-a15)sinnτ3+a9ncosnτ3,D2=a10ncosnτ4-a16sinnτ4,E2=a11ncosnτ5-a17sinnτ5,F2=a12ncosnτ6-a18sinnτ6.
Adding up the squares of both equations into (30), we have
(32)n6+S1n4+S2n2+S3=2{(a4n2-a13)×[(a3-a1n2)cosnτ1-(-n3+a2n)sinnτ1]+a7n[(a3-a1n2)sinnτ1+(-n3+a2n)cosnτ1]+a19(a12nsinnτ1+a18cosnτ1)+(a5n2-a14)×[(a3-a1n2)cosnτ2-(-n3+a2n)sinnτ2]+a19(a11nsinnτ2+a17cosnτ2)+a8n[(a3-a1n2)sinnτ2+(-n3+a2n)cosnτ2]+a19(a10nsinnτ3+a16cosnτ3)+(a6n2-a15)×[(a3-a1n2)cosnτ3-(-n3+a2n)sinnτ3]+a9n[(a3-a1n2)sinnτ3+(-n3+a2n)cosnτ3]+(a5n2-a14)[{n2}a7nsinn(τ1-τ2)+(a5n2-a14)-(a4n2-a13)cosn(τ1-τ2)]+a17[a18cosn(τ1-τ2)-a12nsinn(τ1-τ2)]+a11n[a12ncosn(τ1-τ2)+a18sinn(τ1-τ2)]+(a6n2-a15)[{n2}a7nsinn(τ1-τ3)+(a6n2-a15)-(a4n2-a13)cosn(τ1-τ3)]+a10n[a12ncosn(τ1-τ3)+a18sinn(τ1-τ3)]+a16[a18cosn(τ1-τ3)-a12nsinn(τ1-τ3)]-a9n[{n2}a7ncosn(τ1-τ3)+(a4n2-a13)sinn(τ1-τ3)]+(a6n2-a15)[{n2}a8nsinn(τ2-τ3)+(a6n2-a15)-(a5n2-a14)cosn(τ2-τ3)]+a10n[a11ncosn(τ2-τ3)+a17sinn(τ2-τ3)]-a16[a11nsinn(τ2-τ3)-a17cosn(τ2-τ3)]-a9n[{n2}a8ncosn(τ2-τ3)+(a5n2-a14)sinn(τ2-τ3)]}+4a19{sinn(τ1+τ2)[a10ncosn(τ1+τ2+τ3)sinn(τ1+τ2)-a16sinn(τ1+τ2+τ3)]+4a19+sinn(τ1+τ3)[a11ncosn(τ1+τ2+τ3)+4a19+sinn(τ1+τ3)-a17sinn(τ1+τ2+τ3)]+4a19+sinn(τ2+τ3)[a12ncosn(τ1+τ2+τ3)+sinn(τ2+τ3)+4a19-a18sinn(τ1+τ2+τ3)]},
where
(33)τ4+τ7=2(τ1+τ2)+τ3,τ5+τ6=τ1+τ2+2τ3,τ5+τ7=τ2+2(τ1+τ3),τ6+τ7=τ1+2(τ2+τ3),τ4-τ5=τ2-τ3,τ4-τ6=τ1-τ3,τ5-τ6=τ1-τ2,S1=a12+a42+a52+a62-2a2,S2=a22+a72+a82+a92-a102-a112-a122,S3=a32+a132+a142+a152-a162-a172-a182-a192.
Clearly, the right-hand side of (32) is always less than
(34)2{|n3-a1n2-a2n+a3|×[|a4n2-a13|+|a5n2-a14|+|a6n2-a15|]+a19[|a12n+a18|+|a11n+a17|+|a10n+a16|]+n|-n3-a1n2+a2n+a3|(a7+a8+a9)+n(a10+a11)|a12n+a18|-(a16+a17)|a12n-a18|+|-a4n2+a7n+a13|[|a5n2-a14|+|a6n2-a15|]-n(a8+a9)|a4n2+a7n-a13|-a9n|a5n2+a8n-a14|+|a6n2-a15||-a5n2+a8n+a14||n3-a1n2-a2n+a3|+a10n|a11n+a17|-a16|a11n-a17|}+4a19[|a10n-a16|+|a11n-a17|+|a12n-a18|].
Hence, if the inequality
(35)λ6+S1λ4+S2λ2+S3>2{|λ3-a1λ2-a2λ+a3|×[|a4λ2-a13|+|a5λ2-a14|+|a6λ2-a15|]+a19[|a12λ+a18|+|a11λ+a17|+|a10λ+a16|]+(a7+a8+a9)|-λ3-a1λ2+a2λ+a3|λ+(a10+a11)|a12λ+a18|λ-(a16+a17)|a12λ-a18|+|-a4λ2+a7λ+a13|[|a5λ2-a14|+|a6λ2-a15|]-(a8+a9)|a4λ2+a7λ-a13|λ-a9λ|a5λ2+a8λ-a14|+|a6λ2-a15||-a5λ2+a8λ+a14||λ3-a1λ2-a2λ+a3|+a10λ|a11λ+a17|-a16|a11λ-a17|}+4a19[|a10λ-a16|+|a11λ-a17|+|a12λ-a18|]
has no real solution on 0<λ<n+, then (32) cannot be satisfied. Similar to Section 4, we note that n+ is the positive solution of first equation in (6) (when m=0), which is written as
(36)a1n2=Ψ(n)=[a3-(a4cosnζ+a5cosnξ+a6cosnη)n2+(a7sinnζ+a8sinnξ+a9sinnη+a10sinn(ζ+ξ)+a11sinn(ζ+η)+a12sinn(ξ+η))n+a13cosnζ+a14cosnξ+a15cosnη+a16cosn(ζ+ξ)+a17cosn(ζ+η)+a18cosn(ξ+η)(a4cosnζ+a5cosnξ+a6cosnη)+a19cosn(ζ+ξ+η)]≤a3+|a13|+|a14|+|a15|+|a16|+|a17|+|a18|+|a19|+(|a7|+|a8|+|a9|+|a10|+|a11|+|a12|)n-(|a4|+|a5|+|a6|)n2,i.e.,(a1+|a4|+|a5|+|a6|)n2-(|a7|+|a8|+|a9|+|a10|+|a11|+|a12|)n-(a3+|a13|+⋯+|a19|)=0.
Thus, for n+, we have
(37)n+=12A(B+B2+4AC),
where A=a1+|a4|+|a5|+|a6|≠0, B=|a7|+⋯+|a12| and C=a3+|a13|+⋯+|a19|. It is clear that n≤n+.Rearranging terms, we write (35) as
(38)(|λ3-a1λ2-a2λ+a3|-|a4λ2-a13|)2+(|λ3-a1λ2-a2λ+a3|-|a5λ2-a14|)2+(a11λ-|a12λ+a18|)2+(|λ3-a1λ2-a2λ+a3|-|a6λ2-a15|)2+(|-λ3-a1λ2+a2λ+a3|-λ|a7+a8+a9|)2+(a17+|a12λ-a18|)2+(|a6λ2-a15|-|-a5λ2+a8λ+a14|)2+(a9λ+|a5λ2+a8λ-a14|)2+λ(1+|a1λ4+a2a3|)2+(a16-|a12λ-a18|)2+(|-a4λ2+a7λ+a13|-|a5λ2-a14|)2+(|-a4λ2+a7λ+a13|-|a6λ2-a15|)2+(a10λ-|a12λ+a18|)2+(|a4λ2+a7λ-a13|+λ|a8+a9|)2+λ(|a3+a1a2|-λ2)2+λ(|a1a3+a4a13|+λ)2+(a10λ-|a11λ+a17|)2+(a16+|a11λ-a17|)2+(2a19-|a10λ-a16|)2+(2a19-|a11λ-a17|)2+(a5+a14λ2)2+(a6+a15λ2)2+(a8-a9λ2)2+a72λ2>2(λ3-a1λ2-a2λ+a3)2+(a5λ2-a14)2+2(a6λ2-a15)2+(-λ3-a1λ2+a2λ+a3)2+λ(a7+a8+a9)2+(-a5λ2+a8λ+a14)2+(a5λ2+a8λ-a14)2+λ[(a1λ4+a2a3)2+(a1a2+a3)2+λ4+(a1a3+a4a13)2+λ2{(a1λ4+a2a3)2}]+2(-a4λ2+a7λ+a13)2+(a4λ2+a7λ-a13)2+2(a10λ-a16)2+2(a11λ-a17)2+(a19+|a12λ+a18|)2+(a19+|a11λ+a17|)2+(a19+|a10λ+a16|)2+(2a19+|a12λ-a18|)2+(a92+a142+a152)λ4+(a92+3a122)λ2+a52+a62+a82+a102+a162+2a172+3a182+2a192.
Therefore, the following theorem can be formulated.Theorem 9.
Leta3+|a13|+|a14|+⋯+|a19|≠0 and (38) hold. Then there is no change in stability of E-.Remark 10.
In the special case thatζ=ξ=η, the characteristic equation (4) becomes
(39)χ3+a1χ2+a2χ+a3+ℓ-χη[(a4+a5+a6)χ2+(a7+a8+a9)χ+a13+a14+a15{χ2}]+ℓ-2χη[(a10+a11+a12)χ+a16+a17+a18]+a19ℓ-3χη=0.Hence, this case is a private one of Theorem9.Corollary 11.
If conditions of Theorem9 are not valid and τ3bif is defined as in Theorem 3, then according to Lemma 5 for any τ3∈[0,τb), there exists a τ1bif(τ3)>0(τ2bif(τ3)>0; τ4bif(τ3)>0; τ5bif(τ3)>0; τ6bif(τ3)>0; τ7bif(τ3)>0 resp.) such that the steady state  E- of system (2) is unstable when τ1∈[0,τ1bif(τ3))(τ2∈[0,τ2bif(τ3)); τ4∈[0,τ4bif(τ3)); τ5∈[0,τ5bif(τ3)); τ6∈[0,τ6bif(τ3))τ7∈[0,τ7bif(τ3)) resp.), and an Andronov-Hopf bifurcation take place.


## 6. Application: Numerical Analysis

In this section, we apply the material and results of the previous sections to a nonlinear mathematical model with three discrete time delays, which investigate the dynamics of the gene expression regulated by the miRNA. For a general discussion of the model, we refer the reader to [22, 23]. Hence, we consider the following system:
(40)dy1dt=k1k2+k3y2n1(t-τ1)-γ1y1-k4y1y3,dy2dt=k5y1(t-τ2)-γ2y2,dy3dt=l-k6y3(t-τ3)-k4y1y3,
where the protein, y2, controls its own synthesis through the repression of mRNA, y1,  y3 is the concentration of miRNA, τ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 and γ2 are the kinetic rate constants. According to [2, 24], we consider one representative value for Hill coefficient, n1=2. Model (40) is a modified version of a model considered in [22]; that is, in the present paper our model contains three time delays.The fixed points of the system,E-=(y-1,y-2,y-3), represented by (40) can be analytically estimated and are defined by the following set of algebraic equations, including the rate constants of the model:
(41)y-24+k5(k4l+k6γ1)k4γ1γ2y-23+k2k3y-22+k5k3k4γ1γ2[k2(k4l+k6γ1)-k1k4]y-2-k1k52k6k3k4γ1γ22=0,y-1=γ2k5y-2,y-3=k5lk5k6+γ2k4y-2.According to Descarte’s rule [19, 25], the first equation in (41) has only one real positive root, which ensures that the system has only one physiologically feasible fixed point.Let us consider a small perturbation around the fixed pointE- of the system (40) defined as
(42)yk=y-k+xk(k=1,2,3).
In the case when n1=2, the function  k1/(k2+k3y22(t-τ1)) can be written as a MacLaurin series
(43)k1k2+k3y22(t-τ1)=k1δ+k3φ=k1δ((k3/δφ)+1)=k1δ(1-k3δφ+(k3δ)2φ2-(k3δ)3φ3+⋯),
where δ=k2+k3y-22 and φ=2y-2x2(t-τ1)+x22(t-τ1). If we take only linear term from (43) and after substitution of (42) into differential equation (40), we have,
(44)dx1dt=-c1x1-c2ℓ-χτ1x2-c3x3-c4ℓ-2χτ1x22-k4x1x3,dx2dt=k5ℓ-χτ2x1-γ2x2,dx3dt=-c5x1-(k6ℓ-χτ3+c3)x3-k4x1x3,
where
(45)c1=γ1+k4y-3,c2=2k1k3y-2δ2,c3=k4y-1,c4=k1k3δ2,c5=k4y-3.
It is seen that notations in Section 2 here are
(46)f2=f3=f6=g1=g4=g5=g6=h2=h3=h4=0,f1=-c1,f4=-c2,f5=-c3,g2=k5,g3=-γ2,h1=-c5,h5=-c6,h6=-k6.
In the next few figures we illustrate numerically the validity of the criteria proved Section 5, that is, when the three time delays are different from zero. For our simulations according to [22], the numerical values of the model parameters are
(47)k1=0.3[min-1],k2=k3=1[min-1],k4=0.3[min-1],k5=0.5[min-1],k6∈[0.05,0.3],γ1=0.1[min-1],γ2=0.2[min-1],τ1∈[1,8],τ2∈[12,35],τ3∈[10,300],l=0.1.
In order to compare the predictions with numerical results, the governing equations of the model (40) were solved numerically using MATLAB [26]. In Figure 1(a), 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). 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 criteria written in (7). In Figure 1(b), the unstable solutions of system (40) are shown. It is seen that for larger values of τ3, than bifurcation one, after Andronov-Hopf bifurcation the stable limit cycle (self oscillations) occurs and sustained oscillations take place. In other words, in this case the conditions of Theorem 9 are not satisfied and the steady state of system (40) is unstable.(a) Stable regime of system (40) for τ1=τ2=τ3=0 and (b) unstable regime (sustained oscillations) for τ1=1, τ2=13, τ3=60.
(a)
(b)


## 7. Conclusions

In this paper we have considered a system of three equations with three discrete time delays. 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 behavior is devoted to the use of the bifurcation analysis. Particularly, a Hopf bifurcation theorem was employed. The basic view that the time delays ζ, ξ, and η are a key factor in the dynamical behavior of system (2) has been confirmed by analytical and numerical calculations.



---

*Source: 347071-2013-03-24.xml*