Dynamical Behaviors of a Fractional-Order Three Dimensional Prey-Predator Model

In this paper, the dynamical behavior of a three-dimensional fractional-order prey-predator model is investigated with Holling type III functional response and constant rate harvesting. It is assumed that the middle predator species consumes only the prey species, and the top predator species consumes only the middle predator species. We also prove the boundedness, the non-negativity, the uniqueness, and the existence of the solutions of the proposed model. Then, all possible equilibria are determined, and the dynamical behaviors of the proposed model around the equilibrium points are investigated. Finally, numerical simulations results are presented to confirm the theoretical results and to give a better understanding of the dynamics of our proposed model.


Introduction
The most influential theme in ecology and mathematical modeling is the dynamic of the relationship among species. Many authors extended or modified the work of the Lotka and Voltera model [1,2], and they have investigated thse topics widely by using ordinary differential equations or deference equations, see [3][4][5][6][7][8][9], and references therein. Nowadays, authors formulate their mathematical models by fractional order differential equation due to their ability to give the precise description for various linear and nonlinear phenomena [10][11][12][13][14][15]. The increasing of mathematical models that based on fractional order differential equation has recently obtained popularity in the studying the behavior of biological models. Fractional-order differential equation has been successfully used and applied to model many areas of science, engineering, and phenomena that cannot be formulated by other types of equations [10,16]. In [12,[16][17][18][19][20][21], authors have investigated the effects of the fractional order differential equation on a prey predator model as well as they also discussed the stability analysis of equilibrium points of fractional order model with and without harvest-ing, as well as the existence, uniqueness, and boundedness of the solutions that are proved.
There are several different types of definitions of fractional-order differential equation in the literatures [16,[22][23][24][25], for e.g., Caputo, Riemann-Litouville, Fabrizio, Marechand,Grunwald-Letnikov, wayl, and Riesz fractional-order derivatives. Throughout this work, we used the Caputo fractional-order derivatives because it is not necessary to define the initial conditions of fractional-order and its fractional-order derivative of constant function is zero, as well as the similarity of the initial conditions of fractional order differential and the integer order ones [8]. This work is organized as follows: in Section 2, some useful definitions and concept that concern to the fractional order are presented. In Section 3, a three-dimensional fractional-order prey-predator model is considered. The uniqueness and boundedness as well as the nonnegativity of its solution are proved. In Section 4, all the equilibrium points are determined, and the conditions to achieve its local stability are set. In Section 5, numerical simulations are given to confirm the theoretical results. Finally, conclusions are given in Section 6.

Preliminaries
Definition 1. [24,26]: Caputo's definition of fractional derivatives is given as follows: D α t f ðtÞ = I ½α−α f ½α ðtÞ, α > 0, where ½α is the least integer which is not less than α, and I θ is the Riemann-Liouville integral operator of order θ which is given by Throughout this work, we need some useful results for the fractional-order derivative. We start with the following results which are proved in [13,16,27]. Lemma 1.

Model Formulation
Frist, the following three dimension prey-predator model with functional response of Holling type III is considered as follows: where X, Y, and Z are the density of prey species, middle species, and the top predator species. R 0 , K 0, , c 1 , c 2 , and b 1 are the intrinsic growth rates, the carrying capacity of the prey population, the death rates for middle predator population and top predator population, and capture rate of prey and middle predator, respectively. The parameter a 1 is the conservation rate of prey X, and the parameters a 2 and b 2 are the conservation rate of middle predator Y to the top predator Z, respectively. The constants B 1 and B 2 are the predation rates and half saturation constants, respectively. The parameter h 2 denotes to the constant harvesting rate.
We use the following transformation to reduce the dimension of model (8) and t = R 0 T, and then we obtain the following system: 2 Abstract and Applied Analysis Next, we introduce the fractional-order derivative α in model (2) with the help of fractional order Caputo type derivative. Then, the system (9) becomes as follows: Remark 3. Because of the biological significance, we are interesting solutions that are nonnegative and bounded only, so that we will prove the following theorem.

The Equilibrium Points and Local Stability Analysis
In this section, the existence and the local stability of equilibrium points of the considered system (10) are studied and investigated. Before that, we need the following theorem.
Theorem 6 (see [29]). Consider the fractional-order differential system d α x/dt α = f ðxÞ; xð0Þ = x 0 , with α ∈ ð0, 1Þ and x ∈ ℝ n . The equilibrium points of the above system are solutions to the equation f ðxÞ = 0. An equilibrium point is locally asymptotically stable if all the eigenvalues λ j of the Jacobian matrix J = ∂f /∂x evaluated at the equilibrium satisfy jarg ðλ j Þj > απ/2.
Proof (see [29]). To obtain the equilibrium points of the considered system (10), we solve the following simultaneous equations: Thus, the considered system (10) has the following possible equilibrium points: The Jacobian matrix of system (10) at ðx, y, zÞ is given by where k = ð1 + ax 2 Þ and k 1 = ð1 + dy 2 Þ.
The next theorems give the local stability of E 0 , E 1 , and E 2 , respectively. Theorem 7. The trivial equilibrium point E 0 of the system (10) is always unstable point.
In order to discuss the local stability of the interior point, we need the following definition and lemma.

Theorem 12.
The interior equilibrium point E 3 is locally asymptotically stable if one of the following conditions holds:  Therefore, the characteristic equation of JðE 3 Þ is given by Thus, the discriminant DðpÞ of the cubic polynomial pðλÞ is where Then, all results can be obtained by Lemma 6.

Numerical Simulation
All previous theoretical results in Section 4 are confirmed by giving numerical simulations. At different set of values of parameters, the local stability of E 1 , E 2 , and E 3 of the considered system (10) is investigated numerically. First, for the point E 1 , the following set of values of the parameters is chosen : a 1 = 0:5 ; c 1 = 0:1 ; b 1 = 0:5 ; d 1 = 0:2 ; m 1 = 0:15 ; m 2 = 0:1 ; h = 0:8,α = 0:98, and (0.1,0.5,0.9) is the initial value. Therefore, the condition in Theorem 4 is satisfied. Figure 1 shows that the point E 1 is locally stable.

Conclusions
Fractional-order differential equations have recently been successfully used and applied to model many areas of science, especially, phenomena that cannot be formulated by other types of equations, so that this work concerns a study of a three-dimensional fractional-order prey-predator model with Holling type III functional response and constant rate harvesting. It is shown that the model possesses the existence, nonnegativity, boundedness, and uniqueness of the solutions, as desired in any population dynamics. The local stability of all possible equilibrium points of the considered system is investigated analytically then confirms by numerical verifications.

Data Availability
All simulations in the article are used by the authors to confirm the theoretical results. We have also cited references that were used, and we gave the DOI in the references.

Conflicts of Interest
The authors declare that they have no conflicts of interest.