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–11], fish [12, 13], and spiders [14–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–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 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 Predator–prey 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–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. Predator cannibalism is represented by the last term of the second equation in system (Equation 1). The model can be interpreted as follows: In the absence of predator, prey grows logistically with the intrinsic growth rate r and the environmental carrying capacity K. With the presence of the predator, the prey population density decreases by b1NPk1+N, where b₁ is the maximum predation rate and k₁ is the half-saturation constant. The predation rate follows Holling type-II functional response with the assumption that predators need time to catch and handle the prey. With the prey predation by predator, the predator population density increases by c1NPk1+N, where c₁ is the conversion rate of predation of prey into predator births and c₁ ≤ b₁. Predators die naturally with the death rate e. The term b2(1-m)P2k2+(1-m)P depicts the decrease in predator population density caused by cannibalism with saturated a cannibalism rate, which follows Holling type-II functional response,

The value of Equation (2) monotonically increases with the supremum b₂. (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₂) is assumed to be less than the maximum predator cannibalism rate (b₂).

The model proposed by Rayungsari et al. [21] was constructed in a system of nonlinear differential equations with the first-order 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 12. 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–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 fractional-order 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.

By applying the Caputo fractional-order derivative operator to the left-hand side of system (Equation 1), the model becomes

with α ∈ ℝ, 0 < α ≤ 1, and D*α is the α-order of Caputo fractional derivative operator defined by D*αx(t)=1Γ(1-α)∫0t(t-s)-αx(s) ds.

Since the variables in the system (Equation 3) represent the population densities, the solution of the system must be non-negative. 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 (N(0),P(0))∈ℝ+2.

Proof. Since D*α=N(r(1-NK)-b1Pk1+N), then D*αN(0)=0 if N(0) = 0. D*αN=0 means there is no change of prey population density. Consequently, N(t) = 0, ∀t > 0. Then, we prove that if N(0) > 0 then N(t) ≥ 0 for every t > 0. Suppose that the statement is wrong, so there is t^* > 0 such as

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 region [0, ∞) × Ω, where Ω={X=(N,P)∈ℝ+2:c2<e}. Then, we denote a mapping F(X) = (F₁(X), F₂(X)), where

For all X=(N,P),X¯=(N¯,P¯)∈Ω,

Since in the following discussion, it can be proved that the system solution (Equation 3) is bounded in Ω, there is a positive constant M = max{N, P}, ∀t ≥ 0. Hence, we have

with

By choosing a positive constant L = max {L₁, L₂}, 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

The Caputo derivative α-order of V satisfies,

For any positive constant φ,

If c₂ < e and by choosing 0 < φ < e − c₂, 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,

Eα[-φ(t)α]→0ast→+∞, so that,

Hence, we establish the following theorem.

THEOREM 3. All solutions of Equation (2) with initial values (N(0),P(0))∈{(x,y)∈R+2 : c2<e} are uniformly bounded

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 of the numerical simulations are to confirm the dynamics analysis results and the existence of Hopf bifurcation. Since there are no available data related to our proposed model yet, the parameter values are chosen hypothetically in Table 2.

For parameter values in Simulation 1, E₁ exists, i.e., E₁ = (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₁ 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₁.

For the second simulation, we use the same parameter values but c₁ and e (see Table 2). As a result, the existence condition for E₁ is not satisfied, so the point does not exist. It means that the prey can survive from extinction. For the predator extinction point E₂(1, 0), the stability condition in Theorem 6 is satisfied and E₂ 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₂ go to E₂. 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₃(0.1423, 1.2645), which has the eigenvalues λ^* = 0.0232 + 0.1589i and λ*¯=0.0232-0.1589i. The parameter values satisfy k₁ < K − 2N₃ − 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₃ 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₃. Meanwhile, for α > α^*, the solutions tend to limit cycle around E₃. As confirmation of the bifurcation diagram, two α values satisfying α < α^*, i.e., α = 0.8 and α = 0.89, and two α values satisfying α.α^*, i.e., α = 0.91 and α = 1, are selected to simulate the solutions of N and P with respect to time. For α = 0.8 and α = 0.89, the solutions tend to E₃. The solution with α = 0.89 oscillates longer than α = 0.8 before finally convergent to E₃. 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₂) and refuge coefficient (m). The bifurcation diagrams are shown in Figures 5, 6, respectively.

For bifurcation diagram with parameter b₂, we have three bifurcation points, i.e., b2*=0.2429, b2**=0.306, and b2***=0.372. For b2<b2*, the solutions convergent to prey extinction point E₁. It is in accordance with the analytical result since the stability condition of E₁ is satisfied. When the predator cannibalism rate is increased pass b2*, E₁ is unstable, and the solutions convergent to the coexistence point, which means the predator survive from extinction. The solutions tend to limit cycle when b2**<b2<b2***. 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^**.

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

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

AS and WMK: conceptualization. MR, AS, and ID: methodology. MR: software, data curation, writing—original draft preparation, and visualization. AS, WMK, and ID: validation, writing—reviewing and editing, and supervision. MR and AS: formal analysis. MR, AS, WMK, and ID: investigation. AS: resources and project administration. All authors have read and agreed to the published version of the manuscript.