
Citation: Jing-Feng Tian, Ming-Hu Ha, Hong-Jie Xing. Properties of the power-mean and their applications[J]. AIMS Mathematics, 2020, 5(6): 7285-7300. doi: 10.3934/math.2020466
[1] | Serap Özcan, Saad Ihsan Butt, Sanja Tipurić-Spužević, Bandar Bin Mohsin . Construction of new fractional inequalities via generalized n-fractional polynomial s-type convexity. AIMS Mathematics, 2024, 9(9): 23924-23944. doi: 10.3934/math.20241163 |
[2] | Hasan Kara, Hüseyin Budak, Mehmet Eyüp Kiriş . On Fejer type inequalities for co-ordinated hyperbolic ρ-convex functions. AIMS Mathematics, 2020, 5(5): 4681-4701. doi: 10.3934/math.2020300 |
[3] | Wei Gao, Artion Kashuri, Saad Ihsan Butt, Muhammad Tariq, Adnan Aslam, Muhammad Nadeem . New inequalities via n-polynomial harmonically exponential type convex functions. AIMS Mathematics, 2020, 5(6): 6856-6873. doi: 10.3934/math.2020440 |
[4] | Naila Mehreen, Matloob Anwar . Ostrowski type inequalities via some exponentially convex functions with applications. AIMS Mathematics, 2020, 5(2): 1476-1483. doi: 10.3934/math.2020101 |
[5] | Tekin Toplu, Mahir Kadakal, İmdat İşcan . On n-Polynomial convexity and some related inequalities. AIMS Mathematics, 2020, 5(2): 1304-1318. doi: 10.3934/math.2020089 |
[6] | Ghulam Farid, Saira Bano Akbar, Shafiq Ur Rehman, Josip Pečarić . Boundedness of fractional integral operators containing Mittag-Leffler functions via (s,m)-convexity. AIMS Mathematics, 2020, 5(2): 966-978. doi: 10.3934/math.2020067 |
[7] | M. Emin Özdemir, Saad I. Butt, Bahtiyar Bayraktar, Jamshed Nasir . Several integral inequalities for (α, s,m)-convex functions. AIMS Mathematics, 2020, 5(4): 3906-3921. doi: 10.3934/math.2020253 |
[8] | Anjum Mustafa Khan Abbasi, Matloob Anwar . Ostrowski type inequalities for exponentially s-convex functions on time scale. AIMS Mathematics, 2022, 7(3): 4700-4710. doi: 10.3934/math.2022261 |
[9] | Saad Ihsan Butt, Ahmet Ocak Akdemir, Muhammad Nadeem, Nabil Mlaiki, İşcan İmdat, Thabet Abdeljawad . (m,n)-Harmonically polynomial convex functions and some Hadamard type inequalities on the co-ordinates. AIMS Mathematics, 2021, 6(5): 4677-4690. doi: 10.3934/math.2021275 |
[10] | Andrea Aglić Aljinović, Domagoj Kovačević, Mehmet Kunt, Mate Puljiz . Correction: Quantum Montgomery identity and quantum estimates of Ostrowski type inequalities. AIMS Mathematics, 2021, 6(2): 1880-1888. doi: 10.3934/math.2021114 |
Plankton is divided into two groups: zooplankton and phytoplankton. Phytoplankton are the primary producers in aquatic ecological models, as well as the main supplier of dissolved oxygen to phytoplankton blooms. Phytoplankton opens up the food web of aquatic ecosystems. Zooplankton, as economic aquatic animals, constitute an important feed for fish and other economic animals in the middle and upper waters, which is of great significance to the development of fishery[1,2,3,4].
Plankton has been studied extensively in many ways. In References[1,2,3,4], the researchers found that zooplankton eat plankton and zooplankton smaller than themselves, or feed on algae, bacteria, copepods and other food scraps. Therefore, there is a predator-prey relationship between zooplankton and phytoplankton. In 1939, Fleming[5] presented the first mathematical model of plankton. Since then, researchers have done a lot of work on plankton[6,7,8,9,10], focusing on factors such as nutrients, temperature, light, viral diseases and harvest to understand the bloom and disappearance of algal blooms. However, the researchers found that toxins released by toxic-producing phytoplankton have an effect on the termination of plankton blooms, which means that toxic chemicals can act as biological control for other plankton populations[11].
When the predator captures the prey, the prey seeks shelter because of a survival instinct. Nature provides shelter to the prey, and this behavior keeps the balance of the predator-prey model[12,13,14]. For lake ecosystems, prey refuge can stabilize plankton biomass by preventing phytoplankton from being temporarily eaten by zooplankton. Scholars have found that phytoplankton shelters can be obtained through benthic sediments, which can allow phytoplankton to temporarily escape from zooplankton predation. At the same time, the water layer can also form a temporary shelter for phytoplankton, and the shelter can prevent the extinction of the prey population[15,16,17]. Li et al.[18] proposed a model with the refuge as follows:
{dPdt=rP(1−PK)−β1(P−m)Za1+(P−m),dZdt=β2(P−m)Za1+(P−m)−dZ−θPZa2+P, | (1.1) |
where P and Z represent the number of phytoplankton and zooplankton, respectively. r represents the intrinsic growth rate of phytoplankton, K represents the environmental carrying capacity of phytoplankton, β1 represents the predation rate, β2 represents the conversion rate, d represents the mortality rate of zooplankton, θ represents the toxin release rate, a1 and a2 represent half-full constants, m represents the number of protected phytoplankton when phytoplankton have the ability to shelter and P−m denotes the number of unprotected phytoplankton that can be preyed upon by zooplankton. They[18] studied the effect of refuge of phytoplankton on the phytoplankton - zooplankton model.
In nature, the population diffuse from one area to another in order to survive. The predator-prey model with diffusion can generate complex spatial patterns [19,20,21,22]. There are two types of diffusion: self-diffusion and cross - diffusion. The former refers to the diffusion of one species from the higher - density area to the lower - density area in order to survive, while the latter refers to the diffusion of one species as influenced by other species. Meanwhile, time delay in nature is an important factor affecting the predator - prey model. A dynamical model without time delay can only be an approximation [23]. Generally speaking, there are many kinds of delayed factors in the growth process, such as the digestion of food [24], the maturation of cells [25], pregnancy[26,27]. These processes are not instantaneous and take time to complete. Later, it is found that the existence of time delay would make the positive steady state of a predator - prey model lose stability, resulting in bifurcation or periodic oscillation[21,28]. Zhao et al.[19] proposed a reaction - diffusion model with mature delay:
{∂P∂t=d1ΔP+rP(1−PK)−μPZα+P,∂Z∂t=d2ΔZ+μ1PZα+P−δZ−ρP(t−τ)Zα+P(t−τ), | (1.2) |
where d1 and d2 represent the diffusion coefficients. τ is the time required for phytoplankton to form a mature cell to release toxins. They[19] analyzed the stability of the equilibrium, the existence and properties of Hopf bifurcation. They[19] also found that time delay has an effect on the model (1.2). Recently, Hopf bifurcation has also continued to be investigated in fractional-order dynamical systems [29,30,31,32] and integer-order differential systems [33,34,35,36,37].
In lake ecosystems, fish eat plankton for survival, while zooplankton eat phytoplankton. Thus, phytoplankton, zooplankton and fish form a food chain. There is a lot of work that has done by many researchers[20,38,39].
The functional response reflects the predator's predation on the prey, which is an important part of the predator-prey model. It can be divided into prey-dependent (Holling I-IV type[40], Ivlev type [41], Rosenzweig type [22]) and predator-dependent functional responses (Beddington-DeAngelis type [42,43], Crowley-Martin (C-M) type [44], Hassell-Varley type [45], ratio-dependent type[46]). A predator-prey model that takes into account interactions between predators is more realistic. In [44], Crowley and Martin first proposed the C-M functional response:
F(H,P)=aH(1+abH)(1+cP), | (1.3) |
where F(H,P) represents the predation rate per predator, H represents the density of the prey per unit of area, P represents the density of the predator per unit of area, a represents the attack coefficient, b is the handling time and c is the interference coefficient. This functional response indicates that there is interference between predators as they feed on and handle prey. In [47], the authors proposed a model with the C-M functional response and showed that the system has complex dynamical behavior. In [48], the interaction between mature prey and predator is assumed to be the C-M functional response; the authors analyzed the positivity, boundedness and existence of equilibrium points. They not only analyzed the stability behavior of the delayed and non-delayed system, but also discussed the properties of Hopf bifurcation by choosing delay as the bifurcation parameter. In recent years, the researchers have done a great deal of work on this predator-prey models with a C-M functional response[49,50,51]. Therefore, in this paper, on the basis of Eq (1.3), we will consider that fish predation on zooplankton follows the C-M functional response:
p(Z,F)=γZF(1+aZ)(1+bF), | (1.4) |
where Z and F represent the number of zooplankton and fish, respectively. a represents interference between zooplankton, b is the interference between fish and γ is the maximum rate of fish predation on zooplankton. For different values of a and b, we have that (i) Eq (1.4) has turned into a Holling II functional response if a>0 and b=0; (ii) Eq (1.4) expresses a saturation response with respect to a predator if a=0 and b>0; (iii) Eq (1.4) becomes a linear mass-action functional response if a=0 and b=0.
In this paper, we consider the self-diffusion of populations. Let P(x,t), Z(x,t) and F(x,t) be the densities of phytoplankton, zooplankton and fish at the location x and the time t, respectively. We also give the following assumptions.
(1) The phytoplankton population follows the logistic growth under the condition of no zooplankton.
(2) Zooplankton preyed upon by fish is a C-M functional response.
(3) Phytoplankton preyed upon by zooplankton is a Holling II functional response. Nature provides shelter to the prey, so a constant proportion m∈(0,1) of the phytoplankton take refuge, leaving (1−m)P of the unprotected phytoplankton available for zooplankton grazing, following the Holling II functional response, i.e., f1(P)=(1−m)Pa1+(1−m)P.
(4) The progress of release of toxins takes the Holling II functional response and considers the mature delay of toxins in the cell, i.e., f2(P)=P(t−τ)a2+P(t−τ).
Based on the above assumptions, a schematic diagram that expresses the interactions of phytoplankton, zooplankton and fish is depicted in Figure 1. The corresponding model is
{∂P(x,t)∂t=d1ΔP+r1P(1−PK)−β1(1−m)PZa1+(1−m)P,x∈(0,lπ),t>0,∂Z(x.t)∂t=d2ΔZ+β2(1−m)PZa1+(1−m)P−δP(t−τ)Za2+P(t−τ)−γ1ZF(1+aZ)(1+bF)−g1Z,x∈(0,lπ),t>0,∂F(x,t)∂t=d3ΔF+γ2ZF(1+aZ)(1+bF)−g2F,x∈(0,lπ),t>0,Px(0,t)=Zx(0,t)=Fx(0,t)=0,t>0,Px(lπ,t)=Zx(lπ,t)=Fx(lπ,t)=0,t>0,P(x,t)>0,Z(x,t)>0,F(x,t)>0,x∈[0,lπ],t∈[−τ,0], | (1.5) |
where d1, d2 and d3 represent the diffusion coefficients of each population, respectively; Δ is the Laplace operator, β1 represents the maximum predation rate of phytoplankton by zooplankton, β2 is the conversion rate, K represents the environmental carrying capacity, γ1 represents the predation rate, γ2 represents the conversion rate, m is the refuge proportion of phytoplankton, α1 and α2 are half saturation constants, δ represents the release rate of phytoplankton toxins and g1 and g2 represent the natural mortality rates of zooplankton and fish, respectively. a represents the degree of interference between zooplankton, b represents the degree of interference between fish and τ represents the mature time delay of phytoplankton toxin release. The Neumann boundary condition indicates that the area is closed and no individuals can move across this area.
We rescale the model (1.5) by
˜P=PK,˜Z=Z,˜F=F,˜t=r1t,~d1=d1r1,~β1=β1Kr1,~α1=a1K,˜a=1a,˜b=1b,~d2=d2r1,~β2=β2r1,~α2=a2K,˜δ=δr1,~γ1=γ1abr1,~d3=d3r1,~γ2=γ2abr1,~g1=g1r1,~g2=g2r1. |
For the sake of convenience, omitting the breaking line, the model (1.5) becomes
{∂P∂t=d1ΔP+P(1−P)−β1(1−m)PZα1+(1−m)P,x∈(0,lπ),t>0,∂Z∂t=d2ΔZ+β2(1−m)PZα1+(1−m)P−δP(t−τ)Zα2+P(t−τ)−γ1ZF(a+Z)(b+F)−g1Z,x∈(0,lπ),t>0,∂F∂t=d3ΔF+γ2ZF(a+Z)(b+F)−g2F,x∈(0,lπ),t>0,Px(0,t)=Zx(0,t)=Fx(0,t)=0,t>0,Px(lπ,t)=Zx(lπ,t)=Fx(lπ,t)=0,t>0,P(x,t)>0,Z(x,t)>0,F(x,t)>0,x∈[0,lπ],t∈[−τ,0], | (1.6) |
where 0<P<1, 0<m<1 and all parameters are positive.
Our paper is organized as follows. The existence and stability of equilibrium of the model (2.1) are discussed in Section 2. Meanwhile, the occurrence of Hopf bifurcation is given by choosing m as a bifurcation parameter. The existence and properties of Hopf bifurcation of the model (3.1) are discussed in Section 3. In Section 4, Hopf bifurcation of the reaction-diffusion model (1.6) at the positive equilibrium and its properties are analyzed. In Section 5, a numerical simulation is demonstrated to prove the previous theoretical results by using Matlab software. In the last section, we have a brief discussion.
In this section, we will investigate the dynamical behavior of the model (1.6) with no delay and no diffusion. That is, the model is
{dPdt=P(1−P)−β1(1−m)PZα1+(1−m)P,dZdt=β2(1−m)PZα1+(1−m)P−δPZα2+P−γ1ZF(a+Z)(b+F)−g1Z,dFdt=γ2ZF(a+Z)(b+F)−g2F. | (2.1) |
According to the existence theorem of the solution of ordinary differential equations, we can know that the solution of the model (2.1) exists. By Lemma 2.1 in Reference[52], we give the following lemma to explain the positivity of the solution of the model (2.1).
Lemma 2.1. All solutions of the model (2.1) that start positive remain positive.
Proof. Let (P(t),Z(t),F(t)) be any solution of the model (2.1). Assuming that the initial time is t0 and one solution of the model (2.1) is at least not positive, then we have the following three cases:
(i) there exists time t1 such that P(t0)>0,P(t1)=0,P′(t1)≤0,Z(t)>0,F(t)>0,t0<t<t1;
(ii) there exists time t2 such that Z(t0)>0,Z(t2)=0,Z′(t2)≤0,P(t)>0,F(t)>0,t0<t<t2;
(iii) there exists time t3 such that F(t0)>0,F(t3)=0,F′(t3)≤0,P(t)>0,Z(t)>0,t0<t<t3.
If the first case is true, then we get P′(t1)=0. This contradicts with P′(t1)≤0. Similarly, we have that Z′(t2)=0, which contradicts with Z′(t2)≤0. And, F′(t3)=0, which contradicts with F′(t3)≤0.
Because of the arbitrariness of P(t),Z(t) and F(t), all solutions of the model (2.1) remain positive for all t>t0. Thus, all solutions of the model (2.1) that start positive remain positive.
Obviously, the model (2.1) has three boundary equilibria: one trivial equilibrium E0(0,0,0) and two boundary equilibria E1(1,0,0) and E2(P2,Z2,0) if (H1):β2−δ−g1>0 holds, where
P2=−B+√B2+4(1−m)(β2−δ−g1)g1α1α22(1−m)(β2−δ−g1),Z2=(1−P2)[α1+(1−m)P2](1−m)β1; |
here, B=(β2−g1)(1−m)α2−(δ+g1)α.
Assume that the model (2.1) has the coexistence equilibrium E∗(P∗,Z∗,F∗), where P∗, Z∗ and F∗ satisfy
{P∗(1−P∗)−β1(1−m)P∗Z∗α1+(1−m)P∗=0,β2(1−m)P∗Z∗α1+(1−m)P∗−δP∗Z∗α2+P∗−γ1Z∗F∗(a+Z∗)(b+F∗)−g1Z∗=0,γ2Z∗F∗(a+Z∗)(b+F∗)−g2F∗=0. | (2.2) |
By the first equation of (2.2), we can obtain
Z∗=(1−P∗)[α1+(1−m)P∗](1−m)β1. | (2.3) |
Substituting Eq (2.3) into the third equation of (2.2), we have
F∗=(γ2−bg2)(1−P∗)[α1+(1−m)P∗]−bg2aβ1(1−m)g2[aβ1(1−m)+(1−P∗)(α1+(1−m)P∗)]. | (2.4) |
If (H2):(γ2−bg2)(1−P∗)[α1+(1−m)P∗]−bg2aβ1(1−m)>0 holds, then we have that F∗>0.
Substituting Eqs (2.3) and (2.4) into the second equation of (2.2), P∗ is the positive root of the equation
f(P)=M5P5+M4P4+M3P3+M2P2+M1P+M0=0, | (2.5) |
where
M5=(1−m)2γ2(β2−δ−g1),M4=−γ2(1−m)(1−m−α1)(β2−δ−g1)+δγ2(1−m)(1−m−α1)−γ2(1−m)2(1−α2)(β2−g1)−g1γ2α1(1−m),M3=−[aβ1(1−m)+α1]γ2(1−m)(β2−δ−g1)+(β2−g1γ2(1−m)(1−α2)(1−m−α1)−δγ2(1−m−α1)2+g1γ2α1(1−m−α1)+γ1(1−m)2(γ2−bg2)−α2β2γ2(1−m)2+γ2(1−m)α1(δ+g1)+g1γ2α2(1−m)(1−m−α1),M2=α1(g1γ2α2+γ1β1)(1−m)+α2β2γ2(1−m)(1−m−α2)+(1−m)(1−α2)α1γ2(β2−g1)−γ2α1(1−m−α1)(2δ+g1)+aβ1g1γ2α1(1−m)+g1γ2α2−aβ1δγ2(1−m)(1−m−α1)+β1(1−m)2(1−α2)[aβ2γ2−ag1γ2−γ1(γ2−bg2)],M1=[aβ1(1−m)+α1][γ2α2(β2−g1)+α1g1γ2α2−γ2α1(δ+g1)]−(1−m−α1)g1γ2α1α2+γ1abg2β21(1−m)2−γ1β1(1−m)(γ2−bg2)[α2(1−m−α1)+α1],M0=γ1abg2β21α2(1−m)2−[aβ1(1−m)+α1]g1γ2α1α2−γ1β1(1−m)(γ2−bg2)α1α2. |
Under the condition (H1), we have that M5>0. By the Descarte's rule of signs[53], Eq (2.5) has only one positive root P∗ if and only if one of the following terms is satisfied:
(1)M4>0,M3>0,M2>0,M1>0,M0<0;(2)M4>0,M3>0,M2>0,M1<0,M0<0;(3)M4>0,M3>0,M2<0,M1<0,M0<0;(4)M4>0,M3<0,M2<0,M1<0,M0<0;(5)M4<0,M3<0,M2<0,M1<0,M0<0. |
Here, we give two figures as examples of the first two cases to verify the conclusion that there is only one positive root (see Figure 2). We first give the assumption (H3): one of the conditions (1)–(5) is true. Therefore, the following conclusion can be obtained.
Theorem 2.1. Under the conditions (H1), (H2) and (H3), the model (2.1) has only one positive equilibrium E∗(P∗,Z∗,F∗), which is determined by Eqs (2.3)–(2.5).
The stability of all equilibria will be analyzed in this part.
The Jacobian matrix of the model (2.1) is
A=(a11a120a21a22a230a32a33), | (2.6) |
where
a11=1−2P−α1β1(1−m)Z[α1+(1−m)P]2,a12=−β1(1−m)Pα1+(1−m)P,a21=α1β2(1−m)Z[α1+(1−m)P]2−δα2Z(α2+P)2,a23=−γ1bZ(a+Z)(b+F)2,a22=β2(1−m)Pα1+(1−m)P−δPα2+P−γ1aF(a+Z)2(b+F)−g1,a32=γ2aF(a+Z)2(b+F),a33=γ2bZ(a+Z)(b+F)2−g2. |
The characteristic equation of the model (2.1) is
λ3−(a11+a22+a33)λ2+(a11a22+a11a33+a22a33)λ+a11a32a23−a11a22a33+a12a21a33=0. | (2.7) |
According to Eq (2.7), we have the following conclusions.
Theorem 2.2. The boundary equilibrium E0(0,0,0) of the model (2.1) is always unstable.
Proof. The characteristic equation of the model (2.1) at E0 is
(λ−1)(λ+g1)(λ+g2)=0. |
It has three roots:
λ1=1>0,λ2=−g1<0,λ3=−g2<0. |
Thus, the boundary equilibrium E0 is unstable.
Theorem 2.3. If (H4):β2(1−m)α1+1−m−δα2+1−g1<0 holds, then the boundary equilibrium E1(1,0,0) of the model (2.1) is locally asymptotically stable.
Proof. The characteristic equation of the model (2.1) at E1 is
(λ+1)[λ−(β2(1−m)α1+(1−m)−δα2+1−g1)](λ+g2)=0. |
It has three roots:
λ1=−1<0,λ2=β2(1−m)α1+(1−m)−δα2+1−g1,λ3=−g2<0. |
When (H4) is true, then λ2<0. Therefore, the boundary equilibrium E1 is locally asymptotically stable under the condition (H4).
The characteristic equation of the model (2.1) at E2(P2,Z2,0) is
(λ−s1)(λ2+s2λ+s3)=0, |
where
s1=γ2(1−P2)[α1+(1−m)P2]abβ1(1−m)+b(1−P2)[α1+(1−m)P2]−g2,s2=α1−β1(1−m)Z2+(1−m)P2−2β2(1−m)P2α1+(1−m)P2,s3=−1(α2+P2)[α1+(1−m)P2]2[(1−2P2)(α1+(1−m)P2)−α1(1−P2)][β2(1−m)P2(α2+P2)−(α1+(1−m)P2)((δ+g1)P2+α2g1)]+α1P2(1−(1−m)P2)[α1+(1−m)P2]2−δα2P2(1−P2)(α+P2)2. |
Assume that the condition (H5):s1<0, s22−4s3>0, s2<0 and s3>0 are true; we can get the following result.
Theorem 2.4. Under the assumptions (H1) and (H5), the boundary equilibrium E2(P2,Z2,0) of the model (2.1) is locally asymptotically stable.
The characteristic equation of the model (2.1) at E∗(P∗,Z∗,F∗) is
λ3+M6λ2+M7λ+M8=0, | (2.8) |
where
M6=−(A11+A22+A33),M7=A22A33+A11A33+A11A22−A12A21−A23A32,M8=A11A23A32+A12A21A33−A11A22A33,A11=−P∗+β1(1−m)2Z∗[α1+(1−m)P∗]2,A12=−β1(1−m)P∗α1+(1−m)P∗,A21=α1β2(1−m)Z∗[α1+(1−m)P∗]2−δα2Z∗(α2+P∗)2,A22=γ1F∗Z∗(a+Z∗)2(b+F∗),A23=−γ1bZ∗(a+Z∗)(b+F∗)2,A32=γ2aF∗(a+Z∗)2(b+F∗),A33=−γ2F∗Z∗(a+Z∗)(b+F∗)2. |
From the Routh-Hurwitz criterion[54], if M6>0, M7>0 and M6M7−M8>0, then all solutions of Eq (2.8) have negative real parts. When M6>0, M7>0 and M6M7−M8<0, Eq (2.8) has one negative root and a pair of complex roots with a positive real part. Assume that (H6): M6>0, M7>0 and M6M7−M8>0; then, the stability of the positive equilibrium E∗ will be obtained.
Theorem 2.5. Suppose that the conditions (H1)–(H3) are true. If the condition (H6) holds, then E∗ is locally asymptotically stable. Further, E∗ loses stability when M6M7−M8 passes through 0; in other words, the model (2.1) undergoes a Hopf bifurcation at E∗ when M6M7−M8=0.
Next, we will choose m as the bifurcation parameter to study the occurrence of Hopf bifurcation of the model (2.1) at E∗. By using the results in Reference[55], the following result can be obtained.
Theorem 2.6. If the characteristic equation of the model (2.1) at E∗(P∗,Z∗,F∗) is
λ3+M6(m)λ2+M7(m)λ+M8(m)=0, |
where M6(m), T(m)=M6(m)M7(m)−M8(m) and M8(m) are the smooth functions of m and there exists a positive number m∗ that satisfies
(1) M6(m∗)>0, T(m∗)=0 and M8(m∗)>0;
(2) dTdm|m=m∗≠0,
then Hopf bifurcation occurs at E∗(P∗,Z∗,F∗) when m=m∗.
We used Matlab software for numerical simulations to obtain this result in Section 5. Meanwhile, we can also choose δ as a bifurcation parameter to study the occurrence of Hopf bifurcation of the model (2.1) at E∗(P∗,Z∗,F∗). Since the discussion process is similar, we will only give the bifurcation diagram in Section 5.
Under some conditions, Hopf bifurcation may not take place. Thus, we will discuss the global asymptotical stability of the positive equilibrium E∗ as follows.
Theorem 2.7. Suppose that the conditions (H1)–(H3) are true. The positive equilibrium E∗ of the model (2.1) is globally asymptotically stable if 1−β1(1−m)2Z∗α1[α1+(1−m)P∗]−δ2[α2+P∗]2>0 and [α1+(1−m)P∗]β1P∗δβ2α1(α2+P∗)−2−β21(1−m)P∗γ2α1g1ab>0.
Proof. Let (P,Z,F) be any positive solution of the model (2.1). Define a Lyapunov function
V(t)=P−P∗−P∗lnPP∗+[α1+(1−m)P∗]β1β2α1(Z−Z∗−Z∗lnZZ∗)+F−F∗−F∗lnFF∗. |
Calculating the derivative of V(t) along the solution of the model (2.1), then we have
dV(t)dt=(P−P∗){−(P−P∗)+β1(1−m)(1−m)Z∗(P−P∗)−[α1+(1−m)P∗](Z−Z∗)[α1+(1−m)P∗][α1+(1−m)P]}+[α1+(1−m)P∗]β1β2α1(Z−Z∗){β2(1−m)α1[α1+(1−m)P∗][α1+(1−m)P](P−P∗)+β2(1−m)P∗Z[α1+(1−m)P∗](Z−Z∗)−δα2(α2+P)(α2+P∗)(P−P∗)−P∗δα2+P∗(Z−Z∗)−γ1b(a+Z)(b+F)(b+F∗)(F−F∗)−γ1aF∗Z(a+Z)(a+Z∗)(b+F∗)(Z−Z∗)−g1Z(Z−Z∗)}+(F−F∗){γ2a(a+Z)(b+F)(a+Z∗)(Z−Z∗)−γ2Z∗(b+F)(a+Z∗)(b+F∗)(F−F∗)}≤−{1−β1(1−m)2Z∗[α1+(1−m)P][α1+(1−m)P∗]−δ2α22[α2+P]2[α2+P∗]2}(P−P∗)2−{−2−β1(1−m)P∗α1Z+[α1+(1−m)P∗]β1P∗δβ2α1(α2+P∗)+γ1aF∗Z(a+Z)(a+Z∗)(b+F∗)+g1Z}(Z−Z∗)2−{γ2Z∗(b+F)(a+Z∗)(b+F∗)+γ22a2(a+Z)2(b+F)2(a+Z∗)2}(F−F∗)2≤−{1−β1(1−m)2Z∗α1[α1+(1−m)P∗]−δ2[α2+P∗]2}(P−P∗)2−{−2−β21(1−m)P∗γ2α1g1ab+[α1+(1−m)P∗]β1P∗δβ2α1(α2+P∗)}(Z−Z∗)2−{γ2Z∗(b+F)(a+Z∗)(b+F∗)+γ22a2(a+Z)2(b+F)2(a+Z∗)2}(F−F∗)2. |
Here, it is obvious that γ2Z∗(b+F)(a+Z∗)(b+F∗)+γ22a2(a+Z)2(b+F)2(a+Z∗)2>0. If 1−β1(1−m)2Z∗α1[α1+(1−m)P∗]−δ2[α2+P∗]2>0 and −2−β21(1−m)P∗γ2α1g1ab+[α1+(1−m)P∗]β1P∗δβ2α1(α2+P∗)>0 hold, then we have that the coefficients of (P−P∗)2, (Z−Z∗)2 and (F−F∗)2 are always negative. Thus, dV(t)dt is negative. This completes the proof.
Here, under the conditions (H1)–(H3), we choose delay τ as the bifurcation parameter and study its influence on the stability of the positive equilibrium E∗(P∗,Z∗,F∗). The model (2.1) is
{dPdt=P(1−P)−β1(1−m)PZα1+(1−m)P,dZdt=β2(1−m)PZα1+(1−m)P−δP(t−τ)Zα2+P(t−τ)−γ1ZF(a+Z)(b+F)−g1Z,dFdt=γ2ZF(a+Z)(b+F)−g2F. | (3.1) |
Let u1(t)=P(t)−P∗, u2(t)=Z(t)−Z∗, u3(t)=F(t)−F∗ and u(t)=(u1(t),u2(t),u3(t))T ∈R3. The linearized system of the model (3.1) at E∗ is
{du1(t)dt=b11u1(t)+b12u2(t),du2(t)dt=b21u1(t)+b22u2(t)+b23u3(t)+cu1(t−τ),du3(t)dt=b32u2(t)+b33u3(t), | (3.2) |
where
b11=A11,b12=A12,b21=α1β2Z∗(1−m)[α1+P∗(1−m)]2,b22=A22,b23=A23,b32=A32,b33=A33,c=−δα2Z∗(α2+P∗)2. |
Then, the model (3.2) can also be given by
du(t)dt=L1u(t)+L2u(t−τ), | (3.3) |
where
L1=(b11b120b21b22b230b32b33),L2=(000c00000). | (3.4) |
Thus, O(0, 0, 0) is the zero equilibrium of the model (3.2). The characteristic equation of the model (3.2) at O is
λ3+D2λ2+D1λ+D0+e−λτ(D3λ+D4)=0, | (3.5) |
where
D2=−(b11+b22+b33),D1=b11b22+b11b33+b22b33−b12b21−b23b32,D0=b11b23b32+b12b21b33−b11b22b33,D3=−cb12,D4=cb12b33. |
The roots of Eq (3.5) have been discussed above when τ=0. Next, we will study the effect of delay τ(τ>0) on the model (3.2).
Suppose that λ=iω(ω>0) is a pair of pure imaginary roots of Eq (3.5). Substituting it into Eq (3.5), we can obtain
−iω3−D2ω3+D1iω+D0+(cosωτ−isinωτ)(iD3ω+D4)=0. | (3.6) |
Separating the real and imaginary parts of Eq (3.6), we can obtain
{ω3−D1ω=D3ωcosωτ−D4sinωτ,D2ω2−D0=D3ωsinωτ+D4cosωτ. | (3.7) |
From Eq (3.7), we can get
ω6+ω4(D22−2D1)+ω2(D21−2D0D2−D33)+D20−D24=0. | (3.8) |
Let z=ω2, D5=D22−2D1, D6=D21−2D0D1−D23 and D7=D20−D24. Then, Eq (3.8) can be rewritten as
f(z)=z3+D5z2+D6z+D7=0, | (3.9) |
and we have
f′(z)=3z2+2D5z+D6=0. | (3.10) |
If Eq (3.9) has at least one positive root, then Hopf bifurcation takes place. Assume that (H7):Δ1=D25−3D6≤0 and (H8):Δ1=D25−3D6>0. By Lemmas 2.2 and 4.2 in Reference [56], we can get the following results.
Since limz→+∞f(z)=+∞, Eq (3.9) has at least one positive root when D7<0.
If (H7) holds, then f(z) is monotonically increasing for z∈[0,+∞); so, when D7≥0 and (H7) hold, Eq (3.9) has no positive root for z∈[0,+∞).
When D7≥0 and (H8) hold, Eq (3.10) has two roots, that is, z∗1 and z∗2, where
z∗1=−D5+√Δ13,z∗2=−D5−√Δ13. |
Furthermore, we have
f″(z∗1)=−2D5+2√Δ1+2D5=2√Δ1>0,f″(z∗2)=−2D5−2√Δ1+2D5=−2√Δ1<0. |
Therefore, we can obtain z∗1 and z∗2 as the local minimum and the local maximum of f(z), respectively.
Theorem 3.1. For Eq (3.9), the following conclusions are true.
(1) If D7<0, then Eq (3.9) has at least one positive root.
(2) If the condition (H7) holds and D7≥0, then Eq (3.9) has no positive root.
(3) If the condition (H8) holds and D7≥0, then Eq (3.9) has two positive roots when z∗1>0 and f(z∗1)≤0.
Without loss of generality, suppose that Eq (3.9) has three positive roots defined by z1, z2 and z3, respectively. Thus, Eq (3.8) has three positive roots ω1=√z1, ω2=√z2 and ω3=√z3. From Eq (3.7), we can get
τjk=1ωkarccos[D3ω4k+(D2D4−D1D3)ω2k−D0D4D23ω2k+D24]+2jπωk,k=1,2,3,j=0,1,2,⋯; | (3.11) |
thus, ±iωk is a pair of purely imaginary roots of Eq (3.5) when τ=τjk.
Define
τ0=mink=1,2,3.j∈N0{τjk},ω0=ωk|τ=τ0. |
Let λ(τ)=α(τ)+iβ(τ) be the root of Eq (3.5) near τ=τjk, and let it satisfy α(τjk)=0 and β(τjk)=ωk, j=0,1,2,⋯, k=1,2,3.
Theorem 3.2. Assume that zk=ω2k and f′(zk)≠0; then, dReλ(τ)dτ|λ=iω0≠0.
Proof. Differentiating Eq (3.5) for τ, we have
(dλ(τ)dτ)−1=(3λ2+2D2λ+D1)eλτ+D3λ(D3λ+D4)−τλ. | (3.12) |
Substituting λ(τ)=iω0 into Eq (3.12), we can get
(dλ(τ)dτ)−1|λ=iω0=(−3ω20+2D2iω0+D1)(cosω0τ+isinω0τ)+D3iω0(D3iω0+D4)−τiω0 |
=(D3ω20+iω0D4)(3ω20−D1−2D2ω0i)(cosω0τ+isinω0τ)+D3(D3ω20+iω0D4)D23ω40+D24ω20+τiω0. |
Therefore, we can get
Re(dλ(τ)dτ)−1|λ=iω0=3ω60+(2D22−4D1)ω40+(D21−2D0D2−D23)ω20D23ω40+D24ω20=zkD23ω40+D24ω20[3z2k+(2D22−4D1)zk+(D21−2D0D2−D23)]=zkf′(zk)D23ω40+D24ω20. |
Thus, we have
sign{dλ(τ)dτ|λ=iω0}=sign{Re(dλ(τ)dτ)−1|λ=iω0}=signf′(zk)≠0. |
This ends the proof of Theorem 3.2.
Theorem 3.3. Under the conditions (H1)–(H3), the dynamics of the model (3.1) at the positive equilibrium E∗ can be obtained.
(i) If D7≥0 and the condition (H7) is true, then E∗ is locally asymptotically stable for all τ>0.
(ii) If D7<0 or the condition (H8), D7>0, z∗1>0 and f(z∗1)≤0 are true, then E∗ is locally asymptotically stable when τ∈[0,τ0); but, E∗ is unstable when τ>τ0.
(iii) If the conditions in (ii) are all satisfied and f′(zk)≠0 holds, then a Hopf bifurcation occurs at E∗ when τ=τ0.
For τ=τ0, the existence of Hopf bifurcation has been discussed. The properties of the Hopf bifurcation will be discussed in consideration of the work of Hassard et al. [57].
Let τ=τ0+μ, μ∈R; we can obtain that μ=0 is a Hopf bifurcation value of the model (3.1). Let the phase space be C=C([−1,0],R3) and t→tτ; then, the model (3.1) can be expressed as a functional differential equation (FDE) in C as
˙u(t)=Lμ(ut)+F(μ,ut), | (3.13) |
where Lμ:C→R3 and F:R×C→R3.
Define φ(θ)=(φ1(θ),φ2(θ),φ3(θ))T∈R3,θ∈[−1,0] such that
Lμ(φ)=(τ0+μ)L1φ(0)+(τ0+μ)L2φ(−1) | (3.14) |
and
F(μ,φ)=(τ0+μ)(F1(φ)F2(φ)F3(φ)), | (3.15) |
where L1 and L2 are defined in Eq (3.4), and
F1(φ)=C11φ21(0)+C12φ1(0)φ2(0),F2(φ)=C21φ21(0)+C22φ22(0)+C23φ23(0)+C24φ1(0)φ2(0)+C25φ2(0)φ3(0)+C26φ21(−1)+C27φ1(−1)φ2(0),F3(φ)=C31φ22(0)+C32φ23(0)+C33φ2(0)φ3(0), |
C11=−1+(1−m)2α1β1Z∗[α1+P∗(1−m)]3,C12=−α1β1(1−m)[α1+P∗(1−m)]2,C21=(1−m)2α1β2Z∗[α1+P∗(1−m)]3,C22=γ1aF∗(a+Z∗)3(b+F∗),C23=γ1bZ∗(a+Z∗)(b+F∗)3,C24=α1β2(1−m)[α1+P∗(1−m)]2,C25=−γ1ab(a+Z∗)2(b+F∗)2,C26=δα2Z∗(α2+P∗)2,C27=−δα2(α2+P∗)2,C31=−γ2aF∗(a+Z∗)3(b+F∗),C32=−γ2bZ∗(a+Z∗)(b+F∗)3,C33=γ2ab(a+Z∗)2(b+F∗)2. |
From the Riesz representation theorem[58], there is a matrix function whose elements are functions of bounded variation η(θ,μ)∈C([−1,0],R3); also,
Lμ(φ)=∫0−1dη(θ,μ)φ(θ),φ∈C. | (3.16) |
Let
η(θ,0)=τ0L1δ(θ)−τ0L2δ(θ+1), |
where δ(θ) is a Dirac function, as follows:
δ(θ)={0,θ≠0,1,θ=1. |
For φ∈C=C([−1,0],R3), we define that
A(μ)φ={dφ(θ)dθ,θ∈[−1,0),∫0−1dη(μ,s)φ(s),θ=0,R(μ)φ={0,θ∈[−1,0),F(μ,φ),θ=0. |
Therefore, Eq (3.13) becomes
˙u(t)=Aμ(ut)+Rμut. | (3.17) |
For ψ∈C∗=C([−1,0],R3), let
A∗(μ)ψ(s)={−dψ(s)ds,s∈[0,1),∫0−1ψ(−ξ)dη(ξ,0),s=0, |
and establish the bilinear inner product
⟨ψ(s),φ(θ)⟩=ˉψ(0)φ(0)−∫0−1∫θ0ˉψ(ξ−θ)dη(θ)φ(ξ)dξ,ψ(s)∈C∗,φ(θ)∈C, |
where η(θ)=η(θ,0). Then, A(0) and A∗(0) are adjoint operators. According to Theorem 3.3, we can get that ±iω0τ0 are eigenvalues of A(0). Thus, they are also eigenvalues of A∗(0). Assume that q(θ)=(1,q1,q2)Teiω0τ0θ is the eigenvector of A(0) corresponding to iω0τ0; then, A(0)q(θ)=iω0τ0q(θ). By the definition above, we have
dq(θ)dθ=iωτ0q(θ),θ∈[−1,0), |
L0q(0)=∫0−1dη(0)q(0)=A(0)=iωτ0q(0),θ=0, |
that is,
(τ0+μ)(b11−iω0b120b21+ce−iω0τ0b22−iω0b230b32b33−iω0)(1q1q2)=0. |
Solving it, we can get
(1q1q2)=(1iω0−b11b12b32(iω0−b11)b12(iω0−b33)). |
Similarly, we assume that q∗(s)=D(1,q∗1,q∗2)e−iω0τ0s is the eigenvector of A∗(0) corresponding to −iω0τ0; then, A∗(0)q∗(s)=−iω0τ0q∗(s). By the definition above, we can get
(1q∗1q∗2)=(1−iω0+b11b12+ceiω0τ0b23(iω0+b11)b12(iω0+b33)(b21+ceiω0τ0)). |
By Eq (3.2), ⟨q∗,q⟩ can be expressed as
⟨q∗,q⟩=¯q∗(0)q(0)−∫0−1∫θ0¯q∗(ξ−θ)dη(θ)q(ξ)dξ=ˉD[(1,q∗1,q∗2)(1,q1,q2)T−∫0−1∫θ0(1,¯q∗1,¯q∗2)e−iω0(ξ−θ)τ0dη(θ)(1,q1,q2)Teiξω0τ0dξ]=ˉD[1+q1¯q∗1+q2¯q∗2+τ0c¯q∗1eiω0τ0]. |
Thus, we can choose
ˉD=[1+q1¯q∗1+q2¯q∗2+τ0c¯q∗1eiω0τ0]−1; |
then, ⟨q∗,q⟩=1.
Next, we adopt the ideas of Hassard et al.[57] to compute the coordinates describing the center manifold C0 at μ=0. We assume that ut is the solution of Eq (3.13) when μ=0. Let
z(t)=⟨q∗,ut⟩,W(t,θ)=ut(θ)−2Re{z(t)q(θ)}. | (3.18) |
On the center manifold C0, we can get
W(t,θ)=W(z(t),ˉz(t),θ), |
where
W(z(t),ˉz(t),θ)=W21(θ)z22+W11zˉz+W02(θ)ˉz22+⋯; | (3.19) |
z and ˉz express the local coordinates for the center manifold C0 in the direction of q∗ and ˉq∗. We can get that W is real when ut is real. For the real solution ut∈C0 of Eq (3.17), when μ=0, we have
˙z(t)=⟨q∗,˙u(t)⟩=⟨q∗,A(0)ut+R(0)ut⟩=⟨A∗(0)q∗,ut⟩+¯q∗(0)F(0,ut)=iω0τ0z(t)+¯q∗(0)F(0,W(t,0)+2Re{z(t)q(0)})=iω0τ0z(t)+¯q∗(0)F(0,W(z,ˉz,0)+2Re{z(t)q(0)}def=iω0τ0z(t)+¯q∗(0)F0(z,ˉz). | (3.20) |
Eq (3.20) can also be written as
˙z(t)=iω0τ0z(t)+g(z,ˉz), |
where
g(z,ˉz)=¯q∗(0)F0(z,ˉz)=g20(θ)z22+g11(θ)zˉz+g02(θ)ˉz22+g21z2ˉz2+⋯. | (3.21) |
By Eq (3.18), we can get
˙W=˙ut−˙zq−˙ˉzˉq={AW−2Re{¯q∗(0)F0q(θ)},θ∈[−1,0),AW−2Re{¯q∗(0)F0q(0)}+F0,θ=0,def=AW+H(z,ˉz,θ), | (3.22) |
where
H(z,ˉz,θ)=H20(θ)z22+H11(θ)zˉz+H02(θ)ˉz22+H21z2ˉz2+⋯. | (3.23) |
From Eqs (3.22) and (3.23), we obtain
(A−2iω0τ0)W20(θ)=−H20(θ),AW11(θ)=−H11(θ),(A+2iω0τ0)W02(θ)=−H02(θ). | (3.24) |
By Eq (3.18), we have
ut(θ)=W(t,θ)+2Re{z(t)q(θ)}=W20(θ)z22+W11(θ)zˉz+W02(θ)ˉz22+(1,q1,q2)Teiω0τ0θz+(1,¯q1,¯q2)Te−iω0τ0θˉz+⋯. | (3.25) |
Then, we have
g(z,ˉz)=¯q∗(0)F(0,ut)=ˉD(1,¯q1∗,¯q2∗)τ0(F1(ut)F2(ut)F3(ut))=ˉDτ0[F1(ut)+¯q1∗F2(ut)+¯q2∗F3(ut)]=z22[2ˉDτ0(k11+k21¯q1∗+k31¯q∗2)]+zˉz[ˉDτ0(k12+k22¯q1∗+k32¯q∗2)]+ˉz22[2ˉDτ0(k13+k23¯q1∗+k33¯q∗2)]+z2ˉz2[2ˉDτ0(k14+k24¯q1∗+k34¯q∗2)], | (3.26) |
where
k11=C11+C12q1,k12=2C11+C12(ˉq1+q1),k21=C21+C22q21+C23q22+C24q1+C25q1q2+C26e−2iω0τ0+C27q1e−iω0τ0,k31=C31q21+C32q22+C33q1q2,k33=C31ˉq21+C32ˉq22+C33ˉq1ˉq2,k22=2C21+2C22q1ˉq1+2C23q2ˉq2+C24(q1+ˉq1)+C25(q1ˉq2+ˉq1q2)+2C26+C27(q1eiω0τ0+ˉq1e−iω0τ0),k32=2C31q1ˉq1+2C32q2ˉq2+C33(q1ˉq2+ˉq1q2),k13=C11+C12ˉq1,k23=C21+C22ˉq21+C23ˉq22+C24ˉq1+C25ˉq1ˉq2+C26e2iω0τ0+C27eiω0τ0, |
k14=C11(W(1)20(0)+2W(1)11(0))+C12(12W(1)20(0)ˉq1+q1W(1)11(0)+W(2)11(0)+12W(2)20(0)),k24=C21(W(1)20(0)+2W(1)11(0))+C22(W(2)20(0)+2W(2)11(0))+C23(W(3)20(0)+2W(3)11(0))+C24(12W(1)20(0)ˉq1+q1W(1)11(0)+W(2)11(0)+12W(2)20(0))+C25(12W(2)20(0)ˉq2+W(1)11(0)q2+q1W(3)11(0)+12W(3)20(0)ˉq1)+C26(eiω0τ0W(1)20(−1)+2e−iω0τ0W(1)11(−1))+C27(12ˉq1W(1)20(−1)+q1W(1)11(−1)+e−iω0τ0W(2)11(0)+12eiω0τ0W(2)20(0)),k34=C31(W(2)20(0)+2W(1)11(0))+C33(12ˉq2W(2)11(0)+q2W(2)11(0)+q1W(3)11(0)+12W(3)20(0)ˉq1). |
Comparing the coefficients of Eqs (3.21) and (3.26), we can get
g20=2ˉDτ0(k11+ˉq∗1k21+ˉq∗2k31),g11=ˉDτ0(k12+ˉq∗1k22+ˉq∗2k32),g02=2ˉDτ0(k13+ˉq∗1k23+ˉq∗2k33),g21=2ˉDτ0(k14+ˉq∗1k24+ˉq∗2k34). | (3.27) |
Since the expression of g21 contains W20(θ) and W11(θ), we must compute W20(θ) and W11(θ). According to Eq (3.22), when θ∈[−1,0), we get
H(z,ˉz,θ)=−2Re{ˉq∗(0)F0(z,ˉz)q(θ)}=−2Re{g(z,ˉz)q(θ)}=−g(z,ˉz)q(θ)−ˉg(z,ˉz)ˉq(θ). | (3.28) |
Comparing the coefficients of Eqs (3.23) and (3.28), we can receive
H20(θ)=−g20q(θ)−ˉg02ˉq(θ),H11(θ)=−g11q(θ)−ˉg11ˉq(θ). | (3.29) |
By Eq (3.24), we have
˙W20(θ)=2iω0τ0W20(θ)+g20q(θ)+ˉg02ˉq(θ),˙W11(θ)=g11q(θ)+ˉg11ˉq(θ). | (3.30) |
Solving Eq (3.30), we can obtain
W20(θ)=ig20q(0)ω0τ0eiω0τ0θ+iˉg02ˉq(0)3ω0τ0e−iω0τ0θ+E1e2iω0τ0θ,W11(θ)=−ig11q(0)ω0τ0eiω0τ0θ+iˉg11ˉq(0)ω0τ0e−iω0τ0θ+E2, | (3.31) |
where Ei=(E(1)i,E(2)i,E(3)i)∈R3 (i=1,2) is a constant vector.
For θ=0, from Eq (3.22), we have
H(z,ˉz,0)=−2Re{ˉq∗(0)F0(z,ˉz)q(0)}+F0. |
From Eq (3.23), we can get
H20(0)=−g20q(0)−ˉg02ˉq(0)+2τ0(k11k21k31),H11(0)=−g11q(0)−ˉg11ˉq(0)+τ0(k12k22k32). | (3.32) |
According to the meaning of A(0) and Eq (3.24), we have
∫0−1dη(θ)W20(θ)=2iω0τ0W20−H20(0),∫0−1dη(θ)W11(θ)=−H11(0), | (3.33) |
where η(θ)=η(θ,0).
From Eqs (3.14), (3.16) and (3.33), we can get
τ0L1W20(0)+τ0L2W20(−1)=2iω0τ0−H20(0),τ0L1W11(0)+τ0L2W11(−1)=H11(0). | (3.34) |
Substituting Eqs (3.31) and (3.32) into Eq (3.34), we have
E1=2(2iω0−b11−b120−ce−2iω0τ0−b212iω0−b22−b230−b322iω0−b33)−1(k11k21k31), |
E2=(−b11−b120−c−b21−b22−b230−b32−b33)−1(k11k21k31). |
Thus, all expressions of gij can be represented in full. Also, we have
{c1(0)=i2ω0τ0(g11g20−2|g11|2−|g02|23)+g212,μ1=−Re{c1(0)}Re{λ′(τ0)},μ2=2Re{c1(0)},T1=−Im{c1(0)}+μ1Im{λ′(τ0)}ω0τ0, | (3.35) |
which determine the direction of Hopf bifurcation and the stability of bifurcating periodic solutions on the center manifold at τ=τ0.
Theorem 3.4. From Eq (3.35), we have the following conclusions.
(i) μ1 determines the direction of Hopf bifurcation: if μ1>0, then the Hopf bifurcation is supercritical; if μ1<0, then the Hopf bifurcation is subcritical and the bifurcating periodic solutions exist when τ>τ0;
(ii) μ2 determines the stability of the bifurcating periodic solutions: if μ2<0, then the bifurcating periodic solutions are stable; if μ2>0, then the bifurcating periodic solutions are unstable;
(iii) T1 determines the period of the bifurcating periodic solutions: if T1>0, then the period increases; if T1<0, then the period decreases.
Next, we will analyze the existence and properties of Hopf bifurcation of the model (1.6). The model (1.6) is linearized at the positive equilibrium E∗ in the phase space C=C([−τ,0],R3):
du(t)dt=DΔu(t)+L1u(t)+L2u(t−τ), | (4.1) |
where D=diag{d1,d2,d3} and L1 and L2 are defined in Eq (3.4).
We know that Δ has the eigenvalues −(nl)2 and n∈N0 under Neumann boundary conditions in [0,lπ]. Then, the characteristic equation of the model (4.1) is
λ3+q2nλ2+q1nλ+q0n+e−λτ(q3nλ+q4n)=0, | (4.2) |
where
q2n=(d1+d2+d3)(nl)2−(b11+b22+b33),q1n=(d1d2+d2d3+d1d3)(nl)4−(d1b22+d2b11+d2b33+d3b22+d3b11+d1b33)(nl)2+(b11b22+b22b33+b11b33−b12b21−b23b32),q0n=d1d2d3(nl)6−(d1d2b33+d1d3b22+d2d3b11)(nl)4+(d1b22b33+d2b11b33+d3b11b22)(nl)2+b11b23b32+b12b21b33−b11b22b33,q3n=−cb12,q4n=−d3cb12(nl)2+cb12b33. |
When τ=0, Eq (4.2) can be rewritten as
λ3+q2nλ2+(q1n+q3n)λ+(q0n+q4n)=0. | (4.3) |
Furthermore, we have
q0n+q4n=−d3cb12(nl)2+cb12b33+d1d2d3(nl)6−(d1d2b33+d1d3b22+d2d3b11)(nl)4+(d1b22b33+d2b11b33+d3b11b22)(nl)2+b11b23b32+b12b21b33−b11b22b33, |
q2n(q1n+q3n)−(q0n+q4n)=[(d1+d2+d3)(d1d2+d2d3+d1d3)−d1d2d3](nl)6−[(d1+d2+d3)(d1b22+d2b11+d2b33+d3b22+d3b11+d1b33)+(b11+b22+b33)(d1d2+d2d3+d1d3)−(d1d2b33+d1d3b22+d2d3b11)](nl)4+[(d1+d2+d3)(b11b22+b22b33+b11b33−b12b21−b23b32−cb12)+(b11+b22+b33)(d1b22+d2b11+d2b33+d3b22+d3b11+d1b33)−(d1b22b33+d2b11b33+d3b11b22)+d3cb12](nl)2+(b11+b22+b33)(cb12−b11b22−b22b33−b11b33+b12b21+b23b32)−(b11b23b32+b12b21b33−b11b22b33)−cb12b33. |
Assume that the following condition holds true: (H9): q2n>0, q2n(q1n+q3n)−(q0n+q4n)>0, q0n+q4n>0, n∈N0. Using the Routh-Hurwitz criterion[54], we can obtain the next conclusion.
Theorem 4.1. If the conditions (H1)–(H3) and (H9) hold, then all roots of Eq (4.3) have negative real parts, that is, the positive equilibrium E∗ of the model (1.6) is locally asymptotically stable when τ=0.
When τ≠0, the time delay may have some effect on the model (1.6). Therefore, we will analyze the effect of delay τ on the positive equilibrium E∗. Let λ=iω2(ω2>0) be the solution of Eq (4.2); we can get
−iω32−ω22q2n+iω2q1n+q0n+e−iω2τ(iω2q3n+q4n)=0. | (4.4) |
Separating the real and imaginary parts of Eq (4.4), we can receive
{ω32−q1nω2=q3nω2cosω2τ−q4nsinω2τ,q2nω22−q0n=q3nω3nsinω2τ+q4ncosω2τ, | (4.5) |
which follows that
ω62+(q22n−2q1n)ω42+(q21n−2q0nq2n−q23n)ω22+(q20n−q24n)=0. | (4.6) |
Let p2n=q22n−2q1n, p1n=q21n−2q0nq2n−q23n, p0n=q20n−q24n and ω22=y; we can obtain
y3+p2ny2+p1ny+p0n=0. | (4.7) |
Further, if f(y)=y3+p2ny2+p1ny+p0n, then f′(y)=3y2+2p2ny+p1n. Assume that (H10):Δ2=p22n−3p1n≤0 and (H11):Δ2=p22n−3p1n>0 are true, here y∗1=−p1n+√Δ23 and y∗2=−p1n−√Δ23 are the local minimum and the local maximum of Eq (4.7), respectively. Similarly, we have the following conclusion by Lemmas 2.2 and 4.2 in[56].
Theorem 4.2. For Eq (4.7), the following results are true.
(1) If p0n<0, then Eq (4.7) has at least one positive root.
(2) If p0n≥0 and the condition (H10) holds, then Eq (4.7) has no positive root.
(3) If p0n≥0 and the condition (H11) holds, then Eq (4.7) has positive roots when y∗1>0 and f(y∗1)≤0.
Without loss of generality, we suppose that it has three positive roots defined by y1n, y2 and y3n. So, Eq (4.6) has three positive roots:
ω12=y1n,ω22=y2n,ω32=y3n. |
Substituting ωk2(k=1,2,3) into Eq (4.5), we can get
τjkn=1ωk2arccos[q3n(ωk2)4+(q2nq4n−q1nq3n)(ωk2)2−q0nq4nq23n(ωk2)2+q24n]+2jπωk2,k=1,2,3,j=0,1,2,⋯. | (4.8) |
Thus, when τ=τjkn, we get that λ=±iωk2 is a pair of purely imaginary roots of Eq (4.2). Define
τn0=mink=1,2,3.j∈N0{τjkn},ωn0=ωk2|τ=τkn0. |
Assume that λ(τ)=ε1(τ)+iε2(τ) is the root of Eq (4.2) near τ=τjkn, and that it {satisfies} ε1(τjkn)=0 and ε2(τjkn)=ωk2, k=1,2,3, j∈N0, n∈N0.
Theorem 4.3. If ykn=(ωk2)2 and f′(ykn)≠0, then we have that dReλ(τ)dτ|λ=iωn0≠0.
Proof. Differentiating Eq (4.2) for τ, we can get
(dλ(τ)dτ)−1=(3λ2+2q2nλ+q1n)eλτ+q3nλ(q3nλ+q4n)−τλ. | (4.9) |
Substituting λ=iωn0 into Eq (4.9), we have
(dλ(τ)dτ)−1|λ=iωn0=(−3ω2n0+2q2niωn0+q1n)(cosωn0τ+isinωn0τ)+q3n−q3nω2n0+iωn0q4n+τiωn0; |
then,
Re(dλ(τ)dτ)−1|λ=iωn0=xkny′(xkn)[q3n(ω2n)]2+(ωn0q4n)2≠0. |
Thus, we can get
sign{dReλ(τ)dτ|λ=iωn0}=sign{Re(dλ(τ)dτ)−1|λ=iωn0}≠0. |
The proof of Theorem 4.3 is completed.
By computing as shown above, we have the following conclusion.
Theorem 4.4. Under the conditions (H1)–(H3), for the the positive equilibrium E∗ of the model (1.6), we have the following:
(i) if the condition (H10) is true and p0n≥0, then E∗ is locally asymptotically stable for all τ>0;
(ii) if p0n<0 or p0n≥0, y∗1>0, f(y∗1)≤0 and the condition (H11) holds, then E∗ is locally asymptotically stable when τ∈[0,τn0); but, E∗ is unstable when τ>τn0;
(iii) if the conditions in (ii) are all satisfied and f′(ykn)≠0, then the spatially homogeneous Hopf bifurcation occurs at E∗ when τ=τ0 and n=0; and, the spatially inhomogeneous Hopf bifurcation occurs at E∗ when τ=τn0 and n>0.
Let τn=τjkn, ωn=ωk2n and τ=τn+μn, μn∈R. Thus, μn=0 is a Hopf bifurcation value of the model (1.6). Let t→tτ; then, the model (1.6) can be expressed as an FDE in C=C([−1,0],R3), as follows:
˙u(t)=τnDΔu(t)+L(τn)ut+Fn(μn,ut), | (4.10) |
where L(θ):C→X,Fn(μn,ut):C→X satisfies
L(θ)(φ)=θ(b11φ1(0)+b12φ2(0)b21φ1(0)+b22φ2(0)+b23φ3(0)+cφ1(−1)b32φ2(0)+b33φ3(0)) |
and
Fn(μn,φ)=μnDΔφ(0)+L(μn)φ+F(μn,φ), | (4.11) |
where F(μn,φ) is defined in Eq (3.15), L1 and L2 are defined in Eq (3.4) and φ=(φ1,φ2,φ3)T∈C.
The linear equation of Eq (4.10) at O(0, 0, 0) is
˙u(t)=τnDΔu(t)+L(τn)ut. | (4.12) |
Let Λ={iωnτn,−iωnτn} and zt(θ)∈C=C([−1,0],R3); we consider the following FDE:
˙z(t)=L(τn)(zt). | (4.13) |
On the basis of the Riesz representation theorem, there exists a 3×3 matrix function ηn(θ,μ) (−1≤θ≤0)∈C([−1,0],R3), and it satisfies
L(τn)(φ)=∫0−1dηn(θ,μn)φ(θ),φ∈C([−1,0],R3). |
Let
ηn(θ,μn)=(τn+μn)L1δ(θ)−(τn+μn)L2δ(θ+1), |
where δ(θ) is the Dirac delta function.
We set C∗=C([0,1],R3∗), and R3∗ is the three - dimensional vector space of row vectors. The bilinear inner product is
(ψ(s),φ(θ))=ψ(0)φ(0)−∫0−1∫θ0ψ(ξ−θ)dηn(θ)φ(ξ)dξ,ψ(s)∈C∗,φ(θ)∈C. | (4.14) |
An(τn) describes the infinitesimal generator of the semigroup induced by the solutions of Eq (4.13), and A∗n(τn) denotes the formal adjoint generator of An(τn) satisfying Eq (4.14). Let V and V∗ denote the center spaces of the generators An(τn) and A∗n(τn) corresponding to Λ, respectively. Therefore, V∗ is the adjoint space of V, dimV=dimV∗.
Lemma 4.1. Let
V1=iωn−b11b12,V2=b32(iωn−b11)b12(iωn−b33),V∗1=iωn−b11b21+ce−iωnτn,V∗2=b23(iωn−b11)(iωn−b33)(b21+ce−iωnτn); |
then, p1(θ)=(1,V1,V2)Teiωnτnθ and p2(θ)=¯p1(θ), −1≤θ≤0 form the basis of V associated with Λ; p∗1(s)=(1,V∗1,V∗2)e−iωnτns and p∗2(s)=¯p∗1(s), 0≤s≤1 form the basis of V∗ associated with Λ.
Denote Φ=(Φ1,Φ2) and Ψ∗=(Ψ∗1,Ψ∗2)T, where
Φ1(θ)=p1(θ)+p2(θ)2=(Re{eiωnτnθ},Re{V1eiωnτnθ},Re{V2eiωnτnθ})T,θ∈[−1,0],Φ2(θ)=p1(θ)−p2(θ)2i=(Im{eiωnτnθ},Im{V1eiωnτnθ},Im{V2eiωnτnθ})T,θ∈[−1,0],Ψ1(s)=p∗1(s)+p∗2(s)2=(Re{e−iωnτns},Re{V∗1e−iωnτns},Re{V∗2e−iωnτns}),s∈[0,1],Ψ2(s)=p∗1(s)−p∗2(s)2i=(Im{e−iωnτns},Im{V∗1e−iωnτns},Im{V∗2e−iωnτns}),s∈[0,1]. |
Suppose that (Ψ∗,Φ)=(Ψ∗i,Φj)(i,j=1,2.) is the basis Ψ of V∗, which satisfies
Ψ=(Ψ1,Ψ2)T=(Ψ∗,Φ)−1Ψ∗. |
Thus, we have that (Ψ,Φ)=I2×2.
Let fn=(ξ1n,ξ2n,ξ3n), where ξ1n=(cosnlx,0,0)T, ξ2n=(0,cosnlx,0)T and ξ3n=(0,0,cosnlx)T. ξjn (j=1,2,3) denotes the eigenfunctions on R3 of the eigenvalues −(nl)2, n=0,1,2⋯. Define cn⋅fn=c1ξ1n+c2ξ2n+c3ξ3n, cn=(c1,c2,c3)T,cj∈R, j=1,2,3, and the center space of Eq (4.12) is written as
PCNφ=Φ(Ψ,⟨φ,fn⟩)⋅fn, |
where φ∈C, C=PCNC⊕PsC and PsC expresses the complementary subspace of PCNC.
According to [57] and [59], the center space of the linear model of (4.10) with μn=0 is expressed as PCNC, where
PCNC={12(p1(θ)z+p2(θ)ˉz)⋅fn,z∈C}. |
Thus, the solution of the model (4.10) can be written as
ut=12(p1(θ)z+p2(θ)ˉz)⋅fn+Q(z(t),ˉz(t))(θ), |
where Q(z(t),ˉz(t))(θ)=W(z+ˉz2,iz−ˉz2,0), z=x1−ix2.
By Wu[59], z satisfies
˙z=iωnτnz+gn(z,ˉz), | (4.15) |
where
gn(z,ˉz)=(Ψ1(0)−iΨ2(0))⟨Fn(0,ut),fn⟩,Ψ(0)=(Ψ1(0),Ψ2(0))T. |
Let
Q(z,ˉz)=Q20(θ)z22+Q11(θ)zˉz+Q02(θ)ˉz22+⋯, | (4.16) |
gn(z,ˉz)=˜g20(θ)z22+˜g11(θ)zˉz+˜g02(θ)ˉz22+˜g21z2ˉz2+⋯ |
and
Ψ1(0)−iΨ2(0)=(ψ1,ψ2,ψ3). |
By computing and comparing the coefficients, we have
˜g20=τn2⟨[(C11+C12V1)ψ1+(C21+C22V21+C23V22+C24V1+C25V1V2+C26e−2iωnτn+C27V1e−iωnτn)ψ2+(C31V21+C32V22+C33V1V2)ψ3]cos2(nl)x,cos(nl)x⟩, |
˜g11=τn4⟨[(2C11+C12(V1+ˉV1))ψ1+(2C21+2C22V1ˉV1+2C23V2ˉV2+C24(V1+ˉV1)+C25(V1ˉV2+ˉV1V2)+2C26+C27(V1eiωnτn+ˉV1e−iωnτn))ψ2+(2C31V1ˉV1+2C32V2ˉV2+C33(V1ˉV2+ˉV1V2))ψ3]cos2(nl)x,cos(nl)x⟩,˜g21=τn{⟨[C11(Q(1)20(0)+2Q(1)11(0))+C12(12Q(1)20(0)ˉV1+V1Q(1)11(0)+Q(2)11(0)+12Q(2)20(0))]cos(nl)x,cos(nl)x⟩ψ1+⟨[C21(Q(1)20(0)+2Q(1)11(0))+C22(Q(2)20(0)+2Q(2)11(0))+C23(Q(3)20(0)+2Q(3)11(0))+C24(12Q(1)20(0)ˉV1+V1Q(1)11(0)+Q(2)11(0)+12Q(2)20(0))+C25(12Q(2)20(0)ˉV2+Q(1)11(0)V2+V1Q(3)11(0)+12Q(3)20(0)ˉV1)+C26(eiωnτnQ(1)20(−1)+2e−iωnτnQ(1)11(−1))+C27(12ˉV1Q(1)20(−1)+V1Q(1)11(−1)+e−iωnτnQ(2)11(0)+12eiωnτnQ(2)20(0))]cos(nl)x,cos(nl)x⟩ψ2+⟨[C31(Q(2)20(0)+2Q(1)11(0))+C33(12ˉV2Q(2)11(0)+V2Q(2)11(0)+V1Q(3)11(0)+12Q(3)20(0)ˉV1)]cos(nl)x,cos(nl)x⟩ψ3}. |
We know that ∫π0cos3(nl)xdx=0 and ˜g02=¯˜g20. Therefore, we can get that ˜g20=˜g11=˜g02=0 when n=1,2,3⋯. When n=0, we have
˜g20=τn2[(C11+C12V1)ψ1+(C21+C22V21+C23V22+C24V1+C25p1V2+C26e−2iωnτn+C27V1e−iωnτn)ψ2+(C31V21+C32V22+C33V1V2)ψ3],˜g11=τn4[(2C11+C12(V1+ˉV1))ψ1+(2C21+2C22V1ˉV1+2C23V2ˉV2+C24(V1+ˉV1)+C25(V1ˉV2+ˉV1V2)+2C26+C27(V1eiωnτn+ˉV1e−iωnτn))ψ2+(2C31V1ˉV1+2C32V2ˉV2+C33(V1ˉV2+ˉV1V2))ψ3]. | (4.17) |
Considering the expression of ˜g21, it contains Q20(θ) and Q11(θ), so we must compute Q20(θ) and Q11(θ). Seeking the derivative on both sides of Eq (4.16), we have
˙Q(z,ˉz)=Q20z+Q11˙zˉz+Q11z˙ˉz+Q02˙ˉzˉz+⋯, | (4.18) |
AτnQ(z,ˉz)=AτnQ20z22+AτnQ11zˉz+AτnQ02ˉz22+⋯. | (4.19) |
From Wu [59], Q(z,ˉz) satisfies
˙Q(z,ˉz)=AτnQ(z,ˉz)+S(z,ˉz), | (4.20) |
where
S(z,ˉz)=S20z22+S11zˉz+S02ˉz22+⋯=X0Fn(ut,0)−Φ(Ψ,⟨X0Fn(ut,0),fn⟩)⋅fn, |
and
X0(θ)={I,θ=0,0,−1≤θ<0, |
Sij∈PSC,i+j=2. |
From Eq (4.15) and Eqs (4.17)–(4.20), we can get
{(2iωnτn−Aτn)Q20=S20,−AτnQ11=S11. |
Because Aτn has only two characteristic roots with a zero real part, i.e., ±iωnτn, Eq (4.20) has a unique solution Qij(i+j=2) in PSC, which satisfies
{Q20=(2iωnτn−Aτn)−1S20,Q11=−A−1τnS11. | (4.21) |
From Eq (4.20), we can get that, for θ∈[−1,0),
S(z,ˉz)=−Φ(θ)Ψ(0)⟨Fn(ut,0),fn⟩⋅fn=−τn2[(p1(θ)˜g20+¯˜g02p2(θ))⋅fn⋅z22+(p1(θ)˜g11+¯˜g11p2(θ))⋅fn⋅zˉz+(p1(θ)˜g02+¯˜g20p2(θ))⋅fn⋅ˉz22]+⋯. | (4.22) |
Comparing the coefficients in Eqs (4.20) and (4.22), when θ∈[−1,0), we can receive
S20(θ)=−τn2(p1(θ)˜g20+¯˜g02p2(θ))cos(nl)x,S11(θ)=−τn2(p1(θ)˜g11+¯˜g11p2(θ))cos(nl)x. | (4.23) |
When θ=0, we have
S20(0)=τn2(C11+C12p1C21+C22p21+C23p22+C24p1+C25p1p2+C26e−2iωnτn+C27p1e−iωnτnC31p21+C32p22+C33p1p2)cos2(nl)x−τn2(p1(θ)˜g20+¯˜g02p2(θ))cos(nl)x, | (4.24) |
S11(0)=τn4(2C11+C12(p1+ˉp1)2C21+2C22p1ˉp1+2C23p2ˉp2+C24(p1+ˉp1)+C25(p1ˉp2+ˉp1p2)+2C26+C27(p1eiωnτn+ˉp1e−iωnτn)2C31p1ˉp1+2C32p2ˉp2+C33(p1ˉp2+ˉp1p2))cos2(nl)x−τn2(p1(θ)˜g11+¯˜g11p2(θ))cos(nl)x. | (4.25) |
Using Eqs (4.24) and (4.25), we can get Q20(0), Q11(0), Q20(−1) and Q11(−1). Because p1(θ)=p1(0)eiωnτnθ, θ∈[−1,0), from Eqs (4.21)–(4.25), we have
Q20(θ)=i2[˜g20ωnτnp1(0)eiωnτnθ+¯˜g023ωnτnp1(0)e−iωnτnθ]+e2iωnτnθE3,Q11(θ)=i2[˜g11ωnτnp1(0)eiωnτnθ−¯˜g11ωnτnp1(0)e−iωnτnθ]+E4, |
where E3=(E(1)3,E(2)3,E(3)3)∈R3 and E4=(E(1)4,E(2)4,E(3)4)∈R3 satisfy
E3=12(2iωn−b11−b120−ce−2iωnτn−b212iωn−b22−b230−b322iωn−b33)−1 |
×(C11+C12V1C21+C22V21+C23V22+C24V1+C25V1V2+C26e−2iωnτn+C27V1e−iωnτnC31V21+C32V22+C33V1V2)cos2(nl)x, |
E4=(−b11−b120−c−b21−b22−b230−b32−b33)−1 |
×(2C11+C12(V1+ˉV1)2C21+2C22V1ˉV1+2C23V2ˉV2+C24(V1+ˉV1)+C25(V1ˉV2+ˉV1V2)+2C26+C27(V1eiωnτn+ˉV1e−iωnτn)2C31V1ˉV1+2C32V2ˉV2+C33(V1ˉV2+ˉV1V2))cos2(nl)x. |
Therefore, we can get the following values:
{c2(0)=i2ωnτn(˜g11˜g20−2|˜g11|2−|˜g02|23)+˜g212,μ3=−Re{c2(0)}Re{λ′(τn)},μ4=2Re{c2(0)},T2=−Im{c2(0)}+μ3Im{λ′(τn)}ωnτn, | (4.26) |
which determine the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions on the center manifold at τ=τn.
Theorem 4.5. According to Eq (4.26), we have the following conclusions.
(i) μ3 determines the direction of Hopf bifurcation: if μ3>0, then the Hopf bifurcation is supercritical; if μ3<0, then the Hopf bifurcation is subcritical and the bifurcating periodic solutions exist when τ>τn;
(ii) μ4 determines the stability of the bifurcating periodic solutions: if μ4>0, then the bifurcating periodic solutions are unstable; if μ4<0, then the bifurcating periodic solutions are stable;
(iii) T2 determines the period of the bifurcating periodic solutions: if T2>0, then the bifurcating periodic solutions increase; if T2<0, then the bifurcating periodic solutions decrease.
With the help of Matlab software, the stability of the positive equilibrium E∗(P∗,Z∗,F∗) was simulated with the given values of all parameters in order to confirm the previous theoretical results.
First, we assume that P(0)=0.8, Z(0)=40 and F(0)=0.1 for the model (2.1). And, we take the values of all other parameters as follows: β1=0.016, β2=0.7, γ1=0.0875, γ2=0.075, α1=0.1, α2=0.2, a=0.5, b=0.25, g1=0.1, g2=0.2, m=0.8 and δ=0.35. According to Theorem 2.1, we can know that the model (2.1) has one unique positive equilibrium E∗(0.5043,31.1137,0.1191). From Theorem 2.5, E∗ is locally asymptotically stable (see Figure 3).
In the second section, we obtained Theorem 2.6 by referring to [55]. Now, we will verify some conclusions by taking m as the bifurcation parameter. When δ=0.2, the critical value is m∗=0.82, which satisfies M6(m∗)>0, T(m∗)=0, M8(m∗)>0 and dTds|m=m∗≠0. That is, all conditions in Theorem 2.6 are satisfied. Therefore, Hopf bifurcation occurs at E∗ when m=m∗. Meanwhile, we can obtain that P(t) reaches a maximum value and Z(t) and F(t) are always 0 when m≥0.94. Therefore, we can get the bifurcation diagram as m changes (see Figure 4). If we choose δ as the bifurcation parameter when m=0.75, we can get the critical value δ∗=0.33. Therefore, Hopf bifurcation occurs at E∗ when δ=δ∗. We can obtain that P(t) reaches a maximum value of 1 and Z(t) and F(t) are always 0 when δ≥0.49. Therefore, we can get the bifurcation diagram as δ changes (see Figure 5). From Figure 5, when other parameter values are fixed, the density of zooplankton and fish will decrease to 0 whether the refuge capacity of phytoplankton or the probability of toxin release of phytoplankton-produced toxic substances increase to some certain value. Properly increasing the shelter capacity of phytoplankton and the rate of toxin release of by phytoplankton can stabilize the population and reach a stable state. Then, the plankton and fish populations will always exist.
For the model (3.1), we assume that P(0)=0.5, Z(0)=30 and F(0)=0.115. When m=0.8 and δ=0.25, we have that τ0=4.9397, and the model (2.1) has one unique positive equilibrium E∗(0.2671,35.1379,0.1197) according to Theorem 2.1. From Theorem 3.3, E∗ is locally asymptotically stable when τ∈[0,4.9397], but Hopf bifurcation occurs when τ∈[4.9397,+∞). From Eq (3.35), we can know that c1(0)=−512.86−540.47i<0, μ1=1457>0, μ2=−1025.7<0 and T1=153.5416>0. Thus, the Hopf bifurcation is supercritical, the bifurcating periodic solution is stable and the period of the bifurcating periodic solutions is increasing, which can be seen in Figure 6 (τ=1) and Figure 7 (τ=10). Here, we give the delay bifurcation diagram (see Figure 8). This means that, if the mature delay exceeds the critical value, the model transitions to unstable from stable. At this moment, the model has a Hopf bifurcation near the equilibrium and unstable behavior occurs among populations. In other words, the presence of the mature delay can destabilize the plankton-fish population.
For the model (1.6), we choose l=1, that is, x∈(0,π). The values of other parameters as follows: d1=0.2, d2=0.5, d3=0.1, m=0.84 and δ=0.25. The positive equilibrium is E∗(0.3794,38.9575,0.1202). From Eq (4.8), we have that τn0=6.2832. Based on Theorem 4.1, E∗ is also locally asymptotically stable when τ=0 (see Figure 9), and E∗ is locally asymptotically stable when τ=3.5<τn0=6.2832 (see Figure 10). But, E∗ is unstable when τ=25>τn0=6.2832 (see Figure 11). And, we can compute that c2(0)=−1141.6−2543.9i<0, μ3=7257.8>0, μ4=−2283.3<0 and T2=646.8790>0. From Theorem 4.5, the Hopf bifurcation is supercritical, the bifurcating periodic solution is stable and the period of the bifurcating periodic solutions is increasing. For the reaction-diffusion model, we can know that the model will transitions to unstable from stable if the mature delay exceeds the critical value. At this moment, the model has a spatially homogeneous Hopf bifurcation or spatially inhomogeneous Hopf bifurcation near the equilibrium and unstable behavior occurs between the populations. At this time, the presence of the mature delay can destabilize the plankton-fish population.
By the theoretical conclusions and numerical simulation, we not only find that the existence of delay will deteriorate the system stability under some conditions, but also that the refuge of the prey and the release of toxins will cause the stability of system be damaged in a reaction-diffusion model with delay, even causing Hopf bifurcation to occur at the positive equilibrium E∗.
In our paper, we establish a phytoplankton-zooplankton-fish model with mature delay and population diffusion by considering the refuge of phytoplankton, C-M functional response and Holling II functional response. In [18], the authors found that the refuge affects the stability of the positive equilibrium. However, in our paper, we not only analyzed the effect of refuge, but also studied the effects of diffusion and delay on the model; we obtained that the stability of the system may be destroyed due to the existence of delay. In [19], the authors analyzed Hopf bifurcation caused by delay. However, in our paper, we not only obtained the properties of Hopf bifurcation induced by delay, but also the influence of prey refuge on the population. We determined that the existence of prey refuge can also lead to Hopf bifurcation.
After the parameters in the model were selected, the existence and stability of the equilibrium were analyzed. First, we chose m as a bifurcation parameter to study the dynamical behavior as m changes in the model (2.1). Meanwhile, we consider the effect of the parameter δ on the positive equilibrium in the model (2.1). Through analysis, it could be obtained that the model undergoes a Hopf bifurcation when m=m∗ or δ=δ∗. We found that, when other parameter values are fixed, the densities of zooplankton and fish will decrease to 0 regardless whether the refuge capacity of phytoplankton or the probability of toxin release of phytoplankton-produced toxic substances increase to a certain value. And, we chose the time delay τ as the bifurcation parameter and discussed the dynamical behavior of the model without diffusion, or with diffusion, respectively. We give the direction of the Hopf bifurcation and the stability of the bifurcating periodic solution by the center manifold theorem and normal form theory. We found that the model transitions to unstable from stable when the mature delay exceeds the critical value. At this moment, the model has a spatially homogeneous Hopf bifurcation or spatially inhomogeneous Hopf bifurcation near the positive equilibrium and unstable behavior occurs between the populations. In a word, the existence of time delay has a great influence on such a model. Meanwhile, we used Matlab software for numerical simulation to prove our theoretical results.
In this paper, we have discussed the influence of factors such as prey refuge, the disturbance between predators, time delay and diffusion on the model. However, in nature, there are external factors to influence the model, such as changing temperature, environmental pollution, human activities and noise. We did not take these influencing factors into account. Therefore, in the future work, we will introduce the influence of environmental pollution on the model and analyze the dynamical behavior of the phytoplankton-zooplankton model under the influence of environmental pollution. The model is
{∂P∂t=d1ΔP+r1P(1−PK1)−β1(1−m)PZ1+a1(1−m)P+cZ−m1P3,x∈Ω,t>0,∂Z∂t=d2ΔZ+r2Z(1−ZK2)+β2(1−m)PZ1+a1(1−m)P+cZ−δP(t−τ)Za2+P(t−τ)−m2Z2−gZ,x∈Ω,t>0,Px(x,t)=Zx(x,t)=0,x∈∂Ω,t>0,P(x,t)>0,Z(x,t)>0,x∈Ω,t∈[−τ,0], |
where m1 and m2 are the effects of environmental pollution on phytoplankton and zooplankton, respectively. We leave this work for the future.
This work was supported by the National Natural Science Foundation of China (Grant Nos. 12161054, 11661050 and 11861044), the National Natural Science Foundation of Gansu Province (Grant No. 20JR10RA156), the Doctoral Foundation of Lanzhou University of Technology and the HongLiu First-Class Disciplines Development Program of Lanzhou University of Technology.
The authors declare that they have no conflict of interest.
[1] | Z. H. Yang, Y. M. Chu, An optimal inequalities chain for bivariate means, J. Math. Inequal., 9 (2015), 331-343. |
[2] |
K. B. Stolarsky, Generalizations of the Logarithmic Mean, Math. Mag., 48 (1975), 87-92. doi: 10.1080/0025570X.1975.11976447
![]() |
[3] |
E. B. Leach, M. C. Sholander, Extended mean values, Amer. Math. Monthly, 85 (1978), 84-90. doi: 10.1080/00029890.1978.11994526
![]() |
[4] | C. Gini, Diuna formula comprensiva delle media, Metron, 13 (1938), 3-22. |
[5] | Z. H. Yang, On the homogeneous functions with two parameters and its monotonicity, J. Inequal. Pure Appl. Math., 6 (2005), 1-11. |
[6] | Z. H. Yang, On the monotonicity and log-convexity of a four-parameter homogeneous mean, J. Inequal. Appl., 2008 (2008), 1-12. |
[7] | A. Witkowski, Comparison theorem for two-parameter means, Math. Inequal. Appl., 12 (2009), 11-20. |
[8] | F. Qi, Logarithmic convexities of the extended mean values, Proc. Amer. Math. Soc., 130 (2002), 1787-1796. |
[9] | Z. H. Yang, On the log-convexity of two-parameter homogeneous functions, Math. Inequal. Appl., 10 (2007), 499-516. |
[10] |
Zs. Páles, Inequalities for sums of powers, J. Math. Anal. Appl., 131 (1988), 265-270. doi: 10.1016/0022-247X(88)90204-1
![]() |
[11] |
Zs. Páles, Inequalities for differences of powers, J. Math. Anal. Appl., 131 (1988), 271-281. doi: 10.1016/0022-247X(88)90205-3
![]() |
[12] | L. Losonczi, Zs. Páles, Minkowki's inequality for two variable Gini means, Acta Sci. Math. Szeged, 62 (1996), 413-425. |
[13] |
L. Losonczi, Zs. Páles, Minkowki's inequality for two variable difference means, Proc. Amer. Math. Soc., 126 (1998), 779-791. doi: 10.1090/S0002-9939-98-04125-2
![]() |
[14] |
E. Neuman, Zs. Páles, On comparison of Stolarsky and Gini means, J. Math. Anal. Appl., 278 (2003), 274-284. doi: 10.1016/S0022-247X(02)00319-0
![]() |
[15] |
Y. M. Li, B. Y. Long, Y. M. Chu, Sharp bounds by the power mean for the generalized Heronian mean, J. Inequal. Appl., 2012 (2012), 1-9. doi: 10.1186/1029-242X-2012-1
![]() |
[16] |
M. Raïsouli, J. Sándor, Sub-super-stabilizability of certain bivariate means via mean-convexity, J. Inequal. Appl., 2016 (2016), 1-13. doi: 10.1186/s13660-015-0952-5
![]() |
[17] | Z. H. Yang, On converses of some comparison inequalities for homogeneous means, Hacet. J. Math. Stat., 46 (2017), 629-644. |
[18] | Z. H. Yang, New sharp bounds for identric mean in terms of logarithmic mean and arithmetic mean, J. Math. Inequal., 6 (2012), 533-543. |
[19] | A. Begea, J. Bukor, J. T. Tóhb, On (log-) convexity of power mean, Anna. Math. Inform., 42 (2013), 3-7. |
[20] |
L. Matejíčka, Short note on convexity of power mean, Tamkang J. Math., 46 (2015), 423-426. doi: 10.5556/j.tkjm.46.2015.1808
![]() |
[21] | I. Pinelis, L'Hospital type rules for oscillation, with applications, J. Inequal. Pure Appl. Math., 2 (2001), 1-24. |
[22] |
G. D. Anderson, M. Vamanamurthy, M. Vuorinen, Monotonicity rules in calculus, Amer. Math. Monthly, 113 (2006), 805-816. doi: 10.1080/00029890.2006.11920367
![]() |
[23] | Z. H. Yang, A new way to prove L'Hospital Monotone Rules with applications, arXiv:1409.6408, 2014. Available from: https://arxiv.org/abs/1409.6408. |
[24] |
Z. H. Yang, Y. M. Chu, M. K. Wang, Monotonicity criterion for the quotient of power series with applications, J. Math. Anal. Appl., 428 (2015), 587-604. doi: 10.1016/j.jmaa.2015.03.043
![]() |
[25] | Z. H. Yang, Y. M. Chu, A monotonicity property involving the generalized elliptic integral of the first kind, Math. Inequal. Appl., 20 (2017), 729-735. |
[26] | Z. H. Yang, W. Zhang, Y. M. Chu, Sharp Gautschi inequality for parameter 0 < p < 1 with applications, Math. Inequal. Appl., 20 (2017), 1107-1120. |
[27] |
Z. H. Yang, J. F. Tian, The monotonicity rules for the ratio of two Laplace transforms with applications, J. Math. Anal. Appl., 470 (2019), 821-845. doi: 10.1016/j.jmaa.2018.10.034
![]() |
[28] |
Z. H. Yang, K. F. Tin, Q. Gao, The monotonicity of ratios involving arctangent function with applications, Open Math., 17 (2019), 1450-1467. doi: 10.1515/math-2019-0098
![]() |
[29] | Z. H. Yang, Estimates for Neuman-Sándor mean by power means and their relative errors, J. Math. Inequal., 7 (2013), 711-726. |
[30] | Z. H. Yang, Some monotonictiy results for the ratio of two-parameter symmetric homogeneous functions, Int. J. Math. Math. Sci., 2009 (2019), 1-12. |
[31] | Z. H. Yang, Log-convexity of ratio of the two-parameter symmetric homogeneous functions and an application, J. Inequal. Spec. Func., 2010 (2010), 16-29. |
1. | Muhammad Uzair Awan, Sadia Talib, Muhammad Aslam Noor, Yu-Ming Chu, Khalida Inayat Noor, On post quantum estimates of upper bounds involving twice (p,q)-differentiable preinvex function, 2020, 2020, 1029-242X, 10.1186/s13660-020-02496-5 | |
2. | Shu-Bo Chen, Saima Rashid, Muhammad Aslam Noor, Zakia Hammouch, Yu-Ming Chu, New fractional approaches for n-polynomial P-convexity with applications in special function theory, 2020, 2020, 1687-1847, 10.1186/s13662-020-03000-5 | |
3. | Yu-Ming Chu, Muhammad Uzair Awan, Sadia Talib, Muhammad Aslam Noor, Khalida Inayat Noor, New post quantum analogues of Ostrowski-type inequalities using new definitions of left–right (p,q)-derivatives and definite integrals, 2020, 2020, 1687-1847, 10.1186/s13662-020-03094-x | |
4. | Saad Ihsan Butt, Muhammad Umar, Saima Rashid, Ahmet Ocak Akdemir, Yu-Ming Chu, New Hermite–Jensen–Mercer-type inequalities via k-fractional integrals, 2020, 2020, 1687-1847, 10.1186/s13662-020-03093-y | |
5. | Hua Wang, Humaira Kalsoom, Hüseyin Budak, Muhammad Idrees, Ahmet Ocak Akdemir, q -Hermite–Hadamard Inequalities for Generalized Exponentially s , m ; η -Preinvex Functions, 2021, 2021, 2314-4785, 1, 10.1155/2021/5577340 | |
6. | Humaira Kalsoom, Muhammad Aamir Ali, Muhammad Idrees, Praveen Agarwal, Muhammad Arif, Adrian Neagu, New Post Quantum Analogues of Hermite–Hadamard Type Inequalities for Interval-Valued Convex Functions, 2021, 2021, 1563-5147, 1, 10.1155/2021/5529650 | |
7. | Humaira Kalsoom, Miguel Vivas-Cortez, Muhammad Idrees, Praveen Agarwal, New Parameterized Inequalities for η-Quasiconvex Functions via (p, q)-Calculus, 2021, 23, 1099-4300, 1523, 10.3390/e23111523 | |
8. | Saad Ihsan Butt, Hüseyin Budak, Muhammad Tariq, Muhammad Nadeem, Integral inequalities for n ‐polynomial s ‐type preinvex functions with applications , 2021, 44, 0170-4214, 11006, 10.1002/mma.7465 | |
9. | Humaira Kalsoom, Miguel Vivas-Cortez, Muhammad Amer Latif, Trapezoidal-Type Inequalities for Strongly Convex and Quasi-Convex Functions via Post-Quantum Calculus, 2021, 23, 1099-4300, 1238, 10.3390/e23101238 | |
10. | YONGFANG QI, GUOPING LI, FRACTIONAL OSTROWSKI TYPE INEQUALITIES FOR (s,m)-CONVEX FUNCTION WITH APPLICATIONS, 2023, 31, 0218-348X, 10.1142/S0218348X23501281 |