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...
Panigoro, Hasan S.;Suryanto, Agus;Kusumawinahyu, Wuryansari Muharini;Darti, Isnani
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; firstname.lastname@example.org (H.S.P.); email@example.com (W.M.K.); firstname.lastname@example.org (I.D.) Department of Mathematics, Faculty of Mathematics and Natural Sciences, State University of Gorontalo, Bone Bolango 96119, Indonesia * Correspondence: email@example.com 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  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 , the Allee effect [10,11], the intraspeciﬁc competition [12,13], the cannibalism , 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.  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  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.  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 . 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.  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)  and Atangana–Baleanu operator . 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.  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.  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 ¯