# The Dynamics of a Tritrophic Leslie-Gower Food-Web System with the Effect of Fear

Authors: Firas Hussean Maghool; Raid Kamel Naji
Journal: Journal of Applied Mathematics
 (2021)

---

## Abstract

The avoidance strategy of prey to predation and the predation strategy for predators are important topics in evolutionary biology. Both prey and predators adjust their behaviors in order to obtain the maximal benefits and to raise their biomass for each. Therefore, this paper is aimed at studying the impact of prey’s fear and group defense against predation on the dynamics of the food-web model. Consequently, in this paper, a mathematical model that describes a tritrophic Leslie-Gower food-web system is formulated. Sokol-Howell type of function response is adapted to describe the predation process due to the prey’s group defensive capability. The effects of fear due to the predation process are considered in the first two levels. It is assumed that the generalist predator grows logistically using the Leslie-Gower type of growth function. All the solution properties of the model are studied. Local dynamics behaviors are investigated. The basin of attraction for each equilibrium is determined using the Lyapunov function. The conditions of persistence of the model are specified. The study of local bifurcation in the model is done. Numerical simulations are implemented to show the obtained results. It is watched that the system is wealthy in its dynamics including chaos. The fear factor works as a stabilizing factor in the system up to a specific level; otherwise, it leads to the extinction of the predator. However, increasing the prey’s group defense leads to extinction in predator species.


---

## Body

## 1. Introduction

Mathematical modeling is used to understand the interaction of organisms with their surrounding environment, as well as species evolution that is well-known biological evolution. Nonlinear differential equation systems have always played an important role in simulating such mathematical models. Therefore, in the last few decades, there are increasing interests in the study of the dynamics of such systems. After the pioneering prey-predator model by Lotka-Volterra, the prey-predator models have been essential in mathematical ecology and have been extensively investigated. The Lotka-Volterra model was then adjusted by compiling a prey’s logistic growth and a Holling-type II functional response to describe the predation [1].It is completely recognized that a prey-predator is fundamental interaction for drawing the dynamics of food webs in the real-world environment, predation risks can negatively affect prey biomass and growth efficiency, and thus, predators affect the structure of natural communities. Prey can die of fear, and it can also reduce fertility, as a result of which fear can be just as important as outright killing by predators in affecting prey numbers [2–4]. Many prey-predator models have been suggested and investigated thoroughly in which either the predator kills the prey or the presence of the predator affects the behavior of the prey population due to the fear of the predation process [5]. Upadhyay and Mishra [6] modeled a prey-predator system taken into account the fear effect in prey growth. They watched that a relative rise in the scale of fear may minimize the size of the total population. Panday et al. [7] investigated the action of fear in a three-species food chain model, where the evolution rate of lower predator is minimized due to the fear’s action from a higher predator, and the evolution rate of prey is reduced due to the fear’s action from the lower predator. Das and Samanta [8] studied the fear’s impact from the predator on prey in case of providing additional food for the predator. Recently, Kumar and Kumari [9] discussed the action of fear on the dynamics of the food chain model having three species. They obtained that for low fear’s scale, the system stays chaotic while relatively high fear’s scale leads to stability. Moreover, the rise of fear’s scale further causes population extinction.On the other hand, Aziz-Alaoui [10] investigated a dynamic system consisting of a Leslie-Gower food chain of three species. He displayed that for a realistic parameter set, chaotic dynamics could be obtained. Zhang et al. [11] studied a harvested Leslie-Gower prey-predator model. They got with the help of Pontryagin’s maximal principle the optimal harvesting policy for the model. Cai et al. [12] investigated a Leslie-Gower prey-predator model with the diffusion and Allee effect on prey. They gained that the Allee effect basically rises the model spatiotemporal complexity. Later on, Abid et al. [13] were interested in the stability and optimal harvesting of a modified Leslie-Gowe prey-predator model with Holling-type II functional response. Recently, Singh and Bhadauria [14] studied the dynamics of the prey-predator model with weak Allee effect II and modified Leslie-Gower. They explained that Allee effect II greatly impacts the model and can raise the risk of extinction.Keeping the above in mind, in this paper, a tritrophic Leslie-Gower food-web system is constructed mathematically. Sokol-Howell type of function response [15] is used for describing the predation process due to the prey’s group’s defensive capability. The effects of fear due to the predation process are included and studied. The following section is written out as follows: the construction of the model is done. In Section 3, existence and local stability analysis are studied. In Section 4, persistence is discussed. However, in Section 5, the region of attraction for each equilibrium point of system (4) is determined by using the direct method of Lyapunov. Local bifurcation analysis is discussed in Section 6. In Section 7, the numerical simulation of the model is studied. Finally, Section 8 gives the conclusion.


## 2. The Mathematical Construction of the Model

In this section, a real-world tritrophic food-web system is constructed mathematically. It is supposed that the logistic law is applied for both the prey and the generalist-predator so that the environmental carrying capacity of the generalist-predator is not a constant but a function of the available food quantity so that there are upper limits to the rates of increase of both prey and generalist-predator. Therefore, the Leslie-Gower model is used in the construction of the proposed model, which assumes that both the prey and the generalist-predator grow according to the logistic law. Moreover, both the prey and specialist-predator have a defensive property against the predation when they are available in abundance. Therefore, a simplified Holling-type IV functional response or Sokol-Howell type of functional response is used in the proposed model. This is because the per capita predation rate in the Sokol-Howell functional response rises with prey density to a maximum at a critical prey density beyond which it decreases.LetXT represents the prey density at the time T, YT is the specialist-predator density at the time T, while ZT is the generalist-predator density at time T. Accordingly, the dynamics of such a food-web system can be described as follows:
(1)dXdT=XgX−YPX−Zq1X,Y,dYdT=a3YPX−Zq2X,Y−dY,dZdT=c1Z2−Z2c2+c3X+c4Y,where gX=r1−bX is the logistic growth rate of the prey species, while PX=a0X/1+m0X2 and q1X,Y=a1X/1+m1X2+m2Y2, with q2X,Y=a2X/1+m1X2+m2Y2, are the Sokol-Howell functional responses for the specialist-predator and generalist-predator, respectively. The parameters of system (1) are supposed to be positive, and they are described as the following: the parameter r>0 represents the intrinsic growth rate, 1/b>0 is the carrying capacity of the environment for the prey, the positive parameters a0,a1, and a2 are the maximum attack rates, the positive parameters m0,m1, and m2 are the preference rates of the predator to their prey species, a3>0 represents the conversion rate from individuals of the prey to specialist-predator, and c1>0 is the intrinsic growth rate of a generalist-predator, while the positive parameters c3andc4 are the respective food preference of generalist-predator from the prey and specialist-predator, respectively. In addition to the above, the square term c1Z2 in the third equation that represents the modified Leslie-Gower term signifies the fact that mating frequency is directly proportional to the number of males as well as to that of females, see [16, 17]. Finally, the parameter c2>0 is a constant added in the denominator to normalize the residual reduction in the generalist-predator due to a heavy lack of the preferable food X and Y, while d>0 is the natural death rate of the specialist-predator. Now, to add the factor of the fear for each of the prey and the specialist-predator from the predation by an upper-level predator, the growth of each of them should be reduced by a rate proportional with the abundance of their predators. Therefore, the dynamics of the food-web system given by system (1), in case of the existence of fear, can be represented as
(2)dXdT=rX1+n1Y+n2z−bX2−a0XY1+m0X2−a1XZ1+m1X2+m2Y2,dYdT=a3XY1+m0X211+n3Z−a2YZ1+m1X2+m2Y2−dY,dZdT=c1Z2−Z2c2+c3X+c4Y,where the positive parameters n1,n2, and n3 represent the fear coefficients of the prey and specialist-predator from their direct predator. Clearly, system (2) has 16 parameters, which leads to analysis difficulty. So, in order to simplify the analysis of system (2) and minimize the number of parameters, the following nondimensional variables and parameters are utilized.
(3)x=bXr,y=aoYr,z=aoZr,t=rT,α1=n1rao,α2=n2ra1,α3=mor2b2,α4=m1r2b2,α5=m2r2ao2,α6=a3b,α7=n3ra1,α8=aoa1,α9=dr,α10=ca1,α11=c2,α12=c3rb,α13=c4rao.Therefore, the nondimensional system corresponding to system (2) is given by
(4)dxdt=x11+α1y+α2z−x−y1+α3x2−z1+α4x2+α5y2=xf1,dydt=yα6x1+α3x211+α7z−α8z1+α4x2+α5y2−α9=yf2,dzdt=zα10z−α10zα11+α12x+α13y=zf3.Here, the interaction functions are defined onℝ+3=x,y,z:xt≥0,yt≥0,zt≥0. Moreover, since the right-hand side functions of system (4) are continuous and have continuous partial derivatives, hence they are Lipschitz functions. Thus, the solution of system (4) exists and is unique.Theorem 1.
The solutions of system (4), which start in ℝ+3 are uniformly bounded.Proof.
Letxt,yt,zt be a solution of system (4) with nonnegative initial condition x0,y0,z0∈ℝ+3. From the first equation in system (4), it is obtained that
(5)dxdt≤x−x2,hence, by solving this differential inequality, it is obtained that x≤1 as t⟶∞. Define the function w1=x+y/α6; then, the following is resulting
(6)dw1dt≤14+α9−α9x+yα6,or dw1/dt+α9w1≤1/4+α9.
Then, solving this differential inequality gives fort⟶∞:
(7)w1=x+yα6≤14α9+1.
Assume thatw2=x+y/α6+z, then by differentiating it with respect to time gives
(8)dw2dt≤14+α9−α9x+yα6+z+α9z+α10Z2−α10θZ2,where θ=α11+α12+α13α6/4α9+α6.
Therefore, rearranging and substituting the upper bound of the logistic term give(9)dw2dt+α9w2≤14+α9+β1,where β1=α92θ/4α101−θ.
Then, solving the last differential inequality ast⟶∞ gives
(10)w2≤14α9+1+β1α9=β2.
Hence, all solutions of system (4) that initiate in ℝ+3 are confined in the region Λ=x,y,z∈R+3,w2t≤β2, that is, all the solution are uniformly bounded.☐


## 3. Equilibria and Their Local Stability

System (4) has at most six equilibrium points, which can be described as follows:
(1)
The trivial equilibrium point (TP) denoted byE0=0,0,0 always exist(2)
The first axial equilibrium point (FAP) denoted byE1x=1,0,0 always exist(3)
The second axial equilibrium point (SAP) denoted byE1z=0,0,z^, where z^ is any positive real number, exists if and only if(11)α11=1(4)
The generalist-predator free equilibrium point (GPFP) that denoted byE2=x¯,y¯,0, where(12a)x¯=α6+α62−4α92α32α9α3,y¯=−σ1x¯+σ1x¯2−4α1σ0x¯2α1,with σ1x¯=α1x¯+α1α3x¯3+1>0, and σ0x¯=x¯−11+α3x¯2, exists uniquely in the positive quadrant of xy-plane if and only if the following conditions hold:
(12b)4α92α3≤α62,x¯<1(5)
The specialist-predator free equilibrium point (SPFP) denoted byE3=x̿,0,z̿ can be represented as(13)x̿=1−α11α12,z̿=−β1x̿+β1x̿2−4α2β0x̿2α2,where β1x̿=α2x̿+α2α4x̿3+1>0, and β0x̿=x̿−11+α4x̿2. Obviously, SPFP exists uniquely in the xz-plane if and only if the following condition holds:
(14)α11+α12>1(6)
The coexistence equilibrium point (CP) denoted byE4=x∗,y∗,z∗ can be represented as follows:(15)x∗=1−α11+α13y∗α12.The point y∗,z∗ represents the intersection point of the following two isoclines in the positive quadrant of the yz-plane.
(16a)g1y,z=11+α1y+α2z−x∗−y1+α3x∗2−z1+α4x∗2+α5y2=0,(16b)g2y,z=α6x∗1+α3x∗211+α7z−α8z1+α4x∗2+α5y2−α9=0.Simple computation displays that the two isoclines given by (16a) and (16b) have a unique positive intersection point, denoted by y∗,z∗, provided that the following sufficient conditions hold:
(17a)α11+α13y∗<1,(17b)α11+α12<1,(17c)α9α122+α9α31−α112<α12α61−α11,(17d)y1<y2,(17e)dydz=−∂g1/∂z∂g1/∂y>0,(17f)dydz=−∂g2/∂z∂g2/∂y<0.Here,y1 and y2 represent the positive root of the following two polynomials, respectively:
(18)ς1y4+ς2y3+ς3y2+ς4y+ς5=0,ς6y2+ς7y+ς8=0, where
(19)ς1=α1α3α133>0,ς2=α3α132α13−3α11−α11>0,ς3=α1α122α13−α12+3α13α31−α11α11−α11−α13+α3α132α12>0,ς4=α122α1α11+α13−α1+α31−α112α1α11+α13−α1−2α3α131−α11α12+α11−1−α123,ς5=α12−1−α11α122+1−α112<0,ς6=α9α3α122>0,ς7=α12α9α13−2α9α3α131−α11,ς8=α9α122+α9α31−α112+α11α12α6−α12α6<0.Now, the local stability of the above equilibrium points is studied by computing the Jacobian matrix (JM) and then determining their eigenvalues. It is simply to verify that the JM of system (4) at the point x,y,z can be represented as
(20)J=x∂f1∂x+f1x∂f1∂yx∂f1∂zy∂f2∂xy∂f2∂y+f2y∂f2∂zz2∂f3∂xz2∂f3∂yz2∂f3∂z+2zf3,where
(21)∂f1∂x=−1+2α3xyΓ22+2α4xzΓ12,∂f1∂y=−α1Γ32−1Γ2+2α5yzΓ12,∂f1∂z=−α2Γ32−1Γ1,∂f2∂x=α6−α6α3x2Γ22Γ4+2α8α4xzΓ12,∂f2∂y=2α5α8yzΓ12,∂f2∂z=−α6α7xΓ2Γ42−α8Γ1,∂f3∂x=α10α12Γ52,∂f3∂y=α10α13Γ52,∂f3∂z=0,with Γ1=1+α4x2+α5y2, Γ2=1+α3x2, Γ3=1+α1y+α2z, Γ4=1+α7z, and Γ5=α11+α12x+α13y.Consequently, for the TP, the JM becomes(22)JE0=1000−α90000.Then, the eigenvalues ofJE0 are given by λ10=1>0,λ20=−α9<0,andλ30=0. Since there is a positive equilibrium point, then E0 is an unstable nonhyperbolic point.For the FAP, the JM can be represented as(23)JE1x=−1000α61+α3−α90000.Then, the eigenvalues ofJE1x are clearly given by λ11=−1<0,λ21=α6/1+α3−α9,andλ31=0. Hence, the FAP is a nonhyperbolic point due to the existence of zero eigenvalue.For the SAP, the JM can be represented as(24)JE1z=11+α2z^−z^000−α8z^+α90α10α12z^2α112α10α13z^2α1120.Then, the eigenvalues ofJE1z are clearly given by λ21=1/1+α2z^−z^, λ22=−α8z^+α9<0, and λ23=0. Hence, the SAP is a nonhyperbolic point due to the existence of zero eigenvalue.For the GPFP, the JM can be represented as(25)JE2=aij3×3,where
(26)a11=x¯1−1+2α3x¯1y¯1Γ¯22,a12=x¯1−α1Γ¯32−1Γ¯2,a13=x¯1−α2Γ¯32−1Γ¯1,a21=α6y¯11−α3x¯12Γ¯22,a22=0,a23=x¯1−α6α7y¯1Γ¯22−α8Γ¯1,a31=0,a32=0,a33=0,with Γ¯1=1+α4x¯12+α5y¯12, Γ¯2=1+α3x¯12, and Γ¯3=1+α1y¯1. Clearly, the characteristic equation of JE2 can be represented as
(27)λλ2−a11λ−a12a21=0.Because the JM has a zero eigenvalue (λ=0), then E2 is a nonhyperbolic equilibrium point.For the SPEP, the JM can be determined as(28)JE3=bij3×3,where
(29)b11=x¯¯−1+2α4x¯¯z̿Γ̿12,b12=x¯¯−α1Γ̿32−1Γ̿2,b13=x¯¯−α2Γ̿32−1Γ̿1,b21=0,b22=α6x¯¯Γ̿12Γ̿4−α8z̿Γ̿1−α9,b23=0,b31=α10α12z̿2Γ̿52,b32=α10α13z̿2Γ̿52,b33=0,with Γ̿1=1+α4x¯¯2, Γ̿2=1+α3x¯¯2, Γ̿3=1+α2z̿, Γ̿4=1+α7z̿, and Γ̿5=α11+α12x¯¯. Therefore, the characteristic equation of JE3 can be represented as
(30)b22−λλ2−b11λ−b13b31=0.Clearly, the roots (eigenvalues) of equation (30) can be represented as
(31)λ13=b11+b112−4b13b312,λ23=α6x¯¯Γ̿12Γ̿4−α8z̿Γ̿1−α9,λ33=b11−b112−4b13b312.It is easy to verify that all the eigenvalues ofJE3 are negative, and hence, the SPFP is locally asymptotically stable (L.A.S), presuming that the next conditions hold:
(32a)2α4x¯¯z̿Γ̿12<1,(32b)α6x¯¯Γ̿12Γ̿4−α8z̿Γ̿1−α9<0.The JM of system (4) at the CP can be represented as
(33)JE4=cij3×3,where
(34)c11=x∗−1+2α3x∗y∗Γ2∗2+2α4x∗z∗Γ1∗2,c12=x∗−α1Γ3∗2−1Γ2∗+2α5y∗z∗Γ1∗2,c13=x∗−α2Γ3∗2−1Γ1∗,c21=y∗α61−α3x∗2Γ2∗2Γ4∗+2α8α4x∗z∗Γ1∗2,c22=2α5α8y∗2z∗Γ1∗2,c23=y∗−α6α7x∗Γ2∗Γ4∗2−α8Γ1∗,c31=α10α12z∗2Γ5∗2,c32=α10α13z∗2Γ5∗2,c33=0,with Γ1∗=1+α4x∗2+α5y∗2, Γ2∗=1+α3x∗2, Γ3∗=1+α1y∗+α2z∗, Γ4∗=1+α7z∗, and Γ5∗=α11+α12x∗+α13y∗. Consequently, the characteristic equation of JE4 can be represented as
(35)λ3+A1λ2+A2λ+A3=0,where A1=−c11+c22,
(36)A2=−c13c31−c11c22+c23c32+c12c21,A3=−c13c21c32−c22c31+c23c12c31−c11c32, with
(37)∆=A1A2−A3=c11+c22c11c22−c12c21+c13c11c31+c21c32+c23c22c32+c12c31.A reminder that according to the Routh-Hurwitz criterion, the characteristic equation (35) has negative real part eigenvalues, and then, the CP becomes L.A.S if and only if A1>0, A3>0, and ∆>0. Accordingly, the following theorem can be proved simply.Theorem 2.
The CP is L.A.S if the below sufficient conditions hold.(38a)2α3x∗2y∗Γ1∗2+2α4x∗2z∗Γ2∗2+2α5α8y∗2z∗Γ2∗2<x∗Γ1∗2Γ2∗2,(38b)2α5y∗z∗Γ2∗Γ3∗2<α1Γ1∗2Γ2∗+Γ1∗2Γ3∗2,(38c)α3x∗2<1,(38d)c22c21<α13α12<maxc12c11,−c11c21,−c12c22,(38e)α1Γ1∗2Γ2∗+Γ1∗2Γ3∗2−2α5y∗z∗Γ2∗Γ3∗2α61−α3x∗2Γ1∗2+2α8α4x∗z∗Γ2∗2Γ4∗Γ2∗Γ3∗2Γ4∗Γ1∗2Γ2∗2−2α3x∗y∗Γ1∗2−2α4x∗z∗Γ2∗2α5α8y∗z∗<2.


## 4. Persistence

In this part, the persistence of system (4) is investigated; it is fully recognized that the system is said to persist if and only if none of their species extinct; this means that system (4) persists if the trajectory of the system that initiates at a positive initial point does not have omega limit set on the boundary planes of its domain.System (4) has two subsystems lying in xy-plane and xz-plane, respectively, which can be represented as follows:
(39)dxdt=x11+α1y−x−y1+α3x2=h1x,y,dydt=yα6x1+α3x2−α9=h2x,y,and
(40)dxdt=x11+α2z−x−z1+α4x2=G1x,z,dzdt=Z2α10−α10α11+α12x=G2x,z.It can be simply confirmed that the above subsystems (39) and (40) have positive equilibrium points that coincide with those of system (4) in the interior of boundary planes xy-plane and xz-plane, respectively. Now, to discover the possibility of the existence of periodic dynamics in the boundary planes, the Dulac function approach is used.Consider the following functionB1x,y=1/xy, clearly this function B1x,y>0 and C1 function in the interior of ℝ+2 of xy-plane. Moreover, it is obvious that
(41)∆x,y=∂∂xB1∙h1+∂∂yB1∙h2=−1y+2α3x1+α3x22.Therefore,∆x,y does not alter sign and not vanish under the following condition:
(42)2α3x1+α3x22<1y.Hence, if condition (42) holds, there are no periodic dynamics in the interior of a positive quadrant of xy-plane for subsystem (39).Similarly, with the help of Dulac functionB2x,z=1/xz2, which is B2x,z>0, and C1 function in the interior of ℝ+2 of xz-plane, it is observed that there are no periodic dynamics in the interior of a positive quadrant of xz-plane for subsystem (40) provided that the following condition holds:
(43)2α4x1+α4x22<1z.Theorem 3.
Suppose that the boundary planes have no periodic dynamics; then, system (4) is uniformly persistent as long as that the next conditions hold:
(44a)α9<α61+α3,(44b)1α11+α12<1,(44c)z^1+α2z^<1,(44d)1α11+α12x¯+α13y¯<1,(44e)α8x̿1+α4x̿2+α9<α6x̿1+α3x̿211+α7z̿.Proof.
Define a functionϑx,y,z=xp1yp2zp3, where p1,p2,andp3 are positive constants, and ϑx,y,z>0 for all x,y,z∈Intℝ+3 with ϑx,y,z=0 if any one of x,y, or z approaches zero. Therefore, direct computation gives
(45)Ωx,y,z=ϑ′x,y,zϑx,y,z=p1f1+p2f2+p3f3,where the functions fi;i=1,2,3 are given in system (4). Now, according to the average Lyapunov method, the proof is satisfied provided that Ωx,y,z>0 for all boundary equilibrium points. Therefore,
(46)Ωx,y,z=p11Γ3−x−yΓ2−zΓ1+p2α6xΓ2Γ4−α8zΓ1−α9+p3α10−α10Γ5,where Γi;i=1,2,3,4,5 are given in equation (20). Clearly, we have that
(47)ΩE0=p11+p2−α9.
Obviously, by choosing the arbitrary positive value ofp1 sufficiently larger than that of p2, it is obtained that ΩE0>0.
(48)ΩE1x=p2α61+α3−α9+p3α10−α10α11+α12.
Note that conditions (44a) and (44b) guarantee that ΩE1x>0.
(49)ΩE1z=p111+α2z^−z^+p2−α8z^+α9.
Again, by choosing the arbitrary positive value ofp1 sufficiently larger than that of p2, it is obtained that ΩE1z>0 provided that condition (44c) holds.
(50)ΩE2=p3α10−α10α11+α12x¯+α13y¯.
Clearly, condition (44d) guarantees that ΩE2>0.
(51)ΩE3=p2α6x̿1+α3x̿211+α7z̿−α8x̿1+α4x̿2−α9.
Finally,ΩE3>0 under condition (44e). Hence, system (4) is uniformly persistent, and the proof is done.☐


## 5. Region of Attraction

In this section, the region of attraction for each equilibrium point of system (4) is determined using the direct method of Lyapunov. These regions of attraction for each of their equilibrium point are known as the domain or basin of attraction.Theorem 4.
The FAP is asymptotically stable (A.S) in the subregion ofℝ+3 that satisfies the following sufficient conditions:
(52a)1<α9,(52b)α1<1−α61+α3,(52c)α2<11+α4+α5α6/4α9+α62,(52d)14α10Gα9−1<y,where the symbol G is given in the proof.Proof.
DefineL1x,y,z=∫1xu−1/udu+y+z. Clearly, the function L1x,y,z is a positive definite function that is L11,0,0=0, while L1x,y,z>0, for all values in the region x,y,z∈ℝ+3:x>0,y≥0,z≥0;x,y,z≠1,0,0. Then, using some algebraic manipulation gives that
(53)dL1dt≤−x−121+α1y+α2z+xα1y+α2z1+α1y+α2z−xy1+α3x2+y1+α3x2−xz1+α4x2+α5y2+z1+α4x2+α5y2+α6xy1+α3x211+α7z−α9y+α10z2−Z2α11+α12x+α13y.
Consequently, by using additional computation, the following is obtained:(54)dL1dt≤−x−121+α1y+α2z+α1xy+α2xz−α9−1y−1−α6xy1+α3x2−xz1+α4x2+α5y2+z−α10z21α11+α12x+α13y−1.
Now, according to the upper bound forx and y, so as t⟶∞, we have x≤1 and y≤α6/4α9+α6. Therefore, it is obtained that
(55)dL1dt≤−x−121+α1y+α2z−α9−1y−1−α61+α3−α1xy−11+α4+α5α6/4α9+α62−α2xz+z1−α10z1α11+α12+α13α6/4α9+α6−1.
Consequently, using conditions (52b) and (52c) with the bound of logistic term yields
(56)dL1dt≤−x−121+α1y+α2z−α9−1y−1−α61+α3−α1xy−11+α4+α5α6/4α9+α62−α2xz+14α10G,where G=1/α11+α12+α13α6/4α9+α6−1.
Accordingly, condition (52d) guarantees that dL1/dt<0, which means it is a negative definite. Then, L1 is a Lyapunov function, and hence, the FAP is A.S for any trajectory starting from a point that belongs to the region that satisfies the above conditions.☐Note that the subregion described in Theorem4 represents the basin of attraction of the FAP.Theorem 5.
The SAP is A.S in the subregion ofℝ+3 that satisfies the following sufficient conditions:
(57a)z^<minα11Γ5α10α12Γ1,α8α11Γ5α10α13Γ1,(57b)x+α10Nz−z∧2<α9y,where the symbol N>0 is given in the proof.Proof.
DefineL2x,y,z=x+y+∫z^zm−z^/zdm. Clearly, the function L2x,y,z is a positive definite function that is L20,0,z^=0, while L2x,y,z>0, for all values in the region x,y,z∈ℝ+3:x≥0,y≥0,z>0;x,y,z≠0,0,z^. Then, using some algebraic manipulation gives that
(58)dL2dt≤x−xz1Γ1−α10α12z^α11Γ5−yzα8Γ1−α10α13z^α11Γ5−α9y+α10Nz−z∧2,where N=1−1/α11+α12x+α13y; clearly, N is positive at the SAP.
Consequently, by using conditions (57a) and (57b), it is obtained that dL2/dt<0. Hence, the derivative of the function L2 is a negative definite, and then, L2 is a Lyapunov function. Hence, the SAP becomes A.S point for any point that belongs to the subregion that satisfies the given conditions and the proof is done.☐Theorem 6.
The GPFP is A.S in the subregion ofℝ+3 that satisfies the following sufficient conditions:
(59a)α3y¯x+x¯Γ2Γ¯2<1,(59b)Γ5<1,(59c)ρ122<ρ11,(59d)ρ132<ρ11ρ33,(59e)ρ232<ρ33,(59f)y−y¯2<12y−y¯+ρ33z2,where all the symbols are given in the proof.Proof.
DefineL3x,y,z=∫x¯xm−x/mdm+∫y¯yw−y/wdw+z. Clearly, the function L3x,y,z is a positive definite function that is L3x¯,y¯,0=0, while L3x,y,z>0, for all values in the region x,y,z∈ℝ+3:x>0,y>0,z≥0;x,y,z≠x¯,y¯,0. Then, using some algebraic manipulation gives that
(60)dL3dt=x−x¯−α1y−y¯Γ3Γ¯3−α2zΓ3Γ¯3−x−x¯−y−y¯Γ2+α3y¯Γ2Γ¯2x+x¯x−x¯−zΓ1+y−y¯α6Γ2Γ¯2Γ41−α3x¯xx−x¯−α6α7x¯zΓ¯2Γ4−α8zΓ1+α10z21−1Γ5.
Furthermore, computation gives(61)dL3dt=−1−α3y¯Γ2Γ¯2x+x¯x−x¯2−α1Γ3Γ¯3+1Γ2−α61−α3x¯xΓ2Γ¯2Γ4x−x¯y−y¯−α2Γ3Γ¯3+1Γ1x−x¯z−α6α7x¯Γ¯2Γ4+α8Γ1y−y¯z−α10z21Γ5−1.
Therefore, using conditions (59a), (59b), (59c), and (59d) gives that
(62)dL3dt≤−12ρ11x−x¯+y−y¯2−12ρ11x−x¯+ρ33z2−12y−y¯+ρ33z2+y−y¯2,where ρ11=1−α3y¯/Γ2Γ¯2x+x¯, ρ12=α1/Γ3Γ¯3+1/Γ2−α61−α3x¯x/Γ2Γ¯2Γ4, ρ13=α2/Γ3Γ¯3+1/Γ1, ρ23=α6α7x¯/Γ¯2Γ4+α8/Γ1, ρ33=α101/Γ5−1.
Accordingly, using condition (59f), it is observed that dL3/dt<0, which means that L3 is a Lyapunov function and hence the GPFP is an A.S for any trajectory starting from a point that belongs to the region that satisfies the above conditions.☐Theorem 7.
Assume that the SPFP is L.A.S; then, it has a basin of attraction that belongs to theℝ+3 that satisfies the following conditions:
(63a)α4z̿x+x¯Γ1Γ1̿<1,(63b)F5<1,(63c)q132<4q11q33,(63d)α5z̿Γ1Γ1̿y<α1Γ3Γ3̿,(63e)1Γ2+α1Γ3Γ3̿x̿<α9,(63f)α10α13z̿Γ5Γ5̿<α8Γ1,where the symbols q11,q33, and q13 are given in the proof.Proof.
DefineL4x,y,z=∫x̿xm−x̿/mdm+y+∫z̿zw−z̿/wdw. Clearly, the function L4x,y,z is a positive definite function that is L4x̿,0,z̿=0, while L4x,y,z>0, for all values in the region x,y,z∈ℝ+3:x>0,y≥0,z>0;x,y,z≠x̿,0,z̿. Then, using some algebraic manipulation gives that
(64)dL4dt≤−1−α4z̿x+x̿Γ1Γ1̿x−x̿2−α101Γ5−1z−z̿2−α2Γ3Γ̿3+1Γ1−α10α12z̿Γ5Γ̿5x−x̿z−z̿−α9−1Γ2+α1Γ3Γ3̿x̿y−1−α6Γ4xyΓ2−α8Γ1−α10α13z̿Γ5Γ5̿yz−α1Γ3Γ3̿−α5z̿Γ1Γ1̿yxy.
Consequently, using conditions (63a), (63b), (63d), (63e), and (63f) gives
(65)dL4dt≤−q11x−x̿2−q33z−z̿2−q13x−x̿z−z̿−α9−1Γ2+α1Γ3Γ3̿x̿y.
Moreover, using condition (63c) leads to
(66)dL4dt≤−q11x−x̿+q33z−z̿2−α9−1Γ2+α1Γ3Γ3̿x̿y.
Accordingly,dL4/dt<0, which means that L4 is a Lyapunov function and hence the SPFP is A.S for any trajectory starting from a point that belongs to the region that satisfies the above conditions.☐Theorem 8.
Assume that the CP is L.A.S; then, it has a basin of attraction that belongs to theℝ+3 that satisfies the following conditions:
(67a)d122<d11d22,(67b)d132<d11d33,(67c)d232<d22d33,(67d)α3y∗x+x∗Γ2Γ2∗+α4z∗x+x∗Γ1Γ1∗<1,(67e)α3xx∗<1,(67f)α5y∗2<α4x∗2,where the symbols dij,i,j=1,2,3 are given in the proof.Proof.
Define(68)L5x,y,z=∫x∗xm−x∗mdm+∫y∗yw−y∗wdw+∫z∗zn−z∗ndn.
Clearly, the functionL5x,y,z is a positive definite function that is L5x∗,y∗,z∗=0, while L5x,y,z>0, for all values in the region x,y,z∈ℝ+3:x>0,y>0,x>0;x,y,z≠x∗,y∗,z∗. Thus, after some algebraic manipulation, it is obtained that
(69)dL5dt=−d112x−x∗2−d12x−x∗y−y∗−d222y−y∗2−d332z−z∗2−d13x−x∗z−z∗−d23y−y∗z−z∗,where
(70)d11=1−α3y∗x+x∗Γ2Γ2∗+α4z∗x+x∗Γ1Γ1∗,d12=1−α1Γ3Γ3∗+1+α3x∗Γ2Γ2∗+α3+α7zα6xx∗−α61+α7z∗Γ2Γ2∗Γ4Γ4∗+α5z∗y+y∗−α4α8x+x∗Γ2Γ2∗,d22=α5α8z∗y+y∗Γ1Γ1∗,d33=α10+α10α11+α13x∗Γ5Γ5∗+α10α13y∗Γ2Γ2∗Γ4Γ4∗,d13=α2Γ3Γ3∗+1Γ1−α10α12z∗Γ5Γ5∗,d23=α6α7x∗1−α3x∗Γ2Γ2∗Γ4Γ4∗+α8Γ1Γ1∗+α4α8x∗2Γ1Γ1∗+α8α5y∗2Γ1Γ1∗−α10α13z∗Γ5Γ5∗.
Consequently, using the above conditions gives(71)dL5dt≤−12d11x−x∗+d22y−y∗2−12d22y−y∗+d33z−z∗2−12d11x−x∗+d33z−z∗2.
Obviously,dL5/dt<0, which gives that L5 is a Lyapunov function and hence the CP is A.S for any trajectory starting from a point that belongs to the region that satisfies the above conditions.☐


## 6. Local Bifurcation

An investigation of the impact of varying the parameter values on the dynamics of system (4) around each equilibrium point is done using Sotomayor’s theorem for local bifurcation. Because the existence of a nonhyperbolic equilibrium point represents a necessary but not sufficient condition for bifurcation to occur, then the candidate bifurcation parameter is choosing so that the equilibrium point becomes a nonhyperbolic point. Rewrite system (4) in the following general form:
(72)dXdt=FX,X=x,y,zT,F=xf1,yf2,zf3T.Then, the second directional derivative of system (4) can become as
(73)D2FXΦ,Φ=φi13×1,where
(74)φ11=−2+2α3xy3−α3x2Γ23+2α4xz3−α4x2+3α5y2Γ13v12+−2α1Γ32−21−α3x2Γ22+α5yz4+4α5y2−10α4x2Γ13v1v2−2α2Γ32+21−α4x2+α5y2Γ1v1v3+2α12xΓ33+2α5xz1+α4x2−2α5y2Γ13v22+4xα1α2Γ33+α5yΓ12v2v3+2α22xΓ33v32,φ21=2α3α6xy−3+α3x2Γ23Γ4+2α4α8yz1+α5y2−2α4x2Γ13v12+α4α8xz4+4α4x2−10α5y2Γ13+2α61−α3x2Γ22Γ4v1v2+−2α6α7y1−α3x2Γ22Γ42+4α4α8xyΓ12v1v3+2α5α8yz3−α5y2+3α4x2Γ13v22+α6α72xyΓ2Γ43v32−2α6α7xΓ2Γ42+2α81+α4x2−α5yzΓ12v2v3,φ31=−α10α122z2Γ53v12−2α10α12α13z2Γ53v1v2+4zα10α12Γ52v1v3−α10α132z2Γ53v22+4zα10α13Γ52v2v3+2α10−α10Γ5v32,with Φ=v1,v2,v3T be any nonzero vector, and the symbols Γi;i=1,2,3,4,5 are given in equation (20). Accordingly, the local bifurcation near the equilibrium points is studied in the next theorems.Theorem 9.
Suppose that condition (32a) is satisfied; then, system (4) at the SPFP undergoes a transcritical bifurcation when the parameter α9 passes through the value α9∗=α6x̿/Γ̿2Γ̿4−α8z̿/Γ̿1, provided that the following condition holds:
(75)4α4α8x̿z̿1+α4x̿2Γ̿13+2α61−α3x̿2Γ̿22Γ̿4ρ1−2α6α7x̿Γ̿2Γ̿42+2α81+α4x̿2Γ̿12ρ2≠0,where the symbols ρ1 and ρ2 are given in the proof.Proof.
Direct computation shows that the JM of system (4) at E3,α9∗ can be represented as
(76)J1=x¯¯−1+2α4x¯¯z̿Γ̿12x¯¯−α1Γ̿32−1Γ̿2x¯¯−α2Γ̿32−1Γ̿1000α10α12z̿2Γ̿52α10α13z̿2Γ̿520=bij3×3,where Γ̿i;i=1,2,3,4,5 are given in equation (28). Obviously, the matrix J1 has two eigenvalues having negative real parts due to condition (32a) and the third one is zero, say λ32∗=0. Hence, E3 is a nonhyperbolic point.
LetΘ1=θ11,θ12,θ13T be the eigenvector corresponding to the eigenvalue λ32∗=0. Thus, J1Θ1=0 gives that Θ1=−b32/b31θ12,θ12,b11b32−b31b12/b31b13θ12T=ρ1θ12,θ12,ρ2θ12T, where θ12≠0 be any real number.
Now, letY1=ϵ11,ϵ12,ϵ13T represent the eigenvector corresponding to the eigenvalue λ32∗=0 of the matrix J1T. Thus, J1TY1=0 gives that Y1=0,ϵ12,0T, where ϵ12≠0 be any real number. Now, according to Sotomayor’s theorem, it is obtained that
(77)∂F∂α9=Fα9X,α9=0,−y,0T⇛∂F∂α9=Fα9E3,α9∗=0,0,0T.
Therefore,Y1TFα9E3,α9∗=0; hence, system (4) has no saddle-node bifurcation. Moreover, since
(78)DFα9X,α9=0000−10000⇛DFα9E3,α9∗Θ1=0,−θ12,0T.
Then,Y1TDFα9E3,α9∗Θ1=−θ12ϵ12≠0.
Also, by using equation (73), it is obtained that
(79)D2FE3,α9∗Θ1,Θ1=φ̿i13×1,where
(80)φ̿11=−2+2α4x̿z̿3−α4x̿2Γ̿13ρ1θ122−2α1Γ̿32+21−α3x̿2Γ̿22ρ1θ122+4x̿α1α2Γ̿33ρ2θ122−2α2Γ̿32+21−α4x̿2Γ̿1ρ1ρ2θ122+2α12x̿Γ̿33+2α5x̿z̿1+α4x̿2Γ̿13θ122+2α22x̿Γ̿33ρ2θ122,φ̿21=4α4α8x̿z̿1+α4x̿2Γ̿13+2α61−α3x̿2Γ̿22Γ̿4ρ1θ122−2α6α7x̿Γ̿2Γ̿42+2α81+α4x̿2Γ̿12ρ2θ122,φ̿31=−α10α122z̿2Γ̿53ρ1θ122−2α10α12α13z̿2Γ̿53ρ1θ122+4z̿α10α12Γ̿52ρ1ρ2θ122−α10α132z̿2Γ̿53θ122+4z̿α10α13Γ̿52ρ2θ122+2α10−α10Γ̿5ρ2θ122.
Accordingly, it is easy to verify that(81)Y1TD2FE3,α9∗Θ1,Θ1=4α4α8x̿z̿1+α4x̿2Γ̿13+2α61−α3x̿2Γ̿22Γ̿4ρ1−2α6α7x̿Γ̿2Γ̿42+2α81+α4x̿2Γ̿12ρ2θ122ϵ12.
Clearly,Y1TD2FE3,α9∗Θ1,Θ1≠0, under condition (75), and then, a transcritical bifurcation take place in the sense of Sotomayor.☐Theorem 10.
Suppose that conditions (38a)–(38c) and (38e) are satisfied; then, system (4) at CP undergoes a saddle-node bifurcation when the parameter α12 passes through the value α12∗=α13c11c23−c13c21/c12c23−c13c22, provided that the following condition holds:
(82)ρ5φ11∗+ρ6φ21∗+φ31∗≠0,where φi1∗;i=1,2,3 are given in the proof.Proof.
Direct computation shows that the JM of system (4) at E4,α12∗ can be represented as
(83)J2=c11c12c13c21c22c23c31α12∗c32α12∗0,where all the coefficients cij;i,j=1,2,3 are given in equation (33). Clearly, the determinant of J2 that is given by A3 in equation (35) is zero at α12∗. Therefore, the characteristic equation of J2 has zero root, say λ∗=0, and then, the coexistence equilibrium point becomes a nonhyperbolic point.
LetΘ2=θ21,θ22,θ23T be the eigenvector corresponding to the eigenvalue λ∗=0. Thus, J2Θ2=0 gives that Θ2=ρ4θ23,ρ3θ23,θ23T, where θ23≠0 be any real number, and ρ3=c13c21−c11c23/c11c22−c12c21, while ρ4=c23c12−c13c22/c11c22−c12c21. Clearly, ρ3>0, and ρ4<0, due to conditions (38a)–(38c) with (38e).
Now, letY2=ϵ21,ϵ22,ϵ23T represent the eigenvector corresponding to the eigenvalue λ∗=0 of the matrix J2T. Thus, J2TY2=0 gives that Y2=ρ5ϵ23,ρ6ϵ23,ϵ23T, where ϵ23≠0 be any real number, and ρ5=−c23c31/c11c23−c13c21, while ρ6=c13c31/c11c23−c13c21. Clearly, ρ5>0, and ρ6<0, due to conditions (38a)–(38c) with (38e).
Now, according to Sotomayor’s theorem, it is obtained that(84)∂F∂α12=Fα12X,α12=0,0,α10xz2Γ52T,where Γ5 is given in equation (20). Therefore, it is obtained that
(85)Y2∗Fα12E4,α12∗=α10x∗z∗2Γ5∗2ϵ23≠0,where Γ5∗ is given in equation (33). Now, by using equation (73), it is obtained that
(86)D2FE4,α12∗Θ2,Θ2=φi1∗3×1,where
(87)φ11∗=−2−2α3x∗y∗3−α3x∗2Γ2∗3−2α4x∗z∗3−α4x∗2+3α5y∗2Γ1∗3ρ4θ232−2α1Γ32+21−α3x∗2Γ2∗2−α5y∗z∗4+4α5y∗2−10α4x∗2Γ1∗3ρ3ρ4θ232−2α2Γ3∗2+21−α4x∗2+α5y∗2Γ1∗ρ4θ232+4x∗α1α2Γ3∗3+α5y∗Γ1∗2ρ3θ232+2α12x∗Γ33+2α5x∗z∗1+α4x∗2−2α5y∗2Γ1∗3ρ3θ232+2α22x∗Γ3∗3θ232,φ21∗=2α3α6x∗y∗−3+α3x∗2Γ2∗3Γ4∗+2α4α8y∗z∗1+α5y∗2−2α4x∗2Γ1∗3ρ4θ232+α4α8x∗z∗4+4α4x∗2−10α5y∗2Γ1∗3+2α61−α3x∗2Γ2∗2Γ4∗ρ3ρ4θ232−2α6α7y∗1−α3x∗2Γ2∗2Γ4∗2−4α4α8x∗y∗Γ1∗2ρ4θ232+α6α72x∗y∗Γ2∗Γ4∗3θ232+2α5α8y∗z∗3−α5y∗2+3α4x∗2Γ1∗3ρ3θ232−2α6α7x∗Γ2∗Γ4∗2+2α81+α4x∗2−α5y∗z∗Γ1∗2ρ3θ232,φ31∗=−α10α12∗2z∗2Γ5∗3ρ4θ232−2α10α12∗α13z∗2Γ5∗3ρ3ρ4θ232+4z∗α10α12∗Γ5∗2ρ4θ232−α10α132z∗2Γ5∗3ρ3θ232+4z∗α10α13Γ5∗2ρ4θ232+2α10−α10Γ5∗θ232.
Accordingly, it is simple to verify that(88)Y2TD2FE4,α12∗Θ2,Θ2=ϵ23ρ5φ11∗+ρ6φ21∗+φ31∗.
Clearly,Y2TD2FE4,α12∗Θ2,Θ2≠0, under condition (82), and then, a saddle-node bifurcation takes place in the sense of Sotomayor.☐


## 7. Numerical Simulation

In this part, the dynamics of system (4) are studied numerically. The objective is to complete the vision of the global dynamics of the system and detect the effect of varying the parameter values on the system’s dynamical behavior. Two sets of hypothetical parameter values are used to satisfy our objective. For the following set of parameter values, other sets can be used to study the system:
(89)α1=0.2,α2=0.2,α3=0.1,α4=0.1,α5=0.1,α6=0.25,α7=0.2,α8=0.75,α9=0.1,α10=0.75,α11=0.1,α12=0.9,α13=1.5.It is watched that system (4) has a globally asymptotically stable CP as shown in Figure 1.Figure 1
(a) The trajectories of system (4) approach asymptotically to CP, starting from different initial points with the data given by (89). (b) The trajectory of x versus t. (c) The trajectory of y versus t. (d) The trajectory of z versus t.
(a)(b)(c)(d)According to Figure1, system (4) is globally asymptotically stable given by E4=0.68,0.18,0.08 which means system (4) is persistent at the CP. However, for the data (89) with increasing α1, system (4) loses its persistence and approaches asymptotically to GPFP as shown in Figure 2 for the typical value of α1.Figure 2
Forα1=2, with the data given by (89). (a) The trajectory of x versus t. (b) The trajectory of y versus t. (c) The trajectory of z versus t. (d) The trajectory of system (4) approaches asymptotically to GPFP.
(a)(b)(c)(d)On the other hand, varying the other fear’s parameters leads to a quantitative change in the position of the CP and system (4) still persistent. Now, the effects of varying the parameters α3 and α5 are presented in Figures 3 and 4, respectively.Figure 3
(a) The trajectories of system (4) for the parameters given by (89) with different values of α3=0.5,1,and5 represented by the colors blue, green, and red, respectively, approach asymptotically to CP for the first two values, while the trajectory approaches to SPFP for the third value. (b) The trajectories of x versus t. (c) The trajectories of y versus t. (d) The trajectories of z versus t.
(a)(b)(c)(d)Figure 4
(a) The trajectory of system (4) approaches asymptotically to periodic dynamics in the positive octant for α5=20, with the other parameters as given by (89). (b) Time series for the trajectory given in (a).
(a)(b)According to Figure3, increasing the value of α3 causes decreasing in the predator’s populations and increasing in prey population gradually; then, the solution of system (4) settled at the SPFP first and then at the FAP. However, increasing α5 destabilizes system (4) and the solution persists at periodic dynamics. Moreover, it is observed that the effect of the parameter α4 is similar to that of α3.Now, the effects of varying the parametersα6 and α9 on system (4) dynamical behavior are shown in Figures 5 and 6, respectively.Figure 5
The trajectories of system (4) using the data given by (89) with different values of α6=0.1,0.5,and0.75, those represented by dotted line, dashed line, and solid line, respectively.Figure 6
The trajectories of system (4) using the data given by (89) with different values of α9=0.01,0.15,and0.25, those represented by dotted line, dashed line, and solid line, respectively.Obviously, Figures5 and 6 show that the parameters α6 and α9 have opposite effects on the dynamic behavior of system (4) so that as increases α6 (decreasing α6), the population of the prey decreases while the populations of the predators increase and hence system (4) loses their persistence and approaches to FAP as decreasing α6 (increasing α9). Furthermore, it is observed that varying the parameter α10 has a quantitative effect on the dynamic behavior of system (4) and system (4) persists at the CP.The effects of varying the parametersα11 and α12 on system (4) dynamical behavior are explaining in Figures 7 and 8, respectively.Figure 7
(a) The trajectories of system (4) using the data given by (89) with α11=0.24 approaches asymptotically to SPFP. (b) Time series for the trajectory given in (a). (c) The trajectories of system (4) using the data given by (89) with α11=1 approaches asymptotically to SAP. (d) Time series for the trajectory given in (c).
(a)(b)(c)(d)Figure 8
(a) The trajectories of system (4) using the data given by (89) with α12=0.1 approaches asymptotically to GPFP. (b) Time series for the trajectory given in (a). (c) The trajectories of system (4) using the data given by (89) with α12=1.25 approaches asymptotically to SPFP. (d) Time series for the trajectory given in (c).
(a)(b)(c)(d)According to Figures7 and 8, it is watched that increasing the parameter α11 leads to losing the persistence of the system by approaching the SPFP first and then to SAP. However, decreasing (or increasing) the parameter α12 leads to losing the persistence of the system and the solution approaches to GPFP (SPFP). It is observed that the parameter α13 has similar effects on the dynamic behavior of system (4) as that for α12.Keeping the above in mind, the second set of parameter values, which leads system (4) to chaotic dynamics, is given in the set of data (90):
(90)α1=0.8,α2=0.8,α3=4,α4=8,α5=4,α6=1,α7=0.8,α8=0.25,α9=0.05,α10=0.5,α11=0.1,α12=4,α13=4.The dynamic behavior of system (4) is shown in Figures 9 and 10.Figure 9
(a) The trajectory of system (4) approaches a chaotic attractor using the dataset given by (90), while their projection on xy,xz, and yz planes are shown in (b), (c), and (d), respectively.
(a)(b)(c)(d)Figure 10
Time series of the chaotic attractors given in Figure9.According to Figures9 and 10, system (4) is rich in its dynamics including point attractors, periodic attractors, and chaos. It has many bifurcation points when varying its parameter values. It is well known that the most property that characterizes the chaotic attractor is the sensitivity to the varying initial points. So, Figure 11 shows clearly the sensitivity of the strange attractor given in Figure 9 to a small change in their initial point.Figure 11
Sensitivity of the trajectory of the chaotic attractor that is given in Figure9(a) to the initial values. Solid line starts at (0.8, 0.7, 0.6) while the dotted line starts at (0.81, 0.71, 0.61). (a) Trajectories of x versus t. (b) Trajectories of y versus t. (c) Trajectories of z versus t.
(a)(b)(c)Now, the impacts of varying some parameter values on the chaotic regions and the dynamic behavior of system (4) are studied using the bifurcation diagrams. It is known that the bifurcation diagram shows the birth and death of the attractors as a function of varying the parameter values. The bifurcation diagrams (BD) as a function of the parameters α1,α7,α9,α12, and α13 are shown in Figures 12–16.Figure 12
The BD of system (4) as a function of α1 in the range 0.7<α1<1.1 using the parameter values in (90).Figure 13
The BD of system (4) as a function of α7 in the range 0.1<α7<2 using the parameter values in (90).Figure 14
The BD of system (4) as a function of α9 in the range 0.01<α9<0.1 using the parameter values in (90).Figure 15
The BD of system (4) as a function of α12 in the range 0.01<α12<6 using the parameter values in (90).Figure 16
The BD of system (4) as a function of α13 in the range 1<α13<8 using the parameter values in (90).According to Figures12 and 13, there are ranges around α1=0.8 and α7=0.8 in which the system is chaotic; however, it faces extinction in the specialist predator as these parameters increase and the trajectories approach periodic dynamics in xz-plane. Similar effects on the dynamics behavior of system (4) are observed for the parameters α2 as those of α1 and α7. Finally, the other BD are clearly shown the existence of chaotic regions; however, it is watched that increasing these parameters leads to extinction in one predator species and the solution approaches to periodic dynamics in the plane.


## 8. Discussion and Conclusion

In this paper, a food-web model having specialist and generalist predators is proposed and studied. It is assumed that the prey has antipredator properties against predation and hence Sokol-Howell functional responses are used to satisfy this assumption. The generalist predator has an alternative food source, and hence, the Leslie-Gower type generalist predator is also included in the system. The impact of fear on the growth rate of the prey due to their fear of predation is also considered. All the properties of the solution of the model are discussed. The equilibria and their stability analysis are studied. The persistence conditions are established. The regions of asymptotic stability of these equilibria are determined with the help of the Lyapunov method. Local bifurcation is investigated whenever it exists. Finally, extensive numerical simulations are performed to understand the effects of varying the parameter values on the dynamics of the system using two sets of parameters. It is watched that system (4) is wealthy in its dynamical behavior including all possible attractors’ types such as point attractor, periodic attractor, and chaos. The fear coefficients have a clear control on the dynamics of the system so that as the fear increases, the chaotic regions decrease and the solution settled to CP or periodic dynamics in the boundary planes. Increasing the preference rates, which are equivalent to increasing the prey’s group defense, leads to loss of the system persistence and the solution approaches asymptotically to FAP. Increasing the death rate of the specialist predator causes extinction in both the predators and the system approaches to FAP. Finally, system (4) is very sensitive to changes in the abundance of alternative food sources.



---

*Source: 2112814-2021-09-01.xml*