Dynamical Analysis of a Predator-Prey Model Incorporating Predator Cannibalism and Refuge

Abstract

We consider a mathematical model to describe the interaction between predator and prey that includes predator cannibalism and refuge. We aim to study the dynamics and its long-term behavior of the proposed model, as well as to discuss the effects of crucial parameters associated with the model. We first show the boundedness and positivity of the solution of the model. Then, we study the existence and stability of all possible equilibrium points. The local stability of the model around each equilibrium point is studied via the linearized system, while the global stability is performed by defining a Lyapunov function. The model has four equilibrium points. It is found that the equilibrium point representing the extinction of both prey and predator populations is always unstable, while the other equilibrium points are conditionally stable. In addition, there is forward bifurcation phenomena that occur under certain condition. To support our analytical findings, we perform some numerical simulations.

1. Introduction

Predator–prey interaction is one of the most important issues in ecology, as it is the basis of the food chain. Many mathematical models have been proposed in the literature to understand the dynamics of predator–prey interaction. So far, the development of the predator–prey model is still continuing to accommodate more realistic phenomena. One of the interesting things is cannibalism in the predator–prey system. Cannibalism, or intraspecific predation, is the consumption of the whole or part of another individual of the same species [1,2]. The occurrence of cannibalism is influenced by several important factors, which are temperature, population density, size, stage of development, and so on [3]. The mathematical model of cannibalism has been studied by some researchers [4,5,6]. A predator–prey model with cannibalism is enticing to study, because, in fact, many species in nature have cannibalistic traits. Numerous groups of animals have been found to have cannibalism, i.e., insects [7], primates [8], frogs [9], fish [10], carnivore mammals [11], spiders [12], etc.

Kang et al. [4] studied a single-species cannibalism model with stage structure. The model studied is a dynamical system of one population with an age structure that divides the population into two classes, namely eggs and adult class consisting of larvae, pupae, queen insects, worker insects, and other types. Zhang et al. [5] analyzed predator–prey models with cannibalism and stage structure in predators so that the model studied was a three-dimensional dynamical model. In Zhang’s model, the predator population is divided into two subpopulations, i.e., juvenile predators and adult predators. The birth rate of juvenile predators is assumed to be proportional to the number of adult predators, and follows Malthus’ growth model. Predation of prey and juvenile predators by adult predators follows the type I Holling functional response. This research resulted in local and global stability analysis of the equilibrium point, as well as the forward bifurcation. Meanwhile, Deng et al. [6] studied a two-dimensional predator–prey model with predator cannibalism:

$$\begin{array}{ccc}{\displaystyle \frac{dN}{dt}}& =& {\displaystyle rN\left(1-\frac{N}{K}\right)-b_{1}NP,}\hfill \\ {\displaystyle \frac{dP}{dt}}& =& {\displaystyle c_{1}NP+c_{2}P-eP-\frac{b_2{P}^2}{k_2+P},}\hfill \end{array}$$

(1)

where $N\ge 0$ and $P\ge 0$ represent prey density and predator density. The parameters r, K, $b_1$, $c_1$, $c_2$, e, $b_2$, and $k_2$ are positive constants that, respectively, represent the intrinsic growth rate prey, environmental carrying capacity for prey, the rate of prey predation, conversion rate of prey biomass into predator birth, conversion rate of cannibalism into predator birth, predator death rate, the rate of cannibalism in predators, and the half-saturation constant of cannibalism. The phenomenon of cannibalism is represented by the last term and the second term in the second equation of the system (1).

Besides cannibalism, the phenomenon of predator–prey that is also appealing to study is the prey hiding behavior from predator’s captures and attacks. In the concept of ecology, this behavior is called refuge. Many prey species adopt the technique of refuge to avoid the predation. For example, Daphnia against fish predation in Mediterranean shallow lakes by hiding in the sediments [13], and sea-urchins hide the juvenile ones in articulated coralline algae from predatory crabs [14]. Apart from the natural behavior of prey, the prey refuge can be done with humans’ help, such as creating conservation forests [15], nature reserves, wildlife reserves, or even a simple protection. The mathematical model of predator–prey with prey refuge has also been widely studied [16,17,18]. The technique of refuge is adopted by various species to avoid the cannibalism, such as the komodo dragon and wolf spider. The females of komodo dragons sometimes build decoy chambers in their nests to protect their eggs from other dragons [19]. Then, the young dragons climb trees to stay safe from adult dragons and stay in the trees for two to four years [20]. The wolf spiders hide from stronger ones in leaf litter [21], especially the newly hatched wolf spiders, which take refuge from cannibalism by riding on their mother’s back for several days [22]. Therefore, in the predators–prey model with cannibalistic predators, it can be assumed that there is refuge for cannibalized predators. The predator–prey model with predator refuge is still very rarely studied, so it becomes a very enticing thing to study.

As far as we know, the interaction of predator and prey with predator cannibalism and predator refugee has never been studied mathematically. Thus, we propose a model describing the predator–prey interaction incorporating predator cannibalism and refuge and then perform a dynamical analysis for the proposed model. The proposed model is a development of Deng’s model [6], namely by implementing the Holling type II functional response instead of the Holling type I and assuming that there is predator refuge from cannibalization. The Holling type II functional response is basically implemented to describe the saturated predation mechanism.

The sections of this paper are organized as follows. The development of the model is described in Section 2, followed by the verification on the existence, uniqueness, nonnegativity, and boundedness of solutions of the developed model in Section 3. The existence and analysis of the local stability of the equilibrium points of the model are discussed in Section 4, while the analysis of the global stability of the equilibrium points and the analysis of the forward bifurcation are described, respectively, in Section 5 and Section 6. The numerical simulations and the intrepretations are carried out in Section 7. Finally, we draw some concluding remarks in Section 8.

2. Model Development

If in model (1) it is assumed that there is protection in the cannibalized predator population, called predator refuge, as much as $mP$, then the number of predators available for cannibalization is $(1-m)P$. Hence, the model (1) is modified into

$$\begin{array}{ccc}{\displaystyle \frac{dN}{dt}}& =& {\displaystyle rN\left(1-\frac{N}{K}\right)-b_{1}NP,}\hfill \\ {\displaystyle \frac{dP}{dt}}& =& {\displaystyle c_{1}NP+c_{2}P-eP-\frac{b_2(1-m)P^2}{k_2+(1-m)P}.}\hfill \end{array}$$

(2)

By adding the assumption that predators need time to catch and handle the prey, then the rate of prey predation follows the Holling type II functional response, so that the model (2) becomes

$$\begin{array}{ccc}{\displaystyle \frac{dN}{dt}}& =& {\displaystyle rN\left(1-\frac{N}{K}\right)-\frac{b_{1}NP}{k_1+N},}\hfill \\ {\displaystyle \frac{dP}{dt}}& =& {\displaystyle \frac{c_{1}NP}{k_1+N}+c_{2}P-eP-\frac{b_2(1-m)P^2}{k_2+(1-m)P}.}\hfill \end{array}$$

(3)

The parameter $k_1$ represents half-saturation constant of predation. Holling type II functional response is more realistic than Holling type I because the rate of the prey predation is saturated [23]. This is in accordance with the real conditions, that it is impossible to predators to eat the prey continuously. Even less when the number of prey is abundant, the predators will experience satiety or reach the saturation point.

3. Preliminaries Results

3.1. Existence and Uniqueness

In this section, we show the existence and uniqueness of solution of the system (3) in $\mathsf{\Omega }\times [0,T]$ where $\mathsf{\Omega }=\{X=(N,P)\in {\mathbb{R}}_{+}^2\}$ and $T<\infty$. For this aim, we denote $F\left(X\right)=(F_1\left(X\right),F_2\left(X\right))$, where

$$\begin{array}{ccc}{F}_1\left(X\right)& =& {\displaystyle rN\left(1-\frac{N}{K}\right)-\frac{b_{1}NP}{k_1+N},}\hfill \\ F_2\left(X\right)& =& {\displaystyle \frac{c_{1}NP}{k_1+N}+c_{2}P-eP-\frac{b_2(1-m)P^2}{k_2+(1-m)P},}\hfill \end{array}$$

(4)

such that the system (3) can be written as

$$\frac{dX}{dt}=F\left(X\right).$$

It can be verified that $F_i,\frac{\partial F_i}{\partial N}$, and $\frac{\partial F_i}{\partial P},i=1,2$ are continuous in $\mathsf{\Omega }$. Based on a lemma in [24] (p. 71), $F\left(X\right)$ is locally Lipschitz on $\mathsf{\Omega }$. Consequently, using the fundamental existence-uniqueness theorem (see [24] (p. 74)), we obtain the following theorem.

Theorem 1.

For the system (3) with any non-negative initial condition $N\left(0\right)\ge 0$ and $P\left(0\right)\ge 0$, there exists $T>0$ so that the system (3) has a unique solution defined in Ω.

3.2. Nonnegativity

Since the variables in the system (3) represent the population densities, the solution of the system must be non-negative. The solution of system (3) is guaranteed to be non-negative by the following theorem.

Theorem 2.

All solutions of (3) with initial values $(N\left(0\right),P\left(0\right))\in {\mathbb{R}}_{+}^2$ are non-negative.

Proof. 

We will first prove that if $N\left(0\right)\ge 0$ and $P\left(0\right)\ge 0$ then $N\left(t\right)\ge 0$ and $P\left(t\right)\ge 0$ for every $t>0$. If $N\left(0\right)=0$ then

$$\frac{dN}{dt}=0,t=0.$$

It means that the prey population density N doesn’t change from the beginning to the next. Subsequently, it is assumed that $N\left(0\right)>0$. If $N\left(t\right)\ge 0$ for every $t\ge 0$ is not true, then there is $t_1>0$ such that $N\left(t\right)>0$ for $0<t<t_1$, $N\left(t\right)=0$ for $t=t_1$ and $N\left(t\right)<0$ for $t>t_1$. From (3), we obtain:

$$\frac{dN}{dt}=0,t=t_1.$$

Thus, there is no change in the population density of N when $t=t_1$. This contradicts the statement that $N\left(t\right)<0$ for $t>t_1$. Therefore, the previous assumption is false, which means $N\left(t\right)\ge 0$ for every $t>0$. In the same way, it can be proved that $P\left(t\right)\ge 0$ for every $t>0$. □

3.3. Boundedness

Predator and prey populations in the system (3) must be limited due to the limited carrying capacity of the prey and predator resources.

Theorem 3.

All solutions of (3) with initial values $(N\left(0\right),P\left(0\right))\in {\mathbb{R}}_{+}^2$ are uniformly bounded.

Proof. 

Consider a function defined by

$$V\left(t\right)=p_{1}N\left(t\right)+p_{2}P\left(t\right),$$

where $p_1,p_2>0$. $V\left(t\right)$ has the first derivative:

$$\begin{array}{ccc}{\displaystyle \frac{dV}{dt}}& =& {\displaystyle p_1\frac{dN}{dt}+p_2\frac{dP}{dt}}\hfill \\ & =& {\displaystyle p_1\left[rN\left(1-\frac{N}{K}\right)-\frac{b_{1}NP}{k_1+N}\right]+p_2\left[\frac{c_{1}NP}{k+N}+c_{2}P-eP-\frac{b_2(1-m)P^2}{k_2+(1-m)P}\right]}\hfill \\ & <& {\displaystyle p_{1}rN-\frac{p_{1}r}{K}{N}^2+\frac{p_2{c}_1-p_1{b}_1}{k_1+N}NP+p_2(c_2-e)P.}\hfill \end{array}$$

If $p_1=c_1$ and $p_2=b_1$, then

$${\displaystyle \frac{dV}{dt}<c_{1}rN-\frac{c_{1}r}{K}{N}^2+b_1(c_2-e)P.}$$

For any positive constant $\phi$, it holds that

$$\begin{array}{ccc}{\displaystyle \frac{dV}{dt}+\phi V}& <& {\displaystyle c_{1}rN-\frac{c_{1}r}{K}{N}^2+b_1(c_2-e)P+\phi (c_{1}N+b_{1}P)}\hfill \\ & =& {\displaystyle c_1(r+\phi )N-\frac{c_{1}r}{K}{N}^2+b_1(c_2-e+\phi )P.}\hfill \end{array}$$

By choosing $\phi <e-c_2$, we have

$$\begin{array}{ccc}{\displaystyle \frac{dV}{dt}+\phi V}& <& {\displaystyle c_1(r+\phi )N-\frac{c_{1}r}{K}{N}^2}\hfill \\ & =& {\displaystyle -\frac{c_{1}r}{K}\left[{\left(N-\frac{(r+\phi )K}{2r}\right)}^2-{\left(\frac{(r+\phi )K}{2r}\right)}^2\right]}\hfill \\ & \le & {\displaystyle \frac{c_{1}r}{K}{\left(\frac{(r+\phi )K}{2r}\right)}^2.}\hfill \end{array}$$

Thus, ${\displaystyle \frac{dV}{dt}+\phi V\left(t\right)<w}$ where ${\displaystyle w=\frac{c_{1}r}{K}{\left(\frac{(r+\phi )K}{2r}\right)}^2}$. It is easy to show that the solution of the first order differential inequality satisfies ${\displaystyle V\left(t\right)<\frac{w}{\phi }+\left(V\left(0\right)-\frac{w}{\phi }\right)e^{-\phi t}}$. Since ${\displaystyle \underset{t\to \infty }{lim}{e}^{-\phi t}=0}$, it is clear that $V\left(t\right)$ is uniformly bounded, which also means that all solutions of (3) are uniformly bounded. □

4. Existence and Local Stability of Equilibrium Points

4.1. The Existence of Equilibrium Points

Equilibrium points of the system (3) are determined by setting $dN/dt=0$ and $dP/dt=0$, namely

$$\begin{array}{ccc}{\displaystyle rN\left(1-\frac{N}{K}\right)-\frac{b_{1}NP}{k_1+N}}\hfill & =& 0,\hfill \\ {\displaystyle \frac{c_{1}NP}{k_1+N}+c_{2}P-eP-\frac{b_2(1-m)P^2}{k_2+(1-m)P}}\hfill & =& 0.\hfill \end{array}$$

(5)

By solving the system (5), we obtain the following four equilibrium points:

1.

The point of extinction of both populations, $E_0=(0,0)$, that always exists in ${\mathbb{R}}_{+}^2$.

2.

The prey extinction point ${\displaystyle E_1=\left(0,P_1\right)}$, where ${\displaystyle P_1=\frac{k_2(e-c_2)}{(c_2-e-b_2)(1-m)}}$. $E_1$ exists in ${\mathbb{R}}_{+}^2$ if $0<c_2-e<b_2$. This condition shows that even though prey is extinct, predator still survives as long as the rate of cannibalism is greater than the difference between the birth rate of predator due to cannibalism and the death rate of predator.

3.

The predator extinction point $E_2=(K,0)$, that always exists in ${\mathbb{R}}_{+}^2$ since $K>0$.

4.

The coexistence point $E_3=(N^{*},P^{*})$, where

$$\begin{array}{ccc}{N}^{*}\hfill & =& {\displaystyle \frac{\sqrt[3]{Q_2\pm \sqrt{Q_2^2+\frac{4}{27}{Q}_1^3}}}{\sqrt[3]{2}}-\frac{Q_1\sqrt[3]{2}}{3\sqrt[3]{Q_2\pm \sqrt{Q_2^2+\frac{4}{27}{Q}_1^3}}}+\frac{B}{3A},}\hfill \\ P^{*}\hfill & =& {\displaystyle \frac{r}{b_1}\left(1-\frac{N^{*}}{K}\right)(k_1+N^{*}),}\hfill \\ Q_1\hfill & =& {\displaystyle \frac{3AC-B^2}{3A^2},}\hfill \\ Q_2\hfill & =& {\displaystyle \frac{9ABC-2B^3-27A^{2}D}{27A^3},}\hfill \\ A\hfill & =& {\displaystyle \frac{r}{b_{1}K}(1-m)(b_2-c_1-c_2+e),}\hfill \\ B\hfill & =& {\displaystyle \frac{r}{b_1}(1-m)\left[(c_1+c_2-e-b_2)-\frac{k_1}{K}(c_1+2(c_2-e-b_2))\right],}\hfill \\ C\hfill & =& {\displaystyle (c_1+c_2-e)k_2+\frac{r{k}_1}{b_1}(1-m)\left[c_1+(2-k_1)(c_2-e)-2b_2+\frac{b_2{k}_1}{K}\right],}\hfill \\ D\hfill & =& {\displaystyle k_1\left[k_2(c_2-e)+\frac{r{k}_1}{b_1}(1-m)(c_2-e-b_2)\right],}\hfill \end{array}$$

(6)

if $b_2+e\ne c_1+c_2$. $E_3$ in (6) is obtained using Cardanos’s formula [25,26] and exists in ${\mathbb{R}}_{+}^2$ if the following conditions are met.

$$\begin{array}{c}\hfill Q_2^2+\frac{4}{27}{Q}_1^3\ge 0,\\ \hfill 0<N^{*}<K.\end{array}$$

If $b_2+e=c_1+c_2$, we have the value of $N^{*}$ and $P^{*}$ at the point of coexistence as follows.

$${\displaystyle N^{*}=\frac{-C\pm \sqrt{C^2-4BD}}{2B},P^{*}=\frac{r}{b_1}\left(1-\frac{N^{*}}{K}\right)(k_1+N^{*}),}$$

with

$$\begin{array}{ccc}B\hfill & =& {\displaystyle \frac{c_{1}r{k}_1}{b_{1}K}(1-m),}\hfill \\ C\hfill & =& {\displaystyle b_2{k}_2+\frac{r{k}_1}{b_1}(1-m)\left(k_1(c_1-b_2)-c_1+\frac{b_2{k}_1}{K}\right),}\hfill \\ D\hfill & =& {\displaystyle k_1\left[k_2(b_2-c_1)-\frac{r{c}_1{k}_1}{b_1}(1-m)\right].}\hfill \end{array}$$

The point $E_3$ exists in ${\mathbb{R}}_{+}^2$ if

$$\begin{array}{c}\hfill C^2-4BD\ge 0,\\ \hfill 0<N^{*}<K.\end{array}$$

4.2. Local Stability

Linearization around the equilibrium point is carried out so that the Jacobian matrix is obtained:

$$J\left(E\right)=\left[\begin{array}{cc}{\displaystyle r\left(1-\frac{2N}{K}\right)-\frac{b_1{k}_{1}P}{{(k_1+N)}^2}}& {\displaystyle -\frac{b_{1}N}{k_1+N}}\\ {\displaystyle \frac{c_1{k}_{1}P}{{(k_1+N)}^2}}& {\displaystyle \frac{c_{1}N}{k_1+N}+c_2-e-\frac{b_2(1-m)P\left[2k_2+(1-m)P\right]}{{(k_2+(1-m)P)}^2}}\end{array}\right].$$

(7)

The stability of the equilibrium points of the system (3) are determined by the eigenvalues of the Jacobian matrix (7) and the result is obtained in the following theorem.

Theorem 4.

The local stability of the equilibrium points of the system (3) is as follows.

(i) 

Equilibrium point $E_0(0,0)$ is always unstable (saddle node).

(ii) 

$E_1=\left(0,P_1\right)$ is locally asymtotically stable if ${\displaystyle r<\frac{b_1{P}_1}{k_1}}$ and unstable (saddle node) if ${\displaystyle r>\frac{b_1{P}_1}{k_1}}$.

(iii) 

$E_2=(K,0)$ is locally asymtotically stable if ${\displaystyle c_1<\frac{(e-c_2)(k_1+K)}{K}}$ and unstable if ${\displaystyle c_1>\frac{(e-c_2)(k_1+K)}{K}}$.

(iv) 

$E_3=(N^{*},P^{*})$ is locally asymtotically stable if ${\displaystyle k_1>K-2N^{*}}$.

Proof. 

 

(i)

By substituting $E_0=(0,0)$ to (7),

$$J\left(E_0\right)=\left[\begin{array}{cc}{\displaystyle r}& {\displaystyle 0}\\ {\displaystyle 0}& c_2-e\end{array}\right]$$

and we get the eigen values ${\lambda }_1=r$ and ${\lambda }_2=c_2-e$. Since ${\lambda }_1$ is positive, the equilibrium point of $E_0$ is always unstable. If $c_2>e$ then $E_0$ is source, while if $c_2<e$ then $E_0$ is saddle.

(ii)

The Jacobian matrix (7) for $E_1$,

$$J\left(E_1\right)=\left[\begin{array}{cc}{\displaystyle r-\frac{b_1{P}_1}{k_1}}& {\displaystyle 0}\\ {\displaystyle \frac{c_1{P}_1}{k_1}}& {\displaystyle \frac{(c_2-e)(c_2-e-b_2)}{b_2}}\end{array}\right],$$

has the eigenvalues ${\displaystyle {\lambda }_1=r-\frac{b_1{P}_1}{k_1}}$ and ${\displaystyle {\lambda }_2=\frac{(c_2-e)(c_2-e-b_2)}{b_2}}$. Based on the conditions of existence of $E_1$, then ${\lambda }_2$ is negative. $E_1$ is locally asymptotically stable if ${\displaystyle r<\frac{b_1{P}_1}{k_1}}$. While, if ${\displaystyle r>\frac{b_1{P}_1}{k_1}}$ then ${\lambda }_1>0$ and $E_1$ is unstable (saddle node).

(iii)

The Jacobian matrix for predator extinction point is

$$J\left(E_2\right)=\left[\begin{array}{cc}-r& {\displaystyle -\frac{b_{1}K}{k_1+K}}\\ 0& {\displaystyle \frac{c_{1}K}{k_1+K}+c_2-e}\end{array}\right]$$

so that the eigenvalues are ${\lambda }_1=-r$ and ${\displaystyle {\lambda }_2=\frac{c_{1}K}{k_1+K}+c_2-e}$. If ${\displaystyle c_1<\frac{(e-c_2)(k_1+K)}{K}}$, then ${\lambda }_2<0$ and causes $E_2$ locally asymptotically stable. Otherwise, ${\lambda }_2>0$ if ${\displaystyle c_1>\frac{(e-c_2)(k_1+K)}{K}}$ and $E_2$ become unstable (saddle node). Besides, by the existence condition of $E_1$, it is clear that $E_2$ is unstable if $E_1$ exists.

(iv)

By substituting ${\displaystyle E_3=(N^{*},P^{*})}$ to Jacobian matrix (7), we have the Jacobian matrix for $E_3$,

$$J\left(E_3\right)=\left[\begin{array}{cc}{\displaystyle r\left(1-\frac{2N^{*}}{K}\right)-\frac{b_1{k}_1{P}^{*}}{{(k_1+N^{*})}^2}}& {\displaystyle -\frac{b_1{N}^{*}}{k_1+N^{*}}}\\ {\displaystyle \frac{c_1{k}_1{P}^{*}}{{(k_1+N^{*})}^2}}& {\displaystyle -\frac{b_2{k}_2(1-m)P^{*}}{{(k_2+(1-m)P^{*})}^2}}\end{array}\right].$$

(8)

Since ${\displaystyle P^{*}=\frac{b_1}{k_1}\left(1-\frac{N^{*}}{K}\right)(k_1+N^{*})}$, the Jacobian (8) can be written as

$$J\left(E_3\right)=\left[\begin{array}{cc}{\displaystyle J_{11}}& {\displaystyle J_{12}}\\ {\displaystyle J_{21}}& {\displaystyle J_{22}}\end{array}\right],$$

where

$$\begin{array}{ccc}{J}_{11}\hfill & =\hfill & {\displaystyle \frac{r{N}^{*}}{k_1+N^{*}}\left(1-\frac{k_1+2N^{*}}{K}\right),}\hfill \\ J_{12}\hfill & =\hfill & {\displaystyle -\frac{b_1{N}^{*}}{k_1+N^{*}},}\hfill \\ J_{21}\hfill & =\hfill & {\displaystyle \frac{c_1{k}_{1}r}{b_1(k_1+N^{*})}\left(1-\frac{N^{*}}{K}\right),}\hfill \\ J_{22}\hfill & =\hfill & {\displaystyle -\frac{b_1{b}_2{k}_{2}r(1-m)\left(1-\frac{N^{*}}{K}\right)(k_1+N^{*})}{{\left(b_1{k}_2+r(1-m)\left(1-\frac{N^{*}}{K}\right)(k_1+N^{*})\right)}^2}.}\hfill \end{array}$$

The determinant and the trace of the matrix $J\left(E_3\right)$ are, respectively, given by

$$\begin{array}{ccc}det\left(J\right)\hfill & =\hfill & J_{11}{J}_{22}-J_{12}{J}_{21}\hfill \\ \hfill & =\hfill & {\displaystyle -\frac{b_1{b}_2{k}_2{r}^2(1-m)N^{*}\left(1-\frac{N^{*}}{K}\right)}{{\left(b_1{k}_2+r(1-m)\left(1-\frac{N^{*}}{K}\right)(k_1+N^{*})\right)}^2}\left(1-\frac{k_1+2N^{*}}{K}\right)}\hfill \\ \hfill & \hfill & {\displaystyle +\frac{c_1{k}_{1}r{N}^{*}}{{(k_1+N^{*})}^2}\left(1-\frac{N^{*}}{K}\right),}\hfill \end{array}$$

and

$$\begin{array}{ccc}trace\left(J\right)\hfill & =\hfill & {\displaystyle J_{11}+J_{22}=\frac{r{N}^{*}}{k_1+N^{*}}\left(1-\frac{k_1+2N^{*}}{K}\right)}\hfill \\ \hfill & \hfill & {\displaystyle -\frac{b_1{b}_2{k}_{2}r(1-m)\left(1-\frac{N^{*}}{K}\right)(k_1+N^{*})}{{\left(b_1{k}_2+r(1-m)\left(1-\frac{N^{*}}{K}\right)(k_1+N^{*})\right)}^2}.}\hfill \end{array}$$

Since $0<N^{*}<K$, then if $k_1>K-2N^{*}$, we get $det\left(J\right)>0$ and $trace\left(J\right)<0$, then the coexistence point is locally asymptotically stable. □

5. Global Stability

5.1. Global Stability of $E_1$

Theorem 5.

Assume that $E_1=(0,P_1)$ exists, namely if $0<c_2-e<b_2$. Then, $E_1$ is globally asymtotically stable if ${\displaystyle r<\frac{b_1{P}_1}{k_1}}$ and $K\le k_1$.

Proof. 

We consider a Lyapunov function as follows:

$$V_1(N,P)=N+\frac{b_1}{c_1}\left(P-P_1-P_{1}ln\frac{P}{P_1}\right).$$

The first order derivative of $V_1$ is given by

$$\begin{array}{ccc}{\displaystyle \frac{d{V}_1}{dt}}\hfill & =\hfill & {\displaystyle \frac{\partial V_1}{\partial N}\frac{dN}{dt}+\frac{\partial V_1}{\partial P}\frac{dP}{dt}}\hfill \\ \hfill & =\hfill & {\displaystyle rN\left(1-\frac{N}{K}\right)-\frac{b_{1}NP}{k_1+N}+\frac{b_1}{c_1}\left(\frac{P-P_1}{P}\right)\left(\frac{c_{1}NP}{k_1+N}+c_{2}P-eP-\frac{b_2(1-m)P^2}{k_2+(1-m)P}\right)}\hfill \\ \hfill & =\hfill & {\displaystyle rN\left(1-\frac{N}{K}\right)-\frac{b_{1}N{P}_1}{k_1+N}+\frac{b_1}{c_1{P}_1}\left(P-P_1\right)\left(\frac{k_2(c_2-e)P_1+k_2(e-c_2)P}{k_2+(1-m)P}\right)}\hfill \\ \hfill & =\hfill & {\displaystyle rN\left(1-\frac{N}{K}\right)-\frac{b_{1}N{P}_1}{k_1+N}-\frac{b_1}{c_1{P}_1}{\left(P-P_1\right)}^2\left(\frac{k_2(c_2-e)}{k_2+(1-m)P}\right)}\hfill \\ \hfill & \le \hfill & {\displaystyle rN\left(1-\frac{N}{K}\right)-\frac{b_{1}N{P}_1}{k_1+N}.}\hfill \end{array}$$

If ${\displaystyle r<\frac{b_1{P}_1}{k_1}}$, then ${\displaystyle b_1{P}_1>r{k}_1}$, so

$$\begin{array}{ccc}{\displaystyle \frac{d{V}_1}{dt}}\hfill & \le \hfill & {\displaystyle rN\left(1-\frac{N}{K}\right)-\frac{r{k}_{1}N}{k_1+N}}\hfill \\ \hfill & =\hfill & {\displaystyle \frac{rN}{K(k_1+N)}\left(KN-k_{1}N-N^2\right).}\hfill \end{array}$$

If $K\le k_1$, we have ${\displaystyle \frac{d{V}_1}{dt}\le 0}$. Moreover, ${\displaystyle \frac{d{V}_1}{dt}=0}$ only if $N=0$ and $P=P_1$. According to the LaSalle’s invariance principle, $E_1$ is globally asymptotically stable. □

5.2. Global Stability of $E_2$

The global stability condition $E_2$ is given by the following theorem.

Theorem 6.

$E_2$ is globally asymtotically stable if ${\displaystyle c_1<\frac{(e-c_2)(k_1+K)}{K}}$ with $0<|K-N|<\rho <k_1+K$.

Proof. 

Define a Lyapunov function as

$$V_2(N,P)=\frac{c_1}{b_1}\left(N-N_2-N_{2}ln\frac{N}{N_2}\right)+P$$

(9)

where

$$N_2=K.$$

The first order derivative of the Lyapunov function (9) is

$$\begin{array}{ccc}{\displaystyle \frac{d{V}_2}{dt}}\hfill & =\hfill & {\displaystyle \frac{\partial V_2}{\partial N}\frac{dN}{dt}+\frac{\partial V_2}{\partial P}\frac{dP}{dt}}\hfill \\ \hfill & =\hfill & {\displaystyle \frac{c_1}{b_1}\left(\frac{N-N_2}{N}\right)\left(rN\left(1-\frac{N}{K}\right)-\frac{b_{1}NP}{k_1+N}\right)+\frac{c_{1}NP}{k_1+N}+c_{2}P-eP}\hfill \\ \hfill & \hfill & {\displaystyle -\frac{b_2(1-m)P^2}{k_2+(1-m)P}}\hfill \\ \hfill & =\hfill & {\displaystyle \frac{c_1}{b_1}\left(N-N_2\right)\left(r\left(1-\frac{N}{K}\right)-\frac{b_{1}P}{k_1+N}\right)+\frac{c_{1}NP}{k_1+N}+c_{2}P-eP-\frac{b_2(1-m)P^2}{k_2+(1-m)P}}\hfill \\ \hfill & =\hfill & {\displaystyle \frac{c_{1}r}{b_1}\left(N-K\right)\left(\frac{K-N}{K}\right)+P\left(\frac{c_1{N}_2}{k_1+N}+c_2-e-\frac{b_2(1-m)P}{k_2+(1-m)P}\right)}\hfill \\ \hfill & =\hfill & {\displaystyle -\frac{c_{1}r}{b_{1}K}{\left(N-K\right)}^2+P\left(\frac{c_{1}K}{k_1+N}+c_2-e-\frac{b_2(1-m)P}{k_2+(1-m)P}\right)}\hfill \\ \hfill & \le \hfill & {\displaystyle P\left(\frac{c_{1}K}{k_1+N}+c_2-e-\frac{b_2(1-m)P}{k_2+(1-m)P}\right)}\hfill \\ \hfill & \le \hfill & {\displaystyle P\left(\frac{c_{1}K}{k_1+N}+c_2-e\right).}\hfill \end{array}$$

Assume ${\displaystyle c_1<\frac{(e-c_2)(k_1+K)}{K}}$ with $0<|K-N|<\rho <k_1+K$. Thus, we have

$$\begin{array}{ccc}{\displaystyle \frac{d{V}_2}{dt}}\hfill & \le \hfill & {\displaystyle P\left(\frac{c_{1}K}{k_1+N}-\frac{c_{1}K}{k_1+K-\rho }\right)}\hfill \\ \hfill & =\hfill & {\displaystyle c_{1}KP\left[\frac{k_1+K-\rho -(k_1+N)}{(k_1+N)(k_1+K-\rho )}\right].}\hfill \end{array}$$

(10)

In the case of $N\ge K$, $|K-N|=N-K$ and the inequality (10) can be expressed as

$$\begin{array}{ccc}{\displaystyle \frac{d{V}_2}{dt}}\hfill & <\hfill & {\displaystyle c_{1}KP\left[\frac{k_1+K-(N-K)-(k_1+N)}{(k_1+N)(k_1+K-\rho )}\right]}\hfill \\ \hfill & =\hfill & {\displaystyle c_{1}KP\left[\frac{2(K-N)}{(k_1+N)(k_1+K-\rho )}\right]\le 0.}\hfill \end{array}$$

(11)

For $N<K$, the inequality (10) becomes

$$\begin{array}{ccc}{\displaystyle \frac{d{V}_2}{dt}}\hfill & <\hfill & {\displaystyle c_{1}KP\left[\frac{k_1+K-(K-N)-(k_1+N)}{(k_1+N)(k_1+K-\rho )}\right]=0.}\hfill \end{array}$$

(12)

Based on the (11) and (12), we get ${\displaystyle \frac{d{V}_2}{dt}\le 0}$. ${\displaystyle \frac{d{V}_2}{dt}=0}$ is achieved only if $(N,P)=E_2$. Thus, $E_2$ is globally asymptotically stable. □

5.3. Global Stability of $E_3$

The global stability of the coexistence point is explained by the following theorem.

Theorem 7.

$E_3$ is globally asymptotically stable in ${\displaystyle {\mathsf{\Omega }}^{*}=\left\{(N,P)|\frac{P}{P^{*}}>\frac{N}{N^{*}}>1\right\}.}$

Proof. 

Consider a Lyapunov function

$$V_3(N,P)=N-N^{*}-N^{*}ln\frac{N}{N^{*}}+\frac{b_1}{c_1}\left(P-P^{*}-P^{*}ln\frac{P}{P^{*}}\right)$$

where $N^{*}$ and $P^{*}$ as in Equation (6). The first order derivative of $V_3$ is

$$\begin{array}{ccc}{\displaystyle \frac{d{V}_3}{dt}}\hfill & =\hfill & {\displaystyle \frac{\partial V_3}{\partial N}\frac{dN}{dt}+\frac{\partial V_3}{\partial P}\frac{dP}{dt}}\hfill \\ \hfill & =\hfill & {\displaystyle \left(1-\frac{N^{*}}{N}\right)\left(rN\left(1-\frac{N}{K}\right)-\frac{b_{1}NP}{k_1+N}\right)}\hfill \\ \hfill & \hfill & {\displaystyle +\frac{b_1}{c_1}\left(1-\frac{P^{*}}{P}\right)\left(\frac{c_{1}NP}{k_1+N}+c_{2}P-eP-\frac{b_2(1-m)P^2}{k_2+(1-m)P}\right)}\hfill \\ \hfill & =\hfill & {\displaystyle \left(N-N^{*}\right)\left(r\left(1-\frac{N}{K}\right)-\frac{b_{1}P}{k_1+N}\right)}\hfill \\ \hfill & \hfill & {\displaystyle +\frac{b_1}{c_1}\left(P-P^{*}\right)\left(\frac{c_{1}N}{k_1+N}+c_2-e-\frac{b_2(1-m)P}{k_2+(1-m)P}\right)}\hfill \\ \hfill & =\hfill & {\displaystyle \left(N-N^{*}\right)\left[r\left(1-\frac{N}{K}\right)-\frac{b_{1}P}{k_1+N}-\left(r\left(1-\frac{N^{*}}{K}\right)+\frac{b_1{P}^{*}}{k_1+N^{*}}\right)\right]}\hfill \\ \hfill & \hfill & {\displaystyle +\frac{b_1}{c_1}\left(P-P^{*}\right)\left[\frac{c_{1}N}{k_1+N}-\frac{b_2(1-m)P}{k_2+(1-m)P}-\left(\frac{c_1{N}^{*}}{k_1+N^{*}}-\frac{b_2(1-m)P^{*}}{k_2+(1-m)P^{*}}\right)\right]}\hfill \\ \hfill & =\hfill & {\displaystyle \left(N-N^{*}\right)\left[r\left(\frac{N^{*}-N}{K}\right)-\frac{b_1{k}_1(P-P^{*})}{(k_1+N)(k_1+N^{*})}\right]}\hfill \\ \hfill & \hfill & {\displaystyle +\frac{b_1}{c_1}\left(P-P^{*}\right)\left[\frac{c_1{k}_1(N-N^{*})}{(k_1+N)(k_1+N^{*})}-\frac{b_2{k}_2(1-m)(P-P^{*})}{(k_2+(1-m)P)(k_2+(1-m)P^{*})}\right]}\hfill \\ \hfill & =\hfill & {\displaystyle -\frac{r}{K}{\left(N-N^{*}\right)}^2-\frac{b_1(N-N^{*})(N^{*}P-N{P}^{*})}{(k_1+N)(k_1+N^{*})}}\hfill \\ \hfill & \hfill & {\displaystyle -\frac{b_1{b}_2{k}_2(1-m){(P-P^{*})}^2}{c_1(k_2+(1-m)P)(k_2+(1-m)P^{*})}.}\hfill \end{array}$$

In ${\displaystyle {\mathsf{\Omega }}^{*}=\left\{(N,P)|\frac{P}{P^{*}}>\frac{N}{N^{*}}>1\right\}}$, ${\displaystyle \frac{d{V}_3}{dt}<0}$ so that $E_3$ is globally asymptotically stable. □

6. Existence of Forward Bifurcation

The following theorem can be used to investigate the occurrence of the forward and backward bifurcations in model (3). The proof of the theorem can be studied in [27].

Theorem 8.

Consider the the following system with a parameter β.

$$\frac{d\overrightarrow{X}}{dt}=\overrightarrow{F}(\overrightarrow{X},\beta ),\overrightarrow{F}:{\mathbb{R}}^n\times \mathbb{R}\to {\mathbb{R}}^n.$$

(13)

Without losing generality, it is assumed that $\overrightarrow{0}$ is an equilibrium point for system (13), such that $\overrightarrow{F}(\overrightarrow{0},\beta )\equiv \overrightarrow{0}$ for all β. Assume that

• $A1$: $JM=D_{X}F(\overrightarrow{0},0)$ is the linearization matrix of system (13) around the equilibrium point $\overrightarrow{0}$ with β evaluated at 0. Zero is an eigenvalue of $JM$ and the real parts of the other eigenvalues are negative
• $A2$: Matrix $JM$ has a non-negative right eigenvector $\overrightarrow{v}$ and a non-negative left eigenvector $\overrightarrow{w}$ corresponding to the zero eigenvalue. Let $F_k$ be the k-th component of $\overrightarrow{F}$ and

  $$\begin{array}{ccc}a\hfill & =\hfill & {\displaystyle \sum _{k,i,j=1}^n{w}_k{v}_i{v}_j\frac{{\partial }^2{F}_k(0,0)}{\partial X_i\partial X_j},}\hfill \\ b\hfill & =\hfill & {\displaystyle \sum _{k,i}^n{w}_k{v}_i\frac{{\partial }^2{F}_k(0,0)}{\partial X_i\partial \beta }.}\hfill \end{array}$$

  The local dynamics of (13) around $\overrightarrow{0}$ are totally determined by a and b based on the following 4 cases.

1.

$a>0$, $b>0$. When $-1<\beta <0$, $\overrightarrow{0}$ is locally asymptotically stable, and there exists a positive unstable equilibrium. When $0<\beta <1$, $\overrightarrow{0}$ is unstable and there exists a negative and local asymptotically stable equilibrium;

2.

$a<0$, $b<0$. When $-1<\beta <0$, $\overrightarrow{0}$ is unstable. When $0<\beta <1$, $\overrightarrow{0}$ is locally asymptotically stable and there exists a positive unstable equilibrium;

3.

$a>0$, $b<0$. When $-1<\beta <0$, $\overrightarrow{0}$ is unstable, and there exists a negative unstable equilibrium. When $0<\beta <1$, $\overrightarrow{0}$ is locally asymptotically stable and a positive unstable equilibrium appears;

4.

$a<0$, $b>0$. When β changes from negative to positive, $\overrightarrow{0}$ changes its stability from stable to unstable. Correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable.

On the basis of Theorem 8, forward bifurcation occurs at $\beta =0$ if $a<0$ and $b>0$. Moreover, backward bifurcation occurs at $\beta =0$ if $a>0$ and $b>0$.

Theorem 9.

The system (3) undergoes forward bifurcation of $E_2$ when $c_1$ passes through ${\displaystyle c_1^{*}=\frac{(e-c_2)(k_1+K)}{K}}$.

Proof. 

Assume that ${\displaystyle c_1^{*}=\frac{(e-c_2)(k_1+K)}{K}}$. The Jacobian matrix of $E_2$ has eigenvalues, one of which is negative and the other one is zero if $c_1=c_1^{*}$. The right eigenvector for zero eigenvalue is

$$\overrightarrow{v}=\left[\begin{array}{c}{v}_1\\ v_2\end{array}\right]=\left[\begin{array}{c}1\\ -\frac{r(k_1+K)}{b_{1}K}\end{array}\right]v_1,$$

where $v_1$ is an arbitrary negative real number. Furthermore, the corresponding left eigenvector is

$$\overrightarrow{w}=\left[\begin{array}{c}{w}_1\\ w_2\end{array}\right]=\left[\begin{array}{c}0\\ w_2\end{array}\right],$$

where ${\displaystyle w_2=-\frac{b_{1}K}{r(k_1+K)v_1}}$ is obtained from $\overrightarrow{v}·\overrightarrow{w}=1$. Since $v_1$ is negative, $w_2$ is positive. From Theorem 8, it can be shown that

$$a=w_2{v}_1{v}_2\left[\frac{c_1{k}_1}{{(k_1+K)}^2}+\frac{c_1{k}_1}{{(k_1+K)}^2}\right]=\frac{2w_2{v}_1{v}_2{c}_1{k}_1}{{(k_1+K)}^2},$$

and

$${\displaystyle b=w_2{v}_1\frac{{\partial }^2{F}_2(E_2,c_1^{*})}{\partial N\partial c_1}+w_2{v}_2\frac{{\partial }^2{F}_2(E_2,c_1^{*})}{\partial P\partial c_1}=w_2{v}_2\frac{K}{k_1+K}.}$$

Since $a<0$ and $b>0$, based on Theorem 8, system (3) undergoes forward bifurcation at $c_1=c_1^{*}$. □

7. Numerical Simulations

In this section, numerical simulations of the model (3) are carried out using Matlab software and the 4th order Runge–Kutta method. The purposes of the numerical simulations are to show the results of the dynamic analysis that has been done regarding the stability of the equilibrium states and the forward bifurcation. We can also see the impacts of prey predation by predator, conversion of predator rate, predator cannibalism rate, and predator refuge to the system behavior. To the best of our knowledge, there are no available data related to our proposed model. Hence, the following numerical simulations are performed using hypothetical parameters.

7.1. The Impacts of Prey Predation by Predator

To see the impacts of prey predation by predators, parameter values are used as shown in

Table 1. Values of parameter to see the impacts of predation rate.

Parameter	r	K	${\mathit{b}}_1$	${\mathit{k}}_1$	${\mathit{c}}_1$	${\mathit{c}}_2$	e	${\mathit{b}}_2$	m	${\mathit{k}}_2$
Value	1	1	0.3/0.4/0.5	0.3	0.2	0.2	0.1	0.3	0.3	1

and the predation rate, $b_1$, is in the range $0.2$ to 1. By using those parameter values, we can see the effect of the increasing of $b_1$ to the solution convergence for a sufficient time in the bifurcation diagram, see Figure 1. From this figure, the solution undergoes a Hopf bifurcation at $b_1^{*}=0.3268$ in which the value is obtained numerically, and forward bifurcation at $b_1^{**}=0.4200$ in which the value is obtained analytically from the local stability condition of $E_1$. From the numerical simulation, the solution convergence is divided into three areas, i.e., the convergence to a limit cycle around $E_3$ when $0.2\le b_1<b_1^{*}$, the convergence to the existence point $E_3$ when $b_1^{*}<b_1<b_1^{**}$, and the convergence to the prey existence point $E_1$ when $b_1>b_1^{**}$. To illustrate this, phase portraits with a $b_1$ value for each area are given in Figure 2.

In the first result (Figure 2a), by choosing $b_1=0.3$, the four equilibrium points of the system exist, i.e., $E_0=(0,0)$, $E_1=(0,0.7143)$, $E_2=(1,0)$, and $E_3=(0.0571,1.1224)$. With four initial values, all of the equilibrium points are unstable and the solution leads to a limit cycle around the coexistence point $E_3$. By increasing the predation rate to $b_1=0.4$ in Figure 2b, all equilibrium points exist with $E_3=(0.0066,0.7614)$ while the others remain. The coexistence point is stable, whereas the others are unstable. In the third area, $b_1=0.5$ is selected, and the coexistence point does not exist. The solutions are stable to the prey extinction point as can be seen in Figure 2c. This makes sense because the greater the rate of prey predation by predators, the less the prey population density, and with a large enough value of $b_1$, the prey can become extinct, while the predator can still survive with the presence of conversion from cannibalism.

7.2. The Impacts of Conversion Rate of Prey Predation

To see the impacts of conversion rate of prey predation into predator birth, $c_1$, parameter values in

Table 2. Values of parameter to see the impacts of conversion rate of prey predation.

Parameter	r	K	${\mathit{b}}_1$	${\mathit{k}}_1$	${\mathit{c}}_1$	${\mathit{c}}_2$	e	${\mathit{b}}_2$	m	${\mathit{k}}_2$
Value	1	1	0.5	0.3	0.1/0.3/0.5	0.2	0.3	0.3	0.3	1

are used. The value of the parameter $c_1$ is in the interval $0\le c_1\le 0.5$ for $c_1\le b_1=0.5$. Since $c_1$ is one of the parameters that affect the stability of $E_2$, and $E_2$ will never be stable if $E_1$ exists, then the predator death rate of $e=0.3$ is chosen so that the existence condition of $E_1$ is not met. As can be seen in Figure 3, the forward bifurcation point of $E_2$ occurs at ${\displaystyle c_1=c_1^{*}=0.13}$ and the Hopf bifurcation point occurs ${\displaystyle c_1=c_1^{**}=0.4385}$. The value of $c_1^{*}$ is obtained analytically from the stability condition of $E_2$ while $c_1^{**}$ is obtained numerically.

In range $0\le c_1\le c_1^{*}$, the solution tends to predator extinction point $E_2$. It is shown in the phase portrait in Figure 4a, with the value $c_1=0.1$, the coexistence point doesn’t exist and the solution towards $E_2=(1,0)$. In other words, if the rate of predation conversion is very low, the predator can’t survive and eventually becomes extinct. Predators can preserve their population existence when $c_1\ge c_1^{*}$. The prey also still exists even though the density is decreasing. It is illustrated by the portrait phase with $c_1=0.3$ in Figure 4b, the solutions lead to $E_3=(0.6026,0.7174)$. Hereafter, when $c_1=0.5$ as in Figure 4c, $c_1>c_1^{**}$ and $E_3(0.2144,0.8082)$ is not in ${\mathsf{\Omega }}^{*}$ anymore so that the solutions goes to a limit cycle around $E_3$.

7.3. The Impacts of Predator Cannibalism Rate

The impacts of predator cannibalism rate to the system, parameter values are used as in

Table 3. Values of parameter to see the impacts of predator cannibalism rate.

Parameter	r	K	${\mathit{b}}_1$	${\mathit{k}}_1$	${\mathit{c}}_1$	${\mathit{c}}_2$	e	${\mathit{b}}_2$	m	${\mathit{k}}_2$
Value	1	1	0.3	0.3	0.2	0.12	0.02	0.2/0.28/0.35/0.5	0.3	1

. The changes in the dynamic of prey and predator population density under the changes in the predator cannibalism rate, $0.2\le b_2\le 1$, are as in Figure 5. The stable equilibrium point if $0.2\le b_2<b_2^{*}=0.2436$ is the prey extinction point $E_1$. When $b_2^{*}<b_2<b_2^{**}=0.2835$, $E_3$ is stable while the others are not. Then, as $b_2^{**}<b_2<b_2^{***}=0.3867$, all of the equilibrium points are unstable while the solutions tend to a limit cycle around $E_3$. $E_3$ returns stable for $b_1>b_2^{***}$. $b_2^{*}$ is obtained analytically from the stability condition of $E_1$ while the other bifurcation points are obtained numerically. An example of a phase portrait with different values of $b_2$ representing each of four regions bounded by the three bifurcation points can be seen in Figure 6.

In Figure 6a, with $b_2=0.2$, the equilibrium points of the system are $E_0=(0,0)$, $E_1=(0,1.4286)$, and $E_2=(1,0)$, while the coexistence point doesn’t exist. Because the rate of predator cannibalism is small, then most of the food sources for predators come from prey. Hence, the solution leads to the prey extinction point $E_1$. If $b_2$ is greater than $b_2^{*}$, then the number of prey that is predated by predators will be less so that the prey can maintain the population from extinction. In other words, if $b_2>b_2^{*}$, then prey and predator populations will always exist. When $b_2=0.28$ is chosen in Figure 6b, then $b_2^{*}<b_2<b_2^{**}$ and the coexistence point $E_3=(0.0338,1.0750)$ is stable while the others are unstable. By increasing the value of $b_2$ to 0.35, the solutions tend to a limit cycle around $E_3(0.1423,1.2645)$ (see Figure 6c). Then, $E_3$ is stable again in Figure 6d by selecting $b_2=0.5$ with the increasing prey density and the decreasing predator density.

7.4. The Impacts of Predator Refuge from Cannibalism

To see the impact of predator refuge from cannibalism (m), parameter values are used as shown in

Table 4. Values of parameter to see the impacts of predator refuge.

Parameter	r	K	${\mathit{b}}_1$	${\mathit{k}}_1$	${\mathit{c}}_1$	${\mathit{c}}_2$	e	${\mathit{b}}_2$	m	${\mathit{k}}_2$
Value	1	1	0.3	0.3	0.2	0.2	0.1	0.3	0.2/0.4/0.6	1

. Figure 7 can be interpreted that with a fairly small coefficient of refuge, there are quite a lot of available predators to be cannibalized, so that predators can eat both prey and predators. Prey and predators coexist with oscillating population density values. Then, after m throughs $m^{*}=0.3550$, they stop oscillating and the more predators take refuge from cannibalism, the more predator food sources are diverted to prey. Thus, the prey population density is getting smaller and eventually extinct when m passes $m^{**}=0.5$. The predator population density continues to increase and reaches a value of 500 when $m=1$.

To illustrate the previously described dynamics, we give the portrait phases in Figure 8a–c. $m^{*}$ is determined numerically while $m^{**}$ is determined analytically from the stability condition of $E_1$. In Figure 8a, with a small value of m, i.e., $m=0.2$, the solutions tend to a limit cycle around the coexistence point $E_3=(0.0890,1.1813)$. Then, by increasing the proportion of predators taking refuge to $m=0.4$ both of prey and predator still exist so that $E_3=(0.0274,1.0613)$ is stable with a very small value of $N^{*}$. In the last illustration in Figure 8c, the selected coefficient of refuge is quite large, i.e., $m=0.6$. The prey is extinct, while predators can survive. In this case, the coexistence point does not exist, and the prey extinction point is stable.

8. Conclusions

This paper studies the dynamical analysis of a predator–prey model incorporating predator cannibalism and refuge. For preliminaries results, we have proven the existence, uniqueness, non-negativity, and boundedness properties of the solutions. Then, we determined four types of equilibrium points with their existence, local stability, and global stability conditions. The extinction of both prey and predator point is always unstable. The prey extinction point, predator extinction point, and coexistence point are conditionally stable. Forward bifurcations occur at the prey extinction point and predator extinction point, and have been investigated. The numerical simulations support our findings and also show the occurrence of Hopf bifurcations.