A Rosenzweig–MacArthur Model with Continuous Threshold Harvesting in Predator Involving Fractional Derivatives with Power Law and Mittag–Leffler Kernel

Hasan  S. Panigoro; Agus Suryanto; Wuryansari  Muharini Kusumawinahyu; Isnani Darti

doi:10.3390/axioms9040122

A Rosenzweig–MacArthur Model with Continuous Threshold Harvesting in Predator Involving Fractional Derivatives with Power Law and Mittag–Leffler Kernel

Panigoro, Hasan S.;Suryanto, Agus;Kusumawinahyu, Wuryansari Muharini;Darti, Isnani 2020-10-22 00:00:00 axioms Article A Rosenzweig–MacArthur Model with Continuous Threshold Harvesting in Predator Involving Fractional Derivatives with Power Law and Mittag–Lefﬂer Kernel 1,2 1, 1 Hasan S. Panigoro , Agus Suryanto * , Wuryansari Muharini Kusumawinahyu and Isnani Darti Department of Mathematics, Faculty of Mathematics and Natural Sciences, University of Brawijaya, Malang 65145, Indonesia; hspanigoro@ung.ac.id (H.S.P.); wmuharini@ub.ac.id (W.M.K.); isnanidarti@ub.ac.id (I.D.) Department of Mathematics, Faculty of Mathematics and Natural Sciences, State University of Gorontalo, Bone Bolango 96119, Indonesia * Correspondence: suryanto@ub.ac.id Received: 17 September 2020; Accepted: 19 October 2020; Published: 22 October 2020 Abstract: The harvesting management is developed to protect the biological resources from over-exploitation such as harvesting and trapping. In this article, we consider a predator–prey interaction that follows the fractional-order Rosenzweig–MacArthur model where the predator is harvested obeying a threshold harvesting policy (THP). The THP is applied to maintain the existence of the population in the prey–predator mechanism. We ﬁrst consider the Rosenzweig–MacArthur model using the Caputo fractional-order derivative (that is, the operator with the power-law kernel) and perform some dynamical analysis such as the existence and uniqueness, non-negativity, boundedness, local stability, global stability, and the existence of Hopf bifurcation. We then reconsider the same model involving the Atangana–Baleanu fractional derivative with the Mittag–Lefﬂer kernel in the Caputo sense (ABC). The existence and uniqueness of the solution of the model with ABC operator are established. We also explore the dynamics of the model with both fractional derivative operators numerically and conﬁrm the theoretical ﬁndings. In particular, it is shown that models with both Caputo operator and ABC operator undergo a Hopf bifurcation that can be controlled by the conversion rate of consumed prey into the predator birth rate or by the order of fractional derivative. However, the bifurcation point of the model with the Caputo operator is different from that of the model with the ABC operator. Keywords: Rosenzweig–MacArthur model; fractional derivatives; threshold harvesting 1. Introduction More than 50 years after the model has been proposed, the Rosenzweig–MacArthur predator–prey model [1] has been consistently developed by many scholars to approach the real world phenomena with more realistic mathematical models. The commonsensical modiﬁed Rosenzweig–MacArthur models are accomplishable by considering the biological perspectives of ecosystem conditions, for instance the stage structure [2,3], the refuge effect [4–8], the fear effect [9], the Allee effect [10,11], the intraspeciﬁc competition [12,13], the cannibalism [14], the infectious disease [15–17], and so forth. On the other hand, the modeling also contemplates the optimal management of bioeconomic resources as in ﬁshery and pest management. Some researchers put an intervention into the predator–prey model such as the harvesting to one or more population [8,18–22]. To protect Axioms 2020, 9, 122; doi:10.3390/axioms9040122 www.mdpi.com/journal/axioms Axioms 2020, 9, 122 2 of 23 the population from over-exploitation during the harvesting, some management techniques have been established. One of the famous technique is a continuous threshold harvesting policy (THP) (see [23–29]). THP is regulated as follows: when the population density above the threshold level, harvesting is permitted; when the population density drops below the threshold level, harvesting is prohibited. In 2013, Lv et al. [26] proposed the following Rosenzweig–MacArthur model with THP in predator dx x mxy =rx 1 , dt K a + x (1) dy nxy = dy H(y), dt a + x where < 0 , if y < T, H(y) = h(y T) : , if y T. c + (y T) We describe the biological interpretation of variables and parameters of model (1) in Table 1. Model (1) represents an interaction between two populations with a prey–predator relationship, where THP is only applied for the predator to preserving its populations. Some appealing examples of the ecological model (1) are given by the interaction between Sycanus sp. and Setothosea asigna and between Rhinocoris sp. and Spodoptera litura. Shepard [30] reported that Sycanus sp. and Rhinocoris sp. are the natural predators of the pests such as Setothosea asigna and Spodoptera litura which exist in agricultural lands and plantations. The worrying problem is: How if the density of insects such as Sycanus sp. and Rhinocoris sp. uncontrolled? One solution is applying THP as in model (1). Table 1. Description of variables and parameters of the model (1). Variables and Parameters Description x(t) The density of prey y(t) The density of predator r The intrinsic growth rate of prey K The environmental carrying capacity of prey m The maximum uptake rate for prey n The conversion rate of consumed prey into predator birth a The environment protection for prey d The natural death rate of predator h The harvesting rate c The half saturation constant for harvesting T The threshold level of harvesting Lv et.al. [26] successfully explored the dynamics of the model (1) including the local stability and the existence of various phenomena (saddle-node, Hopf, cusp, and Bogdanov–Taken bifurcations). Despite their success works, the model with the first-order derivative is limited to its capability of involving all previous conditions to the growth rates of both predator and prey. The growth rates of both populations in the model (1) depend only on the current state. Biologically, the growth rates must be dependent on all of the previous conditions which are known as the memory effects. To account for such memory effects, some researchers proposed to apply the fractional-order derivative instead of the first-order derivative when expressing the growth rate of the population. The fractional-order models are naturally related to systems with memory which exists in most biological models [7,31]. The fractional-order models are also well-liked due to their capability in providing an exact description of different nonlinear phenomena [32]. In recent years, the development of fractional-order models grows rapidly and becomes popular in studying the dynamical behavior of predator–prey interaction [17,33–38]. It has been shown that the order of the fractional derivate signiﬁcantly affects the dynamical behavior of Axioms 2020, 9, 122 3 of 23 the models, which is in contrast to the ﬁrst-order derivative models that depend only on the values of parameters. In this paper, we modify the model of Lv et.al. [26] by applying the fractional-order derivative to the left-hand sides of the first-order differential Equation (1). We use two types of fractional-order derivatives, namely the Caputo operator (that is, the operator with the power law kernel) [39] and Atangana–Baleanu operator [40]. The basic difference between these two operators lies on their kernel. Atangana–Baleanu operator has a non-singular and non-local (that is, Mittag–Leffler) kernel while the Caputo operator does not [41,42]. From the biological meanings, a model with Atangana–Baleanu operator gives better results in applying memory effects [43–45]. Nevertheless, the Caputo operator has more complex analytical tools in investigating the dynamics of the model such as the local stability [46–50], the global stability [51,52], bifurcation theory [52–54], and so on. By considering the deﬁciencies and advantages, the model with Caputo and Atangana–Baleanu operator are employed in our work. Based on our literature review, the dynamics of the model (1) with Caputo and Atangana–Baleanu operator have never been studied. For this reason, we are interested in investigating the dynamical behavior of model (1) using both Caputo and Atangana–Baleanu fractional-order operators. If the ﬁrst-order derivatives at the left hand sides of model (1) are replaced by the dt fractional-order derivatives D , then we obtain x m ˆ xy D x =r ˆx 1 , K a + x (2) n ˆ xy D y = dy H(y). a + x Note that the left hand sides of model (2) have the dimension of (time) , while the parameters at ˆ ˆ the right hand sides of model (2) such as r ˆ, m ˆ , n ˆ , d, and h have the dimension of (time) ; this shows the inconsistency of physical dimensions in the model (2) (see [55,56]). To overcome such inconsistency, we rescale the parameters in the model (2) to get the following model x m ˆ xy a a D x =r ˆ x 1 , K a + x (3) n ˆ xy a a D y = d y H(y), a + x where 0 , if y < T, H(y) = h (y T) , if y T. c + (y T) a a a a a ˆ ˆ ˆ ˆ ˆ By applying new parameters r = r , m = m , n = n , d = d , and h = h , we obtain x mxy D x =rx 1 , K a + x (4) xy D y = dy H(y), a + x where 0 , if y < T, H(y) = h(y T) : , if y T. c + (y T) This paper is organized as follows. In Section 2, we investigate the dynamics of model (4) with the Caputo operator. We identify the existence, uniqueness, non-negativity, and boundedness of solutions. Furthermore, we explore the dynamics of the model by examining the existence of the equilibrium points, their local and global stability, and the existence of Hopf bifurcation. In Section 3, the existence and uniqueness of solutions of the model with Atangana–Baleanu operator are veriﬁed. In Section 4, Axioms 2020, 9, 122 4 of 23 we explore the dynamics of the model using both operators numerically. We demonstrate numerically the stability of the equilibrium point, and the occurrence of forward and Hopf bifurcations. We end our works with a brief conclusion in Section 5. 2. The Caputo Fractional-Order Rosenzweig–MacArthur Model with THP in Predator 2.1. Model Formulation The operator of Caputo fractional-order derivative is deﬁned as follows Deﬁnition 1. [48] Let a 2 (0, 1], f 2 C ([0, +¥),R), and G() is the Gamma function. The Caputo fractional derivative of order-a is deﬁned by C a a 0 D f (t) = (t s) f (s)ds, t 0. (5) G(1 a) The kernel of Caputo operator is known as the power law kernel. By applying Deﬁnition 1 to model (4), we get the Caputo fractional order Rosenzweig–MacArthur model with THP in predator x mxy C a D x =rx 1 F , K a + x (6) nxy C a D y = dy H(y) F . a + x 2.2. Existence and Uniqueness In this part, we study the existence and uniqueness of model (6). Lemma 1. [57] Consider a Caputo fractional-order system C a D x(t) = f (t, x(t)), t > 0, x(0) 0, a 2 (0, 1], (7) n n where f : (0, ¥) Q ! R , Q R . A unique solution of Equation (7) on (0, ¥) Q exists if f (t, x(t)) satisﬁes the locally Lipschitz condition with respect to x. Since the right hand-side of model (6) is a piecewise function which is switched when the number of predators passes through the threshold level, we divide the existence and uniqueness of the solution into two cases, namely y T and y < T. We start from y T. Consider the region Q [0, T ] where Q := (x, y) 2 R : max (jxj,jyj) g, y T , T < +¥, and a mapping F(L) = (F (L), F (L)). + 1 2 For any L = (x y) 2 Q and L = (x ¯, y ¯) 2 Q, we obtain , Axioms 2020, 9, 122 5 of 23 ¯ ¯ ¯ F(L) F(L) = F (L) F (L) + F (L) F (L) 1 1 2 2 x mxy x ¯ mx ¯y ¯ = rx 1 rx ¯ 1 + K a + x K a + x nxy nx ¯y ¯ ¯ ¯ dy H(y) dy H(y) a + x a + x ¯ r xy x ¯y ¯ 2 2 ¯ x x ¯ m = r(x x) + K a + x a + x ¯ xy x ¯y ¯ n d (y y ¯) ( H(y) H(y ¯)) a + x a + x ¯ ¯ ¯ r xy xy =rjx x ¯j + jx + x ¯jjx x ¯j + (m + n) + K a + x a + x ¯ h(y T) h(y T) djy y ¯j + c + (y T) c + (y ¯ T) ¯ ¯ ¯ ¯ r ay(x x) + (ax + xx)(y y) =r x x ¯ + x + x ¯ x x ¯ + (m + n) + j j j jj j K (a + x)(a + x ¯) ch(y y) d y y ¯ + j j (c + y T)(c + y ¯ T) 2gr m + n r x x ¯ + x x ¯ + ay(x x ¯) + (ax ¯ + x ¯ x)(y y ¯) + j j j j j j K a ¯ ¯ djy yj + jy yj 2gr (m + n)g rjx x ¯j + jx x ¯j + jx x ¯j + K a (m + n)(ag + g ) h y y ¯ + d y y ¯ + y y ¯ j j j j j j a c 2gr (m + n)g (m + n)(ag + g ) h = r + + jx x ¯j + + d + jy y ¯j K a a c M L L , 2gr (m + n)g (m + n)(ag + g ) h where M = max r + + , + d + . Hence, F(L) satisﬁes the K a a c Lipschitz condition for y T. By using similar approaches, when y < T, it is easy to check that 2gr (m + n)g (m + n)(ag + g ) ¯ ¯ F(L) F(L) M L L , where M = max r + + , + d 2 2 K a a and hence the Lipschitz condition is also satisﬁed. According to Lemma 1, model (6) with non-negative initial condition has a unique solution L(t) = (x(t), y(t)) 2 Q. Thus, we establish the following theorem. Theorem 1. For each non-negative initial condition (x , y ) 2 Q, there exists a unique solution 0 0 (x(t), y(t)) 2 Q of model (6), which is deﬁned for all t 0. 2.3. Non-Negativity and Boundedness The solution of model (6) is required to be nonnegative and bounded to establish a biologically well-behaved model. To determine the non-negativity and boundedness of the solution of model (6), the following lemmas are needed. C a Lemma 2. [58] Let 0 < a 1. Suppose that f (t) 2 C[a, b] and D f (t) 2 C[a, b]. C a If D f (t) 0,8 t 2 (a, b), then f (t) is a non-decreasing function for each t 2 [a, b]. C a If D f (t) 0,8t 2 (a, b), then f (t) is a non-increasing function for each t 2 [a, b]. t Axioms 2020, 9, 122 6 of 23 Lemma 3. (Standard comparison theorem for Caputo fractional-order derivative [31]). Let x(t) 2 C [0, +¥) . ( ) C a 2 If x(t) satisﬁes D x(t) + lx(t) m, x(0) = x , where a 2 (0, 1], (l, m) 2 R and l 6= 0, then x(t) m m x E [lt ] + . 0 a l l In the following theorem, we prove the non-negativity and boundedness of solutions using the above lemmas. Theorem 2. All solutions of model (6) with non-negative initial conditions are non-negative and uniformly bounded. Proof. We start by proving that, if the initial condition is non-negative, then x(t) 0 for all t ! ¥. Suppose that it is incorrect; then, we can ﬁnd t > 0 such that x(t) > 0, 0 t < t , < 1 x(t ) = 0, (8) : + x(t ) < 0 By employing (8) and the ﬁrst equation of model (6), we obtain C a D x(t ) = 0. (9) t 1 x(t )=0 + + Based on Lemma 2, we have x(t ) = 0, which contradicts the fact that x(t ) < 0. Thus, x(t) 0 1 1 for all t 0. Using a similar procedure, we conclude y(t) 0 for all t > 0. Now, we adopt the similar manner as in [34] to show the boundedness of solutions. By setting up my a function V(t) = x + , we get m dmy C a C a C a D V(t) + dV(t) = D x + D y + dx + t t t n n x mxy m nxy dmy =rx 1 + dy H(y) + dx + K a + x n a + x n rx m H(y) =(d + r)x K n r (d + r)K (d + r) K m H(y) = x + K 2r 4r n (d + r) K 4r According to the standard comparison theorem for Caputo fractional-order derivative in Lemma 3, we achieve the following inequality 2 2 (d + r) K (d + r) K V(t) V(0) E [d(t) ] + , 4r 4r where E is the one-parameter Mittag–Lefﬂer function. Since E [d(t) ] ! 0 as t ! 0, we acquire a a (d + r) K V(t) ! for t ! ¥. Hence, with non-negative initial condition, all solutions are restricted to 4r the region Q where my (d + r) K Q := (x, y) 2 R : x + + #, # > 0 . n 4r Consequently, all solutions of model (6) are uniformly bounded. Axioms 2020, 9, 122 7 of 23 2.4. The Existence of Equilibrium Point The ﬁrst commonplace technique in studying the dynamical behavior of a fractional-order system is investigating the existence of the equilibrium point. Due to the biological nature, we give the following deﬁnition. Deﬁnition 2. Consider a Caputo fractional-order system C a D ~ x = f (~ x);~ x(0) = ~ x ; a 2 (0, 1]. (10) A point ~ x is called an equilibrium point of Equation (10) if it satisﬁes f (~ x ) = 0. Particularly, it is called the biological equilibrium point of Equation (10) if it satisﬁes ~ x 0. Based on Deﬁnition 2, the equilibrium point of model (6) is obtained by solving rx my r x =0, K a + x (11) nxy dy H(y) =0. a + x Thus, we get four feasible biological equilibrium points as follows. (i) The equilibrium points when y < T are (i.a) the origin point E = (0, 0) which always exists; (i.b) the predator extinction point E = (K, 0) which always exists; and (K x ˆ) (a + x ˆ)r ad (i.c) the ﬁrst co-existence point E = x ˆ, , with x ˆ = , which exists if mK n d ad TmK n > + d and (K x ˆ) (a + x ˆ) < . K r (ii) The equilibrium point when y T is the second co-existence point E = (x , y ) where y = (K x ) (a + x )r 4 3 2 and x is the positive roots of polynomial b x + b x + b x + b x + b = 0 1 2 3 4 5 mK where b =(n d)r , b = [(an + 2dK) 2(ad + nK)] r , b =(nrK + 4adr + cdm + mnT + hm)rK ((drK + 2anr + cmn + dmT)K + a dr)r, b =((anr + cmn + dmT)K + (2adr + hm + cdm)a)rK ((2adr + hm + cdm + mnT)K + admT)rK, b = (adr + hm)mT (adr + hm + cdm)ar K . [ ] TmK E exists if 0 < x < K and (K x ) (a + x ) . 2.5. Local Asymptotic Stability In this part, we discuss the local stability of E , E , E, and E . For this aim, we need the 0 1 following theorem. Theorem 3. (Matignon condition [48,59]) The equilibrium point ~ x of system (10) is locally asymptotically stable if all eigenvalues l of the Jacobian matrix J = ¶ f /¶~ x evaluated at ~ x satisfy j arg(l )j > ap/2. j j If there exists at least one eigenvalue satisfying j arg(l )j > ap/2 and j arg(l )j < ap/2, k 6= l, then ~ x is k l a saddle-point. Axioms 2020, 9, 122 8 of 23 Now, we present Theorems 4–7 to show the local stability properties of E , E , E, and E . 0 1 Theorem 4. The origin point E = (0, 0) is always a saddle point. " # r 0 Proof. When E = (0, 0), the model (6) has Jacobian matrix J(E ) = , where its eigenvalues 0 0 0 d are l = r > 0 and l = d < 0. It is clear that jarg(l )j = 0 < ap/2 and jarg(l )j = p > ap/2. 1 2 1 2 Therefore, based on Theorem 3, E is always a saddle point. ad Theorem 5. The predator extinction point E = (K, 0) is locally asymptotically stable if n < + d. ad Otherwise, if n > + d, then E = (K, 0) is a saddle point. " # mK a+K Proof. The Jacobian matrix of model (6) evaluated at E is J(E ) = . The eigenvalues 1 1 nK 0 d a+K nK of J(E ) are l = r < 0 and l = d. Clearly, jarg(l )j = p > ap/2 and jarg(l )j = p > 1 1 2 1 2 a + K ad ad ap/2 if n < + d and arg(l ) = 0 < ap/2 if n > + d. Hence, we have the theorem. j j K K Remark 1. It is noted that the existence condition for the ﬁrst co-existence point E contradicts the stability condition of E . Consequently, if E is locally asymptotically stable, then E does not exist. This condition also 1 1 indicates the existence of forward bifurcation, which is conﬁrmed numerically in the next section. 2 2 ˆ ˆ ˆ ˆ ˆ 4 (K x) anrx K x r x 2 D(a + x)K Theorem 6. Let D = 1 and a ˆ = tan . 2 2 (a + x ˆ) K a + x ˆ K p (K a 2x ˆ) rx ˆ Suppose that one of the following statements is obeyed. K a (i) x ˆ > ; or K a (ii) x ˆ < , D > 0, and a < a ˆ . ˆ ˆ (K x) (a + x)r Then, the ﬁrst co-existence point E = x ˆ, is locally asymptotically stable. mK Proof. We ﬁrst observe that the Jacobian matrix of model (6) evaluated at E is 2 3 ˆ ˆ ˆ K x rx mx 6 7 a + x ˆ K a + x ˆ ˆ 6 7 J(E) = . (12) 4 5 (K x ˆ) anr (a + x)mK The eigenvalues of the Jacobian matrix (12) are the solutions of the characteristic equation K x ˆ rx ˆ K x ˆ anrx ˆ ( ) l 1 l + = 0, a + x ˆ K (a + x ˆ) K namely K x ˆ rx ˆ i D l = 1 + , a + x ˆ 2K 2 (13) ˆ ˆ K x rx i D l = 1 . a + x ˆ 2K 2 Axioms 2020, 9, 122 9 of 23 K a When x ˆ > , the real parts of l are negative. Thus, the eigenvalues (13) always satisfy 1,2 K a (K x ˆ) anrx ˆ arg(l ) > ap/2 for any D. If x ˆ < and D 0, then l l = > 0 and l + j j 1,2 1 2 1 2 (a + x ˆ) K K x ˆ rx ˆ K a l = 1 > 0, meaning that arg(l ) = 0 < ap/2. If x ˆ < and D > 0, j j 2 1,2 a + x ˆ K 2 then jarg(l )j > ap/2 for a < a ˆ . Hence, we prove the theorem. 1,2 Theorem 7. Let ((y T) cT)h my r x = + x , 2 2 y (c + y T) (a + x ) K my r ((y T) cT)h x mnx y q = x + 1 . 2 2 2 (a + x ) K y (c + y T) a + x (a + x ) If one of the following statements is satisﬁed,then the second co-existence point E = (x , y ) is locally asymptotically stable: (i) q > 0 and x < 0; or (ii) x < 4q, x > 0, and a < a . Proof. The Jacobian matrix of model (6) evaluated at E is given by 2 3 my r mx 6 7 K a + x (a + x ) 6 7 J(E ) = . (14) 6 7 4 5 x ny ((y T) cT)h a + x a + x y (c + y T) The eigenvalues of (14) are obtained by solving the characteristic equation l xl + q = 0. Thus, x x 4q we have l = . If q > 0 and x < 0, then jarg(l )j > ap/2. If x < 4q and x > 0, 1,2 1,2 2 2 thenjarg(l )j > ap/2 for a < a . Using Theorem 3, the local stability of E is completely proven. 1,2 2.6. Global Asymptotic stability To study the global stability of equilibrium points, we need the following lemmas. Lemma 4. [51] Let x(t) 2 C (R ), x 2 R , and its Caputo fractional derivative of order-a exists for any + + x(t) x a C a C a 2 (0, 1]. Then, for any t > 0, we have D x(t) x x ln 1 D x(t). t t x x(t) Lemma 5. (Generalized Lasalle Invariance Principle [52]). Suppose W is a bounded closed set and every solution of system C a D x(t) = f (x(t)), (15) which starts from a point in W remains in W for all time. Let V(x) : W ! R be a continuously differentiable function such that C a D Vj 0. (15) C a Let E := xj D Vj = 0 and M be the largest invariant set of E. Then, every solution x(t) originating (15) in W tends to M as t ! ¥. Since model (6) is basically a piecewise fractional-order differential equations that depends on T, the analysis of the global stability is split into two regions deﬁned by W := f(x, y) : x 0, y < Tg and W := (x, y) : x 0, y T . Therefore, the global stabilities of E , E, and E are investigated f g 2 1 as follows. Axioms 2020, 9, 122 10 of 23 ad Theorem 8. If n < , then the predator extinction point E = (K, 0) is globally asymptotically stable in the region W . x my Proof. By specifying a positive Lyapunov function L (x, y) = x K K ln + , K n and conforming to Lemma 4, we obtain x K m C a C a C a D L (x, y) D x + D y t t t x n rx my m nxy =(x K) r + dy K a + x n a + x rx mKy dmy =2rx rK + K a + x n r mKy dmy = (x K) + K a + x n r mKy dmy (x K) + K a n d K my. n a ad C a Owing to the fact that n < , we have D L (x, y) 0. In consequence of Lemma 5, E is 1 1 globally asymptotically stable in the region W . ad Remark 2. Notice E is locally asymptotically stable if n < + d and is globally asymptotically stable if ad n < . Hence, if the global stability condition is fulﬁlled, then the local stability is also achieved but not vice versa. This fact reinforces that the global stability condition has a larger attracting region than that of the local stability condition. ˆ ˆ (K x) (a + x)r Theorem 9. If (K x ˆ) (a + x ˆ) < a , then the ﬁrst co-existence point E = x ˆ, is globally mK asymptotically stable in the region W . K x ˆ (a + x ˆ)r ( ) Proof. Let E = (x ˆ, y ˆ) where y ˆ = and j is a positive constant. By considering a mK h i x y ˆ ˆ ˆ ˆ Lyapunov function L (x, y) = x x x ln + j y y y ln , and using Lemma 4, we get x ˆ y ˆ x x ˆ y y ˆ C a C a C a D L (x, y) D x + j D y t t t x y rx my nx =(x x ˆ) r + j(y y ˆ) d K a + x a + x rx ˆ my ˆ rx my nx nx ˆ =(x x ˆ) + + j(y y ˆ) K a + x ˆ K a + x a + x a + x ˆ r y y ˆ x x ˆ = (x x ˆ) (x x ˆ) + m + nj(y y ˆ) K a + x a + x ˆ a + x a + x ˆ r (y y ˆ)(a + x ˆ) + (x ˆ x)y ˆ (x x ˆ)(y y ˆ) = (x x ˆ) (x x ˆ) m + anj K (a + x)(a + x ˆ) (a + x)(a + x ˆ) ˆ ˆ ˆ ˆ ˆ r (x x)(y y)(a + x) (x x) y = (x x ˆ) m + m+ K (a + x)(a + x ˆ) (a + x)(a + x ˆ) (x x ˆ)(y y ˆ) anj (a + x)(a + x ˆ) r my ˆ (x x ˆ)(y y ˆ) (x x ˆ) + (anj (a + x ˆ)m). K a (a + x)(a + x ˆ) Axioms 2020, 9, 122 11 of 23 (a + x ˆ)m By choosing j = and substituting the value of y ˆ, we get an C a 2 D L (x, y) (x x ˆ) a (K x ˆ) (a + x ˆ) . t 2 a K 2 C a Therefore, if (K x ˆ) (a + x ˆ) < a , then D L (x, y) 0. Using Lemma 5, we conclude that E is globally asymptotically stable in the region W . Based on Theorem 2, x(t) and y(t) are bounded. Let Y be the upper bound of y(t) such that 0 < T y(t) Y. The global stability of E is stated in the following theorem. a r cT Theorem 10. If y < min , + T , then the second co-existence point E = (x , y ) is globally mK Y T asymptotically stable in the region W . h i Proof. We consider a positive Lyapunov function L (E ) = x x x ln + j y y y ln . According to Lemma 4, the fractional derivative of L (E ) satisﬁes x x y y C a C a C a D L (x, y) D x + j D y t t t x y rx my nx h(y T) =(x x ) r + j (y y ) d K a + x a + x (c + y T)y rx my rx my =(x x ) + K a + x K a + x nx nx h(y T) h(y T) + j (y y ) + a + x a + x (c + y T)y (c + y T)y r (y y )(a + x )m (x x )my = (x x ) (x x ) K (a + x )(a + x) (x x )an (y y )chT (y y )(y T)(y T)h + j (y y ) (a + x )(a + x) (c + y T)(c + y T)y y r (x x )(y y )(a + x )m (x x ) my = (x x ) + K (a + x )(a + x) (a + x )(a + x) (x x )(y y )anj (cT (y T)(y T))(y y ) hj (a + x )(a + x) (c + y T)(c + y T)y y r (x x )(y y )(a + x )m my 2 2 (x x ) + (x x ) K (a + x )(a + x) a (x x )(y y )anj (cT (y T)(y T))(y y ) hj (a + x )(a + x) (c + y T)(c + y T)y y (a + x )m By taking j = and remembering that y(t) < Y, we obtain an r my (cT (y T)(Y T))(y y ) (a + x )mh C a 2 D L (E ) (x x ) . K a (c + y T)(c + y T)any y a r cT C a It is easily conﬁrmed that, if y < min , + T , then D L (x, y) 0. Based on mK Y T Lemma 5, E is globally asymptotically stable in the region W . 2 Axioms 2020, 9, 122 12 of 23 2.7. The Existence of Hopf Bifurcation The Hopf bifurcation is a local phenomenon when a stable equilibrium point loses its stability and all nearby solutions converge to a periodic solution namely limit-cycle if a parameter is varied [54,60,61]. It is shown that many fractional-order models involving the Caputo operator undergo a Hopf bifurcation which is driven by the order of the derivative (see [2,17,34,53,62]). The difference between the Hopf bifurcation in the integer-order model and that in the fractional-order model lies on the property of the limit-cycle. In the integer-order model, the limit-cycle is a periodic orbit which does not exist in the fractional-order model [63]. In the fractional-order model, the limit-cycle is not a periodic solution, but all nearby solutions converge to a limit-cycle [56,62]. Adapted from Theorem 3 in [62], for two dimensional Caputo fractional-order system, a Hopf bifurcation occurs when the eigenvalues l of the Jacobian matrix evaluated at the equilibrium point 1,2 satisfy the following conditions: (i) l = y wi where y > 0; 1,2 (ii) m(a ) = a p/2 min jarg(l )j = 0; and 1i2 i dm(a) (iii) 6= 0. da a=a Now, consider the stability condition in Theorems 6 and 7. For y < T, the Jacobian matrix of model (6) evaluated at E has a pair of complex eigenvalues if D > 0. We can easily conﬁrm that, K a if x ˆ < , then the real part of the eigenvalues are positive. We also have that m(a ˆ) = 0 and dm(a) 6= 0. Therefore, E undergoes a Hopf bifurcation when a crosses a ˆ . A similar circumstance da a=a ˆ also occurs when y T. When x < 4q, the Jacobian matrix of model (6) has a pair of complex eigenvalues. The real part of the eigenvalues are also positive when x > 0. We can also check that dm(a) m(a ) = 0 and 6= 0. This means the Hopf bifurcation also occurs when y T. Therefore, da a=a we have the following theorem. K a Theorem 11. (i) Let D > 0 and x ˆ < . The ﬁrst co-existence point E undergoes a Hopf bifurcation when a passes through a ˆ in the region W . (ii) Let x < 4q and x > 0. The second co-existence point E undergoes a Hopf bifurcation when a passes through a in the region W . 3. The Atangana–Baleanu Fractional-Order Rosenzweig–MacArthur Model with THP in Predator 3.1. Model Formulation In this section, we consider a fractional-order Rosenzweig–MacArthur model with THP in predator involving the Atangana–Baleanu operator. Speciﬁcally, we consider the Atangana–Baleanu operator in Caputo sense of order-a which is deﬁned as follows. Deﬁnition 3. [40] Suppose 0 < a 1. The Atangana–Baleanu fractional integral and derivative in Caputo sense of order-a (ABC) is deﬁned by 1 a a ABC a a1 I f (t) = f (t) + (t s) f (s)ds, B(a) G(a)B(a) B(a) a ABC a a 0 D f (t) = E (t s) f (s)ds, 1 a 1 a n ¥ where t 0, f 2 C ([0, +¥),R), E (t) = å is the Mittag–Lefﬂer function and B(a) is a k=0 G(ak + 1) normalization function with B(0) = B(1) = 1. In this paper, we deﬁne B(a) = 1 a + . G(a) Axioms 2020, 9, 122 13 of 23 By applying Deﬁnition 3 to model (4), we get the following fractional-order model with ABC operator x mxy ABC a D x =rx 1 G , K a + x (16) nxy ABC a D y = dy H(y) G . t 2 a + x Due to the lack of analytical theory, model (16) is investigated numerically. However, we ﬁrst show the existence and uniqueness of the solution of the model (16). 3.2. Existence and Uniqueness We start this work by representing the Lipschitz condition of the kernels of model (16). Since the harvesting is performed by obeying threshold harvesting policy, we give the proof into two cases i.e., y < T and y T. We start for y T. Let x, x ¯, y, and y ¯ be functions satisfying kxk a , kx ¯k a , kyk b , 1 1 and ky ¯k b . Suppose that r my g =r + (a + a ) + , 1 1 2 K a g =n + d + . Therefore, we get x mxy x ¯ mx ¯y kG (x) G (x ¯)k = rx 1 rx ¯ 1 1 1 K a + x K a + x 2 2 rx mxy rx ¯ mx ¯y = rx rx ¯ + + K a + x K a + x ¯ r x x ¯ 2 2 ¯ ¯ = r(x x) (x x ) my K a + x a + x ¯ (17) r amy(x x ¯) ¯ ¯ ¯ = r(x x) (x + x)(x x) K (a + x)(a + x ¯) r my ¯ ¯ ¯ rkx xk + (a + a ) kx xk + kx xk 1 2 K a r my = r + (a + a ) + kx x ¯k 1 2 K a =g x x ¯ , k k and nxy h(y T) nxy ¯ h(y ¯ T) kG (y) G (y ¯)k = dy dy ¯ 2 2 a + x c + (y T) a + x c + (y T) nxy h(y T) nxy ¯ h(y ¯ T) = dy + dy ¯ + a + x c + (y T) a + x c + (y ¯ T) nx(y y ¯) ch(y y ¯) = d(y y ¯) a + x (c + y T)(c + y ¯ T) (18) ¯ ¯ ¯ nky yk + dky yk + ky yk = n + d + ky y ¯k =g ky y ¯k . 2 Axioms 2020, 9, 122 14 of 23 When y < T, by utilizing the similar manner, we achieve G (y) G (y ¯) g y y ¯ where k k k k 2 2 3 g = n + d. Accordingly, the kernel of (16) satisﬁes the Lipschitz condition. Furthermore, if 0 g < 1 3 1 and 0 g < g < 1, then G and G are contracted. 3 2 1 2 Now, by employing the ﬁxed-point theorem, the solution of model (16) is investigated. By utilizing the Atangana–Baleanu fractional integral operator, model (16) is transformed into the following Volterra-type integral equation. 1 a a a1 x(t) x(0) = G (t, x) + (t s) G (s, x)ds, 1 1 B(a) B(a)G(a) (19) 1 a a a1 y(t) y(0) = G (t, y) + (t s) G (s, y)ds. 2 2 B(a) B(a)G(a) In a recursive form, Equation (19) is written as 1 a a a1 x (t) = G (t, x ) + (t s) G (s, x )ds, 1 n1 1 n1 B(a) B(a)G(a) (20) 1 a a a1 y (t) = G (t, y ) + (t s) G (s, y )ds. n 2 n1 2 n1 B(a) B(a)G(a) The associated initial conditions along with Equation (20) are x (t) = x(0) and y (t) = y(0). 0 0 By considering Equation (20), we have the difference expression of successive terms as follows. F (t) =x (t) x (t) 1,n n1 1 a a a1 = (G (t, x ) G (t, x )) + (t s) (G (t, x ) G (t, x )) ds, 1 n1 1 n2 1 n1 1 n2 B(a) B(a)G(a) (21) F (t) =y (t) y (t) 2,n n n1 1 a a a1 = (G (t, y ) G (t, y )) + (t s) (G (t, y ) G (t, y )) ds. 2 n1 2 n2 2 n1 2 n2 B(a) B(a)G(a) 0 n n According to Equation (21), we get x (t) = F (t) and y (t) = F (t). By applying n å n å 1,i 2,i i=1 i=1 Equations (17), (18) and (21), we have that 1 a a a1 kF (t)k g kF k + g kF (s)k (t s) ds, 1,n 1 1,n1 1 1,n1 B(a) B(a)G(a) (22) 1 a a a1 kF (t)k g kF k + g kF (s)k (t s) ds. 2,n 2 2,n1 2 2,n1 B(a) B(a)G(a) Therefore, by using (22), the existence and uniqueness of model (16) is presented as follows. Theorem 12. Model (16) has a unique solution if we can ﬁnd t such that t g (1 a)g + < 1, i = 1, 2, 3. (23) B(a) B(a)G(a) Proof. Let x(t) and y(t) be bounded functions, and hence the Lipschitz condition is satisﬁed. Thus, according to Equation (22), we obtain the following inequalities. (1 a)g t g 1 1 kF (t)k kx k + , 1,n B(a) B(a)G(a) (24) (1 a)g t g 2 2 F (t) y + . k k k k 2,n 0 B(a) B(a)G(a) Axioms 2020, 9, 122 15 of 23 Therefore, the continuity and existence of solution are proved since F (t) ! 0 and k k 1,n kF (t)k ! 0 as n ! ¥ and t = t . Now, suppose that 2,n 0 x(t) x(0) =x (t) U (t), n 1,n (25) y(t) y(0) =y (t) U (t). n 2,n We conﬁrm that n+1 (1 a) t n+1 kU (t)k + g (26) 1,n B(a) B(a)G(a) It is clear that kU (t)k ! 0 when n ! ¥. By using a similar manner, we acquire 1,n kU (t)k ! 0, n ! ¥. Finally, the uniqueness of solution for the model is proven. Suppose that there 2,n exists different set of solutions denote by x ˜(t) and y ˜(t); then, we have 1 a a a1 x(t) x ˜(t) = (G (t, x) G (t, x ˜)) + (G (s, x) G (s, x ˜))(t s) ds. (27) 1 1 1 1 B(a) B(a)G(a) Taking the norm for both sides and using a simpliﬁcation as in (22) and (24), we obtain (1 a)g t g 1 1 kx(t) x ˜(t)k 1 0. (28) B(a) B(a)G(a) For t = t , we have (23), thus kx(t) x ˜(t)k = 0 and hence x(t) = x ˜(t). Applying the same algebraic procedures, we can show that y(t) = y ˜(t). Therefore, the solution is unique. 4. Numerical Simulations In this section, the numerical simulations of Caputo model (6) and ABC model (16) are performed to illustrate the previous theoretical results. In the literature, there exist many numerical methods to solve a system of fractional differential equations, such as the Monte Carlo method [64], the Grünwald–Letnikov method [65,66] and the predictor–corrector method [67–69]. We apply the predictor–corrector scheme proposed by Diethelm et al. [67] to solve the Caputo fractional-order model and the predictor–corrector scheme proposed by Baleanu et al. [69] to solve the Atangana–Baleanu in Caputo sense model (ABC). Due to the limitation of ﬁeld data, we use hypothetical parameter values that correspond to the theoretical results. The initial parameter values are given as follows. r = 0.5, K = 1, m = 0.3, a = 0.2, d = 0.1, h = 0.1, T = 0.9, c = 0.1, anda = 0.9. (29) In Figure 1, we plot a bifurcation diagram by varying the conversion rate of consumed prey into predator birth n in interval [0.08, 0.2]. We notice that the parameter values (29) and the interval 0.08 n 0.2 lead to the non-existence of equilibrium point in W . Therefore, the ﬁrst numerical simulations are focused on the dynamics in W . For 0.08 n < n = 0.12, Theorem 5 states that the predator extinction point E = (1, 0) is the only equilibrium point which is asymptotically stable. To see this behavior, we take n = 0.1 and plot the phase portrait and the time series as in Figure 2. It is seen that all solutions with initial values in both W and W are convergent to E . When the initial value is in W , then the solution is oscillating 1 2 1 2 when it crosses the threshold harvesting level and eventually converges to E . When n passes through n , the equilibrium point E = (1, 0) undergoes a forward bifurcation. In this case, there appear two equilibrium points, namely the unstable E and an asymptotically stable ˆ ˆ E. Figure 1 shows that E is asymptotically stable if 0.12 < n . n = 0.1557. In Figure 3, we show the phase portrait and time series for the case of n = 0.14. We see that E = (0, 0) and E = (1, 0) 0 1 are saddle points, while E = (0.5, 0.5833) is asymptotically stable. This circumstance corresponds to Theorems 4–6 and 9. Axioms 2020, 9, 122 16 of 23 1.0 Caputo Upper bound of limit-cycle 0.8 ABC n = 0.14 E-stable E-unstable 0.6 n = 0.2 0.4 Lower bound of limit-cycle 0.2 E1-stable E1-unstable n = 0.1 0.0 1.2 E -stable E -unstable 1 1 1.0 n = 0.1 Upper bound of limit-cycle 0.8 0.6 0.4 n = 0.14 E-stable E-unstable 0.2 Caputo n = 0.2 Lower bound of limit-cycle ABC 0.0 0.08 0.10 0.12 0.14 0.16 0.18 0.20 Figure 1. Bifurcation diagram of Caputo model (6) and ABC model (16) driven by the conversion rate of consumed prey into predator birth (n) with constant parameter values (29). There exists two bifurcations namely a forward bifurcation which occurs when n passes through n 0.12, and a Hopf bifurcation which occurs when n passes through n 0.1557. 1.2 1.0 1.00 Caputo Caputo 0.95 0.8 0.90 ABC ABC 0.85 THP 1.0 0.80 2 3 4 5 6 7 8 0.6 0.4 0.8 0.2 0.0 0.6 1.2 1.0 0.4 0.8 0.2 0.6 Caputo 0.4 ABC 0.0 0.2 0.0 0.5 1.0 1.5 2.0 0 250 500 750 1000 1250 1500 1750 2000 x(t) t (a) (b) Figure 2. Numerical simulations of Caputo model (6) and ABC model (16) with parameter values (29) and n = 0.1: Figure (a) phase portrait; and Figure (b) time series. Furthermore, if we increase the value of n such that n > n , then the equilibrium point E loses its stability, and all solutions converge to a limit-cycle. This situation conﬁrms the occurrence of a Hopf bifurcation driven by n. For example, if we select n = 0.2 then all equilibrium points are unstable and all solutions eventually converge to limit cycle (see Figure 4). Now, we perform simulation using the same parameter values as in Figure 4, but with a lower threshold value, namely T = 0.5. In this case, there is no equilibrium point E in W , and E = (0.5954, 0.5364) occurs in W . According to Theorem 7, E is asymptotically stable. Such dynamics can be clearly seen in Figure 5. This simulation shows that by applying the THP when the interior equilibrium point is stable, we can choose a suitable constant of threshold so that the existence of both prey and predator are maintained. y(t) x(t) y(t) x(t) y(t) y(t) Axioms 2020, 9, 122 17 of 23 1.2 Caputo Caputo 0.8 ABC ABC THP 0.7 1.0 0.6 0.8 0.5 0.4 0.6 Caputo 1.0 ABC 0.4 0.8 0.6 0.2 0.4 0.2 0.0 E 0.00 0.25 0.50 0.75 1.00 1.25 1.50 0 100 200 300 400 500 x(t) t (a) (b) Figure 3. Numerical simulations of Caputo model (6) and ABC model (16) with parameter values (29) and n = 0.14: Figure (a) phase portrait; and Figure (b) time series. Notice that, in Figures 2–5, we see that both model with Caputo operator (6) and Atangana–Baleanu operator (16) have similar dynamical behavior. The noticeable difference between them is the orbit of solutions and the diameter of the limit-cycle. In Figure 4, the diameter of the limit-cycle of the model with ABC operator looks shorter than that of the Caputo operator, which may give different dynamics when a Hopf bifurcation occurs. To get more detail view, we plot a bifurcation diagram by varying the order of the fractional derivative (a) (see Figure 6). In this simulation, we use parameter values as in Figure 4 and vary the order-a in the interval [0.6, 0.9]. It is clearly seen that, besides the diameter of the limit-cycle, the bifurcation points of Caputo model and ABC model are also different. The Caputo model has an earlier bifurcation point than that of the ABC model. To show this situation, we show some numerical simulations using different values of a (see Figure 7). For a = 0.7, the equilibrium point E of both model are asymptotically stable. For a = 0.772, the equilibrium point E of the Caputo model loses its stability, while that of the ABC model remains asymtotically stable. For a = 0.83, the equilibrium point E of both models are unstable. 1.2 1.0 Caputo LC-ABC Caputo 0.8 ABC THP ABC LC-Caputo 1.0 0.6 0.4 0.8 0.2 0.0 0.6 1.0 Caputo 0.8 0.4 ABC 0.6 0.2 0.4 E 0.2 0.0 0.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0 250 500 750 1000 1250 1500 1750 2000 x(t) t (a) (b) Figure 4. Numerical simulations of Caputo model (6) and ABC model (16) with parameter values (29) and n = 0.2: Figure (a) phase portrait; and Figure (b) time series. y(t) y(t) x(t) y(t) x(t) y(t) Axioms 2020, 9, 122 18 of 23 1.1 1.0 Caputo ABC THP Caputo ABC 1.0 0.8 0.9 0.6 0.8 0.4 0.7 0.8 0.6 * 0.6 0.5 0.4 0.4 0.2 Caputo ABC 0.3 0.0 0.0 0.2 0.4 0.6 0.8 1.0 0 100 200 300 400 500 x(t) t (a) Phase portrait (b) n = 0.14 Figure 5. Numerical simulations of Caputo model (6) and ABC model (16) with parameter values (29), n = 0.2 and T = 0.5. Caputo 0.8 Upper bound of limit-cycle ABC 0.7 The difference in stability of Caputo and ABC 0.6 E-stable = 0.7 = 0.772 = 0.83 E-unstable 0.5 0.4 Lower bound of limit-cycle 0.3 Caputo 0.6 Upper bound of limit-cycle ABC 0.5 The difference in stability of Caputo and ABC 0.4 0.3 E-stable = 0.7 = 0.772 = 0.83 E-unstable 0.2 0.1 Lower bound of limit-cycle 0.60 0.65 0.70 0.75 0.80 0.85 0.90 Figure 6. Bifurcation diagram of Caputo model (6) and ABC model (16) driven by the order of the fractional-derivative (a) with constant parameter values (29) and n = 0.2. There exists a Hopf bifurcation where the bifurcation points of the Caputo model and ABC model are different. 0.8 0.6 E E Caputo Caputo 0.4 = 0.772 = 0.83 Caputo = 0.7 Limit-cycle Limit-cycle 0.2 0.8 0.6 ABC 0.4 = 0.83 ABC ABC = 0.7 = 0.772 Limit-cycle 0.2 0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 x(t) x(t) x(t) (a) Figure 7. Cont. y(t) y(t) y(t) x(t) y(t) x(t) y(t) Axioms 2020, 9, 122 19 of 23 0.8 0.6 0.4 Caputo Caputo Caputo = 0.7 = 0.772 = 0.83 0.2 x(t) x(t) x(t) y(t) y(t) y(t) 0.0 0.8 0.6 0.4 ABC ABC ABC = 0.7 = 0.772 = 0.83 0.2 x(t) x(t) x(t) y(t) y(t) y(t) 0.0 0 500 1000 1500 2000 0 500 1000 1500 2000 0 500 1000 1500 2000 t t t (b) Figure 7. Numerical simulations of Caputo model (6) and ABC model (16) with parameter values (29), n = 0.2, and a = f0.7, 0.772, 0.83g: Figure (a) phase portrait; and Figure (b) time series. From the biological point of view, all previous numerical simulations show that the dynamical properties of both Caputo model and ABC model are similar except when the eigenvalues of the Jacobian matrix evaluated at the interior equilibrium point E are a pair of complex conjugate with positive real part. There is a biological condition such that the prey and predator densities are eventually periodic for the Caputo model, while for ABC model, the densities of both predator and prey are eventually constant. 5. Conclusions The dynamics of a Rosenzweig–MacArthur model with continuous threshold harvesting in predator involving the Caputo fractional-order derivative and ABC fractional-order derivative are studied. We prove the existence and uniqueness of solutions of both Caputo and ABC models. Particularly, we completely investigate the dynamics of the Caputo model including the non-negativity, boundedness, local stability, global stability, and the existence of Hopf bifurcation. From the biological meanings, the extinction of both populations never occurs since the origin point (E ) is a saddle point. Some of the situations that might occur are: (1) the predator goes extinct while the prey still survives, which is indicated by the stability of E ; (2) both predator and prey co-exist and converge to a constant population density, which happens if the interior point E or E are asymptotically stable; and (3) both predator and prey co-exist where both population densities change periodically, namely when a Hopf bifurcation occurs. We show numerically that our model may undergo a forward bifurcation or a Hopf bifurcation. The Hopf bifurcation in models with both Caputo operator and ABC operator can be driven by the conversion rate of consumed prey into the predator birth rate or by the order of fractional derivative. Our numerical simulations show that the Hopf bifurcation point of both models are different. Author Contributions: All authors equally contributed to this article. All authors read and approved the ﬁnal manuscript. Funding: This research was funded by the Directorate of Research and Community Service, The Directorate General of Strengthening Research and Development, the Ministry of Research, Technology, and Higher Education (Brawijaya University), Indonesia, via Doctoral Dissertation Research, in accordance with the Research Contract No. 037/SP2H/LT/DRPM/2020, dated 9 March 2020. Acknowledgments: The authors would like to thank the reviewers for all useful and helpful comments to improve the manuscript. Conﬂicts of Interest: All authors report no conﬂict of interest relevant to this article. x(t), y(t) x(t), y(t) Axioms 2020, 9, 122 20 of 23 References 1. Rosenzweig, M.L.; MacArthur, R.H. Graphical representation and stability conditions of predator-prey interactions. Am. Nat. 1963, 97, 209–223. doi:10.1086/282272. 2. Moustafa, M.; Mohd, M.H.; Ismail, A.I.; Abdullah, F.A. Stage structure and refuge effects in the dynamical analysis of a fractional order Rosenzweig-MacArthur prey-predator model. Prog. Fract. Differ. Appl. 2019, 5, 49–64. doi:10.18576/pfda/050106. 3. Beay, L.K.; Suryanto, A.; Darti, I.; Trisilowati, T. Hopf bifurcation and stability analysis of the Rosenzweig-MacArthur predator-prey model with stage-structure in prey. Math. Biosci. Eng. 2020, 17, 4080–4097. doi:10.3934/mbe.2020226. 4. González-Olivares, E.; Ramos-Jiliberto, R. Dynamic consequences of prey refuges in a simple model system: More prey, fewer predators and enhanced stability. Ecol. Model. 2003, 166, 135–146. doi:10.1016/S0304-3800(03)00131-5. 5. Chen, L.; Chen, F.; Chen, L. Qualitative analysis of a predator–prey model with Holling type II functional response incorporating a constant prey refuge. Nonlinear Anal. Real World Appl. 2010, 11, 246–252. doi:10.1016/j.nonrwa.2008.10.056. 6. Almanza-Vasquez, E.; Ortiz-Ortiz, R.D.; Marin-Ramirez, A.M. Bifurcations in the dynamics of Rosenzweig-Macarthur predator-prey model considering saturated refuge for the preys. Appl. Math. Sci. 2015, 9, 7475–7482. doi:10.12988/ams.2015.510640. 7. Das, M.; Maiti, A.; Samanta, G.P. Stability analysis of a prey-predator fractional order model incorporating prey refuge. Ecol. Genet. Genom. 2018, 7-8, 33–46. doi:10.1016/j.egg.2018.05.001. 8. Moustafa, M.; Mohd, M.H.; Ismail, A.I.; Abdullah, F.A. Dynamical analysis of a fractional-order Rosenzweig–MacArthur model incorporating a prey refuge. Chaos Solitons Fractals 2018, 109, 1–13. doi:10.1016/j.chaos.2018.02.008. 9. Sarkar, K.; Khajanchi, S. Impact of fear effect on the growth of prey in a predator-prey interaction model. Ecol. Complex. 2020, 42, 100826. doi:10.1016/j.ecocom.2020.100826. 10. Zu, J.; Mimura, M. The impact of Allee effect on a predator-prey system with Holling type II functional response. Appl. Math. Comput. 2010, 217, 3542–3556. doi:10.1016/j.amc.2010.09.029. 11. Pal, P.J.; Saha, T. Qualitative analysis of a predator-prey system with double Allee effect in prey. Chaos Solitons Fractals 2015, 73, 36–63. doi:10.1016/j.chaos.2014.12.007. 12. Bodine, E.N.; Yust, A.E. Predator–prey dynamics with intraspeciﬁc competition and an Allee effect in the predator population. Lett. Biomath. 2017, 4, 23–38. doi:10.1080/23737867.2017.1282843. 13. Mukherjee, D. Role of fear in predator–prey system with intraspeciﬁc competition. Math. Comput. Simul. 2020, 177, 263–275. doi:10.1016/j.matcom.2020.04.025. 14. Van Den Bosch, F. Cannibalism in an age-structured predator-prey system. Bull. Math. Biol. 1997, 59, 551–567. doi:10.1016/s0092-8240(96)00107-3. 15. Mondal, S.; Lahiri, A.; Bairagi, N. Analysis of a fractional order eco-epidemiological model with prey infection and type 2 functional response. Math. Methods Appl. Sci. 2017, 40, 6776–6789. doi:10.1002/mma.4490. 16. Suryanto, A.; Darti, I.; Anam, S. Stability analysis of pest-predator interaction model with infectious disease in prey. AIP Conf. Proc. 2018, 1937, 020018. doi:10.1063/1.5026090. 17. Panigoro, H.S.; Suryanto, A.; Kusumawinahyu, W.M.; Darti, I. Dynamics of a fractional-order predator-prey model with infectious diseases in prey. Commun. Biomath. Sci. 2019, 2, 105–117. doi:10.5614/cbms.2019.2.2.4. 18. Kumar, T.; And, K.; Chakraborty, K. Effort dynamics in a prey-predator model with harvesting. Int. J. Inf. Syst. Sci. 2000, 6, 318–332. 19. Javidi, M.; Nyamoradi, N. Dynamic analysis of a fractional order prey-predator interaction with harvesting. Appl. Math. Model. 2013, 37, 8946–8956. doi:10.1016/j.apm.2013.04.024. 20. Zhu, C.; Kong, L. Bifurcations analysis of Leslie-Gower predator-prey models with nonlinear predator-harvesting. Discret. Contin. Dyn. Syst. 2017, 10, 1187–1206. doi:10.3934/dcdss.2017065. 21. Suryanto, A.; Darti, I. Dynamics of Leslie-Gower pest-predator model with disease in pest including pest-harvesting and optimal implementation of pesticide. Int. J. Math. Math. Sci. 2019, 2019, 1–9. doi:10.1155/2019/5079171. Axioms 2020, 9, 122 21 of 23 22. Ang, T.K.; Safuan, H.M. Harvesting in a toxicated intraguild predator–prey ﬁshery model with variable carrying capacity. Chaos Solitons Fractals 2019, 126, 158–168. doi:10.1016/j.chaos.2019.06.004. 23. Leard, B.; Rebaza, J. Analysis of predator-prey models with continuous threshold harvesting. Appl. Math. Comput. 2011, 217, 5265–5278. doi:10.1016/j.amc.2010.11.050. 24. Bohn, J.; Rebaza, J.; Speer, K. Continuous threshold prey harvesting in predator-prey models. Int. J. Math. Comput. Sci. 2011, 5, 964–971. 25. Rebaza, J. Dynamics of prey threshold harvesting and refuge. J. Comput. Appl. Math. 2012, 236, 1743–1752. doi:10.1016/j.cam.2011.10.005. 26. Lv, Y.; Yuan, R.; Pei, Y. Dynamics in two nonsmooth predator–prey models with threshold harvesting. Nonlinear Dyn. 2013, 74, 107–132. doi:10.1007/s11071-013-0952-2. 27. Wu, D.; Zhao, H.; Yuan, Y. Complex dynamics of a diffusive predator–prey model with strong Allee effect and threshold harvesting. J. Math. Anal. Appl. 2019, 469, 982–1014. doi:10.1016/j.jmaa.2018.09.047. 28. Toaha, S. The effect of harvesting with threshold on the dynamics of prey predator model. J. Phys. Conf. Ser. 2019, 1341. doi:10.1088/1742-6596/1341/6/062021. 29. Panigoro, H.S.; Suryanto, A.; Kusumawinahyu, W.M.; Darti, I. Continuous threshold harvesting in a gause-type predator-prey model with fractional-order. AIP Conf. Proc. 2020, 2264, 040001. doi:10.1063/5.0023513. 30. Shepard, B.M.; Carner, G.R.; Barrion, A.T.; Ooi, P.A.C.; Van den Berg, H. Insects and Their Natural Enemies Associated with Vegetables and Soybean in Southeast Asia; Clemson Univ Coastal Research &: Orangeburg, USA, 31. Li, H.L.; Zhang, L.; Hu, C.; Jiang, Y.L.; Teng, Z. Dynamical analysis of a fractional-order predator-prey model incorporating a prey refuge. J. Appl. Math. Comput. 2017, 54, 435–449. doi:10.1007/s12190-016-1017-8. 32. Das, S.; Gupta, P.K. A mathematical model on fractional Lotka-Volterra equations. J. Theor. Biol. 2011, 277, 1–6. doi:10.1016/j.jtbi.2011.01.034. 33. Panja, P. Dynamics of a fractional order predator-prey model with intraguild predation. Int. J. Model. Simul. 2019, 39, 256–268. doi:10.1080/02286203.2019.1611311. 34. Suryanto, A.; Darti, I.; S. Panigoro, H.; Kilicman, A. A fractional-order predator–prey model with ratio-dependent functional response and linear harvesting. Mathematics 2019, 7, 1100. doi:10.3390/math7111100. 35. Alidousti, J.; Ghafari, E. Dynamic behavior of a fractional order prey-predator model with group defense. Chaos Solitons Fractals 2020, 134, 109688. doi:10.1016/j.chaos.2020.109688. 36. Ghanbari, B.; Djilali, S. Mathematical analysis of a fractional-order predator-prey model with prey social behavior and infection developed in predator population. Chaos Solitons Fractals 2020, 138, 109960. doi:10.1016/j.chaos.2020.109960. 37. Xie, Y.; Wang, Z.; Meng, B.; Huang, X. Dynamical analysis for a fractional-order prey–predator model with Holling III type functional response and discontinuous harvest. Appl. Math. Lett. 2020, 106, 106342. doi:10.1016/j.aml.2020.106342. 38. Sekerci, Y. Climate change effects on fractional order prey-predator model. Chaos Solitons Fractals 2020, 134, 109690. doi:10.1016/j.chaos.2020.109690. 39. Caputo, M. Linear models of dissipation whose Q is almost fFrequency independent–II. Geophys. J. Int. 1967, 13, 529–539. doi:10.1111/j.1365-246X.1967.tb02303.x. 40. Atangana, A.; Baleanu, D. New fractional derivatives with nonlocal and non-singular kernel: Theory and application to heat transfer model. Therm. Sci. 2016, 20, 763–769. doi:10.2298/TSCI160111018A. 41. Atangana, A.; Koca, I. Chaos in a simple nonlinear system with Atangana–Baleanu derivatives with fractional order. Chaos Solitons Fractals 2016, 89, 447–454. doi:10.1016/j.chaos.2016.02.012. 42. Tajadodi, H. A Numerical approach of fractional advection-diffusion equation with Atangana–Baleanu derivative. Chaos Solitons Fractals 2020, 130, 109527. doi:10.1016/j.chaos.2019.109527. 43. Ghanbari, B.; Kumar, D. Numerical solution of predator-prey model with Beddington-DeAngelis functional response and fractional derivatives with Mittag-Lefﬂer kernel. Chaos Interdiscip. J. Nonlinear Sci. 2019, 29, 063103. doi:10.1063/1.5094546. 44. Morales-Delgado, V.F.; Gómez-Aguilar, J.F.; Saad, K.; Escobar Jiménez, R.F. Application of the Caputo-Fabrizio and Atangana-Baleanu fractional derivatives to mathematical model of cancer chemotherapy effect. Math. Methods Appl. Sci. 2019, 42, 1167–1193. doi:10.1002/mma.5421. Axioms 2020, 9, 122 22 of 23 45. Shah, S.A.A.; Khan, M.A.; Farooq, M.; Ullah, S.; Alzahrani, E.O. A fractional order model for Hepatitis B virus with treatment via Atangana–Baleanu derivative. Phys. A Stat. Mech. Its Appl. 2020, 538, 122636. doi:10.1016/j.physa.2019.122636. 46. Bourafa, S.; Abdelouahab, M.S.; Moussaoui, A. On some extended Routh–Hurwitz conditions for fractional-order autonomous systems of order a 2 (0, 2) and their applications to some population dynamic models. Chaos Solitons Fractals 2020, 133. doi:10.1016/j.chaos.2020.109623. 47. Ahmed, E.; El-Sayed, A.; El-Saka, H.A. On some Routh–Hurwitz conditions for fractional order differential equations and their applications in Lorenz, Rössler, Chua and Chen systems. Phys. Lett. A 2006, 358, 1–4. doi:10.1016/j.physleta.2006.04.087. 48. Petras, I. Fractional-Order Nonlinear Systems: Modeling, Analysis and Simulation; Springer: London, UK; Beijing, China, 2011. 49. Diethelm, K. The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using Differential Operators of Caputo Type; Springer: Braunschweig, Germany, 2010. 50. Podlubny, I. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications; Academic Press: San Diego CA, USA, 1999. 51. Vargas-De-León, C. Volterra-type Lyapunov functions for fractional-order epidemic systems. Commun. Nonlinear Sci. Numer. Simul. 2015, 24, 75–85. doi:10.1016/j.cnsns.2014.12.013. 52. Huo, J.; Zhao, H.; Zhu, L. The effect of vaccines on backward bifurcation in a fractional order HIV model. Nonlinear Anal. Real World Appl. 2015, 26, 289–305. doi:10.1016/j.nonrwa.2015.05.014. 53. Abdelouahab, M.S.; Hamri, N.E.; Wang, J. Hopf bifurcation and caos in fractional-order modiﬁed hybrid optical system. Nonlinear Dyn. 2012, 69, 275–284. doi:10.1007/s11071-011-0263-4. 54. Deshpande, A.S.; Daftardar-Gejji, V.; Sukale, Y.V. On Hopf bifurcation in fractional dynamical systems. Chaos Solitons Fractals 2017, 98, 189–198. doi:10.1016/j.chaos.2017.03.034. 55. Diethelm, K. A fractional calculus based model for the simulation of an outbreak of dengue fever. Nonlinear Dyn. 2013, 71, 613–619. doi:10.1007/s11071-012-0475-2. 56. Moustafa, M.; Mohd, M.H.; Ismail, A.I.; Abdullah, F.A. Dynamical analysis of a fractional-order eco-epidemiological model with disease in prey population. Adv. Differ. Equ. 2020, 2020, 48. doi:10.1186/s13662-020-2522-5. 57. Li, Y.; Chen, Y.; Podlubny, I. Stability of fractional-order nonlinear dynamic systems: Lyapunov direct method and generalized Mittag–Lefﬂer stability. Comput. Math. Appl. 2010, 59, 1810–1821. doi:10.1016/j.camwa.2009.08.019. 58. Odibat, Z.M.; Shawagfeh, N.T. Generalized Taylor ’s formula. Appl. Math. Comput. 2007, 186, 286–293. doi:10.1016/j.amc.2006.07.102. 59. Matignon, D. Stability results for fractional differential equations with applications to control processing. Comput. Eng. Syst. Appl. 1996, pp. 963–968. 60. Kuznetsov, Y.A. Elements of Applied Bifurcation Theory, 3rd ed.; Springer: New York, NY, USA, 2004. 61. Baisad, K.; Moonchai, S. Analysis of stability and Hopf bifurcation in a fractional Gauss-type predator–prey model with Allee effect and Holling type-III functional response. Adv. Differ. Equ. 2018, 2018. doi:10.1186/s13662-018-1535-9. 62. Li, X.; Wu, R. Hopf bifurcation analysis of a new commensurate fractional-order hyperchaotic system. Nonlinear Dyn. 2014, 78, 279–288. doi:10.1007/s11071-014-1439-5. 63. Tavazoei, M.S.; Haeri, M. A proof for non existence of periodic solutions in time invariant fractional order systems. Automatica 2009, 45, 1886–1890. doi:10.1016/j.automatica.2009.04.001. 64. Fulger, D.; Scalas, E.; Germano, G. Monte Carlo simulation of uncoupled continuous-time random walks yielding a stochastic solution of the space-time fractional diffusion equation. Phys. Rev. E 2008, 77, 021122. doi:10.1103/PhysRevE.77.021122. 65. Scherer, R.; Kalla, S.L.; Tang, Y.; Huang, J. The Grünwald–Letnikov method for fractional differential equations. Comput. Math. Appl. 2011, 62, 902–917. doi:10.1016/j.camwa.2011.03.054. 66. Suryanto, A.; Darti, I. Stability analysis and nonstandard Grünwald-Letnikov scheme for a fractional order predator-prey model with ratio-dependent functional response. AIP Conf. Proc. 2017, 1913, 020011. doi:10.1063/1.5016645. 67. Diethelm, K.; Ford, N.J.; Freed, A.D. A predictor-corrector approach for the numerical solution of fractional differential equations. Nonlinear Dyn. 2002, 29, 3–22. doi:10.1023/A:1016592219341. Axioms 2020, 9, 122 23 of 23 68. Wang, Z. A numerical method for delayed fractional-order differential equations. J. Appl. Math. 2013, 2013, 256071. doi:10.1155/2013/256071. 69. Baleanu, D.; Jajarmi, A.; Hajipour, M. On the nonlinear dynamical systems within the generalized fractional derivatives with Mittag–Lefﬂer kernel. Nonlinear Dyn. 2018, 94, 397–414. doi:10.1007/s11071-018-4367-y. Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional afﬁliations. © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Axioms Multidisciplinary Digital Publishing Institute http://www.deepdyve.com/lp/multidisciplinary-digital-publishing-institute/a-rosenzweig-macarthur-model-with-continuous-threshold-harvesting-in-bnyCmdp5KY

Loading next page...

References (71)

John Bohn, J. Rebaza, K. Speer (2011)
Continuous Threshold Prey Harvesting in Predator-Prey Models
World Academy of Science, Engineering and Technology, International Journal of Mathematical, Computational, Physical, Electrical and Computer Engineering, 5
B. Ghanbari, Devendra Kumar (2019)
Numerical solution of predator-prey model with Beddington-DeAngelis functional response and fractional derivatives with Mittag-Leffler kernel.
Chaos, 29 6
Yan Li, Y. Chen, I. Podlubny (2010)
Stability of fractional-order nonlinear dynamic systems: Lyapunov direct method and generalized Mittag-Leffler stability
Comput. Math. Appl., 59
I. Podlubny (1999)
Fractional differential equations : an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications
M. Tavazoei, M. Haeri (2009)
A proof for non existence of periodic solutions in time invariant fractional order systems
Autom., 45
M. Rosenzweig, R. Macarthur (1963)
Graphical Representation and Stability Conditions of Predator-Prey Interactions
The American Naturalist, 97
Tau Ang, H. Safuan (2019)
Harvesting in a toxicated intraguild predator–prey fishery model with variable carrying capacity
Chaos, Solitons & Fractals
A. Suryanto, I. Darti, Hasan Panigoro, A. Kılıçman (2019)
A Fractional-Order Predator–Prey Model with Ratio-Dependent Functional Response and Linear Harvesting
Mathematics
Y. Kuznetsov (1998)
Elements of applied bifurcation theory (2nd ed.)
Meghadri Das, A. Maiti, G. Samanta (2018)
Stability analysis of a prey-predator fractional order model incorporating prey refuge
Ecological Genetics and Genomics
C. Vargas‐De‐León (2015)
Volterra-type Lyapunov functions for fractional-order epidemic systems
Communications in Nonlinear Science and Numerical Simulation, 24
H. Panigoro, A. Suryanto, Wuryansari Kusumahwinahyu, I. Darti (2019)
Dynamics of a Fractional-Order Predator-Prey Model with Infectious Diseases in Prey
, 2
P. Pal, Tapan Saha (2015)
Qualitative analysis of a predator–prey system with double Allee effect in prey
Chaos Solitons & Fractals, 73
Van den Berg, H. Insects and Their Natural Enemies Associated with Vegetables and Soybean
M. Moustafa, M. Mohd, A. Ismail, F. Abdullah (2020)
Dynamical analysis of a fractional-order eco-epidemiological model with disease in prey population
Advances in Difference Equations, 2020
Subir Das, P. Gupta (2011)
A mathematical model on fractional Lotka-Volterra equations.
Journal of theoretical biology, 277 1
Liujuan Chen, Fengde Chen, Lijuan Chen (2010)
Qualitative analysis of a predator–prey model with Holling type II functional response incorporating a constant prey refuge
Nonlinear Analysis-real World Applications, 11
A. Suryanto, I. Darti (2019)
Dynamics of Leslie-Gower Pest-Predator Model with Disease in Pest Including Pest-Harvesting and Optimal Implementation of Pesticide
Int. J. Math. Math. Sci., 2019
Shuvojit Mondal, A. Lahiri, N. Bairagi (2017)
Analysis of a fractional order eco‐epidemiological model with prey infection and type 2 functional response
Mathematical Methods in the Applied Sciences, 40
D. Matignon (1996)
Stability results for fractional differential equations with applications to control processing
D. Baleanu, A. Jajarmi, M. Hajipour (2018)
On the nonlinear dynamical systems within the generalized fractional derivatives with Mittag–Leffler kernel
Nonlinear Dynamics, 94
K. Chakraborty (2010)
EFFORT DYNAMICS IN A PREY-PREDATOR MODEL WITH HARVESTING
J. Kuznecov (1998)
Elements of applied bifurcation theory
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional afﬁliations. c (cid:13) 2020 by the authors. Licensee MDPI, Basel, Switzerland
A. Atangana, I. Koca (2016)
Chaos in a simple nonlinear system with Atangana-Baleanu derivatives with fractional order
Chaos Solitons & Fractals, 89
Benjamin Leard, J. Rebaza (2011)
Analysis of predator-prey models with continuous threshold harvesting
Appl. Math. Comput., 217
E. Ahmed, A. El-Sayed, H. El-Saka (2006)
On some Routh–Hurwitz conditions for fractional order differential equations and their applications in Lorenz, Rössler, Chua and Chen systems
Physics Letters A, 358
Jingjing Huo, Hongyong Zhao, Linhe Zhu (2015)
The effect of vaccines on backward bifurcation in a fractional order HIV model
Nonlinear Analysis-real World Applications, 26
J. Zu, M. Mimura (2010)
The impact of Allee effect on a predator-prey system with Holling type II functional response
Appl. Math. Comput., 217
Syed Shah, M. Khan, Muhammad Farooq, Saif Ullah, E. Alzahrani (2020)
A fractional order model for Hepatitis B virus with treatment via Atangana–Baleanu derivative
Physica A-statistical Mechanics and Its Applications, 538
M. Moustafa, M. Mohd, A. Ismail, F. Abdullah (2019)
Stage Structure and Refuge Effects in the Dynamical Analysis of a Fractional Order Rosenzweig-MacArthur Prey-Predator Model
Progress in Fractional Differentiation and Applications
K. Diethelm, N. Ford, A. Freed (2013)
A Predictor-Corrector Approach for the Numerical Solution of Fractional Differential Equations
Nonlinear Dynamics, 29
Xiang Li, R. Wu (2014)
Hopf bifurcation analysis of a new commensurate fractional-order hyperchaotic system
Nonlinear Dynamics, 78
S. Bourafa, Mohammed Abdelouahab, A. Moussaoui (2020)
On some extended Routh–Hurwitz conditions for fractional-order autonomous systems of order α ∈ (0, 2) and their applications to some population dynamic models
Chaos Solitons & Fractals, 133
P. Panja (2019)
Dynamics of a fractional order predator-prey model with intraguild predation
International Journal of Modelling and Simulation, 39
D. Fulger, E. Scalas, G. Germano (2007)
Monte Carlo simulation of uncoupled continuous-time random walks yielding a stochastic solution of the space-time fractional diffusion equation.
Physical review. E, Statistical, nonlinear, and soft matter physics, 77 2 Pt 1
K. Diethelm (2013)
A fractional calculus based model for the simulation of an outbreak of dengue fever
Nonlinear Dynamics, 71
C. Zhu, Lei Kong (2017)
Bifurcations analysis of Leslie-Gower predator-prey models with nonlinear predator-harvesting
Discrete and Continuous Dynamical Systems - Series S, 10
Zhen Wang (2013)
A Numerical Method for Delayed Fractional-Order Differential Equations
J. Appl. Math., 2013
H. Panigoro, A. Suryanto, W. Kusumawinahyu, I. Darti (2020)
Continuous threshold harvesting in a gause-type predator-prey model with fractional-order
B. Ghanbari, S. Djilali (2020)
Mathematical analysis of a fractional-order predator-prey model with prey social behavior and infection developed in predator population
Chaos Solitons & Fractals, 138
D. Mukherjee (2020)
Role of fear in predator-prey system with intraspecific competition
Math. Comput. Simul., 177
Yunfei Lv, R. Yuan, Yongzhen Pei (2013)
Dynamics in two nonsmooth predator–prey models with threshold harvesting
Nonlinear Dynamics, 74
R. Scherer, S. Kalla, Yifa Tang, Jianfei Huang (2011)
The Grünwald-Letnikov method for fractional differential equations
Comput. Math. Appl., 62
M. Moustafa, M. Mohd, A. Ismail, F. Abdullah (2018)
Dynamical analysis of a fractional-order Rosenzweig–MacArthur model incorporating a prey refuge
Chaos Solitons & Fractals, 109
Yingkang Xie, Zhen Wang, Bo Meng, Xia Huang (2020)
Dynamical analysis for a fractional-order prey-predator model with Holling III type functional response and discontinuous harvest
Appl. Math. Lett., 106
Daiyong Wu, Hongyong Zhao, Yuan Yuan (2019)
Complex dynamics of a diffusive predator–prey model with strong Allee effect and threshold harvesting
Journal of Mathematical Analysis and Applications
Edilber Almanza-Vasquez, Ruben-Dario Ortiz-Ortiz, Ana-Magnolia Marin-Ramirez (2015)
Bifurcations in the dynamics of Rosenzweig-MacArthur predator-prey model considering saturated refuge for the preys
Applied mathematical sciences, 9
I. Petráš (2011)
Fractional-Order Nonlinear Systems: Modeling, Analysis and Simulation
K. Diethelm (2010)
The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using Differential Operators of Caputo Type
E. González‐Olivares, R. Ramos‐Jiliberto (2003)
Dynamic consequences of prey refuges in a simple model system: more prey, fewer predators and enhanced stability
Ecological Modelling, 166
J. Rebaza (2012)
Dynamics of prey threshold harvesting and refuge
J. Comput. Appl. Math., 236
F. Bosch, W. Gabriel (1997)
Cannibalism in an age-structured predator-prey system.
Bulletin of Mathematical Biology, 59
Erin Bodine, Anne Yust (2017)
Predator–prey dynamics with intraspecific competition and an Allee effect in the predator population
Letters in Biomathematics, 4
Hongli Li, Long Zhang, Cheng Hu, Yaolin Jiang, Z. Teng (2016)
Dynamical analysis of a fractional-order predator-prey model incorporating a prey refuge
Journal of Applied Mathematics and Computing, 54
Mohammed Abdelouahab, N. Hamri, Junwei Wang (2011)
Hopf bifurcation and chaos in fractional-order modified hybrid optical system
Nonlinear Dynamics, 69
M. Caputo (1967)
Linear Models of Dissipation whose Q is almost Frequency Independent-II
Geophysical Journal International, 13
Y. Sekerci (2020)
Climate change effects on fractional order prey-predator model
Chaos Solitons & Fractals, 134
S. Toaha (2019)
The effect of harvesting with threshold on the dynamics of prey predator model
Journal of Physics: Conference Series, 1341
A. Suryanto, I. Darti, S. Anam (2018)
Stability analysis of pest-predator interaction model with infectious disease in prey
, 1937
A. Atangana, D. Baleanu (2016)
New Fractional Derivatives with Nonlocal and Non-Singular Kernel: Theory and Application to Heat Transfer Model
arXiv: General Mathematics
Lazarus Beay, A. Suryanto, I. Darti, T. Trisilowati (2020)
Hopf bifurcation and stability analysis of the Rosenzweig-MacArthur predator-prey model with stage-structure in prey.
Mathematical biosciences and engineering : MBE, 17 4
K. Sarkar, S. Khajanchi (2020)
Impact of fear effect on the growth of prey in a predator-prey interaction model
Ecological Complexity, 42
Amey Deshpande, Varsha Daftardar-Gejji, Yogita Sukale (2017)
On Hopf bifurcation in fractional dynamical systems
Chaos Solitons & Fractals, 98
H. Tajadodi (2020)
A Numerical approach of fractional advection-diffusion equation with Atangana–Baleanu derivative
Chaos Solitons & Fractals, 130
M. Javidi, N. Nyamoradi (2013)
Dynamic analysis of a fractional order prey–predator interaction with harvesting
Applied Mathematical Modelling, 37
A. Suryanto, I. Darti (2017)
Stability analysis and nonstandard Grünwald-Letnikov scheme for a fractional order predator-prey model with ratio-dependent functional response
, 1913
Z. Odibat, N. Shawagfeh (2007)
Generalized Taylor's formula
Appl. Math. Comput., 186
J. Alidousti, Elham Ghafari (2020)
Dynamic behavior of a fractional order prey-predator model with group defense
Chaos Solitons & Fractals, 134
V. Morales-Delgado, J. Gómez‐Aguilar, K. Saad, Ricardo Jimenez (2019)
Application of the Caputo‐Fabrizio and Atangana‐Baleanu fractional derivatives to mathematical model of cancer chemotherapy effect
Mathematical Methods in the Applied Sciences, 42
Kanokrat Baisad, Sompop Moonchai (2018)
Analysis of stability and Hopf bifurcation in a fractional Gauss-type predator–prey model with Allee effect and Holling type-III functional response
Advances in Difference Equations, 2018

Publisher: Multidisciplinary Digital Publishing Institute
Copyright: © 1996-2020 MDPI (Basel, Switzerland) unless otherwise stated Disclaimer The statements, opinions and data contained in the journals are solely those of the individual authors and contributors and not of the publisher and the editor(s). MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Terms and Conditions Privacy Policy
ISSN: 2075-1680
DOI: 10.3390/axioms9040122
Publisher site: See Article on Publisher Site

Abstract

axioms Article A Rosenzweig–MacArthur Model with Continuous Threshold Harvesting in Predator Involving Fractional Derivatives with Power Law and Mittag–Lefﬂer Kernel 1,2 1, 1 Hasan S. Panigoro , Agus Suryanto * , Wuryansari Muharini Kusumawinahyu and Isnani Darti Department of Mathematics, Faculty of Mathematics and Natural Sciences, University of Brawijaya, Malang 65145, Indonesia; hspanigoro@ung.ac.id (H.S.P.); wmuharini@ub.ac.id (W.M.K.); isnanidarti@ub.ac.id (I.D.) Department of Mathematics, Faculty of Mathematics and Natural Sciences, State University of Gorontalo, Bone Bolango 96119, Indonesia * Correspondence: suryanto@ub.ac.id Received: 17 September 2020; Accepted: 19 October 2020; Published: 22 October 2020 Abstract: The harvesting management is developed to protect the biological resources from over-exploitation such as harvesting and trapping. In this article, we consider a predator–prey interaction that follows the fractional-order Rosenzweig–MacArthur model where the predator is harvested obeying a threshold harvesting policy (THP). The THP is applied to maintain the existence of the population in the prey–predator mechanism. We ﬁrst consider the Rosenzweig–MacArthur model using the Caputo fractional-order derivative (that is, the operator with the power-law kernel) and perform some dynamical analysis such as the existence and uniqueness, non-negativity, boundedness, local stability, global stability, and the existence of Hopf bifurcation. We then reconsider the same model involving the Atangana–Baleanu fractional derivative with the Mittag–Lefﬂer kernel in the Caputo sense (ABC). The existence and uniqueness of the solution of the model with ABC operator are established. We also explore the dynamics of the model with both fractional derivative operators numerically and conﬁrm the theoretical ﬁndings. In particular, it is shown that models with both Caputo operator and ABC operator undergo a Hopf bifurcation that can be controlled by the conversion rate of consumed prey into the predator birth rate or by the order of fractional derivative. However, the bifurcation point of the model with the Caputo operator is different from that of the model with the ABC operator. Keywords: Rosenzweig–MacArthur model; fractional derivatives; threshold harvesting 1. Introduction More than 50 years after the model has been proposed, the Rosenzweig–MacArthur predator–prey model [1] has been consistently developed by many scholars to approach the real world phenomena with more realistic mathematical models. The commonsensical modiﬁed Rosenzweig–MacArthur models are accomplishable by considering the biological perspectives of ecosystem conditions, for instance the stage structure [2,3], the refuge effect [4–8], the fear effect [9], the Allee effect [10,11], the intraspeciﬁc competition [12,13], the cannibalism [14], the infectious disease [15–17], and so forth. On the other hand, the modeling also contemplates the optimal management of bioeconomic resources as in ﬁshery and pest management. Some researchers put an intervention into the predator–prey model such as the harvesting to one or more population [8,18–22]. To protect Axioms 2020, 9, 122; doi:10.3390/axioms9040122 www.mdpi.com/journal/axioms Axioms 2020, 9, 122 2 of 23 the population from over-exploitation during the harvesting, some management techniques have been established. One of the famous technique is a continuous threshold harvesting policy (THP) (see [23–29]). THP is regulated as follows: when the population density above the threshold level, harvesting is permitted; when the population density drops below the threshold level, harvesting is prohibited. In 2013, Lv et al. [26] proposed the following Rosenzweig–MacArthur model with THP in predator dx x mxy =rx 1 , dt K a + x (1) dy nxy = dy H(y), dt a + x where < 0 , if y < T, H(y) = h(y T) : , if y T. c + (y T) We describe the biological interpretation of variables and parameters of model (1) in Table 1. Model (1) represents an interaction between two populations with a prey–predator relationship, where THP is only applied for the predator to preserving its populations. Some appealing examples of the ecological model (1) are given by the interaction between Sycanus sp. and Setothosea asigna and between Rhinocoris sp. and Spodoptera litura. Shepard [30] reported that Sycanus sp. and Rhinocoris sp. are the natural predators of the pests such as Setothosea asigna and Spodoptera litura which exist in agricultural lands and plantations. The worrying problem is: How if the density of insects such as Sycanus sp. and Rhinocoris sp. uncontrolled? One solution is applying THP as in model (1). Table 1. Description of variables and parameters of the model (1). Variables and Parameters Description x(t) The density of prey y(t) The density of predator r The intrinsic growth rate of prey K The environmental carrying capacity of prey m The maximum uptake rate for prey n The conversion rate of consumed prey into predator birth a The environment protection for prey d The natural death rate of predator h The harvesting rate c The half saturation constant for harvesting T The threshold level of harvesting Lv et.al. [26] successfully explored the dynamics of the model (1) including the local stability and the existence of various phenomena (saddle-node, Hopf, cusp, and Bogdanov–Taken bifurcations). Despite their success works, the model with the first-order derivative is limited to its capability of involving all previous conditions to the growth rates of both predator and prey. The growth rates of both populations in the model (1) depend only on the current state. Biologically, the growth rates must be dependent on all of the previous conditions which are known as the memory effects. To account for such memory effects, some researchers proposed to apply the fractional-order derivative instead of the first-order derivative when expressing the growth rate of the population. The fractional-order models are naturally related to systems with memory which exists in most biological models [7,31]. The fractional-order models are also well-liked due to their capability in providing an exact description of different nonlinear phenomena [32]. In recent years, the development of fractional-order models grows rapidly and becomes popular in studying the dynamical behavior of predator–prey interaction [17,33–38]. It has been shown that the order of the fractional derivate signiﬁcantly affects the dynamical behavior of Axioms 2020, 9, 122 3 of 23 the models, which is in contrast to the ﬁrst-order derivative models that depend only on the values of parameters. In this paper, we modify the model of Lv et.al. [26] by applying the fractional-order derivative to the left-hand sides of the first-order differential Equation (1). We use two types of fractional-order derivatives, namely the Caputo operator (that is, the operator with the power law kernel) [39] and Atangana–Baleanu operator [40]. The basic difference between these two operators lies on their kernel. Atangana–Baleanu operator has a non-singular and non-local (that is, Mittag–Leffler) kernel while the Caputo operator does not [41,42]. From the biological meanings, a model with Atangana–Baleanu operator gives better results in applying memory effects [43–45]. Nevertheless, the Caputo operator has more complex analytical tools in investigating the dynamics of the model such as the local stability [46–50], the global stability [51,52], bifurcation theory [52–54], and so on. By considering the deﬁciencies and advantages, the model with Caputo and Atangana–Baleanu operator are employed in our work. Based on our literature review, the dynamics of the model (1) with Caputo and Atangana–Baleanu operator have never been studied. For this reason, we are interested in investigating the dynamical behavior of model (1) using both Caputo and Atangana–Baleanu fractional-order operators. If the ﬁrst-order derivatives at the left hand sides of model (1) are replaced by the dt fractional-order derivatives D , then we obtain x m ˆ xy D x =r ˆx 1 , K a + x (2) n ˆ xy D y = dy H(y). a + x Note that the left hand sides of model (2) have the dimension of (time) , while the parameters at ˆ ˆ the right hand sides of model (2) such as r ˆ, m ˆ , n ˆ , d, and h have the dimension of (time) ; this shows the inconsistency of physical dimensions in the model (2) (see [55,56]). To overcome such inconsistency, we rescale the parameters in the model (2) to get the following model x m ˆ xy a a D x =r ˆ x 1 , K a + x (3) n ˆ xy a a D y = d y H(y), a + x where 0 , if y < T, H(y) = h (y T) , if y T. c + (y T) a a a a a ˆ ˆ ˆ ˆ ˆ By applying new parameters r = r , m = m , n = n , d = d , and h = h , we obtain x mxy D x =rx 1 , K a + x (4) xy D y = dy H(y), a + x where 0 , if y < T, H(y) = h(y T) : , if y T. c + (y T) This paper is organized as follows. In Section 2, we investigate the dynamics of model (4) with the Caputo operator. We identify the existence, uniqueness, non-negativity, and boundedness of solutions. Furthermore, we explore the dynamics of the model by examining the existence of the equilibrium points, their local and global stability, and the existence of Hopf bifurcation. In Section 3, the existence and uniqueness of solutions of the model with Atangana–Baleanu operator are veriﬁed. In Section 4, Axioms 2020, 9, 122 4 of 23 we explore the dynamics of the model using both operators numerically. We demonstrate numerically the stability of the equilibrium point, and the occurrence of forward and Hopf bifurcations. We end our works with a brief conclusion in Section 5. 2. The Caputo Fractional-Order Rosenzweig–MacArthur Model with THP in Predator 2.1. Model Formulation The operator of Caputo fractional-order derivative is deﬁned as follows Deﬁnition 1. [48] Let a 2 (0, 1], f 2 C ([0, +¥),R), and G() is the Gamma function. The Caputo fractional derivative of order-a is deﬁned by C a a 0 D f (t) = (t s) f (s)ds, t 0. (5) G(1 a) The kernel of Caputo operator is known as the power law kernel. By applying Deﬁnition 1 to model (4), we get the Caputo fractional order Rosenzweig–MacArthur model with THP in predator x mxy C a D x =rx 1 F , K a + x (6) nxy C a D y = dy H(y) F . a + x 2.2. Existence and Uniqueness In this part, we study the existence and uniqueness of model (6). Lemma 1. [57] Consider a Caputo fractional-order system C a D x(t) = f (t, x(t)), t > 0, x(0) 0, a 2 (0, 1], (7) n n where f : (0, ¥) Q ! R , Q R . A unique solution of Equation (7) on (0, ¥) Q exists if f (t, x(t)) satisﬁes the locally Lipschitz condition with respect to x. Since the right hand-side of model (6) is a piecewise function which is switched when the number of predators passes through the threshold level, we divide the existence and uniqueness of the solution into two cases, namely y T and y < T. We start from y T. Consider the region Q [0, T ] where Q := (x, y) 2 R : max (jxj,jyj) g, y T , T < +¥, and a mapping F(L) = (F (L), F (L)). + 1 2 For any L = (x y) 2 Q and L = (x ¯, y ¯) 2 Q, we obtain , Axioms 2020, 9, 122 5 of 23 ¯ ¯ ¯ F(L) F(L) = F (L) F (L) + F (L) F (L) 1 1 2 2 x mxy x ¯ mx ¯y ¯ = rx 1 rx ¯ 1 + K a + x K a + x nxy nx ¯y ¯ ¯ ¯ dy H(y) dy H(y) a + x a + x ¯ r xy x ¯y ¯ 2 2 ¯ x x ¯ m = r(x x) + K a + x a + x ¯ xy x ¯y ¯ n d (y y ¯) ( H(y) H(y ¯)) a + x a + x ¯ ¯ ¯ r xy xy =rjx x ¯j + jx + x ¯jjx x ¯j + (m + n) + K a + x a + x ¯ h(y T) h(y T) djy y ¯j + c + (y T) c + (y ¯ T) ¯ ¯ ¯ ¯ r ay(x x) + (ax + xx)(y y) =r x x ¯ + x + x ¯ x x ¯ + (m + n) + j j j jj j K (a + x)(a + x ¯) ch(y y) d y y ¯ + j j (c + y T)(c + y ¯ T) 2gr m + n r x x ¯ + x x ¯ + ay(x x ¯) + (ax ¯ + x ¯ x)(y y ¯) + j j j j j j K a ¯ ¯ djy yj + jy yj 2gr (m + n)g rjx x ¯j + jx x ¯j + jx x ¯j + K a (m + n)(ag + g ) h y y ¯ + d y y ¯ + y y ¯ j j j j j j a c 2gr (m + n)g (m + n)(ag + g ) h = r + + jx x ¯j + + d + jy y ¯j K a a c M L L , 2gr (m + n)g (m + n)(ag + g ) h where M = max r + + , + d + . Hence, F(L) satisﬁes the K a a c Lipschitz condition for y T. By using similar approaches, when y < T, it is easy to check that 2gr (m + n)g (m + n)(ag + g ) ¯ ¯ F(L) F(L) M L L , where M = max r + + , + d 2 2 K a a and hence the Lipschitz condition is also satisﬁed. According to Lemma 1, model (6) with non-negative initial condition has a unique solution L(t) = (x(t), y(t)) 2 Q. Thus, we establish the following theorem. Theorem 1. For each non-negative initial condition (x , y ) 2 Q, there exists a unique solution 0 0 (x(t), y(t)) 2 Q of model (6), which is deﬁned for all t 0. 2.3. Non-Negativity and Boundedness The solution of model (6) is required to be nonnegative and bounded to establish a biologically well-behaved model. To determine the non-negativity and boundedness of the solution of model (6), the following lemmas are needed. C a Lemma 2. [58] Let 0 < a 1. Suppose that f (t) 2 C[a, b] and D f (t) 2 C[a, b]. C a If D f (t) 0,8 t 2 (a, b), then f (t) is a non-decreasing function for each t 2 [a, b]. C a If D f (t) 0,8t 2 (a, b), then f (t) is a non-increasing function for each t 2 [a, b]. t Axioms 2020, 9, 122 6 of 23 Lemma 3. (Standard comparison theorem for Caputo fractional-order derivative [31]). Let x(t) 2 C [0, +¥) . ( ) C a 2 If x(t) satisﬁes D x(t) + lx(t) m, x(0) = x , where a 2 (0, 1], (l, m) 2 R and l 6= 0, then x(t) m m x E [lt ] + . 0 a l l In the following theorem, we prove the non-negativity and boundedness of solutions using the above lemmas. Theorem 2. All solutions of model (6) with non-negative initial conditions are non-negative and uniformly bounded. Proof. We start by proving that, if the initial condition is non-negative, then x(t) 0 for all t ! ¥. Suppose that it is incorrect; then, we can ﬁnd t > 0 such that x(t) > 0, 0 t < t , < 1 x(t ) = 0, (8) : + x(t ) < 0 By employing (8) and the ﬁrst equation of model (6), we obtain C a D x(t ) = 0. (9) t 1 x(t )=0 + + Based on Lemma 2, we have x(t ) = 0, which contradicts the fact that x(t ) < 0. Thus, x(t) 0 1 1 for all t 0. Using a similar procedure, we conclude y(t) 0 for all t > 0. Now, we adopt the similar manner as in [34] to show the boundedness of solutions. By setting up my a function V(t) = x + , we get m dmy C a C a C a D V(t) + dV(t) = D x + D y + dx + t t t n n x mxy m nxy dmy =rx 1 + dy H(y) + dx + K a + x n a + x n rx m H(y) =(d + r)x K n r (d + r)K (d + r) K m H(y) = x + K 2r 4r n (d + r) K 4r According to the standard comparison theorem for Caputo fractional-order derivative in Lemma 3, we achieve the following inequality 2 2 (d + r) K (d + r) K V(t) V(0) E [d(t) ] + , 4r 4r where E is the one-parameter Mittag–Lefﬂer function. Since E [d(t) ] ! 0 as t ! 0, we acquire a a (d + r) K V(t) ! for t ! ¥. Hence, with non-negative initial condition, all solutions are restricted to 4r the region Q where my (d + r) K Q := (x, y) 2 R : x + + #, # > 0 . n 4r Consequently, all solutions of model (6) are uniformly bounded. Axioms 2020, 9, 122 7 of 23 2.4. The Existence of Equilibrium Point The ﬁrst commonplace technique in studying the dynamical behavior of a fractional-order system is investigating the existence of the equilibrium point. Due to the biological nature, we give the following deﬁnition. Deﬁnition 2. Consider a Caputo fractional-order system C a D ~ x = f (~ x);~ x(0) = ~ x ; a 2 (0, 1]. (10) A point ~ x is called an equilibrium point of Equation (10) if it satisﬁes f (~ x ) = 0. Particularly, it is called the biological equilibrium point of Equation (10) if it satisﬁes ~ x 0. Based on Deﬁnition 2, the equilibrium point of model (6) is obtained by solving rx my r x =0, K a + x (11) nxy dy H(y) =0. a + x Thus, we get four feasible biological equilibrium points as follows. (i) The equilibrium points when y < T are (i.a) the origin point E = (0, 0) which always exists; (i.b) the predator extinction point E = (K, 0) which always exists; and (K x ˆ) (a + x ˆ)r ad (i.c) the ﬁrst co-existence point E = x ˆ, , with x ˆ = , which exists if mK n d ad TmK n > + d and (K x ˆ) (a + x ˆ) < . K r (ii) The equilibrium point when y T is the second co-existence point E = (x , y ) where y = (K x ) (a + x )r 4 3 2 and x is the positive roots of polynomial b x + b x + b x + b x + b = 0 1 2 3 4 5 mK where b =(n d)r , b = [(an + 2dK) 2(ad + nK)] r , b =(nrK + 4adr + cdm + mnT + hm)rK ((drK + 2anr + cmn + dmT)K + a dr)r, b =((anr + cmn + dmT)K + (2adr + hm + cdm)a)rK ((2adr + hm + cdm + mnT)K + admT)rK, b = (adr + hm)mT (adr + hm + cdm)ar K . [ ] TmK E exists if 0 < x < K and (K x ) (a + x ) . 2.5. Local Asymptotic Stability In this part, we discuss the local stability of E , E , E, and E . For this aim, we need the 0 1 following theorem. Theorem 3. (Matignon condition [48,59]) The equilibrium point ~ x of system (10) is locally asymptotically stable if all eigenvalues l of the Jacobian matrix J = ¶ f /¶~ x evaluated at ~ x satisfy j arg(l )j > ap/2. j j If there exists at least one eigenvalue satisfying j arg(l )j > ap/2 and j arg(l )j < ap/2, k 6= l, then ~ x is k l a saddle-point. Axioms 2020, 9, 122 8 of 23 Now, we present Theorems 4–7 to show the local stability properties of E , E , E, and E . 0 1 Theorem 4. The origin point E = (0, 0) is always a saddle point. " # r 0 Proof. When E = (0, 0), the model (6) has Jacobian matrix J(E ) = , where its eigenvalues 0 0 0 d are l = r > 0 and l = d < 0. It is clear that jarg(l )j = 0 < ap/2 and jarg(l )j = p > ap/2. 1 2 1 2 Therefore, based on Theorem 3, E is always a saddle point. ad Theorem 5. The predator extinction point E = (K, 0) is locally asymptotically stable if n < + d. ad Otherwise, if n > + d, then E = (K, 0) is a saddle point. " # mK a+K Proof. The Jacobian matrix of model (6) evaluated at E is J(E ) = . The eigenvalues 1 1 nK 0 d a+K nK of J(E ) are l = r < 0 and l = d. Clearly, jarg(l )j = p > ap/2 and jarg(l )j = p > 1 1 2 1 2 a + K ad ad ap/2 if n < + d and arg(l ) = 0 < ap/2 if n > + d. Hence, we have the theorem. j j K K Remark 1. It is noted that the existence condition for the ﬁrst co-existence point E contradicts the stability condition of E . Consequently, if E is locally asymptotically stable, then E does not exist. This condition also 1 1 indicates the existence of forward bifurcation, which is conﬁrmed numerically in the next section. 2 2 ˆ ˆ ˆ ˆ ˆ 4 (K x) anrx K x r x 2 D(a + x)K Theorem 6. Let D = 1 and a ˆ = tan . 2 2 (a + x ˆ) K a + x ˆ K p (K a 2x ˆ) rx ˆ Suppose that one of the following statements is obeyed. K a (i) x ˆ > ; or K a (ii) x ˆ < , D > 0, and a < a ˆ . ˆ ˆ (K x) (a + x)r Then, the ﬁrst co-existence point E = x ˆ, is locally asymptotically stable. mK Proof. We ﬁrst observe that the Jacobian matrix of model (6) evaluated at E is 2 3 ˆ ˆ ˆ K x rx mx 6 7 a + x ˆ K a + x ˆ ˆ 6 7 J(E) = . (12) 4 5 (K x ˆ) anr (a + x)mK The eigenvalues of the Jacobian matrix (12) are the solutions of the characteristic equation K x ˆ rx ˆ K x ˆ anrx ˆ ( ) l 1 l + = 0, a + x ˆ K (a + x ˆ) K namely K x ˆ rx ˆ i D l = 1 + , a + x ˆ 2K 2 (13) ˆ ˆ K x rx i D l = 1 . a + x ˆ 2K 2 Axioms 2020, 9, 122 9 of 23 K a When x ˆ > , the real parts of l are negative. Thus, the eigenvalues (13) always satisfy 1,2 K a (K x ˆ) anrx ˆ arg(l ) > ap/2 for any D. If x ˆ < and D 0, then l l = > 0 and l + j j 1,2 1 2 1 2 (a + x ˆ) K K x ˆ rx ˆ K a l = 1 > 0, meaning that arg(l ) = 0 < ap/2. If x ˆ < and D > 0, j j 2 1,2 a + x ˆ K 2 then jarg(l )j > ap/2 for a < a ˆ . Hence, we prove the theorem. 1,2 Theorem 7. Let ((y T) cT)h my r x = + x , 2 2 y (c + y T) (a + x ) K my r ((y T) cT)h x mnx y q = x + 1 . 2 2 2 (a + x ) K y (c + y T) a + x (a + x ) If one of the following statements is satisﬁed,then the second co-existence point E = (x , y ) is locally asymptotically stable: (i) q > 0 and x < 0; or (ii) x < 4q, x > 0, and a < a . Proof. The Jacobian matrix of model (6) evaluated at E is given by 2 3 my r mx 6 7 K a + x (a + x ) 6 7 J(E ) = . (14) 6 7 4 5 x ny ((y T) cT)h a + x a + x y (c + y T) The eigenvalues of (14) are obtained by solving the characteristic equation l xl + q = 0. Thus, x x 4q we have l = . If q > 0 and x < 0, then jarg(l )j > ap/2. If x < 4q and x > 0, 1,2 1,2 2 2 thenjarg(l )j > ap/2 for a < a . Using Theorem 3, the local stability of E is completely proven. 1,2 2.6. Global Asymptotic stability To study the global stability of equilibrium points, we need the following lemmas. Lemma 4. [51] Let x(t) 2 C (R ), x 2 R , and its Caputo fractional derivative of order-a exists for any + + x(t) x a C a C a 2 (0, 1]. Then, for any t > 0, we have D x(t) x x ln 1 D x(t). t t x x(t) Lemma 5. (Generalized Lasalle Invariance Principle [52]). Suppose W is a bounded closed set and every solution of system C a D x(t) = f (x(t)), (15) which starts from a point in W remains in W for all time. Let V(x) : W ! R be a continuously differentiable function such that C a D Vj 0. (15) C a Let E := xj D Vj = 0 and M be the largest invariant set of E. Then, every solution x(t) originating (15) in W tends to M as t ! ¥. Since model (6) is basically a piecewise fractional-order differential equations that depends on T, the analysis of the global stability is split into two regions deﬁned by W := f(x, y) : x 0, y < Tg and W := (x, y) : x 0, y T . Therefore, the global stabilities of E , E, and E are investigated f g 2 1 as follows. Axioms 2020, 9, 122 10 of 23 ad Theorem 8. If n < , then the predator extinction point E = (K, 0) is globally asymptotically stable in the region W . x my Proof. By specifying a positive Lyapunov function L (x, y) = x K K ln + , K n and conforming to Lemma 4, we obtain x K m C a C a C a D L (x, y) D x + D y t t t x n rx my m nxy =(x K) r + dy K a + x n a + x rx mKy dmy =2rx rK + K a + x n r mKy dmy = (x K) + K a + x n r mKy dmy (x K) + K a n d K my. n a ad C a Owing to the fact that n < , we have D L (x, y) 0. In consequence of Lemma 5, E is 1 1 globally asymptotically stable in the region W . ad Remark 2. Notice E is locally asymptotically stable if n < + d and is globally asymptotically stable if ad n < . Hence, if the global stability condition is fulﬁlled, then the local stability is also achieved but not vice versa. This fact reinforces that the global stability condition has a larger attracting region than that of the local stability condition. ˆ ˆ (K x) (a + x)r Theorem 9. If (K x ˆ) (a + x ˆ) < a , then the ﬁrst co-existence point E = x ˆ, is globally mK asymptotically stable in the region W . K x ˆ (a + x ˆ)r ( ) Proof. Let E = (x ˆ, y ˆ) where y ˆ = and j is a positive constant. By considering a mK h i x y ˆ ˆ ˆ ˆ Lyapunov function L (x, y) = x x x ln + j y y y ln , and using Lemma 4, we get x ˆ y ˆ x x ˆ y y ˆ C a C a C a D L (x, y) D x + j D y t t t x y rx my nx =(x x ˆ) r + j(y y ˆ) d K a + x a + x rx ˆ my ˆ rx my nx nx ˆ =(x x ˆ) + + j(y y ˆ) K a + x ˆ K a + x a + x a + x ˆ r y y ˆ x x ˆ = (x x ˆ) (x x ˆ) + m + nj(y y ˆ) K a + x a + x ˆ a + x a + x ˆ r (y y ˆ)(a + x ˆ) + (x ˆ x)y ˆ (x x ˆ)(y y ˆ) = (x x ˆ) (x x ˆ) m + anj K (a + x)(a + x ˆ) (a + x)(a + x ˆ) ˆ ˆ ˆ ˆ ˆ r (x x)(y y)(a + x) (x x) y = (x x ˆ) m + m+ K (a + x)(a + x ˆ) (a + x)(a + x ˆ) (x x ˆ)(y y ˆ) anj (a + x)(a + x ˆ) r my ˆ (x x ˆ)(y y ˆ) (x x ˆ) + (anj (a + x ˆ)m). K a (a + x)(a + x ˆ) Axioms 2020, 9, 122 11 of 23 (a + x ˆ)m By choosing j = and substituting the value of y ˆ, we get an C a 2 D L (x, y) (x x ˆ) a (K x ˆ) (a + x ˆ) . t 2 a K 2 C a Therefore, if (K x ˆ) (a + x ˆ) < a , then D L (x, y) 0. Using Lemma 5, we conclude that E is globally asymptotically stable in the region W . Based on Theorem 2, x(t) and y(t) are bounded. Let Y be the upper bound of y(t) such that 0 < T y(t) Y. The global stability of E is stated in the following theorem. a r cT Theorem 10. If y < min , + T , then the second co-existence point E = (x , y ) is globally mK Y T asymptotically stable in the region W . h i Proof. We consider a positive Lyapunov function L (E ) = x x x ln + j y y y ln . According to Lemma 4, the fractional derivative of L (E ) satisﬁes x x y y C a C a C a D L (x, y) D x + j D y t t t x y rx my nx h(y T) =(x x ) r + j (y y ) d K a + x a + x (c + y T)y rx my rx my =(x x ) + K a + x K a + x nx nx h(y T) h(y T) + j (y y ) + a + x a + x (c + y T)y (c + y T)y r (y y )(a + x )m (x x )my = (x x ) (x x ) K (a + x )(a + x) (x x )an (y y )chT (y y )(y T)(y T)h + j (y y ) (a + x )(a + x) (c + y T)(c + y T)y y r (x x )(y y )(a + x )m (x x ) my = (x x ) + K (a + x )(a + x) (a + x )(a + x) (x x )(y y )anj (cT (y T)(y T))(y y ) hj (a + x )(a + x) (c + y T)(c + y T)y y r (x x )(y y )(a + x )m my 2 2 (x x ) + (x x ) K (a + x )(a + x) a (x x )(y y )anj (cT (y T)(y T))(y y ) hj (a + x )(a + x) (c + y T)(c + y T)y y (a + x )m By taking j = and remembering that y(t) < Y, we obtain an r my (cT (y T)(Y T))(y y ) (a + x )mh C a 2 D L (E ) (x x ) . K a (c + y T)(c + y T)any y a r cT C a It is easily conﬁrmed that, if y < min , + T , then D L (x, y) 0. Based on mK Y T Lemma 5, E is globally asymptotically stable in the region W . 2 Axioms 2020, 9, 122 12 of 23 2.7. The Existence of Hopf Bifurcation The Hopf bifurcation is a local phenomenon when a stable equilibrium point loses its stability and all nearby solutions converge to a periodic solution namely limit-cycle if a parameter is varied [54,60,61]. It is shown that many fractional-order models involving the Caputo operator undergo a Hopf bifurcation which is driven by the order of the derivative (see [2,17,34,53,62]). The difference between the Hopf bifurcation in the integer-order model and that in the fractional-order model lies on the property of the limit-cycle. In the integer-order model, the limit-cycle is a periodic orbit which does not exist in the fractional-order model [63]. In the fractional-order model, the limit-cycle is not a periodic solution, but all nearby solutions converge to a limit-cycle [56,62]. Adapted from Theorem 3 in [62], for two dimensional Caputo fractional-order system, a Hopf bifurcation occurs when the eigenvalues l of the Jacobian matrix evaluated at the equilibrium point 1,2 satisfy the following conditions: (i) l = y wi where y > 0; 1,2 (ii) m(a ) = a p/2 min jarg(l )j = 0; and 1i2 i dm(a) (iii) 6= 0. da a=a Now, consider the stability condition in Theorems 6 and 7. For y < T, the Jacobian matrix of model (6) evaluated at E has a pair of complex eigenvalues if D > 0. We can easily conﬁrm that, K a if x ˆ < , then the real part of the eigenvalues are positive. We also have that m(a ˆ) = 0 and dm(a) 6= 0. Therefore, E undergoes a Hopf bifurcation when a crosses a ˆ . A similar circumstance da a=a ˆ also occurs when y T. When x < 4q, the Jacobian matrix of model (6) has a pair of complex eigenvalues. The real part of the eigenvalues are also positive when x > 0. We can also check that dm(a) m(a ) = 0 and 6= 0. This means the Hopf bifurcation also occurs when y T. Therefore, da a=a we have the following theorem. K a Theorem 11. (i) Let D > 0 and x ˆ < . The ﬁrst co-existence point E undergoes a Hopf bifurcation when a passes through a ˆ in the region W . (ii) Let x < 4q and x > 0. The second co-existence point E undergoes a Hopf bifurcation when a passes through a in the region W . 3. The Atangana–Baleanu Fractional-Order Rosenzweig–MacArthur Model with THP in Predator 3.1. Model Formulation In this section, we consider a fractional-order Rosenzweig–MacArthur model with THP in predator involving the Atangana–Baleanu operator. Speciﬁcally, we consider the Atangana–Baleanu operator in Caputo sense of order-a which is deﬁned as follows. Deﬁnition 3. [40] Suppose 0 < a 1. The Atangana–Baleanu fractional integral and derivative in Caputo sense of order-a (ABC) is deﬁned by 1 a a ABC a a1 I f (t) = f (t) + (t s) f (s)ds, B(a) G(a)B(a) B(a) a ABC a a 0 D f (t) = E (t s) f (s)ds, 1 a 1 a n ¥ where t 0, f 2 C ([0, +¥),R), E (t) = å is the Mittag–Lefﬂer function and B(a) is a k=0 G(ak + 1) normalization function with B(0) = B(1) = 1. In this paper, we deﬁne B(a) = 1 a + . G(a) Axioms 2020, 9, 122 13 of 23 By applying Deﬁnition 3 to model (4), we get the following fractional-order model with ABC operator x mxy ABC a D x =rx 1 G , K a + x (16) nxy ABC a D y = dy H(y) G . t 2 a + x Due to the lack of analytical theory, model (16) is investigated numerically. However, we ﬁrst show the existence and uniqueness of the solution of the model (16). 3.2. Existence and Uniqueness We start this work by representing the Lipschitz condition of the kernels of model (16). Since the harvesting is performed by obeying threshold harvesting policy, we give the proof into two cases i.e., y < T and y T. We start for y T. Let x, x ¯, y, and y ¯ be functions satisfying kxk a , kx ¯k a , kyk b , 1 1 and ky ¯k b . Suppose that r my g =r + (a + a ) + , 1 1 2 K a g =n + d + . Therefore, we get x mxy x ¯ mx ¯y kG (x) G (x ¯)k = rx 1 rx ¯ 1 1 1 K a + x K a + x 2 2 rx mxy rx ¯ mx ¯y = rx rx ¯ + + K a + x K a + x ¯ r x x ¯ 2 2 ¯ ¯ = r(x x) (x x ) my K a + x a + x ¯ (17) r amy(x x ¯) ¯ ¯ ¯ = r(x x) (x + x)(x x) K (a + x)(a + x ¯) r my ¯ ¯ ¯ rkx xk + (a + a ) kx xk + kx xk 1 2 K a r my = r + (a + a ) + kx x ¯k 1 2 K a =g x x ¯ , k k and nxy h(y T) nxy ¯ h(y ¯ T) kG (y) G (y ¯)k = dy dy ¯ 2 2 a + x c + (y T) a + x c + (y T) nxy h(y T) nxy ¯ h(y ¯ T) = dy + dy ¯ + a + x c + (y T) a + x c + (y ¯ T) nx(y y ¯) ch(y y ¯) = d(y y ¯) a + x (c + y T)(c + y ¯ T) (18) ¯ ¯ ¯ nky yk + dky yk + ky yk = n + d + ky y ¯k =g ky y ¯k . 2 Axioms 2020, 9, 122 14 of 23 When y < T, by utilizing the similar manner, we achieve G (y) G (y ¯) g y y ¯ where k k k k 2 2 3 g = n + d. Accordingly, the kernel of (16) satisﬁes the Lipschitz condition. Furthermore, if 0 g < 1 3 1 and 0 g < g < 1, then G and G are contracted. 3 2 1 2 Now, by employing the ﬁxed-point theorem, the solution of model (16) is investigated. By utilizing the Atangana–Baleanu fractional integral operator, model (16) is transformed into the following Volterra-type integral equation. 1 a a a1 x(t) x(0) = G (t, x) + (t s) G (s, x)ds, 1 1 B(a) B(a)G(a) (19) 1 a a a1 y(t) y(0) = G (t, y) + (t s) G (s, y)ds. 2 2 B(a) B(a)G(a) In a recursive form, Equation (19) is written as 1 a a a1 x (t) = G (t, x ) + (t s) G (s, x )ds, 1 n1 1 n1 B(a) B(a)G(a) (20) 1 a a a1 y (t) = G (t, y ) + (t s) G (s, y )ds. n 2 n1 2 n1 B(a) B(a)G(a) The associated initial conditions along with Equation (20) are x (t) = x(0) and y (t) = y(0). 0 0 By considering Equation (20), we have the difference expression of successive terms as follows. F (t) =x (t) x (t) 1,n n1 1 a a a1 = (G (t, x ) G (t, x )) + (t s) (G (t, x ) G (t, x )) ds, 1 n1 1 n2 1 n1 1 n2 B(a) B(a)G(a) (21) F (t) =y (t) y (t) 2,n n n1 1 a a a1 = (G (t, y ) G (t, y )) + (t s) (G (t, y ) G (t, y )) ds. 2 n1 2 n2 2 n1 2 n2 B(a) B(a)G(a) 0 n n According to Equation (21), we get x (t) = F (t) and y (t) = F (t). By applying n å n å 1,i 2,i i=1 i=1 Equations (17), (18) and (21), we have that 1 a a a1 kF (t)k g kF k + g kF (s)k (t s) ds, 1,n 1 1,n1 1 1,n1 B(a) B(a)G(a) (22) 1 a a a1 kF (t)k g kF k + g kF (s)k (t s) ds. 2,n 2 2,n1 2 2,n1 B(a) B(a)G(a) Therefore, by using (22), the existence and uniqueness of model (16) is presented as follows. Theorem 12. Model (16) has a unique solution if we can ﬁnd t such that t g (1 a)g + < 1, i = 1, 2, 3. (23) B(a) B(a)G(a) Proof. Let x(t) and y(t) be bounded functions, and hence the Lipschitz condition is satisﬁed. Thus, according to Equation (22), we obtain the following inequalities. (1 a)g t g 1 1 kF (t)k kx k + , 1,n B(a) B(a)G(a) (24) (1 a)g t g 2 2 F (t) y + . k k k k 2,n 0 B(a) B(a)G(a) Axioms 2020, 9, 122 15 of 23 Therefore, the continuity and existence of solution are proved since F (t) ! 0 and k k 1,n kF (t)k ! 0 as n ! ¥ and t = t . Now, suppose that 2,n 0 x(t) x(0) =x (t) U (t), n 1,n (25) y(t) y(0) =y (t) U (t). n 2,n We conﬁrm that n+1 (1 a) t n+1 kU (t)k + g (26) 1,n B(a) B(a)G(a) It is clear that kU (t)k ! 0 when n ! ¥. By using a similar manner, we acquire 1,n kU (t)k ! 0, n ! ¥. Finally, the uniqueness of solution for the model is proven. Suppose that there 2,n exists different set of solutions denote by x ˜(t) and y ˜(t); then, we have 1 a a a1 x(t) x ˜(t) = (G (t, x) G (t, x ˜)) + (G (s, x) G (s, x ˜))(t s) ds. (27) 1 1 1 1 B(a) B(a)G(a) Taking the norm for both sides and using a simpliﬁcation as in (22) and (24), we obtain (1 a)g t g 1 1 kx(t) x ˜(t)k 1 0. (28) B(a) B(a)G(a) For t = t , we have (23), thus kx(t) x ˜(t)k = 0 and hence x(t) = x ˜(t). Applying the same algebraic procedures, we can show that y(t) = y ˜(t). Therefore, the solution is unique. 4. Numerical Simulations In this section, the numerical simulations of Caputo model (6) and ABC model (16) are performed to illustrate the previous theoretical results. In the literature, there exist many numerical methods to solve a system of fractional differential equations, such as the Monte Carlo method [64], the Grünwald–Letnikov method [65,66] and the predictor–corrector method [67–69]. We apply the predictor–corrector scheme proposed by Diethelm et al. [67] to solve the Caputo fractional-order model and the predictor–corrector scheme proposed by Baleanu et al. [69] to solve the Atangana–Baleanu in Caputo sense model (ABC). Due to the limitation of ﬁeld data, we use hypothetical parameter values that correspond to the theoretical results. The initial parameter values are given as follows. r = 0.5, K = 1, m = 0.3, a = 0.2, d = 0.1, h = 0.1, T = 0.9, c = 0.1, anda = 0.9. (29) In Figure 1, we plot a bifurcation diagram by varying the conversion rate of consumed prey into predator birth n in interval [0.08, 0.2]. We notice that the parameter values (29) and the interval 0.08 n 0.2 lead to the non-existence of equilibrium point in W . Therefore, the ﬁrst numerical simulations are focused on the dynamics in W . For 0.08 n < n = 0.12, Theorem 5 states that the predator extinction point E = (1, 0) is the only equilibrium point which is asymptotically stable. To see this behavior, we take n = 0.1 and plot the phase portrait and the time series as in Figure 2. It is seen that all solutions with initial values in both W and W are convergent to E . When the initial value is in W , then the solution is oscillating 1 2 1 2 when it crosses the threshold harvesting level and eventually converges to E . When n passes through n , the equilibrium point E = (1, 0) undergoes a forward bifurcation. In this case, there appear two equilibrium points, namely the unstable E and an asymptotically stable ˆ ˆ E. Figure 1 shows that E is asymptotically stable if 0.12 < n . n = 0.1557. In Figure 3, we show the phase portrait and time series for the case of n = 0.14. We see that E = (0, 0) and E = (1, 0) 0 1 are saddle points, while E = (0.5, 0.5833) is asymptotically stable. This circumstance corresponds to Theorems 4–6 and 9. Axioms 2020, 9, 122 16 of 23 1.0 Caputo Upper bound of limit-cycle 0.8 ABC n = 0.14 E-stable E-unstable 0.6 n = 0.2 0.4 Lower bound of limit-cycle 0.2 E1-stable E1-unstable n = 0.1 0.0 1.2 E -stable E -unstable 1 1 1.0 n = 0.1 Upper bound of limit-cycle 0.8 0.6 0.4 n = 0.14 E-stable E-unstable 0.2 Caputo n = 0.2 Lower bound of limit-cycle ABC 0.0 0.08 0.10 0.12 0.14 0.16 0.18 0.20 Figure 1. Bifurcation diagram of Caputo model (6) and ABC model (16) driven by the conversion rate of consumed prey into predator birth (n) with constant parameter values (29). There exists two bifurcations namely a forward bifurcation which occurs when n passes through n 0.12, and a Hopf bifurcation which occurs when n passes through n 0.1557. 1.2 1.0 1.00 Caputo Caputo 0.95 0.8 0.90 ABC ABC 0.85 THP 1.0 0.80 2 3 4 5 6 7 8 0.6 0.4 0.8 0.2 0.0 0.6 1.2 1.0 0.4 0.8 0.2 0.6 Caputo 0.4 ABC 0.0 0.2 0.0 0.5 1.0 1.5 2.0 0 250 500 750 1000 1250 1500 1750 2000 x(t) t (a) (b) Figure 2. Numerical simulations of Caputo model (6) and ABC model (16) with parameter values (29) and n = 0.1: Figure (a) phase portrait; and Figure (b) time series. Furthermore, if we increase the value of n such that n > n , then the equilibrium point E loses its stability, and all solutions converge to a limit-cycle. This situation conﬁrms the occurrence of a Hopf bifurcation driven by n. For example, if we select n = 0.2 then all equilibrium points are unstable and all solutions eventually converge to limit cycle (see Figure 4). Now, we perform simulation using the same parameter values as in Figure 4, but with a lower threshold value, namely T = 0.5. In this case, there is no equilibrium point E in W , and E = (0.5954, 0.5364) occurs in W . According to Theorem 7, E is asymptotically stable. Such dynamics can be clearly seen in Figure 5. This simulation shows that by applying the THP when the interior equilibrium point is stable, we can choose a suitable constant of threshold so that the existence of both prey and predator are maintained. y(t) x(t) y(t) x(t) y(t) y(t) Axioms 2020, 9, 122 17 of 23 1.2 Caputo Caputo 0.8 ABC ABC THP 0.7 1.0 0.6 0.8 0.5 0.4 0.6 Caputo 1.0 ABC 0.4 0.8 0.6 0.2 0.4 0.2 0.0 E 0.00 0.25 0.50 0.75 1.00 1.25 1.50 0 100 200 300 400 500 x(t) t (a) (b) Figure 3. Numerical simulations of Caputo model (6) and ABC model (16) with parameter values (29) and n = 0.14: Figure (a) phase portrait; and Figure (b) time series. Notice that, in Figures 2–5, we see that both model with Caputo operator (6) and Atangana–Baleanu operator (16) have similar dynamical behavior. The noticeable difference between them is the orbit of solutions and the diameter of the limit-cycle. In Figure 4, the diameter of the limit-cycle of the model with ABC operator looks shorter than that of the Caputo operator, which may give different dynamics when a Hopf bifurcation occurs. To get more detail view, we plot a bifurcation diagram by varying the order of the fractional derivative (a) (see Figure 6). In this simulation, we use parameter values as in Figure 4 and vary the order-a in the interval [0.6, 0.9]. It is clearly seen that, besides the diameter of the limit-cycle, the bifurcation points of Caputo model and ABC model are also different. The Caputo model has an earlier bifurcation point than that of the ABC model. To show this situation, we show some numerical simulations using different values of a (see Figure 7). For a = 0.7, the equilibrium point E of both model are asymptotically stable. For a = 0.772, the equilibrium point E of the Caputo model loses its stability, while that of the ABC model remains asymtotically stable. For a = 0.83, the equilibrium point E of both models are unstable. 1.2 1.0 Caputo LC-ABC Caputo 0.8 ABC THP ABC LC-Caputo 1.0 0.6 0.4 0.8 0.2 0.0 0.6 1.0 Caputo 0.8 0.4 ABC 0.6 0.2 0.4 E 0.2 0.0 0.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0 250 500 750 1000 1250 1500 1750 2000 x(t) t (a) (b) Figure 4. Numerical simulations of Caputo model (6) and ABC model (16) with parameter values (29) and n = 0.2: Figure (a) phase portrait; and Figure (b) time series. y(t) y(t) x(t) y(t) x(t) y(t) Axioms 2020, 9, 122 18 of 23 1.1 1.0 Caputo ABC THP Caputo ABC 1.0 0.8 0.9 0.6 0.8 0.4 0.7 0.8 0.6 * 0.6 0.5 0.4 0.4 0.2 Caputo ABC 0.3 0.0 0.0 0.2 0.4 0.6 0.8 1.0 0 100 200 300 400 500 x(t) t (a) Phase portrait (b) n = 0.14 Figure 5. Numerical simulations of Caputo model (6) and ABC model (16) with parameter values (29), n = 0.2 and T = 0.5. Caputo 0.8 Upper bound of limit-cycle ABC 0.7 The difference in stability of Caputo and ABC 0.6 E-stable = 0.7 = 0.772 = 0.83 E-unstable 0.5 0.4 Lower bound of limit-cycle 0.3 Caputo 0.6 Upper bound of limit-cycle ABC 0.5 The difference in stability of Caputo and ABC 0.4 0.3 E-stable = 0.7 = 0.772 = 0.83 E-unstable 0.2 0.1 Lower bound of limit-cycle 0.60 0.65 0.70 0.75 0.80 0.85 0.90 Figure 6. Bifurcation diagram of Caputo model (6) and ABC model (16) driven by the order of the fractional-derivative (a) with constant parameter values (29) and n = 0.2. There exists a Hopf bifurcation where the bifurcation points of the Caputo model and ABC model are different. 0.8 0.6 E E Caputo Caputo 0.4 = 0.772 = 0.83 Caputo = 0.7 Limit-cycle Limit-cycle 0.2 0.8 0.6 ABC 0.4 = 0.83 ABC ABC = 0.7 = 0.772 Limit-cycle 0.2 0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 x(t) x(t) x(t) (a) Figure 7. Cont. y(t) y(t) y(t) x(t) y(t) x(t) y(t) Axioms 2020, 9, 122 19 of 23 0.8 0.6 0.4 Caputo Caputo Caputo = 0.7 = 0.772 = 0.83 0.2 x(t) x(t) x(t) y(t) y(t) y(t) 0.0 0.8 0.6 0.4 ABC ABC ABC = 0.7 = 0.772 = 0.83 0.2 x(t) x(t) x(t) y(t) y(t) y(t) 0.0 0 500 1000 1500 2000 0 500 1000 1500 2000 0 500 1000 1500 2000 t t t (b) Figure 7. Numerical simulations of Caputo model (6) and ABC model (16) with parameter values (29), n = 0.2, and a = f0.7, 0.772, 0.83g: Figure (a) phase portrait; and Figure (b) time series. From the biological point of view, all previous numerical simulations show that the dynamical properties of both Caputo model and ABC model are similar except when the eigenvalues of the Jacobian matrix evaluated at the interior equilibrium point E are a pair of complex conjugate with positive real part. There is a biological condition such that the prey and predator densities are eventually periodic for the Caputo model, while for ABC model, the densities of both predator and prey are eventually constant. 5. Conclusions The dynamics of a Rosenzweig–MacArthur model with continuous threshold harvesting in predator involving the Caputo fractional-order derivative and ABC fractional-order derivative are studied. We prove the existence and uniqueness of solutions of both Caputo and ABC models. Particularly, we completely investigate the dynamics of the Caputo model including the non-negativity, boundedness, local stability, global stability, and the existence of Hopf bifurcation. From the biological meanings, the extinction of both populations never occurs since the origin point (E ) is a saddle point. Some of the situations that might occur are: (1) the predator goes extinct while the prey still survives, which is indicated by the stability of E ; (2) both predator and prey co-exist and converge to a constant population density, which happens if the interior point E or E are asymptotically stable; and (3) both predator and prey co-exist where both population densities change periodically, namely when a Hopf bifurcation occurs. We show numerically that our model may undergo a forward bifurcation or a Hopf bifurcation. The Hopf bifurcation in models with both Caputo operator and ABC operator can be driven by the conversion rate of consumed prey into the predator birth rate or by the order of fractional derivative. Our numerical simulations show that the Hopf bifurcation point of both models are different. Author Contributions: All authors equally contributed to this article. All authors read and approved the ﬁnal manuscript. Funding: This research was funded by the Directorate of Research and Community Service, The Directorate General of Strengthening Research and Development, the Ministry of Research, Technology, and Higher Education (Brawijaya University), Indonesia, via Doctoral Dissertation Research, in accordance with the Research Contract No. 037/SP2H/LT/DRPM/2020, dated 9 March 2020. Acknowledgments: The authors would like to thank the reviewers for all useful and helpful comments to improve the manuscript. Conﬂicts of Interest: All authors report no conﬂict of interest relevant to this article. x(t), y(t) x(t), y(t) Axioms 2020, 9, 122 20 of 23 References 1. Rosenzweig, M.L.; MacArthur, R.H. Graphical representation and stability conditions of predator-prey interactions. Am. Nat. 1963, 97, 209–223. doi:10.1086/282272. 2. Moustafa, M.; Mohd, M.H.; Ismail, A.I.; Abdullah, F.A. Stage structure and refuge effects in the dynamical analysis of a fractional order Rosenzweig-MacArthur prey-predator model. Prog. Fract. Differ. Appl. 2019, 5, 49–64. doi:10.18576/pfda/050106. 3. Beay, L.K.; Suryanto, A.; Darti, I.; Trisilowati, T. Hopf bifurcation and stability analysis of the Rosenzweig-MacArthur predator-prey model with stage-structure in prey. Math. Biosci. Eng. 2020, 17, 4080–4097. doi:10.3934/mbe.2020226. 4. González-Olivares, E.; Ramos-Jiliberto, R. Dynamic consequences of prey refuges in a simple model system: More prey, fewer predators and enhanced stability. Ecol. Model. 2003, 166, 135–146. doi:10.1016/S0304-3800(03)00131-5. 5. Chen, L.; Chen, F.; Chen, L. Qualitative analysis of a predator–prey model with Holling type II functional response incorporating a constant prey refuge. Nonlinear Anal. Real World Appl. 2010, 11, 246–252. doi:10.1016/j.nonrwa.2008.10.056. 6. Almanza-Vasquez, E.; Ortiz-Ortiz, R.D.; Marin-Ramirez, A.M. Bifurcations in the dynamics of Rosenzweig-Macarthur predator-prey model considering saturated refuge for the preys. Appl. Math. Sci. 2015, 9, 7475–7482. doi:10.12988/ams.2015.510640. 7. Das, M.; Maiti, A.; Samanta, G.P. Stability analysis of a prey-predator fractional order model incorporating prey refuge. Ecol. Genet. Genom. 2018, 7-8, 33–46. doi:10.1016/j.egg.2018.05.001. 8. Moustafa, M.; Mohd, M.H.; Ismail, A.I.; Abdullah, F.A. Dynamical analysis of a fractional-order Rosenzweig–MacArthur model incorporating a prey refuge. Chaos Solitons Fractals 2018, 109, 1–13. doi:10.1016/j.chaos.2018.02.008. 9. Sarkar, K.; Khajanchi, S. Impact of fear effect on the growth of prey in a predator-prey interaction model. Ecol. Complex. 2020, 42, 100826. doi:10.1016/j.ecocom.2020.100826. 10. Zu, J.; Mimura, M. The impact of Allee effect on a predator-prey system with Holling type II functional response. Appl. Math. Comput. 2010, 217, 3542–3556. doi:10.1016/j.amc.2010.09.029. 11. Pal, P.J.; Saha, T. Qualitative analysis of a predator-prey system with double Allee effect in prey. Chaos Solitons Fractals 2015, 73, 36–63. doi:10.1016/j.chaos.2014.12.007. 12. Bodine, E.N.; Yust, A.E. Predator–prey dynamics with intraspeciﬁc competition and an Allee effect in the predator population. Lett. Biomath. 2017, 4, 23–38. doi:10.1080/23737867.2017.1282843. 13. Mukherjee, D. Role of fear in predator–prey system with intraspeciﬁc competition. Math. Comput. Simul. 2020, 177, 263–275. doi:10.1016/j.matcom.2020.04.025. 14. Van Den Bosch, F. Cannibalism in an age-structured predator-prey system. Bull. Math. Biol. 1997, 59, 551–567. doi:10.1016/s0092-8240(96)00107-3. 15. Mondal, S.; Lahiri, A.; Bairagi, N. Analysis of a fractional order eco-epidemiological model with prey infection and type 2 functional response. Math. Methods Appl. Sci. 2017, 40, 6776–6789. doi:10.1002/mma.4490. 16. Suryanto, A.; Darti, I.; Anam, S. Stability analysis of pest-predator interaction model with infectious disease in prey. AIP Conf. Proc. 2018, 1937, 020018. doi:10.1063/1.5026090. 17. Panigoro, H.S.; Suryanto, A.; Kusumawinahyu, W.M.; Darti, I. Dynamics of a fractional-order predator-prey model with infectious diseases in prey. Commun. Biomath. Sci. 2019, 2, 105–117. doi:10.5614/cbms.2019.2.2.4. 18. Kumar, T.; And, K.; Chakraborty, K. Effort dynamics in a prey-predator model with harvesting. Int. J. Inf. Syst. Sci. 2000, 6, 318–332. 19. Javidi, M.; Nyamoradi, N. Dynamic analysis of a fractional order prey-predator interaction with harvesting. Appl. Math. Model. 2013, 37, 8946–8956. doi:10.1016/j.apm.2013.04.024. 20. Zhu, C.; Kong, L. Bifurcations analysis of Leslie-Gower predator-prey models with nonlinear predator-harvesting. Discret. Contin. Dyn. Syst. 2017, 10, 1187–1206. doi:10.3934/dcdss.2017065. 21. Suryanto, A.; Darti, I. Dynamics of Leslie-Gower pest-predator model with disease in pest including pest-harvesting and optimal implementation of pesticide. Int. J. Math. Math. Sci. 2019, 2019, 1–9. doi:10.1155/2019/5079171. Axioms 2020, 9, 122 21 of 23 22. Ang, T.K.; Safuan, H.M. Harvesting in a toxicated intraguild predator–prey ﬁshery model with variable carrying capacity. Chaos Solitons Fractals 2019, 126, 158–168. doi:10.1016/j.chaos.2019.06.004. 23. Leard, B.; Rebaza, J. Analysis of predator-prey models with continuous threshold harvesting. Appl. Math. Comput. 2011, 217, 5265–5278. doi:10.1016/j.amc.2010.11.050. 24. Bohn, J.; Rebaza, J.; Speer, K. Continuous threshold prey harvesting in predator-prey models. Int. J. Math. Comput. Sci. 2011, 5, 964–971. 25. Rebaza, J. Dynamics of prey threshold harvesting and refuge. J. Comput. Appl. Math. 2012, 236, 1743–1752. doi:10.1016/j.cam.2011.10.005. 26. Lv, Y.; Yuan, R.; Pei, Y. Dynamics in two nonsmooth predator–prey models with threshold harvesting. Nonlinear Dyn. 2013, 74, 107–132. doi:10.1007/s11071-013-0952-2. 27. Wu, D.; Zhao, H.; Yuan, Y. Complex dynamics of a diffusive predator–prey model with strong Allee effect and threshold harvesting. J. Math. Anal. Appl. 2019, 469, 982–1014. doi:10.1016/j.jmaa.2018.09.047. 28. Toaha, S. The effect of harvesting with threshold on the dynamics of prey predator model. J. Phys. Conf. Ser. 2019, 1341. doi:10.1088/1742-6596/1341/6/062021. 29. Panigoro, H.S.; Suryanto, A.; Kusumawinahyu, W.M.; Darti, I. Continuous threshold harvesting in a gause-type predator-prey model with fractional-order. AIP Conf. Proc. 2020, 2264, 040001. doi:10.1063/5.0023513. 30. Shepard, B.M.; Carner, G.R.; Barrion, A.T.; Ooi, P.A.C.; Van den Berg, H. Insects and Their Natural Enemies Associated with Vegetables and Soybean in Southeast Asia; Clemson Univ Coastal Research &: Orangeburg, USA, 31. Li, H.L.; Zhang, L.; Hu, C.; Jiang, Y.L.; Teng, Z. Dynamical analysis of a fractional-order predator-prey model incorporating a prey refuge. J. Appl. Math. Comput. 2017, 54, 435–449. doi:10.1007/s12190-016-1017-8. 32. Das, S.; Gupta, P.K. A mathematical model on fractional Lotka-Volterra equations. J. Theor. Biol. 2011, 277, 1–6. doi:10.1016/j.jtbi.2011.01.034. 33. Panja, P. Dynamics of a fractional order predator-prey model with intraguild predation. Int. J. Model. Simul. 2019, 39, 256–268. doi:10.1080/02286203.2019.1611311. 34. Suryanto, A.; Darti, I.; S. Panigoro, H.; Kilicman, A. A fractional-order predator–prey model with ratio-dependent functional response and linear harvesting. Mathematics 2019, 7, 1100. doi:10.3390/math7111100. 35. Alidousti, J.; Ghafari, E. Dynamic behavior of a fractional order prey-predator model with group defense. Chaos Solitons Fractals 2020, 134, 109688. doi:10.1016/j.chaos.2020.109688. 36. Ghanbari, B.; Djilali, S. Mathematical analysis of a fractional-order predator-prey model with prey social behavior and infection developed in predator population. Chaos Solitons Fractals 2020, 138, 109960. doi:10.1016/j.chaos.2020.109960. 37. Xie, Y.; Wang, Z.; Meng, B.; Huang, X. Dynamical analysis for a fractional-order prey–predator model with Holling III type functional response and discontinuous harvest. Appl. Math. Lett. 2020, 106, 106342. doi:10.1016/j.aml.2020.106342. 38. Sekerci, Y. Climate change effects on fractional order prey-predator model. Chaos Solitons Fractals 2020, 134, 109690. doi:10.1016/j.chaos.2020.109690. 39. Caputo, M. Linear models of dissipation whose Q is almost fFrequency independent–II. Geophys. J. Int. 1967, 13, 529–539. doi:10.1111/j.1365-246X.1967.tb02303.x. 40. Atangana, A.; Baleanu, D. New fractional derivatives with nonlocal and non-singular kernel: Theory and application to heat transfer model. Therm. Sci. 2016, 20, 763–769. doi:10.2298/TSCI160111018A. 41. Atangana, A.; Koca, I. Chaos in a simple nonlinear system with Atangana–Baleanu derivatives with fractional order. Chaos Solitons Fractals 2016, 89, 447–454. doi:10.1016/j.chaos.2016.02.012. 42. Tajadodi, H. A Numerical approach of fractional advection-diffusion equation with Atangana–Baleanu derivative. Chaos Solitons Fractals 2020, 130, 109527. doi:10.1016/j.chaos.2019.109527. 43. Ghanbari, B.; Kumar, D. Numerical solution of predator-prey model with Beddington-DeAngelis functional response and fractional derivatives with Mittag-Lefﬂer kernel. Chaos Interdiscip. J. Nonlinear Sci. 2019, 29, 063103. doi:10.1063/1.5094546. 44. Morales-Delgado, V.F.; Gómez-Aguilar, J.F.; Saad, K.; Escobar Jiménez, R.F. Application of the Caputo-Fabrizio and Atangana-Baleanu fractional derivatives to mathematical model of cancer chemotherapy effect. Math. Methods Appl. Sci. 2019, 42, 1167–1193. doi:10.1002/mma.5421. Axioms 2020, 9, 122 22 of 23 45. Shah, S.A.A.; Khan, M.A.; Farooq, M.; Ullah, S.; Alzahrani, E.O. A fractional order model for Hepatitis B virus with treatment via Atangana–Baleanu derivative. Phys. A Stat. Mech. Its Appl. 2020, 538, 122636. doi:10.1016/j.physa.2019.122636. 46. Bourafa, S.; Abdelouahab, M.S.; Moussaoui, A. On some extended Routh–Hurwitz conditions for fractional-order autonomous systems of order a 2 (0, 2) and their applications to some population dynamic models. Chaos Solitons Fractals 2020, 133. doi:10.1016/j.chaos.2020.109623. 47. Ahmed, E.; El-Sayed, A.; El-Saka, H.A. On some Routh–Hurwitz conditions for fractional order differential equations and their applications in Lorenz, Rössler, Chua and Chen systems. Phys. Lett. A 2006, 358, 1–4. doi:10.1016/j.physleta.2006.04.087. 48. Petras, I. Fractional-Order Nonlinear Systems: Modeling, Analysis and Simulation; Springer: London, UK; Beijing, China, 2011. 49. Diethelm, K. The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using Differential Operators of Caputo Type; Springer: Braunschweig, Germany, 2010. 50. Podlubny, I. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications; Academic Press: San Diego CA, USA, 1999. 51. Vargas-De-León, C. Volterra-type Lyapunov functions for fractional-order epidemic systems. Commun. Nonlinear Sci. Numer. Simul. 2015, 24, 75–85. doi:10.1016/j.cnsns.2014.12.013. 52. Huo, J.; Zhao, H.; Zhu, L. The effect of vaccines on backward bifurcation in a fractional order HIV model. Nonlinear Anal. Real World Appl. 2015, 26, 289–305. doi:10.1016/j.nonrwa.2015.05.014. 53. Abdelouahab, M.S.; Hamri, N.E.; Wang, J. Hopf bifurcation and caos in fractional-order modiﬁed hybrid optical system. Nonlinear Dyn. 2012, 69, 275–284. doi:10.1007/s11071-011-0263-4. 54. Deshpande, A.S.; Daftardar-Gejji, V.; Sukale, Y.V. On Hopf bifurcation in fractional dynamical systems. Chaos Solitons Fractals 2017, 98, 189–198. doi:10.1016/j.chaos.2017.03.034. 55. Diethelm, K. A fractional calculus based model for the simulation of an outbreak of dengue fever. Nonlinear Dyn. 2013, 71, 613–619. doi:10.1007/s11071-012-0475-2. 56. Moustafa, M.; Mohd, M.H.; Ismail, A.I.; Abdullah, F.A. Dynamical analysis of a fractional-order eco-epidemiological model with disease in prey population. Adv. Differ. Equ. 2020, 2020, 48. doi:10.1186/s13662-020-2522-5. 57. Li, Y.; Chen, Y.; Podlubny, I. Stability of fractional-order nonlinear dynamic systems: Lyapunov direct method and generalized Mittag–Lefﬂer stability. Comput. Math. Appl. 2010, 59, 1810–1821. doi:10.1016/j.camwa.2009.08.019. 58. Odibat, Z.M.; Shawagfeh, N.T. Generalized Taylor ’s formula. Appl. Math. Comput. 2007, 186, 286–293. doi:10.1016/j.amc.2006.07.102. 59. Matignon, D. Stability results for fractional differential equations with applications to control processing. Comput. Eng. Syst. Appl. 1996, pp. 963–968. 60. Kuznetsov, Y.A. Elements of Applied Bifurcation Theory, 3rd ed.; Springer: New York, NY, USA, 2004. 61. Baisad, K.; Moonchai, S. Analysis of stability and Hopf bifurcation in a fractional Gauss-type predator–prey model with Allee effect and Holling type-III functional response. Adv. Differ. Equ. 2018, 2018. doi:10.1186/s13662-018-1535-9. 62. Li, X.; Wu, R. Hopf bifurcation analysis of a new commensurate fractional-order hyperchaotic system. Nonlinear Dyn. 2014, 78, 279–288. doi:10.1007/s11071-014-1439-5. 63. Tavazoei, M.S.; Haeri, M. A proof for non existence of periodic solutions in time invariant fractional order systems. Automatica 2009, 45, 1886–1890. doi:10.1016/j.automatica.2009.04.001. 64. Fulger, D.; Scalas, E.; Germano, G. Monte Carlo simulation of uncoupled continuous-time random walks yielding a stochastic solution of the space-time fractional diffusion equation. Phys. Rev. E 2008, 77, 021122. doi:10.1103/PhysRevE.77.021122. 65. Scherer, R.; Kalla, S.L.; Tang, Y.; Huang, J. The Grünwald–Letnikov method for fractional differential equations. Comput. Math. Appl. 2011, 62, 902–917. doi:10.1016/j.camwa.2011.03.054. 66. Suryanto, A.; Darti, I. Stability analysis and nonstandard Grünwald-Letnikov scheme for a fractional order predator-prey model with ratio-dependent functional response. AIP Conf. Proc. 2017, 1913, 020011. doi:10.1063/1.5016645. 67. Diethelm, K.; Ford, N.J.; Freed, A.D. A predictor-corrector approach for the numerical solution of fractional differential equations. Nonlinear Dyn. 2002, 29, 3–22. doi:10.1023/A:1016592219341. Axioms 2020, 9, 122 23 of 23 68. Wang, Z. A numerical method for delayed fractional-order differential equations. J. Appl. Math. 2013, 2013, 256071. doi:10.1155/2013/256071. 69. Baleanu, D.; Jajarmi, A.; Hajipour, M. On the nonlinear dynamical systems within the generalized fractional derivatives with Mittag–Lefﬂer kernel. Nonlinear Dyn. 2018, 94, 397–414. doi:10.1007/s11071-018-4367-y. Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional afﬁliations. © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Journal

Axioms – Multidisciplinary Digital Publishing Institute

Published: Oct 22, 2020

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

A Rosenzweig–MacArthur Model with Continuous Threshold Harvesting in Predator Involving Fractional Derivatives with Power Law and Mittag–Leffler Kernel

A Rosenzweig–MacArthur Model with Continuous Threshold Harvesting in Predator Involving Fractional Derivatives with Power Law and Mittag–Leffler Kernel

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

A Rosenzweig–MacArthur Model with Continuous Threshold Harvesting in Predator Involving Fractional Derivatives with Power Law and Mittag–Leffler Kernel

A Rosenzweig–MacArthur Model with Continuous Threshold Harvesting in Predator Involving Fractional Derivatives with Power Law and Mittag–Leffler Kernel

References (71)

Abstract

Journal

Recommended Articles

There are no references for this article.

Our policy towards the use of cookies