# Dynamic Analysis and Hopf Bifurcation of a Lengyel–Epstein System with Two Delays

Authors: Long Li; Yanxia Zhang
Journal: Journal of Mathematics
 (2021)

---

## Abstract

In this paper, a Lengyel–Epstein model with two delays is proposed and considered. By choosing the different delay as a parameter, the stability and Hopf bifurcation of the system under different situations are investigated in detail by using the linear stability method. Furthermore, the sufficient conditions for the stability of the equilibrium and the Hopf conditions are obtained. In addition, the explicit formula determining the direction of Hopf bifurcation and the stability of bifurcating periodic solutions are obtained with the normal form theory and the center manifold theorem to delay differential equations. Some numerical examples and simulation results are also conducted at the end of this paper to validate the developed theories.


---

## Body

## 1. Introduction

It is of great significance to study the dynamical behaviors of chemical reaction models to understand the reaction mechanism and evolution law of the reaction process of reactants. Experiments show that the description of the internal reaction mechanism of the chemical reaction system by delay reaction-diffusion equation is more realistic. Therefore, many scholars have conducted in-depth studies on the equation and obtained a large number of practical conclusions. In 1990, De Kepper et al. discovered Turing patterns in the CIMA reaction experiment, which also verified Turing’s theoretical work over 40 years ago. In 1991 and 1992, Lengyel and Epstein in [1, 2] proposed the following mathematical model to depict the experimental process of the CIMA:(1)∂u∂t=Δu+a−u−4uv1+u2,t>0,x∈Ω,∂v∂t=σcΔv+bu−uv1+u2,t>0,x∈Ω,∂u∂n=∂v∂n=0,t>0,x∈∂Ω,ux,0=u0x≥0,vx,0=v0x≥0,x∈Ω,where ux,t and vx,t are the concentrations of two chemical reactants, respectively, a and b are positive constants representing the quantities related to the feed concentrations, the ratio of diffusion coefficients and the rescaling parameter are denoted as c and σ>0, respectively, Ω is a bounded open domain in ℝ with smooth boundary ∂Ω, and Δ is the Laplace operator.System (1) has been extensively studied by many scholars (see [3–7]). In [5], the researchers mainly analyzed some basic properties of system (1). It can be found through the study that when the parameter Ω (size of the reactor), d=c/b(effective diffusion rate), and a (initial concentration of the reactants) are relatively small, the system has no nonconstant steady state. In addition, when the initial concentration of the reactant a is within an appropriate range, the system has nonconstant steady states for a large diffusion rate, d=c/b. In 2005, on the basis of [3], Jang et al. studied the global bifurcation problem of the nonconstant positive steady-state solutions in one-dimensional case and considered its limiting behavior by using the shadow system method. In [6], Yi et al. discussed the conditions of homogeneous equilibrium solution and periodic solution and gave the detailed Hopf bifurcation analysis for the ODE and PDE models. In 2009, Yi et al. [7] proved the global asymptotical behavior of the constant positive equilibrium solution and convergence of all solutions of the system by constructing the Lyapunov function. In [4], Jin et al. proved the bifurcations of steady-state solutions and spatially nonhomogeneous periodic solutions of (1).Obviously, in the initial study of the model, most scholars did not consider the influence of time-delay factors. It is well known that the ordinary differential equation reflects that the development of things only depends on the current state, while the delay differential equation is used to describe the development system which depends on both the current state and the past state. In practical problems, if the delay factor is not considered, the properties and states of the system may change, even lead to wrong conclusions. Usually, a long delay will destroy the stability of the equilibrium point of the system. Therefore, it is more realistic to consider the delay factor in the system (see [8–15]). Based on this fact, Celik et al. [16] considered a coupled delayed-Lengyel–Epstein model as follows:(2)u˙t=a−ut−4utvt−τ1+u2t,t>0,v˙t=σbut−utvt−τ1+u2t,t>0.By choosingτ as the varying parameter and analyzing the related characteristic equations, Celik et al. [16] investigated the stability of the constant equilibrium point and the existence of Hopf bifurcation of system (2) and determined the necessary conditions for the parameters.In [17, 18], Merdan and Kayan proposed the following delayed models:(3)u˙t=a−ut−4ut−τvt1+u2t,t>0,v˙t=σbut−ut−τvt1+u2t,t>0,and(4)u˙t=a−ut−τ−4ut−τvt1+u2t−τ,t>0,v˙t=σbut−τ−ut−τvt1+u2t−τ,t>0.They investigated the existence of Hopf bifurcation of the systems and the properties of Hopf bifurcation.Recently, Zhang and He [19] assumed that the delay only occurred in the self-decomposition of the activator and further studied the delayed differential equation model as the following form:(5)u˙t=a−ut−τ−4utvt1+u2t,t>0,v˙t=σbut−utvt1+u2t,t>0.They analyzed the influence of the change of delayτ on the dynamical behaviors and found that the equilibrium of system (5) would eventually become unstable after passing through several stable switches and Hopf bifurcations at some certain critical values of τ.For a long time, most scholars are more concerned about biological mathematics. Therefore, in a period of time, biological mathematics has made a rapid development, and some practical conclusions have been obtained. However, there are many chemical reaction functional differential equations in nature, so it is of great practical significance to study the existence, stability, and Hopf bifurcation of periodic solutions of these functional differential equations in the field of chemistry. Of course, researchers did not realize the influence of time delay on the properties of periodic solutions of differential equations in the early stages. With the maturity of the biological mathematics theory, researchers began to realize the importance of time delay in chemical reaction systems. Apparently, these models (2)–(5) discussed above are based on the hypothesis that the delay effect is a single activator or inhibitor alone. In the field of chemistry, scholars are concerned about the range of time delay, when the system produces Hopf bifurcation and the equilibrium point is stable, which means that the parameters of the model should be selected in an appropriate range to save energy to the greatest extent and control environmental pollution for the model of the system. Many models represent the dynamic systems of chemical reactions through ordinary differential equations with no delay or differential equations with the same delay. However, the effects of time delay in the activator and inhibitor are completely different. Based on the fact, a more reasonable mathematical model to describe CIMA reaction should be depicted by the Lengyel–Epstein model with two delays as follows:(6)u˙t=a−ut−τ1−4ut−τ1vt−τ21+u2t−τ1,t>0,v˙t=σbut−τ2−ut−τ2vt−τ11+u2t−τ2,t>0,where τi≥0i=1,2 are the time lags of the activator and the inhibitor.We know that the most widely studied delay differential equations are the cases with a single delay or multiple identical delays. When the delay is taken as a parameter, the unique positive equilibrium of the system often loses stability at the critical value after Hopf bifurcation. However, there are multiple stable intervals and unstable intervals in the system when it has many different time delays because of their different critical values. We pay more attention to the changes happened to the dynamic properties of the system in this case. If we fix one of the time delays to its stable interval and take the other as a dynamic parameter, we can study the rich dynamic properties of the system. Since the time lags between catalysts and inhibitors are different in actual chemical reactions, it would be interesting and realistic to consider the two different delays in the system. To the best knowledge of the authors, the Lengyel–Epstein model with two different delays is seldom studied by scholars. By choosing different delays as the bifurcation parameter, the main purpose of this article is to analyze the effect of the two different delaysτi≥0i=1,2 on the dynamical behaviors of system (6).This article is organized as follows. In Section2, we investigate the stability of the positive equilibrium and the existence of local Hopf bifurcation of system (6) in three cases according to different symbols of τi. In Section 3, by choosing τ1 as a parameter and fixing τ2 in its stable interval, the direction and stability of bifurcating periodic solutions are obtained by using the normal form method and the center manifold reduction for delayed differential equations developed by Hassard et al. [20]. In Section 4, some numerical examples and simulations results are included to illustrate the validity of the main results.


## 2. Stability and Hopf Bifurcation

In this section, we use the method in [21] to discuss the system’s positivity and uniform boundedness.Lemma 1.
Letu0>0andv0>0. Then, the solutions utandvt of system (6) are nonnegative, for all t>0.Proof.
From the first equation, we have(7)dudt≥−u−4uv1+u2=−u1+4v1+u2≥−u1+4v,that is,(8)ut≥u0e−∫0t1+4vdω>0.
Therefore, one can show thatut>0, for all t>0.
Similarly, we can apply the same method to provevt>0, for all t>0. In fact, from the second equation, we can obtain(9)dvdt≥−σbuv1+u2,that is,(10)vt≥v0e−∫0tσbu/1+u2dω>0.
Therefore, one can show thatut>0,vt>0, for all t>0.
For the boundedness of all solutions of system (6), according to [5], we have the following result.Lemma 2.
Letu0>0andv0>0. Then, all solutions of system (6) are uniformly bounded, that is,(11)limsuputt⟶+∞≤a,limsupvtt⟶+∞≤1+a2.Now, we will analyze the stability of the positive equilibrium point and the existence of local Hopf bifurcation of system (6) in detail in three cases τ1=τ2=0,τ1=0andτ2>0, and τ1>0andτ2>0.Letα^=a/5; it is easy to get that system (6) has a unique positive equilibrium point E∗u∗,v∗=α^,1+α^2. Then, the linearized system of (6) is given as follows:(12)u˙t=P1ut−τ1−Q1vt−τ2,v˙t=P2ut−τ1−Q2vt−τ2,where(13)P1=3α^2−51+α^2,P2=2σbα^21+α^2,Q1=−4α^1+α^2,Q2=−σbα^1+α^2.The associated characteristic equation of the linear system (12) is given by the following form:(14)λ−P1e−λτ1−Q1e−λτ2−P2e−λτ2λ−Q2e−λτ1=0.That is,(15)λ2−P1+Q2λe−λτ1+P1Q2e−2λτ1−P2Q1e−2λτ2=0.Case 1.
τ1=τ2=0.
Then, the characteristic equation (15) can be simplified to(16)λ2−P1+Q2λ+P1Q2−P2Q1=0.
By the Routh–Hurwitz criterion, we can see that all roots of (16) have always negative real parts if the condition,(17)H1:P1+Q2<0,P1Q2−P2Q1>0,holds.
Hence, the positive equilibrium of system (6) is locally asymptotically stable if H1 holds.Case 2.
τ1=0,τ2>0.
Then, the associated characteristic equation (15) can be transformed to(18)λ2−P1+Q2λ+P1Q2−P2Q1e−2λτ2=0.
Letλ=iωω>0 be a root of the characteristic equation (18); then, we have(19)−ω2−iP1+Q2ω+P1Q2−P2Q1cos2ωτ2−isin2ωτ2=0.
Separating the real and imaginary parts of equation (19), we obtain(20)ω2−P1Q2=−P2Q1cos2ωτ2,P1+Q2ω=P2Q1sin2ωτ2.
It follows from (20) that(21)ω4+P12+Q22ω2+P12Q22−P22Q12=0.
Letz=ω2. Then, (21) becomes(22)z2+P12+Q22z+P12Q22−P22Q12=0.
Notice thatP12+Q22>0 and(23)Δ=P12+Q222−4P12Q22−P22Q12=P12−Q222+4P22Q12>0.
Hence, one can see that if the condition(24)H2:P1Q22−P2Q12>0,holds, then equation (22) has no positive solution. Therefore, all solutions of equation (18) have negative real parts when τ2>0, and condition H2 holds.Theorem 1.
Assume thatτ1=0 and H2 holds; then, the positive equilibrium E∗ of system (6) is asymptotically stable, for all τ2>0.Now, we suppose that the following condition is given:(25)H3:P1Q22−P2Q12<0.Then, one can obtain thatz0 is a unique positive root of equation (22) and ±iω0 is a pair of purely imaginary roots of the characteristic equation (18). Here,(26)ω0=z0=12−P12+Q22+P12−Q222+4P22Q12.By substitutingω0 into (20), one can solve a series of τ2 as(27)τ2j=12ω0arccosω02−P1Q2−P2Q1+2jπω0,j=0,1,2,….Letλτ2=ατ2+iωτ2 be a root of (18) near τ2=τ2j such that ατ2j=0 and ωτ2j=ω0. Differentiating both sides of equation (18) with respect to τ2, the following can be performed:(28)2λdλdτ2−P1+Q2dλdτ2−P2Q1e−2λτ2−2τ2dλdτ2−2λ=0.It follows that(29)dλdτ2−1=2λ−P1+Q2+2τ2P2Q1e−2λτ2−2λP2Q1e−2λτ2=2λ−P1+Q2e2λτ2−2λP2Q1−τ2λ.Furthermore, we can know that(30)dReλτ2dτ2τ2=τ2j−1=Re2λ−P1+Q2e2λτ2−2λP2Q1|τ2=τ2j=2ω02+P12+Q222P2Q12>0.According to the above discussion and Corollary 2.4 of [3], we can present the following results.Theorem 2.
Assume thatτ1=0 and (H3) hold, τ2jj=0,1,2,… can be defined by (27), and ω0 can be given by (26). Then,(i)
The positive equilibriumE∗ of system (6) is asymptotically stable, for all τ2∈0,τ20.(ii)
The positive equilibriumE∗ of system (6) is unstable for τ2∈τ20,+∞(iii)
System (6) undergoes Hopf bifurcation at the positive equilibrium E∗ for τ2=τ20Case 3.
τ1>0,τ2>0.
In this case, we consider the characteristic equation (15) by choosing τ1 as a parameter and fixing τ2 in its stable interval 0,τ20. Without loss of generality, we analyze system (6) under condition H3.
Multiplying both sides of equation (15) by eλτ1, then one can get the following equation:(31)λ2eλτ1−P1+Q2λ+P1Q2e−λτ1−P2Q1e−2λτ2e−λτ1=0.
Letλ=iωω>0 be a root of the characteristic equation (15); then, equation (31) becomes the following form:(32)−ω2cosωτ1+isinωτ1−iP1+Q2ω+P1Q2cosωτ1−isinωτ1−P2Q1cos2ωτ2−isin2ωτ2cosωτ1+isinωτ1=0.
Separating the real and imaginary parts of equation (32), we obtain(33)−ω2+P1Q2−P2Q1cos2ωτ2cosωτ1−P2Q1sin2ωτ2sinωτ1=0,−ω2−P1Q2−P2Q1cos2ωτ2sinωτ1+P2Q1sin2ωτ2cosωτ1=P1+Q2ω.
It follows from (33) that(34)sinωτ1=P1+Q2ω−ω2+P1Q2−P2Q1cos2ωτ2−ω2−P2Q1cos2ωτ22−P1Q22+P2Q1sin2ωτ22,cosωτ1=P1+Q2ωP2Q1sin2ωτ2−ω2−P2Q1cos2ωτ22−P1Q22+P2Q1sin2ωτ22.
According to (34), we can obtain(35)ω8+l1ω6+l2ω4+l3ω2+l4=0,where(36)l1=4P2Q1cos2ωτ2−P1+Q22,l2=4P2Q1cos2ωτ22+2P2Q12−P1Q22+2P1+Q22P2Q1cos2ωτ2,l3=4P2Q1P2Q12−P1Q22cos2ωτ2−P1+Q22P1Q22+P2Q12+2P1+Q22P1Q2P2Q1cos2ωτ2,l4=P2Q12−P1Q222.
Now, we give the following assumption for equation (35).
H4: assume that equation (35) has at least one positive root.
If conditionH4 holds, then equation (35) has a positive root ω∗ such that equation (15) has a pair of purely imaginary roots ±ω∗. Thus, we have(37)τ1j=1ω∗arccosP1+Q2ω∗2P2Q1sin2ω∗τ2−ω∗2−P2Q1cos2ω∗τ22−P1Q22+P2Q1sin2ω∗τ22+2jπω∗,j∈ℕ0.
Differentiating both sides of equation (15) with respect to τ1, the following equation can be easily obtained:(38)2λdλdτ1−P1+Q2dλdτ1e−λτ1+λe−λτ1−dλdτ1τ1−λ+P1Q2e−2λτ1−2dλdτ1τ1−2λ−P2Q1e−2λτ2−2τ2dλdτ1=0.
It follows that(39)dλdτ1−1=2λ−P1+Q2e−λτ1+P1+Q2λτ1e−λτ1−2P1Q2τ1e−2λτ1+2P2Q1τ2e−2λτ2−P1+Q2λ2e−λτ1−2λP1Q2e−2λτ1=2λ−P1+Q2e−λτ1+2P2Q1τ2e−2λτ2−P1+Q2λ2e−λτ1+2λP1Q2e−2λτ1−τ1λ=2λeλτ1−P1+Q2+2P2Q1τ2e−2λτ2eλτ1−P1+Q2λ2+2λP1Q2e−λτ1−τ1λ.
Letλ=iω∗, and we can know that(40)dReλτ1dτ1τ1=τ1j−1=M1ω∗N1ω∗+M2ω∗N2ω∗M12ω∗+M22ω∗,where(41)M1ω∗=P1+Q2ω∗2+2P1Q2ω∗sinω∗τ1,M2ω∗=2P1Q2ω∗cosω∗τ1,N1ω∗=2P2Q1τ2cos2ω∗τ2cosω∗τ1−2ω∗−2P2Q1τ2sin2ω∗τ2sinω∗τ1−P1+Q2,N2ω∗=2P2Q1τ2cos2ω∗τ2sinω∗τ1+2ω∗−2P2Q1τ2sin2ω∗τ2cosω∗τ1.
Next, we assume that the following condition holds:(42)H5:M1ω∗N1ω∗+M2ω∗N2ω∗≠0.
Therefore, from the above discussions and Hopf bifurcation Theorem2, the following results can be directly deduced.Theorem 3.
Assume that (H3), (H4), and (H5) hold, τ2∈0,τ20 and τ1jj∈ℕ0 can be defined by (37), and ω∗ can be given by (35). Then,(i)
The positive equilibriumE∗ of system (6) is asymptotically stable for all τ1∈0,τ10(ii)
The positive equilibriumE∗ of system (6) is unstable for τ1∈τ10,+∞(iii)
System (6) undergoes Hopf bifurcation at the positive equilibrium E∗ for τ1=τ10


## 3. Direction of Hopf Bifurcation and Stability of Bifurcating Periodic Solution

In the previous section, the stability of the positive equilibriumE∗ and the existence of Hopf bifurcation of system (6) by regarding different delays τii=1,2 as the bifurcation parameter were discussed. In this section, by applying the normal form theory and the center manifold theorem of Hassard et al. [20], we derive the explicit formulas, determining the direction of Hopf bifurcation and the stability of bifurcation periodic solutions. We always assume that system (6) undergoes Hopf bifurcation at the positive equilibrium E∗ for τ1=τ10, and the corresponding purely imaginary roots of the characteristic equation at E∗ are denoted by ±iω∗.Let us assume thatτ2∗<τ10, where τ2∈0,τ20 and τ20 is defined by (27). For the convenience of discussion, let u¯t=ut−u∗, v¯t=vt−v∗, t=τt¯, and τ1=τ10+μ,μ∈ℝ; then, μ=0 is the Hopf bifurcation value of (6). We drop the bars of these notations to simplify the next discussion, and system (6) can be concisely regarded as the abstract FDE in C=C−1,0,ℝ2:(43)u˙t=Lμut+Fμ,ut,where ut=u1t,u2tT, utθ=ut+θ, Lμ:C⟶ℝ2, and Fμ,⋅:ℝ×C⟶ℝ2 are, respectively, the linear operator and the nonlinear operator with the following forms:(44)Lμϕ=τ10+μAϕ1−τ2∗τ10ϕ2−τ2∗τ10+Bϕ1−1ϕ2−1,Fμ,ϕ=τ10+μf1,f2T,where ϕθ=ϕ1θ,ϕ2θT∈C, A=0Q1P20, B=P100Q2, and(45)f1=4α3−α21+α22ϕ12−1+−41−α21+α22ϕ1−1ϕ2−τ2∗τ10+⋯,f2=σbα3−α21+α22ϕ12−τ2∗τ10+−σb1−α21+α22ϕ1−τ2∗τ10ϕ2−1+⋯.According to Riesz representation theorem, there exists a matrix functionηθ,μ,θ∈−1,0 with bounded variation components such that(46)Lμϕ=∫−10dηθ,μϕθ,ϕθ∈C−1,0,ℝ2.In fact, we can take(47)ηθ,μ=τ10+μA+B,θ=0,τ10+μA+B,θ∈−τ2∗τ10,0,τ10+μB,θ∈−1,−τ2∗τ10,0.θ=−1.In the following discussion, we define the linear differential operatorsAμ and Rμ by(48)Aμϕ=dϕθdθ,θ∈−1,0,∫−10dηs,μϕs,θ=0,Rμϕ=0,θ∈−1,0,Fμ,ϕ,θ=0.Then, system (43) can be further represented as(49)u˙t=Aμut+Rμut,where utθ=ut+θ,θ∈−1,0. For ψ∈C10,1,ℝ2∗, the linear differential operator A∗ is defined by(50)A∗ψs=−dψsds,s∈0,1,∫−10dηTt,0ψ−t,s=0.Furthermore, forϕ∈C1−1,0,ℝ2∗ and ψ∈C10,1,ℝ2∗, we give the definition of the bilinear inner product:(51)ψs,ϕθ=ψ¯0ϕ0−∫−10∫ξ=0θψ¯ξ−θdηθϕξdξ,where ηθ=ηθ,0. Obviously, A0 and A∗ are adjoint operators. It is easy to see that they are eigenvalues of the linear operator A0. It follows that they are also eigenvalues of the linear operator A∗.Assume thatqθ=1,αTeiω∗τ10θ is the eigenvector of the operator A0 corresponding to the eigenvalues iω∗τ10; then, A0qθ=iω∗τ10qθ. From the definitions of A, Lμϕ, and ηθ,μ, we can obtain(52)iω∗−P1e−iω∗τ10−Q1e−iω∗τ2∗−P2e−iω∗τ2∗iω∗−Q2e−iω∗τ101α=00.By simple calculation, we have(53)α=iω∗−P1e−iω∗τ10Q1e−iω∗τ2∗.By analogy, we can get the eigenvectorq∗s=D1,α∗eiω∗τ10s of the operator A∗ corresponding to the eigenvalues −iω∗τ10, where(54)α∗=iω∗+P1e−iω∗τ10P2e−iω∗τ2∗.Now, we evaluate the value ofD such that q∗s,qθ=1. From the bilinear inner product of (51), it follows that(55)q∗s,qθ=q¯∗0q0−∫−10∫ξ=0θq¯∗ξ−θdηθqξdξ=D¯1,α¯∗1,αT−∫−10∫ξ=0θD¯1,α¯∗e−iω∗τ10ξ−θdηθ1,αTeiω∗τ10ξdξ=D¯1,α¯∗1,αT−D¯1,α¯∗∫−10∫ξ=0θdηθ1,αTeiω∗τ10θdξ=D¯1,α¯∗1,αT−D¯1,α¯∗∫−10θeiω∗τ10θdηθ1,αT=D¯1,α¯∗1,αT−D¯1,α¯∗τ10Aϕ−τ2∗τ10+τ10Bϕ−11,αT=D¯1,α¯∗1,αT−D¯1,α¯∗τ100Q1P20−τ2∗τ10e−iω∗τ2∗−τ10P100Q2e−iω∗τ101,αT=D¯1+αα¯∗+D¯τ10P1e−iω∗τ10+P2α¯∗τ2∗τ10e−iω∗τ2∗+αQ1τ2∗τ10e−iω∗τ2∗+Q2αα¯∗e−iω∗τ10=D¯1+αα¯∗+D¯P1τ10e−iω∗τ10+P2α¯∗τ2∗e−iω∗τ2∗+αQ1τ2∗e−iω∗τ2∗+Q2τ10αα¯∗e−iω∗τ10.Thus, we have(56)D=11+α¯α∗+P1τ10eiω∗τ10+P2α∗τ2∗eiω∗τ2∗+αQ1τ2∗eiω∗τ2∗+Q2τ10α¯α∗eiω∗τ10.In addition, fromψ,Aϕ=A∗ψ,ϕ and Aq¯θ=−iω∗q¯θ, we can obtain(57)−iω∗q∗,q¯=q∗,Aq¯=A∗q∗,q¯=−iω∗q∗,q¯=iω∗q∗,q¯.Hence,q∗θ,q¯θ=0.In the next, we calculate the coordinates describing the center manifoldC0 at μ=0 by using the method of Hassard et al. [20]. Let ut be the solution of equation (43) and zt=q∗,ut; then, from (49) and (51), we have(58)z˙t=q∗,u˙t=q∗,A0ut+R0ut=q∗,A0ut+q∗,R0ut=A∗0q∗,ut+q¯∗0F0,ut=iω∗τ10z+gz,z¯,where(59)gz,z¯=q¯∗0F0,ut=g20θz22+g11θzz¯+g02θz¯22+g21θz2z¯2+,….Assume(60)Wt,θ=utθ−ztqθ−z¯tq¯θ=utθ−2Reztqθ.Then, we have(61)Wt,θ=Wzt,z¯t,θ=W20z22+W11zz¯+W02z¯22+W30z36+⋯,where z and z¯ refer to local coordinates on the center manifold C0 in the directions of q∗ and q¯∗. As W is real, if ut is real, we only discuss real solutions. From (60), it follows that(62)utθ=u1tθ,u2tθT=Wt,θ+2Reztqθ=Wt,θ+ztqθ+z¯tq¯θ=1,αTeiω∗τ10θz+1,α¯Te−iω∗τ10θz¯+W20z22+W11zz¯+W02z¯22+⋯.Therefore,(63)u1t−1=ze−iω∗τ10+z¯eiω∗τ10+W201−1z22+W111−1zz¯+W021−1z¯22+⋯,u2t−1=zαe−iω∗τ10+z¯α¯eiω∗τ10+W202−1z22+W112−1zz¯+W022−1z¯22+⋯,u1t−τ2∗τ10=ze−iω∗τ2∗+z¯eiω∗τ2∗+W201−τ2∗τ10z22+W111−τ2∗τ10zz¯+W021−τ2∗τ10z¯22+⋯,u2t−τ2∗τ10=zαe−iω∗τ2∗+z¯α¯eiω∗τ2∗+W202−τ2∗τ10z22+W112−τ2∗τ10zz¯+W022−τ2∗τ10z¯22+⋯.Furthermore, we have(64)f1=4α3−α21+α22ϕ12−1+−41−α21+α22ϕ1−1ϕ2−τ2∗τ10+⋯,f2=σbα3−α21+α22ϕ12−τ2∗τ10+−σb1−α21+α22ϕ1−τ2∗τ10ϕ2−1+⋯,gz,z¯=q¯∗0F0,ut=D¯1,α¯∗F0,ut=D¯1,α¯∗τ10α1u1t2−1+α2u1t−1u2t−τ2∗τ10β1u1t2−τ2∗τ10+β2u1t−τ2∗τ10u2t−1=D¯τ10α1ze−iω∗τ10+z¯eiω∗τ10+W201−1z22+W111−1zz¯+W021−1z¯22+⋯2+α2ze−iω∗τ10+z¯eiω∗τ10+W201−1z22+W111−1zz¯+W021−1z¯22+⋯⋅zαe−iω∗τ2∗+z¯α¯eiω∗τ2∗+W202−τ2∗τ10z22+W112−τ2∗τ10zz¯+W022−τ2∗τ10z¯22+⋯+α¯∗β1ze−iω∗τ2∗+z¯eiω∗τ2∗+W201−τ2∗τ10z22+W111−τ2∗τ10zz¯+W021−τ2∗τ10z¯22+⋯2+α¯∗β2ze−iω∗τ2∗+z¯eiω∗τ2∗+W201−τ2∗τ10z22+W111−τ2∗τ10z¯+W021−τ2∗τ10z¯22+⋯⋅zαe−iω∗τ10+z¯α¯eiω∗τ10+W202−1z22+W112−1zz¯+W022−1z¯22+⋯=D¯τ10α1e−2iω∗τ10+α2αe−iω∗τ10e−iω∗τ2∗+α¯∗β1e−2iω∗τ2∗+α¯∗β2αe−iω∗τ10e−iω∗τ2∗z2+2α1+α2α¯e−iω∗τ10eiω∗τ2∗+α2αeiω∗τ10e−iω∗τ2∗+2α¯∗β1+α¯∗β2α¯eiω∗τ10e−iω∗τ2∗+αe−iω∗τ10eiω∗τ2∗zz¯+α1e2iω∗τ10+α2α¯eiω∗τ10eiω∗τ2∗+α¯∗β1e2iω∗τ2∗+α¯∗β2α¯eiω∗τ10eiω∗τ2∗z¯2+α12W111−1+W201−1e−iω∗τ10+α2W112−τ2∗τ10e−iω∗τ10+12W202−τ2∗τ10α¯eiω∗τ10+12W201−1α¯eiω∗τ2∗+W111−1αe−iω∗τ2∗+α¯∗β12W111−τ2∗τ10e−iω∗τ2∗+W201−τ2∗τ10eiω∗τ2∗+α¯∗β2W112−1e−iω∗τ2∗+12W202−1eiω∗τ2∗+12W201−τ2∗τ10α¯eiω∗τ10+W111−τ2∗τ10αe−iω∗τ10z2z¯+⋯,where(65)α1=4α3−α21+α22,α2=−41−α21+α22,β1=σbα3−α21+α22,β2=−σb1−α21+α22.According to the above analysis, we can obtain(66)g20=2D¯τ10α1e−2iω∗τ10+α2αe−iω∗τ10e−iω∗τ2∗+α¯∗β1e−2iω∗τ2∗+α¯∗β2αe−iω∗τ10e−iω∗τ2∗,g11=2D¯τ10α1+α2Reαeiω∗τ10e−iω∗τ2∗+α¯∗β1+α¯∗β2Reαe−iω∗τ10eiω∗τ2∗;g02=2D¯τ10α1e2iω∗τ10+α2α¯eiω∗τ10eiω∗τ2∗+α¯∗β1e2iω∗τ2∗+α¯∗β2α¯eiω∗τ10eiω∗τ2∗,g21=2D¯τ10α12W111−1+W201−1e−iω∗τ10+α2W112−τ2∗τ10e−iω∗τ10+12W202−τ2∗τ10α¯eiω∗τ10+12W201−1α¯eiω∗τ2∗+W111−1αe−iω∗τ2∗+α¯∗β12W111−τ2∗τ10e−iω∗τ2∗+W201−τ2∗τ10eiω∗τ2∗+α¯∗β2W112−1e−iω∗τ2∗+12W202−1eiω∗τ2∗+12W201−τ2∗τ10α¯eiω∗τ10+W111−τ2∗τ10αe−iω∗τ10.Afterg20, g21, and g02 are obtained, W20θ and W11θ in g21 can be determined by using method in [22–24]. From (49) and (60), we can obtain(67)W˙=u˙t−z˙q−z¯˙q¯=AW−2Req¯∗0F0,utqθθ∈−1,0AW−2Req¯∗0F0,utq0+F0,utθ=0≜AW+Hz,z¯,θ,where(68)Hz,z¯,θ=H20θz22+H11θzz¯+H02θz¯22+⋯.Thus, we have(69)AWt,θ−W˙=−Hz,z¯,θ=−H20θz22−H11θzz¯−H02θz¯22−⋯.From (61), one can obtain(70)AWt,θ=AW20θz22+AW11θzz¯+AW02θz¯22+AW30θz36+⋯,W˙=Wzz˙+Wz¯z¯˙=W20θzz˙+W11θz˙z¯+zz¯˙+⋯=2iω∗τ10W20θz22+⋯.Therefore, we can further obtain(71)A−2iω∗τ10W20θ=−H20θ,AW11θ=−H11θ.Forθ∈−1,0, we have(72)Hz,z¯,θ=−q¯∗0F0,utqθ−q∗0F¯0,utq¯θ=−gz,z¯qθ−g¯z,z¯q¯θ.A comparison of coefficients for (72) and (68) leads to(73)H20θ=−g20qθ+g¯02q¯θ,(74)H11θ=−g11qθ+g¯11q¯θ.By (71) and (73) and the definition of A, one can obtain(75)W˙20θ=2iω∗τ10W20θ+g20qθ+g¯02q¯θ.Notice thatqθ=q0eiω∗τ10θ=1,αTeiω∗τ10θ. By computing (75), we have(76)W20θ=ig20ω∗τ10q0eiω∗τ10θ+ig¯023ω∗τ10q¯0e−iω∗τ10θ+E1e2iω∗τ10θ,where E1=E11,E12∈ℝ2 is also a two-dimensional constant vector.Similarly, by (71) and (74) and the definition of A, one can obtain(77)W˙11θ=g11qθ+g¯11q¯θ,(78)W11θ=−ig11ω∗τ10q0eiω∗τ10θ+ig¯11ω∗τ10q¯0e−iω∗τ10θ+E2,where E2=E21,E22∈ℝ2 is a two-dimensional constant vector.Next, we need to computeE1 and E2 in (76) and (78). By the definition of A and (71), we obtain that(79)∫−10dηθW20θ=2iω∗τ10W200−H200,(80)∫−10dηθW11θ=−H110.Thus, we have(81)H200=−g20q0−g¯02q¯0+2τ10α1e−2iω∗τ10+α2αe−iω∗τ10e−iω∗τ2∗β1e−2iω∗τ2∗+β2αe−iω∗τ10e−iω∗τ2∗,(82)H110=−g11q0−g¯11q¯0+2τ10α1+α2Reαeiω∗τ10e−iω∗τ2∗β1+β2Reαe−iω∗τ10eiω∗τ2∗.Noting that(83)−iω∗τ10I−∫−10e−iω∗τ10θdηθq¯0=0.Substituting (76) and (81) into equation (80) yields(84)2iω∗τ10I−∫−10e2iω∗τ10θdηθE1=2τ10α1e−2iω∗τ10+α2αe−iω∗τ10e−iω∗τ2∗β1e−2iω∗τ2∗+β2αe−iω∗τ10e−iω∗τ2∗.That is,(85)2iω∗−P1e−2iω∗τ10−Q1e−2iω∗τ2∗−P2e−2iω∗τ2∗2iω∗−Q2e−2iω∗τ10E1=2α1e−2iω∗τ10+α2αe−iω∗τ10e−iω∗τ2∗β1e−2iω∗τ2∗+β2αe−iω∗τ10e−iω∗τ2∗.It follows that(86)E11=2Δ1α1e−2iω∗τ10+α2αe−iω∗τ10e−iω∗τ2∗−Q1e−2iω∗τ2∗β1e−2iω∗τ2∗+β2αe−iω∗τ10e−iω∗τ2∗2iω∗−Q2e−2iω∗τ10,E12=2Δ12iω∗−P1e−2iω∗τ10α1e−2iω∗τ10+α2αe−iω∗τ10e−iω∗τ2∗−P2e−2iω∗τ2∗β1e−2iω∗τ2∗+β2αe−iω∗τ10e−iω∗τ2∗,where(87)Δ1=2iω∗−P1e−2iω∗τ10−Q1e−2iω∗τ2∗−P2e−2iω∗τ2∗2iω∗−Q2e−2iω∗τ10.Similarly, substituting (78) and (82) into (80) yields(88)∫−10dηθE2=−2τ10α1+α2Reαeiω∗τ10e−iω∗τ2∗β1+β2Reαe−iω∗τ10eiω∗τ2∗.That is,(89)−P1−Q1−P2−Q2E2=2α1+α2Reαeiω∗τ10e−iω∗τ2∗β1+β2Reαe−iω∗τ10eiω∗τ2∗.It follows that(90)E21=2Δ2α1+α2Reαeiω∗τ10e−iω∗τ2∗−Q1β1+β2Reαe−iω∗τ10eiω∗τ2∗−Q2,E22=2Δ2−P1α1+α2Reαeiω∗τ10e−iω∗τ2∗−P2β1+β2Reαe−iω∗τ10eiω∗τ2∗,where(91)Δ2=P1Q1P2Q2.Thus, from (76), (78), (86), and (90), we can determine g21 and compute the following values:(92)c10=i2ω∗τ10g20g11−2g112−13g022+12g21,μ2=−Rec10Reλ′τ10,β2=2Rec10,T2=−1ω∗τ10Imc10+μ2Imλ′τ10.Now, the main results in this section are given.Theorem 4.
Ifμ2>0 (resp.μ2<0), then the periodic solution is supercritical (resp.subcritical); if β2<0 (resp.β2>0), then the bifurcating periodic solutions are orbitally asymptotically stable with an asymptotical phase (resp. unstable); if T2>0 (resp.T2<0), then the period of the bifurcating periodic solutions increases (resp. decreases).


## 4. Examples and Numerical Simulations

In this section, by applying the Matlab software package, we shall give some examples and numerical simulations to verify the results developed in the previous sections.Example 1.
Leta=3.3andσb=2; then, system (6) becomes(93)u˙t=3.3−ut−τ1−4ut−τ1vt−τ21+u2t−τ1,v˙t=2ut−τ2−ut−τ2vt−τ11+u2t−τ2.
By calculating, system (93) has a unique positive equilibrium u∗,v∗=0.66,1.4356. When τ1=0, it is easily obtained that the condition H2:P1Q22−P2Q12=0.6136>0 is satisfied. By Theorem 1, the positive equilibrium E∗ of (93) is asymptotically stable for any τ2>0, as shown in Figures 1 and 2.Figure 1
The trajectories graphs and phase graphs of system (93) with τ1=0andτ2=1 at the initial value 1,1.2.Figure 2
The trajectories graphs and phase graphs of system (93) with τ1=0andτ2=20 at the initial value 1,1.2.Example 2.
Leta=5.5andσb=2; we have the following system:(94)u˙t=5.5−ut−τ1−4ut−τ1vt−τ21+u2t−τ1,v˙t=2ut−τ2−ut−τ2vt−τ11+u2t−τ2.
By calculation with Matlab, system (94) has a unique positive equilibrium u∗,v∗=1.1,2.21. When τ1=τ2=0, we obtain P1+Q2=−1.6154<0andP1Q2−P2Q1=4.9774>0. That is, condition H1 holds. Thus, the positive equilibrium E∗ of (94) is asymptotically stable, as depicted in Figure 3. When τ1=0,τ2>0, we obtain that condition H3:P1Q22−P2Q12=−18.6311<0 holds and ω0=1.9192,τ20=0.2061,dReλτ2/dτ2λ=iw20,τ2=τ20−1=0.2299>0. It is known from Theorem 2, the positive equilibrium E∗ of system (94) is asymptotically stable when τ2∈0,τ20. However, τ2 passes through the critical value τ20, the positive equilibrium E∗ loses its stability, and a Hopf bifurcation occurs, as shown in Figures 4 and 5.
Letτ2=0.15∈0,0.2061, and choose τ1>0 as a parameter. The numerical results show that ω∗=2.3361andτ10=0.2107 and the conditions H3:P1Q22−P2Q12=−18.6311<0,H4:ω∗=2.3361>0,andH5:M1ω∗N1ω∗−M2ω∗N2ω∗=0.4023>0 are satisfied. By Theorem 3, the positive equilibrium E∗ of system (94) is asymptotically stable. And, the transversality condition dReλτ1/dτ1λ=iw∗,τ1=τ10−1=0.0065>0 is also satisfied. Therefore, when τ1 passes through the critical value τ10, the positive equilibrium E∗ loses its stability and a family of periodic solutions bifurcates from the positive equilibrium E∗, as shown in Figures 6 and 7. From formula (92) in Section 3, it follows that c10=−1.1543+7.05i,μ2=1.7611,β2=−2.3087, and T2=−22.7964. Since μ2>0andβ2<0, the direction of Hopf bifurcation of system (94) at τ10=0.2107 is supercritical, and these bifurcating periodic solutions are stable, as shown in Figure 7.Figure 3
The trajectories graphs and phase graphs of system (94) with τ1=τ2=0 at the initial value 0,0.Figure 4
The trajectories graphs and phase graphs of system (94) with τ1=0,τ2=0.15<τ20=0.2061 at the initial value 1,1.Figure 5
The trajectories graphs and phase graphs of system (94) with τ1=0,τ2=0.23>τ20=0.2061 at the initial value 1,1.Figure 6
The trajectories graphs and phase graphs of system (94) with τ1=0.2<τ10=0.2107,τ2=0.15∈0,0.2061 at the initial value 1,1.Figure 7
The trajectories graphs and phase graphs of system (94) with τ1=0.23>τ10=0.2107,τ2=0.15∈0,0.2061 at the initial value 1,1.


## 5. Conclusions

Researchers often apply differential equations with parametric-to-dynamic system modeling because the solution of the equations depends on the behaviors of the parameters, which can effectively help researchers understand the system. In the past few decades, the dynamic behavior of the time-delay model has been widely concerned in many fields (chemistry, economics, medicine, ecology, etc.), which is inseparable from its biological and physical significance. The research on the dynamics of delay differential equations has always been the focus of attention in various fields. Of course, it is especially widely used in chemistry. The main reason is that the chemical reaction of the system with a particular input is often not immediate but delayed, and the system has different time delays for a chemical reaction with different inputs. Therefore, it is necessary to analyze the dynamical behaviors of the system with multiple delays.In this article, we investigated the dynamical behaviors of the Lengyel–Epstein model with two delays which are different from the one studied in [5, 12, 14, 16]. By choosing the delay τii=1,2 as the bifurcation parameter, the stability of the positive equilibrium and the existence of Hopf bifurcation of the system are discussed in detail in three different cases. When τ1=0andτ2>0 and condition H2 holds, the positive equilibrium E∗ of system (6) is asymptotically stable for any τ2>0. For the same τi, when condition H3 holds, the positive equilibrium of system (6) is asymptotically stable for τ2∈0,τ20, and with the increase of τ2, the positive equilibrium loses its stability and a sequence of Hopf bifurcations occur. More generally, we further analyze the case, where τ2 is fixed to the stable interval and τ1 is taken as the bifurcation parameter. When τ1>0andτ2>0 and conditions H3−H5 hold, the positive equilibrium of system (6) is asymptotically stable for τ1∈0,τ10. Therefore, under certain conditions, if the value of the time-delay factor is lower than the critical value, the system is in an ideal stable state, which is convenient for the effective control of the chemical reaction in the system (See Figure 6). With the increase of τ1, the positive equilibrium loses its stability and a sequence of Hopf bifurcations occurs. In this case, the phenomenon of small amplitude periodic solution is generated near the equilibrium point of the system, which is not conducive to the effective control of the reaction in the system (see Figure 7). In addition, by applying the normal form and the center manifold theorem, the explicit formula determining the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions are obtained. Moreover, some numerical examples also support our theoretical results. The above results indicate that the stability of chemical reactions can be controlled by setting parameters, which plays a positive role in the safety and effectiveness of chemical experiments. Secondly, it is worth noting that there are many factors that affect the dynamic properties of the chemical reaction system. In this paper, we only consider the influence of delay parameters on the system, and there are many works worth study, such as considering the influence of the concentration of reactants in the chemical reaction process, diffusion, and other factors on the dynamic properties of the chemical reaction system and analyzing the global Hopf bifurcation of the model. Furthermore, due to the discrete nature of chemical reactions, the study of discrete systems is meaningful and interesting, and there are opportunities to further consider the dynamic properties of discrete systems with multiple delays in the future.



---

*Source: 5554562-2021-09-15.xml*