On approximate solutions for a fractional prey–predator model involving the Atangana–Baleanu derivative

Mathematical modeling has always been one of the most potent tools in predicting the behavior of dynamic systems in biology. In this regard, we aim to study a three-species prey–predator model in the context of fractional operator. The model includes two competing species with logistic growing. It is considered that one of the competitors is being predated by the third group with Holling type II functional response. Moreover, one another competitor is in a commensal relationship with the third category acting as its host. In this model, the Atangana–Baleanu fractional derivative is used to describe the rate of evolution of functions in the model. Using a creative numerical trick, an iterative method for determining the numerical solution of fractional systems has been developed. This method provides an implicit form for determining solution approximations that can be solved by standard methods in solving nonlinear systems such as Newton’s method. Using this numerical technique, approximate answers for this system are provided, assuming several categories of possible choices for the model parameters. In the continuation of the simulations, the sensitivity analysis of the solutions to some parameters is examined. Some other theoretical features related to the model, such as expressing the necessary conditions on the stability of equilibrium points as well as the existence and uniqueness of solutions, are also examined in this article. It is found that utilizing the concept of fractional derivative order the flexibility of the model in justifying different situations for the system has increased. The use of fractional operators in the study of other models in computational biology is recommended.


Introduction
The study of the ecosystem using mathematical models has been one of the most crucial joint research fields between applied mathematicians, engineers, economists, ecologists, biologists, and epidemiologists for decades. It should be noted that differential equations have always been used in modeling and studying phenomena. The reason for this efficiency is that the future of systems is influenced by the past and present situations of understudied phenomena. In recent years there has been a strong tendency to use differential equations to model real-world problems. Mathematical biology is one of the most important and growing branches whose studies rely more on the use of differential equations than other fields [1][2][3][4][5][6][7][8][9][10][11][12]. One of the widespread applications of dynamic systems in ecology is to model the growth of monogamous and multicultural populations. The systems proposed in ecology determine the interaction of species with each other and with their environment, as well as the distribution and structure of populations. Such interactions may take place over a wider range of spatial and temporal scales. One of the most important types of interactions that affect the qualitative behavior of all species is hunting. For this reason, predator and prey models have been the focus of all ecological studies. In any ecosystem, different species of animals are usually in constant interaction with each other. In this case, the population dynamics of each species is affected by other species. In fact, there is a general set of interacting species called the food web [13][14][15][16]. One of the most important models in these prey and predator species is Lotka-Volterra systems [17]. In conventional models, in order to simplify the surveys, the interaction between a limited number of populations is usually examined. In any ecosystem, there are usually many species of animals in constant interaction with each other. In this case, the population dynamics of each species is affected by other species. The diversity of animals in the area in question always complicates the dynamic systems. Therefore, in conventional models, in order to simplify the analysis, the interaction between a limited number of populations is usually considered.
As we know, fractional calculus can be considered as an advanced generalization for classic differential calculus of integer order [18][19][20][21][22][23][24][25][26][27][28][29][30][31]. Nowadays, fractional differential calculus is one of the most widely used aspects of differential calculus, which is less used in practical aspects of real-life problems compared to the classical concepts. Perhaps the most important reason for this shortcoming can be considered the lack of objective background to describe practical phenomena, as well as the lack of tools and definitions needed in this area. However, this trend has accelerated in recent years [32][33][34][35][36][37]. And these days, its use in various branches of science, including mathematics, physics, and engineering has been highlighted [29,38,39]. One of the most important factors in the efficiency of such derivatives concerning derivatives is the correct equipping of these definitions with the concept of memory, which causes the basic information related to the phenomenon to be preserved and used during the study time. This fundamental advantage is very important when applying them to the structure of dynamic systems.
A review of the literature reveals that the dynamic models in biology are one of the most fundamental applications of fractional differential calculus with a wide range of potential theoretical and practical aspects. For instance, in [40] the authors have employed the Caputo fractional derivative to propose a new fractional-order predator-prey model with group defense. In [41], the authors have proposed an efficient numerical method to handle a predator-prey model with Beddington-DeAngelis functional response and fractional derivatives with Mittag-Leffler kernel. In [42] one of the Kolmogorov model variants called a Gauss-type predator-prey model involving the Allee effect and Holling type-III functional response in a fractional perspective has been studied. The authors of [43] have utilized a novel discretization based on the numerical discretization of the Riemann-Liouville integral to manipulate a fractional predator-prey model with the harvesting rate. To see more models, please refer to [44][45][46][47][48][49].
These research activities give us the idea of using the tools in fractional differential calculus in computational biological models. In this paper, we aim to incorporate the Atangana-Baleanu fractional derivative [50] in the structure of a prey and predator model. One of the key reasons for this choice is that this fractional order derivative has the ability to store and use the necessary system information from the beginning to the desired time [29,38,39]. This useful feature can play a key role in the study of dynamic systems related to computational biology. The paper is structured as follows. Some mathematical background, mainly about recent definitions on fractional calculus model, is formulated in Sect. 2. The proposed system is analyzed in Sect. 3. In Sect. 4, we investigate some mathematical analysis aspects of the model, including the calculation of equilibrium points of the system, the existence and uniqueness of the solution for the model different equilibrium points, and their stability. Then, we present the corresponding numerical simulations in Sect. 5. Finally, some concluding remarks are stated in the last section of the paper.

Mathematical backgrounds
In this section, we have an overview of some definitions and basic features of fractional differential calculus.

Definition 1
The Caputo (C) fractional operators of order ν for n ∈ N are respectively defined as [19] By taking the Laplace transform from Eq. (1), one achieves Definition 2 The Caputo-Fabrizio (CF) fractional operators of order ν are respectively defined as [51] where M(ν) = 2 2-ν .
The following are some of the basic and widely used features of AB fractional operators: • Both AB derivative and integral operators have linear properties as follows: • The fractional AB derivative operator satisfies the Lipschitz condition where p(t) = max 0≤t≤t f |p(t)|. • For a given function p(t) and the corresponding norm • The combination of the integral and derivative of AB-type operators yields • Utilizing the Laplace transform on fractional definition stated in Eq. (8) yields

Structural expression of the system
To describe the model, let us assume a three-species food-web system of X 1 (t), X 2 (t), and X 3 (t) species. In this model, X 1 (t) and X 2 (t) are the two competing species bearing intrinsic growth rates. Also, it is assumed that their intrinsic growth rates and carrying capacities are r i and K i (i = 1, 2), respectively. Another symbol used in this model is X 3 (t) which preserves the predating over X 2 (t) with Holling type II response. The species X 1 (t) is assumed to benefit from the presence of X 3 (t), as such X 1 (t) is commensal of X 3 (t). The proportion p of species X 2 (t) is refuged from predation. The coefficients α j (i, j = 1, 2; i = j) stand for the inter-species competition coefficients. Finally, the commensal coefficient is denoted by δ 13 .
The following nonlinear differential equation system can be defined to express the interactions and features of this described ecological model: subject to the following nonnegative initial conditions X 1 (0), X 2 (0), X 3 (0) > 0. Now, in order to make system (16) dimensionless, we consider the following new variables in the system: Taking variables (17) and system (16) into account, the following dimensionless form of the problem is obtained [53]: Replacing the standard derivative with the AB fractional derivative D ν AB , we modify model (18) as the following fractional-order version: 4 Some theoretical features of model (19) Some theoretical features of the fractional model are explored in this section.

The equilibrium points of the model
Determining the equilibrium points of a dynamic system is of great importance in better analyzing and understanding the behavior of the model. Finding equilibrium points is equivalent to finding positions for which all derivatives in model (19) are equal to zero. By following this procedure, the corresponding equilibrium points of the system are determined as follows: In order to analyze the local stability of the possible biological equilibrium points, the Jacobian matrix corresponding to model (19) at the equilibrium point B i = (X * 1,i , X * 2,i , X * 3,i ) is calculated as follows: Necessary and sufficient conditions for the existence and stability of all the equilibrium points of B 1 -B 6 have been discussed in reference [53].

Existence of the solution
In this section, we aim to prove that the understudied fractional model in the paper possesses at least one solution. To this end, we first apply the β integral operator defined in Eq. (9) on system (19) and get where (t) = [X 1 (t), X 2 (t), X 3 (t)]. Now, let us define the operator (t) = [Q 1 ( (t)), Q 2 ( (t)), Q 3 ( (t))], where Moreover, 0 = [X 1 (0), X 2 (0), X 3 (0)] is assumed. Using the above notations, Eq. (21) is reformulated as follows: Now, inspired by (23) and starting from 0 (t) = (0), we define the following iterative scheme: Subtracting two consecutive terms gives Now, let us define ς n (t) = n (t)n-1 (t). So, one gets As a consequence, one obtains Assuming that operator K satisfies the Lipschitz condition, it reads Finally, we have the following inequality: Now, replacing ς n-1 (t) by its value, we get And finally the following result is derived: By choosing the general structure for (t) is suggested: where μ n (t) → 0 when n(t) → ∞. Thus, Utilizing the norm on both sides of the above equation, we obtain If n → ∞, then right-hand side vanishes, we reach the following conclusion: This result ensures the existence of solution (t) in the system.

The uniqueness of the solution
Now, let us assume that the system possesses two possible distinct solutions of (t) and Y(t). Hence we must have . . .

A numerical method
Designing numerical algorithms to solve problems is always one of our most urgent needs in the face of presenting new concepts in differential calculus of fractional order. In addition, these methods can provide an accurate description of the theoretical behaviors of problems, which is important from different aspects.
In this section, with an idea of what has been studied in article [54], an efficient numerical approximation form for solving fractional problems will be presented to the following form: By applying the integral operator defined in Eq. (9), the Volterra integral equation is resulted as follows: Setting t = t n = t 0 + n t in (36) yields (37) In this step, it is necessary to approximate the integrand function W (θ , X(θ )) using linear interpolation to the following form: (38) where X i = X(t i ). Furthermore, substituting this linear approximation for W (θ , X(θ )) in integral (37) along with performing symbolic calculations with the help of computational software, the following structure in calculating the approximate solution of the problem is introduced [41,[55][56][57]: where ς n = (-1 + n) ν+1n(-1 + nν) (ν + 2) , , i = 1, 2, . . . , n -1.
After determining the general structure for the approximate method in (39) and (40), the method can be employed in solving problem (19). Hence the following iterative structures will provide the desired numerical results for the problem: By running these iterative schemes, approximate solutions to the main fractional problem of the paper are obtained. By observing the obtained relations, it can be seen that these relations are implicit formulas in determining approximate solutions to the problem. Therefore, it seems necessary to use an additional method in solving these nonlinear algebraic systems during calculations. One of the main suggestions in this regard is to use Newton's method [54].

Figure 1
Dynamics for system (19) for different parameters in Set 1

Figure 2
Dynamics for system (19) for different parameters in Set 1

Simulation results and discussion
The fractional-order system (19) is solved using the proposed iterative scheme outlined in (41). The model has been solved for different values of ν, which are 0.88, 0.878, 0.906, 0.934, 0.962, and 0.99. Now, let us consider some different groups for the parameters in the model.  Figures 1 and 2 show the numerical approximate solution of system (19) using the initial condition (X 1 (0), X 2 (0), X 3 ) = (0.75, 0.5, 0.25). It is seen that the system is locally asymptotically stable around the interior equilibrium point Figures 5 and 6 present the numerical simulations of system (19) with the initial condition (X 1 (0), X 2 (0), X 3 ) = (0.75, 0.5, 0.25). It is seen that the system is locally asymptotically stable around the X 1 -free equilibrium point   Figures 7 and 8 show the numerical approximate solution of system (19) using the initial condition (X 1 (0), X 2 (0), X 3 ) = (0.75, 0.5, 0.25). It is seen that the system is locally asymp- In this part, we seek to study the sensitivity analysis of some parameters in the model. In this case, we fix all the other parameters, and we examine the effects of changing a parameter on the system behavior. The results show that for smaller values of m 1 , the system solution shows more oscillating behaviors and the point equilibrium will be unstable. But as the parameter increases, the system solution quickly converges to its equilibrium point.
Further, let us examine the effect of different values for parameter p on the solution. In this case, we take the parameters in  The results show that, for smaller values of p, the system solution shows more oscillating behaviors and the point equilibrium is unstable. But as the parameter increases, the system solution quickly converges to its corresponding interior equilibrium point.

Conclusion
Employing mathematical modeling using differential equations helps us better understand the behavior of dynamic biological systems. In addition, the use of fractional order differential equations can address some of the shortcomings observed in the modeling of biological systems related to the concept of memory. This has made the importance of using this tool even more obvious for researchers. From a numerical point of view, the use of fast supercomputers with high processing speeds allows us to provide more advanced numerical methods in solving various types of differential equations, including equations of fractional order. In this paper, the AB fractional derivative is employed to study some computational aspects of a three-species prey-predator model in mathematical biology. One of the most prominent features of taking this fractional derivative operator is using a nonsingular kernel in the structure of the operator. Some theoretical requirements for the model, such as the existence and uniqueness of the solutions, are also discussed in this paper. Preservation of basic information in biological systems is one of the most important features necessary for modeling these phenomena. Since the concept of memory in fractional derivatives can provide us with such facili- ties, the use of this important tool in biological modeling will be of great importance.
In addition to the constant coefficients in the model, the fractional-order parameter ν is a useful tool for calibrating the results with real problem data. As we observe in the numerical simulations, the employed fractional operator provides all the expected theoretical aspects of the studied model. From obtained theoretical and numerical investigations one concludes that the inter-species competition coefficients, prey refuge rate, half saturation constant, the death rate of predator species, and conservation rate of prey species have a significant role in the stability of the bio-diversity of an ecological system.