#### Abstract

We propose a fractional-order model of the interaction within two-prey and one-predator system. We prove the existence and the uniqueness of the solutions of this model. We investigate in detail the local asymptotic stability of the equilibrium solutions of this model. Also, we illustrate the analytical results by some numerical simulations. Finally, we give an example of an equilibrium solution that is centre for the integer order system, while it is locally asymptotically stable for its fractional-order counterpart.

#### 1. Introduction

The branch of mathematics that concerned biology is called mathematical biology. Mathematical biology tries to model, study, analyze, and interpret biological phenomenon such as the interactions, coexistence, and evolution of different species [1–4]. These interactions may be among the individuals of the same species and among the individuals of different species or interactions against the environment, disease, and food supply. The most important one is the interaction among the individuals of different species which can be predatory, competitive, or mutualistic. A large number of simple mathematical models have been suggested to understand the predator-prey interaction. The work done independently by Lotka [5] and Volterra [6], which is known as Lotka-Volterra model, was the first stone in this field. Later, the model was extended to include density-dependent prey growth and a functional response [7]. After that, a huge number of variants of this model were suggested. Initially, the authors proposed a system of two-dimensional coupled differential equations [8, 9]. Then, the scenario changed to discrete mathematical models which included a lot of complex dynamical behaviours [10, 11]. In the last two decades, the fractional-order differential equations appeared and began to study the predator-prey models in the fractional-order form [12–19].

Fractional calculus generalizes the concept of ordinary differentiation and integration to noninteger order. Fractional calculus is a fertile field for researchers to study very important real phenomena in many fields like physics, engineering, biology, and so forth [20–30]. The fractional differential equations are naturally related to systems with memory. Also, they are closely related to fractals which are numerous in biological system. The definition of fractional derivative involves an integration which is nonlocal operator. This means that the fractional derivative is a nonlocal operator. Studies and results of solutions obtained using the fractional differential equations are more general and are as stable as their integer order counterpart.

There are a lot of approaches to define the fractional differential operator such as Grunwald-Letnikov, Riemann-Liouville, Caputo, and Hadamard. The Riemann-Liouville and Caputo approaches are the most widely used in applications [21, 22, 29, 31]. The fractional derivative is defined via the fractional integral operator. So, we will start by the definition of the fractional integral operator. The fractional integral operator of order (>0) of a function such that is given byThe Riemann-Liouville derivative of order (>0) is given bywhere and

In this paper, we used Caputo approach to define the fractional derivative. It is a modification to the Riemann-Liouville definition. The Caputo fractional derivative of order (>0) is denoted by and is given in the following form:where and .

The following are some important properties of the fractional derivatives and integrals [18]. Let and ; then we have the following:(1)If and , then (2) uniformly, where and (3) weakly.(4)If is constant and , then (5)If is absolutely continuous on , then

In this paper, we proposed a fractional-order model to study the interaction of a system consists of two-prey and one-predator species. In Section 2, we proved the existence and the uniqueness of the solutions of our model. In Section 3, we studied the local asymptotic stability of the equilibrium points of the system. The numerical solution of the fractional-order two-prey one-predator model is given in Section 4.

#### 2. The Fractional-Order Prey-Predator Model

Let represent the first prey herds (gazelles) and be the second prey herds (buffalos) densities, respectively. Suppose represent the predator (lions) density. It is known that the logistic scenario is the most appropriate to describe the growth of the preys. So, the terms and are the growth of the two preys, where the positive parameters and are the intrinsic growth rate of them. The nature requires cooperations between preys against the predator to avoid the predation and to facilitate getting food. This cooperation is represented by the term . Also, the predation to the preys is represented by the terms and . Considering these assumptions, we get the following fractional-order two-prey one-predator system:with the initial values where is the death rate of the predator, , and , , . The constants , and are all positive. In the following, we studied the above model to understand the long time behaviour prey-predator interaction.

Lemma 1. *The initial value problem (1, 2) can be written as the following matrix form: where*

*Definition 2. *Let be the class of continuous column vector where is the class of continuous functions defined on the interval and ,

Lemma 3. *The initial value problem (1, 2) has a unique solution in the region , where and *

*Proof. *Let be a continuous function where is a continuous column vector. Suppose that and are two distinct continuous column vectors solutions of the initial value problem (1, 2) such that , Then, Since for each

Then, we have Then where So, the continuous function satisfies Lipschitz condition and has a unique solution [32].

#### 3. Equilibrium Solutions and Stability Analysis

So far there is no known method for solving nonlinear equations. Therefore, it is difficult and even impossible to find an analytical solution to system (4). So, we will need the qualitative study. Finding the equilibrium points and studying their stability are the most important. To evaluate the equilibrium points of system (4), let where is the equilibrium point of system (4).

Solving the equations of system (12), we get the following equilibrium points. The dormant state , the boundary states , , , , , the two coexistence states , and which exists under the following conditions:To study the stability of these equilibrium points, we have to linearize system (4) and compute its Jacobian matrix:

Substituting by the equilibrium point in the above Jacobian matrix, we getwhich has the eigenvalues . Since we have two positive eigenvalues, then the equilibrium point is unstable. This means that the state of extinction of all species is not acceptable. There will be no life. The Jacobian matrix for the is and their eigenvalues are . Then the equilibrium point is unstable. Similarly, the eigenvalues of the equilibrium point are . So, it is unstable. Similarly, the eigenvalues for the equilibrium point are . So, point is unstable. This means that the state of existence of one prey alone cannot be continuous forever.

For equilibrium point , the Jacobian matrix iswhich has the following characteristic polynomial: where and

A sufficient condition to say that an equilibrium point is a locally asymptotically stable is that all eigenvalues satisfy [33]. This condition implies that the characteristic polynomial of that point should satisfy Routh-Hurwitz conditions [17]. Since and are both positive, then the stability of point depends on the first bracket in (18). Thus is locally asymptotically stable if

Similarly, equilibrium point has the Jacobian matrix: which has the characteristic polynomial: where and Since and are both positive, then is locally asymptotically stable if

The study of the coexistence points is more important. This study gives us the conditions that lead to coexistence between the prey and the predator. The Jacobian matrix for the first coexistence equilibrium point has the characteristic polynomial: where and Since and are both positive and . Thus is locally asymptotically stable.

For the equilibrium point , the Jacobian matrix is given bywhere The above matrix has the following characteristic polynomial: where , , Since , , and , then, from the above, , , and .

*Definition 4. *The discriminant of the polynomial is defined by , where is the derivative of the function and is the determinant of the corresponding Sylvester matrix; then,

Proposition 5. *If the discriminant is positive, then the Routh-Hurwitz conditions (, and ) are the necessary and sufficient conditions for that all eigenvalues satisfy *

Proposition 6. *From the above definition of the discriminant and the coefficients of the characteristic equation, we find that (with long calculations) all the four necessary and sufficient conditions are satisfied. Then equilibrium point is locally asymptotically stable.*

#### 4. Numerical Results

The Adams-type predictor-corrector method for the numerical solution of the fractional differential equations was discussed in [10]. This method can be used for both linear and nonlinear problems. It may be extended to multiterm equations (involving more than one differential operator) too [10].

In this paper, we used Adams-type predictor-corrector method for the numerical solution of our system. First, we will give the Adams-type predictor-corrector method for solving general initial value problem with Caputo derivative: with the initial condition and . We assumed that a set of points , where , are the points used for our approximation and , , The general formula for Adams-type predictor-corrector method iswhereApplying the above algorithm for system (4), we have the following:where Therefore, for ,

We used the step size to be in Figures 1, 2, and 3. In Figure 1, we used the initial values , , and and the parameter values , , , , , and , such that equilibrium point is locally asymptotically stable. We have the approximate solutions between and in Figure 1(a), between and in Figure 1(b), and between and in Figure 1(c) and we showed that the solutions are stable.

**(a)**The approximate solutions between and

**(b)**The approximate solutions between and

**(c)**The approximate solutions between and

**(d)**The stability of approximate solutions , and

**(a)**The approximate solutions between and

**(b)**The approximate solutions between and

**(c)**The approximate solutions between and

**(d)**The stability of approximate solutions , and

**(a)**The approximate solutions between and

**(b)**The approximate solutions between and

**(c)**The approximate solutions between and

**(d)**The stability of approximate solutions , , andIn Figures 2 and 3, we used the initial values , , and . In Figure 2, the constant values are , , , , , and . In Figure 3, the constant values are , , , , , and which satisfied the existence conditions of the locally asymptotically stable equilibrium point We plot the approximate solutions between and in Figures 2(a) and 3(a), between and in Figures 2(b) and 3(b), and between and in Figures 2(c) and 3(c) and we showed that the solutions are stable.

We used the trapezoidal method [34] to compare the results and plot some new figures; see Figures 4, 5, and 6 to show this comparison.

**(a)**The approximate solutions between and

**(b)**The approximate solutions between and

**(c)**The approximate solutions between and

**(d)**The stability of approximate solutions , , and

**(a)**The approximate solutions between and

**(b)**The approximate solutions between and

**(c)**The approximate solutions between and

**(d)**The stability of approximate solutions , , and

**(a)**The approximate solutions between and

**(b)**The approximate solutions between and

**(c)**The approximate solutions between and

**(d)**The stability of approximate solutions , , andWe used the step size to be in Figures 4, 5, and 6. In Figure 4, we used the initial values , , and and the parameter values , , , , , and , such that equilibrium point is locally asymptotically stable. We have the approximate solutions between and in Figure 4(a), between and in Figure 4(b), and between and in Figure 4(c) and we showed that the solutions are stable.

In Figures 5 and 6, we used the initial values , , and . In Figure 5, the constant values are , , , , , and . In Figure 6, the constant values are , , , , , and which satisfied the existence conditions of the locally asymptotically stable equilibrium point We plot the approximate solutions between and in Figures 5(a) and 6(a), between and in Figures 5(b) and 6(b), and between and in Figures 5(c) and 6(c) and we showed that the solutions are stable.

#### 5. Conclusion

In this paper, we studied the existence, the uniqueness, the stability of the equilibrium points, and numerical solutions of a fractional-order two-prey one-predator model. The coexistence equilibrium points and were stable equilibrium points under some conditions at the ordinary differential equation form of the model. But, in our fractional form, we found that the same points are stable without any conditions on the parameters. This means that the predators (lions) can live stably together with the two preys (gazelles and buffalos) without extinction of any of them. This is an example of the equilibrium point which is centre for the integer order system but locally asymptotically stable for its fractional-order counterpart. This means that the fractional-order differential equations are, at least, as stable as their integer order counterpart. So, we recommend to restudy most of the complex models (biological, medical, etc.) in the fractional-order form.

#### Conflicts of Interest

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

#### Acknowledgments

The authors would like to express their gratitude to King Khalid University, Saudi Arabia, for providing administrative and technical support.