A Fractional-in-Time Prey–Predator Model with Hunting Cooperation: Qualitative Analysis, Stability and Numerical Approximations

A prey–predator system with logistic growth of prey and hunting cooperation of predators is studied. The introduction of fractional time derivatives and the related persistent memory strongly characterize the model behavior, as many dynamical systems in the applied sciences are well described by such fractional-order models. Mathematical analysis and numerical simulations are performed to highlight the characteristics of the proposed model. The existence, uniqueness and boundedness of solutions is proved; the stability of the coexistence equilibrium and the occurrence of Hopf bifurcation is investigated. Some numerical approximations of the solution are finally considered; the obtained trajectories confirm the theoretical findings. It is observed that the fractional-order derivative has a stabilizing effect and can be useful to control the coexistence between species.


Introduction
Population dynamics are regulated by several factors: availability of resources, predation, diseases, etc. Among these factors, the interaction between prey and predators is probably the most studied one in ecology, tracing back to the works of Lotka and Volterra in the early 20th century. Since then, a number of predator-prey models have been proposed and studied (see the excellent reviews [1,2]), considering different extensions of the original one: the realistic assumption that prey populations are limited by food resources and not just by predation leads to the inclusion of terms representing carrying capacity for the prey population; the characterization of specific behaviors of the predator population results in the introduction of several functional responses for them; the complexity and dishomogeneity in the environment often requires a spatial description [3]. It is well-known that species diffusing at different rates can generate spatial patterns, observed in several biological contexts. Such Turing patterns affecting spatial predator-prey models have been deeply investigated in the literature. In addition, the inclusion of group defense has a significant impact on the dynamics of the predator-prey system. Of course, this movement is conditioned by the abundance or scarcity of the other species, so that spatio-temporal population models can include cross-diffusion terms in addition to self-diffusion ones.
Cooperation is a common behavior in many biological groups and can sensibly affect the growth rate of populations in an ecological community. Limiting our interest to the predator-prey dynamics, many mechanisms have been identified, all capable of facilitating reproduction, breeding, foraging and defense in the prey population. All of them can induce a demographic Allee effect [4]. In predators, hunting cooperation, among other interactions, can also result in Allee effects; different intensities of such a cooperative behavior impact on the survival of both species and modify the stability of the ecological system (see [5] and references therein).
Even if the interest on fractional calculus traces back at least to the 1970s (without considering the seminal suggestions given by Leibniz and Euler more than 300 years ago), in the last decades there has been an explosion of research activities on its application to several areas. Fractional order systems are not only an extension of conventional integerorder systems in mathematics but also have memory and hereditary properties that integer order systems do not have. In the last decades, there has been a huge interest in such models, due to their specific feature of describing memory effects in dynamical systems: many nonlinear models with time fractional derivatives have been considered, not only in population dynamics [6,7], but also in chemistry and biochemistry, medicine, mechanics, engineering, finance, psychology (see, among the recent surveys, [8,9]). In some situations, indeed, fractional-order models, enabling the description of the memory and hereditary properties inherent in various materials, processes and biological systems, seem more consistent with the real phenomena than integer-order models. Models with fractionalorder derivatives can take different forms depending on the system considered; among them, fractional differential equations and fractional partial differential equations are greatly applied to describe continuous systems, both deterministic and stochastic [10]. In addition, heterogeneous non-ergodic diffusion processes [11] and the effect of anomalous diffusion on a population survival have been investigated [12].
The authors considered, in previous studies, some interacting population models that explored the mentioned characteristics (limited resources, nonlinear growth, fear effect, cooperative behavior), also in the presence of spatial dishomogeneity [13][14][15][16][17]. In the above-mentioned research, systems of both ordinary and partial differential equations have been studied. Precisely, ODE descriptions for the dynamics of an intraguild predator-prey model [13] and for the spreading of waterborne diseases [15,17] have been considered and the stability of the model solutions discussed. Furthermore, in order to highlight how the spatial diffusion can both play an important role in the population evolution and lead to the formation of spatial patterns, a reaction-diffusion system modeling hunting cooperation [14] and a predator-prey system with fear and group defense [16] have been analyzed. In such researches, the effect of spatial diffusion on the stability of the equilibria has been highlighted and the conditions on the parameters ensuring the formation of patterns (stripes, spots, mixed, etc.) have been found.
In the present paper, the authors generalize the model introduced in [5], by replacing the ordinary time derivative with a fractional one so to investigate how the fractionalin-time derivative impacts the system dynamics. It is worth underlining, for the sake of completeness, that such a model has been already extended in [14] to include spatiotemporal dynamics and the related occurrence of Turing patterns. Here, instead, the authors, to better highlight the fractional derivative impact on the populations dynamics, take into account the original ODE model and consider the corresponding fractional-intime system.
Such a modeling provides challenges and ideas in many other fields of applied mathematics in which nonlinear mathematical models having a similar structure are considered [18][19][20][21][22]. After reviewing in Section 2 the main prerequisites on fractional calculus, the model is formulated in Section 3 and existence and boundedness of solutions is proved. Section 4 discusses the stability of the coexistence equilibrium, while Section 5 shows some numerical approximations by different schemes. Section 6 concludes the manuscript.

Preliminaries
In this section, we recall the fundamental definitions, concepts and results that we will use throughout the paper. For further details, we refer the reader to [23][24][25]. Definition 1. The fractional integral of order θ > 0 for a Lebesgue integrable function f : R + → R is defined by where Γ(θ) is the Euler's Gamma function.
We underline that, under natural conditions on f (t), the Caputo derivative coincides with the classical derivative whenever θ ∈ N.
Then for x(t) it holds true where E θ (·)is the Mittag-Leffler function defined by . (2) , with t > t 0 , be the system with the initial condition x(t 0 ) and x) satisfies the Lipschitz condition with respect to x, then there exists a unique solution of the above system on [t 0 , ∞) × Ω.
, with x(t 0 ) = x 0 and 0 < θ < 1, be an autonomous nonlinear fractional-order system. A point x * is called an equilibrium point of the system if it satisfies φ(x * ) = 0. This equilibrium point is locally asymptotically stable if all eigenvalues λ i of the Jacobian matrix J = ∂φ ∂x evaluated at x * satisfy the Matignon conditions | arg(λ i )| > θπ 2 .

Fractional Model with Caputo Derivative: Statement of the Problem, Boundedness, Existence and Uniqueness
We present the following fractional-order prey-predator model corresponding to the non-dimensional model introduced in [ where D θ t represents the Caputo fractional derivative, given in Definition 1, with 0 < θ < 1 and n, p are the non-dimensional variables corresponding to the prey and predator densities. All the non-dimensional parameters are positive constant. Precisely, k comprises the dimensional carrying capacity, the conversion efficiency as well as the per-capita predator mortality and attack rate, while σ is linked to the per-capita intrinsic growth rate of prey and per-capita mortality rate of predator and α is linked to the predator cooperation in hunting rate, the attack rate and the per-capita predator mortality rate. Details on both the derivation of model and the biological meaning of parameters can be found in [5,14]. To (3) we append the initial conditions

Boundedness
In this subsection, we prove the positivity and boundedness of the solution of the Theorem 1. The solution of the fractional-order system (3) and (4) is bounded in R 2 + . Moreover, the density of the population remains in a nonnegative region.
Proof. Let us define the function and ξ be a positive constant. Applying the Caputo fractional derivative on both sides in (5) and using (3), it follows that where ξ < 1. Using Lemma 1, it follows that . Hence, all the solutions to (3) starting from R 2 + are confined in the following domain A ⊆ R 2 which is positively invariant.

Existence and Uniqueness
In this subsection, we find the conditions for the existence and uniqueness of the solution to fractional-order prey-predator model (3) Denoting by F 1 and F 2 the following kernels the following theorem holds.
, then the kernels F i (t, n, p), (i = 1, 2) agree with the contraction and Lipschitz conditions in the region Θ × (t 0 , T]. Proof. Let n(t 1 ) = n 1 and n(t 2 ) = n 2 be two functions for the kernel F 1 and p(t 1 ) = p 1 and p(t 2 ) = p 2 be two functions for the kernel F 2 . Then we have and where M 1 = σ + (1 + αK)K + 2 σ k K and M 2 = |K − 1 + 2αK 2 |. Therefore, the Lipschitz conditions are satisfied for kernels F 1 and F 2 , and if 0 < M i < 1, (i = 1, 2) then M 1 and M 2 are also contractions for F 1 and F 2 , respectively. Assume that the conditions (9) and (10) hold and let us consider the kernels F 1 and F 2 . Then (7) can be written The initial conditions and the recurrence form of the model (11) are, respectively n 0 (t) = n(t 0 ), p 0 (t) = p(t 0 ) (12) and The successive difference between the terms is defined as where Taking the norm of (14), it follows from the conditions (9)-(10) that The following theorem holds.
Theorem 3. Assume that the conditions (9)-(10) hold. If then the solution of the fractional model given in (3)-(4) exists and is unique.
Proof. Let us consider n(t), p(t) bounded functions and the kernels F 1 and F 2 satisfy the Lipschitz condition. From (15) and (16), it follows that This shows the existence for the solutions. Moreover, in order to prove that (18) are solutions to (3) and (4), we consider where r 1m , r 2m are the remaining terms. We show that the terms in (19) satisfy r 1∞ (t) → 0 and r 2∞ (t) → 0. Now, we consider the conditions On using recursive techniques, we get For m → ∞, it follows that r 1m (t) → 0. In a similar way, we conclude that r 2m (t) → 0. In order to prove the uniqueness of the solution to the model (3) and (4), let us suppose there exists another solution of the systemn(t) andp(t). Then In view of the Lipschitz condition, it follows that and hence In view of (17), it follows that n(t) −n(t) = 0 and p(t) −p(t) = 0.

Stability Analysis
The equilibrium points of system (3) are obtained by considering According to [5,14], besides the trivial equilibrium E 0 ≡ (0, 0), the system admits the predator-free equilibrium E b ≡ (k, 0) and the coexistence equilibrium E * = (n * , p * ) < k ≤ 1 and σ > 1 α there exist two coexistence equilibria. In all other cases, no coexistence equilibria are admissible.
In the following, we will discuss the stability properties of the coexistence equilibrium by using the linearization method.
The Jacobian matrix of system (3) evaluated at the equilibrium point E * = (n * , p * ) is given by Now, the characteristic equation of the Jacobian matrix is λ 2 − Iλ + A = 0 and the roots of the characteristic equation are where I and A are the trace and determinant of the Jacobian matrix J(E * ) given by The stability analysis of E * will be divided in the following four cases: (1) If I 2 − 4A ≥ 0 and I < 0, then λ i (i = 1, 2) are negative real numbers. Hence, The coexistence equilibrium is asymptotically stable.
(4) When I = 0 and In this case, the equilibrium is asymptotically stable.
A Hopf bifurcation occurs when a pair of complex eigenvalues of the Jacobian matrix at an equilibrium point exists and the stability of the equilibrium changes from stable to unstable when a bifurcation parameter crosses a critical value. We can choose the order of derivation θ to be the bifurcation parameter and by using the following well-known result we find the conditions for the Hopf bifurcation to appear. Lemma 4. [26] When bifurcation parameter θ passes through the critical value θ * ∈ (0, 1), fractional-order system (3) undergoes a Hopf bifurcation at the equilibrium point, if the following conditions hold 1.
the corresponding characteristic equation of system (3) has a pair of complex conjugate roots We find the conditions for the model (3) tp undergo a Hopf bifurcation at equilibrium point E * ≡ (n * , p * ) when the order of derivation passes through the critical value θ * = 2 Let us assume that I 2 − 4A < 0 and I > 0; let the critical value θ * = 2 Let us define γ = I 2 and ω = |I 2 − 4A| 2 . Then, it follows that the eigenvalues are a pair of complex conjugate λ 1,2 = γ ± iω with γ > 0. In addition, let m(θ) = θ π 2 − min 1≤i≤2 | arg λ i |, then, it follows that Therefore, we can conclude that a Hopf bifurcation of (3) will appear at E * .

Numerical Solution
Even if the interest on fractional calculus traces back at least to the 1970s (without considering the seminal suggestions given by Leibniz and Euler more than 300 years ago), in the last decades there has been an explosion of research activities on its application to several areas. Such a growing interest for fractional-order models had led to the development of specific algorithms devoted to the numerical approximation of the solution to fractional-order differential problems. Several solvers have been proposed in the last years [27][28][29], all trying to balance efficiency and accuracy while guaranteeing reliability of the numerical approximation. It is well-known, indeed, that the persistent memory of fractional-order operators reflects on the numerical evaluation of the solution: as a natural consequence, at each new timestep all the past history of the solution is to be considered. The number of involved steps increases with time, and so does the computational burden.
In a sense, numerical methods for fractional-order differential systems are then naturally multi-step. They generalize some of the ODE methods, but with significant differences in both complexity and accuracy. We follow here the excellent survey [30] and consider some of the schemes reported there. The simplest schemes are derived by approximating the integrand function x(t) = (n(t), p(t)) in (7) by a piecewise polynomial and proceeding to its exact integration. This leads to "rectangular" (explicit or implicit) product integration formulas when a zero-order approximation of the integrand function in each subinterval is assumed, or "trapezoidal" formulas where first-order approximation is chosen. Of course, implicit formulas are more stable, but they require solving the nonlinear equations that generally arise. To avoid this additional burden, a predictor-corrector approach can be useful.
The convergence order of the rectangular formulas is one: the distance between the numerical approximation and the exact solution in any time point t n decreases linearly with the timestep, and this result replicates the well-known one for the ODE formula. Unfortunately, this is not always true for the trapezoidal rule: its convergence order is limited by the minimum between 1 + θ and 2 [28], so that for 0 < θ < 1 the rate of the corresponding ODE method cannot be reached. Similarly, higher order formulas for fractional systems do not lead to significant improvements in accuracy and convergence rate, and they are not worthy to be considered. As an alternative to product integration formulas, Fractional Linear Multistep methods have been proposed to generalize the standard multistep methods for ODEs. They are very robust and can reach higher convergence order, but their convolution weights are not known explicitly in advance and have to be evaluated. This computation can however be performed very efficiently by FFT derived algorithms.
Numerical approximations of the solution of (3) can then be obtained by applying several different schemes. Even developed from the same ODE methods, they present distinguishing characteristics that deserve to be investigated. For this reason, we considered and compared the following schemes: • Implicit Rectangular Product Integration rule (PI1) • Implicit Trapezoidal Product Integration rule (PI2) with a (θ) • Predictor-Corrector Product Integration rule (PI12) • Fractional Backward Differentiation Formula (FLMM2) All these numerical schemes have been implemented in Matlab routines as given in [30,31]. For the problem at hand, the performance of these schemes can be evaluated only qualitatively, due to the lack of an exact (closed form) solution, by comparing their behavior for decreasing timesteps; however, we refer the reader to the above cited literature for an exhaustive assessment of their results on several test cases.

Preliminary Assessment
As a preliminary test, we checked the results of all methods in reproducing trajectories when θ = 1, because these results can be directly compared with the reference solutions given by classical ODE solvers. Even if the considered methods are all devised for fractionalorder systems, they are expected to reproduce a good numerical approximation of the solution also for the integer order case. When the parameter settings (σ = 3, α = 0.3, k = 5) lead to a stable coexistence equilibrium, all methods are able to reproduce the expected trajectories: the left panel of Figure 1 shows the convergence of all numerical approximations of n(t) towards the equilibrium value n * ≈ 0.66. On the other hand, when Hopf instability occurs (σ = 3, α = 10, k = 0.8), the simplest PI1 method shows a sensible damping of the oscillations around the equilibrium n * ≈ 0.188, as shown in the right panel of the same Figure, while the other methods' trajectories are practically indistinguishable from the reference one (as reported in [14]).

Stabilizing Effect of the Persistent Memory
As a second test case, we consider a parameter setting for which the corresponding ODE system shows instability: if we choose as before σ = 3, α = 10, k = 0.8, simply considering a slightly lower value for θ (θ = 0.9) has a powerful stabilizing effect. As Figure 2 shows, all the numerical approximations agree in a fast damping of the oscillations for both the variables. Again, the damping is more evident in the lower order method PI1, while the other considered schemes agree in accurately reconstructing the system trajectories. To further analyze the different performance of the considered numerical methods, we present in the following Figure 3 Figure 1), the trajectories of the fractional system rapidly stabilizes for both variables (n(t) shown in the left panel, p(t) in the right one). The lowest order approximation PI1 clearly shows more damped oscillations, while the higher order approximations are all in excellent agreement.

Hopf Instability for the Fractional System
As a third example, we show how the instability can still occur for the fractional system with a suitable choice of the parameters. We start by considering a parameter setting for which, as stated in Section 4, the conditions for the Hopf instability to appear are met: we set σ = 3, α = 3.5, k = 5. In this case, it is I > 0 and I 2 − 4A < 0. The eigenvalues of the Jacobian are a couple of complex conjugate numbers with a positive real part so that for θ > θ * = arctan 2/π √ 4A − I 2 /I ≈ 0.89 the coexistence equilibrium E * ≈ (0.27, 0.77) loses its stability. The following Figure 5 shows the trajectories of (n(t), p(t)) when θ = 0.92 as reconstructed by the more accurate numerical algorithms and the corresponding phase plan portrait, clearly showing the appearance of a limit cycle.

Conclusions
Fractional calculus, implying non-locality and memory effects, allows the description of numerous phenomena in a wide variety of scientific domains. Fractional-order operators have been proved to be a very powerful modeling instrument to represent a variety of processes and biological systems. In this framework, we studied in this paper the dynamic behavior of a fractional-in-time prey-predator model with hunting cooperation. The existence and uniqueness of a non-negative solution has been proved and the local stability of the coexistence equilibrium point has been analyzed by using the linearization technique. The conditions for the occurrence of a Hopf bifurcation have been found. Finally, numerical simulations have been presented to confirm by some selected examples how the order of derivation θ affects the dynamical behavior of prey and predator density. The findings of the present study, in our opinion, can provide hints for the investigation of the dynamics of other processes describing real-world phenomena. As a final consideration, we outline that several authors investigated pattern formation in the related spatial fractional-order systems. Recently, Yin and Wen [32] have shown that fractional-derivative systems can result in persistent spatial patterns even though their ODE counterpart does not induce any steady pattern. Such mechanisms of pattern formation along with the ones induced by anomalous diffusion suggest further directions of research in spatial generalizations of the considered model.