Dynamics analysis of a predator–prey fractional-order model incorporating predator cannibalism and refuge

In this article, we consider a predator–prey interaction incorporating cannibalism, refuge, and memory effect. To involve the memory effect, we apply Caputo fractional-order derivative operator. We verify the non-negativity, existence, uniqueness, and boundedness of the model solution. We then analyze the local and global stability of the equilibrium points. We also investigate the existence of Hopf bifurcation. The model has four equilibrium points, i.e., the origin point, prey extinction point, predator extinction point, and coexistence point. The origin point is always unstable, while the other equilibrium points are conditionally locally asymptotically stable. The stability of the coexistence point depends on the order of the Caputo derivative, α. The prey extinction point, predator extinction point, and coexistence point are conditionally globally and asymptotically stable. There exists Hopf bifurcation of coexistence point with parameter α. The analytic results of stability properties and Hopf bifurcations are confirmed by numerical simulations.


. Introduction
Predator-prey interaction, as the basis of the food chain, is among the most essential ecological issues. In numerous published research, mathematical models have been developed to explain the dynamics of Predator-prey interaction, such as by incorporating social behavior [1,2], age structure [3,4], ratio-dependent functional response [5,6], harvesting [7,8], and so on. The Predator-prey model is still being developed by considering many factors that occur in nature. Cannibalism, the consuming of the same species in whole or in part, is one of its most intriguing aspects since many animals in nature exhibit cannibalistic behaviors, such as carnivore mammals [9][10][11], fish [12,13], and spiders [14][15][16]. Cannibalism may provide adaptive advantages such as exploiting conspecifics as a food source or eliminating possible competitors [17].
Some researchers have investigated the mathematical model involving cannibalism [18][19][20][21]. Kang et al. [18] studied a single-species cannibalism model with stage structure. The model studied is a dynamic system of one population such an age structure that divides the population into two classes, i.e., eggs and an adult class consisting of larvae, pupae, queen insects, worker insects, and other types. Zhang et al. [19] analyzed predator-prey .
/fams. . 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 and adult predators. The birth rate of juvenile predators is assumed to be proportional to the number of adult predators and follows the Malthus growth model. Predation of prey and juvenile predators by adult predators follows the type-I Holling functional response. Meanwhile, Deng et al. [20] studied a two-dimensional predator-prey model with predator cannibalism. Aside from cannibalism, another interesting Predatorprey phenomenon to investigate is prey hiding behavior from predator captures and attacks. This is known as refuge behavior in the context of ecology. The mathematical model of Predator-prey interaction with prey refuge has also piqued the interest of researchers [21][22][23][24][25]. Rayungsari et al. [21] modified model proposed by Deng et al. [20] by adding the assumption that there is a refuge in the cannibalized predator population, as much as mP. Moreover, it is also assumed that predators need time to catch and handle the prey, so that the rate of prey predation follows the Holling type-II functional response. The Predator-prey model incorporating predator cannibalism and refuge proposed by Rayungsari et al. [21] is as follows: where N ≥ 0 and P ≥ 0 represent prey density and predator density, respectively. The parameters of system (Equation 1) are positive constants described in Table 1 The value of Equation (2) monotonically increases with the supremum b 2 . (1 − m)P is the amount of the available predator to be cannibalized, as m is the coefficient of refuge. The conversion rate of cannibalism into predator birth (c 2 ) is assumed to be less than the maximum predator cannibalism rate (b 2 ). The model proposed by Rayungsari et al. [21] was constructed in a system of nonlinear differential equations with the firstorder derivative, where the change of population density at any time depends on the current population density instantaneously. Whereas in reality, the current condition is also affected by the history of all previous conditions, which is called the memory effect [26]. The phenomenon or systems that have memory and genetic characteristics can be described by a fractional differential system [27]. The definition of fractional-order derivative was first introduced by Liouville [28] motivated by L'Hôpital and Leibniz's critical thinking on derivatives of order 1 2 . Liouville's definition was modified by Riemann by applying a direct generalization of the Cauchy formula and named Riemann-Liouville fractional derivative operator [29]. The fractional-order derivative concept by Liouville and Riemann utilizes Euler's study of fractional integration, which led him to construct the Gamma function as generalization of the factorial concept for fractional numbers [30]. In 1967, Michele Caputo modified the Riemann-Liouville operator so that when solving differential equations, no initial conditions are required. The definition of the modified operator is named by Caputo fractional-order derivative operator. Predator-prey models using Caputo-type fractional derivatives have been widely studied recently [24,[31][32][33]. Hence, in this article, we modify and analyze the Predator-prey model incorporating predator cannibalism and refuge in Rayungsari et al. [21] by applying the Caputo fractionalorder derivative operator.
This article is organized as follows. In Section 2, model development and basic properties are described. The basic properties consist of verification of the non-negativity, existence, uniqueness, and boundedness of solutions of the model. In Section 3, the results of dynamical analysis are presented. The results consist of the existence and stability of equilibrium points. Both local and global stability are investigated, while the analyzed bifurcation is the Hopf bifurcation. In Section 4, the numerical simulations and intrepretations are carried out to confirm the analytical results. Finally, in Section 5, we draw some concluding remarks.

. Model development and basic properties
By applying the Caputo fractional-order derivative operator to the left-hand side of system (Equation 1), the model becomes with α ∈ R, 0 < α ≤ 1, and D α * is the α-order of Caputo fractional derivative operator defined by D α * x(t) = 1 Since the variables in the system (Equation 3) represent the population densities, the solution of the system must be nonnegative. The solution of system (Equation 3) is guaranteed to be non-negative by the following theorem.

THEOREM 1. All solutions of Equation (3) are non-negative for any initial values
From Equations (3), (4), we get that D α * N = 0, t = t * . Thus, there is no change in the population density of N when t = t * . From the prior statement, N(t) = 0, t = t * , so that N(t) = 0, t > t * . This contradicts the statement that N(t) < 0 for t > t * . Therefore, N(t) ≥ 0 for all t > 0 is correct. In the same way, it can be proved that P(t) ≥ 0 for every t > 0.
Next, we show the existence and uniqueness of solution of the system (Equation 3) using Theorem 3.7 in Li et al. [34]. Consider a For all X = (N, P),X = (N,P) ∈ , Since in the following discussion, it can be proved that the system solution ( By choosing a positive constant L = max {L 1 , L 2 }, we get Based on Equation (6), the function F(X) satisfies the Lipschitz condition so that there exist a unique solution X(t) of the system (Equation 3) with any initial value of X(0) = (N(0), P(0)). Thus, we derive the following theorem. THEOREM 2. For the system (Equation 3) with any non-negative initial condition (N(0), P(0)) ∈ , there exist a unique solution X(t) ∈ .
Next, due to the limited carrying capacity of the prey and predator resources, the size of both populations in the system (Equation 3) must be limited. Consider a function defined by . /fams. .
The Caputo derivative α-order of V satisfies, For any positive constant ϕ, If c 2 < e and by choosing 0 < ϕ < e − c 2 , we get Based on Equation (7), Generalized Mean Value Theorem in Odibat and Shawagfeh [35], and Lemma 6.1 (Fractional Comparison Principle) in Li et al. [34], we get that, Hence, we establish the following theorem. In the similar way as in Rayungsari et al. [21], the system (Equation 3) has four equilibrium points, namely E 0 = (0, 0), . If b 2 + e = c 1 + c 2 , then N 3 and P 3 in E 3 is obtained from the solution of a cubic equation using the Cardano's formula [36,37], i.e., with Whereas, if b 2 + e = c 1 + c 2 , we have the value of N 3 and P 3 as follows: Two of the equilibrium points need existence conditions. E 1 exists in R 2

. . Local stability
Local stability of the equilibrium points of Equation (3) are determined by the arguments of the eigenvalues of Jacobian matrix and applying Matignon Local Stability Theorem in Petras [38]. Suppose that E * is an equilibrium point of system (Equation 3). Based on Matignon Local Stability Theorem, E * is local asymptotically stable if all of the eigenvalues λ j of the Jacobian matrix, that satisfies | arg(λ j )| > απ 2 .
Proof. With the same way, we get the Jacobian matrix for E 2 as follows: The eigenvalues are λ 1 = −r and λ 2 = | arg(λ 2 )| = 0 < απ 2 , and E 2 is an unstable saddle node.
For existence point E 3 (N 3 , P 3 ), the Jacobian matrix is as follows: where Thus, the eigenvalues are obtained from the following quadratic equation. where then Suppose that d is the discriminant of Equation (16), i.e., The cases are divided into two parts, those are for d ≥ 0 and for d < 0.
Hence, we establish the following theorem.  (14). E 3 = (N 3 , P 3 ) is locally asymptotically stable if one of the following conditions are satisfied.
with a is as in Equation (19).

. . Global stability
Next, we investigate the global stability of E 1 , E 2 , and E 3 . For this aim, we use the help of Lemma 3.1 in Vargas-De-Leon [39] and Generalized Lasalle Invariance Principle in Huo et al. [40].
For prey extinction point E 1 (0, P 1 ), we consider a Lyapunov function, The Caputo derivative α-order of V 1 is as follows: For N > 0, if K ≤ k 1 , then D α * V 1 ≤ 0 and according to Generalized Lasalle Invariance Principle [40], E 1 is globally asymptotically stable. We write the global stability conditions of E 1 in the following theorem. THEOREM 8. If E 1 = (0, P 1 ) exists, then E 1 is globally asymtotically stable if r < b 1 P 1 k 1 and K ≤ k 1 .
Then, we construct a Lyapunov function as follows: for E 2 (K, 0). We have, Suppose that e > c 1 K k 1 + c 2 . Thus, we have, We get that D α * V 2 ≤ 0, ∀(N, P) ∈ R 2 + . Hence, E 2 is globally asymptotically stable with the condition as in the following theorem. THEOREM 9. E 2 is globally asymtotically stable if e > c 1 K k 1 + c 2 .
To investigate the global stability of coexistence point, we consider a Lyapunov function where N 3 and P 3 as in Equation (9). The α-order derivative of V 3 satisfies .
Consider a domain * = (N, P) P P 3 > N N 3 > 1 . Then, D α * V 3 < 0 and E 3 is globally asymptotically stable in * . Hence, we derive the following theorem. THEOREM 10. E 3 is globally asymptotically stable in the domain * = (N, P) . . Existence of Hopf bifurcation THEOREM 11. If d < 0 and k 1 < K − 2N 3 − aK with a is given in Equation (19), then E 3 undergoes Hopf bifurcation when the order of Caputo derivative, α, pass α * with and λ * is an eigenvalue of E 3 .

. Numerical simulations
In this section, numerical simulations of the model (Equation 3) are carried out using Matlab software and the predictor-corrector scheme, which is introduced by Diethelm et al. [42]. The purposes  Table 2.
For parameter values in Simulation 1, E 1 exists, i.e., E 1 = (0, 0.7143) and the local stability condition in Theorem 5 is satisfied. We conduct numerical simulations with several values of α. The results in Figure 1 show that the solutions tend to the prey extinction point for all α values chosen. This is consistent with the analytical results since the Jacobi matrix eigenvalues are negative real numbers, which involve E 1 always asymptotically stable with the selected parameter values for any order derivative of the α ∈ (0, 1]. However, we can see a difference in the solution's behavior for each α. The lower the α value, the slower the solution reaches E 1 .
For the second simulation, we use the same parameter values but c 1 and e (see Table 2). As a result, the existence condition for E 1 is not satisfied, so the point does not exist. It means that the prey can survive from extinction. For the predator extinction point E 2 (1, 0), the stability condition in Theorem 6 is satisfied and E 2 is asymptotically stable for any fractional order of α ∈ (0, 1]. It fits the numerial simulation results in Figure 2. Represented by some values of α, we can see that the solutions with initial value close to E 2 go to E 2 . With a greater α value, the solution will reach the predator extinction point faster. Next, we demonstrate the effect of the derivative order on the behavior of the solution, with 0.8 ≤ α ≤ 1. The parameter values in the last column of Table 2 were chosen. With those parameter values, the coexistence point exists, i.e., E 3 (0.1423, 1.2645), which has the eigenvalues λ * = 0.0232 + 0.1589i andλ * = 0.0232 − 0.1589i. The parameter values satisfy k 1 < K − 2N 3 − aK and the discriminant of the quadratic equation of the eigenvalues is negative, i.e., d = −0.1010. Based on the Theorem 7, the stability of E 3 is determined by the argument of the order derivative α. The threshold is α * = 0.9077, which satisfies α * < 2 π λ * −λ * λ * +λ * . From the bifurcation diagram in Figure 3, we can see that for α < α * , the solutions tend to E 3 . Meanwhile, for α > α * , the finally convergent to E 3 . Meanwhile, for α = 0.91 and α = 1, each solution convergent to a limit cycle. The amplitude of the limit cycle solution with α = 1 is greater than α = 0.91. Numerical simulations in Figures 3, 4 show the existence of Hopf bifurcation in system (3) with α as bifurcation parameter. In addition, the system also undergoes one-parameter Hopf bifurcation with other bifurcation parameters such as cannibalism rate (b 2 ) and refuge coefficient (m). The bifurcation diagrams are shown in Figures 5, 6, respectively.
For bifurcation diagram with parameter b 2 , we have three bifurcation points, i.e., b * 2 = 0.2429, b * * 2 = 0.306, and b * * * 2 = 0.372. For b 2 < b * 2 , the solutions convergent to prey extinction point E 1 . It is in accordance with the analytical result since the stability condition of E 1 is satisfied. When the predator cannibalism rate is increased pass b * 2 , E 1 is unstable, and the solutions convergent to the coexistence point, which means the predator survive from extinction. The solutions tend to limit cycle when b * * 2 < b 2 < b * * * 2 . For bifurcation diagram with parameter m, we have two bifurcation points, i.e., m * , m * * . For m < m * , the solutions convergent to coexistence point. The solutions tend to limit cycle in the refuge coefficient range m * < m < m * * .

. Conclusion
A first-order system of Predator-prey interaction incorporating predator cannibalism and refuge is modified by applying Caputo fractional-order derivative operator. We verify the non-negativity, existence, uniqueness, and boundedness of the model solution. The local and global stability of equilibrium points are then examined. In addition, the existence of Hopf bifurcation is investigated. There are four equilibrium points in the model: the origin point, the prey extinction point, the predator extinction point, and the coexistence point. Since the eigenvalues are real numbers, the first three equilibrium points have the same local stability as the first-order system. However, the local stability of the coexistence point differs from that of the first-order system. The coexistence point with positive real-part eigenvalues is locally .
/fams. . asymptotically stable in the modified system as long as the absolute of the eigenvalue arguments are bigger than απ 2 . Even though it is based on different theories, the global stability properties of all equilibrium points are identical to those in the firstorder system. Under certain conditions, the Hopf bifurcation exists for the coexistence point. Numerical simulations confirmed the analytical results of stability properties and the existence of Hopf bifurcation.

Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.