Next Article in Journal / Special Issue
An Information-Theoretic Framework for Optimal Design: Analysis of Protocols for Estimating Soft Tissue Parameters in Biaxial Experiments
Previous Article in Journal
Trees in Positive Entropy Subshifts
Previous Article in Special Issue
A Quadratic Mean Field Games Model for the Langevin Equation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Maria Francesca Carfora
* and
Isabella Torcicollo
Istituto per le Applicazioni del Calcolo “Mauro Picone” CNR, 80131 Napoli, Italy
*
Author to whom correspondence should be addressed.
Submission received: 31 March 2021 / Revised: 23 April 2021 / Accepted: 27 April 2021 / Published: 30 April 2021
(This article belongs to the Special Issue Differential Models, Numerical Simulations and Applications)

Abstract

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

1. 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 integer-order 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 fractional-order 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 fractional-in-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 spatio-temporal 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-in-time 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.

2. 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
I θ f ( t ) = 1 Γ ( θ ) t 0 t ( t τ ) θ 1 f ( τ ) d τ
and the Caputo fractional derivative of order θ ( m 1 , m ) of a sufficiently differentiable function f ( t ) is defined by
D t θ f ( t ) = 1 Γ ( m θ ) t 0 t ( t τ ) m θ 1 f ( m ) ( τ ) d τ ,
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 .
Lemma 1.
Let x ( t ) be a continuous function on [ t 0 , ) satisfying D θ x ( t ) A x ( t ) + B , x ( t 0 ) = x 0 , with 0 < θ < 1 and A , B R , A > 0 , t 0 0 . Then for x ( t ) it holds true
x ( t ) x 0 B A E θ [ A ( t t 0 ) θ ] + B A ,
where E θ ( · ) is the Mittag–Leffler function defined by
E θ = k = 0 z k Γ ( θ k + 1 ) .
Lemma 2.
Let D θ x ( t ) = ϕ ( t , x ) , with t > t 0 , be the system with the initial condition x ( t 0 ) and 0 < θ 1 , ϕ : [ t 0 , ) × Ω R m , Ω R m . If ϕ ( t , x ) satisfies the Lipschitz condition with respect to x, then there exists a unique solution of the above system on [ t 0 , ) × Ω .
Lemma 3.
Let D t θ x = ϕ ( x ) , 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 .

3. 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 [5]
D t θ n = σ n 1 n k ( 1 + α p ) n p , D t θ p = ( 1 + α p ) n p p ,
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
n ( t 0 ) = n 0 , p ( t 0 ) = p 0 .

3.1. Boundedness

In this subsection, we prove the positivity and boundedness of the solution of the system (3). Let R + 2 = { x ( t ) R 2 : x ( t ) 0 } and x ( t ) = n ( t ) p ( t ) .
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
U ( t ) = n ( t ) + p ( t )
and ξ be a positive constant. Applying the Caputo fractional derivative on both sides in (5) and using (3), it follows that
D t θ U ( t ) + ξ U ( t ) = D t θ n ( t ) + D t θ p ( t ) + ξ n ( t ) + ξ p ( t ) = σ n 1 n k p + ξ n + ξ p = σ n 2 k + ( σ + ξ ) n + ( ξ 1 ) p σ k n k ( σ + ξ ) 2 σ 2 + k ( σ + ξ ) 2 2 σ k ( σ + ξ ) 2 2 σ ,
where ξ < 1 . Using Lemma 1, it follows that
U ( t ) U ( t 0 ) k ( σ + ξ ) 2 2 σ ξ E θ [ ξ t θ ] + k ( σ + ξ ) 2 2 σ ξ U ( t 0 ) E θ [ ξ t θ ] + k ( σ + ξ ) 2 2 σ ξ 1 E θ [ ξ t θ ] .
Then, for t it follows that U ( t ) k ( σ + ξ ) 2 2 σ ξ . Hence, all the solutions to (3) starting from R + 2 are confined in the following domain A R + 2
A = { ( n ( t ) , p ( t ) ) R + 2 : n ( t ) + p ( t ) k ( σ + ξ ) 2 2 σ ξ + ε , ε > 0 }
which is positively invariant. □

3.2. 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) in the region Θ × ( t 0 , T ] , T < , where
Θ = { ( n , p ) R 2 : max ( | n | , | p | ) K } ,
by using the fixed point technique. The existence of K is guaranteed by the boundedness of the solution. The model can be reformulated in the fractional integral form which gives
n ( t ) n ( t 0 ) = 1 Γ ( θ ) t 0 t ( t τ ) θ 1 σ n ( τ ) 1 n ( τ ) k ( 1 + α p ( τ ) ) n ( τ ) p ( τ ) d τ , p ( t ) p ( t 0 ) = 1 Γ ( θ ) t 0 t ( t τ ) θ 1 ( 1 + α p ( τ ) ) n ( τ ) p ( τ ) p ( τ ) d τ .
Denoting by F 1 and F 2 the following kernels
F 1 ( t , n ( t ) , p ( t ) ) = σ n ( t ) 1 n ( t ) k [ 1 + α p ( t ) ] n ( t ) p ( t ) , F 2 ( t , n ( t ) , p ( t ) ) = [ 1 + α p ( t ) ] n ( t ) p ( t ) p ( t ) ,
the following theorem holds.
Theorem 2.
Let M 1 = σ + ( 1 + α K ) K + 2 σ k K and M 2 = | K 1 + 2 α K 2 | . If 0 < M i < 1 , ( i = 1 , 2 ) , 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
F 1 ( t , n 1 , p ) F 1 ( t , n 2 , p = σ 1 n 1 k n 1 ( 1 + α p ) p n 1 σ 1 n 2 k n 2 + ( 1 + α p ) p n 2 = [ σ ( 1 + α p ) p ] n 1 [ σ ( 1 + α p ) p ] n 2 σ k ( n 1 2 n 2 2 ) = [ σ ( 1 + α p ) p ] ( n 1 n 2 ) σ k ( n 1 n 2 ) ( n 1 + n 2 ) σ + ( 1 + α K ) K + 2 σ k K n 1 n 2 M 1 n 1 n 2
and
F 2 ( t , n , p 1 ) F 2 ( t , n , p 2 = ( 1 + α p 1 ) n p 1 p 1 ( 1 + α p 2 ) n p 2 + p 2 = ( n 1 ) ( p 1 p 2 ) + α n ( p 1 2 p 2 2 ) = [ ( n 1 ) + α n ( p 1 + p 2 ) ] ( p 1 p 2 ) | K 1 + 2 α K 2 | p 1 p 2 M 2 p 1 p 2 ,
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
n ( t ) = n ( t 0 ) + 1 Γ ( θ ) t 0 t ( t τ ) θ 1 F 1 ( τ , n , p ) d τ , p ( t ) = p ( t 0 ) + 1 Γ ( θ ) t 0 t ( t τ ) θ 1 F 2 ( τ , n , p ) d τ .
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 )
and
n m ( t ) = n ( t 0 ) + 1 Γ ( θ ) t 0 t ( t τ ) θ 1 F 1 ( τ , n m 1 , p ) d τ , p m ( t ) = p ( t 0 ) + 1 Γ ( θ ) t 0 t ( t τ ) θ 1 F 2 ( τ , n , p m 1 ) d τ .
The successive difference between the terms is defined as
Φ 1 m ( t ) = n m ( t ) n m 1 ( t ) = 1 Γ ( θ ) t 0 t ( t τ ) θ 1 [ F 1 ( τ , n m 1 , p ) F 1 ( τ , n m 2 , p ) ] d τ , Φ 2 m ( t ) = p m ( t ) p m 1 ( t ) = 1 Γ ( θ ) t 0 t ( t τ ) θ 1 [ F 2 ( τ , n , p m 1 ) F 2 ( τ , n , p m 2 ) ] d τ ,
where
n m ( t ) = i = 1 m Φ 1 i ( t ) , p m ( t ) = i = 1 m Φ 2 i ( t ) .
Taking the norm of (14), it follows from the conditions (9)–(10) that
Φ 1 m ( t ) = n m ( t ) n m 1 ( t ) M 1 Γ ( θ ) t 0 t ( t τ ) θ 1 Φ 1 ( m 1 ) ( τ ) d τ , Φ 2 m ( t ) = p m ( t ) p m 1 ( t ) M 2 Γ ( θ ) t 0 t ( t τ ) θ 1 Φ 2 ( m 1 ) ( τ ) d τ .
The following theorem holds.
Theorem 3. 
Assume that the conditions (9)–(10) hold. If
M i T θ Γ ( θ + 1 ) < 1 , i = 1 , 2 ,
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
Φ 1 m ( t )     n 0 ( t ) M 1 T θ Γ ( θ + 1 ) m , Φ 2 m ( t )     p 0 ( t ) M 2 T θ Γ ( θ + 1 ) m .
This shows the existence for the solutions. Moreover, in order to prove that (18) are solutions to (3) and (4), we consider
n ( t ) n ( t 0 ) = n m ( t ) r 1 m ( t ) , p ( t ) p ( t 0 ) = p m ( t ) r 2 m ( t ) ,
where r 1 m , r 2 m 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
r 1 m ( t )     1 Γ ( θ ) t 0 t ( t τ ) θ 1 | F 1 ( τ , n , p ) F 1 ( τ , n m 1 , p ) | d τ 1 Γ ( θ ) t 0 t ( t τ ) θ 1 F 1 ( τ , n , p ) F 1 ( τ , n m 1 , p ) d τ M 1 T θ Γ ( θ + 1 ) n n m 1 .
On using recursive techniques, we get
r 1 m ( t )     T θ Γ ( θ + 1 ) m + 1 M 1 m + 1 .
For m , it follows that r 1 m ( t ) 0 . In a similar way, we conclude that r 2 m ( 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 system n ¯ ( t ) and p ¯ ( t ) . Then
n ( t ) n ¯ ( t )     1 Γ ( θ ) t 0 t ( t τ ) θ 1 F 1 ( τ , n , p ) F 1 ( τ , n ¯ , p ) d τ , p ( t ) p ¯ ( t )     1 Γ ( θ ) t 0 t ( t τ ) θ 1 F 2 ( τ , n , p ) F 2 ( τ , n , p ¯ ) d τ .
In view of the Lipschitz condition, it follows that
n ( t ) n ¯ ( t )     M 1 t θ Γ ( θ + 1 ) n ( t ) n ¯ ( t ) , p ( t ) p ¯ ( t )     M 2 t θ Γ ( θ + 1 ) p ( t ) p ¯ ( t ) ,
and hence
n ( t ) n ¯ ( t ) 1 M 1 t θ Γ ( θ + 1 )     0 , p ( t ) p ¯ ( t ) 1 M 2 t θ Γ ( θ + 1 )     0 .
In view of (17), it follows that n ( t ) n ¯ ( t )   =   0 and p ( t ) p ¯ ( t )   =   0 .  □

4. Stability Analysis

The equilibrium points of system (3) are obtained by considering
D t θ n ( t ) = 0 , D t θ p ( t ) = 0 .
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 * ) 1 ( 1 + α p * ) , p * where p * satisfies α 2 k p 3 + 2 α k p 2 + k ( 1 α σ ) p + σ ( 1 k ) = 0 . As shown in [14], if k > 1 , E * is unique, while if 1 + 1 + 3 σ α α σ < 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
J ( E * ) = σ k ( 1 + α p * ) 1 + 2 α p * 1 + α p * ( 1 + α p * ) p * α p * 1 + α p * .
Now, the characteristic equation of the Jacobian matrix is λ 2 I λ + A = 0 and the roots of the characteristic equation are
λ 1 , 2 = I ± I 2 4 A 2 ,
where I and A are the trace and determinant of the Jacobian matrix J ( E * ) given by
I = k α p * σ k ( 1 + α p * ) , A = p * k ( 1 + α p * ) 2 [ k ( 1 + α p * ) 2 ( 1 + 2 α p * ) α σ ] .
The stability analysis of E * will be divided in the following four cases:
(1)
If I 2 4 A 0 and I < 0 , then λ i ( i = 1 , 2 ) are negative real numbers. Hence, | arg ( λ i ) | = π > θ π 2 . The coexistence equilibrium is asymptotically stable.
(2)
If I 2 4 A 0 and I 0 then at least one of λ 1 , λ 2 will be positive. Therefore, at least one of λ i will be such that | arg λ i | = 0 < θ π 2 and hence the coexistence point is unstable.
(3)
When I 2 4 A < 0 , then the eigenvalues are a pair of complex conjugate λ 1 , 2 = I 2 ± i 4 A I 2 2 .
(a)
If I > 0 then R e ( λ 1 ) = R e ( λ 2 ) > 0 and I m ( λ 1 ) = I m ( λ 2 ) = 4 A I 2 2 > 0 . Then, when 4 A I 2 I > tan θ π 2 , the coexistence equilibrium is asymptotically stable, otherwise it is unstable.
(b)
If I < 0 , then R e ( λ 1 ) = R e ( λ 2 ) < 0 . Then, | arg ( λ i ) | > π 2 and the equilibrium is asymptotically stable.
(4)
When I = 0 and I 2 4 A = 4 A < 0 , then R e ( λ 1 ) = R e ( λ 2 ) = 0 . Hence, | arg ( λ i ) | = π 2 > θ π 2 . 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
  • the corresponding characteristic equation of system (3) has a pair of complex conjugate roots λ 1 , 2 = γ ± i ω with γ > 0 ;
  • m ( θ * ) = θ * π 2 min 1 i 2 | arg λ i | = 0 ;
  • d m ( θ ) d θ | θ = θ * 0 .
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 π arctan [ | I 2 4 A | I ] .
Let us assume that I 2 4 A < 0 and I > 0 ; let the critical value θ * = 2 π arctan [ | I 2 4 A | I ] . Let us define γ = I 2 and ω = | I 2 4 A | 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
m ( θ * ) = θ * π 2 min 1 i 2 | arg λ i | = θ * π 2 arctan ω γ = arctan ω γ arctan ω γ = 0 .
Finally,
d m ( θ ) d θ | θ = θ * = π 2 0 .
Therefore, we can conclude that a Hopf bifurcation of (3) will appear at E * .

5. 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)
    x n = x 0 + h θ i = 1 n b n i ( θ ) F ( t i , x i ) , b n ( θ ) = ( n + 1 ) θ n θ Γ ( θ + 1 ) ;
  • Implicit Trapezoidal Product Integration rule (PI2)
    x n = x 0 + h θ 1 Γ ( θ + 2 ) a n ˜ ( θ ) F ( t 0 , x 0 ) + i = 1 n a n i ( θ ) F ( t i , x i ) ,
    with a 0 ( θ ) = 1 and
    a n ˜ ( θ ) = ( n 1 ) θ + 1 n θ ( n θ 1 ) , a n ( θ ) = ( n 1 ) θ + 1 2 n θ + 1 + ( n + 1 ) θ + 1 ;
  • Predictor–Corrector Product Integration rule (PI12)
    x n P = x 0 + h θ i = 0 n 1 b n i 1 ( θ ) F ( t i , x i ) ,
    x n = x 0 + h θ 1 Γ ( θ + 2 ) a n ˜ ( θ ) F ( t 0 , x 0 ) + i = 1 n 1 a n i ( θ ) F ( t i , x i ) + a 0 ( θ ) F ( t n , x n P ) ;
  • Fractional Backward Differentiation Formula (FLMM2)
    x n = x 0 + h θ i = 0 n ( 1 ) n i θ n i F ( t i , x i ) , where θ n i = Γ ( 1 θ ) Γ ( n + 1 ) Γ ( θ n + 1 ) .
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.

5.1. 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 fractional-order 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]).

5.2. 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 a detail of the first part of the same trajectories as reconstructed by any of the numerical schemes for decreasing time steps h 1 = 2 4 , h 2 = 2 5 , h 3 = 2 6 , h 4 = 2 7 . The less accurate results of PI1 with respect to the other schemes can be clearly seen. Of course, lower values for the fractional-order θ result in an even faster damping of the oscillations in the trajectories, confirming the strong stabilizing effect of the persistent memory in the model. For the same test case, Figure 4 compares the numerical approximations by PI2 of the populations’ trajectories for different values of the order θ (0.95, 0.9, 0.85, 0.8), confirming the strong stabilizing effect of the fractional derivative.

5.3. 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 4 A < 0. The eigenvalues of the Jacobian are a couple of complex conjugate numbers with a positive real part so that for θ > θ * = arctan 2 / π 4 A 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.

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

Author Contributions

Conceptualization, M.F.C. and I.T.; formal analysis, I.T.; numerical simulations, M.F.C.; writing (original draft preparation, review and editing), M.F.C. and I.T.; funding acquisition, M.F.C. and I.T. All authors have read and agreed to the published version of the manuscript.

Funding

This research was partially supported by Regione Campania Projects REMIAM, ADVISE and MEDIA.

Acknowledgments

This paper has been performed under the auspices of the G.N.F.M. and G.N.C.S. of INdAM.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Abrams, P.A. The Evolution of Predator-Prey Interactions: Theory and Evidence. Annu. Rev. Ecol. Syst. 2000, 31, 79–105. [Google Scholar] [CrossRef]
  2. Berryman, A.A. The Origins and Evolution of Predator-Prey Theory. Ecology 1992, 73, 1530–1535. [Google Scholar] [CrossRef] [Green Version]
  3. Briggs, C.J.; Hoopes, M.F. Stabilizing effects in spatial parasitoid–host and predator–prey models: A review. Theor. Popul. Biol. 2004, 65, 299–315. [Google Scholar] [CrossRef] [PubMed]
  4. Courchamp, F.; Berec, L.; Gascoigne, J. Allee Effects in Ecology and Conservation; Oxford University Press: Oxford, UK, 2008. [Google Scholar] [CrossRef]
  5. Alves, M.T.; Hilker, F. Hunting cooperation and Allee effects in predators. J. Theor. Biol. 2017, 419, 13–22. [Google Scholar] [CrossRef]
  6. Tang, B. Dynamics for a fractional-order predator-prey model with group defense. Sci. Rep. 2020, 10, 4906. [Google Scholar] [CrossRef]
  7. Amirian, M.M.; Towers, I.; Jovanoski, Z.; Irwin, A.J. Memory and mutualism in species sustainability: A time-fractional Lotka-Volterra model with harvesting. Heliyon 2020, 6, e04816. [Google Scholar] [CrossRef]
  8. Patnaik, S.; Hollkamp, J.P.; Semperlotti, F. Applications of variable-order fractional operators: A review. Proc. R. Soc. Math. Phys. Eng. Sci. 2020, 476, 20190498. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  9. Sun, H.; Zhang, Y.; Baleanu, D.; Chen, W.; Chen, Y. A new collection of real world applications of fractional calculus in science and engineering. Commun. Nonlinear Sci. Numer. Simul. 2018, 64, 213–231. [Google Scholar] [CrossRef]
  10. Abundo, M.; Ascione, G.; Carfora, M.F.; Pirozzi, E. A fractional PDE for first passage time of time-changed Brownian motion and its numerical solution. Appl. Numer. Math. 2020, 155, 103–118. [Google Scholar] [CrossRef]
  11. Cherstvy, A.G.; Chechkin, A.V.; Metzler, R. Ageing and confinement in non-ergodic heterogeneous diffusion processes. J. Phys. Math. Theor. 2014, 47, 485002. [Google Scholar] [CrossRef]
  12. Wang, W.; Cherstvy, A.G.; Liu, X.; Metzler, R. Anomalous diffusion and nonergodicity for heterogeneous diffusion processes with fractional Gaussian noise. Phys. Rev. E 2020, 102, 012146. [Google Scholar] [CrossRef]
  13. Capone, F.; Carfora, M.F.; De Luca, R.; Torcicollo, I. On the dynamics of an intraguild predator–prey model. Math. Comput. Simul. 2018, 149, 17–31. [Google Scholar] [CrossRef]
  14. Capone, F.; Carfora, M.F.; De Luca, R.; Torcicollo, I. Turing patterns in a reaction-diffusion system modeling hunting cooperation. Math. Comput. Simul. 2019, 165, 172–180. [Google Scholar] [CrossRef]
  15. Capone, F.; Carfora, M.F.; De Luca, R.; Torcicollo, I. Analysis of a model for waterborne diseases with Allee effect on bacteria. Nonlinear Anal. Model. Control 2020, 25, 1035–1058. [Google Scholar] [CrossRef]
  16. Carfora, M.F.; Torcicollo, I. Cross-diffusion-driven instability in a predator-prey system with fear and group defense. Mathematics 2020, 8, 1244. [Google Scholar] [CrossRef]
  17. Carfora, M.F.; Torcicollo, I. Identification of epidemiological models: The case study of Yemen cholera outbreak. Appl. Anal. 2020. [Google Scholar] [CrossRef]
  18. Torcicollo, I. On the nonlinear stability of a continuous duopoly model with constant conjectural variation. Int. J. Non-Linear Mech. 2016, 81, 268–273. [Google Scholar] [CrossRef] [Green Version]
  19. Rionero, S.; Torcicollo, I. On the dynamics of a nonlinear reaction-diffusion duopoly model. Int. J. Non-Linear Mech. 2018, 99, 105–111. [Google Scholar] [CrossRef]
  20. Rionero, S.; Torcicollo, I. Stability of a Continuous Reaction-Diffusion Cournot-Kopel Duopoly Game Model. Acta Appl. Math. 2014, 132, 505–513. [Google Scholar] [CrossRef] [Green Version]
  21. De Angelis, F.; De Angelis, M. On solutions to a FitzHugh-Rinzel type model. Ric. Mat. 2020. [Google Scholar] [CrossRef] [Green Version]
  22. Capone, F.; De Luca, R.; Rionero, S. On the stability of non-autonomous perturbed Lotka-Volterra models. Appl. Math. Comput. 2013, 219, 6868–6881. [Google Scholar] [CrossRef]
  23. Podlubny, I. Fractional Differential Equations. An Introduction to Fractional Derivatives, Fractional Differential Equations, Some Methods of Their Solution and Some of Their Applications; Academic Press: San Diego, CA, USA; London, UK, 1999. [Google Scholar]
  24. Diethelm, K. The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using Differential Operators of Caputo Type; Springer: Berlin/Heidelberg, Germany, 2010. [Google Scholar] [CrossRef] [Green Version]
  25. Kilbas, A.A.; Srivastava, H.M.; Trujillo, J.J. Theory and Applications of Fractional Differential Equations, Volume 204 (North-Holland Mathematics Studies); Elsevier Science Inc.: Amsterdam, The Netherlands, 2006. [Google Scholar]
  26. Li, X.; Wu, R. Hopf bifurcation analysis of a new commensurate fractional-order hyperchaotic system. Nonlinear Dyn. 2014, 78, 279–288. [Google Scholar] [CrossRef]
  27. Diethelm, K.; Ford, N.J.; Freed, A.D.; Luchko, Y. Algorithms for the fractional calculus: A selection of numerical methods. Comput. Methods Appl. Mech. Eng. 2005, 194, 743–773. [Google Scholar] [CrossRef] [Green Version]
  28. Garrappa, R. Trapezoidal methods for fractional differential equations: Theoretical and computational aspects. Math. Comput. Simul. 2015, 110, 96–112. [Google Scholar] [CrossRef] [Green Version]
  29. Cai, M.; Li, C. Numerical Approaches to Fractional Integrals and Derivatives: A Review. Mathematics 2020, 8, 43. [Google Scholar] [CrossRef] [Green Version]
  30. Garrappa, R. Numerical Solution of Fractional Differential Equations: A Survey and a Software Tutorial. Mathematics 2018, 6, 16. [Google Scholar] [CrossRef] [Green Version]
  31. Garrappa, R. FLMM2. MATLAB Central File Exchange. 2021. Available online: https://www.mathworks.com/matlabcentral/fileexchange/47081-flmm2 (accessed on 26 March 2021).
  32. Yin, H.; Wen, X. Pattern Formation through Temporal Fractional Derivatives. Sci. Rep. 2018, 8, 5070. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Figure 1. Trajectories for the numerical approximations PI1 and PI2 of n ( t ) solution of system (3) for θ = 1 in case of a stable equilibrium (left panel) and in the presence of Hopf instability (right panel). The lowest order approximation PI1 clearly shows damped oscillations, while the trajectories of the higher order approximations are in excellent agreement with the ones obtained by classical ODE solvers.
Figure 1. Trajectories for the numerical approximations PI1 and PI2 of n ( t ) solution of system (3) for θ = 1 in case of a stable equilibrium (left panel) and in the presence of Hopf instability (right panel). The lowest order approximation PI1 clearly shows damped oscillations, while the trajectories of the higher order approximations are in excellent agreement with the ones obtained by classical ODE solvers.
Axioms 10 00078 g001
Figure 2. Trajectories for all the numerical approximations of the solutions of system (3) in case θ = 0.9. While the corresponding ODE systems show instability (see the right panel of 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.
Figure 2. Trajectories for all the numerical approximations of the solutions of system (3) in case θ = 0.9. While the corresponding ODE systems show instability (see the right panel of 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.
Axioms 10 00078 g002
Figure 3. Details of the initial part of the trajectories shown in the left panel of Figure 2 as reconstructed by all the numerical methods for decreasing timesteps h 1 , h 2 , h 3 , h 4 . The lowest order approximation PI1 clearly shows more damped oscillations, and only for the smallest timestep it agrees with the other methods, whose results are practically identical.
Figure 3. Details of the initial part of the trajectories shown in the left panel of Figure 2 as reconstructed by all the numerical methods for decreasing timesteps h 1 , h 2 , h 3 , h 4 . The lowest order approximation PI1 clearly shows more damped oscillations, and only for the smallest timestep it agrees with the other methods, whose results are practically identical.
Axioms 10 00078 g003
Figure 4. Trajectories for the numerical approximations PI2 of the solutions of system (3) for different θ values (0.95, 0.9, 0.85, 0.8) for both populations ( n ( t ) shown in the left panel, p ( t ) in the right one). The stabilizing effect of the fractional operators can be clearly seen.
Figure 4. Trajectories for the numerical approximations PI2 of the solutions of system (3) for different θ values (0.95, 0.9, 0.85, 0.8) for both populations ( n ( t ) shown in the left panel, p ( t ) in the right one). The stabilizing effect of the fractional operators can be clearly seen.
Axioms 10 00078 g004
Figure 5. Numerical approximation PI12 of the solutions of system (3) confirming the Hopf instability predicted by the theory in case θ = 0.92 with σ = 3, α = 3.5, k = 5. In the left panel, the trajectories for both variables are reported as functions of time; the right panel shows the limit cycle in the phase plan.
Figure 5. Numerical approximation PI12 of the solutions of system (3) confirming the Hopf instability predicted by the theory in case θ = 0.92 with σ = 3, α = 3.5, k = 5. In the left panel, the trajectories for both variables are reported as functions of time; the right panel shows the limit cycle in the phase plan.
Axioms 10 00078 g005
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Carfora, M.F.; Torcicollo, I. A Fractional-in-Time Prey–Predator Model with Hunting Cooperation: Qualitative Analysis, Stability and Numerical Approximations. Axioms 2021, 10, 78. https://0-doi-org.brum.beds.ac.uk/10.3390/axioms10020078

AMA Style

Carfora MF, Torcicollo I. A Fractional-in-Time Prey–Predator Model with Hunting Cooperation: Qualitative Analysis, Stability and Numerical Approximations. Axioms. 2021; 10(2):78. https://0-doi-org.brum.beds.ac.uk/10.3390/axioms10020078

Chicago/Turabian Style

Carfora, Maria Francesca, and Isabella Torcicollo. 2021. "A Fractional-in-Time Prey–Predator Model with Hunting Cooperation: Qualitative Analysis, Stability and Numerical Approximations" Axioms 10, no. 2: 78. https://0-doi-org.brum.beds.ac.uk/10.3390/axioms10020078

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop