Stability Analysis and Numerical Computation of the Fractional Predator–Prey Model with the Harvesting Rate

In this work, a fractional predator-prey model with the harvesting rate is considered. Besides the existence and uniqueness of the solution to the model, local stability and global stability are experienced. A novel discretization depending on the numerical discretization of the Riemann–Liouville integral was introduced and the corresponding numerical discretization of the predator–prey fractional model was obtained. The net reproduction number R 0 was obtained for the prediction and persistence of the disease. The dynamical behavior of the equilibria was examined by using the stability criteria. Furthermore, numerical simulations of the model were performed and their graphical representations are shown to support the numerical discretizations, to visualize the effectiveness of our theoretical results and to monitor the effect of arbitrary order derivative. In our investigations, the fractional operator is understood in the Caputo sense.


Introduction
Fractional calculus (FC) is a common field trying to understand the real-world phenomena that are modeled with non-integer-order derivatives and it is a field wherein the differentiation and the integrations are done with non-integer order derivatives as well. A fractional derivative, one can understand, is a type of derivative in which the order is non-integer-based but satisfies certain conditions: when the order of the derivative is zero we have the primary function, and when the order is one we converge to the first order integer order derivative [1]. The advantages of the fractional derivatives are the memory impact and the illustrative physical properties that are conserved. Using these types of operators, more effective and up-to-date studies have been revealing over time. In this context, fractional calculus theory and its illustrative applications are attracting attention all over the world day by day. New fractional operators that have different features have been defined and have been used extensively to model real-life problems. The emergence of the new operators in the literature can be considered as a result of the reproduction of new problems that model different types of real-life events. Fractional derivative operators that address the kind of nonlinear differential equations can be stated as non-local. There exist nowadays, many types of fractional derivatives with and without singular kernels. The fractional derivative begins with Leibniz's question in 1695. The list of the existing fractional derivatives is very long. With singular kernels, we have the Caputo-derivative [2], the Riemann-Liouville derivative [2] and the Katugampola derivative [3].

Some Preliminaries
In this section, we give the fundamental definitions that can be used throughout the paper. These definitions generally explain the fractional derivative in the power kernel sense. Definition 1 ([1]). The Riemann-Liouville (R-L) representation of fractional integral operator of order γ > 0 of a function ϕ : (0, ∞) → R is given by  (2) Definition 3 ([1]). The Caputo fractional operator of order γ > 0 of a function ϕ : (0, ∞) → R is given by where the operator C 0 D γ t satisfies:

Definition 4 ([1]
). The Laplace of the Caputo fractional operator of a function ϕ(t) of order γ > 0 is presented with We introduce two lemmas which will be used to establish local stability and global stability, respectively. The first lemma is called the Matignon criterion and the second one is the Lyapunov characterization for global stability.

Lemma 1 ([68]
). The fractional differential equation C 0 D γ t x = Px, with P ∈ R n×n , x(t 0 ) = x 0 , 0 < γ < 1 and x ∈ R n , is local asymptotically stable if only if where spc(P) is considered as the spectrum of the matrix P.

Lemma 2 ([69]
). We assume the vectors w ∈ R n which are differentiables. Under the assumption t ≥ t 0 , we have the following condition

Fractional Predator-Prey Model with Caputo Derivative
In this paper, we present the predator-prey model in the context of the fractional operator. Note that the fractional-order derivatives are more accurate at modeling biological processes since these types of operators consider the memory impact which gives more realistic results in real-life models. It is the deterministic property of the dynamical system. The following dynamics represent the fractional equation proposed in this section [53,70,71] where C 0 D γ t represents the Caputo fractional derivative which is given in Definition 3. We consider the initial conditions defined by where x denotes the prey species; y represents the predator species; the prey increases logistically with the growth rate denoted by r and the carrying capacity k; a represents the rate of predation; c denotes the efficiency of predation; d is the mortality rate of the predator species; and H (y) represents the harvesting function. Now, we consider the model of Equations (7) and (8). Firstly, we describe harvesting function H (y) of the predators in our model, Equations (7) and (8), which has the following form Using the model in Equations (7) and (8) as our baseline model, we suppose that harvesting takes place, but only the predator population is under harvesting, and introduce harvesting function H (y) of the predator to prey-predator model in Equations (7) and (8) for discussing its dynamical features. We assume that the harvesting rate is proportional to the predator population size until it reaches a threshold value due to limited facilities of harvesting or resource protection. Let us show the harvesting threshold value as h = my 0 ; thus Equations (7) and (8) can be written as the following When y 0 < y, then the system in Equations (7) and (8) turns to the following

Positivity and Boundedness
In this subsection, the positivity and boundedness of the solution for the proposed model (Equations (7) and (8)) are given. Let R 2 + = χ(t) ∈ R 2 : χ(t) ≥ 0} and χ(t) = [x (t) , y (t)] T . Theorem 1. The solution of the proposed fractional-order model (Equations (7) and (8)) along initial conditions (9) is bounded in R 2 + . Moreover, the density of the population remains in a nonnegative region.
Proof. Let the function W(t) = x + 1 c y and λ be a positive constant. Applying the Caputo derivative, we have the following relationship where λ < d + m. From the property in [72] and using comparison principle, since t converges to infinity, we have the following relationship where E γ,β (.) is the Mittag-Leffer function of two parameters. Finally, the following set in the domain R 2 + is positively invariant We have replaced the classical derivative by the Caputo derivative in the study. It is important to justify the replacement and to prove the physical meanings of the new model. Another point is also to show that the solution to the new model exists and is unique. All these points will be discussed in the next sections.

Qualitative Properties of the P-PM
In this Part, we give the qualitative properties of the solutions of the predator-prey model which is given in the system (7) and (8). Firstly we start by taking the Riemann-Liouville integral which is given in Definition 1 of both sides the mentioned system and we get which gives the following Volterra-type integral equations: Let us define the following kernels as Then the following theorem arises: The kernels ϕ and φ satisfy the Lipschitz assumptions and contractions if the following inequality is verified: where x ≤ q, y ≤ l, z 1 = r + al + 2qr/k, q, l ≥ 0 and z 2 = acq + d, or z 2 = acq + d + m, for the constant harvesting rate or depending on the predator population, respectively.
Proof. Let x 1 and x 2 be two functions for the kernel ϕ; and y 1 and y 2 be two functions for the kernel φ. Then we have and or where . is the Euclidean norm, x ≤ q, y ≤ l, z 1 = r + al + 2qr/k and z 2 = acq + d, or z 2 = acq + d + m, for the constant harvesting rate or depending on the predator population, respectively. Therefore, the Lipschitz conditions are satisfied for kernels ϕ and φ, and if 0 ≤ z 1 , z 2 < 1, then z 1 and z 2 are also contractions for ϕ and φ, respectively. This proofs the theorem.
By considering the kernels ϕ and φ, we can rewrite the system which is given in Equation (16) as follows: We can proceed with the following recursive formula where x 0 (t) = x (0) and y 0 (t) = y (0) . Then we can write where x n (t) = ∑ n j=1 Ψ n (t) and y n (t) = ∑ n j=1 Φ n (t). By taking the norm of both sides of Equation (24), we have (25) Since the kernels satisfy the Lipschitz condition (see Theorem 2), we get Then we achieve from the last inequality These results give us the following theorem: The predator-prey model defined by the fractional operator with the power kernel has a solution under the condition that we are able to find t max holding: Proof. Considering the functions x (t) and y (t) are bounded and their kernels ϕ and φ hold the Lipschitz condition, we can give the following by taking Equation (27) into account, Now we show that the functions in Equation (29) are the solutions of the given predator-prey model. We suppose where p n and q n are remaining terms. Then we will demonstrate that the terms which are given in Equation (30) hold that p ∞ (t) → 0 and q ∞ (t) → 0. Since we have repeating this process recursively, we get Considering these last two inequalities at t max point, we have For the last step, after applying the limit to both sides of the last inequalities as n → ∞, and by taking into account the results of Theorem 2, we get p ∞ (t) → 0 and q ∞ (t) → 0.

Theorem 4.
The predator-prey model defined by the fractional operator with the power kernel in Equations (7) and (8) has a unique solution.
Proof. Let say there exists another solution of the system, namely, x 1 (t) and y 1 (t). Then we can write If we apply the norm to both sides of Equation (33), we obtain Since the Lipschitz condition is satisfied by the kernels ϕ and φ, we can write which gives Hence, we have x (t) − x 1 (t) = 0 and y (t) − y 1 (t) = 0 which gives x (t) = x 1 (t) and y (t) = y 1 (t) . This concludes that the model has a unique solution and proofs the theorem.

Stability Analysis of the Predator-Prey Model
In this section, we determine the equilibrium points and study their stability. The local stability of the equilibrium points will be examined by using the Jacobian matrix, and the global asymptotic stability will be studied by constructing a Lyapunov function.

Existence of Equilibria
In this subsection, we examine the existence of all nonnegative equilibria and present our results of the existence of positive equilibria as follows. The equilibrium points of the P-PM described by Equations (7) and (8) are obtained by solving the equations represented by Before going further, we proceed with the net reproduction number for the predator population as the expected number of predator individuals producing as the predator population is introduced into a stable prey population [53]. Thus we adopt the following procedure. Fractional differential Equations (7) and (8) can be written in the form where ξ = (y, x) and the matrix F and V are described as follows The Jacobian matrixes of the functions F and V at the predator-free point (k, 0) gives the following matrixes F and V. We have the following Finally, the reproduction number is obtained by determining the spectral radius of the matrix FV −1 : The reproduction number is significant in the classification of the biological models. This number, in general, tells us the number of species that can be infected by a single infected person. For the control of the biological model, notably, its usage in the stability analysis is essential. As we will notice, for the rest of the paper, the stability analysis of the extinction point, the predator-free equilibrium and the non-trivial equilibrium will be focused according to the reproduction number R 0 . Many other interpretations of the net reproduction number exist as well.
Now we proceed with the system in Equation (11) which considers the linear predator harvesting strategy to present its equilibria. The resolution of Equation (11) gives three different equilibrium points, which are the extinction point (0, 0), the predator-free equilibrium (k, 0) and the non-trivial predator-prey equilibrium point d+m ac , r a 1 − d+m akc ; it is straightforward to reach to the first two equilibria. In the subregion A with the linear predator harvest strategy, i.e., when 0 ≤ y ≤ y 0 , a positive equilibrium, satisfies the following system which follows which means that there exists a positive coexistence equilibrium in the region when 0 ≤ y ≤ y 0 if r − ay 0 r ≤ 0, or r − ay 0 r > 0 and R 0 ≤ 1 1 − ay 0 r = r r − ay 0 .

Stability of Equilibria
For the investigation related to the stability of the equilibrium points, we give the form of the Jacobian matrix; that is, For the extinction point (0, 0), we fix x = 0 and y = 0, substituting them into the Equation (43), andwe obtain the Jacobian matrix computed at the (0, 0) given by the following matrix The eigenvalues of the Jacobian matrix are given by λ 1 = r and λ 2 = −d − m. In the context of fractional order derivative, we evaluate arg (λ 1 ) = 0 < γπ/2 and arg (λ 2 ) = π > γπ/2, for all γ ∈ (0, 1). Therefore, the condition described in Lemma 1 is not satisfied; that is, the extinction point (0, 0) is unstable.
The local stability of the predator-free equilibrium is described in the following procedure. The Jacobian matrix defined in Equation (43) evaluated at the point (k, 0) is determined by the matrix Thus, the eigenvalues of the Jacobian matrix are given by λ 1 = −r and . Now, in the context of fractional order derivative, we evaluate arg (λ 1 ) = π > γπ/2 and arg (λ 2 ) = π > γπ/2, for all γ ∈ (0, 1), and when the reproduction number satisfies the condition R 0 < 1. Thus, the predator-free equilibrium is locally asymptotically stable if the condition R 0 < 1 is held. We also notice when R 0 > 1, then arg (λ 2 ) = 0 > γπ/2, which in turn is impossible. Thus, the predator free-equilibrium should be unstable for all fractional time order satisfying the condition γ ∈ (0, 1).

Theorem 5.
The coexistence equilibrium E * (x * , y * ) of the proposed fractional predator-prey model is locally asymptotically stable if R 0 > 1; otherwise, it is unstable.

Proof.
To study the stability criterion of coexistence equilibrium, the Jacobian matrix which has been calculated in Equation (43) is considered in the equilibrium point E * (x * , y * ) . Thus, it follows We can evaluate the eigenvalues by solving the following corresponding characteristic equation which gives the equation of the form and we then get from the last equation the roots as . This means that all two eigenvalues have negative real parts as long as R 0 > 1. Therefore, we conclude from this fact that the coexistence equilibrium E * (x * , y * ) is locally asymptotically stable if R 0 > 1. Moreover, if R 0 < 1, one of the eigenvalues is positive which means that the coexistence equilibrium is unstable. Now, we prove that the coexistence equilibrium point is globally asymptotically stable. For that proof we take into account the Lyapunov direct technique. To arrive at our end, we propose the following Lyapunov function Using Lemma 2, the fractional derivative along the trajectories of the Lyapunov function yields that Using the fact that x * = d+m ac and y * = r a 1 − d+m akc , Equation (50) can be written in the following form: which implies that the non-trivial equilibrium point d+m ac , r a 1 − d+m akc is globally asymptotically stable.

Numerical Scheme of the Predator-Prey Model
This section addresses the numerical discretizations of the fractional predator-prey equations represented by Equations (7) and (8). The algorithm adopted in this section begins with the solutions of the fractional differential equation of Equations (7) and (8) in terms of the R-L integral. We mainly use the implicit discretization method to construct the scheme. The solutions of the fractional predator-prey model described by Equations (7) and (8) are represented as follows: At the point t n , the solutions represented by Equations (52) and (53) can be discretized as the following form We set the step-size as z and the grid step as t n = nz. Thus the Riemann-Liouville integral can be discretized in an implicit sense to the following form.
where the parameters are represented by the expressions enumerated as follows and for n = 1, 2, . . ., we set the following parameters in other cases as the following form: and a (γ) Substituting Equations (56) and (57) into Equations (54) and (55), respectively, the final numerical discretization of the P-PM in the context of the Caputo derivative is presented by the following form [73].
where the following relationships describe the numerical discretizations of the intermediary functions Ψ and Φ.
We set x(t n ) and y(t n ), the numerical approximations of the fractional predator-prey model, and x n and y n , their associated exact solutions. Note that our numerical discretization follow the following errors terms: which converge to zero as the step-size z converge to zero too. In the next section, we give the graphical representations to support our implicit numerical discretization for the fractional predator-prey model constructed with the Caputo derivative.

Numerical Simulation of the Implicit Scheme
We illustrate the numerical discretizations of the predator-prey model described in the previous section. For the graphical representations, we consider different contexts by considering the different values of the parameters of the model. We, firstly, analyze the solution according to the fixed harvesting rate (case 1). In first cases, we fix the following assumptions [53]: the growth rate r = 0.03 and the carrying capacity k = 0.25, the rate of predation is a = 0.5, the efficiency of predation is c = 0.4, d = 0.01 is death rate of the predator species and γ = 0.95 is the order. In the first subcase, we fix the harvesting rate as h = 0.00033; in Figure 1, we depict the evolutions of the predator and the prey species in time. We notice the trajectory describes a cycle and converges to the non-trivial equilibrium point. In Figure 2, we depict the dynamics of the P-PM under no harvesting rate (h = 0). We notice the same dynamics as in the presence of harvesting rate; the solutions describe a cycle and converge to a non-trivial equilibrium point. We notice when the harvesting rate varies and exceeds h = 0.00035, in Figure 3, we see the line instead of the cycle. It corresponds that the non-trivial equilibrium point is not stable. Furthermore, the obtained solutions for the predator-prey model are unrealistic, since the species can not be negative; they do not make sense. In conclusion, the harvesting rate has a significant impact on the dynamics of the P-PMs. We note after certain values that the behaviors of the solutions are in contradiction with the real possible phenomena. Thus, it is crucial to control the harvesting rate into the predator-prey models. There exist many methods to control it, such as the maximization principle using the Hamiltonian approach. This problem is not addressed in this paper.
Let us give illustrative results to support our results. We change the values of the parameters, and we suppose the following assumptions [53]: the growth rate r = 0.004 and the carrying capacity k = 0.25, the rate of predation a = 0.1, the efficiency of predation c = 0.4, the death rate of the predator species d = 0.001, the harvesting rate h = 0.000005 and the order is maintained at γ = 0.95.
We notice in Figure 4 that the evolutions of the predator and the prey species have been depicted in time. We see the trajectory describes a cycle and converges to the non-trivial equilibrium point.
In conclusion, after variation in the parameters, the solutions respect the cycles; there are not many changes in the behaviors of the solutions. Our fractional model is, in general, stable but strongly depends on the values of the harvesting rate, because when the harvesting rate value exceeds certain values, the cycle is displayed.  To support our numerical scheme, we consider the second case where the harvesting rate is expressed as the linear representation of the predator population, i.e., h = my. We set the following assumptions [53]: the growth rate r = 0.1 and the carrying capacity k = 0.5, a = 0.25 represents the rate of predation, c = 0.8 denotes the efficiency of predation, d = 0.01 is death rate of the predator species, the order γ = 0.95 and m = 0.0003. In Figure 5, we represent graphically the dynamics of the predator versus prey. In Figure 6, we represent graphically the dynamics of the predator in time. In Figure 7, we represent graphically the dynamics of the prey in time. In the figures of the first case, the order converges to the classical order, i.e., to 1. We now depict the figures when the order is arbitrary in the interval (0, 1) which we consider it as γ = 0.85. In Figure 8, we represent graphically the dynamics of the predator versus the prey.
In Figure 9, we represent graphically the dynamics of the predator in time.
In Figure 10, we represent graphically the dynamics of the prey in time.
We notice by comparing the plots in Figures 6, 7, 9 and 10 that the fractional order derivative has a significant impact on the dynamics of the suggested predator-prey model represented by Equation (11). In general, it has an acceleration effect on the model process. It is remarkable that the model has not been previously considered in terms of fractional derivatives. For this reason, making a comparison of it with those existing in the literature is a bit difficult. A possible comparison can be done with the prey-predator model which was constructed by integer-order derivative in [53]. In terms of comparison we notice that the fractional order derivative generates an acceleration effect in the dynamics. These mentioned differences can be seen when looking at the comparisons of

Concluding Remarks
In this paper, modeling and analysis of the fractional order P-PM which contains the harvesting rate have been provided. Since the harvesting rate has a significant impact on the dynamics of the predator-prey models, we have demonstrated after certain values, the behaviors of the solutions are in contradiction with the real possible phenomena. Moreover, the basic reproduction number R 0 has been computed by the next generation matrix method which performs as a threshold parameter in the disease transmission and determines whether the disease persists or vanishes from the population. The existence and uniqueness of the solutions of the proposed fractional system have been examined. Additionally, the stability conditions of the equilibrium points for the fractional system have been discussed. Meanwhile, the global dynamics of the equilibria have been obtained by the Lyapunov functional approach method. It is well-known that the infection spreads in the population when R 0 > 1. A new numerical algorithm based on the numerical discretization of the Riemann-Liouville integral has been introduced and it has been successfully applied for the corresponding numerical solution of the proposed predator-prey fractional-order model to carry out the numerical simulations for different values of the fractional order γ. The model introduced in the paper is new in the context of fractional order derivatives; therefore, to analyze the dynamics of the predator and prey, we need to get the solutions to the proposed model. Since it is straightforward that the analytical solution to the proposed model is not possible, an alternative way to deal with this issue is to construct the numerical scheme. This is because we used an effective and accurate numerical scheme including the discretization of the Riemann-Liouville integral, and by using this scheme, we have managed to obtain numerical solutions to the aforementioned problem. Finally, it has been demonstrated that physical processes are better described using the derivative of fractional order which is more accurate and reliable in comparison with the classical order case. Hence, we replace the integer order time derivative with the Caputo type fractional order derivative.