Dynamical analysis of a new fractional-order predator–prey system with Holling type-III functional

In this paper, we consider a new fractional-order predator–prey model with Holling type-III functional response and stage structure. Based on the Lyapunov stability theory and by constructing a suitable Lyapunov functional, we obtain some sufficient conditions for the existence and uniqueness of positive solutions and the asymptotic stability of the positive equilibrium to the system. Finally, we give some numerical examples to illustrate the feasibility of our results.


Introduction
It is well known that the existence and uniqueness of positive solutions for predator-prey models with Holling type-III functional response and stage structure have been widely investigated by many researchers [1][2][3][4][5][6][7][8][9]. Stage structure is an important natural phenomenon and represents, for example, the division of a population into immature and mature individuals. In [4], the authors studied the global properties of a predator-prey model with nonlinear functional response and stage structure for the predator, and the conditions of the existence and the global stability of the positive steady state were established. In [9], the authors discussed the existence and local stability of equilibrium points, and in order to protect the stability of this kind of systems, they proposed a hybrid control method such that the Hopf bifurcation can be controlled.
Fractional order calculus is the extension of integer order calculus to arbitrary real number order; both appeared almost at the same time. Most of the present works were focused on fractional differential equations, see [10][11][12][13][14][15][16] and the references cited therein. The basic theory of fractional calculus can be found in the monographs of Miller and Ross [10]. With the improvement of fractional calculus theory, the fractional differential equations are widely used in various fields, such as physics [17][18][19], economics [20][21][22], medicine [23,24], and biology [25,26], etc. Because of the memory and heritage properties of fractional calculus, it is more suitable for describing the population dynamics system than integer calculus. Therefore, it is more in line with the laws of nature and has practical sig-nificance to study predator-prey models by using fractional calculus. There are many definitions of fractional derivative and integral, and those of Grünwald-Letnikov, Riemann-Liouville, and Caputo are commonly used. To use the Riemann-Liouville definition one must specify the value of a fractional derivative of the unknown solution at the initial value t = 0, but this can simplify the calculation of the fractional derivative. For complex multiscale analysis, the Riemann-Liouville definition can simplify the calculation process. The shortcomings of the Riemann-Liouville fractional derivative can be compensated by modifying the definition. Therefore, in this paper we will adopt the modified Riemann-Liouville-type definition given in [12].
Recently, a fractional predator-prey model to describe the ecosystem was considered in [27], where it performed well on a practical problem. The author considered the fractional predator-prey model with Holling type-II functional response. Later, the research on the dynamics of fractional-order predator-prey model with Holling type-II or Holling type-III functional response has become a hot topic [28][29][30][31].
To the best of our knowledge, up to now, few results are available for fractional-order predator-prey systems with Holling type-III functional response and stage structure. Therefore, it is a challenging and important problem in theory and applications.
Motivated by the above discussion, in this paper, we consider the following fractionalorder predator-prey system with Holling type-III functional response and stage structure: where x 1 (t) and x 2 (t) represent the densities of the immature and mature prey at time t, respectively; y(t) represents the density of the predator at time t; the parameters a, r 1 , b, r 2 , b 1 , a 1 , m, a 2 , and r are positive constants. Through calculation, it is easy to get that system (1.1) always has a trivial equilibrium E 0 (0, 0, 0) and if the condition ab > r 2 (b + r 1 ) holds, then (1.1) has a predator-extinction equilibrium E 1 (x 0 1 , x 0 2 , 0), where If the conditions ab > r 2 (b + r 1 ) and (a 2mr)(abr 2 (b + r 1 )) 2 > rb 2 1 (b + r 1 ) 2 hold, then system (1.1) has a unique coexistence equilibrium E * (x * 1 , x * 2 , y * ), where Our main purpose of this paper is to study the existence and asymptotical stability of equilibria for system (1.1). Especially, we will focus on the considerations that system (1.1) has a unique coexistence equilibrium.
The main contributions of this paper are listed as follows. Firstly, fractional-order predator-prey system with Holling type-III functional response and stage structure defined by modified fractional derivative is proposed. Secondly, we consider the existence and stability of the equilibrium point of the system, which has important practical significance for the sustainable development of the ecosystem. Finally, examples and numerical simulations are given to verify the effectiveness of the conclusion.
The organization of the rest of this paper is as follows. In Sect. 2, we introduce some good properties on the fractional derivative to study fractional systems. We show that there is a unique nonnegative solution of (1.1) in Sect. 3. In Sect. 4, we use the Lyapunov stability theory [19] of the fractional system to prove that the positive equilibrium is asymptotically stable.

Preliminaries
In this section, we introduce notations, some basic definitions, and preliminaries that will be used in this paper.
, be a continuous (but not necessarily differentiable) function, and let h > 0 be a constant discretization stepsize. Define the forward operator FW (h), i.e., (the symbol := means that the left side is defined by the right one) and its derivative of fractional order is defined by the expression

Definition 2.2 ([12]
, Riemann-Liouville definition revisited) Refer to the function of Definition 2.1. Then its fractional derivative of order α is defined by the expression For positive α, one will set Now, we give some properties of the modified Riemann-Liouville derivative [14] which are used further in this paper:

Lemma 2.1 ([12])
The following equalities hold: The solution of the equation is defined by the following result: In [19], the authors consider the equation 3) has an asymptotically stable null solution;

Existence and uniqueness of the nonnegative solution
Throughout this paper, let R 3 + denote the positive cone of R 3 , namely, R 3

Theorem 3.1 For any initial value
Proof Since the right-hand side coefficients of the system (1.1) are locally Lipschitz continuous, referring to the proof of Theorem 4 of [27], we omit the proof of existence of a local solution in this paper. Hence, for any given initial value (x 1 (0), x 2 (0), y(0)) ∈ R 3 + , there is a unique maximal local solution (x 1 (t), x 2 (t), y(t)) on t ∈ [0, τ e ), where τ e is the explosion time. To show this solution is global, we need to show that τ e = ∞. Let k 0 ≥ 1 be sufficiently large so that x 1 (0), x 2 (0) and y(0) all lie within the interval [ 1 k 0 , k 0 ]. For each integer k > k 0 , define the stopping time by where, throughout this paper, we set inf ∅ = ∞ (∅ denotes the empty set). Obviously, τ k is increasing as k → ∞. Set τ ∞ = lim k→∞ τ k , where τ ∞ < τ e . If we can show that τ ∞ = ∞, then τ e = ∞ and (x 1 (t), x 2 (t), y(t)) ∈ R 3 + for all t ≥ 0. Namely, to complete the proof, it is sufficient to show that τ ∞ = ∞. If this statement is false, then there is a constant T > 0 such that Hence, there is an integer k 1 > k 0 such that the nonnegativity of this function can be seen from u -1 -log u ≥ 0, ∀u.
Using the chain rule, we get Now, we pay attention to the term 1 2 ( 1 √ x 1 -1 x 1 )ax 2 . If x 1 < 1, then 1 2 ( 1 √ x 1 -1 x 1 )ax 2 < 0, hence, this term can be omitted from the right-hand side of the inequality. If x 1 > 1, then 1 2 So we get the following inequality: where M is a positive constant. Then d α V ≤ M dt α . Therefore, on the one hand, Using equality (2.2), we get On the other hand, using the differential relation Note that if at least one of x 1 (τ k ), x 2 (τ k ) and y(τ k ) equals either k or 1 k , we know that V (x 1 (τ k ∧ T), x 2 (τ k ∧ T), y(τ k ∧ T)) is no less than min{( Letting k → ∞ leads to the contradiction This contradiction shows that τ ∞ = ∞, which completes the proof.

The stability of the solution
In this section, by constructing an appropriate Lyapunov functional, we shall study the stability of the solution of system (1.1).

Define a Lyapunov function
where c i (i = 1, 2, 3) are positive numbers to be determined; the nonnegativity of this function can be seen from the property that lim t→+∞ x 1 (t) = x * 1 , lim t→+∞ x 2 (t) = x * 2 , lim t→+∞ y(t) = y * .
Then using the chain rule, we obtain .
. Then we have Then, from Lemma 2.3 and the condition of Theorem 4.1, the above implies d α V dt α < 0. Thus system (1.1) is asymptotically stable. The proof is complete.
Hence, system (5.1) has a unique coexistence equilibrium which is asymptotically stable. The results are verified by the numerical simulations in Figs. 1-3.
It is easy to check that the condition of Theorem 4.1 is satisfied. So, system (5.2) has a unique coexistence equilibrium