Existence theory and numerical analysis of three species prey–predator model under Mittag-Leffler power law

In this manuscript, the fractional Atangana–Baleanu–Caputo model of prey and predator is studied theoretically and numerically. The existence and Ulam–Hyers stability results are obtained by applying fixed point theory and nonlinear analysis. The approximation solutions for the considered model are discussed via the fractional Adams Bashforth method. Moreover, the behavior of the solution to the given model is explained by graphical representations through the numerical simulations. The obtained results play an important role in developing the theory of fractional analytical dynamic of many biological systems.


Introduction
A predator-prey model is a two-component system, where one of them lives at the expense of the other. A diversity of mathematical techniques is applied at modeling a predator-prey system due to numerous factors that may affect its evolution. In this regard there have been introduced some models in [1][2][3][4][5][6][7] in which the first model, which regards in a specified way only substantial phenomena (gluttony and fertility), is of the type where p(t) and q(t) are the number of prey, the number of predators, respectively, and a 1 , a 2 and a 3 are the average of death of predators, the measurement of the tendency of prey to predation, and the predatory capability, respectively. The model (1) has a unique solution. However, the solutions of (1) are not structurally stable w.r.t. perturbation of the initial conditions. Inside the restricted scope of quadratic differential equations, those which cover competition and also predation, must be slightly more realistic. A second model with competition within preys is formulated as ⎧ ⎨ ⎩ p (t) = p(t)[a 1a 2 q(t) + a 5 p(t)], where a 1 a 3 > a 4 a 5 , a 5 > 0 describes the competition of the prey. According to biologically sensible hypotheses, there exists a unique positive solution of model (2) which is asymptotically stable.
Dai, and Zhao investigated the dynamic complexities of a predator-prey model with state dependent on impulsive influences as The authors used the analogue of the Poincaré norm to obtain the existence and stability of the model (3). For details, see [8]. The following dynamic model, which addresses the case of predatory prey with disease, was analyzed by Das et al. [9]: dy dt = (1 -y+x k )r 2 y + a 2 xyα 2 x(z + w)my, where x(0), y(0), z(0), w(0) > 0. The authors showed that the model is globally stable on every side of the internal equilibrium point according to certain standard conditions. So, their analysis shows that the force of infection and predation average are the main parameters on the dynamics of the model. Fractional calculus deals with differentiation and integration involving fractional order, which is advantageous over the ordinary integer order in the explanation of real-world problems, as also in the modeling of real phenomena due to characterization of memory and hereditary properties [10,11]. Further, the integer-order derivative does not describe the dynamics between two various points. Various types of fractional-order or nonlocal derivatives were proposed in the present literature to deal with the reduction of a traditional derivative. For instance, based on a power-law, Riemann-Liouville introduced the idea of a fractional derivative. Afterwards Caputo-Fabrizio in [12] have proposed a new fractional derivative utilizing the exponential kernel. This derivative has a few problems related to the locality of the kernel. Newly, to overcome Caputo-Fabrizio's problem, Atangana and Baleanu (AB) in [13] have proposed a new modified version of a fractional derivative with the aid of a generalized Mittag-Leffler function (MLF) as a nonsingular kernel and being nonlocal. Since the generalized MLF is used as the kernel it is guaranteed to have no singularity. Furthermore, the AB fractional derivative supplies a description of memory as discussed in [14][15][16][17][18][19][20].
Most of the published work describes the mathematical system of predators and prey as a problem of Cauchy type of a system of classical differential equations [21][22][23][24][25]. However, recently, there has been great interest in studying the behavior of the solution for some biological systems using fractional differential equations involving the Atangana-Baleanu operator by several authors for the purpose of investigating several real-world systems and modeling infectious diseases; see [26][27][28][29][30][31][32][33][34][35][36]. Some fractional-order models have been investigated via the new operators recently. For instance its use has been suggested for the dynamics of smoking in [32]. Along the same line, the transference model for the Ebola virus together with AB operator was studied in [31]. A fractional-order model of leptospirosis infection was considered in [26]. The dynamical behavior of coronavirus (COVID-19) epidemic infection model through the ABC derivative has been studied in [33]. Also, the existence results and analytic solutions of fractional-order dynamics of COVID-19 with ABC derivative has been obtained in [34]. There is no literature available on prey-predator fractional models with three species under the aforesaid derivative. Just some fractional models have been found in the previous years; however, they have been confined to a standard fractional derivative. Furthermore, in the presence of the mentioned derivatives, recently some fruitful results have been published in [37][38][39].
Due to the success of this operator in modeling the biological systems and infectious diseases, we have studied the dynamical behavior of the mathematical model which describes three prey-predator species by a nonlocal Atangana-Baleanu-Caputo (ABC) derivative operator with 0 < α ≤ 1 as with the initial conditions where ABC D α 0 + (·) is the ABC fractional derivative of order α, P 0 is the initial population density of prey, S 0 is the initial population density of susceptible predator, and I 0 initial population density infected predator. Here a denotes the saturation constant whereas susceptible predators threaten the prey, b is a search rate of the prey across a susceptible predator, c is the conversion rate of the susceptible predator due to prey, and d is the disease transmission coefficient. The symbol k represents the carrying capacity of the prey population, the proportionality constant is denoted by b 1 , the growth rate of the prey population is represented as r 1 . In the proposed model, m and n define the death rate of sensitive predator and death rate of the infected predator, respectively. Further, we ramark that the right hand sides of our considered model (4) under ABC fractional derivtive are assumed to vanish at zero, (for details, see Theorem 3.1 in [28]).
The main aim of the paper is to demonstrate the existence, uniqueness and Ulam stability of the solution for the model (4)-(5) by using the Picard and fixed point techniques. Moreover, the numerical simulations via the fractional version of the Adams Bashforth technique to approximate the ABC fractional operator are performed. Graphical presentations are also given of the numerical results. This paper is organized as follows: Sect. 1 presents an introduction which contains a survey of the literature. Section 2 consists of some foundational preliminaries related to fractional calculus and nonlinear analysis. The existence and Ulam stability results on a proposed model are obtained in Sects. 3, 4. The numerical solution and numerical simulations of the model at hand are presented in Sect. 5.

Lemma 1 (see Proposition 3 in [40])
The solution of the proposed problem for α ∈ (0, 1] is given by here κ is the Lipschitz constant for Π . If κ < 1 we say that Π is a contraction.

Existence of solutions for the proposed model (4)-(5)
Now, we address the existence and uniqueness results of the model (4)-(5) by utilizing the fixed point technique. Let us reformulate model (4) in the appropriate form where Utilizing Lemma 1, the model (8) can be turned to the fractional integral equation in the sense of AB fractional integral as follows: Theorem 2 The kernels W ( = 1, 2, 3) agree with the contraction and Lipchitz conditions if there exists a constant L such that 0 ≤ L < 1, = 1, 2, 3.

Theorem 3 Assume that the conditions
Then the solution of the fractional model given in (4)-(5) exists and is unique.

Proof 2
The initial conditions and the recurrence form of the model (10) are, respectively, The successive difference between the terms is defined as Clearly Taking the norm of Eqs. (15), it follows from the conditions (11)-(13) that Let us consider P, S and I as bounded functions that comply with the Lipschitz condition. It follows from Eqs. (16) and (17) that This shows the existence for the solutions. Moreover, to prove that Eqs. (18) are solutions for the model (4)-(5), we consider

Now, we consider the conditions
On using recessive techniques, we get As n → ∞, M 1n (t) → 0. In a similar way, we conclude that M 2n (t) and M 3n (t) tends to 0. Next, we address the uniqueness of the solution to the proposed mode (4)- (5). To this end, let P * (t), S * (t) and I * (t) be other solutions. Then It means that From our hypothesis It follows that P(t)-P * (t) = 0. Likewise, we conclude that S(t)-S * (t) = 0 and I(t)-I * (t) = 0.

Ulam-Hyers stability
For the notion of Ulam stability, see [42,43]. The aforesaid stability has been scrutinized for classical fractional derivatives in many of the research articles; we refer to some of them like [44][45][46][47]. Additionally, since stability is a prerequisite in respect of approximate solution, we endeavor on Ulam type stability for the model (4) via using nonlinear functional analysis.
2 Furthermore, one has Note that we will only discuss the first equation from the proposed system and the rest of the equations are similar in technique, i.e. P -P E ≤ λ 1 1 .
Proof 3 Thanks to Remark 1, and Lemma 1, the solution of (21) is given by Also, we have It follows from Remark 1 that (11), the system (4)-(5) will be Ulam-Hyers stable in Ω.

Numerical interpretation and discussion
Now, to present the numerical simulations of the ABC fractional model (4)-(5), we apply the iterative solution contained in (28)- (30). Take the time range up to 100 units. The numerical values of the parameters applied in the simulations are specified in Table 1.
The graphical representations of numerical solution for species P, S, I at various fractional orders, α = 0.4, 0.6, 0.8, 1.0, of the considered model (4) are given in Figs. 1-3, respectively.  From Figs. 1-3, we observe that species I depends on species P and S. Therefore the papulation density of specie P and S gradually go on decreasing with different rate due to the fractional order in the first 50 days. The lower the fractional order, the faster the decay rate and hence the more rapidly the system becomes stable and vice versa. On the other hand, the species I is going on increasing with different rate, the lower the order the slower is the growth rate until it becomes stable and vice versa. The fractional order greatly Figure 3 Graphical representation of numerical solution for specie I at various fractional orders of the considered model (4) affects the stability of the system and also provides the global nature of the dynamics of the considered model.
Here we claim that the established numerical technique is powerful and converges for the ABC fractional derivative. Meanwhile the iterative techniques like perturbation and decomposition methods do not show the perfect behavior for the said derivatives for approximate solutions in many cases.

Conclusion
In this paper, the population density model of prey and sensitive predatory and infected predatory has been studied theoretically and numerically. Theoretically, the existence and stability results in the sense of Ulam-Hyers have been obtained through the help of fixed point theory and nonlinear analysis. Numerically, the approximation solution of the ABC fractional model (4) has been discussed via the use of a fractional Adam Bashforth method. Moreover, the behavior of the solutions of the model (4) has also been explained through graphs using some numerical values for the parameter. The obtained results play an important role in developing the theory of fractional analytical dynamics of various phenomena of real-world problems.