STABILITY ANALYSIS OF A FRACTIONAL PREDATOR-PREY SYSTEM WITH TWO DELAYS AND INCOMMENSURATE ORDERS

In this paper, we consider a fractional predator-prey system with two delays and incommensurate orders. Firstly, the local stability of positive equilibrium of the system without delay is discussed. Secondly, we calculate the critical value of Hopf bifurcation by taking one delay as bifurcation parameter. Then, as two nonidentical delays change simultaneously, the stability switching curves, the directions of crossing and the existence of Hopf bifurcation are obtained. Finally, numerical simulations are presented to verify the given theoretical results.


Introduction
Fractional calculus is a kind of generalization of integer-order calculus. Due to the effects of memory and hereditary properties, fractional calculus can more accurately describe the complex and rich dynamic behavior of the system rather than integerorder calculus [12,15]. Because of the complexity of the calculation, fractional calculus has not been widely used for a long time. It was not until fractal theory was introduced in [18] that fractional calculus has progressively become one of the research hotspots. In recent years, fractional calculus has been successfully introduced into various fields, such as physics, chemistry, electricity, biology, economics, etc [3,5,7,8,21,22].
It takes a certain amount of time to complete biological evolution and physical process, therefore delay is widely found in nature. The appearance of time delay means that the development of the system is not only dependent on the current state, but also related to the state in the previous period. On the other hand, time delay can cause Hopf bifurcation. When the parameter changes slightly near a critical value, the stability of the equilibrium of system changes, and there is the phenomenon of periodic solution near equilibrium. For fractional-order delay differential systems, Hopf bifurcation is also a common bifurcation phenomenon.
In [2,17,25,30], the authors chose the delay as bifurcation parameter, and studied Hopf bifurcation of fractional-order delay differential systems.
The predator-prey model is an important model in population dynamics models [31]. Authors in [1] analyzed the dynamic behavior of the prey-predator system with time delays, and it was shown that the stability of the system can be changed with the change of harvest. In [26], authors studied a delayed predator-prey system with Beddington-DeAngelis functional response and got that the phenomenon of Hopf bifurcation is the main factor that switching the stability of the system to unstability with respect to delays. In [28], Wang and Tang took the Holling-III functional response, the complex diversity of biological environment and efforts into account, and discussed the Hopf bifurcation of a time-delayed predator-prey system which is as follows: where x(t) and y(t) are the densities of prey and predator, respectively. The growth rate and the environmental capacity of prey are represented by r and K, respectively. α denotes the attack coefficient; c (0 < c < 1) represents a dimensionless parameter that measures habitat complexity; d is the death rate of the predator; q is the coefficient of catchability; E 0 denotes the harvesting effort; and α(1−c)x 2 y is Holling-III functional response. All parameters of system (1.1) are positive. For integer-order delay differential equations, a lot of researchers have paid too much attention to predator-prey systems with single delay, such as [4,14]. Even though the integer-order predator-prey models with two delays have been discussed in [9,29], the approaches are to make the two delays equal or fix one delay and choose the other one as bifurcation parameter. There is little work to consider differential equations with two delays varying simultaneously. In [10,16], integerorder systems with two time delays varying simultaneously were discussed. By analyzing the characteristic equations of the following two forms D(s; τ 1 , τ 2 ) = U 0 (s) + U 1 (s)e −sτ1 + U 2 (s)e −sτ2 and D(s; τ 1 , τ 2 ) = U 0 (s) + U 1 (s)e −sτ1 + U 2 (s)e −sτ2 + U 3 (s)e −s (τ1+τ2) respectively in [10,16], calculating the explicit expression of the stability switching curves, and giving the judgment method of the directions of change in stability, then the stability of the system was obtained. The method of [10,16] has been applied to study the stability of the integer-order systems with two time delays varying simultaneously, such as [11,20,23].
For fractional-order delay differential equations, most of systems concerned by scholars are systems with one or two delays, but the ways to study the systems with two delays are either to make two delays equal or only fix one delay and choose the other delay as bifurcation parameter, such as [13,25,27,32]. As far as we are aware, few authors have studied the stability of fractional-order systems when two delays change simultaneously. Therefore, it is significant to extend the method of [10,16] to fractional-order systems with two delays.
Motived by the method of stability switching curves of [10,16] and based on the system of [28], we consider a fractional-order predator-prey system with two delays and incommensurate orders: where γ 1 , γ 2 ∈ (0, 1], γ 1 ̸ = γ 2 . The other parameters have the same biological significance as system (1.1), and initial conditions are x(t) > 0 and y(t) > 0, t ∈ [− max{τ i }, 0] (i = 1, 2). The highlights of this paper are generalized as follows: (i) We generalize integer-order delayed predator-prey system to a new fractional predator-prey system with two delays and incommensurate orders.
(ii) The stability and the Hopf bifurcation of a fractional predator-prey system with two delays and incommensurate orders are obtained. It is shown that the delay and the order play an important role in the stability and the existence of Hopf bifurcation of the corresponding fractional-order system.
(iii) To the best of our knowledge, there are not many results on the stability of fractional-order system and the existence of Hopf bifurcation with two delays varying simultaneously. Discriminating from the general ways that making two delays equal or fixing one delay and choosing another delay as bifurcation parameter, this paper is the case that the stability of system and the existence of Hopf bifurcation are obtained by taking two delays as bifurcation parameters, and the stable region is a two-dimensional region about τ 1 and τ 2 . (iv) The method of determining the stability of the system by calculating the stability switching curves is first applied to the fractional differential equations with two delays. It is meaningful to study fractional-order differential equations with two delays.
The paper is organized as folllows. In Section 2, some basic knowledge and necessary lemmas on fractional calculus are presented. In Section 3, using the method in [10,16], the stability and the existence of Hopf bifurcation of system (1.2) are obtained by calculating the stability switching curves and the directions of change in stability, taking the two delays as bifurcation parameters and considering the simultaneous change of the two delays. In Section 4, we perform numerical simulations to confirm the theoretical results. This paper ends with a conclusion.

Preliminaries
In fractional derivatives, the Grunwald-Letnikov(G-L) definition, the Riemann-Liouville(R-L) definition and the Caputo definition are commonly used. Then the Caputo definition can not only simplify the Laplace transform properly, but also allow that the initial conditions of the corresponding fractional-order equation can be expressed in integer order, which is more suitable for practical mathematical problems. In this paper, all the fractional derivatives are the Caputo definition. Next, we will introduce the Caputo definition and two fundamental lemmas for the following theoretical analysis. Definition 2.1 ( [24]). The Caputo fractional-order derivative is defined as follows: When 0 < γ < 1, we have where γ i ∈ (0, 1] (i = 1, 2, · · · , n), and the initial conditions where ∆(s) is the characteristic matrix of system (2.2). The zero solution of system (2.2) is locally asymptotically stable if all the roots of det(∆(s)) = 0 possess negative real parts.

Main results
In this section, we discuss the stability of system (1.2) and the existence of Hopf bifurcation. We first analyze the existence of positive equilibrium of system (1.2).
Next, we discuss the local stability of the positive equilibrium of the system (1.2) without delay by Routh-Hurwitz criterion. Then, we address the existence of Hopf bifurcation at the positive equilibrium of the system (1.2) with one delay by taking τ 1 and τ 2 as bifurcation parameter, respectively. Finally, applying the method of [10,16], we study the change in stability of the positive equilibrium and the existence of Hopf bifurcation of the system (1.2) when two delays change simultaneously.

Existence of positive equilibrium of system (1.2)
The equilibria J 0 , J K , J * of system (1.2) can be obtained by the following equations: When the following assumptions (H 1 ) and (H 2 ) hold: J * is a unique positive equilibrium of system (1.2).
In conclusion, all common eigenvalues of Eq.

Stability analysis of positive equilibrium J * of system (1.2) with one delay
In this subsection, we analyze the stability of positive equilibrium J * of system (1.2) with one delay by taking τ 1 and τ 2 as bifurcation parameter, respectively.
To further present our main results, the following assumption is needed: Proof. Taking the derivative of Eq.(3.7) with respect to τ 1 , then it is deduced that Thus, we can get that It acquires from Eq.(3.12) that (3.13) Hence, (H 4 ) implies that transversality condition holds. This ends the proof of Lemma 3.1.
According to Lemma 2.1, Lemma 2.2 and Lemma 3.1, the following theorem can be concluded. (ii) Hopf bifurcation will happen around the positive equilibrium J * of system (1.2) for τ 1 = τ 10 .
In order to present our main results, the additional hypothesis is useful and necessary: where M 1 and M 2 can be obtained by Eq. (3.21). Proof. Taking the derivative of Eq.(3.14) with respect to τ 2 , it is calculated that Clearly, we can get that  (3.21) Therefore, (H 5 ) implies that transversality condition holds. This fulfills the proof of Lemma 3.2.
According to Lemma 2.1, Lemma 2.2 and Lemma 3.2, we can get the following theorem. To analyze the stability of system (1.2) when two delays change at the same time, stability switching curves, directions of crossing and Hopf bifurcation are discussed in this subsection. Let R, R + , N 0 , Z and C be the sets of real numbers, nonnegative real numbers, nonnegative integer numbers, integer numbers and complex numbers, respectively. Similary, R 2 and R 2 + denotes the sets of 2-dimensional vectors with components in R and R + .
Based on [16], in order to ensure that the equation of the form of Eq.
(H 9 ): The polynomials U i (s) (i = 0, 1, 2, 3) satisfy the condition, which is To make sure that the Eq.
That means that (H 6 ) is true.

Stability switching curves
Applying the method of [10,16], the feasible region of the system (1.2) is found in this subsubsection. Then all pairs of points (τ 1 , τ 2 ) in the feasible region such that the characteristic equation (3.3) at least has one pair of pure imaginary roots, which constitute the stability switching curves T . We need to find a series of points (τ 1 , τ 2 ) ∈ R 2 + , such that the charisteristic equation (3.3) has a root s = iω = ω(cos π 2 + i sin π 2 ) (ω > 0). It's obvious by (H 7 ) that s ̸ = 0. Therefore, the characteristic equation (3.3) can be written as Since |e −iωτ2 | = 1, we can get (3.23) Square both sides of Eq.(3.23), and it could be equal to By calculation, we can obtain that where There is a continnous function such that which means that Therefore, Eq.(3.24) can be equivalently written as Since | cos(θ 1 (ω) + ωτ 1 )| ≤ 1, then for any τ 1 ∈ R + , there's a necessary condition that satisfies the above equation: Based on (H 7 ), we can get that F (0) ̸ = 0. When ω → +∞, then F (ω) → +∞ by (H 6 ) and (H 9 ). Therefore, F (ω) has a finite number of roots on R + . The range of ω satisfying Eq.(3.26) is denoted by Ω, which is Another way to calculate τ 2 is to analyze τ 2 in the same way that we analyzed for τ 1 , which gives where Similarly to Eq.(3.26), we have By comparing Eq.(3.26) and Eq.(3.32), it is not difficult to find that they are equivalent. The range of ω that satisfies Eq.(3.26) is the same as the range of ω that satisfies Eq.(3.32), and they are denoted as Ω. Ω is also known as the feasible region and it is consistent with the range of all frequencies corresponding to all points on the stability switching curves.
From the above derivation, we can get that the stability switching curves are (3.33)

Directions of crossing
With the previous discussion, in this subsubsection, we will focus on the directions of change in the stability of system (1.2).
In the neighborhood of (τ * 1 , τ * 2 ), η(τ * 1 , τ * 2 ) = 0 and ω(τ * 1 , τ * 2 ) = ω * are satisfied. We call the increasing direction of ω ∈ Ω as the positive direction of the stability switching curves T , and as moving along the positive direction of the curves T , the left-hand (right-hand) side is called as the region on the left (right).
Due to the tangent vector of T along the positive direction is m = ( . Also, as η increases from negative to positive through 0, the direction of a pair of pure imaginary roots of characteristic equation (3.3) across the imaginary axis to the right on the complex pane is determined by the sign of the inner product of m and n, which is (3.34) If δ(ω * ) > 0 (δ(ω * ) < 0), then the region on the right (left) has characteristic roots with positive real parts when moving along the positive direction of stability switching curves T .

Hopf bifurcation
Taking the derivstive of τ 1 of Eq.(3.3), we can get that 0, 1, 2, 3). Eq.(3.22) can be written equivalently as follows: By Eq.(3.40) and Eq.(3.41), it can get that Let If the following condition holds: then the transversal condition is true.
Clearly, we can obtain the following theorem.

Conclusion
In this paper,we have investigated a fractional predator-prey system with two delays and incommensurate orders. Taking time delay as bifurcation parameter, the  (i) When τ 1 = τ 2 = 0, the local stability of the positive equilibrium of the system (1.2) is analyzed by Routh-Hurwitz criterion.
(iii) When τ 1 > 0, τ 2 > 0, and τ 1 ̸ = τ 2 , applying the method of [10,16], we can calculate the stability switching curves and the directions of crossing, and obtain the change in the stability of the positive equilibrium of the system (1.2) and the existence of Hopf bifurcation as two delays change simultaneously.
From the discussion of the above cases, we can get the following results: (i) Delay has an important effect on the stability of system. When the delay crosses a critical value, Hopf bifurcation occurs in the system (1.2).
(ii) Contrasted the system (1.2) with integer order and fractional order, the latter has widely stability region. The order has an effect on the stability of the system.