A predator–prey model involving variable-order fractional differential equations with Mittag-Leffler kernel

*Correspondence: jgomez@cenidet.edu.mx; dr.zareenkhan@ymail.com 3CONACyT-Tecnológico Nacional de México/CENIDET, Interior Internado Palmira S/N, Col. Palmira, C.P. 62490, Cuernavaca, Morelos, México 4College of Science, Department of Mathematical Sciences, Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia Full list of author information is available at the end of the article Abstract


Introduction
During the past decades, several mathematical models have been investigated to model prey-predator systems [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. The understanding of the relationship between herbivores and plants is extremely important in the behavior of the ecosystems. Those studies have been implemented in dynamical systems of organic models incorporating integer-order differential equations. However, in these models, the effects of long-range memories are neglected. Recently, mathematical systems with fractional order (FO) have become more worthy than classical systems as FO models provide the description of the memory effects [17][18][19][20][21][22][23][24][25][26][27][28][29]. Owolabi et al. in [30] studied local and global stability of the fractional dynamical system of prey-predator with Holling type-II involving a time delay. Javidi and Nyamoradi provided a detailed study of the local stability of the FO predator-prey model [31]. Alkahtani in [32] investigated the FO triadic predator-prey model. This system was produced utilizing the latest FO differentiation related to the function of generalized Mittag-Leffler (GML). In [33], the authors proposed harvesting terms of a FO delayed predator-prey system. The stability of the model was studied and numerical results were presented to verify the theoretical outcomes. In [34], a FO singular predator-prey model with Holling type-II functional response was studied. The authors investigated the local stability and the solvability condition of this model. Numerical simulations were obtained for different values of the order of the FO derivative. Owolabi in [35] analyzed a FO predator-prey model. The Riemann-Liouville operator based on the pseudo-spectral technique was considered.
In this work, we develop a FO predator-prey system with constant and variable-order Atangana-Baleanu FO operator in the Liouville-Caputo sense. We assume the deterministic integer-order predator-prey model analyzed and studied by Li in [36]. Our system assumes a more realistic prey-predator model with two properties with and without time delay. The mathematical model is built using this FO derivative with the generalized Mittag-Leffler law as a kernel caused by nonlocality of the model and broad usability of the GML function. The Mittag-Leffler kernel applied in Atangana-Baleanu FO differentiation ordinarily exists in various physical models as a generalized exponential decay and as power-law asymptotic for a very large time.
The paper is organized as follows: In Sect. 2, the FO model involving Atangana-Baleanu-Caputo (ABC) operators is disputed. Results are displayed in Sect. 3 numerically, and conclusion is given in Sect. 4.

Predator-prey model
In this paper, we study the mathematical system proposed by Li in [36] and assume Atangana-Baleanu FO operator with constant and variable order. We aim to manufacture a superior estimate model to the real dynamics of a predator-prey system involving a prey refuge among heterogeneous populations. In this investigation, two FO models are treated. The first one depends on the Atangana-Baleanu-Caputo FO differential equations with constant order. We present the existence of positive solutions for the unique proposed design. We study the uniqueness and existence of the solutions. Furthermore, we consider this model involving FO differential equations of Atangana-Baleanu-Caputo type with variable order. Variable order is fixed as a smooth function to specify time differentiation bounded in (0;1]. To get simulation solutions for the suggested system, we employ the Adams process. Finally, the second one assumes a FO model with time delay; in this scenario, we portray the computative technique focused on the Adams-Bashforth-Moulton method to clarify FO delay differential equations relative to variable order through the Atangana-Baleanu-Caputo FO derivative.
In the first case, the system considered is defined as [36] ABC 0 with the initial conditions x(τ 0 ) = x τ (0) and y(τ 0 ) = y τ (0) . Fractional-order operator used in the proposed model (2.1) is in the ABC sense, ABC 0 D τ . FO derivative with constant order is defined as [37] ABC 0 where the function ( ) represents normalization, 1 = (0) = (1). FO operator is using the Mittag-Leffler rule as a nonsingular and nonlocal kernel, where E is a Mittag-Leffler function of one parameter 3) The Atangana-Baleanu FO integral of function f (τ ) with order is characterized as [37] AB System (2.1) can be made more realistic as the prey and the predators should not follow the same FO dynamics. For this reason, we provide two various orders of the fractionaldifferential derivatives ∈ (0, 1] and ς ∈ (0, 1]. In our model, x(τ ) and y(τ ) represent prey and predator populations at time τ , respectively. The environmental carrying capacity of prey and intrinsic growth rate and population are denoted by k and r, respectively; c represents predators attack rate on prey population; e denotes recently conceived predators for each prey catch; d represents the predator mortality rate. This system involves a constant shelter covering Mx(τ ) of the prey such that m ∈ [0, 1), allowing (1 -M)x(τ ) of the prey be accessible to the predator [36]. These constants are all positive.
On the other hand, compared with constant FO systems, the research on variable FO differential equations is relatively new, and a numerical approximation of these equations is still at an early stage of development [37][38][39][40][41][42][43][44][45][46][47][48][49]. The variable FO operator of Atangana-Baleanu type in the sense of Liouville-Caputo is given as follows: where the normalized function ( (τ )) = 1 -(τ ) + (τ ) ( (τ )) . A variable FO integral in the sense of Atangana-Baleanu is listed as follows: In the second case, we consider that the predatory-prey species is not an instantaneous operation; for this case, the model given by Eq.
where the variable FO operator ABC 0 Dτ (τ ) denotes Atangana-Baleanu-Caputo for each operator, and time delay κ represents the gestation delay.

Predator-prey model involving constant-order fractional derivative
First case This portion is to investigate the existence and uniqueness of the solutions under constant order assuming the Atangana-Baleanu-Caputo FO operator.

Existence and uniqueness of the solution For every real-valued continuous function, let
, containing the subnorm and shaft ={x, y ∈ y, x(τ ) > 0 and y(τ Let be a real Banach space with a cone H. H initiates a restricted order in the succeeding approach [45] A cone speaks to typical if one can locate a positive constant α with the end goal that r, l, ∈ , ϒ < r < l → r ≤ α l , where ϒ is the zero component of . [45]. To study the existence of solutions, let us assume the system given by Eq. (2.1).

Theorem 1 Let a Banach space Z have a closed set subspace H and T : H → H with Lipschitz constant υ < 1 be a contraction mapping so that T carries a fixed point τ
Model (2.1) is identical to that of Volterra when it is integral in the Atangana-Baleanu sense. We apply the AB fractional integral in Eq. (2.1) to get the system Now we assume the following lemmas.

Lemma 1 Suppose that T : H → H is a mapping and gives
To handle this more easily, let us assume the following:

Lemma 2 Let M ⊂ H be bounded and constants
.
Due to Eqs. (3.6) and (3.7), we assume that T is bounded. Now observe x(τ ) ∈ M, τ 1 , τ 2 ; τ 1 < τ 2 and if |τ 2τ 1 | < μ for given φ > 0, hence Now we assume an integral part Now we will study the following: Equation (3.10) by taking the Lipschitz condition of the operator becomes  Taking the Lipschitz condition of the operator and Eq. (3.13), we consider the following: is equicontinuous and compact [45] due to the Arzela-Ascoli principle. This completes the proof.

. Then the concluding equation has a positive solution.
Proof We require a fixed point operator of T to research the existence of a positive solution. With Lemma 1 structure, assume that the operator T : H → H is entirely continuous. Consider two arbitrary predator population densities of x 1 and x 2 in H which meet x 1 ≤ x 2 as well as the densities of prey population S 1 and S 2 in H which satisfy S 1 ≤ S 2 , then S is a positive function and are satisfied, thus the mapping T is increasing, we have Tψ 2 ≥ψ 2 and Tψ 2 ≤ψ 2 . So T : ψ 1 ,ψ 2 → ψ 1 ,ψ 2 is a compact operator by the guideline of Lemma 2 and continuous by the perspective of Lemma 1. As H is a normal cone of T [45], this completes the proof.

Uniqueness of the special solution Consider
and Therefore, if the following conditions hold: then the system produces a unique positive solution because mapping T is a contraction with fixed point.

Numerical approximation The Atangana-Baleanu FO integral numerical approximation by imposing the Adams-Moulton rule is granted by [24]
AB where b r = (r + 1) 1--(r) 1-. With a similar previous numerical method, we have Now consider FO differential equations with variable-order ABC of model Eq. (2.7), in which, to justify temporal changes, a FO derivative as a variable is used instead of integerorder derivatives.

Predator-prey model involving variable-order fractional derivative
Adams-Moulton method with variable order Now, we show an alteration of variable order (τ ) of the Adams-Moulton approach. Along these lines, the algorithm computes numerical results of a FO differential equation of the kind ABC 0 Equation (3.23) has one solution categorized in τ ∈ [0, T], which can be reformulated by applying the FO Atangana-Baleanu integral as follows: (3.24) where ( (τ )) = 1 -(τ ) + (τ ) ( (τ )) represents a normalization function. The solution plot by applying the trapezoidal quadrature formulation appears as described: , j = 0, 1, 2, . . . , i.

Theorem 3 We show the existence and uniqueness of the solution.
Proof Assuming that the space we definē (3.32) Find g f (τ )y 0 < y max to establish the condition of well-posedness for which Taking the triangular inequality, we obtain For the above, we need to have . (3.35) Presently, use an operator of Picard to contract on T max under some condition.
Let f 1 and f 2 be two functions in C[A T max (τ 0 ), D y max (y 0 )], then let us calculate the following: The condition for contraction is Under the prior setting, the derived Picard's operation on a Banach space with the measurement included by norm uniformly is a contraction, hence has the feature that an exceptional function x exists means x = x, which is the special solution of Eq. (3.24).

Numerical simulations
In this section, we show the dynamics of the fractional-order operator on the proposed model in the sense of constant fractional-order and variable fractional-order operator (τ ). By assuming the parameter values and choosing a numerical method, we illustrate the mathematical results of the system by Eqs. (2.1), (2.7), and (2.8) with the numerical simulations (3.22), (3.26), and (4.5). We observed in the Figs. 2a-2f and 3a-3f that convergency of fractional-order derivative is faster than that of noninteger operator. We analyze  tions of ( -ς), ( (τ )-ς(τ )) or delay (κ), arbitrary selected. It is significant that expectation relies upon the estimation of alpha not just close to the hypothetical boundaries.

Conclusions
A model of predator-prey consisting of constant and variable-order time FO with and without delay (gestation delay) involving operators related to a Mittag-Leffler kernel of type Liouville-Caputo is investigated in this paper. Because of this system's nonlocality, the predator-prey model was manufactured utilizing another FO differentiation developed along with the generalized Mittag-Leffler function as a kernel. With finality to make our model more realistic, preys and predators did not follow the same fractional dynamics and the orders of FO derivatives of prey and predator population concerning time ∈ (0, 1] and ς ∈ (0, 1] were not equal. The fixed point theorem is beneficial to explain the existence of a novel model having a special set of solutions. We offered a numerical pattern dependent on the Adams-Bashforth-Moulton method for variable order. The existence and uniqueness of this numerical framework and complex dynamics of the structure were studied with modification of fractional orders (ς), ( (τ )ς(τ )) or delay (κ).
The presented findings emphasize that the variable-order fractional methodology was more appropriate to explain complex dynamics configurations under investigation. Eventually, the predator-prey model may show rich dynamical conduct. An Intel Core i7, 2.6 GHz CPU, 16.0-GB RAM (Matlab R.2013a) computer program was used to get the results in this article.