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

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.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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,$$\begin{array}{}\text{(1)}& \begin{array}{c}\frac{\text{d}{X}_1}{\text{d}T}=X_{1}r1-\frac{X_1}{K}-A_1{X}_2,\hfill \\ \frac{\text{d}{X}_2}{\text{d}T}=X_2-D+M_1+A_2{X}_1-\frac{G{X}_2}{H+X_2},\hfill \end{array}\end{array}$$where $X_1$ and $X_2$ are the sizes of prey and predator population at the time $T$; $r$ and $K$ denote the growth rate and carrying capacity of prey; $A_1$ denotes the reduction rate of prey; $D$ represents the death rate of predator; $A_2$ denotes the conversion efficiency of the predator; $G\gamma X_2/H+X_2$ is the cannibalism term, $\gamma$ denotes the rate of cannibalism; $M_1<G$, $M_1$ 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,$$\begin{array}{}\text{(2)}& \begin{array}{c}\frac{\text{d}{X}_1}{\text{d}T}=Xr+M_1-\frac{r{X}_1}{K}-A_1{X}_1-\frac{G{X}_1}{H+X_1},\hfill \\ \frac{\text{d}{X}_2}{\text{d}T}=X_2-d+A_2{X}_1,\hfill \end{array}\end{array}$$where $X_1$ and $X_2$ are the size of prey and predator at time $T$. $r$, $k$, $A_1$, $G$, $H$, $D$, $M_1$, and $A_2$ 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$$\begin{array}{}\text{(3)}& \begin{array}{c}\frac{\text{d}{X}_1}{\text{d}T}=X_{1}r1-\frac{X_1}{K}-\frac{A_1{X}_2}{X+B_1},\hfill \\ \frac{\text{d}{X}_2}{\text{d}T}=X_2-d+\frac{A_2{X}_1}{X_1+B_2},\hfill \end{array}\end{array}$$where $X_1$ and $X_2$ are the sizes of prey and predator at time $T$. $B_1$ and $B_2$ are the half-saturation constants for the prey and predator. In this study, we assumed $B_1=B_2$ 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:$$\begin{array}{}\text{(4)}& \begin{array}{c}\frac{\text{d}{X}_1}{\text{d}T}=X_{1}r1-\frac{X_1}{K}-\frac{A_1{X}_2}{X_1+B},\hfill \\ \frac{\text{d}{X}_2}{\text{d}T}=-D+M_1{X}_2+\frac{A_2{X}_{1}T-\tau X_{2}T-\tau }{X_{1}T-\tau +B}-\frac{G{X}_2^2}{H+X_2},\hfill \end{array}\end{array}$$where $X_1$ and $X_2$ are the sizes of prey and predator populations at time $T$, with $r$, $K$, $A_1$, $A_2$, $B$, $D$, $M_1$, $G$, and $H$ are all positive constants; $X_1\mathrm{\Phi }={\mathrm{\Theta }}_1\mathrm{\Phi }$, $X_2\mathrm{\Phi }={\mathrm{\Theta }}_2\mathrm{\Phi }$, $\mathrm{\Phi }=-\tau ,0$, and ${\mathrm{\Theta }}_{\mathrm{1,2}}:-\tau ,0⟶{\mathrm{ℛ}}^2$, see Reference [22].

This paper is organized as follows: in Section 2, 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 $\tau =0$. By using the transformation $x_1=X_1/K$, $x_2=X_2$, and $t=rT$, the model (4) with reduced parameters takes the form$$\begin{array}{}\text{(5)}& \begin{array}{c}\frac{\text{d}{x}_1}{\text{d}t}=x_{1}1-x_1-\frac{a_1{x}_2}{x_1+b},\hfill \\ \frac{\text{d}{x}_2}{\text{d}t}=-d+m_1{x}_2+\frac{a_2{x}_1{x}_2}{x_1+b}-\frac{\gamma x_2^2}{\delta +x_2},\hfill \end{array}\end{array}$$where $a_1=A_1/rK$, $d=D/r$, $m_1=M_1/r$, $a_2=A_2/r$, $b=B/K$, $\gamma =G/r$, $\delta =H$, $x_1\ge 0$, and $x_2\ge 0$. $a_1$, $a_2$, $b$, $d$, $m_1$, $\gamma$, and $\delta$ are all positive constants. In biological point of view, it is necessary to study in the biologically feasible region $x_1\ge 0$ and $x_2\ge 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.

Since $x_{1}0\ge 0$ and $x_{2}0\ge 0$ then from model (5), we have$$\begin{array}{c}\text{(6)}& x_{1}t=x_{1}0\exp {\displaystyle {\int }_0^{t}1-x_{1}s-\frac{a_1{x}_{2}s}{x_{1}s+b}}\text{d}s\ge 0,x_{2}t=x_{2}0\exp {\displaystyle {\int }_0^t-d+m_1+\frac{a_2{x}_{1}s}{x_{1}s+b}-\frac{\gamma x_{2}s}{\delta +x_{2}s}}\text{d}s\ge 0.\hfill \end{array}$$

Hence, all solutions of model (5) are non-negative.

Theorem 2.

Let $\mathcal{A}$ be the set defined by$$\begin{array}{}\text{(7)}& \mathcal{A}=x_1,x_2\in {\mathrm{ℛ}}_{+}^2:0\le x_1\le \mathrm{1,0}\le x_1+\frac{a_1}{a_2}{x}_2\le 1+\frac{1}{4d-m_1}.\end{array}$$

Then,

(i) $\mathcal{A}$is positive invariant.

(ii) All non-negative solutions of (5) are uniformly bounded forward in time and eventually enter the set $\mathcal{A}$.

(iii) Model (5) is dissipative.

Proof.

From the first equation of model (5), we have$$\begin{array}{c}\text{(8)}& \frac{d{x}_1}{dt}=x_{1}1-x_1-\frac{a_1{x}_2}{x_1+b},\frac{d{x}_1}{dt}\le x_{1}1-x_1.\hfill \end{array}$$

Then, solving above equation, it is bounded that$$\begin{array}{}\text{(9)}& x_{1}t\le \frac{c}{c+e^{-t}},\end{array}$$where $c$ is a constant. Then, using the concept of differential inequality [23], as $t⟶\infty$, we have$$\begin{array}{}\text{(10)}& 0\le x_{1}t\le 1.\end{array}$$

Now, we define a function$$\begin{array}{}\text{(11)}& \sigma t=x_{1}t+\frac{a_1}{a_2}{x}_{2}t.\end{array}$$

Then, differentiating the above equation with respect to $t$, we have$$\begin{array}{c}\text{(12)}& \frac{\text{d}\sigma }{\text{d}t}=\frac{\text{d}{x}_1}{\text{d}t}+\frac{a_1}{a_2}\frac{\text{d}{x}_2}{\text{d}t}=x_{1}1-x_1-\frac{a_1{x}_1{x}_2}{x_1+b}+\frac{{-d+m}_1{a}_1}{a_2}{x}_2\\ +\frac{a_1{x}_1{x}_2}{x_1+b}-\frac{\gamma a_1{x}_2^2}{a_2\delta +x_2}\\ \le x_{1}1-x_1-\frac{d-m_1{a}_1}{a_2}{x}_2\frac{d\sigma }{dt}+d-m_1\sigma t\\ \le x_{1}1-x_1+d-m_1{x}_1\le d-m_1+\frac{1}{4}.\end{array}$$

Since, in $\mathcal{A}$, $0\le x_1\le 1$ and ${\max }_{\mathrm{0,1}}{x}_{1}1-x_1=1/4$. Using Lemma 2 as in [24], we get, for all $t\ge \overline{T}\ge 0$,$$\begin{array}{}\text{(13)}& \sigma t\le 1+\frac{1}{4d-m_1}-1+\frac{1}{4d-m_1}-\sigma \overline{T}{e}^{-d-m_{1}t-\overline{T}}.\end{array}$$

Then, if $\overline{T}=0$,$$\begin{array}{}\text{(14)}& \sigma t\le 1+\frac{1}{4d-m_1}-1+\frac{1}{4d-m_1}-x_{1}0+\frac{a_1}{a_2}{x}_{2}0e^{-d-m_{1}t}.\end{array}$$

Hence, since $x_{1}0,x_{2}0\in \mathcal{A}$,$$\begin{array}{}\text{(15)}& x_{1}t+\frac{a_1}{a_2}{x}_{2}t\le 1+\frac{1}{4d-m_1}\text{for\hspace{0.17em}all}\hspace{1em}t\ge 0.\end{array}$$

It follows from the above results. Since solutions of the initial value problem $\text{d}{x}_1/\text{d}t=x_{1}1-x_1$, $x_{1}0\ge 0$ satisfies ${\text{lim}}_{t⟶+\infty }\text{sup\hspace{0.17em}}{x}_{1}t\le 1$.

Let $\epsilon >0$ be given. Then, there exists a $T_1>0$ such that $x_{1}t\le 1+\epsilon /2$ for all $t\ge T_1$. From () with $\overline{T}=T_1$, we get for all $t\ge T_1\ge 0$.$$\begin{array}{}\text{(16)}& \sigma t=x_{1}t+\frac{a_1}{a_2}{x}_{2}t\le 1+\frac{1}{4d-m_1}-1+\frac{1}{4d-m_1}-x_1{T}_1+\frac{a_1}{a_2}{x}_2{T}_1{e}^{d-m_1{T}_1}{e}^{-d-m_{1}t}.\end{array}$$

Then,$$\begin{array}{}\text{(17)}& x_{1}t+\frac{a_1}{a_2}{x}_{2}t\le 1+\frac{1}{4d-m_1}+\frac{\epsilon }{2}-1+\frac{1}{4d-m_1}+\frac{\epsilon }{2}-x{T}_1+\frac{a_1}{a_2}{x}_2{T}_1{e}^{d-m_1{T}_1}{e}^{-d-m_{1}t}.\end{array}$$

For all, $t\ge T_1$. Let $T_2\ge T_1$ be such that$$\begin{array}{}\text{(18)}& 1+\frac{1}{4d-m_1}+\frac{\epsilon }{2}-x{T}_1+\frac{a_1}{a_2}{x}_2{T}_1{e}^{d-m_1{T}_1}{e}^{-d-m_{1}t}\le \frac{\epsilon }{2}.\end{array}$$

For all, $t\ge T_2$. Then,$$\begin{array}{}\text{(19)}& x_{1}t+\frac{a_1}{a_2}{x}_{2}t\le 1+\frac{1}{4d-m_1}+\epsilon \text{for\hspace{0.17em}all}\hspace{1em}t\ge T_2.\end{array}$$

Hence,$$\begin{array}{}\text{(20)}& {\lim }_{t⟶+\infty }\sup x_{1}t+\frac{a_1}{a_2}{x}_{2}t\le 1+\frac{1}{4d-m_1}.\end{array}$$Model (5) is then obviously dissipative in ${\mathrm{ℛ}}_{+}^2$.

2.2. Existence of Equilibria

We find all the positive equilibrium points of the model (5) by solving the following equations:$$\begin{array}{}\text{(21)}& 1-x_1-\frac{a_1{x}_2}{x_1+b}=0,\text{(22)}& -d+m_1+\frac{a_2{x}_1}{x_1+b}-\frac{\gamma x_2}{\delta +x_2}=0.\end{array}$$

Then, we have the following equilibria for the model (5).

(i) The trivial equilibrium point $E_0\mathrm{0,0}$.

(ii) The first axial equilibrium point $E_1\mathrm{1,0}$.

(iii) The second axial equilibrium point $E_{2}0,\overline{x_2}$, where $\overline{x_2}=--d+m_1\delta /-d+m_1-\gamma$.

(iv) The interior equilibrium point $E^{\ast }{x}_1^{\ast },x_2^{\ast }$, where $x_2^{\ast }=-x_1^{\ast }-1b+x_1^{\ast }/a_1$ and $x_1^{\ast }$ is the positive root of the following equation.$$\begin{array}{}\text{(23)}& {\alpha }_1{x}_1^{\ast 3}+{\alpha }_2{x}_1^{\ast 2}+{\alpha }_3{x}_1^{\ast }+{\alpha }_4=0,\end{array}$$where$$\begin{array}{c}\text{(24)}& {\alpha }_1=--d+m_1-a_2+\gamma ,{\alpha }_2=-d+m_{1}1-b-b-d+m_1+a_{2}1-b-\gamma 1-b+\gamma b,\\ {\alpha }_3=a_1\delta -d+m_1+b-d+m_1+b1-b-d+m_1+a_1{a}_2\delta +b{a}_2-\gamma b-\gamma b1-b,\\ {\alpha }_4=a_1\delta b-d+m_1+b^2-d+m_1-\gamma b^2.\end{array}$$

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^{\ast }{x}_1^{\ast },x_2^{\ast }$ 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), $E_0\mathrm{0,0}$ and $E_1\mathrm{1,0}$ are always exists. $E_{2}0,\overline{x_2}$ is exist, if $-d+m_1>0$ and $-d+m_1<\gamma$. The number of sign changes in ${\alpha }_i,i=\mathrm{1,2,3,4},$ determines the positive roots of (23). Hence, if $x_1^{\ast }<1$, then assume that $E^{\ast }{x}_1^{\ast },x_2^{\ast }$ be the arbitrary interior equilibrium exist for the model (5).

In order to have a quick glance on existence of equilibrium points, let $a_1=0.5,b=0.4,d=0.2,a_2=0.5,\delta =0.2$, $m_1=0.3>d$, and $\gamma =0.2$, we have $E_0\mathrm{0,0}$, $E_1\mathrm{1,0}$, and $E_2\mathrm{0,0.2}$, and from (23), we have$$\begin{array}{}\text{(25)}& -0.4x_1^{\ast 3}+0.28x_1^{\ast 2}+0.196x_1^{\ast }-0.012=0.\end{array}$$

Since two sign change occurs in the above equation, we have two positive roots, that is, $x_1^{\ast }=0.0569659<1$and $x_1^{\ast }=1.11525>1$. Then, from Lemma 1, we have the interior equilibrium $E^{\ast }{x}_1^{\ast },x_2^{\ast }=0.056965,0.861869$for the model (5), which is shown in Figure 1(a). Similarly, if $m_1=0.2=d$and $m_1=0.1<d$, then there is only $E_0$, $E_1$, and $E^{\ast }$, which exist. Moreover, we can say if the death rate of predator $d$ is less than the birth rate due to cannibalism $m_1$, then the predator-only equilibrium $E_2$ can exist. Next, we show the equilibria of the model (5) with the effect of cannibalism graphically. For the different $m_1$ and $\gamma$, we plot the nullcline plots in Figures 1(a)–1(c).

[figure(s) omitted; refer to PDF]

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 equilibrium $E{x}_1,x_2$, which is$$\begin{array}{}\text{(26)}& J=\begin{array}{cc}1-2x_1-\frac{b{a}_1{x}_2}{{x_1+b}^2}& -\frac{a_1{x}_1}{x_1+b}\\ \frac{b{a}_2{x}_2}{{x_1+b}^2}& -d+m_1+\frac{a_2{x}_1}{x_1+b}-\frac{2\gamma \delta x_2+\gamma x_2^2}{{\delta +x_2}^2}\end{array}.\end{array}$$

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) $E_0$ is unstable if $m_1<d$ and saddle if $m_1>d$.

(ii) $E_1$ is stable if $m_1+a_2/1+b<d$ and saddle if $m_1+a_2/1+b>d$.

(iii) $E_2$ is stable only if $1<a_1\overline{x_2}/b$ and $-d+m_1<2\gamma \delta \overline{x_2}+\gamma {\overline{x_2}}^2/{\delta +\overline{x_2}}^2$, otherwise unstable.

Proof.

The Jacobian matrices at $E_0\mathrm{0,0}$ and $E_1\mathrm{1,0}$ are$$\begin{array}{c}\text{(27)}& J_{E_0}=\begin{array}{cc}1\hfill & 0\hfill \\ 0\hfill & -d+m_1\hfill \end{array},J_{E_1}=\begin{array}{cc}-1\hfill & -\frac{a_1}{1+b}\hfill \\ 0\hfill & -d+m_1+\frac{a_2}{1+b}\hfill \end{array}.\end{array}$$

The eigenvalues of $J_{E_0}$ are ${\lambda }_1=1$ and ${\lambda }_2=-d+m_1$. The eigenvalues of $J_{E_1}$ are ${\lambda }_1=-1$ and ${\lambda }_2=-d+m_1+a_2/1+b$.

The Jacobian matrix at $E_{2}0,\overline{x_2}$$$\begin{array}{}\text{(28)}& J_{E_2}=\begin{array}{cc}1-\frac{a_1\overline{x_2}}{b}\hfill & 0\hfill \\ \frac{a_2\overline{x_2}}{b}\hfill & -d+m_1-\frac{2\gamma \delta \overline{x_2}+\gamma {\overline{x_2}}^2}{{\delta +\overline{x_2}}^2}\hfill \end{array}.\end{array}$$

The eigenvalues of $J_{E_2}$ are ${\lambda }_1=1-a_1\overline{x_2}/b$ and ${\lambda }_2=-d+m_1-2\gamma \delta \overline{x_2}+\gamma {\overline{x_2}}^2/{\delta +\overline{x_2}}^2$.

Theorem 4.

For model (5), if$$\begin{array}{}\text{(29)}& {\gamma }^{\ast }<\gamma ,\frac{b{a}_1{a}_2{x}_1^{\ast }{x}_2^{\ast }}{{x_1^{\ast }+b}^3}<-x_1^{\ast }+\frac{a_1{x}_1^{\ast }{x}_2^{\ast }}{{x_1^{\ast }+b}^2}\frac{\delta \gamma x_2^{\ast }}{{\delta +x_2^{\ast }}^2},\end{array}$$then the interior equilibrium $E^{\ast }{x}_1^{\ast },x_2^{\ast }$ is locally asymptotically stable.

Proof.

The Jacobian matrix at $E^{\ast }{x}_1^{\ast },x_2^{\ast }$$$\begin{array}{}\text{(30)}& J_{E^{\ast }}=\begin{array}{cc}-x_1^{\ast }+\frac{a_1{x}_1^{\ast }{x}_2^{\ast }}{{x_1^{\ast }+b}^2}\hfill & -\frac{a_1{x}_1^{\ast }}{x_1^{\ast }+b}\hfill \\ \frac{b{a}_2{x}_2^{\ast }}{{x_1^{\ast }+b}^2}\hfill & -\frac{\delta \gamma x_2^{\ast }}{{\delta +x_2^{\ast }}^2}\hfill \end{array}.\end{array}$$

The characteristic equation of the above matrix is$$\begin{array}{}\text{(31)}& {\lambda }^2-{\beta }_1\lambda +{\beta }_2=0,\end{array}$$where$$\begin{array}{c}\text{(32)}& {\beta }_1=-x_1^{\ast }+\frac{a_1{x}_1^{\ast }{x}_2^{\ast }}{{x_1^{\ast }+b}^2}-\frac{\delta \gamma x_2^{\ast }}{{\delta +x_2^{\ast }}^2}=-x_1^{\ast }+\frac{a_1{x}_1^{\ast }{x}_2^{\ast }}{{x_1^{\ast }+b}^2}\frac{{\delta +x_2^{\ast }}^2}{\delta x_2^{\ast }}-\gamma \\ ={\gamma }^{\ast }-\gamma ,\\ {\beta }_2=-x_1^{\ast }+\frac{a_1{x}_1^{\ast }{x}_2^{\ast }}{{x_1^{\ast }+b}^2}-\frac{\delta \gamma x_2^{\ast }}{{\delta +x_2^{\ast }}^2}+\frac{b{a}_1{a}_2{x}_1^{\ast }{x}_2^{\ast }}{{x_1^{\ast }+b}^3}.\end{array}$$

By Routh–Hurwitz criterion, the roots of (31) has negative real parts if ${\beta }_1<0$ and ${\beta }_2>0$. Hence, if (29) holds, we achieve ${\beta }_1<0$ and ${\beta }_2>0$.

2.4. Hopf Bifurcation

Theorem 5.

Assume that $\gamma ={\gamma }^{\ast }$ and$$\begin{array}{}\text{(33)}& \frac{b{a}_1{a}_2{x}_1^{\ast }{x}_2^{\ast }}{{x_1^{\ast }+b}^3}<-x_1^{\ast }+\frac{a_1{x}_1^{\ast }{x}_2^{\ast }}{{x_1^{\ast }+b}^2}\frac{\delta \gamma x_2^{\ast }}{{\delta +x_2^{\ast }}^2}.\end{array}$$

Then, the model (5) undergoes Hopf bifurcation near $E^{\ast }{x}_1^{\ast },x_2^{\ast }$.

Proof.

It is known that if (i) ${\beta }_1=0$ and (ii) ${\beta }_2>0$, then the roots of (31) have a pair of imaginary roots. Since we aim to study the effect of cannibalism, we take $\gamma$ as a bifurcation parameter. When the parameter $\gamma$ crosses the critical value $\gamma ={\gamma }_c$, (i) and (ii) hold if it satisfies $\gamma ={\gamma }^{\ast }$ and (33), then the roots of (31) has a pair of imaginary roots. Moreover, the transversality condition is given by (iii) $d/d\gamma {\beta }_1=-\delta x_2^{\ast }/{\delta +x_2^{\ast }}^2\ne 0$. These conditions ensure that the model (5) undergoes Hopf bifurcation at $E^{\ast }{x}_1^{\ast },x_2^{\ast }$. Note $E^{\ast }{x}_1^{\ast },x_2^{\ast }$ is in terms of $\gamma$, and it is calculated at the critical value $\gamma ={\gamma }_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 $\tau \ne 0$ with reduced parameters, which is of the form$$\begin{array}{}\text{(34)}& \begin{array}{c}\frac{\text{d}{x}_1}{\text{d}t}=x_{1}1-x_1-\frac{a_1{x}_2}{x_1+b},\hfill \\ \frac{\text{d}{x}_2}{\text{d}t}=-d+m_1{x}_2+\frac{a_2{x}_{1}t-\tau x_{2}t-\tau }{x_{1}t-\tau +b}-\frac{\gamma x_2^2}{\delta +x_2}.\hfill \end{array}\end{array}$$

With $x_1\varphi ={\theta }_1\varphi$, $x_2\varphi ={\theta }_2\varphi$, $\varphi =-\tau ,0$, and ${\theta }_{\mathrm{1,2}}:-\tau ,0⟶{\mathrm{ℛ}}^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^{\ast }{x}_1^{\ast },x_2^{\ast }$.

3.1. Positivity and Boundedness

One can write from the first equation of model (34) as follows:$$\begin{array}{}\text{(35)}& \frac{\text{d}{x}_1}{x_1}=1-x_1-\frac{a_1{x}_2}{x_1+b}.\end{array}$$

Integrating between 0 and $t$, we get$$\begin{array}{}\text{(36)}& x_{1}t=x_{1}0\exp {\displaystyle {\int }_0^{t}1-x_1-\frac{a_1{x}_2}{x_1+b}}ds\hspace{1em}\ge 0.\end{array}$$

Similarly from second equation of model (34), we have$$\begin{array}{}\text{(37)}& x_{2}t=x_{2}0\exp {\displaystyle {\int }_0^t-d+m_1+\frac{a_2{x}_{1}s-\tau }{x_{1}s-\tau +b}-\frac{\gamma x_{2}s}{\delta +x_{2}s}}ds\hspace{1em}\ge 0.\end{array}$$

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^{\ast }=x_1,x_2\in {\mathrm{ℛ}}_{+}^2:\zeta t\le 1+a_1{1-d+m}^2/4\gamma a_2$ as $t⟶\infty$ for all $x_1\varphi ,x_2\varphi \in {\mathrm{ℛ}}_{+}^2$, where $\zeta t=x_{1}t-\tau +a_1/a_2{x}_{2}t$.

Proof.

define a function $\zeta t=x_{1}t-\tau +a_1/a_2{x}_{2}t$. Then, we obtain$$\begin{array}{c}\text{(38)}& \frac{\text{d}\zeta t}{\text{d}t}=\frac{\text{d}{x}_{1}t-\tau }{\text{d}t}+\frac{a_1}{a_2}\frac{\text{d}{x}_{2}t}{\text{d}t}=x_{1}t-\tau 1-x_{1}t-\tau -\frac{a_1{x}_{1}t-\tau x_{2}t-\tau }{x_{1}t-tau+b}\frac{a_1}{a_2}-d+m_1{x}_2+\frac{a_1{x}_{1}t-\tau x_{2}t-\tau }{x_{1}t-\tau +b}-\frac{a_1}{a_2}\frac{\gamma x_2^{2}t}{\delta +x_{2}t}\\ =-x_{1}t-\tau +\frac{a_1}{a_2}{x}_2+2x_{1}t-\tau -x_1^{2}t-\tau +\frac{a_1\gamma }{a_2}{x}_{2}t\frac{-d+m+1}{\gamma }-\frac{x_{2}t}{\delta +x_{2}t},\\ \frac{\text{d}\zeta t}{\text{d}t}+\zeta t=2x_{1}t-\tau -x_1^{2}t-\tau +\frac{a_1\gamma }{a_2}{x}_{2}t\frac{-d+m+1}{\gamma }-\frac{x_{2}t}{\delta +x_{2}t}\\ \le 1+\frac{a_1{1-d+m}^2}{4\gamma a_2}.\end{array}$$

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 $x_1=x_1-x_1^{\ast }$ and $x_2=x_2-x_2^{\ast }$ is given as follows:$$\begin{array}{c}\text{(39)}& \begin{array}{c}{\dot{x}}_1\hfill \\ {\dot{x}}_2\hfill \end{array}=\begin{array}{cc}-x_1^{\ast }+\frac{a_1{x}_1^{\ast }{x}_2^{\ast }}{{x_1^{\ast }+b}^2}\hfill & -\frac{a_1{x}_1^{\ast }}{x_1^{\ast }+b}\hfill \\ 0\hfill & -d+m_1-\frac{2\gamma \delta x_2^{\ast }-\gamma x_2^{\ast 2}}{{\delta +x_2^{\ast }}^2}\hfill \end{array}\begin{array}{c}{x}_1\hfill \\ x_2\hfill \end{array}+\begin{array}{cc}0\hfill & 0\hfill \\ \frac{b{a}_2{x}_2^{\ast }}{{x_1^{\ast }+b}^2}\hfill & \frac{a_2{x}_1^{\ast }}{x_1^{\ast }+b}\hfill \end{array}\begin{array}{c}{x}_1-\tau \hfill \\ x_2-\tau \hfill \end{array},\\ =\begin{array}{cc}{\vartheta }_1\hfill & {\vartheta }_2\hfill \\ 0\hfill & {\vartheta }_3\hfill \end{array}\begin{array}{c}{x}_1\hfill \\ x_2\hfill \end{array}+\begin{array}{cc}0\hfill & 0\hfill \\ {\vartheta }_4\hfill & {\vartheta }_5\hfill \end{array}\begin{array}{c}{x}_1-\tau \hfill \\ x_2-\tau \hfill \end{array}.\end{array}$$

We removed the bars for our convenience, and then, the characteristic equation for (39) is given by$$\begin{array}{}\text{(40)}& \begin{array}{cc}{\vartheta }_1-\lambda \hfill & {\vartheta }_2\hfill \\ {\vartheta }_4{e}^{-\lambda \tau }\hfill & {\vartheta }_3+{\vartheta }_5{e}^{-\lambda \tau }-\lambda \hfill \end{array}=0.\end{array}$$

That is,$$\begin{array}{}\text{(41)}& {\lambda }^2-{\vartheta }_1+{\vartheta }_3\lambda +{\vartheta }_1{\vartheta }_3+-{\vartheta }_5\lambda +{\vartheta }_1{\vartheta }_5-{\vartheta }_2{\vartheta }_4{e}^{-\lambda \tau }=0.\end{array}$$

Substitute $\lambda =i\varphi$ in (41) and separate the real and imaginary parts:$$\begin{array}{c}\text{(42)}& {\vartheta }_2{\vartheta }_4-{\vartheta }_1{\vartheta }_5\sin \varphi \tau -\varphi {\vartheta }_5\cos \varphi \tau =\varphi {\vartheta }_1+{\vartheta }_3,{\vartheta }_5\varphi \sin \varphi \tau +{\vartheta }_2{\vartheta }_4-{\vartheta }_1{\vartheta }_5\cos \varphi \tau =-{\varphi }^2+{\vartheta }_1{\vartheta }_3.\end{array}$$

After simplification, we get$$\begin{array}{c}\text{(43)}& \sin \varphi \tau =\frac{\varphi {\vartheta }_1+{\vartheta }_3{\vartheta }_2{\vartheta }_4-{\vartheta }_1{\vartheta }_5+{\vartheta }_5\varphi -{\varphi }^2+{\vartheta }_1{\vartheta }_3}{{\varphi }^2{\vartheta }_5^2+{{\vartheta }_2{\vartheta }_4-{\vartheta }_1{\vartheta }_5}^2},\cos \varphi \tau =\frac{-{\varphi }^2+{\vartheta }_1{\vartheta }_3{\vartheta }_2{\vartheta }_4-{\vartheta }_1{\vartheta }_5-{\vartheta }_5{\varphi }^2{\vartheta }_1+{\vartheta }_3}{{\varphi }^2{\vartheta }_5^2+{{\vartheta }_2{\vartheta }_4-{\vartheta }_1{\vartheta }_5}^2}.\end{array}$$

Since ${\sin }^2\varphi \tau +{\cos }^2\varphi \tau =1$, we have$$\begin{array}{}\text{(44)}& {\phi }_1{\varphi }^6+{\phi }_2{\varphi }^4+{\phi }_3{\varphi }^2+{\phi }_4=0,\end{array}$$where$$\begin{array}{c}\text{(45)}& {\phi }_1={\vartheta }_5^2,{\phi }_2=-2{\vartheta }_1{\vartheta }_3{\vartheta }_5^2+{{\vartheta }_2{\vartheta }_4-{\vartheta }_1{\vartheta }_5}^2+{\vartheta }_5^2{{\vartheta }_1+{\vartheta }_3}^2-{\vartheta }_5^4,\\ {\phi }_3={{\vartheta }_2{\vartheta }_4-{\vartheta }_1{\vartheta }_5}^2{{\vartheta }_1+{\vartheta }_3}^2+{\vartheta }_1^2{\vartheta }_3^2{\vartheta }_5^2-2{{\vartheta }_2{\vartheta }_4-{\vartheta }_1{\vartheta }_5}^2{\vartheta }_1{\vartheta }_3-2{{\vartheta }_2{\vartheta }_4-{\vartheta }_1{\vartheta }_5}^2{\vartheta }_5^2,\\ {\phi }_4={\vartheta }_1^2{\vartheta }_3^2{{\vartheta }_2{\vartheta }_4-{\vartheta }_1{\vartheta }_5}^2-{{\vartheta }_2{\vartheta }_4-{\vartheta }_1{\vartheta }_5}^4.\end{array}$$

Let us assume that (44) has at least one positive root, which is ${\varphi }^{\ast }$; then, we get the critical value ${\tau }_n$ as follows:$$\begin{array}{}\text{(46)}& {\tau }_n=\frac{1}{{\varphi }^{\ast }}arc\sin \frac{{\varphi }^{\ast }{\vartheta }_1+{\vartheta }_3{\vartheta }_2{\vartheta }_4-{\vartheta }_1{\vartheta }_5+{\vartheta }_5{\varphi }^{\ast }-{\varphi }^{\ast 2}+{\vartheta }_1{\vartheta }_3}{{\varphi }^{\ast 2}{\vartheta }_5^2+{{\vartheta }_2{\vartheta }_4-{\vartheta }_1{\vartheta }_5}^2}+\frac{2n\pi }{{\varphi }^{\ast }},\hspace{1em}n=\mathrm{0,1,2},\dots .\end{array}$$

Theorem 7.

The following transversality condition is holds:$$\begin{array}{}\text{(47)}& \frac{dR\lambda \varphi }{\text{d}\tau }{|}_{\tau ={\tau }^{\ast }}\ne 0.\end{array}$$

Proof.

Let us substitute $\lambda \tau =\theta \tau +i\varphi \tau$ in () and differentiate $w$$r$$t$$\tau$, we get$$\begin{array}{}\text{(48)}& {\frac{\text{d}\lambda }{\text{d}\tau }}^{-1}=\frac{-2\lambda +{\vartheta }_1+{\vartheta }_3{e}^{\lambda \tau }+{\vartheta }_5}{\lambda {\vartheta }_5\lambda -{\vartheta }_1{\vartheta }_5-{\vartheta }_2{\vartheta }_4}-\frac{\tau }{\lambda }.\end{array}$$

Then,$$\begin{array}{c}\text{(49)}& R{\frac{\text{d}\lambda }{\text{d}\tau }}^{-1}{|}_{\lambda =i{\varphi }^{\ast }}=\frac{-2\lambda +{\vartheta }_1+{\vartheta }_3{e}^{\lambda \tau }+{\vartheta }_5}{\lambda {\vartheta }_5\lambda -{\vartheta }_1{\vartheta }_5-{\vartheta }_2{\vartheta }_4}{|}_{\lambda =i\varphi }=\frac{1}{{\vartheta }_5^2{\varphi }^{\ast 3}+{\varphi }^{\ast }{{\vartheta }_2{\vartheta }_4-{\vartheta }_1{\vartheta }_5}^2}\\ -{\vartheta }_5^2+{\vartheta }_1{\vartheta }_5-{\vartheta }_2{\vartheta }_4{\vartheta }_1+{\vartheta }_3-2{\vartheta }_5{\varphi }^{\ast }{}^2\sin {\varphi }^{\ast }\tau -{\vartheta }_5{\varphi }^{\ast }{\vartheta }_1+{\vartheta }_3-2{\varphi }^{\ast }{\vartheta }_1{\vartheta }_5-{\vartheta }_2{\vartheta }_4\cos {\varphi }^{\ast }\tau .\end{array}$$

Theorem 8.

If (29) holds then for the model (34), we have the following: (i) the interior equilibrium point $E^{\ast }{x}_1^{\ast },x_2^{\ast }$ is locally asymptotically stable for $\tau \in 0,{\tau }^{\ast }$ and (ii) undergoes Hopf bifurcation near $E^{\ast }{x}_1^{\ast },x_2^{\ast }$ at $\tau ={\tau }_n$, $n=\mathrm{0,1,2},\dots$. Also, each ${\tau }_n$ satisfies (47), which is the critical points that were the stability switches of $E^{\ast }$ 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$$\begin{array}{c}\text{(50)}& a_1=0.5,b=0.4,\\ d=0.1,\\ a_2=\mathrm{0.5.}\end{array}$$

Case 1.

The nondelayed model is as follows:

With the parameter values in (50), the model (5) becomes$$\begin{array}{}\text{(51)}& \begin{array}{c}\frac{\text{d}{x}_1}{\text{d}t}=x_{1}1-x_1-\frac{0.5x_2}{x_1+0.4},\hfill \\ \frac{\text{d}{x}_2}{\text{d}t}=-0.1+m_1{x}_2+\frac{0.5x_1{x}_2}{x_1+0.4}-\frac{\gamma x_2^2}{\delta +x_2}.\hfill \end{array}\end{array}$$

First, we showed that for $m_1=0.3>d$, $\gamma =0.2$, and $\delta =0.2$, then the model (51) has the equilibrium points $E_0\mathrm{0,0}$, $E_1\mathrm{1,0}$, $E_2\mathrm{0,0.2}$, and $E^{\ast }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 $E_2$ is exist only if $m_1>d$, which is shown in Figure 1(a). If $E_2$ is stable, then the predator may alone survive without prey by getting energy from its own species $m_1>d$; that is, a death rate of the predator is greater than birth due to cannibalism. If $E_2$ is unstable, then the trajectories near $E_2$ either approach $E_0$ or $E^{\ast }$. From Theorem 3, we have $E_0$ is a saddle point, $E_1$ is also a saddle point, since $m_1+a_2/1+b=0.657143>d$. And, $a_1\overline{x_2}/b<1$, then $E_2$ is unstable for the model (51). Also, the predator-only equilibrium $E_2$ does not exist for $m_1\le d$ (see Figures 1(b) and 1(c)). The local stability of $E^{\ast }$ is achieved in Theorem 4, which ensures the long-time survival of both species. The model (51) at $\gamma =0.36$, $d=0.1$, and $m_1=0.2$ with other parameters in (50) is locally asymptotically stable near the interior equilibrium $E^{\ast }$, see Figure 2. Furthermore, if $\gamma ={\gamma }^{\ast }=0.34324$, also, then the model (51) has interior equilibrium $E^{\ast }{x}_1^{\ast },x_2^{\ast }=0.234183,0.971336$. Then, the following conditions hold: (i) ${\beta }_1=0$; (ii) ${\beta }_2=0.0868205$; (iii) $d/d\gamma {\beta }_1=-0.141591\ne 0$. Then, from Theorem 5, the model (51) undergoes Hopf bifurcation at $\gamma ={\gamma }^{\ast }=0.34324$ at critical value ${\gamma }_c=0.34324$ and is shown in Figure 3 for $\gamma =0.34$. Further, on increasing $\delta =0.8$, then the model (51) is locally asymptotically stable (see Figure 4). For clear representation, one-parameter bifurcation diagram with respect to $\gamma \in 0.2,0.5$ and $\delta \in \mathrm{0,0.8}$ is shown in Figures 5 and 6, correspondingly. In Figure 5, it is shown that for $\gamma \in \mathrm{0,0.34324}$, the periodic solution arises near the $E^{\ast }$, and for $\gamma \in 0.34324,0.5$, the $E^{\ast }$ is locally asymptotically stable. Also in Figure 5, we noticed on increasing the $\gamma$ value, the size of prey population increases, and the predator population size declines. Figure 6 shows that $E^{\ast }$ is locally asymptotically stable for $\delta \in \mathrm{0,0.1905}$, the periodic solution arises near the $E^{\ast }$ for $\delta \in 0.1905,0.593$, and the $E^{\ast }$ is locally asymptotically stable $\gamma \in 0.593,0.8$. Hence, we notice from Figure 6, that on increasing the values of $\delta$, both population sizes reduce. Moreover, the existence of exchange of stability (unstable to stable) on varying the parameters $\gamma$ and $\delta$ occurs; hence, we can conclude that the parameters $\gamma$ and $\delta$ has stabilizing effect in the nondelayed model (5).

Next, to demonstrate the impact of cannibalism parameters $\gamma$, $\delta$, and $m_1$ on the model (51), we plotted the two-parameter bifurcation diagram with respect to the parameters $\delta$ and $\gamma$ in Figure 7(a), and $m_1$ and $\delta$ in Figure 7(b). It is shown that less $\gamma$ value causes the extinction of species. Also, the stable, unstable, and extinction regions are clearly depicted. Further, we showed that the $\gamma$ value should be above some threshold value in order to ensure the survival of the species. Similarly, $m_1$ 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(s) omitted; refer to PDF]

Case 2.

The delayed model.

Let $d=0.1$, $m_1=0.2$, $\gamma =0.4$, $\delta =0.4$, and other parameter values in (50); then, the model (34) becomes$$\begin{array}{}\text{(52)}& \begin{array}{c}\frac{\text{d}{x}_1}{\text{d}t}=x_{1}1-x_1-\frac{0.5x_2}{x_1+0.4},\hfill \\ \frac{\text{d}{x}_2}{\text{d}t}=-0.1+m_1{x}_2+\frac{0.5x_{1}t-\tau x_{2}t-\tau }{x_{1}t-\tau +0.4}-\frac{\gamma x_2^2}{\text{delta}+x_2}.\hfill \end{array}\end{array}$$

Then, the model (52) has interior equilibrium point $E^{\ast }{x}_1^{\ast },x_2^{\ast }=0.346151,0.97574$. Moreover, from (44), we have$$\begin{array}{}\text{(53)}& 0.0538044{\varphi }^6+0.00677379{\varphi }^4-0.0001068487{\varphi }^2-0.0000251728=0.\end{array}$$

If it has a positive root ${\varphi }^{\ast }=0.2374849$, then we get the critical value $\tau =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 $\tau \in \mathrm{0,1.69709}$, model (52) is locally asymptotically stable (see Figure 8) and undergoes Hopf bifurcation at ${\tau }^{\ast }=1.69709$ (see Figure 9). For clear presentation, the bifurcation diagram with respect to $\tau$ 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 $\tau$ with the cannibalism parameters $\gamma$ and $\delta$, we plotted the two-parameter bifurcation diagram with respect to $\delta$ and $\tau$ in Figure 11(a), and $\gamma$ and $\tau$ 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 $\delta$ value, the interval of stability with respect to $\tau$ decreases. For instance, if $\delta =0.2$, the interior equilibrium $E^{\ast }$ for the model (52) is locally asymptotically stable for $\tau \in \mathrm{0,1.69704}$ and Hopf bifurcation occurs when $\tau$ crosses the critical value $\tau =1.69704$, and if $\delta =0.5$, the interior equilibrium $E^{\ast }$ for the model (52) is locally asymptotically stable for $\tau \in \mathrm{0,0.20826}$ and Hopf bifurcation occurs when $\tau$ crosses the critical value $\tau =0.20826$. Similarly, for larger values of $\gamma$, the interval for local stability increases for $\tau$.

[figure(s) omitted; refer to PDF]

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 $E_2$ 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.

References

[1] H. Zhang, Y. Cai, S. Fu, W. Wang, "Impact of the fear effect in a prey-predator model incorporating a prey refuge," Applied Mathematics and Computation, vol. 356, pp. 328-337, DOI: 10.1016/j.amc.2019.03.034, 2019.

[2] S. Vinoth, R. Sivasamy, K. Sathiyanathan, B. Unyong, G. Rajchakit, R. Vadivel, N. Gunasekaran, "The dynamics of a Leslie type predator-prey model with fear and allee effect," Advances in Difference Equations, vol. 2021 no. 1,DOI: 10.1186/s13662-021-03490-x, 2021.

[3] S. Vinoth, R. Sivasamy, K. Sathiyanathan, G. Rajchakit, P. Hammachukiattikul, R. Vadivel, N. Gunasekaran, "Dynamical analysis of a delayed food chain model with additive allee effect," Advances in Difference Equations, vol. 2021 no. 1,DOI: 10.1186/s13662-021-03216-z, 2021.

[4] S. Vinoth, R. Sivasamy, K. Sathiyanathan, B. Unyong, R. Vadivel, N. Gunasekaran, "A novel discrete-time Leslie-Gower model with the impact of allee effect in predator population," Complexity, vol. 2022,DOI: 10.1155/2022/6931354, 2022.

[5] H. Deng, F. Chen, Z. Zhu, Z. Li, "Dynamic behaviors of Lotka-Volterra predator-prey model incorporating predator cannibalism," Advances in Difference Equations, vol. 2019 no. 1,DOI: 10.1186/s13662-019-2289-8, 2019.

[6] Q. Lin, C. Liu, X. Xie, Y. Xue, "Global attractivity of Leslie-Gower predator-prey model incorporating prey cannibalism," Advances in Difference Equations, vol. 2020 no. 1,DOI: 10.1186/s13662-020-02609-w, 2020.

[7] H. J. Alsakaji, S. Kundu, F. A. Rihan, "Delay differential model of one-predator two-prey system with Monod-Haldane and holling type II functional responses," Applied Mathematics and Computation, vol. 397,DOI: 10.1016/j.amc.2020.125919, 2021.

[8] J. D. Murray, Mathematical Biology I. An Introduction, 2002.

[9] T. Peter, Complex Population Dynamics, 2013.

[10] M. L. Rosenzweig, R. H. MacArthur, "Graphical representation and stability conditions of predator-prey interactions," The American Naturalist, vol. 97 no. 895, pp. 209-223, DOI: 10.1086/282272, 1963.

[11] D. Claessen, A. M. De Roos, L. Persson, "Population dynamic theory of size-dependent cannibalism," Proceedings of the Royal Society of London, Series B: Biological Sciences, vol. 271 no. 1537, pp. 333-340, DOI: 10.1098/rspb.2003.2555, 2004.

[12] Ph Getto, O. Diekmann, A. de Roos, "On the (dis) advantages of cannibalism," Journal of Mathematical Biology, vol. 51 no. 6, pp. 695-712, DOI: 10.1007/s00285-005-0342-6, 2005.

[13] A. Al Basheer, R. D. Parshad, E. Quansah, S. Yu, R. K. Upadhyay, "Exploring the dynamics of a Holling-Tanner model with cannibalism in both predator and prey population," International Journal of Biomathematics, vol. 11 no. 1,DOI: 10.1142/s1793524518500109, 2018.

[14] F. A. Rihan, H. J. Alsakaji, "Persistence and extinction for stochastic delay differential model of prey predator system with hunting cooperation in predators," Advances in Difference Equations, vol. 2020 no. 1,DOI: 10.1186/s13662-020-02579-z, 2020.

[15] F. A. Rihan, H. J. Alsakaji, C. Rajivganthi, "Stability and hopf bifurcation of three-species prey-predator system with time delays and allee effect," Complexity, vol. 2020,DOI: 10.1155/2020/7306412, 2020.

[16] K. Yang, Delay Differential Equations, 2012.

[17] F. A. Rihan, Delay Differential Equations and Applications to Biology, 2021.

[18] F. A. Rihan, H. J. Alsakaji, "Stochastic delay differential equations of three-species prey-predator system with cooperation among prey species," Discrete & Continuous Dynamical Systems-S, vol. 15 no. 2,DOI: 10.3934/dcdss.2020468, 2022.

[19] Z. G. Song, B. Zhen, J. Xu, "Species coexistence and chaotic behavior induced by multiple delays in a food chain system," Ecological Complexity, vol. 19,DOI: 10.1016/j.ecocom.2014.01.004, 2014.

[20] R. Vadivel, P. Hammachukiattikul, N. Gunasekaran, R. Saravanakumar, H. Dutta, "Strict dissipativity synchronization for delayed static neural networks: an event-triggered scheme," Chaos, Solitons & Fractals, vol. 150,DOI: 10.1016/j.chaos.2021.111212, 2021.

[21] R. Vadivel, S. Saravanan, B. Unyong, P. Hammachukiattikul, K.-S. Hong, G. M. Lee, "Stabilization of delayed fuzzy neutral-type systems under intermittent control," International Journal of Control, Automation and Systems, vol. 19 no. 3, pp. 1408-1425, DOI: 10.1007/s12555-020-0526-2, 2021.

[22] J. K. Hale, S. M. V. Lunel, "Introduction to functional differential equations," Springer Science & Business Media, vol. 99, 2013.

[23] G. Birkhoff, G. C Rota, "Ordinray Differential Equations," Ginn Boston, 1982.

[24] M. A. Aziz-Alaoui, "Study of a Leslie-Gower-type tritrophic population model," Chaos, Solitons & Fractals, vol. 14 no. 8, pp. 1275-1293, DOI: 10.1016/s0960-0779(02)00079-6, 2002.

[25] N. Pal, S. Samanta, S. Biswas, M. Alquran, K. Al-Khaled, J. Chattopadhyay, "Stability and bifurcation analysis of a three-species food chain model with delay," International Journal of Bifurcation and Chaos, vol. 25 no. 9,DOI: 10.1142/s0218127415501230, 2015.

[26] K. Manna, M. Banerjee, "Stability of hopf-bifurcating limit cycles in a diffusion-driven prey-predator system with allee effect and time delay," Mathematical Biosciences and Engineering, vol. 16 no. 4, pp. 2411-2446, DOI: 10.3934/mbe.2019121, 2019.

[27] M. Banerjee, L. Zhang, "Time delay can enhance spatio-temporal chaos in a prey-predator model," Ecological Complexity, vol. 27, pp. 17-28, DOI: 10.1016/j.ecocom.2015.12.001, 2016.

[28] S. Magudeeswaran, S. Vinoth, K. Sathiyanathan, M. Sivabalan, "Impact of fear on delayed three species food-web model with holling type-II functional response," International Journal of Biomathematics, vol. 15 no. 4,DOI: 10.1142/s1793524522500140, 2022.

[29] A. Basheer, E. Quansah, S. Bhowmick, R. D. Parshad, "Prey cannibalism alters the dynamics of Holling-Tanner-type predator-prey models," Nonlinear Dynamics, vol. 85 no. 4, pp. 2549-2567, DOI: 10.1007/s11071-016-2844-8, 2016.