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

: 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 ﬁrst 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 deﬁning 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 ﬁndings, we perform some numerical simulations.


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 predatorprey 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: where N ≥ 0 and P ≥ 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 Sections 5 and 6. The numerical simulations and the intrepretations are carried out in Section 7. Finally, we draw some concluding remarks in Section 8.

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

Existence and Uniqueness
In this section, we show the existence and uniqueness of solution of the system (3) in Ω × [0, T] where Ω = {X = (N, P) ∈ R 2 + } and T < ∞. For this aim, we denote F(X) = (F 1 (X), F 2 (X)), where such that the system (3) can be written as It can be verified that F i , ∂F i ∂N , and ∂F i ∂P , i = 1, 2 are continuous in Ω. Based on a lemma in [24] (p. 71), F(X) is locally Lipschitz on Ω. Consequently, using the fundamental existenceuniqueness theorem (see [24] (p. 74)), we obtain the following theorem. Theorem 1. For the system (3) with any non-negative initial condition N(0) ≥ 0 and P(0) ≥ 0, there exists T > 0 so that the system (3) has a unique solution defined in Ω.

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 nonnegative by the following theorem. Theorem 2. All solutions of (3) with initial values (N(0), P(0)) ∈ R 2 + are non-negative.
It means that the prey population density N doesn't change from the beginning to the next. Subsequently, it is assumed that N(0) > 0. If N(t) ≥ 0 for every t ≥ 0 is not true, then there is t 1 > 0 such that N(t) > 0 for 0 < t < t 1 , N(t) = 0 for t = t 1 and N(t) < 0 for t > t 1 . From (3), we obtain: Thus, there is no change in the population density of N when t = t 1 . This contradicts the statement that N(t) < 0 for t > t 1 . Therefore, the previous assumption is false, which means N(t) ≥ 0 for every t > 0. In the same way, it can be proved that P(t) ≥ 0 for every t > 0.

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(0), P(0)) ∈ R 2 + are uniformly bounded.

Proof. Consider a function defined by
where p 1 , p 2 > 0. V(t) has the first derivative: If p 1 = c 1 and p 2 = b 1 , then For any positive constant ϕ, it holds that By choosing ϕ < e − c 2 , we have It is easy to show that the solution of the first order differential inequality satisfies V(t) < w is uniformly bounded, which also means that all solutions of (3) are uniformly bounded.

The Existence of Equilibrium Points
Equilibrium points of the system (3) are determined by setting dN/dt = 0 and dP/dt = 0, namely 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 R 2 + .

2.
The prey extinction point 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 R 2 (6) is obtained using Cardanos's formula [25,26] and exists in R 2 + if the following conditions are met.
If b 2 + e = c 1 + c 2 , we have the value of N * and P * at the point of coexistence as follows.

Local Stability
Linearization around the equilibrium point is carried out so that the Jacobian matrix is obtained: 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. (i) Equilibrium point E 0 (0, 0) is always unstable (saddle node). (ii) Proof.
(i) By substituting E 0 = (0, 0) to (7), and we get the eigen values λ 1 = r and λ 2 = c 2 − e. Since λ 1 is positive, the equilibrium point of E 0 is always unstable. If c 2 > e then E 0 is source, while if The Jacobian matrix (7) for E 1 , Based on the conditions of existence of E 1 , then λ 2 is negative. E 1 is locally asymptotically stable The Jacobian matrix for predator extinction point is so that the eigenvalues are λ 1 = −r and λ 2 = , then λ 2 < 0 and causes E 2 locally asymptotically stable.
Besides, by the existence condition of By substituting E 3 = (N * , P * ) to Jacobian matrix (7), we have the Jacobian matrix for E 3 , (8) can be written as The determinant and the trace of the matrix J(E 3 ) are, respectively, given by Since 0 < N * < K, then if k 1 > K − 2N * , we get det(J) > 0 and trace(J) < 0, then the coexistence point is locally asymptotically stable.

Global Stability
Then, E 1 is globally asymtotically stable if r < b 1 P 1 k 1 and K ≤ k 1 .

Proof.
We consider a Lyapunov function as follows: The first order derivative of V 1 is given by If K ≤ k 1 , we have dV 1 dt ≤ 0. Moreover, dV 1 dt = 0 only if N = 0 and P = P 1 . According to the LaSalle's invariance principle, E 1 is globally asymptotically stable.

Global Stability of E 2
The global stability condition E 2 is given by the following theorem.
Proof. Define a Lyapunov function as where N 2 = K.
The first order derivative of the Lyapunov function (9) is In the case of N ≥ K, |K − N| = N − K and the inequality (10) can be expressed as For N < K, the inequality (10) becomes Based on the (11) and (12), we get Thus, E 2 is globally asymptotically stable.

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 Ω * = (N, P)| P P * > N N * > 1 .

Proof. Consider a Lyapunov function
where N * and P * as in Equation (6). The first order derivative of V 3 is In Ω * = (N, P)| P P * > N N * > 1 , dV 3 dt < 0 so that E 3 is globally asymptotically stable.

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 β.
Without losing generality, it is assumed that 0 is an equilibrium point for system (13), such that F( 0, β) ≡ 0 for all β.
The local dynamics of (13) around 0 are totally determined by a and b based on the following 4 cases.
On the basis of Theorem 8, forward bifurcation occurs at β = 0 if a < 0 and b > 0. Moreover, backward bifurcation occurs at β = 0 if a > 0 and b > 0.

Theorem 9. The system (3) undergoes forward bifurcation of E 2 when c 1 passes through c
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 where v 1 is an arbitrary negative real number. Furthermore, the corresponding left eigenvector is Since v 1 is negative, w 2 is positive.
From Theorem 8, it can be shown that Since a < 0 and b > 0, based on Theorem 8, system (3) undergoes forward bifurcation at c 1 = c * 1 .

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.

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 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 ≤ 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.

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 are used. The value of the parameter c 1 is in the interval 0 ≤ c 1 ≤ 0.5 for c 1 ≤ 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 c 1 = c * 1 = 0.13 and the Hopf bifurcation point occurs 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 ≤ c 1 ≤ 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 ≥ 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 Ω * anymore so that the solutions goes to a limit cycle around E 3 . Table 2. Values of parameter to see the impacts of conversion rate of prey predation.

The Impacts of Predator Cannibalism Rate
The impacts of predator cannibalism rate to the system, parameter values are used as in Table 3. The changes in the dynamic of prey and predator population density under the changes in the predator cannibalism rate, 0.2 ≤ b 2 ≤ 1, are as in Figure 5. The stable equilibrium point if 0.2 ≤ 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.   Table 3: (a) N state and (b) P state.

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. 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.    Table 4 and (a) m = 0.2, (b) m = 0.5, and (c) m = 0.6.

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.