# Dynamical Behavior of a Delayed Holling Type-II Predator-Prey Model with Predator Cannibalism

Authors: R. Lavanya; S. Vinoth; K. Sathiyanathan; Zeric Njitacke Tabekoueng; P. Hammachukiattikul; R. Vadivel
Journal: Journal of Mathematics
 (2022)

---

## Abstract

In this article, we propose the Holling type-II predator-prey model involving cannibalism and gestation delay in predators. We study the existence of all possible equilibrium points of the proposed model. We give the condition for the local stability and Hopf bifurcation analysis for the nondelayed model. Next, we also establish the local stability and Hopf bifurcation analysis for the corresponding delayed model. Finally, we discuss how cannibalism and delay play an important role in stabilizing and destabilizing the proposed system both theoretically and numerically.


---

## Body

## 1. Introduction

Studies of predator-prey interactions have been considered among the most challenging problems in population ecology. It is a study project that is relevant not just to ecologists and biologists but also to mathematicians because it considers the interactions of different species. Its importance can be seen through the large number of proposed models, which describe the interaction between the prey and predator in different cases. Increasing interest has been noticed in modeling those types of interactions. The basic model for the interaction of two species is modeled by Lotka in 1925 and Volterra in 1926. The predator-prey model is influenced by many biological factors such as prey refuge [1], fear [2], and the Allee effect [3, 4]. In such models, the predator consumes its prey is given in the form of functional response, which is the function of its lone prey or both predator and its prey. Also, the predator dies out exponentially in the absence of the predator. The other important biological phenomenon is called cannibalism, which refers to the act of killing and at least partially consuming its species. The energy from its species helps to grow its population size. Recently, there has been much attention on the study of the behavior of the prey-predator model involving cannibalism. In Reference [5], the authors extended the classical Lotka–Volterra model with cannibalism in predators, which is of the form,(1)dX1dT=X1r1−X1K−A1X2,dX2dT=X2−D+M1+A2X1−GX2H+X2,where X1 and X2 are the sizes of prey and predator population at the time T; r and K denote the growth rate and carrying capacity of prey; A1 denotes the reduction rate of prey; D represents the death rate of predator; A2 denotes the conversion efficiency of the predator; GγX2/H+X2 is the cannibalism term, γ denotes the rate of cannibalism; M1<G, M1 denotes the birth rate from the predator cannibalism; and H is the positive constant.Lin et al. [6] studied the predator-prey model with cannibalism in prey population, which is of the form,(2)dX1dT=Xr+M1−rX1K−A1X1−GX1H+X1,dX2dT=X2−d+A2X1,where X1 and X2 are the size of prey and predator at time T. r, k, A1, G, H, D, M1, and A2 are all positive constants.One of the key terms for studying the predator-prey model is the interaction term. The most common type of interaction between prey and predator is given by the Holling types [7–9]. In Reference [10], the authors studied the predator-prey model with Holling type-II functional response, which is of the form(3)dX1dT=X1r1−X1K−A1X2X+B1,dX2dT=X2−d+A2X1X1+B2,where X1 and X2 are the sizes of prey and predator at time T. B1 and B2 are the half-saturation constants for the prey and predator. In this study, we assumed B1=B2 and considered as B. Model (3) has been assumed that the prey interacts with a predator in the form of Holling type-II function, and in the absence of prey species, the predator population dies exponentially. We, therefore, formulate the model by considering cannibalism, which is the assumption that the predator population gets energy from its species; that is, the predator can survive without its lone prey. Recent empirical and theoretical work has shown that cannibalistic interactions can have important consequences for species interactions and the structure of whole communities and can result in dynamics that are not predictable from unstructured systems [11]. Mathematical biologists have spent a lot of time investigating cannibalism because of the strange consequences it can have on population dynamics. Cannibalism can be used as a survival mechanism in a population on the verge of extinction through the employment of the lifeboat mechanism [12]. Cannibalism in both prey and predator population has been considered in Reference [13].On the other hand, some researchers have studied the presence of time delay in dynamic systems. The time delay in the biological system is unavoidable. In ecology, there are many reasons for incorporating delays, such as delay in maturation and gestation time. The energy attained by consuming food will not take place immediately, and there is a time delay during the gestation process. The study of the dynamical behavior of the predator-prey model with delay is very significant. The presence of delay in the dynamical system can affect the stability of the system, that is, stabilizing and destabilizing effects. For instance, the persistence and extinction for delayed stochastic prey-predator system with hunting cooperation in predators have been studied in Reference [14]. Also, the authors derived the Lyapunov functional to provide adequate conditions for persistence and extinction. In order to perform the numerical simulations, they used Milstein’s method. They found that the smaller white noise can assist the survival of both species, but higher the strength of noise can lead to the extinction of the predator. The authors in Reference [15] studied a two-prey, one-predator food chain model with the Allee effects in each species and two unique delays. Sufficient conditions for the local stability of coexisting equilibrium and occurrence of Hopf bifurcations in terms of both delays are established. The influence of the Allee effect and time delays in the model raises the complexity of the model and enriches the system dynamics. The presence of time delay can cause complex behavior in the predator-prey model, and other dynamical systems have been found in the literature [16–21].From the existing literature, the study of the Holling type-II predator-prey model with cannibalism and gestation delay in predators has not been considered in the existing works. In this article, we extend the work in Reference [10] and study the dynamics of a predator-prey model involving cannibalism in predators, where the interaction between prey and predator is in the form of a Holling type-II functional response. And, to make a more realistic model, the time delay is considered due to the delay in the gestation process of a predator. Furthermore, under specific parametric conditions, we describe the positivity and boundedness of solutions, as well as the existence and local stability of the equilibria. It also exhibits rich dynamics, such as the extinction of populations and occurrence of the Hopf bifurcation for both nondelayed and delayed models. The model considered in this study is given as follows:(4)dX1dT=X1r1−X1K−A1X2X1+B,dX2dT=−D+M1X2+A2X1T−τX2T−τX1T−τ+B−GX22H+X2,where X1 and X2 are the sizes of prey and predator populations at time T, with r, K, A1, A2, B, D, M1, G, and H are all positive constants; X1Φ=Θ1Φ, X2Φ=Θ2Φ, Φ=−τ,0, and Θ1,2:−τ,0⟶ℛ2, see Reference [22].This paper is organized as follows: in Section2, we provide the condition for positivity and boundedness of solutions, existence, and local stability of all positive equilibria, and also, Hopf bifurcation analysis is carried out near the interior equilibrium point for the model without time delay. The positivity and boundedness of solutions, local stability, and Hopf bifurcation analysis for the delayed model are given in Section 3. In Section 4, the numerical simulation is carried out to ensure our analytical findings and concluded in Section 5.


## 2. The Nondelayed Model

In this section, we study the existence of equilibria and its local stability for the model (4) without time delay τ=0. By using the transformation x1=X1/K, x2=X2, and t=rT, the model (4) with reduced parameters takes the form(5)dx1dt=x11−x1−a1x2x1+b,dx2dt=−d+m1x2+a2x1x2x1+b−γx22δ+x2,where a1=A1/rK, d=D/r, m1=M1/r, a2=A2/r, b=B/K, γ=G/r, δ=H, x1≥0, and x2≥0. a1, a2, b, d, m1, γ, and δ are all positive constants. In biological point of view, it is necessary to study in the biologically feasible region x1≥0 and x2≥0.


### 2.1. Positivity and Boundedness

In this section, positivity and boundedness of solutions of the proposed model (5) have been investigated.Theorem 1.
All solutions of model (5) are non-negative.Proof.
Sincex10≥0 and x20≥0 then from model (5), we have(6)x1t=x10exp∫0t1−x1s−a1x2sx1s+bds≥0,x2t=x20exp∫0t−d+m1+a2x1sx1s+b−γx2sδ+x2sds≥0.
Hence, all solutions of model (5) are non-negative.Theorem 2.
LetA be the set defined by(7)A=x1,x2∈ℛ+2:0≤x1≤1,0≤x1+a1a2x2≤1+14d−m1.
Then,(i)
Ais positive invariant.(ii)
All non-negative solutions of (5) are uniformly bounded forward in time and eventually enter the set A.(iii)
Model (5) is dissipative.Proof.
From the first equation of model (5), we have(8)dx1dt=x11−x1−a1x2x1+b,dx1dt≤x11−x1.
Then, solving above equation, it is bounded that(9)x1t≤cc+e−t,where c is a constant. Then, using the concept of differential inequality [23], as t⟶∞, we have(10)0≤x1t≤1.
Now, we define a function(11)σt=x1t+a1a2x2t.
Then, differentiating the above equation with respect tot, we have(12)dσdt=dx1dt+a1a2dx2dt=x11−x1−a1x1x2x1+b+−d+m1a1a2x2+a1x1x2x1+b−γa1x22a2δ+x2≤x11−x1−d−m1a1a2x2dσdt+d−m1σt≤x11−x1+d−m1x1≤d−m1+14.
Since, inA, 0≤x1≤1 and max0,1x11−x1=1/4. Using Lemma 2 as in [24], we get, for all t≥T¯≥0,(13)σt≤1+14d−m1−1+14d−m1−σT¯e−d−m1t−T¯.
Then, ifT¯=0,(14)σt≤1+14d−m1−1+14d−m1−x10+a1a2x20e−d−m1t.
Hence, sincex10,x20∈A,(15)x1t+a1a2x2t≤1+14d−m1for allt≥0.
It follows from the above results. Since solutions of the initial value problemdx1/dt=x11−x1, x10≥0 satisfies limt⟶+∞supx1t≤1.
Letε>0 be given. Then, there exists a T1>0 such that x1t≤1+ε/2 for all t≥T1. From () with T¯=T1, we get for all t≥T1≥0.(16)σt=x1t+a1a2x2t≤1+14d−m1−1+14d−m1−x1T1+a1a2x2T1ed−m1T1e−d−m1t.
Then,(17)x1t+a1a2x2t≤1+14d−m1+ε2−1+14d−m1+ε2−xT1+a1a2x2T1ed−m1T1e−d−m1t.
For all,t≥T1. Let T2≥T1 be such that(18)1+14d−m1+ε2−xT1+a1a2x2T1ed−m1T1e−d−m1t≤ε2.
For all,t≥T2. Then,(19)x1t+a1a2x2t≤1+14d−m1+εfor allt≥T2.
Hence,(20)limt⟶+∞supx1t+a1a2x2t≤1+14d−m1.Model (5) is then obviously dissipative in ℛ+2.


### 2.2. Existence of Equilibria

We find all the positive equilibrium points of the model (5) by solving the following equations:(21)1−x1−a1x2x1+b=0,(22)−d+m1+a2x1x1+b−γx2δ+x2=0.Then, we have the following equilibria for the model (5).(i)
The trivial equilibrium pointE00,0.(ii)
The first axial equilibrium pointE11,0.(iii)
The second axial equilibrium pointE20,x2¯, where x2¯=−−d+m1δ/−d+m1−γ.(iv)
The interior equilibrium pointE∗x1∗,x2∗, where x2∗=−x1∗−1b+x1∗/a1 and x1∗ is the positive root of the following equation.(23)α1x1∗3+α2x1∗2+α3x1∗+α4=0,where(24)α1=−−d+m1−a2+γ,α2=−d+m11−b−b−d+m1+a21−b−γ1−b+γb,α3=a1δ−d+m1+b−d+m1+b1−b−d+m1+a1a2δ+ba2−γb−γb1−b,α4=a1δb−d+m1+b2−d+m1−γb2.By using Descartes’ rule of sign change, we can say the number of positive roots of (23). Note it is difficult to say about all possible positive roots of the (23) analytically. Hence, the existence of coexisting equilibrium E∗x1∗,x2∗ can be discussed by plotting the nullclines; that is, the intersection points of the prey and predator nullclines are the interior equilibrium points of the model (5).Lemma 1.
For the model (5), E00,0 and E11,0 are always exists. E20,x2¯ is exist, if −d+m1>0 and −d+m1<γ. The number of sign changes in αi,i=1,2,3,4, determines the positive roots of (23). Hence, if x1∗<1, then assume that E∗x1∗,x2∗ be the arbitrary interior equilibrium exist for the model (5).
In order to have a quick glance on existence of equilibrium points, leta1=0.5,b=0.4,d=0.2,a2=0.5,δ=0.2, m1=0.3>d, and γ=0.2, we have E00,0, E11,0, and E20,0.2, and from (23), we have(25)−0.4x1∗3+0.28x1∗2+0.196x1∗−0.012=0.
Since two sign change occurs in the above equation, we have two positive roots, that is,x1∗=0.0569659<1and x1∗=1.11525>1. Then, from Lemma 1, we have the interior equilibrium E∗x1∗,x2∗=0.056965,0.861869for the model (5), which is shown in Figure 1(a). Similarly, if m1=0.2=dand m1=0.1<d, then there is only E0, E1, and E∗, which exist. Moreover, we can say if the death rate of predator d is less than the birth rate due to cannibalism m1, then the predator-only equilibrium E2 can exist. Next, we show the equilibria of the model (5) with the effect of cannibalism graphically. For the different m1 and γ, we plot the nullcline plots in Figures 1(a)–1(c).Figure 1
Nullcline plots for different values ofm1 and γ, where black line represents the prey nullcline and the colored lines represent the predator nullclines. Also, cyan dot implies E0, black dot implies E1, colored dots implies E2, and cyan boxed dots implies E∗. The parameter values are chosen as a1=0.5,b=0.4,d=0.2,a2=0.5,δ=0.2, and m1=0.3>d.
(a)(b)(c)


### 2.3. Local Stability

In order to study the local stability properties of the equilibria, we use the following Jacobian matrix at some arbitrary interior equilibriumEx1,x2, which is(26)J=1−2x1−ba1x2x1+b2−a1x1x1+bba2x2x1+b2−d+m1+a2x1x1+b−2γδx2+γx22δ+x22.The eigenvalues of the Jacobian matrices are calculated at each equilibria, in order to say the local stability properties. Then, we have the following results:Theorem 3.
For model (5),(i)
E0 is unstable if m1<d and saddle if m1>d.(ii)
E1 is stable if m1+a2/1+b<d and saddle if m1+a2/1+b>d.(iii)
E2 is stable only if 1<a1x2¯/b and −d+m1<2γδx2¯+γx2¯2/δ+x2¯2, otherwise unstable.Proof.
The Jacobian matrices atE00,0 and E11,0 are(27)JE0=100−d+m1,JE1=−1−a11+b0−d+m1+a21+b.
The eigenvalues ofJE0 are λ1=1 and λ2=−d+m1. The eigenvalues of JE1 are λ1=−1 and λ2=−d+m1+a2/1+b.
The Jacobian matrix atE20,x2¯(28)JE2=1−a1x2¯b0a2x2¯b−d+m1−2γδx2¯+γx2¯2δ+x2¯2.
The eigenvalues ofJE2 are λ1=1−a1x2¯/b and λ2=−d+m1−2γδx2¯+γx2¯2/δ+x2¯2.Theorem 4.
For model (5), if(29)γ∗<γ,ba1a2x1∗x2∗x1∗+b3<−x1∗+a1x1∗x2∗x1∗+b2δγx2∗δ+x2∗2,then the interior equilibrium E∗x1∗,x2∗ is locally asymptotically stable.Proof.
The Jacobian matrix atE∗x1∗,x2∗(30)JE∗=−x1∗+a1x1∗x2∗x1∗+b2−a1x1∗x1∗+bba2x2∗x1∗+b2−δγx2∗δ+x2∗2.
The characteristic equation of the above matrix is(31)λ2−β1λ+β2=0,where(32)β1=−x1∗+a1x1∗x2∗x1∗+b2−δγx2∗δ+x2∗2=−x1∗+a1x1∗x2∗x1∗+b2δ+x2∗2δx2∗−γ=γ∗−γ,β2=−x1∗+a1x1∗x2∗x1∗+b2−δγx2∗δ+x2∗2+ba1a2x1∗x2∗x1∗+b3.
By Routh–Hurwitz criterion, the roots of (31) has negative real parts if β1<0 and β2>0. Hence, if (29) holds, we achieve β1<0 and β2>0.


### 2.4. Hopf Bifurcation

Theorem 5.
Assume thatγ=γ∗ and(33)ba1a2x1∗x2∗x1∗+b3<−x1∗+a1x1∗x2∗x1∗+b2δγx2∗δ+x2∗2.
Then, the model (5) undergoes Hopf bifurcation near E∗x1∗,x2∗.Proof.
It is known that if (i)β1=0 and (ii) β2>0, then the roots of (31) have a pair of imaginary roots. Since we aim to study the effect of cannibalism, we take γ as a bifurcation parameter. When the parameter γ crosses the critical value γ=γc, (i) and (ii) hold if it satisfies γ=γ∗ and (33), then the roots of (31) has a pair of imaginary roots. Moreover, the transversality condition is given by (iii) d/dγβ1=−δx2∗/δ+x2∗2≠0. These conditions ensure that the model (5) undergoes Hopf bifurcation at E∗x1∗,x2∗. Note E∗x1∗,x2∗ is in terms of γ, and it is calculated at the critical value γ=γc.


## 2.1. Positivity and Boundedness

In this section, positivity and boundedness of solutions of the proposed model (5) have been investigated.Theorem 1.
All solutions of model (5) are non-negative.Proof.
Sincex10≥0 and x20≥0 then from model (5), we have(6)x1t=x10exp∫0t1−x1s−a1x2sx1s+bds≥0,x2t=x20exp∫0t−d+m1+a2x1sx1s+b−γx2sδ+x2sds≥0.
Hence, all solutions of model (5) are non-negative.Theorem 2.
LetA be the set defined by(7)A=x1,x2∈ℛ+2:0≤x1≤1,0≤x1+a1a2x2≤1+14d−m1.
Then,(i)
Ais positive invariant.(ii)
All non-negative solutions of (5) are uniformly bounded forward in time and eventually enter the set A.(iii)
Model (5) is dissipative.Proof.
From the first equation of model (5), we have(8)dx1dt=x11−x1−a1x2x1+b,dx1dt≤x11−x1.
Then, solving above equation, it is bounded that(9)x1t≤cc+e−t,where c is a constant. Then, using the concept of differential inequality [23], as t⟶∞, we have(10)0≤x1t≤1.
Now, we define a function(11)σt=x1t+a1a2x2t.
Then, differentiating the above equation with respect tot, we have(12)dσdt=dx1dt+a1a2dx2dt=x11−x1−a1x1x2x1+b+−d+m1a1a2x2+a1x1x2x1+b−γa1x22a2δ+x2≤x11−x1−d−m1a1a2x2dσdt+d−m1σt≤x11−x1+d−m1x1≤d−m1+14.
Since, inA, 0≤x1≤1 and max0,1x11−x1=1/4. Using Lemma 2 as in [24], we get, for all t≥T¯≥0,(13)σt≤1+14d−m1−1+14d−m1−σT¯e−d−m1t−T¯.
Then, ifT¯=0,(14)σt≤1+14d−m1−1+14d−m1−x10+a1a2x20e−d−m1t.
Hence, sincex10,x20∈A,(15)x1t+a1a2x2t≤1+14d−m1for allt≥0.
It follows from the above results. Since solutions of the initial value problemdx1/dt=x11−x1, x10≥0 satisfies limt⟶+∞supx1t≤1.
Letε>0 be given. Then, there exists a T1>0 such that x1t≤1+ε/2 for all t≥T1. From () with T¯=T1, we get for all t≥T1≥0.(16)σt=x1t+a1a2x2t≤1+14d−m1−1+14d−m1−x1T1+a1a2x2T1ed−m1T1e−d−m1t.
Then,(17)x1t+a1a2x2t≤1+14d−m1+ε2−1+14d−m1+ε2−xT1+a1a2x2T1ed−m1T1e−d−m1t.
For all,t≥T1. Let T2≥T1 be such that(18)1+14d−m1+ε2−xT1+a1a2x2T1ed−m1T1e−d−m1t≤ε2.
For all,t≥T2. Then,(19)x1t+a1a2x2t≤1+14d−m1+εfor allt≥T2.
Hence,(20)limt⟶+∞supx1t+a1a2x2t≤1+14d−m1.Model (5) is then obviously dissipative in ℛ+2.


## 2.2. Existence of Equilibria

We find all the positive equilibrium points of the model (5) by solving the following equations:(21)1−x1−a1x2x1+b=0,(22)−d+m1+a2x1x1+b−γx2δ+x2=0.Then, we have the following equilibria for the model (5).(i)
The trivial equilibrium pointE00,0.(ii)
The first axial equilibrium pointE11,0.(iii)
The second axial equilibrium pointE20,x2¯, where x2¯=−−d+m1δ/−d+m1−γ.(iv)
The interior equilibrium pointE∗x1∗,x2∗, where x2∗=−x1∗−1b+x1∗/a1 and x1∗ is the positive root of the following equation.(23)α1x1∗3+α2x1∗2+α3x1∗+α4=0,where(24)α1=−−d+m1−a2+γ,α2=−d+m11−b−b−d+m1+a21−b−γ1−b+γb,α3=a1δ−d+m1+b−d+m1+b1−b−d+m1+a1a2δ+ba2−γb−γb1−b,α4=a1δb−d+m1+b2−d+m1−γb2.By using Descartes’ rule of sign change, we can say the number of positive roots of (23). Note it is difficult to say about all possible positive roots of the (23) analytically. Hence, the existence of coexisting equilibrium E∗x1∗,x2∗ can be discussed by plotting the nullclines; that is, the intersection points of the prey and predator nullclines are the interior equilibrium points of the model (5).Lemma 1.
For the model (5), E00,0 and E11,0 are always exists. E20,x2¯ is exist, if −d+m1>0 and −d+m1<γ. The number of sign changes in αi,i=1,2,3,4, determines the positive roots of (23). Hence, if x1∗<1, then assume that E∗x1∗,x2∗ be the arbitrary interior equilibrium exist for the model (5).
In order to have a quick glance on existence of equilibrium points, leta1=0.5,b=0.4,d=0.2,a2=0.5,δ=0.2, m1=0.3>d, and γ=0.2, we have E00,0, E11,0, and E20,0.2, and from (23), we have(25)−0.4x1∗3+0.28x1∗2+0.196x1∗−0.012=0.
Since two sign change occurs in the above equation, we have two positive roots, that is,x1∗=0.0569659<1and x1∗=1.11525>1. Then, from Lemma 1, we have the interior equilibrium E∗x1∗,x2∗=0.056965,0.861869for the model (5), which is shown in Figure 1(a). Similarly, if m1=0.2=dand m1=0.1<d, then there is only E0, E1, and E∗, which exist. Moreover, we can say if the death rate of predator d is less than the birth rate due to cannibalism m1, then the predator-only equilibrium E2 can exist. Next, we show the equilibria of the model (5) with the effect of cannibalism graphically. For the different m1 and γ, we plot the nullcline plots in Figures 1(a)–1(c).Figure 1
Nullcline plots for different values ofm1 and γ, where black line represents the prey nullcline and the colored lines represent the predator nullclines. Also, cyan dot implies E0, black dot implies E1, colored dots implies E2, and cyan boxed dots implies E∗. The parameter values are chosen as a1=0.5,b=0.4,d=0.2,a2=0.5,δ=0.2, and m1=0.3>d.
(a)(b)(c)


## 2.3. Local Stability

In order to study the local stability properties of the equilibria, we use the following Jacobian matrix at some arbitrary interior equilibriumEx1,x2, which is(26)J=1−2x1−ba1x2x1+b2−a1x1x1+bba2x2x1+b2−d+m1+a2x1x1+b−2γδx2+γx22δ+x22.The eigenvalues of the Jacobian matrices are calculated at each equilibria, in order to say the local stability properties. Then, we have the following results:Theorem 3.
For model (5),(i)
E0 is unstable if m1<d and saddle if m1>d.(ii)
E1 is stable if m1+a2/1+b<d and saddle if m1+a2/1+b>d.(iii)
E2 is stable only if 1<a1x2¯/b and −d+m1<2γδx2¯+γx2¯2/δ+x2¯2, otherwise unstable.Proof.
The Jacobian matrices atE00,0 and E11,0 are(27)JE0=100−d+m1,JE1=−1−a11+b0−d+m1+a21+b.
The eigenvalues ofJE0 are λ1=1 and λ2=−d+m1. The eigenvalues of JE1 are λ1=−1 and λ2=−d+m1+a2/1+b.
The Jacobian matrix atE20,x2¯(28)JE2=1−a1x2¯b0a2x2¯b−d+m1−2γδx2¯+γx2¯2δ+x2¯2.
The eigenvalues ofJE2 are λ1=1−a1x2¯/b and λ2=−d+m1−2γδx2¯+γx2¯2/δ+x2¯2.Theorem 4.
For model (5), if(29)γ∗<γ,ba1a2x1∗x2∗x1∗+b3<−x1∗+a1x1∗x2∗x1∗+b2δγx2∗δ+x2∗2,then the interior equilibrium E∗x1∗,x2∗ is locally asymptotically stable.Proof.
The Jacobian matrix atE∗x1∗,x2∗(30)JE∗=−x1∗+a1x1∗x2∗x1∗+b2−a1x1∗x1∗+bba2x2∗x1∗+b2−δγx2∗δ+x2∗2.
The characteristic equation of the above matrix is(31)λ2−β1λ+β2=0,where(32)β1=−x1∗+a1x1∗x2∗x1∗+b2−δγx2∗δ+x2∗2=−x1∗+a1x1∗x2∗x1∗+b2δ+x2∗2δx2∗−γ=γ∗−γ,β2=−x1∗+a1x1∗x2∗x1∗+b2−δγx2∗δ+x2∗2+ba1a2x1∗x2∗x1∗+b3.
By Routh–Hurwitz criterion, the roots of (31) has negative real parts if β1<0 and β2>0. Hence, if (29) holds, we achieve β1<0 and β2>0.


## 2.4. Hopf Bifurcation

Theorem 5.
Assume thatγ=γ∗ and(33)ba1a2x1∗x2∗x1∗+b3<−x1∗+a1x1∗x2∗x1∗+b2δγx2∗δ+x2∗2.
Then, the model (5) undergoes Hopf bifurcation near E∗x1∗,x2∗.Proof.
It is known that if (i)β1=0 and (ii) β2>0, then the roots of (31) have a pair of imaginary roots. Since we aim to study the effect of cannibalism, we take γ as a bifurcation parameter. When the parameter γ crosses the critical value γ=γc, (i) and (ii) hold if it satisfies γ=γ∗ and (33), then the roots of (31) has a pair of imaginary roots. Moreover, the transversality condition is given by (iii) d/dγβ1=−δx2∗/δ+x2∗2≠0. These conditions ensure that the model (5) undergoes Hopf bifurcation at E∗x1∗,x2∗. Note E∗x1∗,x2∗ is in terms of γ, and it is calculated at the critical value γ=γc.


## 3. The Delayed Model

The time delay in the predator-prey model can cause complex behavior in the dynamics. For example, the authors in [3] considered gestation delay in the food chain model with Crowley–Martin functional response. They showed that the presence of delay helps to stabilize the unstable near the interior equilibrium point to stable. In Reference [25], the authors considered the three species food chain model with the interaction between the species in the form of Holling type-II functional response. Moreover, they considered time delay in the gestation process of the top predator. They showed that the presence of delay exhibits chaos in the considered model with the help of the bifurcation diagram and maximum Lyapunov exponents. In Reference [26], the authors considered the spatiotemporal prey-predator model with additive Allee effect in prey growth, Holling type-II functional response, and gestation delay in predator population. With the increment in time delay, the stationary pattern gets converted into another one which eventually turns into a chaotic pattern for the sufficiently large time delay. Also, they showed that the transition where cold spot pattern turns into a stationary mixture pattern, and finally, the mixture pattern eventually settles into a chaotic pattern through the quasiperiodic one with the increase in the magnitude of time delay. Two interesting scenarios for the temporal model correspond to the spatiotemporal model, where the bistable scenario for an intermediate range of parameter values is chosen as bifurcation parameter [27]. Based upon these two bifurcation scenarios, the authors are interested in understanding the role of time delay on spatiotemporal pattern formation. For other interesting results on time delay, we refer the readers to [16, 19, 28]. In this section, we consider the model (4) in presence delay τ≠0 with reduced parameters, which is of the form(34)dx1dt=x11−x1−a1x2x1+b,dx2dt=−d+m1x2+a2x1t−τx2t−τx1t−τ+b−γx22δ+x2.Withx1ϕ=θ1ϕ, x2ϕ=θ2ϕ, ϕ=−τ,0, and θ1,2:−τ,0⟶ℛ2, we note the equilibrium points for the models (34) and (5) are same. Moreover, in the following, we analyze the local stability and Hopf bifurcation near the interior equilibrium point E∗x1∗,x2∗.


### 3.1. Positivity and Boundedness

One can write from the first equation of model (34) as follows:(35)dx1x1=1−x1−a1x2x1+b.Integrating between 0 andt, we get(36)x1t=x10exp∫0t1−x1−a1x2x1+bds≥0.Similarly from second equation of model (34), we have(37)x2t=x20exp∫0t−d+m1+a2x1s−τx1s−τ+b−γx2sδ+x2sds≥0.Hence, all the solutions of model (34) are non-negative.Theorem 6.
All solutions of model (34) starting in R+2 are confined to the region D∗=x1,x2∈ℛ+2:ζt≤1+a11−d+m2/4γa2 as t⟶∞ for all x1ϕ,x2ϕ∈ℛ+2, where ζt=x1t−τ+a1/a2x2t.Proof.
define a functionζt=x1t−τ+a1/a2x2t. Then, we obtain(38)dζtdt=dx1t−τdt+a1a2dx2tdt=x1t−τ1−x1t−τ−a1x1t−τx2t−τx1t−tau+ba1a2−d+m1x2+a1x1t−τx2t−τx1t−τ+b−a1a2γx22tδ+x2t=−x1t−τ+a1a2x2+2x1t−τ−x12t−τ+a1γa2x2t−d+m+1γ−x2tδ+x2t,dζtdt+ζt=2x1t−τ−x12t−τ+a1γa2x2t−d+m+1γ−x2tδ+x2t≤1+a11−d+m24γa2.
Using Lemma 2 as in [24], thus all the solutions of model (34) are bounded.


### 3.2. Local Stability and Hopf Bifurcation

The delayed model (34) after linearization using the transformation x1=x1−x1∗ and x2=x2−x2∗ is given as follows:(39)x˙1x˙2=−x1∗+a1x1∗x2∗x1∗+b2−a1x1∗x1∗+b0−d+m1−2γδx2∗−γx2∗2δ+x2∗2x1x2+00ba2x2∗x1∗+b2a2x1∗x1∗+bx1−τx2−τ,=ϑ1ϑ20ϑ3x1x2+00ϑ4ϑ5x1−τx2−τ.We removed the bars for our convenience, and then, the characteristic equation for (39) is given by(40)ϑ1−λϑ2ϑ4e−λτϑ3+ϑ5e−λτ−λ=0.That is,(41)λ2−ϑ1+ϑ3λ+ϑ1ϑ3+−ϑ5λ+ϑ1ϑ5−ϑ2ϑ4e−λτ=0.Substituteλ=iϕ in (41) and separate the real and imaginary parts:(42)ϑ2ϑ4−ϑ1ϑ5sinϕτ−ϕϑ5cosϕτ=ϕϑ1+ϑ3,ϑ5ϕsinϕτ+ϑ2ϑ4−ϑ1ϑ5cosϕτ=−ϕ2+ϑ1ϑ3.After simplification, we get(43)sinϕτ=ϕϑ1+ϑ3ϑ2ϑ4−ϑ1ϑ5+ϑ5ϕ−ϕ2+ϑ1ϑ3ϕ2ϑ52+ϑ2ϑ4−ϑ1ϑ52,cosϕτ=−ϕ2+ϑ1ϑ3ϑ2ϑ4−ϑ1ϑ5−ϑ5ϕ2ϑ1+ϑ3ϕ2ϑ52+ϑ2ϑ4−ϑ1ϑ52.Sincesin2ϕτ+cos2ϕτ=1, we have(44)φ1ϕ6+φ2ϕ4+φ3ϕ2+φ4=0,where(45)φ1=ϑ52,φ2=−2ϑ1ϑ3ϑ52+ϑ2ϑ4−ϑ1ϑ52+ϑ52ϑ1+ϑ32−ϑ54,φ3=ϑ2ϑ4−ϑ1ϑ52ϑ1+ϑ32+ϑ12ϑ32ϑ52−2ϑ2ϑ4−ϑ1ϑ52ϑ1ϑ3−2ϑ2ϑ4−ϑ1ϑ52ϑ52,φ4=ϑ12ϑ32ϑ2ϑ4−ϑ1ϑ52−ϑ2ϑ4−ϑ1ϑ54.Let us assume that (44) has at least one positive root, which is ϕ∗; then, we get the critical value τn as follows:(46)τn=1ϕ∗arcsinϕ∗ϑ1+ϑ3ϑ2ϑ4−ϑ1ϑ5+ϑ5ϕ∗−ϕ∗2+ϑ1ϑ3ϕ∗2ϑ52+ϑ2ϑ4−ϑ1ϑ52+2nπϕ∗,n=0,1,2,….Theorem 7.
The following transversality condition is holds:(47)dRλϕdτ|τ=τ∗≠0.Proof.
Let us substituteλτ=θτ+iϕτ in () and differentiate wrtτ, we get(48)dλdτ−1=−2λ+ϑ1+ϑ3eλτ+ϑ5λϑ5λ−ϑ1ϑ5−ϑ2ϑ4−τλ.
Then,(49)Rdλdτ−1|λ=iϕ∗=−2λ+ϑ1+ϑ3eλτ+ϑ5λϑ5λ−ϑ1ϑ5−ϑ2ϑ4|λ=iϕ=1ϑ52ϕ∗3+ϕ∗ϑ2ϑ4−ϑ1ϑ52−ϑ52+ϑ1ϑ5−ϑ2ϑ4ϑ1+ϑ3−2ϑ5ϕ∗2sinϕ∗τ−ϑ5ϕ∗ϑ1+ϑ3−2ϕ∗ϑ1ϑ5−ϑ2ϑ4cosϕ∗τ.Theorem 8.
If (29) holds then for the model (34), we have the following: (i) the interior equilibrium point E∗x1∗,x2∗ is locally asymptotically stable for τ∈0,τ∗ and (ii) undergoes Hopf bifurcation near E∗x1∗,x2∗ at τ=τn, n=0,1,2,…. Also, each τn satisfies (47), which is the critical points that were the stability switches of E∗ that occurs.


## 3.1. Positivity and Boundedness

One can write from the first equation of model (34) as follows:(35)dx1x1=1−x1−a1x2x1+b.Integrating between 0 andt, we get(36)x1t=x10exp∫0t1−x1−a1x2x1+bds≥0.Similarly from second equation of model (34), we have(37)x2t=x20exp∫0t−d+m1+a2x1s−τx1s−τ+b−γx2sδ+x2sds≥0.Hence, all the solutions of model (34) are non-negative.Theorem 6.
All solutions of model (34) starting in R+2 are confined to the region D∗=x1,x2∈ℛ+2:ζt≤1+a11−d+m2/4γa2 as t⟶∞ for all x1ϕ,x2ϕ∈ℛ+2, where ζt=x1t−τ+a1/a2x2t.Proof.
define a functionζt=x1t−τ+a1/a2x2t. Then, we obtain(38)dζtdt=dx1t−τdt+a1a2dx2tdt=x1t−τ1−x1t−τ−a1x1t−τx2t−τx1t−tau+ba1a2−d+m1x2+a1x1t−τx2t−τx1t−τ+b−a1a2γx22tδ+x2t=−x1t−τ+a1a2x2+2x1t−τ−x12t−τ+a1γa2x2t−d+m+1γ−x2tδ+x2t,dζtdt+ζt=2x1t−τ−x12t−τ+a1γa2x2t−d+m+1γ−x2tδ+x2t≤1+a11−d+m24γa2.
Using Lemma 2 as in [24], thus all the solutions of model (34) are bounded.


## 3.2. Local Stability and Hopf Bifurcation

The delayed model (34) after linearization using the transformation x1=x1−x1∗ and x2=x2−x2∗ is given as follows:(39)x˙1x˙2=−x1∗+a1x1∗x2∗x1∗+b2−a1x1∗x1∗+b0−d+m1−2γδx2∗−γx2∗2δ+x2∗2x1x2+00ba2x2∗x1∗+b2a2x1∗x1∗+bx1−τx2−τ,=ϑ1ϑ20ϑ3x1x2+00ϑ4ϑ5x1−τx2−τ.We removed the bars for our convenience, and then, the characteristic equation for (39) is given by(40)ϑ1−λϑ2ϑ4e−λτϑ3+ϑ5e−λτ−λ=0.That is,(41)λ2−ϑ1+ϑ3λ+ϑ1ϑ3+−ϑ5λ+ϑ1ϑ5−ϑ2ϑ4e−λτ=0.Substituteλ=iϕ in (41) and separate the real and imaginary parts:(42)ϑ2ϑ4−ϑ1ϑ5sinϕτ−ϕϑ5cosϕτ=ϕϑ1+ϑ3,ϑ5ϕsinϕτ+ϑ2ϑ4−ϑ1ϑ5cosϕτ=−ϕ2+ϑ1ϑ3.After simplification, we get(43)sinϕτ=ϕϑ1+ϑ3ϑ2ϑ4−ϑ1ϑ5+ϑ5ϕ−ϕ2+ϑ1ϑ3ϕ2ϑ52+ϑ2ϑ4−ϑ1ϑ52,cosϕτ=−ϕ2+ϑ1ϑ3ϑ2ϑ4−ϑ1ϑ5−ϑ5ϕ2ϑ1+ϑ3ϕ2ϑ52+ϑ2ϑ4−ϑ1ϑ52.Sincesin2ϕτ+cos2ϕτ=1, we have(44)φ1ϕ6+φ2ϕ4+φ3ϕ2+φ4=0,where(45)φ1=ϑ52,φ2=−2ϑ1ϑ3ϑ52+ϑ2ϑ4−ϑ1ϑ52+ϑ52ϑ1+ϑ32−ϑ54,φ3=ϑ2ϑ4−ϑ1ϑ52ϑ1+ϑ32+ϑ12ϑ32ϑ52−2ϑ2ϑ4−ϑ1ϑ52ϑ1ϑ3−2ϑ2ϑ4−ϑ1ϑ52ϑ52,φ4=ϑ12ϑ32ϑ2ϑ4−ϑ1ϑ52−ϑ2ϑ4−ϑ1ϑ54.Let us assume that (44) has at least one positive root, which is ϕ∗; then, we get the critical value τn as follows:(46)τn=1ϕ∗arcsinϕ∗ϑ1+ϑ3ϑ2ϑ4−ϑ1ϑ5+ϑ5ϕ∗−ϕ∗2+ϑ1ϑ3ϕ∗2ϑ52+ϑ2ϑ4−ϑ1ϑ52+2nπϕ∗,n=0,1,2,….Theorem 7.
The following transversality condition is holds:(47)dRλϕdτ|τ=τ∗≠0.Proof.
Let us substituteλτ=θτ+iϕτ in () and differentiate wrtτ, we get(48)dλdτ−1=−2λ+ϑ1+ϑ3eλτ+ϑ5λϑ5λ−ϑ1ϑ5−ϑ2ϑ4−τλ.
Then,(49)Rdλdτ−1|λ=iϕ∗=−2λ+ϑ1+ϑ3eλτ+ϑ5λϑ5λ−ϑ1ϑ5−ϑ2ϑ4|λ=iϕ=1ϑ52ϕ∗3+ϕ∗ϑ2ϑ4−ϑ1ϑ52−ϑ52+ϑ1ϑ5−ϑ2ϑ4ϑ1+ϑ3−2ϑ5ϕ∗2sinϕ∗τ−ϑ5ϕ∗ϑ1+ϑ3−2ϕ∗ϑ1ϑ5−ϑ2ϑ4cosϕ∗τ.Theorem 8.
If (29) holds then for the model (34), we have the following: (i) the interior equilibrium point E∗x1∗,x2∗ is locally asymptotically stable for τ∈0,τ∗ and (ii) undergoes Hopf bifurcation near E∗x1∗,x2∗ at τ=τn, n=0,1,2,…. Also, each τn satisfies (47), which is the critical points that were the stability switches of E∗ that occurs.


## 4. Numerical Simulation

In this section, we perform some simulation results to show the local stability and bifurcation behavior of both nondelayed and delayed models.Let us take the fixed parameter values as(50)a1=0.5,b=0.4,d=0.1,a2=0.5.Case 1.
The nondelayed model is as follows:
With the parameter values in (50), the model (5) becomes(51)dx1dt=x11−x1−0.5x2x1+0.4,dx2dt=−0.1+m1x2+0.5x1x2x1+0.4−γx22δ+x2.
First, we showed that form1=0.3>d, γ=0.2, and δ=0.2, then the model (51) has the equilibrium points E00,0, E11,0, E20,0.2, and E∗0.056965,0.0861869. Further, the existence of equilibria is shown with the help of plotting the nullcline plots in Figures 1(a)–1(c). The predator-only equilibrium E2 is exist only if m1>d, which is shown in Figure 1(a). If E2 is stable, then the predator may alone survive without prey by getting energy from its own species m1>d; that is, a death rate of the predator is greater than birth due to cannibalism. If E2 is unstable, then the trajectories near E2 either approach E0 or E∗. From Theorem 3, we have E0 is a saddle point, E1 is also a saddle point, since m1+a2/1+b=0.657143>d. And, a1x2¯/b<1, then E2 is unstable for the model (51). Also, the predator-only equilibrium E2 does not exist for m1≤d (see Figures 1(b) and 1(c)). The local stability of E∗ is achieved in Theorem 4, which ensures the long-time survival of both species. The model (51) at γ=0.36, d=0.1, and m1=0.2 with other parameters in (50) is locally asymptotically stable near the interior equilibrium E∗, see Figure 2. Furthermore, if γ=γ∗=0.34324, also, then the model (51) has interior equilibrium E∗x1∗,x2∗=0.234183,0.971336. Then, the following conditions hold: (i) β1=0; (ii) β2=0.0868205; (iii) d/dγβ1=−0.141591≠0. Then, from Theorem 5, the model (51) undergoes Hopf bifurcation at γ=γ∗=0.34324 at critical value γc=0.34324 and is shown in Figure 3 for γ=0.34. Further, on increasing δ=0.8, then the model (51) is locally asymptotically stable (see Figure 4). For clear representation, one-parameter bifurcation diagram with respect to γ∈0.2,0.5 and δ∈0,0.8 is shown in Figures 5 and 6, correspondingly. In Figure 5, it is shown that for γ∈0,0.34324, the periodic solution arises near the E∗, and for γ∈0.34324,0.5, the E∗ is locally asymptotically stable. Also in Figure 5, we noticed on increasing the γ value, the size of prey population increases, and the predator population size declines. Figure 6 shows that E∗ is locally asymptotically stable for δ∈0,0.1905, the periodic solution arises near the E∗ for δ∈0.1905,0.593, and the E∗ is locally asymptotically stable γ∈0.593,0.8. Hence, we notice from Figure 6, that on increasing the values of δ, both population sizes reduce. Moreover, the existence of exchange of stability (unstable to stable) on varying the parameters γ and δ occurs; hence, we can conclude that the parameters γ and δ has stabilizing effect in the nondelayed model (5).
Next, to demonstrate the impact of cannibalism parametersγ, δ, and m1 on the model (51), we plotted the two-parameter bifurcation diagram with respect to the parameters δ and γ in Figure 7(a), and m1 and δ in Figure 7(b). It is shown that less γ value causes the extinction of species. Also, the stable, unstable, and extinction regions are clearly depicted. Further, we showed that the γ value should be above some threshold value in order to ensure the survival of the species. Similarly, m1 should be less than some critical value for the existence of species, or else the species will die out. So, cannibalism in the predator population can cause complex behavior in the model (51).Figure 2
Locally asymptotically stable. (a) Time plot; (b) phase portrait for the model (51) with m1=0.2, δ=0.2, and γ=0.36.
(a)(b)Figure 3
Periodic oscillations. (a) Time plot; (b) phase portrait for the model (51) with m1=0.2, δ=0.2, and γ=0.34.
(a)(b)Figure 4
Locally asymptotically stable. (a) Time plot; (b) phase portrait for the model (51) with m1=0.2, δ=0.2, and γ=0.34.
(a)(b)Figure 5
The bifurcation diagram of the model (51) with m1=0.2, δ=0.2, and γ∈0.2,0.5.
(a)(b)Figure 6
The bifurcation diagram of the model (51) with m1=0.2, γ=0.34, and δ∈0,0.8.
(a)(b)Figure 7
(a) The two-parameter bifurcation diagram of the model (51) with m1=0.2, γ∈0,0.5, and δ∈0,0.8. (b) The two-parameter bifurcation diagram of the model with γ=0.34, m1∈0,0.5, and δ∈0,0.8. The cyan color line separates the stable and extinction region. The blue color line (Hopf bifurcation curve) separates the stable and unstable region.
(a)(b)Case 2.
The delayed model.
Letd=0.1, m1=0.2, γ=0.4, δ=0.4, and other parameter values in (50); then, the model (34) becomes(52)dx1dt=x11−x1−0.5x2x1+0.4,dx2dt=−0.1+m1x2+0.5x1t−τx2t−τx1t−τ+0.4−γx22delta+x2.
Then, the model (52) has interior equilibrium point E∗x1∗,x2∗=0.346151,0.97574. Moreover, from (44), we have(53)0.0538044ϕ6+0.00677379ϕ4−0.0001068487ϕ2−0.0000251728=0.
If it has a positive rootϕ∗=0.2374849, then we get the critical value τ=1.69709, where the stability switches occur for the model (52). The local stability and Hopf bifurcation results are given in Theorem 8 that, for τ∈0,1.69709, model (52) is locally asymptotically stable (see Figure 8) and undergoes Hopf bifurcation at τ∗=1.69709 (see Figure 9). For clear presentation, the bifurcation diagram with respect to τ for the existence of Hopf bifurcation is shown in Figure 10. Hence, the time delay in the model (34) has a destabilizing effect. To show the effectiveness of delay parameter τ with the cannibalism parameters γ and δ, we plotted the two-parameter bifurcation diagram with respect to δ and τ in Figure 11(a), and γ and τ in Figure 11(b). The exchange of stability regions (stable to unstable) for the model (52) is clearly depicted. It shows that on increasing the δ value, the interval of stability with respect to τ decreases. For instance, if δ=0.2, the interior equilibrium E∗ for the model (52) is locally asymptotically stable for τ∈0,1.69704 and Hopf bifurcation occurs when τ crosses the critical value τ=1.69704, and if δ=0.5, the interior equilibrium E∗ for the model (52) is locally asymptotically stable for τ∈0,0.20826 and Hopf bifurcation occurs when τ crosses the critical value τ=0.20826. Similarly, for larger values of γ, the interval for local stability increases for τ.Figure 8
Locally asymptotically stable. (a) Time plot; (b) phase portrait for the model (52) with m1=0.2, γ=0.2, δ=0.4, and τ=1.3.
(a)(b)Figure 9
Periodic oscillations. (a) Time plot; (b) phase portrait for the model (52) with m1=0.2, γ=0.2, δ=0.4, and τ=2.
(a)(b)Figure 10
The bifurcation diagram of the model (52) with m1=0.2, γ=0.2, δ=0.4, and τ∈1.6,2.4.
(a)(b)Figure 11
(a) The two-parameter bifurcation diagram of the model (52) with γ=0.4, τ∈0.2,1.8, and δ∈0.2,0.5. (b) The two parameter bifurcation diagram of the model (52) with δ=0.2, τ∈0,6, and γ∈0.34,0.42. The blue color line (Hopf bifurcation curve) separates the stable and unstable region.
(a)(b)


## 5. Conclusion

In this work, we considered the Holling type-II prey-predator model involving predator cannibalism, and also, the delay is considered due to the gestation process in the predator population. We were given a description of what a two-species predation model should be and how its solutions should behave. There have been relatively few attempts to suggest explicit models for cannibalism. While in the model, we have proposed both delay and cannibalism, all the results that we have deduced on the behavior of the model in terms of stability and bifurcation analysis. There is some biological evidence to suggest that complicated population systems have a tendency to be more stable than simple systems. On the other hand, the removal of one species can lead to a collapse of population systems. It is important to know what the predator population involves cannibalism. In the case when the model (4) has no cannibalism, then it follows the well-known Rosenzweig–MacArthur model [10], in the absence of its lone prey, the predator dies out exponentially. But in the presence of cannibalism, if the death of the predator is greater than that of the birth due to cannibalism, then the prey-free equilibrium exists (see Figure 1(a)). Besides, if the birth due to cannibalism is less than or equal to the death rate of a predator, then predators cannot survive alone; that is, the prey-free equilibrium E2 does not exist, which is clearly shown in Figures 1(b) and 1(c).In Reference [5], the authors showed that model (1) without cannibalism has a boundary equilibrium, and it is globally asymptotically stable. For a suitable rate of cannibalism, the model (1) has a unique interior equilibrium, and it is globally asymptotically stable. They showed that the high rate due to cannibalism causes the extinction of the prey population. Further, predator-only equilibrium exists, and it is globally asymptotically stable. The presence of cannibalism has both positive and negative effects. The cannibalism in the prey cannot stabilize the unstable interior equilibrium in the ODE case but can destabilize the stable interior equilibrium, leading to a stable limit cycle [29]. The authors in Reference [13] reported spatial patterning in two-species predator-prey models are driven solely via the joint effect of predator and prey cannibalism. Interestingly, higher levels of equilibrium prey provide stability, while lower levels drive instability. In this study, we derived the condition to undergo Hopf bifurcation for both cannibalism and delay parameters. The local stability conditions prevent both populations from extinction risk. Also, the model may have bifurcation for other model parameters, but we are particularly interested in varying the cannibalism and delay parameters. We found that model (5) in the presence of cannibalism is more stable for high cannibalism rate and unstable for a higher value of delay parameter in the model (34). Additionally, the proposed model can be studied in discrete and stochastic forms, which may result in richer dynamical features than the proposed model. This will also be our future goal.



---

*Source: 4071375-2022-09-17.xml*