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

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.

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 let 𝐶_𝑟 denote the set 𝐶[[−𝑟,0],𝑅^𝑛] with the norm defined by ‖𝜙‖=max{|𝜙(𝑡)∶−𝑟≤𝑡≤0|} where ‖⋅‖ denotes a norm (e.g., the Euclidean norm). Given a function 𝑥(⋅) defined on −𝑟≤𝑡≤0, let 𝑥_𝑡 be the function determined by 𝑥_𝑡(𝑠)=𝑥(𝑡+𝑠) for −𝑟≤𝑠≤0. A retarded functional differential equation (with delay 𝑟) is an equation of the form 𝑥̇(𝑡)=𝐹𝑡,𝑥_𝑡, (1) where 𝐹∈𝐶[Ω,𝑅^𝑛] and Ω is an open set in (𝑅×𝐶_𝑟). The notations of the dynamical system are determined by varying (𝑡₀,𝑥_𝑡0=𝜓) 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) where(^′=𝑑/𝑑𝑡);𝑓,𝑔,ℎ∶𝑅×𝑅×𝑅×𝑅×𝑅×𝑅→𝑅, are such that solutions to initial value problems exist and are continuable.

2. Derivation of the Characteristic Equation

We assume that 𝑓,  𝑔, and ℎ have continuous first partial derivatives with respect to their arguments and that there exists unique 𝑥>0, 𝑦>0, and 𝑧>0, such that 𝑓(𝑥,𝑥,𝑦,𝑦,𝑧,𝑧)=𝑔(𝑥,𝑥,𝑦,𝑦,𝑧,𝑧)=ℎ(𝑥,𝑥,𝑦,𝑦,𝑧,𝑧)=0. We utilize the following notation: 𝑓_𝑖,𝑖=1,…,6 that represents the partial derivative of 𝑓 with respect to its 𝑖th argument evaluated 𝐸(𝑥,𝑥,𝑦,𝑦,𝑧,𝑧). A similar meaning is given to 𝑔_𝑗,𝑗=1,…,6, and ℎ_𝑘,𝑘=1,…,6. Hence, for example, 𝑓₂=(𝜕𝑓(𝑥,𝑥,𝑦,𝑦,𝑧,𝑧)/𝜕𝑥(𝑡−𝜁)) and ℎ₂=(𝜕ℎ(𝑥,𝑥,𝑦,𝑦,𝑧,𝑧)/𝜕𝑥(𝑡−𝜁)). The variational system with respect to 𝐸 is 𝑢(𝑡)𝑣(𝑡)𝜌(𝑡)^′=𝑓₁𝑓₃𝑓₅𝑔₁𝑔₃𝑔₅ℎ₁ℎ₃ℎ₅𝑢(𝑡)𝑣(𝑡)𝜌(𝑡)+𝑓₂𝑓₄𝑓₆𝑔₂𝑔₄𝑔₆ℎ₂ℎ₄ℎ₆𝑢()𝑣()𝜌.𝑡−𝜁𝑡−𝜉𝑡−𝜂 (3) We seek a solution of (3) of the form 𝑢(𝑡)=𝑐₁ℓ^𝜒𝑡, 𝑣(𝑡)=𝑐₂ℓ^𝜒𝑡 and 𝜌(𝑡)=𝑐₃ℓ^𝜒𝑡. For nontrivial solutions, this leads to the following characteristic equation in 𝜒: 𝜒³+𝑎₁𝜒²+𝑎₂𝜒+𝑎₃+ℓ^−𝜒𝜁𝑎₄𝜒²+𝑎₇𝜒+𝑎₁₃+ℓ^−𝜒𝜉𝑎₅𝜒²+𝑎₈𝜒+𝑎₁₄+ℓ^−𝜒𝜂𝑎₆𝜒²+𝑎₉𝜒+𝑎₁₅+ℓ^{−𝜒(𝜁+𝜉)}𝑎₁₀𝜒+𝑎₁₆+ℓ^{−𝜒(𝜁+𝜂)}𝑎₁₁𝜒+𝑎₁₇+ℓ^{−𝜒(𝜉+𝜂)}𝑎₁₂𝜒+𝑎₁₈+𝑎₁₉ℓ^{−𝜒(𝜁+𝜉+𝜂)}=0, (4) where 𝑎₁=−𝑓₁+𝑔₃+ℎ₅,𝑎₂=𝑓₁𝑔₃+ℎ₅−𝑓₃𝑔₁−𝑓₅ℎ₁+𝑔₃ℎ₅−𝑔₅ℎ₃,𝑎₃=𝑓₁𝑔₅ℎ₃−𝑔₃ℎ₅+𝑓₃𝑔₁ℎ₅−𝑔₅ℎ₁+𝑓₅𝑔₃ℎ₁−𝑔₁ℎ₃,𝑎₄=−𝑓₂,𝑎₅=−𝑔₄,𝑎₆=−ℎ₆,𝑎₇=𝑓₂𝑔₃+ℎ₅−𝑓₃𝑔₂−𝑓₅ℎ₂,𝑎₈=𝑔₄𝑓₁+ℎ₅−𝑓₄𝑔₁−𝑔₅ℎ₄,𝑎₉=ℎ₆𝑓₁+𝑔₃−𝑓₆ℎ₁−𝑔₆ℎ₃,𝑎₁₀=𝑓₂𝑔₄−𝑓₄𝑔₂,𝑎₁₁=𝑓₂ℎ₆−𝑓₆ℎ₂,𝑎₁₂=𝑔₄ℎ₆−𝑔₆ℎ₄,𝑎₁₃=𝑓₂𝑔₅ℎ₃−𝑔₃ℎ₅+𝑓₃𝑔₂ℎ₅−𝑔₅ℎ₂+𝑓₅𝑔₃ℎ₂−𝑔₂ℎ₃,𝑎₁₄=𝑓₁𝑔₅ℎ₄−𝑔₄ℎ₅+𝑓₄𝑔₁ℎ₅−𝑔₅ℎ₁+𝑓₅𝑔₄ℎ₁−𝑔₁ℎ₄,𝑎₁₅=𝑓₁𝑔₆ℎ₃−𝑔₃ℎ₆+𝑓₃𝑔₁ℎ₆−𝑔₆ℎ₁+𝑓₆𝑔₃ℎ₁−𝑔₁ℎ₃,𝑎₁₆=𝑓₂𝑔₅ℎ₄−𝑔₄ℎ₅+𝑓₄𝑔₂ℎ₅−𝑔₅ℎ₂+𝑓₅𝑔₄ℎ₂−𝑔₂ℎ₄,𝑎₁₇=𝑓₂𝑔₆ℎ₃−𝑔₃ℎ₆+𝑓₃𝑔₂ℎ₆−𝑔₆ℎ₂+𝑓₆𝑔₃ℎ₂−𝑔₂ℎ₃,𝑎₁₈=𝑓₁𝑔₆ℎ₄−𝑔₄ℎ₆+𝑓₄𝑔₁ℎ₆−𝑔₆ℎ₁+𝑓₆𝑔₄ℎ₁−𝑔₁ℎ₄,𝑎₁₉=𝑓₂𝑔₆ℎ₄−𝑔₄ℎ₆+𝑓₄𝑔₂ℎ₆−𝑔₆ℎ₂+𝑓₆𝑔₄ℎ₂−𝑔₂ℎ₄. (5)

Remark 1. We note that 𝜒=0 is a root of (4) if and only if 𝑎₃=𝑎₁₃=𝑎₁₄ = 𝑎₁₅=𝑎₁₆ = 𝑎₁₇=𝑎₁₈ = 𝑎₁₉=0.

It is well known that the stability of the equilibrium state 𝐸 depends on the sign of the real parts of the roots of (4). If 𝜒=𝑚+𝑖𝑛(𝑚,𝑛∈𝑅) satisfies (4), then 𝑚 and 𝑛 are real solutions; that is, we rewrite (4) in terms of its real and imaginary parts as 𝑚³−3𝑚𝑛²+𝑎₁𝑚²−𝑛²+𝑎₂𝑚+𝑎₃+ℓ^−𝑚𝜁𝑎₄𝑚²−𝑛²cos𝑛𝜁+2𝑎₄𝑚𝑛sin𝑛𝜁+𝑎₇(𝑚cos𝑛𝜁+𝑛sin𝑛𝜁)+𝑎₁₃cos𝑛𝜁+ℓ^−𝑚𝜉𝑎₅𝑚²−𝑛²cos𝑛𝜉+2𝑎₅𝑚𝑛sin𝑛𝜉+𝑎₈(𝑚cos𝑛𝜉+𝑛sin𝑛𝜉)+𝑎₁₄cos𝑛𝜉+ℓ^−𝑚𝜂𝑎₆𝑚²−𝑛²cos𝑛𝜂+2𝑎₆𝑚𝑛sin𝑛𝜂+𝑎₉𝑚cos𝑛𝜂+𝑛sin𝑛𝜂+𝑎₁₅cos𝑛𝜂+ℓ^{−𝑚(𝜁+𝜉)}𝑎₁₀(𝑚cos𝑛(𝜁+𝜉)+𝑛sin𝑛(𝜁+𝜉))+𝑎₁₆cos𝑛(𝜁+𝜉)+ℓ^{−𝑚(𝜁+𝜂)}𝑎₁₁𝑚cos𝑛𝜁+𝜂+𝑛sin𝑛𝜁+𝜂+𝑎₁₇cos𝑛𝜁+𝜂+ℓ^{−𝑚(𝜉+𝜂)}𝑎₁₂𝑚cos𝑛𝜉+𝜂+𝑛sin𝑛𝜉+𝜂+𝑎₁₈cos𝑛𝜉+𝜂+𝑎₁₉ℓ^{−𝑚(𝜁+𝜉+𝜂)}cos𝑛𝜁+𝜉+𝜂=0,3𝑚²𝑛−𝑛³+2𝑎₁𝑚𝑛+𝑎₂𝑛+ℓ^−𝑚𝜁2𝑎₄𝑚𝑛cos𝑛𝜁−𝑎₄𝑚²−𝑛²sin𝑛𝜁+𝑎₇(𝑛cos𝑛𝜁−𝑚sin𝑛𝜁)−𝑎₁₃sin𝑛𝜁+ℓ^−𝑚𝜉2𝑎₅𝑚𝑛cos𝑛𝜉−𝑎₅𝑚²−𝑛²sin𝑛𝜉+𝑎₈(𝑛cos𝑛𝜉−𝑚sin𝑛𝜉)−𝑎₁₄sin𝑛𝜉+ℓ^−𝑚𝜂2𝑎₆𝑚𝑛cos𝑛𝜂−𝑎₆𝑚²−𝑛²sin𝑛𝜂+𝑎₉𝑛cos𝑛𝜂−𝑚sin𝑛𝜂−𝑎₁₅sin𝑛𝜂+ℓ^{−𝑚(𝜁+𝜉)}𝑎₁₀(𝑛cos𝑛(𝜁+𝜉)−𝑚sin𝑛(𝜁+𝜉))−𝑎₁₆sin𝑛(𝜁+𝜉)+ℓ^{−𝑚(𝜁+𝜂)}𝑎₁₁𝑛cos𝑛𝜁+𝜂−𝑚sin𝑛𝜁+𝜂−𝑎₁₇sin𝑛𝜁+𝜂+ℓ^{−𝑚(𝜉+𝜂)}𝑎₁₂𝑛cos𝜉+𝜂−𝑚sin𝑛𝜉+𝜂−𝑎₁₈sin𝑛𝜉+𝜂−𝑎₁₉ℓ^{−𝑚(𝜁+𝜉+𝜂)}sin𝑛𝜁+𝜉+𝜂=0. (6) In the absence of delays (𝜁=𝜉=𝜂=0), 𝐸 is locally asymptotically stable if 𝑝=𝑎₁+𝑎₄+𝑎₅+𝑎₆>0,𝑞=𝑎₂+𝑎₇+𝑎₈+𝑎₉+𝑎₁₀+𝑎₁₁+𝑎₁₂>0,𝑟=𝑎₃+𝑎₁₃+𝑎₁₄+𝑎₁₅+𝑎₁₆+𝑎₁₇+𝑎₁₈+𝑎₁₉>0,𝑅=𝑝𝑞−𝑟>0. (7) 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 𝑚³−3𝑚𝑛²+𝐾₁𝑚²−𝑛²+𝐾₂𝑚+𝐾₃+ℓ^−𝑚𝜂𝑇₁𝑇𝑚cos𝑛𝜂+𝑛sin𝑛𝜂+cos𝑛𝜂₂+𝑎₆𝑚²−𝑛²+2𝑎₆𝑚𝑛sin𝑛𝜂=0,3𝑚²𝑛−𝑛³+2𝐾₄𝑚𝑛+𝐾₂𝑛+ℓ^−𝑚𝜂𝑇₁𝑇𝑛cos𝑛𝜂−𝑚sin𝑛𝜂−sin𝑛𝜂₂+𝑎₆𝑚²−𝑛²=0, (8) where 𝐾₁=𝑎₁+𝑎₄+𝑎₅,𝐾₂=𝑎₂+𝑎₇+𝑎₈+𝑎₁₀,𝐾₃=𝑎₃+𝑎₁₃+𝑎₁₄+𝑎₁₆,𝐾₄=𝑎₁+𝑎₄+𝑎₅+𝑎₆,𝑇₁=𝑎₉+𝑎₁₁+𝑎₁₂,𝑇₂=𝑎₁₅+𝑎₁₇+𝑎₁₈+𝑎₁₉. (9)

To find the first bifurcation point, we look for purely imaginary roots 𝜒=±𝑖𝑛, 𝑛∈𝑅 of (4) (when 𝜁=𝜉=0); that is, we set 𝑚=0. Then, the above two equations are reduced to −𝐾₁𝑛²+𝐾₃=−𝑛𝑇₁sin𝑛𝜂−𝑇₂−𝑎₆𝑛²cos𝑛𝜂,−𝑛³+𝐾₂𝑛=−𝑛𝑇₁cos𝑛𝜂+𝑇₂−𝑎₆𝑛²sin𝑛𝜂 (10) or another one 𝑛cos𝑛𝜂=²−𝐾₂𝑛²𝑇₁+𝑇₂−𝑎₆𝑛²𝐾₁𝑛²−𝐾₃𝑛²𝑇²₁+𝑇₂−𝑎₆𝑛²²,𝐾sin𝑛𝜂=¹𝑛²−𝐾₃𝑛𝑇₁−𝑛𝑇₂−𝑎₆𝑛²𝑛²−𝐾₂𝑛²𝑇²₁+𝑇₂−𝑎₆𝑛²². (11) We note here that it is not possible for (𝑛𝑇₁)² and (𝑇₂−𝑎₆𝑛²)² to be both zero and 𝑛=0 can be a solution of (11) if 𝐾₃=𝑇₂. If the first bifurcation point is (𝑛⁰_𝑏,𝜂⁰_𝑏), then the other bifurcation points (𝑛_𝑏,𝜂_𝑏) satisfy 𝑛_𝑏𝜂_𝑏=𝑛⁰_𝑏𝜂⁰_𝑏+2𝜈𝜋, 𝜈=1,2,…,∞.

One can notice that if 𝑛 is a solution of (10) (or (11)), then so is −𝑛. Hence, in the following we only investigate for positive solutions 𝑛 of (10), or (11) respectively. By squaring the two equations into system (10) and then adding them, it follows that 𝑛⁶+𝐾²₁−2𝐾₂−𝑎²₆𝑛⁴+𝐾²₂−2𝐾₁𝐾₃−𝑇²₁+2𝑎₆𝑇₂𝑛²+𝐾²₃−𝑇²₂=0. (12) As 𝐸 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 𝑛² and the left-hand side is positive for very large values of 𝑛² and also at 𝑛=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 that 𝑛_𝑏 is the least positive simple root of (12). Then, 𝑖𝑛(𝜂_𝑏)=𝑖𝑛_𝑏 is a simple root of (4) (at 𝜁=𝜉=0) and 𝑚(𝜂)+𝑖𝑛(𝜂) is differentiable with respect to 𝜂 in a neighbourhood of 𝜂=𝜂_𝑏.

To establish an Andronov-Hopf bifurcation at 𝜂=𝜂_𝑏, we need to show that the following transversality condition 𝑑𝑚/𝑑𝜂|_{𝜂=𝜂𝑏}≠0 is satisfied.

Hence, if we denote 𝐻𝜒,𝜂=𝜒³+𝐾₁𝜒²+𝐾₂𝜒+𝐾₃+ℓ^−𝜒𝜂𝑎₆𝜒²+𝑇₁𝜒+𝑇₂ (13) then 𝑑𝜒𝑑𝜂=−𝜕𝐻/𝜕𝜂𝜕𝐻/𝜕𝜒=𝜒ℓ^−𝜒𝜂𝑎₆𝜒²+𝑇₁𝜒+𝑇₂×3𝜒²+2𝐾₁𝜒+𝐾₂−𝜂ℓ𝑎₆𝜒²+𝑇₁𝜒+𝑇₂+ℓ^−𝜒𝜂2𝑎₆𝜒+𝑇₁⁻¹. (14)

Evaluating the real part of this equation at 𝜂=𝜂_𝑏 and setting 𝜒=𝑖𝑛_𝑏 yield 𝑑𝑚𝑑𝜂_{𝜂=𝜂𝑏}=Re𝑑𝜒𝑑𝜂_{𝜂=𝜂𝑏}=𝑛²_𝑏3𝑛²_𝑏+2𝐾²₁−2𝐾₂−𝑎²₆𝑛²_𝑏+𝐾²₂−2𝐾₁𝐾₃+2𝑎₆𝑇₂−𝑇²₁×𝐿²+𝐼²⁻¹, (15) where 𝐿=−3𝑛²_𝑏+𝐾₂+𝜂_𝑏(−𝐾₁𝑛²_𝑏+𝐾₃)+𝑇₁cos𝑛_𝑏𝜂_𝑏+2𝑎₆𝑛_𝑏sin𝑛_𝑏𝜂_𝑏 and 𝐼=2𝐾₁𝑛_𝑏+𝜂_𝑏(−𝑛³_𝑏+𝐾₂𝑛_𝑏)+2𝑎₆𝑛_𝑏cos𝑛_𝑏𝜂_𝑏−𝑇₁sin𝑛_𝑏𝜂_𝑏.

Let 𝜃=𝜂²_𝑏, then (12) reduces to 𝑔(𝜃)=𝜃³+𝐾²₁−2𝐾₂−𝑎²₆𝜃²+𝐾²₂−2𝐾₁𝐾₃+2𝑎₆𝑇₂−𝑇²₁𝜃+𝐾²₃−𝑇²₂. (16) Then for 𝑔^′(𝜃), we have 𝑔^′(𝜃)_{𝜂=𝜂𝑏}=𝑑𝑔𝑑𝜃_{𝜂=𝜂𝑏}=3𝜃²+2𝐾²₁−2𝐾₂−𝑎²₆𝜃+𝐾²₂−2𝐾₁𝐾₃+2𝑎𝑇₂−𝑇²₁. (17)

If 𝜂_𝑏 is the least positive simple root of (12), then 𝑑𝑔𝑑𝜂_{𝜃=𝑛2𝑏}>0. (18)

Hence, 𝑑𝑚𝑑𝜂_{𝜂=𝜂𝑏}=Re𝑑𝜒𝑑𝜂_{𝜂=𝜂𝑏}=𝑛²_𝑏𝑔^′𝑛²_𝑏𝐿²+𝐼²>0. (19)

According to the Hopf bifurcation theorem [21], we define the following theorem.

Theorem 3. If 𝑛_𝑏 is the least positive root of (12), then an Andronov-Hopf bifurcation occurs as 𝜂 passes through 𝜂_𝑏.

Corollary 4. When 𝜂<𝜂_𝑏, then the steady state 𝐸 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 𝜒³+𝐾₄₁𝜒²+𝐾₅𝜒+𝐾₆+ℓ^−𝜒𝜏1𝑎₅𝜒²+𝑇₃𝜒+𝑇₄+ℓ^−𝜒𝜏2𝑎₆𝜒²+𝑇₅𝜒+𝑇₆+ℓ^−𝜒𝜏3𝑎₁₂𝜒+𝑇₇=0, (20) where 𝐾₄₁=𝑎₁+𝑎₄, 𝐾₅=𝑎₂+𝑎₇, 𝐾₆=𝑎₃+𝑎₁₃, 𝑇₃=𝑎₈+𝑎₁₀, 𝑇₄=𝑎₁₄+𝑎₁₆, 𝑇₅=𝑎₉+𝑎₁₁, 𝑇₆=𝑎₁₅+𝑎₁₇, 𝑇₇=𝑎₁₈+𝑎₁₉, and 𝜏=[𝜏₁=𝜉,𝜏₂=𝜂,𝜏₃=𝜉+𝜂]^𝑇 denotes a point in the time delay space; that is, 𝜏∈Ω⊂𝑅³₊. Ω is the time delay space and 𝑅³₊ 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 Section 2, we obtain that in the absence of delays, 𝐸 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 polynomial 𝑃(𝜒,ℓ^−𝜒𝜏1,…,ℓ^{−𝜒𝜏𝑚})=𝜒^𝑛+𝑝₁⁽⁰⁾𝜒^𝑛−1+⋯+𝑝⁽⁰⁾_𝑛−1𝜒+𝑝_𝑛⁽⁰⁾ + [𝑝₁⁽¹⁾𝜒^𝑛−1+⋯+𝑝⁽¹⁾_𝑛−1𝜒+𝑝_𝑛⁽¹⁾]ℓ^−𝜒𝜏1+⋯+[𝑝₁^(𝑚)𝜒^𝑛−1+⋯+𝑝^(𝑚)_𝑛−1𝜒+𝑝_𝑛^(𝑚)]ℓ^{−𝜒𝜏𝑚}, where 𝜏_𝑖≥0(𝑖=1,2,…,𝑚) and 𝑝_𝑗^(𝑖)(𝑖=0,1,…,𝑚;𝑗=1,2,…,𝑛) are constants. As (𝜏₁,𝜏₂,…,𝜏_𝑚) vary, the sum of the order of the zeros of 𝑃(𝜒,ℓ^−𝜒𝜏1,…,ℓ^{−𝜒𝜏𝑚}) on the open right half plane can change only if a zero appears there or crosses the imaginary axis.

Obviously, 𝑖𝑛(𝑛>0) is a root of (20) if and only if 𝑛 satisfies −𝑛³𝑖−𝐾₄₁𝑛²+𝐾₅𝑛𝑖+𝐾₆+cos𝑛𝜏₁−𝑖sin𝑛𝜏₁−𝑎₅𝑛²+𝑇₃𝑛𝑖+𝑇₄+cos𝑛𝜏₂−𝑖sin𝑛𝜏₂−𝑎₆𝑛²+𝑇₅𝑛𝑖+𝑇₆+cos𝑛𝜏₃−𝑖sin𝑛𝜏₃𝑎₁₂𝑛𝑖+𝑇₇=0. (21) Separating the real and imaginary parts into (21), we obtain −𝐾₄₁𝑛²+𝐾₆+𝑎₁₂𝑛sin𝑛𝜏₃+𝑇₇cos𝑛𝜏₃=𝑎₅𝑛²−𝑇₄cos𝑛𝜏₁−𝑇₃𝑛sin𝑛𝜏₁+𝑎₆𝑛²−𝑇₆cos𝑛𝜏₂−𝑇₅𝑛sin𝑛𝜏₂,−𝑛³+𝐾₅𝑛+𝑎₁₂𝑛cos𝑛𝜏₃−𝑇₇sin𝑛𝜏₃=𝑇₄−𝑎₅𝑛²sin𝑛𝜏₁−𝑇₃𝑛cos𝑛𝜏₁+𝑇₆−𝑎₆𝑛²sin𝑛𝜏₂−𝑇₅𝑛cos𝑛𝜏₂. (22) We square and add (22), and after simplifying, we get that 𝜏 and 𝑛 must be among the real solutions of 𝑛⁶+𝐾²₄₁−𝑎²₅−𝑎²₆−2𝐾₅𝑛⁴+𝐾²₅+𝑎²₁₂−𝑇²₃−𝑇²₅+2𝑎₅𝑇₄+𝑎₆𝑇₆−𝐾₄𝐾₆𝑛²+𝐾²₆+𝑇²₇=2𝐾₄₁𝑛²−𝐾₆𝑇₇cos𝑛𝜏₃+𝑎₁₂𝑛sin𝑛𝜏₃+𝑛𝑛²−𝐾₅𝑎₁₂𝑛cos𝑛𝜏₃−𝑇₇sin𝑛𝜏₃+𝑎₅𝑛²−𝑇₄𝑎₆𝑛²−𝑇₆cos𝑛𝜏₁−𝜏₂+𝑇₅𝑛𝑎₅𝑛²−𝑇₄sin𝑛𝜏₁−𝜏₂+𝑇₃𝑛𝑇₆−𝑎₆𝑛²sin𝑛𝜏₁−𝜏₂+𝑇₃𝑇₅𝑛²cos𝑛𝜏₁−𝜏₂. (23) We note that the right-hand side of (23) is always less than 2[|𝐾₄₁𝑛²−𝐾₆||𝑇₇+𝑎₁₂𝑛|+|𝑛²−𝐾₅||𝑎₁₂𝑛−𝑇₇|𝑛+|𝑎₅𝑛²−𝑇₄||𝑎₆𝑛²−𝑇₆|+|𝑇₅||𝑎₅𝑛²−𝑇₄|𝑛+|𝑇₃||𝑇₆ − 𝑎₆𝑛²|𝑛 + |𝑇₃𝑇₅|𝑛²]. Hence if the inequality 𝜔⁶+𝐾²₄₁−𝑎²₅−𝑎²₆−2𝐾₅𝜔⁴+𝐾²₅+𝑎²₁₂−𝑇²₃−𝑇²₅+2𝑎₅𝑇₄+𝑎₆𝑇₆−𝐾₄₁𝐾₆𝜔²+𝐾²₆+𝑇²₇>2𝑇₃𝑇₅𝜔²+𝐾₄₁𝜔²−𝐾₆𝑇₇+𝑎₁₂𝜔+𝜔²−𝐾₅𝑎₁₂𝜔−𝑇₇𝜔+𝑎₅𝜔²−𝑇₄𝑎₆𝜔²−𝑇₆+𝑇₅𝑎₅𝜔²−𝑇₄𝜔+𝑇₃𝑇₆−𝑎₆𝜔²𝜔 (24) has no real solution on 0<𝜔<𝑛₊, then (23) cannot be satisfied. Note that 𝑛₊ is the positive solution of first equation in (6) (when 𝑚=0 and 𝜁=0), which we write as 𝑎₁+𝑎₄𝑛²(𝑛)=𝑎=𝜓₃+𝑎₁₃−𝑎₅cos𝑛𝜉+𝑎₆cos𝑛𝜂𝑛²+𝑎₈sin𝑛𝜉+𝑎₉sin𝑛𝜂+𝑎₁₀sin𝑛𝜉+𝑎₁₁sin𝑛𝜂+𝑎₁₂sin𝑛𝜉+𝜂𝑛+𝑎₁₄cos𝑛𝜉+𝑎₁₅cos𝑛𝜂+𝑎₁₆cos𝑛𝜉+𝑎₁₇cos𝑛𝜂+𝑎₁₈cos𝑛𝜉+𝜂+𝑎₁₉cos𝑛𝜉+𝜂≤𝑎₃+𝑎₁₃+𝑎₁₄+𝑎₁₅+𝑎₁₆+𝑎₁₇+𝑎₁₈+𝑎₁₉+𝑎₈+𝑎₉+𝑎₁₀+𝑎₁₁+𝑎₁₂𝑛−𝑎₅+𝑎₆𝑛²𝑎,i.e.,₁+𝑎₄+𝑎₅+𝑎₆𝑛²−𝑎₈+𝑎₉+𝑎₁₀+𝑎₁₁+𝑎₁₂𝑛−𝑎₃+𝑎₁₃+⋯+𝑎₁₉=0. (25) Thus, for 𝑛₊, we have 𝑛₊=12𝑎𝑏+𝑏²+4𝑎𝑐, (26) where 𝑎=𝑎₁+𝑎₄+|𝑎₅|+|𝑎₆|≠0, 𝑏=|𝑎₈|+⋯+|𝑎₁₂| and 𝑐=𝑎₃+𝑎₁₃+|𝑎₁₄|+⋯+|𝑎₁₉|. It is clear that 𝑛≤𝑛₊.

Rearranging terms, we write (24) as 𝐾₄₁𝜔²−𝐾₆−𝑇₇+𝑎₁₂𝜔²+𝜔𝜔²−𝐾₅−𝑎₁₂𝜔−𝑇₇²+𝑎₅𝜔²−𝑇₄²+𝑎₆𝜔²−𝑇₆²+𝑇²₃+𝑇²₅𝜔²>𝑎₅𝜔²−𝑇₄+𝑎₆𝜔²−𝑇₆²+𝑎₅𝜔²−𝑇₄+𝑇₅𝜔²+𝑇₃𝜔+𝑇₆−𝑎₆𝜔²²+𝑇₃+𝑇₅²𝜔²+𝑎²₁₂𝜔²+𝑇²₄+𝑇²₆+𝑇²₇. (27) Hence, the following theorem can be formulated.

Theorem 6. Let 𝑎₃+𝑎₁₃+|𝑎₁₄|+⋯+|𝑎₁₉|≠0 and (27) hold. Then there is no change in stability of 𝐸.

Remark 7. In the special case that 𝜉=𝜂, the characteristic equation (20) becomes 𝜒³+𝐾₄₁𝜒²+𝐾₅𝜒+𝐾₆+ℓ^−𝜒𝜏1𝑇₅₆𝜒²+𝑇₃₅𝜒+𝑇₄₆+ℓ^{−2𝜒𝜏1}𝑎₁₂𝜒+𝑇₇=0, (28) where 𝑇₅₆=𝑎₅+𝑎₆, 𝑇₃₅=𝑇₃+𝑇₅ and 𝑇₄₆=𝑇₄+𝑇₆. Therefore, this case is a private one of Theorem 6.

Corollary 8. If conditions of Theorem 6 are not valid and 𝜏₂^𝑏𝑖𝑓 is defined as in Theorem 3, then according to Lemma 5 for any 𝜏₂∈[0,𝜏_𝑏), there exists a 𝜏₃^𝑏𝑖𝑓(𝜏₂)>0(𝜏₁^𝑏𝑖𝑓(𝜏₂)>0 resp.) such that the steady state 𝐸 of system (2) is unstable when 𝜏₃∈[0,𝜏₃^𝑏𝑖𝑓(𝜏₂))(𝜏₁∈[0,𝜏₁^𝑏𝑖𝑓(𝜏₂)) resp.), and an Andronov-Hopf bifurcation takes place.

5. The General Case 𝜁,𝜉,𝜂>0

Similar to Section 4, we set that 𝜒=𝑖𝑛(𝑛>0) is a root of (4) if and only if 𝑛 satisfies −𝑛³𝑖−𝑎₁𝑛²+𝑎₂𝑛𝑖+𝑎₃+cos𝑛𝜏₁−𝑖sin𝑛𝜏₁−𝑎₄𝑛²+𝑎₇𝑛𝑖+𝑎₁₃+cos𝑛𝜏₂−𝑖sin𝑛𝜏₂−𝑎₅𝑛²+𝑎₈𝑛𝑖+𝑎₁₄+cos𝑛𝜏₃−𝑖sin𝑛𝜏₃−𝑎₆𝑛²+𝑎₉𝑛𝑖+𝑎₁₅+cos𝑛𝜏₄−𝑖sin𝑛𝜏₄𝑎₆𝑛𝑖+𝑎₁₆+cos𝑛𝜏₅−𝑖sin𝑛𝜏₅𝑎₁₁𝑛𝑖+𝑎₁₇+cos𝑛𝜏₆−𝑖sin𝑛𝜏₆𝑎₁₂𝑛𝑖+𝑎₁₈+𝑎₁₉cos𝑛𝜏₇−𝑖sin𝑛𝜏₇=0, (29) where 𝜏=[𝜏₁=𝜁, 𝜏₂=𝜉, 𝜏₃=𝜂, 𝜏₄=𝜁+𝜉, 𝜏₅=𝜁+𝜂, 𝜏₆=𝜉+𝜂, 𝜏₇=𝜁+𝜉+𝜂]^𝑇 denotes a point in the time delay space; that is, 𝜏∈Ω⊂𝑅⁷₊. Ω is the time delay space and 𝑅⁷₊ denotes the set of nonnegative real numbers.

Separating the real and imaginary parts into (29), we have −𝑎₁𝑛²+𝑎₃−𝐴₁−𝐵₁−𝐶₁=−𝐷₁−𝐸₁−𝐹₁−𝑎₁₉cos𝑛𝜏₇,−𝑛³+𝑎₂𝑛+𝐴₂+𝐵₂+𝐶₂=𝐷₂+𝐸₂+𝐹₂+𝑎₁₉sin𝑛𝜏₇, (30) where 𝐴₁=𝑎₄𝑛²−𝑎₁₃cos𝑛𝜏₁−𝑎₇𝑛sin𝑛𝜏₁,𝐵₁=𝑎₅𝑛²−𝑎₁₄cos𝑛𝜏₂−𝑎₈𝑛sin𝑛𝜏₂,𝐶₁=𝑎₆𝑛²−𝑎₁₅cos𝑛𝜏₃−𝑎₉𝑛sin𝑛𝜏₃,𝐷₁=𝑎₁₀𝑛sin𝑛𝜏₄+𝑎₁₆cos𝑛𝜏₄,𝐸₁=𝑎₁₁𝑛sin𝑛𝜏₅+𝑎₁₇cos𝑛𝜏₅,𝐹₁=𝑎₁₂𝑛sin𝑛𝜏₆+𝑎₁₈cos𝑛𝜏₆,𝐴₂=𝑎₄𝑛²−𝑎₁₃sin𝑛𝜏₁+𝑎₇𝑛cos𝑛𝜏₁,𝐵₂=𝑎₅𝑛²−𝑎₁₄sin𝑛𝜏₂+𝑎₈𝑛cos𝑛𝜏₂,𝐶₂=𝑎₆𝑛²−𝑎₁₅sin𝑛𝜏₃+𝑎₉𝑛cos𝑛𝜏₃,𝐷₂=𝑎₁₀𝑛cos𝑛𝜏₄−𝑎₁₆sin𝑛𝜏₄,𝐸₂=𝑎₁₁𝑛cos𝑛𝜏₅−𝑎₁₇sin𝑛𝜏₅,𝐹₂=𝑎₁₂𝑛cos𝑛𝜏₆−𝑎₁₈sin𝑛𝜏₆. (31) Adding up the squares of both equations into (30), we have 𝑛⁶+𝑆₁𝑛⁴+𝑆₂𝑛²+𝑆₃=2𝑎₄𝑛²−𝑎₁₃×𝑎₃−𝑎₁𝑛²cos𝑛𝜏₁−−𝑛³+𝑎₂𝑛sin𝑛𝜏₁+𝑎₇𝑛𝑎₃−𝑎₁𝑛²sin𝑛𝜏₁+−𝑛³+𝑎₂𝑛cos𝑛𝜏₁+𝑎₁₉𝑎₁₂𝑛sin𝑛𝜏₁+𝑎₁₈cos𝑛𝜏₁+𝑎₅𝑛²−𝑎₁₄×𝑎₃−𝑎₁𝑛²cos𝑛𝜏₂−−𝑛³+𝑎₂𝑛sin𝑛𝜏₂+𝑎₁₉𝑎₁₁𝑛sin𝑛𝜏₂+𝑎₁₇cos𝑛𝜏₂+𝑎₈𝑛𝑎₃−𝑎₁𝑛²sin𝑛𝜏₂+−𝑛³+𝑎₂𝑛cos𝑛𝜏₂+𝑎₁₉𝑎₁₀𝑛sin𝑛𝜏₃+𝑎₁₆cos𝑛𝜏₃+𝑎₆𝑛²−𝑎₁₅×𝑎₃−𝑎₁𝑛²cos𝑛𝜏₃−−𝑛³+𝑎₂𝑛sin𝑛𝜏₃+𝑎₉𝑛𝑎₃−𝑎₁𝑛²sin𝑛𝜏₃+−𝑛³+𝑎₂𝑛cos𝑛𝜏₃+𝑎₅𝑛²−𝑎₁₄𝑎₇𝑛sin𝑛𝜏₁−𝜏₂−𝑎₄𝑛²−𝑎₁₃cos𝑛𝜏₁−𝜏₂+𝑎₁₇𝑎₁₈cos𝑛𝜏₁−𝜏₂−𝑎₁₂𝑛sin𝑛𝜏₁−𝜏₂+𝑎₁₁𝑛𝑎₁₂𝑛cos𝑛𝜏₁−𝜏₂+𝑎₁₈sin𝑛𝜏₁−𝜏₂+𝑎₆𝑛²−𝑎₁₅𝑎₇𝑛sin𝑛𝜏₁−𝜏₃−𝑎₄𝑛²−𝑎₁₃cos𝑛𝜏₁−𝜏₃+𝑎₁₀𝑛𝑎₁₂𝑛cos𝑛𝜏₁−𝜏₃+𝑎₁₈sin𝑛𝜏₁−𝜏₃+𝑎₁₆𝑎₁₈cos𝑛𝜏₁−𝜏₃−𝑎₁₂𝑛sin𝑛𝜏₁−𝜏₃−𝑎₉𝑛𝑎₇𝑛cos𝑛𝜏₁−𝜏₃+𝑎₄𝑛²−𝑎₁₃sin𝑛𝜏₁−𝜏₃+𝑎₆𝑛²−𝑎₁₅𝑎₈𝑛sin𝑛𝜏₂−𝜏₃−𝑎₅𝑛²−𝑎₁₄cos𝑛𝜏₂−𝜏₃+𝑎₁₀𝑛𝑎₁₁𝑛cos𝑛𝜏₂−𝜏₃+𝑎₁₇sin𝑛𝜏₂−𝜏₃−𝑎₁₆𝑎₁₁𝑛sin𝑛𝜏₂−𝜏₃−𝑎₁₇cos𝑛𝜏₂−𝜏₃−𝑎₉𝑛𝑎₈𝑛cos𝑛𝜏₂−𝜏₃+𝑎₅𝑛²−𝑎₁₄sin𝑛𝜏₂−𝜏₃+4𝑎₁₉sin𝑛𝜏₁+𝜏₂𝑎₁₀𝑛cos𝑛𝜏₁+𝜏₂+𝜏₃−𝑎₁₆sin𝑛𝜏₁+𝜏₂+𝜏₃+sin𝑛𝜏₁+𝜏₃𝑎₁₁𝑛cos𝑛𝜏₁+𝜏₂+𝜏₃−𝑎₁₇sin𝑛𝜏₁+𝜏₂+𝜏₃+sin𝑛𝜏₂+𝜏₃𝑎₁₂𝑛cos𝑛𝜏₁+𝜏₂+𝜏₃−𝑎₁₈sin𝑛𝜏₁+𝜏₂+𝜏₃, (32) where 𝜏₄+𝜏₇=2𝜏₁+𝜏₂+𝜏₃,𝜏₅+𝜏₆=𝜏₁+𝜏₂+2𝜏₃,𝜏₅+𝜏₇=𝜏₂+2𝜏₁+𝜏₃,𝜏₆+𝜏₇=𝜏₁+2𝜏₂+𝜏₃,𝜏₄−𝜏₅=𝜏₂−𝜏₃,𝜏₄−𝜏₆=𝜏₁−𝜏₃,𝜏₅−𝜏₆=𝜏₁−𝜏₂,𝑆₁=𝑎²₁+𝑎²₄+𝑎²₅+𝑎²₆−2𝑎₂,𝑆₂=𝑎²₂+𝑎²₇+𝑎²₈+𝑎²₉−𝑎²₁₀−𝑎²₁₁−𝑎²₁₂,𝑆₃=𝑎²₃+𝑎²₁₃+𝑎²₁₄+𝑎²₁₅−𝑎²₁₆−𝑎²₁₇−𝑎²₁₈−𝑎²₁₉. (33) Clearly, the right-hand side of (32) is always less than 2𝑛³−𝑎₁𝑛²−𝑎₂𝑛+𝑎₃×𝑎₄𝑛²−𝑎₁₃+𝑎₅𝑛²−𝑎₁₄+𝑎₆𝑛²−𝑎₁₅+𝑎₁₉𝑎₁₂𝑛+𝑎₁₈+𝑎₁₁𝑛+𝑎₁₇+𝑎₁₀𝑛+𝑎₁₆+𝑛−𝑛³−𝑎₁𝑛²+𝑎₂𝑛+𝑎₃𝑎₇+𝑎₈+𝑎₉+𝑛𝑎₁₀+𝑎₁₁𝑎₁₂𝑛+𝑎₁₈−𝑎₁₆+𝑎₁₇𝑎₁₂𝑛−𝑎₁₈+−𝑎₄𝑛²+𝑎₇𝑛+𝑎₁₃𝑎₅𝑛²−𝑎₁₄+𝑎₆𝑛²−𝑎₁₅−𝑛𝑎₈+𝑎₉𝑎₄𝑛²+𝑎₇𝑛−𝑎₁₃−𝑎₉𝑛𝑎₅𝑛²+𝑎₈𝑛−𝑎₁₄+𝑎₆𝑛²−𝑎₁₅−𝑎₅𝑛²+𝑎₈𝑛+𝑎₁₄+𝑎₁₀𝑛𝑎₁₁𝑛+𝑎₁₇−𝑎₁₆𝑎₁₁𝑛−𝑎₁₇+4𝑎₁₉𝑎₁₀𝑛−𝑎₁₆+𝑎₁₁𝑛−𝑎₁₇+𝑎₁₂𝑛−𝑎₁₈. (34) Hence, if the inequality 𝜆⁶+𝑆₁𝜆⁴+𝑆₂𝜆²+𝑆₃>2𝜆³−𝑎₁𝜆²−𝑎₂𝜆+𝑎₃×𝑎₄𝜆²−𝑎₁₃+𝑎₅𝜆²−𝑎₁₄+𝑎₆𝜆²−𝑎₁₅+𝑎₁₉𝑎₁₂𝜆+𝑎₁₈+𝑎₁₁𝜆+𝑎₁₇+𝑎₁₀𝜆+𝑎₁₆+𝑎₇+𝑎₈+𝑎₉−𝜆³−𝑎₁𝜆²+𝑎₂𝜆+𝑎₃𝜆+𝑎₁₀+𝑎₁₁𝑎₁₂𝜆+𝑎₁₈𝜆−𝑎₁₆+𝑎₁₇𝑎₁₂𝜆−𝑎₁₈+−𝑎₄𝜆²+𝑎₇𝜆+𝑎₁₃𝑎₅𝜆²−𝑎₁₄+𝑎₆𝜆²−𝑎₁₅−𝑎₈+𝑎₉𝑎₄𝜆²+𝑎₇𝜆−𝑎₁₃𝜆−𝑎₉𝜆𝑎₅𝜆²+𝑎₈𝜆−𝑎₁₄+𝑎₆𝜆²−𝑎₁₅−𝑎₅𝜆²+𝑎₈𝜆+𝑎₁₄+𝑎₁₀𝜆𝑎₁₁𝜆+𝑎₁₇−𝑎₁₆𝑎₁₁𝜆−𝑎₁₇+4𝑎₁₉𝑎₁₀𝜆−𝑎₁₆+𝑎₁₁𝜆−𝑎₁₇+𝑎₁₂𝜆−𝑎₁₈ (35) has no real solution on 0<𝜆<𝑛₊, then (32) cannot be satisfied. Similar to Section 4, we note that 𝑛₊ is the positive solution of first equation in (6) (when 𝑚=0), which is written as 𝑎₁𝑛²(𝑛)=𝑎=Ψ₃−𝑎₄cos𝑛𝜁+𝑎₅cos𝑛𝜉+𝑎₆cos𝑛𝜂𝑛²+𝑎₇sin𝑛𝜁+𝑎₈sin𝑛𝜉+𝑎₉sin𝑛𝜂+𝑎₁₀sin𝑛(𝜁+𝜉)+𝑎₁₁sin𝑛𝜁+𝜂+𝑎₁₂sin𝑛𝜉+𝜂𝑛+𝑎₁₃cos𝑛𝜁+𝑎₁₄cos𝑛𝜉+𝑎₁₅cos𝑛𝜂+𝑎₁₆cos𝑛(𝜁+𝜉)+𝑎₁₇cos𝑛𝜁+𝜂+𝑎₁₈cos𝑛𝜉+𝜂+𝑎₁₉cos𝑛𝜁+𝜉+𝜂≤𝑎₃+𝑎₁₃+𝑎₁₄+𝑎₁₅+𝑎₁₆+𝑎₁₇+𝑎₁₈+𝑎₁₉+𝑎₇+𝑎₈+𝑎₉+𝑎₁₀+𝑎₁₁+𝑎₁₂𝑛−𝑎₄+𝑎₅+𝑎₆𝑛²𝑎,i.e.,₁+𝑎₄+𝑎₅+𝑎₆𝑛²−𝑎₇+𝑎₈+𝑎₉+𝑎₁₀+𝑎₁₁+𝑎₁₂𝑛−𝑎₃+𝑎₁₃+⋯+𝑎₁₉=0. (36) Thus, for 𝑛₊, we have 𝑛₊=12𝐴𝐵+𝐵²+4𝐴𝐶, (37) where 𝐴=𝑎₁+|𝑎₄|+|𝑎₅|+|𝑎₆|≠0, 𝐵=|𝑎₇|+⋯+|𝑎₁₂| and 𝐶=𝑎₃+|𝑎₁₃|+⋯+|𝑎₁₉|. It is clear that 𝑛≤𝑛₊.

Rearranging terms, we write (35) as 𝜆³−𝑎₁𝜆²−𝑎₂𝜆+𝑎₃−𝑎₄𝜆²−𝑎₁₃²+𝜆³−𝑎₁𝜆²−𝑎₂𝜆+𝑎₃−𝑎₅𝜆²−𝑎₁₄²+𝑎₁₁𝜆−𝑎₁₂𝜆+𝑎₁₈²+𝜆³−𝑎₁𝜆²−𝑎₂𝜆+𝑎₃−𝑎₆𝜆²−𝑎₁₅²+−𝜆³−𝑎₁𝜆²+𝑎₂𝜆+𝑎₃−𝜆𝑎₇+𝑎₈+𝑎₉²+𝑎₁₇+𝑎₁₂𝜆−𝑎₁₈²+𝑎₆𝜆²−𝑎₁₅−−𝑎₅𝜆²+𝑎₈𝜆+𝑎₁₄²+𝑎₉𝜆+𝑎₅𝜆²+𝑎₈𝜆−𝑎₁₄²+𝜆1+𝑎₁𝜆⁴+𝑎₂𝑎₃²+𝑎₁₆−𝑎₁₂𝜆−𝑎₁₈²+−𝑎₄𝜆²+𝑎₇𝜆+𝑎₁₃−𝑎₅𝜆²−𝑎₁₄²+−𝑎₄𝜆²+𝑎₇𝜆+𝑎₁₃−𝑎₆𝜆²−𝑎₁₅²+𝑎₁₀𝜆−𝑎₁₂𝜆+𝑎₁₈²+𝑎₄𝜆²+𝑎₇𝜆−𝑎₁₃+𝜆𝑎₈+𝑎₉²+𝜆𝑎₃+𝑎₁𝑎₂−𝜆²²+𝜆𝑎₁𝑎₃+𝑎₄𝑎₁₃+𝜆²+𝑎₁₀𝜆−𝑎₁₁𝜆+𝑎₁₇²+𝑎₁₆+𝑎₁₁𝜆−𝑎₁₇²+2𝑎₁₉−𝑎₁₀𝜆−𝑎₁₆²+2𝑎₁₉−𝑎₁₁𝜆−𝑎₁₇²+𝑎₅+𝑎₁₄𝜆²²+𝑎₆+𝑎₁₅𝜆²²+𝑎₈−𝑎₉𝜆²²+𝑎²₇𝜆²>2𝜆³−𝑎₁𝜆²−𝑎₂𝜆+𝑎₃²+𝑎₅𝜆²−𝑎₁₄²+2𝑎₆𝜆²−𝑎₁₅²+−𝜆³−𝑎₁𝜆²+𝑎₂𝜆+𝑎₃²+𝜆𝑎₇+𝑎₈+𝑎₉²+−𝑎₅𝜆²+𝑎₈𝜆+𝑎₁₄²+𝑎₅𝜆²+𝑎₈𝜆−𝑎₁₄²+𝜆𝑎₁𝜆⁴+𝑎₂𝑎₃²+𝑎₁𝑎₂+𝑎₃²+𝜆⁴+𝑎₁𝑎₃+𝑎₄𝑎₁₃²+𝜆²+2−𝑎₄𝜆²+𝑎₇𝜆+𝑎₁₃²+𝑎₄𝜆²+𝑎₇𝜆−𝑎₁₃²+2𝑎₁₀𝜆−𝑎₁₆²+2𝑎₁₁𝜆−𝑎₁₇²+𝑎₁₉+𝑎₁₂𝜆+𝑎₁₈²+𝑎₁₉+𝑎₁₁𝜆+𝑎₁₇²+𝑎₁₉+𝑎₁₀𝜆+𝑎₁₆²+2𝑎₁₉+𝑎₁₂𝜆−𝑎₁₈²+𝑎²₉+𝑎²₁₄+𝑎²₁₅𝜆⁴+𝑎²₉+3𝑎²₁₂𝜆²+𝑎²₅+𝑎²₆+𝑎²₈+𝑎²₁₀+𝑎²₁₆+2𝑎²₁₇+3𝑎²₁₈+2𝑎²₁₉. (38) Therefore, the following theorem can be formulated.

Theorem 9. Let 𝑎₃+|𝑎₁₃|+|𝑎₁₄|+⋯+|𝑎₁₉|≠0 and (38) hold. Then there is no change in stability of 𝐸.

Remark 10. In the special case that 𝜁=𝜉=𝜂, the characteristic equation (4) becomes 𝜒³+𝑎₁𝜒²+𝑎₂𝜒+𝑎₃+ℓ^−𝜒𝜂𝑎₄+𝑎₅+𝑎₆𝜒²+𝑎₇+𝑎₈+𝑎₉𝜒+𝑎₁₃+𝑎₁₄+𝑎₁₅+ℓ^−2𝜒𝜂𝑎₁₀+𝑎₁₁+𝑎₁₂𝜒+𝑎₁₆+𝑎₁₇+𝑎₁₈+𝑎₁₉ℓ^−3𝜒𝜂=0. (39)

Hence, this case is a private one of Theorem 9.

Corollary 11. If conditions of Theorem 9 are not valid and 𝜏₃^𝑏𝑖𝑓 is defined as in Theorem 3, then according to Lemma 5 for any 𝜏₃∈[0,𝜏_𝑏), there exists a 𝜏₁^𝑏𝑖𝑓(𝜏₃)>0(𝜏₂^𝑏𝑖𝑓(𝜏₃)>0; 𝜏₄^𝑏𝑖𝑓(𝜏₃)>0; 𝜏₅^𝑏𝑖𝑓(𝜏₃)>0; 𝜏₆^𝑏𝑖𝑓(𝜏₃)>0; 𝜏₇^𝑏𝑖𝑓(𝜏₃)>0 resp.) such that the steady state 𝐸 of system (2) is unstable when 𝜏₁∈[0,𝜏₁^𝑏𝑖𝑓(𝜏₃))(𝜏₂∈[0,𝜏₂^𝑏𝑖𝑓(𝜏₃)); 𝜏₄∈[0,𝜏₄^𝑏𝑖𝑓(𝜏₃)); 𝜏₅∈[0,𝜏₅^𝑏𝑖𝑓(𝜏₃)); 𝜏₆∈[0,𝜏₆^𝑏𝑖𝑓(𝜏₃))𝜏₇∈[0,𝜏₇^𝑏𝑖𝑓(𝜏₃)) 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: 𝑑𝑦₁=𝑘𝑑𝑡𝑘₂+𝑘₃𝑦^1𝑛1₂𝑡−𝜏₁−𝛾₁𝑦₁−𝑘₄𝑦₁𝑦₃,𝑑𝑦₂𝑑𝑡=𝑘₅𝑦₁𝑡−𝜏₂−𝛾₂𝑦₂,𝑑𝑦₃𝑑𝑡=𝑙−𝑘₆𝑦₃𝑡−𝜏₃−𝑘₄𝑦₁𝑦₃, (40) where the protein, 𝑦₂, controls its own synthesis through the repression of mRNA, 𝑦₁,  𝑦₃ is the concentration of miRNA, 𝜏₁ is the time delay for translation, 𝜏₂ is the time delay for transcription, 𝜏₃ is the average time delay for degradation of miRNA, 𝑛₁ is often referred as a Hill coefficient or a cooperativity coefficient, 𝑘_𝑖(𝑖=1,…,6), and 𝛾₁ and 𝛾₂ are the kinetic rate constants. According to [2, 24], we consider one representative value for Hill coefficient, 𝑛₁=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, 𝐸=(𝑦₁,𝑦₂,𝑦₃), represented by (40) can be analytically estimated and are defined by the following set of algebraic equations, including the rate constants of the model: 𝑦⁴₂+𝑘₅𝑘₄𝑙+𝑘₆𝛾₁𝑘₄𝛾₁𝛾₂𝑦³₂+𝑘₂𝑘₃𝑦²₂+𝑘₅𝑘₃𝑘₄𝛾₁𝛾₂𝑘₂𝑘₄𝑙+𝑘₆𝛾₁−𝑘₁𝑘₄𝑦₂−𝑘₁𝑘²₅𝑘₆𝑘₃𝑘₄𝛾₁𝛾²₂=0,𝑦₁=𝛾₂𝑘₅𝑦₂,𝑦₃=𝑘₅𝑙𝑘₅𝑘₆+𝛾₂𝑘₄𝑦₂. (41)

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 point 𝐸 of the system (40) defined as 𝑦_𝑘=𝑦_𝑘+𝑥_𝑘(𝑘=1,2,3). (42) In the case when 𝑛₁=2, the function 𝑘₁/(𝑘₂+𝑘₃𝑦²₂(𝑡−𝜏₁)) can be written as a MacLaurin series 𝑘₁𝑘₂+𝑘₃𝑦²₂𝑡−𝜏₁=𝑘₁𝛿+𝑘₃𝜑=𝑘₁𝛿𝑘₃=𝑘/𝛿𝜑+1₁𝛿𝑘1−³𝛿𝑘𝜑+³𝛿²𝜑²−𝑘₃𝛿³𝜑³+⋯, (43) where 𝛿=𝑘₂+𝑘₃𝑦²₂ and 𝜑=2𝑦₂𝑥₂(𝑡−𝜏₁)+𝑥²₂(𝑡−𝜏₁). If we take only linear term from (43) and after substitution of (42) into differential equation (40), we have, 𝑑𝑥₁𝑑𝑡=−𝑐₁𝑥₁−𝑐₂ℓ^−𝜒𝜏1𝑥₂−𝑐₃𝑥₃−𝑐₄ℓ^{−2𝜒𝜏1}𝑥²₂−𝑘₄𝑥₁𝑥₃,𝑑𝑥₂𝑑𝑡=𝑘₅ℓ^−𝜒𝜏2𝑥₁−𝛾₂𝑥₂,𝑑𝑥₃𝑑𝑡=−𝑐₅𝑥₁−𝑘₆ℓ^−𝜒𝜏3+𝑐₃𝑥₃−𝑘₄𝑥₁𝑥₃, (44) where 𝑐₁=𝛾₁+𝑘₄𝑦₃,𝑐₂=2𝑘₁𝑘₃𝑦₂𝛿²,𝑐₃=𝑘₄𝑦₁,𝑐₄=𝑘₁𝑘₃𝛿²,𝑐₅=𝑘₄𝑦₃. (45) It is seen that notations in Section 2 here are 𝑓₂=𝑓₃=𝑓₆=𝑔₁=𝑔₄=𝑔₅=𝑔₆=ℎ₂=ℎ₃=ℎ₄𝑓=0,₁=−𝑐₁,𝑓₄=−𝑐₂,𝑓₅=−𝑐₃,𝑔₂=𝑘₅,𝑔₃=−𝛾₂,ℎ₁=−𝑐₅,ℎ₅=−𝑐₆,ℎ₆=−𝑘₆. (46) 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 𝑘₁=0.3min⁻¹,𝑘₂=𝑘₃=1min⁻¹,𝑘₄=0.3min⁻¹,𝑘₅=0.5min⁻¹,𝑘₆∈[],𝛾0.05,0.3₁=0.1min⁻¹,𝛾₂=0.2min⁻¹,𝜏₁∈[1,8],𝜏₂∈[],𝜏12,35₃∈[10,300],𝑙=0.1. (47) 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 (𝑦₁), the concentration of protein (𝑦₂), and the concentration of miRNA (𝑦₃) are shown for absence of time delay (i.e., 𝜏₁=𝜏₂=𝜏₃=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 𝜏₃, 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)

(b)

Figure 1 

(a) Stable regime of system (40) for 𝜏1=𝜏2=𝜏3=0 and (b) unstable regime (sustained oscillations) for 𝜏1=1, 𝜏2=13, 𝜏3=60.

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.