Dynamical behavior of two predators–one prey model with generalized functional response and time-fractional derivative

The behavior of any complex dynamic system is a natural result of the interaction between the components of that system. Important examples of these systems are biological models that describe the characteristics of complex interactions between certain organisms in a biological environment. The study of these systems requires the use of precise and advanced computational methods in mathematics. In this paper, we discuss a prey–predator interaction model that includes two competitive predators and one prey with a generalized interaction functional. The primary presumption in the model construction is the competition between two predators on the only prey, which gives a strong implication of the real-world situation. We successfully establish the existence and stability of the equilibria. Further, we investigate the impact of the memory measured by fractional time derivative on the temporal behavior. We test the obtained mathematical results numerically by a proper numerical scheme built using the Caputo fractional-derivative operator and the trapezoidal product-integration rule.


Introduction
Mathematical modeling of the real-world phenomenon is a potent tool for predicting some ecological and biological components. The validity of this mathematical approximation depends on the model itself. The crucial component that describes the interaction between different species in a certain environment is the interaction functional. There are many types of these functionals in the literature. Each one describes a specific manner of intermingling between two species. More precisely, let us assume that x(t) and y(t) are the densities of the prey and predator populations at the time t, respectively. One of the first interaction functionals is the Holling I interaction functional ax(t)y(t), where a is the hunting (predation) rate of the prey population by a predator. This functional response is one of the most used tools in the literature. For instance, we refer to [27, 29-38, 42, 43, 49-51]. The issue with this interaction functional is the unboundedness of the consumption of the prey by a predator, which is inconsistent with the actual conditions in nature. To fix this defect, Holling [26] constructed the saturated functional response ax(t)y (t) 1+at h x(t) , where t h is the average handling time of prey by a predator. This functional resolves the unboundedness of the prey consumption by a predator. There are many other functional responses in the literature. In each case a specific manner on the intermingling between two species has been assumed.

• Hassel-Varley interaction functional ax(t)y(t)
1+bx(t)+cy(t) [6]. Here all the parameters are assumed to be positive. The reason for this great diversity in functionals is due to the variety of environmental conditions in the problem. Some of the factors that influence the selection of these parameters are the behavior of the prey and predator, and the studied area. For the last factor, many components play a crucial role such as rivers (water availability), food (for the prey), and the density of prey and predator. Overall, the functional selection depends on many factors.
The predator-prey models with three species have been attracted many researchers. In the environment the intermingling is not limited to just two populations, but interactions can be defined between more than two species in one single place. The scientists interested in this point of view have put efforts to model such complex interactions in the last few decades. We can take as an example two types of prey and one predator [15], where the predator has the capability of hunting both prey populations. Moreover, in prey-predator-superpredator models [35] the predator feeds the prey only, and the superpredator feeds both prey and predator. In some models, we study the interaction between two predators and one prey model where two types of predators are fed the same prey. Due to the intrinsic nature of the predators, there will always be a constant struggle to capture this one prey. In real situations, it is seen that one predator determines its own hunting territory. The presence of other predators in such territories is entirely unacceptable. This situation is called competition. The models in which competition is found have also received much attention in many research papers such as [2,38]. In this paper, we are interested in studying the intermingling and competition between two competitive predators on one prey with a generalized class of interaction functionals in the presence of the time-fractional derivative. In [14], it is highlighted that the fractional-time derivative explains the memory effect of a dynamical system, where the order of the derivative is called the memory rate, and the kernel of the factional derivative is called the memory function. This point of view has been applied in many disciplines such as mathematics, engineering, signal proceeding, mechanics, and especially in mathematical biology [12, 20-24, 39-41, 44]. By summarizing all the previously mentioned components let us focus on the following time-fractional formulation with a generalized consumption functional: where d δ dt ε represents Caputo's derivative in terms of time, which is defined by The conditions on the functionals f and g are defined as ( , and z(t) are the densities of prey, first predator, and other predator populations at time t, respectively. We assume that the prey population reproduces logistically with the increasing rate r and the carrying capacity of the space k, e 1 and e 2 are respectively the conversion rate of the prey biomass into the first predator population and the diversion of the prey biomass into the second predator biomass, μ 1 and μ 2 are the mortality rates of the first and second predators, respectively, β (resp., γ ) is the competition rate of the first predator with the second one (resp., of the second predator with the first one). The functionals f and g are respectively the interaction functionals for the first and second predator populations with the prey population. In the literature, there are a few papers that deal with a generalization of an interaction functional in a three-species model; we refer, for instance, to [29,37,45,[49][50][51], which give an additional motivation to our research. Furthermore, it is been also applied in understanding some epidemiological interactions; we refer, for example, to [3,4,16,17,31,36,[46][47][48]. In nature the interaction between animals is affected by many factors, such as the weather, animal nature, environmental structure, natural resources (water, food), which can affect the interaction between the three studied populations. Hence it is wise to consider a wide class of interaction functionals, which provides a wide choice of applications of the obtained results in predicting the evolution of the species. For more reading about some recent methods of modeling ecological interactions, we refer the readers to [52,53]. The applicability of model (1) is an additional motivation for us to present this paper. There are many functionals compatible with conditions (A 1 ) and (A 2 ), including Holling I-III interaction functional. Also, numerous other functionals fit with these functionals (we consider only the functional f , and the same can be assumed for the functional g): 1+at h x α [9].
This demonstrates the broad class of interaction functionality that we can consider in this study. For further details, e see [7,10,11,13]. Based on the above-mentioned mathematical and biological backgrounds, the paper is structured as follows: • In Sect. 2, we offer some tools that will be useful in dealing with the fractional operator. • Sect. 3 is devoted to the mathematical investigation of model (1), where the asymptotic behavior of the solutions is investigated. • Sect. 4 offers a numerical scheme for the fractional system (1). The graphical representations of the solution for different parameters are also involved. • The concluding section of the paper is intended to highlight the biological meanings of the acquired numerical results.

Mathematical analysis and asymptotic behavior of the solution 2.1 Equilibria of the model
In this subsection, we determine the local behavior of system (1). First, we determine the equilibria of system (1), which are the solutions of the following system: As a first remark, we deduce that system (2) has the following particular cases: (i) O = (0, 0, 0), which represents the extinction of the three populations.
(ii) 0 = (k, 0, 0), which implies the extinction of two types of predators. The point is called the predator-free equilibrium (PFE). (iii) Searching for the first predator-free equilibrium (FPFE) as 1 = (x 1 , 0, z 1 ), we insert y = 0. By replacing this result in the third equation of system (2) we get x 1 = g -1 ( μ 2 e 2 ), where g -1 is the reciprocal function of g, which exists since g is a bijective function. Substituting this last result into the first equation of (2) yields Summarizing all the results, we can conclude that FPFE Seeking for the second predator-free equilibrium (SPFE) 1 = (x 2 , y 2 , 0) by replacing z = 0 in (2). By substituting this result into the second equation of system (2) we get where f -1 is the inverse function of f , which exists since f is a bijective function. Taking this last result along with the first equation of (2), we get which is biologically relevant if x 2 < k. Summarizing all the results, we can deduce that SPFE 2 = (x 2 , y 2 , 0) exists if x 2 < k.
Remark 1 It is assumed that both functional f and g are increasing in x. From x 2 and x 1 , if lim x→+∞ f (x) = a (resp., lim x→+∞ g(x) = b), then another condition on the parameters arises, μ 1 e 1 < a (resp., μ 2 e 2 < b), which is a necessary condition for having a solution for the equation f (x) = μ 1 e 1 (resp., g(x) = μ 2 e 2 ). This feature can be seen in the Holling II-III interaction functional.
(v) Now we are in a position to seek the coexistence equilibrium, which is the positive solution of the following system: (3) Substituting (4) and (5) into the first equation of (3), we get Some straightforward calculations suggest that To guarantee at least one nontrivial intersection between two curves of the functionals F 1 and F 2 , we introduce the following assumption: which it can be rewritten as .

Asymptotic behavior of (1)
In this part, we are interested in determining the asymptotic stability of the equilibria obtained in the previous section.
Remark 3 Consider the following fractional-order system: For the time-fractional-order derivative, the concept of the local stability is very different from the first-order derivative, where in this case, we have an expansion of the stability region in comparison with the first-order derivative.
Let (x, y, z) be an equilibrium for system (1). The Jacobian matrix of system (1) at (x, y, z) is expressed as At the predator-free equilibrium the Jacobian matrix (7) reduces into The Jacobian matrix (8) has the eigenvalues ϑ 1 = -r < 0, ϑ 2 = e 1 f (k)μ 1 , and ϑ 3 = e 2 g(k)μ 2 . Hence which means that the eigenvalues ϑ i , i = 1, 2, 3, satisfy | arg(ϑ i )| > επ 2 if and only if x <x := min{x 1 , x 2 }; in this case the predator-free equilibrium is locally stable, and it is unstable for x >x. We summarize the obtained results in the following lemma.

Lemma 1
The predator-free equilibrium is locally stable if x <x and unstable if x >x. Now we analyze the linear stability of FPFE point of 1 . The Jacobian matrix corresponding to the equilibrium FPFE is evaluated as As a first look, we can deduce that ϑ 2 = e 1 f (x 1 )μ 1βz 1 is an eigenvalue of the Jacobian matrix (9). By replacing the explicit formula of z 1 we obtain ϑ 2 = e 1 f (x 1 )μ 1 -βrx 1 e 2 μ 2 (1 - βx 1 e 2 (k-x 1 ) , < 0 for r > r 1 .
Under the condition ϑ 2 > 0, we get | arg(ϑ 2 )| < επ 2 . This means that FPFE is an unstable equilibrium point. Besides, from ϑ 2 > 0 we conclude that | arg(ϑ 2 )| < επ 2 . This means that two remaining eigenvalues of the Jacobian matrix (9) determine the stability (resp., instability) of this equilibrium. Note that these significant eigenvalues are the eigenvalues of the matrix To determine the nature of the eigenvalues of the reduced matrix (10), we define the characteristic equation of (10) as To determine the sign of 1 , the following results arise.
Lemma 2 Assuming x 1 < k, we obtain the following results: Proof The proof can be easily deduced from the explicit formula of .
From Lemma 2, we can see that 1 ≤ 0 next to 2 > 0 gives | arg{ϑ 1,2 }| > επ 2 , which means that this equilibrium is stable. Now assume that 1 > 0 (which means that the second assumption of Lemma 2 holds). In this situation, (11) admits two complex roots of ϑ ± = a 1 ± a 2 i, a 1 , a 2 ∈ R. Then these roots satisfy To guarantee the stability of the equilibrium, we must have tan 2 (arg{ϑ i }) > tan 2 ( επ 2 ), i = 1, 2, which implies that Obviously, it is unstable for r > r * . Hence we resume the stability conditions for the equilibrium (x 1 , 0, z 1 ) by the following theorem.

Theorem 1 If x 1 < k, then we have:
(i) If x 1 > x 2 and r < r 1 , then the FPFE is unstable.
For ϑ 3 > 0, we obtain | arg(ϑ 3 )| < επ 2 , which means that SPFE is unstable. In fact, for ϑ 3 > 0, we get | arg(ϑ 3 )| < επ 2 . This means that the other two eigenvalues of the Jacobian matrix (12) determine the stability (resp., instability) of this equilibrium. These eigenvalues are the eigenvalues of the matrix The nature of the eigenvalues of the reduced matrix (13) can be determined through the quadratic equation where 1 = rx 1 1 - To prove the positivity of 1 , we use the following lemma.

Lemma 3 Let x 2 < k. Then:
Proof Form the explicit formula of we can determine its positivity. The proof is completed.
From Lemma 2, we can see that 1 ≤ 0 and 2 > 0 lead to | arg{ϑ 1,2 }| > επ 2 , and hence the SPFE is stable. Now assume that 1 > 0 (which means that the second assumption of Lemma 3 holds). In this case, equation (14) has complex roots ϑ ± = b 1 ± b 2 , b 1 b 2 ∈ R. Then we get To ensure the stability of SPFE, we must have tan 2 (arg{ϑ i }) > tan 2 ( επ 2 ), i = 1, 2, which implies that and it is unstable for r > r * * . We summarize the obtained results in the following theorem.
Theorem 2 Let x 2 < k. Then: (i) If x 1 < x 2 and r < r 2 , then the SPFE is unstable.
Now we are in a position to focus on studying the local behavior of the coexistence equilibrium. For this equilibrium, we suppose that assumption (H 1 ) holds. The Jacobian matrix at this equilibrium is expressed as The characteristic equation associated with Jacobian (16) is Let D = 18 2 1 0 + ( 2 1 ) 2 -4 0 3 2 -4 3 1 -27 3 0 .
Using the Routh-Hurwitz stability criterion for fractional calculus defined in [19], we get the stability conditions for the nontrivial equilibrium.

Theorem 3
The positive equilibrium is stable if one of the following assumptions holds: 3 Numerical analysis of system (1)

A numerical scheme
The main purpose of this section is to numerically solve the following fractal problem: By applying the fundamental theorem of fractional calculus on (1) we get Letting t = t n = n in (18), we arrive at Now we can approximate the function P(t, V (t)) by the following linear approximation: with the notation V i = V (t i ). By substituting Eq. (20) into (19) and applying some algebra (for more detail, see [18]) we get , n = 1, 2, . . . .
Using the numerical method presented in formula (21) to solve problem (1), we obtain the following iterative schemes: where

Example 1
In this section, we numerically investigate system (1) to ensure the obtained results in the previous sections. First, we choose the functional f to be a Holling II interaction functional For these values, we obtain the local stability of FPFE. Indeed, these values satisfy x 1 = x 2 and 1 = -1.309 < 0 (condition (i) in Lemma 2), which means that condition (ii) in Theorem 2 holds. Figure 2 We consider the following parameters: For these values, we get the local stability of the SPFE 2 = (0.0497, 0.0346, 0). More precisely, these values satisfy r < r 1 = 1.12, which means that 1 = (0.3, 0, 0.145) is unstable and satisfy the second condition of Theorem 3 (x1 > x2 and 1 = 0.148 > 0), which means that the SPFE is stable.  and multivalues of the memory rate. In this example, we realize the effect of the memory on the stability of SPFE. More precisely, this figure highlights that the memory has a huge influence on the condition of the stability obtained in (15), which shows a great agreement between the theoretical findings and numerical ones, and thus we conclude the efficiency of the obtained numerical method.

Example 2
In this example, we consider another behavior of the prey population known by herd behavior (see [1,34,56]). In fact, we consider the functionals f (x) = ax α and g(x) = bx α , where 0 < α < 1 represents the herd shape rate.
For the graphical representations, we take the following values of parameters. These values provide the local stability of the SPFE. Figure 6 We considered the same values as in Fig. 5, except α = 0.8. This leads to the instability of the SPFE and the occurrence of oscillations of the solution.   These values provide the instability of the SPFE. Figure 9 We consider the values r = 1.5, k = 10.5, a = 1.5, b = 0.5, e 1 = 0.5, e 2 = 0.5, For these values, we have studied the effect of the memory rate on the stability of the SPFE. Figure 10 The following values are used in numerical simulations r = 5.5, k = 3.5, a = 0.5, b = 1.5, e 1 = 0.5, e 2 = 0.5, μ 1 = 0.5, μ 2 = 0.5, β = 0.5, γ = 0.5, α = 0.7.
Using these values, we can investigate the effect of the memory rate on the stability of the FPFE.

Conclusions
In this research, we studied an ecological model with two predators fighting on one prey with a generalized functional response. The reason behind considering a comprehensive generalized class of functional interaction is to model the diversity in predator-prey interaction with the environment. These interactions can be affected by many factors, such as the environment and the adaptation of the three species. In the first section, we studied the existence of the equilibria of system (1), where we can have many equilibrium points next to the predator-free equilibrium. By analyzing the existence of the equilibria we ob-tained that these populations may have many scenarios. They include the extinction of three populations, two types of predators, the extinction of each population of predators, and finally the coexistence of the three populations. For the coexistence stage, we provided some conditions on the model parameters for the existence of this equilibrium. To determine which scenario will prevail, we have utilized the local asymptotic stability using the Jacobian matrix. Further, we proposed a numerical scheme to find an approximate solution of the fractional-derivative system (1), which is based on the trapezoidal productintegration rule. For the application part, we have utilized two examples. The first example models the interaction of the three populations in the case of the solitary behavior for the prey population. We obtained an excellent agreement between the mathematical results and graphical representation, which proves the accuracy of the obtained approximations by proposed iterative schemes. In fact, in the case of the extinction of the second predator, two sceneries may happen. The first case is to stabilize on the SPFE as it has been highlighted in Fig. 1 for the prey solitary behavior (PSB) of the prey population. Moreover, Fig. 5 displays the case of the prey herd behavior (PHB). The second case corresponds to the instability of the SPFE and the presence of oscillations. The numerical simulations for this case have been plotted in Fig. 7, and also Fig. 8 for PHB. The same remark is obtained for the extinction of the first predator population, presented in Figs. 2, 3, 7, in different cases (the stability and instability of SPFE). Furthermore, the effect of the memory rate measured by the time-fractional derivative order ε is established numerically in Figs. 4,9,10. In these figures, we proved that the memory rate can destabilize the equilibrium states, which demonstrates the importance of considering it in the model construction, where in the real-world, the influence of the time-fractional-order derivative can represent the influence of these two populations on the hunting. Further, we can highlight that system (1) undergoes Hopf bifurcation at FPFE and SPFE at r = r * for FPFE. Moreover, at r = r * * system undergoes Hopf bifurcation for FPFE. The oscillatory behaviors generated by these equilibria are presented in Figs. 3, 4, 6, 7, 8, 9. For more detail about the proof, we refer to the readers to [25]. Based on the obtained results in this contribution, the diversity of the ecological systems can be measured by the choice of the functional response.