# Order and Chaos in a Prey-Predator Model Incorporating Refuge, Disease, and Harvesting

Authors: Dahlia Khaled Bahlool; Huda Abdul Satar; Hiba Abdullah Ibrahim
Journal: Journal of Applied Mathematics
 (2020)

---

## Abstract

In this paper, a mathematical model consisting of a prey-predator system incorporating infectious disease in the prey has been proposed and analyzed. It is assumed that the predator preys upon the nonrefugees prey only according to the modified Holling type-II functional response. There is a harvesting process from the predator. The existence and uniqueness of the solution in addition to their bounded are discussed. The stability analysis of the model around all possible equilibrium points is investigated. The persistence conditions of the system are established. Local bifurcation analysis in view of the Sotomayor theorem is carried out. Numerical simulation has been applied to investigate the global dynamics and specify the effect of varying the parameters. It is observed that the system has a chaotic dynamics.


---

## Body

## 1. Introduction

The interaction between a predator and their prey is one of the most important interaction types in ecology as well as applied mathematics due to its wide existence in our real life. Although such interaction models may appear to be simple mathematically, they are very challenging and complicated [1]. Since the classical Lotka (1925) and Volterra (1926) model on prey-predator interaction, which was developed on the basis of chemical principle of mass action, many research works have been done in literature. Later, the simple growth and decay terms in the Lotka-Volterra model were replaced by more complex functions such as the logistic growth; prey-dependent functional responses such as Holling’s I, II, and III types; ratio-dependent functional response, and predator-dependent type of functional response, see, for example, [2–8] and the references therein. It is well known that the two species models are very restrictive in their dynamics and can give rise to only two basic dynamics approaches to equilibrium or to a limit cycle, while ecological situations demonstrate very complex dynamics in nature. Accordingly, multispecies systems comprising of three species or more was considered to study the complex ecological community [9–13]. Furthermore, due to the existence of different important factors in the environment, such as stage structure, refuge, cannibalism, harvesting, toxicant, and scavenger, different prey-predator models incorporating such factors have been proposed and studied in literature, see, for example, [14–17].On the other hand, the spread of diseases between the individuals of species is an important subject for study, which is known as an epidemiology, due to the existence of a variety of diseases in the environment. It is well known that after the SIR epidemic model of Kermack–Mckendrick [18], the field of epidemiology has received a lot of attention from the researchers, see [19] and the references therein. Recently, Zhou and Yao [20] proposed and studied a host-vector epidemic model with a stage structure for the vector. Further articles have been done also in the field of epidemiology. In fact, recently, there is increasing interest from the scientists by combining ecology and epidemiology in one field that is called ecoepidemiology, which has become a major field of study, see, for example, [21–24] and the references therein. Recently, Wang et al. [25] proposed and studied a generalized ecoepidemiological system with prey refuge. The saturation incidence kinetics and a generalized functional response are used to describe the contact process and the predation process, respectively. They show that the prey refuge can control the spread of disease by the relative level of prey refuge.Keeping the above studies in view, an ecoepidemiological model is proposed consisting of prey-predator incorporating a prey refuge, infectious disease of SI type in prey species only, and the predator preys upon the prey (susceptible and infected) according to the modified Holling type-II functional responses. It is assumed further that there is a harvesting process on the predator species due to external resources.


## 2. Model Formulation

In this section, a prey-predator system incorporating disease in a prey is formulated mathematically. It is assumed that there is a refuge for prey to protect themselves against predation and the harvesting process on predator. The following hypotheses are adopted in constricting the mathematical model that describes the dynamics of the above real-world ecoepidemiological system.(1)
In the presence of infectious disease, the prey population is divided into two compartments, namely, susceptible prey that is denoted to their density at timet by ST and infected prey, which is denoted to their density at time T by IT. Therefore, at time T, the total prey population density is XT=ST+IT(2)
It is assumed that in the absence of predator that is denoted to their density at timeT by YT, the susceptible prey grows logistically with carrying a capacity K>0 and an intrinsic growth rate r>0, while the infected prey cannot reproduce and rather than that, it competes for food and space(3)
The availability of refuge for the prey in the environment keepsm∈0,1 portion of them safe from predation and left 1−m from the total density of the prey facing predator(4)
The disease transfers from infected to susceptible individuals by contact with an infected rate given byβ>0(5)
There is a harvesting process on the predator, due to an external force, represented by a ratel1E/l2E+l3Y, where E>0 is the effort applied to harvest the prey; l1>0 is the catchability coefficient while l2>0 and l3>0 are suitable constants(6)
The predator consumed the prey according to the modified Holling type-II with attack ratesa1>0 and a2>0 for susceptible prey and infected prey, respectively, half saturation constant b>0, and infected prey’s preference rate c>0, while the conversion rate is given by e>0(7)
Finally, there is a disease death rate given byd1>0 and the natural death rate for predator d2>0According to the above hypotheses, the dynamics of the above described ecoepidemiological real-world system can be represented by the following set of differential equations.(1)dSdT=rS1−S+IK−βSI−a11−mSIb+1−mS+c1−mI,dIdT=βSI−a21−mIYb+1−mS+c1−mI−d1I,dYdT=ea1S−a2I1−mYb+1−mS+c1−mI+d2Y−l1EYl2E+l3Y,with initial data S0≥0, I0≥0, and Y0≥0. Clearly, the interaction functions in system (1) are continuously differential functions on the positive octant.Hence, they are Lipschitzian. Therefore, the solution of system (1) exits and is unique. Now, for further simplifying the above system, the following dimensionless variables and constants are used.
(2)rT=t,SK=x,IK=y,a1YrK=z,w1=βKr,w2=bK,w3=a2a1w4=d1r,w5=ea1r,w6=dr,w7=l1a1El3r2K,w8=l2Ea1l3rKThus, system (1) can be turned into the following dimensionless form:
(3)dxdt=x1−x−y−w1y−1−mzw2+1−mx+c1−my=xf1x,y,z,dydt=yw1x−w31−mzw2+1−mx+c1−my−w4=yf2x,y,z,dzdt=zw51−mx+w3yw2+1−mx+c1−my−w6−w7w8+z=zf3x,y,z,with initial data x0,y0,z0 belonging to ℝ+3=x,y,z∈ℝ3:xt>0,yt>0,zt>0. Clearly, according to the biological meaning of the parameters, we have always that w5<1. Moreover, all solutions of system (3) initiated in ℝ+3 are uniformly bounded as shown in the following theorem.Theorem 1.
All the solutions of system (3), which initiate in ℝ+3, are uniformly bounded.Proof.
According to the first equation of system (3) with the help of the Gronwall lemma, we have xt<1. Let xt,yt,zt be any solution of system (3) initiated in ℝ+3, and let W be a function defined by Wt=xt+yt+zt. Then, by differentiating this function with respect to time, it gives
(4)dWdt=dxdt+dydt+dzdt≤2−μ1W,where μ1=min1,w4,w6. Therefore, it obtained
(5)dWdt+μ1W≤2.
Again, using the Gronwall lemma gives that ast⟶∞, Wt≤2/μ1. Thus, all the solutions of system (1) are uniformly bounded.


## 3. Stability and Persistence

In this section, the stability analysis of all possible equilibrium points of system (3) is investigated. The persistence conditions of system (3) are established. Now, straightforward computation shows that the system (1) has at most five nonnegative equilibrium points. The existence conditions for each of them can be summarized as follows: The trivial equilibrium point (TEP), say E0=0,0,0 always exists. The disease-predator-free equilibrium point (DPFEP), namely, E1=x′,0,0=1,0,0 always exists. The predator-free equilibrium point (PFEP), namely, E2=x¯,y¯,0=w4/w1,w1−w4/w11+w1,0 exists in the xy-plane provided that
(6)w4<w1.The disease-free equilibrium point (DFEP) is denoted byE3=x^,0,z^ where
(7)z^=1−x^w2+1−mx^1−m,while x^ is a positive root of the following third-order equation:
(8)α1x^3+α2x^2+α3x^+α4=0,with
(9)α1=N2w6−w5,α2=N2w5−w6+w2N2w6−w5,α3=w8N2w5−w6+w2Nw5−2w6+w22w6−N2w7,α4=−w2Nw6w8+w7−w22w6<0.Hence,E3 exists uniquely in the positive quadrant of xz-plane if one of the following sets of conditions holds.
(10)α1>0andα2>0α1>0andα3<0The coexistence equilibrium point (CEP) is denoted byE4=x∗,y∗,z∗, where
(11)y∗=w3+w4−w1+w3x∗w31+w1,z∗=w1x∗−w4w2+Nx∗+cNy∗Nw3,while x∗ is a positive root of the following third-order equation.
(12)σ1x∗3+σ2x∗2+σ3x∗+σ4=0.Here,
(13)σ1=γ1θ1−θ3w6,σ2=θ1γ2+θ2γ1−θ3γ2w6−θ4γ1w6,σ3=θ1θ6+θ2γ2−θ3θ6w6−θ4γ2w6−θ5θ3,σ4=θ2θ6−θ4θ6w6−θ3θ4,with
(14)θ1=w3w5N1−w3,θ2=w3w5Nw3+w4,θ3=Nw31+w1−cw1+w3,θ4=w2w31+w1+cNw3+w4,θ5=w7w32N1+w1,θ6=θ3+w8Nw321+w1,γ1=w1Nw31+w3−cw1+w3,γ2=w31+w1w1w2−w4N+cNw1w3+w4+cNw4w1+w3.Clearly, the CEP that is denoted byE4 exists uniquely in the interior of the positive octant provided that the following condition holds.
(15)w4w1<x∗<w3+w4w1+w3,with one set of the following sets of conditions
(16)σ1>0,σ4<0,σ2>0σ1>0,σ4<0,σ3<0σ1<0,σ4>0,σ2<0σ1<0,σ4>0,σ3>0In the following, the local stability analysis of the above equilibrium points is carried out. The Jacobian matrix of system (1) at point x,y,z can be written as
(17)Jx,y,z=uij3×3Here,
(18)u11=−x+1−m2xzw2+1−mx+c1−my2+f1x,y,z,u12=−x−w1x+c1−m2xzw2+1−mx+c1−my2,u13=−1−mxw2+1−mx+c1−my,u21=w1y+w31−m2yzw2+1−mx+c1−my2,u22=w31−m2yzw2+1−mx+c1−my2+f2x,y,z,u23=−w31−myw2+1−mx+c1−my,u31=w2w51−mz+w5c−w31−m2yzw2+1−mx+c1−my2,u32=w2w3w51−mz+w5w3−c1−m2xzw2+1−mx+c1−my2,u33=w7zw8+z2+f3x,y,z.Therefore, the eigenvalues of the Jacobian matrix atE0 are
(19)λ0x=1>0,λ0y=−w4<0,λ0z=−w6−w7w8<0.Hence, the TEP is a saddle point with unstable manifold inx-direction and stable manifold in the yz-plane.The Jacobian matrix of system (3) at the DPFEP that is represented by E1=1,0,0 is given by
(20)JE1=−1−1−w1−1−mw2+1−m0w1−w4000w51−mw2+1−m−w6−w7w8Therefore, the eigenvalues are written as(21)λ1x=−1<0,λ1y=w1−w4,λ1z=w51−mw2+1−m−w6−w7w8.It is clear that all the eigenvalues are negative and then, the DPFEP is locally asymptotically stable provided that the following conditions are hold.(22a)w1<w4,(22b)w51−mw2+1−m<w6+w7w8,The Jacobian matrix of system (3) at the PFEP is denoted by E2=x¯,y¯,0.
(23)JE2=−x¯−1+w1x¯−1−mx¯w2+1−mx¯+c1−my¯w1y¯0−w31−my¯w2+1−mx¯+c1−my¯00w51−mx¯+w3y¯w2+1−mx¯+c1−my¯−w6−w7w8Clearly, the characteristic equation ofJE2 can be written as
(24)λ22+x¯λ2+w11+w1x¯y¯w51−mx¯+w3y¯w2+1−mx¯+c1−my¯−w6−w7w8−λ2=0.Accordingly, all the eigenvalues ofJE2 have negative real parts and then, the PFEP is locally asymptotically stable provided that the following condition holds
(25)w51−mx¯+w3y¯w2+1−mx¯+c1−my¯<w6+w7w8.The Jacobian matrix of system (3) at the DFEP that is given by E3 can be written as
(26)JE3=bij3×3,where
(27)b11=−1+1−m2z^w2+1−mx^2x^,b12=−1−w1+1−m2z^w2+1−mx^2x^,b13=−1−mx^w2+1−mx^,b21=b23=0,b22=w1x^−w31−mz^w2+1−mx^−w4,b31=w2w51−mz^w2+1−mx^2,b32=w2w3w51−m+w51−m2w3−cx^w2+1−mx^2z^,b33=w7z^w8+z^2.Hence, the characteristic equation ofJE3 can be written as
(28)λ32−b11+b33λ3+b11b33−b13b31b22−λ3=0.Straightforward computations show that all the roots (eigenvalues) of the characteristic equation (28) have negative real parts, and then, the DFEP is locally asymptotically stable provided that the following conditions hold.
(29a)1−m2x^z^w2+1−mx^2+w7z^w8+z^2<x^,(29b)w7w8+z^2<w71−m2z^w8+z^2w2+1−mx^2+w2w51−m2w2+1−mx^3,(29c)w1x^<w31−mz^w2+1−mx^+w4.Now, the local stability conditions for the CEP are established in the following theorem.Theorem 2.
Suppose that the CEP exists; then, it is locally asymptotically stable in theℝ+3 if the following sufficient conditions hold.
(30a)1−m2x∗z∗+w31−m2y∗z∗w2+1−mx∗+c1−my∗2+w7z∗w8+z∗2<x∗,(30b)c1−m2z∗w2+1−mx∗+c1−my∗2<1+w1,(30c)w3≤c<w2w3+w31−mx∗1−mx∗orw31−my∗−w21−my∗<c≤w3,(30d)a13a31a33<a11<mina12a21a22,a12a31a32,(30e)a22<a21a32a31,(30f)M1>M2,where all symbols are clearly described in the proof.Proof.
According to Equation (17), the Jacobian matrix at E4 is given by
(31)JE4=aij3×3,with aij=uijx∗,y∗,z∗ for all i,j=1,2,3. Therefore, the characteristic equation of JE4 can be written as
(32)λ43+D1λ42+D2λ4+D3=0,where
(33)D1=−a11+a22+a33,D2=a11a22−a12a21+a11a33−a13a31+a22a33−a23a32,D3=−a33a11a22−a12a21+a23a12a31−a11a32+a13a21a32−a22a31.
Moreover, it is easy to verify that(34)D1D2−D3=−a11+a22a11a22−a12a21−a11+a33a11a33−a13a31−a22+a33a22a33−a23a32−2a11a22a33+a12a23a31+a13a21a32=M1−M2,where
(35)M1=−a11+a33a11a33−a13a31−2a11a22a33+a12a23a31>0,M2=a11+a22a11a22−a12a21+a22+a33a22a33−a23a32−a13a21a32>0.
Recall that according to Routh-Hurwitz criterion, Equation (32) has roots (eigenvalues) with negative real parts if and only if D1>0,D3>0 and D1D2−D3>0. Now, straightforward computation shows that all the requirements of the Routh-Hurwitz criterion are satisfied provided that the conditions (30a), (30b), (30c), (30e), and (30f) hold. Hence, the CEP is locally asymptotically stable and the proof is complete.
It is well known that the persistence of an ecological system is an important subject for both environment and ecologists. The coexistence of all the species along the entire time life is known as persistence. Mathematically, this is equivalent to unapproachable of the trajectory to any of the boundary planes of axis. It is well known that system (3) has two subsystems. The first subsystem exists in the xy-plane, in case in the absence of the predator, and can be represented in system (36). However, the second subsystem exists in the xz-plane, in the absence of disease, and can be represented in system (37).
(36)dxdt=x1−x−y−w1y=g1x,y,dydt=yw1x−w4=g2x,y,(37)dxdt=x1−x−1−mzw2+1−mx=g3x,z,dzdt=zw51−mxw2+1−mx−w6−w7w8+z=g4x,z.
Recall that according to Bendixson-Dulac Criteria [26], the planar system X′=GX, where X=x1,x2T and G=g,fT has no periodic dynamics provided that there is a C1 function φx1,x2 (known as Dulac function) so that the expression ∂φg/∂x+∂φf/∂y has the same sign and does not equal zero almost everywhere in a simply connected region of the plane.
Consequently, by choosing the Dulac function1/xy, it is easy to verify that system (36) has no periodic dynamics in the interior of xy-plane, while with the Dulac function 1/xz, system (37) has no periodic dynamics fall in the interior of xz-plane provided that the following condition holds.
(38)1−m2w2+1−mx2+w7xw8+z2≠1z.
Now, in order to established the persistence conditions for system (3), the average Lyapunov function will be used as shown in the following theorem.Theorem 3.
Assuming that there is no periodic dynamics in the boundary planes; then, system (3) is persistent provided that the following conditions hold.
(39a)w1>w4,(39b)w51−mw2+1−m>w6+w7w8,(39c)w51−mx¯+w3y¯w2+1−mx¯+c1−my¯>w6+w7w8,(39d)w1x^>w31−mz^w2+1−mx^+w4.Proof.
Consider the following functionPx,y,z=xq1yq2zq3, where qi,∀i=1,2,3 are positive constants to be determined. Then, we obtain
(40)Λ=P′P=q11−x−y−w1y−1−mzw2+1−mx+c1−my+q2w1x−w31−mzw2+1−mx+c1−my−w4+q3w51−mx+w3yw2+1−mx+c1−my−w6−w7w8+z.
Now, according to the average Lyapunov function, the proof is done if we could prove that this functionΛ is positive at each boundary equilibrium point. Note that at E0, it obtains
(41)ΛE0=q1−w4q2−w6+w7w8q3.
So, for any positive constantsqi,∀i=1,2,3 with q1 sufficiently large compared with q2 and q3, we get ΛE0>0. At the DPFEP that is represented by E1=1,0,0, we have
(42)ΛE1=w1−w4q2+w51−mw2+1−m−w6+w7w8q3.
Clearly,ΛE1>0 under the conditions (39a) and (39b). Now, at the PFEP that is denoted by E2=x¯,y¯,0, it obtains
(43)ΛE2=q3w51−mx¯+w3y¯w2+1−mx¯+c1−my¯−w6+w7w8.
Here,ΛE2>0 under condition (39c). Finally, at the DFEP that is given by E3=x^,0,z^, we have
(44)ΛE3=q2w1x^−w31−mz^w2+1−mx^−w4.
Clearly, condition (39d) guarantees that ΛE3>0. Hence, the system has no omega limit sets on the boundary planes and thus, the system is persists.
Now, the global dynamics of system (3) is investigated using suitable Lyapunov functions as shown in the following theorem.Theorem 4.
Assuming that the DPFEP which is given byE1=1,0,0 is locally asymptotically stable, then, it is a globally asymptotically stable provided that
(45a)1+w1<w4,(45b)1−mw2<w6.Proof.
Consider the positive definite function(46)L1=x−1−lnx+y+z.
It is clear thatL1:ℝ+3⟶ℝ and L11,0,0=0 with L1x,y,z>0, for all x,y,z≠1,0,0. Now the derivative of L1 with respect to time can be written, after some algebraic computation, as
(47)dL1dt<−x−12−w4−1+w1y−w6−1−mw2z.
Therefore, by using conditions (45a) and (45b), it obtains dL1/dt which is negatively defined in ℝ+3. Thus, according to Lyapunov’s second theorem, E1 is globally asymptotically stable.Theorem 5.
Assuming that the PFEP that is given by is locally asymptotically stable, then, it is globally asymptotically stable inℝ+3 provided that
(48)1+w1w1w31−my¯w2+1−mx¯w2<w6w5.Proof.
Consider the positive definite function(49)L2=C1x−x¯−x¯lnxx¯+C2y−y¯−y¯lnyy¯+C2z.
It is clear thatL2:ℝ+3⟶ℝ and L2x¯,y¯,0=0 with L2x,y,z>0, for all x,y,z≠x¯,y¯,0=w4/w1,w1−w4/w11+w1,0. Then, the derivative of L2 with respect to time can be written, after some algebraic steps as
(50)dL2dt=−C1x−x¯2−C11+w1−C2w1x−x¯y−y¯−C1−C2w51−mxzw2+1−mx+c1−my−C2−C3w5w31−myzw2+1−mx+c1−my−C3w6−C2w31−my¯w2−C11−mx¯w2z−C3w7zw8+z.
So by choosing the positive constants asC1=1, C2=1+w1/w1, and C3=1/w5, the following is obtained
(51)dL2dt=−x−x¯2−1w1w31−myzw2+1−mx+c1−my−w6w5−1+w1w1w31−my¯w2−1−mx¯w2z−1w5w7zw8+z,dL2dt=−x−x¯2−1w1w31−myzw2+1−mx+c1−my.
Clearly, under the condition (48), dL2/dt is negatively semidefinite, because dL2/dt=0 if and only if x=x¯ and z=0 with any value of y. Therefore, according to Lyapunov’s second theorem, E2 is a stable point.
Now, form the first equation of system (3), we have dx/dt=x¯1−x¯−1+w1y, which is equal to zero if and only if y=y¯. Hence, the set E2 is the largest invariant set for which dL2/dt=0. Accordingly, the global asymptotic stability of the PFEP follows from LaSalle’s invariance principle. This completes the proof.Theorem 6.
Assuming that the DFEP that is given byE3=x^,0,z^ is locally asymptotically stable, then, it is globally asymptotically stable in ℝ+3 provided that the following sufficient conditions hold
(52a)1−m2z^w2N^1<1,(52b)c1−m2w4z^+1+w1w22N^1x^<1+w1w2w4N^1,(52c)p132<4p11p33,(52d)2p33z−z^2<p11x−x^+p33z−z^2,where all symbols are given in the proof.Proof.
Consider the positive definite function(53)L3=ρ1x−x^−x^lnxx^+ρ2y+ρ3z−z^−z^lnzz^.
It is clear thatL3:ℝ+3⟶ℝ and L3x^,0,z^=0 with L3x,y,z>0, for all x,y,z≠x^,0,z^. Then, the derivative of L3 with respect to time can be written, after some algebraic steps as
(54)dL3dt≤−ρ11−1−m2z^w2N^1x−x^2−1−mN1N^1ρ1N^1−ρ3w2w5x−x^z−z^+ρ3w7w8N^2z−z^2−ρ11+w1w2N^1−c1−m2z^w2N^1−ρ2w2xy−ρ2w4−ρ11+w1x^y−1−m2N1N^1ρ1c−ρ3w5x^yz^−1−mN1ρ2w3+ρ3w51−mx^N^1yz.where N1=w2+1−mx+c1−my, N^1=w2+1−mx^, N2=w8+z, and N^2=w8+z^. So, by choosing the positive constants as ρ1=1, ρ2=1+w1w2N^1−c1−m2z^/w22N^1, which is positive under (52a), and ρ3=c/w5. Therefore, we obtain
(55)dL3dt<−1−1−m2z^w2N^1x−x^2−1−mN1N^1N^1−cw2x−x^z−z^−cw7w5w8N^2z−z^2+2cw7w5w8N^2z−z^2−1+w1w2N^1−c1−m2z^w22N^1w4−1+w1x^y.
Therefore, according to the sufficient conditions (52a)–(52c), it is easy to obtain
(56)dL3dt<−p11x−x^+p33z−z^2+2p33z−z^2−1+w1w2N^1−c1−m2z^w22N^1w4−1+w1x^y.
Here,p11=1−1−m2z^/w2N^1, p13=1−m/N1N^1N^1−cw2, and p33=cw7/w5w8N^2.Theorem 7.
Assuming that the positive equilibrium pointE4=x∗,y∗,z∗ is locally asymptotically stable, then, all the points that satisfy the following conditions belong to the basin of attraction of E4.
(57a)1−m2z∗w2N1∗<1,(57b)q122<2q11q22,(57c)q132<2q11q33,(57d)q232<16q22q33,(57e)2q22y−y∗+q33z−z∗2<12q11x−x∗+2q22y−y∗2+12q11x−x∗+2q33z−z∗2.Proof.
Consider the positive definite function(58)L4=x−x∗−x∗lnxx∗+y−y∗−y∗lnyy∗+z−z∗−z∗lnzz∗.
It is clear thatL4:ℝ+3⟶ℝ and L4x∗,y∗,z∗=0 with L4x,y,z>0, for all x,y,z≠x∗,y∗,z∗ with x>0,y>0 and z>0. Then, the derivative of L4 with respect to time can be written, after some algebraic steps as
(59)dL4dt=x−x∗−x−x∗−1+w1y−y∗−1−mzw2+1−mx+c1−my+1−mz∗w2+1−mx∗+c1−my∗+y−y∗w1x−x∗−w31−mzw2+1−mx+c1−my+w31−mz∗w2+1−mx∗+c1−my∗+z−z∗−w7w8+z+w7w8+z∗+w51−mx+w3yw2+1−mx+c1−my−w51−mx∗+w3y∗w2+1−mx∗+c1−my∗.
Therefore, we obtain(60)dL4dt≤−q112x−x∗2−q12x−x∗y−y∗−q22y−y∗2−q112x−x∗2−q13x−x∗z−z∗−q33z−z∗2+2q22y−y∗2−q23y−y∗z−z∗+2q33z−z∗2,where
(61)q11=1−1−m2z∗w2N1∗,q12=1−c1−m2z∗w2N1∗−w31−m2z∗w2N1∗,q13=1−mN11−w2w5N1∗−cw51−my∗N1∗+w3w51−my∗N1∗,q23=1−mN1w3+cw51−mx∗N1∗−w2w3w5N1∗−w3w51−mx∗N1∗,q22=cw31−m2z∗w2N1∗,q33=w7w8N2∗,with N1=w2+1−mx+c1−my, N1∗=w2+1−mx∗+c1−my∗, and N2∗=w8+z∗.
Clearly according to the condition (57a), then, q11>0, while q22 and q33 are always positive. Therefore, by using the conditions (57b)–(57d), we obtain
(62)dL4dt≤−12q11x−x∗+2q22y−y∗2−12q11x−x∗+2q33z−z∗2+2q22y−y∗+q33z−z∗2.
The first tow terms are clearly negative while the third term is clearly positive; therefore, the sufficient condition (57e) guarantees the negative definite of the derivative. Hence, the solution starting with any point satisfying conditions (57a)–(57e) approaches asymptotically to E4. So these points are belonging to the basin of attraction.


## 4. Bifurcation Analysis

Now, in order to compute the second derivative of the Jacobian matrix, system (3) is rewritten in the vector form as follows:dX/dt=FX, X=x,y,zT,F=F1,F2,F3T=xf1,yf2,zf3T.(50)Therefore, according to the Jacobian matrix of system (3) at the point x,y,z that is given by Equation (17), we obtain
(63)D2Fx,y,zV,V=mij3×1,where V=v1,v2,v3T be any vector and
(64)m11=−21−1−m2w2+c1−myzN13v12−21−mw2+c1−myN12v1v3−21+w1−cz1−m2w2−1−mx+c1−myN13v1v2−2c21−m3xzN13v22+2c1−m2xN12v2v3,m21=−2w31−m3yzN13v12+2w31−m2yN12v1v3+2w1+w3z1−m2w2+1−mx−c1−myN13v1v2+2cw3z1−m2w2+1−mxN13v22−2w31−mw2+1−mxN12v2v3,m31=−2w51−m2zw2+1−mc−w3yN13v12−2w51−m2zw2c+w3−c−w31−mx+cc−w31−myN13v1v2−2cw51−m2zw2w3−1−mc−w3xN13v22+2w2w51−m+w51−m2c−w3yN12v1v3+2w2w3w51−m−w51−m2c−w3xN12v2v3+2w7w8N23v32.Consequently, the conditions guarantee the occurrence of local bifurcations around the nonhyperbolic equilibrium points that are established in the following theorems.Theorem 8.
Assume that condition (22b) holds. Letting the parameter w1 pass through the value w1∗≡w4, then, system (3) undergoes a transcritical bifurcation but neither saddle node nor pitchfork bifurcation can occurProof.
Note that when the parameterw1 passes through the positive value w1∗≡w4, it is easy to verify that the equilibrium point E1 becomes a nonhyperbolic point and the Jacobian matrix of system (3) at E1,w1∗ can be written as
(65)J1=JE1,w1∗=−1−1−w4−1−mw2+1−m00000w51−mw2+1−m−w6−w7w8.
Thus,J1 has the following eigenvalues λ111=−1, λ121=0, and λ131=w51−m/w2+1−m−w6−w7/w8. Clearly, λ131 is negative under the condition (22b).
LetV1=v11,v12,v13T be the eigenvector of J1 corresponding to the eigenvalue λ121=0. Then, straightforward computation gives V1=β1v12,v12,0T, where β1=−1+w4; v12 represents any nonzero real number.
Also, letΨ1=ψ11,ψ12,ψ13T represent the eigenvector of J1T that corresponds to the eigenvalue λ121=0. Then, direct calculation shows that Ψ1=0,ψ12,0T, where ψ12 is any nonzero real number.
Since the partial derivative of vector fieldF with respect to the parameter w1 is given by ∂F/∂w1=Fw1=−xy,xy,0T, hence, we obtain Fw1E1,w1∗=0,0,0T. Therefore,
(66)Ψ1TFw1E1,w1∗=0.
Thus system (3) at E1 with w1=w1∗ does not experience saddle-node bifurcation in view of the Sotomayor theorem [26]. Now, we have
(67)Ψ1TDFw1E1,w1∗V1=v12ψ12≠0,where DFw1 represents the derivative of Fw1 with respect to X. In addition, by substituting the value of E1, w1∗, V1, and Ψ1 in Equation (63), then we get
(68)Ψ1TD2FE1,w1∗V1,V1=2ψ12v122β1w4≠0.
Accordingly, by the Sotomayor theorem, system (3) near the equilibrium point E1 with w1=w1∗ possesses a transcritical bifurcation, while pitchfork bifurcation cannot occur and hence, the proof is complete.Theorem 9.
Assume that the parameterw6 passes through the value w6∗, where
(69)w6∗≡w51−mx¯+w3y¯w2+1−mx¯+c1−my¯−w7w8.
Then, system (3) near the PFEP undergoes a transcritical bifurcation but neither saddle-node bifurcation nor pitchfork bifurcation can occur provided that
(70)w2+1−mc−w3y¯Γ1+w2w3−1−mc−w3x¯Γ2+w7N¯12w51−mw82≠0,where Γ1,Γ2 are given in the proof.Proof.
Note that, straightforward computation shows that as the parameterw6 passes through the value w6∗, the PFEP becomes a nonhyperbolic equilibrium point and the Jacobian matrix of system (3) at E2,w6∗ can be written as
(71)J2=JE2,w6∗=−x¯−1+w1x¯−1−mx¯N¯1w1y¯0−w31−my¯N¯1000where N¯1=w2+1−mx¯+c1−my¯. Hence, J2 has three eigenvalues given by λ212=−x¯/2−x¯2−4w11+w1x¯y¯/2, λ222=−x¯/2+x¯2−4w11+w1x¯y¯/2, and λ232=0.
Obviously, the first two eigenvalues,λ212 and λ222, have a negative real part.
Now, letV2=v21,v22,v23T be the eigenvector of J2 corresponding to the eigenvalue λ232=0. Then, straightforward computation gives V2=Γ1v23,Γ2v23,v23T, where Γ1=w31−m/w1N¯1>0,Γ2=−w1+w31−m/w11+w1N¯1<0,and v23 represents any nonzero real number.
Also, letΨ2=ψ21,ψ22,ψ23T represent the eigenvector of J2T that corresponds to the eigenvalue λ232=0. Then, direct calculation shows that Ψ2=0,0,ψ23T, where ψ23 is any nonzero real number.
Since the partial derivative of vector fieldF with respect to the parameter w6 is given by ∂F/∂w6=Fw6=0,0,−zT, hence, we obtain Fw6E2,w6∗=0,0,0T. Therefore,
(72)Ψ2TFw6E2,w6∗=0.
Thus, system (3) at E2 with w6=w6∗ does not experience saddle-node bifurcation in view of the Sotomayor theorem. Now, we have
(73)Ψ2TDFw6E2,w6∗V2=−v23ψ23≠0,where DFw6 represents the derivative of Fw6 with respect to X. In addition, by substituting the value of E2, w6∗, V2, and Ψ2 in Equation (63), then we get
(74)Ψ2TD2FE2,w6∗V2,V2=2w51−mψ23v232N¯12w2+1−mc−w3y¯Γ+w2w3−1−mc−w3x¯Γ2+w7w51−mw821.
Obviously,Ψ2TD2FE2,w6∗V2,V2≠0 under the condition (70), and hence by the Sotomayor theorem, system (3) near the equilibrium point E2 with w6=w6∗ possesses a transcritical bifurcation, while pitchfork bifurcation cannot occur and hence, the proof is complete.Theorem 10.
Assume that conditions (29a) and (29b) hold. Let the parameter w4 pass through the value w4∗=w1x^−w31−mz^/w2+1−mx^, then, system (3) near the DFEP undergoes a transcritical bifurcation but neither saddle-node bifurcation nor pitchfork bifurcation can occur provided that
(75)w1N^12+w31−m2z^ω1+cw31−m2z^−w31−mN^1ω2≠0,where ω1,ω2 are given in the proof.Proof.
Straightforward computation shows that as the parameterw4 passes through the value w4∗ then, the DFEP becomes a nonhyperbolic equilibrium point and the Jacobian matrix of system (3) at E3,w4∗ can be written as
(76)J3=JE3,w4∗=b^ij3×3,where b^ij=bij for all i,j=1,2,3 are given in Equation (28) with b^22w4∗=0. Hence, J3 has three eigenvalues given by
(77)λ313=−b11+b332−b11+b332−4b11b33−b13b312,λ323=0,λ333=−b11+b332+b11+b332−4b11b33−b13b312.
Clearly, the eigenvaluesλ313 and λ333 have negative real parts under the conditions (29a) and (29b).
Now, letV3=v31,v32,v33T be the eigenvector of J3 corresponding to the eigenvalue λ323=0. Then, straightforward computation gives V3=ω1v32,v32,ω2v32T, where ω1=b13b32−b12b33/b11b33−b13b31andω2=b12b31−b32b11/b11b33−b13b31, while v32 represents any nonzero real number.
Also, letΨ3=ψ31,ψ32,ψ33T represent the eigenvector of J3T that corresponds to the eigenvalue λ323=0. Then, direct calculation shows that Ψ3=0,ψ32,0T, where ψ32 is any nonzero real number.
Since the partial derivative of vector fieldF with respect to the parameter w4 is given by ∂F/∂w4=Fw4=0,−y,0T, hence, we obtain Fw4E3,w4∗=0,0,0T. Therefore,
(78)Ψ3TFw4E3,w4∗=0.
Thus, system (3) at E3 with w4=w4∗ does not experience saddle-node bifurcation in view of the Sotomayor theorem. Now, we have
(79)Ψ3TDFw4E3,w4∗V3=−v32ψ32≠0where DFw4 represents the derivative of Fw4 with respect to X. In addition, by substituting the value of E3, w4∗, V3, and Ψ3 in Equation (63), then, we get
(80)Ψ3TD2FE3,w4∗V3,V3=2ψ32v322N^12w1N^12+w31−m2z^ω1+cw31−m2z^−w31−mN^1ω2.
Obviously,Ψ3TD2FE3,w4∗V3,V3≠0 under the condition (75), and hence, by the Sotomayor theorem, system (3) near the DFEP that is given by E3 with w4=w4∗ possesses a transcritical bifurcation, while pitchfork bifurcation cannot occur and hence, the proof is complete.Theorem 11.
Assume that in addition to the conditions (30a), (30b), (30c), and (30e) the following conditions are satisfied
(81)a12a31a32<a11<a12a21a22.
Then, as the parameterw3 passes through the value w3∗, system (3) near the CEP undergoes a saddle-node bifurcation, but neither transcritical bifurcation nor pitchfork bifurcation can occur provided that the following conditions hold
(82a)w5−ξ4≠0,(82b)ξ3m114+ξ4m214+m314≠0,where all the symbols are clearly described in the proof, while
(83)w3∗≡−N1∗1−my∗a33a11a22−a12a21+a13a21a32−a31a22a11a32−a12a31,with aij,∀i,j=1,2,3 are given in Equation (31) and N1∗=w2+1−mx∗+c1−my∗.Proof.
Straightforward computation shows that the above given conditions guarantee that the parameterw3∗ is positive and the determinant (D3) of the Jacobian matrix J4=JE4,w3∗=aijw3∗3×3 that is given in Equation (32) is zero.
Therefore, the Jacobian matrixJ4 has the following eigenvalues
(84)λ414=−D12+12D12−4D2,λ424=−D12−12D12−4D2,λ434=0.
It is easy to verify thatλ414 and λ424 have negative real parts under the given conditions while the third eigenvalue is zero. So the CEP becomes a nonhyperbolic point.
LetV4=v41,v42,v43T be the eigenvector of J4 corresponding to the zero eigenvalue λ434=0. Then, direct computation shows that V4=ξ1v43,ξ2v43,v43T, where v43 represents any nonzero real number and ξ1=a12a23−a22a13/a11a22−a12a21<0 while ξ2=a21a13−a23a11/a11a22−a12a21>0.
LetΨ4=ψ41,ψ42,ψ43T represent the eigenvector of J4T that is corresponding to the eigenvalue λ434=0. Then, simple calculation shows that Ψ4=ξ3ψ43,ξ4ψ43,ψ43T, where ψ43 is any nonzero real number and ξ3=a21a32−a22a31/a11a22−a12a21<0 while ξ4=a12a31−a11a32/a11a22−a12a21>0.
Now, we have∂F/∂w3=Fw3=0,−1−myz/N1,w51−myz/N1T, where N1=w2+1−mx+c1−my; hence, we obtain Fw3E4,w3∗=0,−1−my∗z∗/N1∗,w51−my∗z∗/N1∗T.
Therefore, by using condition (82a) we get
(85)Ψ4TFw3E4,w3∗=ψ431−my∗z∗N1∗w5−ξ4≠0.
Consequently, the first condition of saddle-node bifurcation in view of the Sotomayor theorem is satisfied. Now, since(86)Ψ4TD2FE4,w3∗V4,V4=ξ3m114+ξ4m214+m314ψ43,where mi14=mi1E4,w3∗,∀i=1,2,3; then, according to the condition (82b), we have Ψ4TD2FE4,w3∗V4,V4≠0. Hence, the system (3) undergoes a saddle-node bifurcation near the coexistence equilibrium point but neither transcritical bifurcation nor pitchfork bifurcation can occur. This completes the proof.


## 5. Numerical Simulation

In this section, system (3) is simulated numerically. The objectives understand the dynamical behavior of the solution of the system as the time increasing with the parameters varied and confirm our obtained analytical results. A variety of mathematical tools is used in this investigation including the bifurcation diagrams, chaotic attractors, phase portrait, and time series. Four steps of the predictor-corrector method are used to solve system (3), and then, the obtained trajectories are drawn using MATLAB version 6. Now, the following set of hypothetical parameters are used to solve system (3) numerically.
(87)w1=0.5,m=0.1,w2=0.2,c=0.55,w3=0.4w4=0.1,w5=0.8,w6=0.1,w7=0.4,W8=0.8In order to investigate the behavior of the system (3) as a function of c, bifurcation diagram is drawn in Figure 1.Figure 1
Bifurcation diagram of system (3) using data (87) with varying c, the trajectory of z versus c in the range 0,0.6.According to the bifurcation diagram given in Figure1, system (3) has a complex dynamics including chaos. Moreover, the system enters to chaos through cascades of periodic doubling and then, after a specific value (c≥0.58), the system loses the persistence, and the trajectory approaches asymptotically to PFEP in the xy-plane as shown in Figure 2 for typical values of c.Figure 2
Trajectory of system (3) using the data given by Equation (87). (a) Period 3 in Int.ℝ+3 where c=0.45. (b) Time series of attractor in (a). (c) Period 5 in Int.ℝ+3 where c=0.47. (d) Time series of attractor in (c). (e) Chaotic attractor in Int.ℝ+3 where c=0.55. (f) Time series of attractor in (e). (g) Asymptotic stable PFEP in the xy-plane given by E2=0.2,0.53,0 where c=0.6.
(a)(b)(c)(d)(e)(f)(g)(h)The dynamical behavior of system (3) as a function of refuge rate m in the range 0,0.15 is investigated with the help of bifurcation diagram as presented in Figure 3.Figure 3
Bifurcation diagram of system (3) using data (87) with varying m, the trajectory of z versus m in the range 0,0.15.Clearly, system (3) has chaotic dynamics where in between there are small periodic regions; moreover, the system approaches asymptotically to PFEP in the xy-plane as m increases above a specific value (m≥0.12), see Figure 4 for typical values of m.Figure 4
Trajectory of system (3) using the data given by Equation (87). (a) Period 3 in Int.ℝ+3 where m=0.024. (b) Time series of attractor in (a). (c) Chaotic attractor in Int.ℝ+3 where m=0.026. (d) Time series of attractor in (c). (e) Chaotic attractor in Int.ℝ+3 where m=0.06. (f) Time series of attractor in (e). (g) Asymptotic stable PFEP in the xy-plane given by E2=0.2,0.53,0 where m=0.15.
(a)(b)(c)(d)(e)(f)(g)(h)The dynamics of system (3) as a function of w1 in the range (0.3≤w1≤0.55) is shown in the bifurcation diagram given in Figure 5 and the typical attractors given in Figure 6.Figure 5
Bifurcation diagram of system (3) using data (87) with varying w1. (a) The trajectory of y versus w1 in the range 0.3,0.55. (b) The trajectory of z versus w1 in the range 0.3,0.55.
(a)(b)Figure 6
Trajectory of system (3) using the data given by Equation (87). (a) Periodic attractor in the xy-plane where w1=0.3. (b) Time series of attractor in (a). (c) Periodic attractor in Int.ℝ+3 where w1=0.4. (d) Time series of attractor in (c). (e) Chaotic attractor in Int.ℝ+3 where w1=0.495. (f) Time series of attractor in (e). (g) Asymptotic stable PFEP in the xy-plane given by E2=0.16,0.52,0 where w1=0.6.
(a)(b)(c)(d)(e)(f)(g)(h)Further investigation for the effect of varyingw2 in the range 0.15,0.21 on the dynamical behavior of system (3) is carried out using the bifurcation diagram and phase portraits with their time series at typical values of w2 as shown in Figures 7 and 8, respectively.Figure 7
Bifurcation diagram of system (3) using data (87) with varying w2, the trajectory of z versus w2 in the range 0.15,0.21.Figure 8
Trajectory of system (3) using the data given by Equation (87). (a) Asymptotic stable PFEP in the xy-plane given by E2=0.2,0.53,0 where w2=0.15. (b) Time series of attractor in (a). (c) Chaotic attractor in Int.ℝ+3 where w2=0.19. (d) Time series of attractor in (c).
(a)(b)(c)(d)According to the above figures, system (3) has a range of complex dynamics including chaos given by 0.16,0.21, while it approaches asymptotically to PFEP otherwise.Other parameters have been investigated too, and the obtained results are summarized in Table1.Table 1
The dynamical behavior of system (3) using the data (87) with varying a specific parameter.
The parameterThe rangeThe dynamical behavior of system (3)w30<w3<0.39Approaches to periodic dynamic inxy-plane0.39≤w3<0.42Approaches to chaotic attractor inInt.ℝ+30.42≤w3<0.64Approaches to periodic attractor inInt.ℝ+30.64≤w3<1Approaches to periodic dynamic inxz-planew40<w4<0.1Approaches to PFEP inxy-plane0.1≤w4<0.11Approaches to chaotic attractor inInt.ℝ+30.11≤w4<0.18Approaches to periodic attractor inInt.ℝ+30.18≤w4<1Approaches to periodic dynamic inxz-planew50<w5<0.8Approaches to PFEP inxy-plane0.8≤w5<0.81Approaches to chaotic attractor inInt.ℝ+30.81≤w5<1Approaches to periodic attractor inInt.ℝ+3w60<w6<0.1Approaches to periodic attractor inInt.ℝ+30.1≤w6<0.11Approaches to chaotic attractor inInt.ℝ+30.11≤w6<1Approaches to PFEP inxy-planew70<w7<0.34Approaches to periodic dynamic inxz-plane0.34≤w7<0.4Approaches to periodic attractor inInt.ℝ+30.4≤w7<0.41Approaches to chaotic attractor inInt.ℝ+30.41≤w7<1Approaches to PFEP inxy-planew80<w8<0.79Approaches to PFEP inxy-plane0.79≤w8<0.81Approaches to chaotic attractor inInt.ℝ+30.81≤w8<1Approaches to periodic attractor inInt.ℝ+3


## 6. Discussion and Conclusion

In this paper, an ecoepidemiological model consisting of the prey-predator system with infectious disease in prey has been proposed and analyzed. As a defensive property against predation, it is assumed that the prey has a refuge and the predator preys upon the available prey individuals only according to the modified Holling type-II functional response. There is a harvesting process applied on the predator species only with a nonlinear harvesting rate. The existence, uniqueness, and boundedness of the solution of the system are studied. The stability analyses around all possible equilibrium points are investigated. The persistence conditions of the system are derived analytically and shown numerically. The occurrences of local bifurcation near the nonhyperbolic equilibrium points are investigated. Finally, numerical simulations are carried out to confirm our obtained analytical results and specify the control set of parameters using biologically feasible set of parameters given in (87). It is observed that the system is very sensitive to variation in the parameter set and it has a complex dynamics for a wide range of parameters including chaos. Different types of attracting sets and bifurcation diagrams are drawn to present the behavior of the system. Finally, in the following, we summarize the obtained numerical results depending on the set of data given by Equation (87).
(1)
The system (3) is sensitive for varying in the preference rate of predation as shown in the bifurcation diagram represented by Figure 1, so that it has periodic dynamics in the interior of the positive octant and then periodic doubling leading to chaos, as shown in Figure 2, for the range c<0.58. Otherwise, the system loses the persistence and approach asymptotically to the PFEP in the interior of the xy-plane(2)
According to the bifurcation diagram as a function of the refuge rate given by Figure3 and the attracting sets drawn in Figure 4, system (3) has a chaotic dynamics for the range m<0.12 in between; there is windows of periodic dynamics too. Moreover, as the refuge rate increases further, the system loses their persistence and approach asymptotically to the PFEP in the interior of the xy-plane(3)
According to the bifurcation diagram as a function of the infected rate given by Figure5 and the attracting sets drawn in Figure 6, system (3) has a periodic dynamics in the interior of the xz-plane for the range w1<0.37 and then, it becomes persistent in the form of periodic dynamics and periodic doubling leading to chaos for the range w1∈0.37,0.53. Otherwise, system (3) loses their persistence and approach asymptotically to the PFEP in the interior of the xy-plane(4)
According to the bifurcation diagram as a function of half saturation constant given by Figure6 and the attracting sets drawn in Figure 7, system (3) has a chaotic dynamics and in between there is periodic windows for the range w2∈0.16,0.21. Otherwise, system (3) loses their persistence and approach asymptotically to the PFEP in the interior of the xy-plane(5)
According to the results given in Table1, system (3) is sensitive to varying in other parameters too and still persists in the form of chaos and periodic in the interior of the positive octant for some ranges of parameters



---

*Source: 5373817-2020-06-03.xml*