Open Access proceedings Journal of Physics: Conference series How to cite: M. Ikbal and Riskawati, “Dynamics of Predator-Prey Model Interaction with Harvesting Effort”, J. Mat. Mantik, vol. 6, no. 2, pp.93-103, October 2020. Jurnal Matematika MANTIK Vol. 6, No. 2, October 2020, pp. 93-103 ISSN: 2527-3159 (print) 2527-3167 (online) Dynamics of Predator-Prey Model Interaction with Harvesting Effort Muhammad Ikbal1, Riskawati2 1Universitas Muslim Maros, ibbal@umma.ac.id 2Universitas Muslim Maros, riskawati@umma.ac.id doi: https://doi.org/10.15642/mantik.2020.6.2.93-103 Abstrak. Di dalam penelitian ini, kami mempelajari dan mengonstruksi suatu dinamika model mangsa-pemangsa. Kami memasukkan unsur kompetisi intraspesifik pada kedua pemangsa. Kami menformulasikan fungsi respon Holling tipe I pada masing-masing pemangsa. Kami menganggap semua populasi bernilai ekonomis sehingga dapat dipanen. Kami menganalisis solusi positifnya, keeksisan titik keseimbangannya, dan kestabilan pada titik-titik keseimbangannya itu. Kondisi kestabilan lokalnya kami peroleh dengan pendekatan kriteria Routh-Hurwitz. Kami juga mensimulasikan model tersebut. Penelitian ini bisa dikembangkan dengan formulasi fungsi respon yang berbeda dan pengoptimalan pemanenan. Kata kunci: Mangsa-pemangsa; Intraspesifik; Pemanenan; Routh-Hurwitz Abstract. In this research, we study and construct a dynamic prey-predator model. We include an element of intraspecific competition in both predators. We formulated the Holling type I response function for each predator. We consider all populations to be of economic value so that they can be harvested. We analyze the positive solution, the existence of the equilibrium points, and the stability of the balance points. We obtained the local stability condition by using the Routh-Hurwitz criterion approach. We also simulate the model. This research can be developed with different response function formulations and harvest optimization. Keywords: Prey-predator; Intraspecific; Harvesting; Routh-Hurwitz http://u.lipi.go.id/1458103791 mailto:ibbal@umma.ac.id mailto:riskawati@umma.ac.id Jurnal Matematika MANTIK Volume 6, No. 2, October 2020, pp. 93-103 94 1. Introduction Mathematical models can be used in observing individual behavior, population dynamics and population linkages in a system. Mathematical models can also be used in determining a policy. Mathematical modeling in the field of ecology is very interesting to study considering the many factors that affect the growth and life of living populations and the balance of organisms. The process of dynamics of organisms can be modeled mathematically by using differential equations involving continuous time or discrete time. One of the mathematical models used to explain this natural phenomenon is the prey- predator population model. Competition between predators and harvesting factors in populations is very important in the discipline of ecology. Many researchers can evoke interesting things from behavioral dynamics in population ecosystems. By combining the two aspects above, namely the aspects of competition between predators and harvesting, population dynamics can be expressed in a model. One of the policies related to the use of living things is harvesting. Intraspecific competition factors are also interesting to study. Intraspecific competition is competition between predators in competing for prey. This is another factor in population dynamics that can affect the stability of a system. There are many researchers who model prey-predator interactions. [1] examined the resistance of predators in the prey-predator model system with non-periodic solutions. [2] discusses the dynamics of the prey-predator with diseased predators. [3] in his journal discussed global dynamics of a prey-predator model with antipredator behavior and two predators. [4] discuss the dynamics of the prey-predator model by quadratic harvesting. Research from [5] discusses global dynamics and control of predator prey models with Holling type III response functions. [6] examined the effect of harvesting and competition between predators in the prey-predator model. [7] discuss a model of interaction of three species in one habitat. [8] discuss the complex dynamics of a three-species food chain model with Holling type III response functions. Many previous researchers have examined prey-predator population models. We examine prey-predator population models with respect to intraspecific competition for predators and considering the economic value of all populations. Our study constructs the factors influencing prey-predator population dynamics as investigated by previous researchers, but we add intraspecific competition and harvesting factors simultaneously to all three populations. 2. Assumptions and model In this model, there is an interference between predators as modeled by other researchers [1], [6], [8], [9]. There are researchers who studied the intraspecific competition coefficient [9]. Researchers frequently use the Holling-type I response function [3], [5], [8]. The response function is used by researchers in their models [10]– [13]. Researchers also use the harvesting rate [4], [6], [14], [15]. The assumptions used are: • The prey growth rate uses the logistical growth rate, • Predators compete with each other for prey, • All of predators uses the Holling type I response function for predation, • There is an intraspecific competition for each predator. • All of population have interest economic values. The model is formulated as follows: 𝑑𝑃 𝑑𝑡 = 𝑟𝑃 (1 − 𝑃 𝐾 ) − 𝛼1𝑃𝐻1 − 𝛼2𝑃𝐻2 − 𝑞1𝐸1𝑃 (1) M. Ikbal, Riskawati Dynamics of Predator-Prey Model Interaction with Harvesting Effort 95 𝑑𝐻1 𝑑𝑡 = 𝑒1𝛼1𝑃𝐻1 − 𝑔1𝐻1 2 − 𝛽1𝐻1𝐻2 − 𝑑1𝐻1 − 𝑞2𝐸2𝐻1 𝑑𝐻2 𝑑𝑡 = 𝑒2𝛼2𝑃𝐻2 − 𝑔2𝐻2 2 − 𝛽2𝐻1𝐻2 − 𝑑2𝐻2 − 𝑞3𝐸3𝐻2 with initial condition P(0) ≥ 0, H1(0) ≥ 0, H2(0) ≥ 0. (2) The first predator (H1) and second predator (H2) are assumed to have direct access to prey (P). The effect of disturbance in the growth rate of competitors is assumed to be proportional to the density of the predator population with β 1 and β 2 respectively given disturbance rates. The parameter α1 and α2 represent predator rates for predator species H1 and H2 respectively. The parameter g1 and g2 represent coefficient intraspecific competition for two predators H1 and H2 respectively, here d1 and d2 are their mortality rate. The parameter e1 and e2 are the predator’s conversion efficiency. However, the predation functions of the two predators were made different - one following a Holling type I response and the other following a Holling type II response. Besides experiencing a reduction due to the predation function, the prey population grew logistically with r as the intrinsic growth rate and K as the holding capacity. The parameter q 1 , q 2 , and q 3 are the catchability coefficient of susceptible prey, first predator, and second predator, respectively. The parameter E1, E2, and E3 are the harvesting effort prey, first predator, and second predator, respectively. 3. Equilibrium points and stability analysis 3.1. Equilibrium points Equilibrium points of the system (1) are given below: • The trivial equilibrium point 𝑇0= (0, 0, 0), • The predator free equilibrium point 𝑇1 = (𝐾 − 𝐾𝑞 1 𝐸1 𝑟 , 0, 0), • The H1-free boundary equilibrium state 𝑇2 = ( 𝐾(𝑞 3 𝐸3𝛼2 + 𝛼2𝑑2 + 𝑟𝑔2 − 𝑞1𝑔2𝐸1) 𝐾𝑒2𝛼2 2 + 𝑟𝑔 2 , 0, 𝐾𝑟𝑒2𝛼2 − 𝑟𝑞3𝐸3 − 𝑑2𝑟 − 𝐾𝑒2𝛼2𝑞1𝐸1 𝐾𝑒2𝛼2 2 + 𝑟𝑔 2 ) • The H2-free boundary equilibrium state 𝑇3 = ( 𝐾(𝑞2𝐸2𝛼1 + 𝛼1𝑑1 + 𝑟𝑔1 − 𝑞1𝑔1𝐸1) 𝐾𝑒1𝛼1 2 + 𝑟𝑔1 , 𝐾𝑟𝑒1𝛼1 − 𝑟𝑞2𝐸2 − 𝑑1𝑟 − 𝐾𝑒1𝛼1𝑞1𝐸1 𝐾𝑒1𝛼1 2 + 𝑟𝑔1 , 0) • The interior equilibrium point 𝑇4 = (𝑃 ∗, 𝐻1 ∗ , 𝐻2 ∗), where P * = 𝐾(𝑞1𝐸1𝑔1𝑔2 + 𝑞2𝐸2𝛼2𝛽2 + 𝑞3𝐸3𝛼1𝛽1 + 𝑑1𝛼2𝛽2 + 𝑑2𝛼1𝛽1 + 𝑟𝛽1𝛽2) 𝐾𝑒1𝛼1𝛼2𝛽2 + 𝐾𝑒2𝛼1𝛼2𝛽1 + 𝑟𝛽1𝛽2 − 𝐾𝑒1𝛼1 2𝑔2 − 𝐾𝑒2𝛼2 2𝑔1 − 𝑟𝑔1𝑔2 − 𝐾(𝑞1𝐸1𝛽1𝛽2 + 𝑞2𝐸2𝛼1𝑔2 + 𝑞3𝐸3𝛼2𝑔1 + 𝑑1𝛼1𝑔2 + 𝑑2𝛼2𝑔1 + 𝑟𝑔1𝑔2) 𝐾𝑒1𝛼1𝛼2𝛽2 + 𝐾𝑒2𝛼1𝛼2𝛽1 + 𝑟𝛽1𝛽2 − 𝐾𝑒1𝛼1 2𝑔2 − 𝐾𝑒2𝛼2 2𝑔1 − 𝑟𝑔1𝑔2 H1 * = 𝐾𝑒1𝛼1𝑞1𝐸1𝑔2 + 𝐾𝑒2𝛼2 2𝑞2𝐸2 + 𝐾𝑑1𝑒2𝛼2 2 + 𝐾𝑟𝑒2𝛼2𝛽1 + 𝑔2𝑟𝑞2𝐸2 + 𝑟𝑑1𝑔2 𝐾𝑒1𝛼1𝛼2𝛽2 + 𝐾𝑒2𝛼1𝛼2𝛽1 + 𝑟𝛽1𝛽2 − 𝐾𝑒1𝛼1 2𝑔2 − 𝐾𝑒2𝛼2 2𝑔1 − 𝑟𝑔1𝑔2 − (𝐾𝑒1𝛼1𝛼2𝑞3𝐸3 + 𝐾𝑒2𝛼2𝛽1𝑞1𝐸1 + 𝐾𝑒1𝛼1𝛼2𝑑2 + 𝐾𝑟𝑒1𝛼1𝑔2 + 𝑟𝛽1 𝑞3𝐸3 + 𝑟𝑑2𝛽1) 𝐾𝑒1𝛼1𝛼2𝛽2 + 𝐾𝑒2𝛼1𝛼2𝛽1 + 𝑟𝛽1𝛽2 − 𝐾𝑒1𝛼1 2𝑔2 − 𝐾𝑒2𝛼2 2𝑔1 − 𝑟𝑔1𝑔2 H2 * = 𝐾𝑒1𝛼1𝑞1𝐸1𝛽2 + 𝐾𝑒2𝛼1𝛼2𝑞2𝐸2 + 𝐾𝑑1𝑒2𝛼1𝛼2 + 𝐾𝑟𝑒2𝛼2𝑔1 + 𝛽2𝑟𝑞2𝐸2 + 𝑟𝑑1𝛽2 𝐾𝑒1𝛼1𝛼2𝛽2 + 𝐾𝑒2𝛼1𝛼2𝛽1 + 𝑟𝛽1𝛽2 − 𝐾𝑒1𝛼1 2𝑔2 − 𝐾𝑒2𝛼2 2𝑔1 − 𝑟𝑔1𝑔2 − (𝐾𝑒1𝛼1 2𝑞3𝐸3 + 𝐾𝑒2𝛼2𝑔1𝑞1𝐸1 + 𝐾𝑒1𝛼1 2𝑑2 + 𝐾𝑟𝑒1𝛼1𝛽2 + 𝑟𝑔1𝑞3𝐸3 + 𝑟𝑑2𝑔1) 𝐾𝑒1𝛼1𝛼2𝛽2 + 𝐾𝑒2𝛼1𝛼2𝛽1 + 𝑟𝛽1𝛽2 − 𝐾𝑒1𝛼1 2𝑔2 − 𝐾𝑒2𝛼2 2𝑔1 − 𝑟𝑔1𝑔2 Jurnal Matematika MANTIK Volume 6, No. 2, October 2020, pp. 93-103 96 3.2. Stability analysis The stability analysis equilibrium point of the system (1) is studied and determined. The point 𝑇0 is trivial equilibrium point. Jacobian matrix of the model system (1) is 𝐽 = [ 𝐽11 𝐽12 𝐽13 𝐽21 𝐽22 𝐽23 𝐽31 𝐽32 𝐽33 ] (3) where, 𝐽11 = 𝑟 − 2𝑟𝑃 𝐾 − 𝛼1𝐻1 − 𝛼2𝐻2 − 𝑞1𝐸1 𝐽12 = −𝛼1𝑃 𝐽13 = −𝛼2𝑃 𝐽21 = 𝑒1𝛼1𝐻1 𝐽22 = 𝑒1𝛼1𝑃 − 2𝑔1𝐻1 − 𝛽1𝐻2 − 𝑑1 − 𝑞2𝐸2 𝐽23 = −𝛽1𝐻1 𝐽31 = 𝑒2𝛼2𝐻2 𝐽32 = −𝛽2𝐻2 𝐽33 = 𝑒2𝛼2𝑃 − 2𝑔2𝐻2 − 𝛽2𝐻1 − 𝑑2 − 𝑞3𝐸3 Theorem 1. Equilibrium point 𝑇1 local stable if 𝑟 > 𝑞1𝐸1, 𝑞2𝐸2 + 𝑑1 > 𝐾𝑒1𝛼1(𝑟−𝑞1𝐸1) 𝑟 , and 𝑞3𝐸3 + 𝑑2 > 𝐾𝑒2𝛼2(𝑟−𝑞1𝐸1) 𝑟 . Proof. The result of substitution equilibrium point T1 to Jacobian Matrix (3) 𝐽(𝑇1) = [ 𝐽11 1 𝐽12 1 𝐽13 1 𝐽21 1 𝐽22 1 𝐽23 1 𝐽31 1 𝐽32 1 𝐽33 1 ] (4) where 𝐽11 1 = 𝑞1𝐸1 − 𝑟 𝐽12 1 = − 𝛼1𝐾(𝑟 − 𝑞1𝐸1) 𝑟 𝐽13 1 = − 𝛼2𝐾(𝑟 − 𝑞1𝐸1) 𝑟 𝐽21 1 = 0 𝐽22 1 = 𝐾𝑒1𝛼1(𝑟 − 𝑞1𝐸1) 𝑟 − 𝑞2𝐸2 − 𝑑1 𝐽23 1 = 0 𝐽31 1 = 0 𝐽32 1 = 0 𝐽33 2 = 𝐾𝑒2𝛼2(𝑟 − 𝑞1𝐸1) 𝑟 − 𝑞3𝐸3 − 𝑑2 Characteristic equation matrix 𝐽(𝑇1) is (𝜆 − 𝑞1𝐸1 + 𝑟) (𝜆 − 𝐾𝑒1𝛼1(𝑟 − 𝑞1𝐸1) 𝑟 + 𝑞2𝐸2 + 𝑑1) (𝜆 − 𝐾𝑒2𝛼2(𝑟 − 𝑞1𝐸1) 𝑟 + 𝑞3𝐸3 + 𝑑2) = 0 (5) The roots of the equation (5) is negative if 𝑟 > 𝑞1𝐸1, 𝑞2𝐸2 + 𝑑1 > 𝐾𝑒1𝛼1(𝑟−𝑞1𝐸1) 𝑟 , and 𝑞3𝐸3 + 𝑑2 > 𝐾𝑒2𝛼2(𝑟−𝑞1𝐸1) 𝑟 . M. Ikbal, Riskawati Dynamics of Predator-Prey Model Interaction with Harvesting Effort 97 Theorem 2. Equilibrium point 𝑇2 local stable if 𝐽11 2 < 0, 𝐽22 2 < 0, 𝐽33 2 < 0, 𝐽11 2 𝐽22 2 + 𝐽11 2 𝐽33 2 + 𝐽22 2 𝐽33 2 > 𝐽13 2 𝐽31 2 and 𝐽13 2 𝐽31 2 𝐽22 2 > 𝐽11 2 𝐽22 2 𝐽33 2 . Proof. The result of substitution equilibrium point T2 to Jacobian Matrix (3) 𝐽(𝑇2) = [ 𝐽11 2 𝐽12 2 𝐽13 2 𝐽21 2 𝐽22 2 𝐽23 2 𝐽31 2 𝐽32 2 𝐽33 2 ] (6) where 𝐽11 2 =𝑟 − 2𝑟(𝑞3𝐸3𝛼2 + 𝛼2𝑑2 + 𝑟𝑔2 − 𝑞1𝑔2𝐸1) 𝐾𝑒2𝛼2 2 + 𝑟𝑔2 −𝛼2 [ 𝐾𝑟𝑒2𝛼2 − 𝑟𝑞3𝐸3 − 𝑑2𝑟 − 𝐾𝑒2𝛼2𝑞1𝐸1 𝐾𝑒2𝛼2 2 + 𝑟𝑔2 ] − 𝑞1𝐸1 𝐽12 2 = −𝛼1 [ 𝐾(𝑞3𝐸3𝛼2 + 𝛼2𝑑2 + 𝑟𝑔2 − 𝑞1𝑔2𝐸1) 𝐾𝑒2𝛼2 2 + 𝑟𝑔2 ] 𝐽13 2 = −𝛼2 [ 𝐾(𝑞3𝐸3𝛼2 + 𝛼2𝑑2 + 𝑟𝑔2 − 𝑞1𝑔2𝐸1) 𝐾𝑒2𝛼2 2 + 𝑟𝑔2 ] 𝐽21 2 = 0 𝐽22 2 = 𝑒1𝛼1 [ 𝐾(𝑞3𝐸3𝛼2 + 𝛼2𝑑2 + 𝑟𝑔2 − 𝑞1𝑔2𝐸1) 𝐾𝑒2𝛼2 2 + 𝑟𝑔2 ] −𝛽1 [ 𝐾𝑟𝑒2𝛼2 − 𝑟𝑞3𝐸3 − 𝑑2𝑟 − 𝐾𝑒2𝛼2𝑞1𝐸1 𝐾𝑒2𝛼2 2 + 𝑟𝑔2 ] − 𝑑1 − 𝑞2𝐸2 𝐽23 2 = 0 𝐽31 2 = 𝑒2𝛼2 [ 𝐾𝑟𝑒2𝛼2 − 𝑟𝑞3𝐸3 − 𝑑2𝑟 − 𝐾𝑒2𝛼2𝑞1𝐸1 𝐾𝑒2𝛼2 2 + 𝑟𝑔2 ] 𝐽32 2 = −𝛽2 [ 𝐾𝑟𝑒2𝛼2 − 𝑟𝑞3𝐸3 − 𝑑2𝑟 − 𝐾𝑒2𝛼2𝑞1𝐸1 𝐾𝑒2𝛼2 2 + 𝑟𝑔2 ] 𝐽33 2 = 𝑒2𝛼2 [ 𝐾(𝑞3𝐸3𝛼2 + 𝛼2𝑑2 + 𝑟𝑔2 − 𝑞1𝑔2𝐸1) 𝐾𝑒2𝛼2 2 + 𝑟𝑔2 ] −2𝑔2 [ 𝐾𝑟𝑒2𝛼2 − 𝑟𝑞3𝐸3 − 𝑑2𝑟 − 𝐾𝑒2𝛼2𝑞1𝐸1 𝐾𝑒2𝛼2 2 + 𝑟𝑔2 ] − 𝑑2 − 𝑞3𝐸3 Characteristic equation matrix 𝐽(𝑇2) is 𝜆3 + 𝐴1𝜆 2 + 𝐴2𝜆 + 𝐴3 = 0 (7) where, 𝐴1 = −(𝐽11 2 + 𝐽22 2 + 𝐽33 2 ) 𝐴2 = 𝐽11 2 𝐽22 2 + 𝐽11 2 𝐽33 2 + 𝐽22 2 𝐽33 2 − 𝐽13 2 𝐽31 2 𝐴3 = 𝐽13 2 𝐽31 2 𝐽22 2 − 𝐽11 2 𝐽22 2 𝐽33 2 The roots of the equation (7) is negative if 𝐽11 2 < 0, 𝐽22 2 < 0, 𝐽33 2 < 0, 𝐽11 2 𝐽22 2 + 𝐽11 2 𝐽33 2 + 𝐽22 2 𝐽33 2 > 𝐽13 2 𝐽31 2 , 𝐽13 2 𝐽31 2 𝐽22 2 > 𝐽11 2 𝐽22 2 𝐽33 2 and A1A2 > A3 Theorem 3. Equilibrium point 𝑇3 local stable if 𝐽11 3 + 𝐽22 3 + 𝐽33 3 < 0, 𝐽11 3 𝐽22 3 + 𝐽11 3 𝐽33 3 + 𝐽22 3 𝐽33 3 > 𝐽12 3 𝐽21 3 , 𝐽12 3 𝐽21 3 𝐽33 3 > 𝐽11 3 𝐽22 3 𝐽33 3 , and 𝐵1𝐵2 > 𝐵3 Proof. The result of substitution equilibrium point 𝑇3 to Jacobian Matrix (3) Jurnal Matematika MANTIK Volume 6, No. 2, October 2020, pp. 93-103 98 𝐽(𝑇3) = [ 𝐽11 3 𝐽12 3 𝐽13 3 𝐽21 3 𝐽22 3 𝐽23 3 𝐽31 3 𝐽32 3 𝐽33 3 ] (8) where 𝐽11 3 = 𝑟 − 2𝑟 [ (𝑞2𝐸2𝛼1+𝛼1𝑑1+𝑟𝑔1−𝑞1 𝑔1𝐸1) 𝐾𝑒1𝛼1 2+𝑟𝑔1 ] − 𝛼1 [ 𝐾𝑟𝑒1𝛼1−𝑟𝑞2𝐸2−𝑑1𝑟−𝐾𝑒1𝛼1𝑞1 𝐸1 𝐾𝑒1𝛼1 2+𝑟𝑔1 ] − 𝑞1𝐸1 𝐽12 3 = − 𝛼1 [ (𝑞2𝐸2𝛼1 + 𝛼1𝑑1 + 𝑟𝑔1 − 𝑞1𝑔1𝐸1) 𝐾𝑒1𝛼1 2 + 𝑟𝑔1 ] 𝐽13 3 = −𝛼2 [ (𝑞2𝐸2𝛼1 + 𝛼1𝑑1 + 𝑟𝑔1 − 𝑞1𝑔1𝐸1) 𝐾𝑒1𝛼1 2 + 𝑟𝑔1 ] 𝐽21 3 = 𝑒1𝛼1 [ 𝐾𝑟𝑒1𝛼1 − 𝑟𝑞2𝐸2 − 𝑑1𝑟 − 𝐾𝑒1𝛼1𝑞1𝐸1 𝐾𝑒1𝛼1 2 + 𝑟𝑔1 ] 𝐽22 3 = 𝑒1𝛼1𝑃 − 2𝑔1 [ 𝐾𝑟𝑒1𝛼1 − 𝑟𝑞2𝐸2 − 𝑑1𝑟 − 𝐾𝑒1𝛼1𝑞1𝐸1 𝐾𝑒1𝛼1 2 + 𝑟𝑔1 ] −𝛽1 [ 𝐾𝑟𝑒1𝛼1 − 𝑟𝑞2𝐸2 − 𝑑1𝑟 − 𝐾𝑒1𝛼1𝑞1𝐸1 𝐾𝑒1𝛼1 2 + 𝑟𝑔1 ] − 𝑑1 − 𝑞2𝐸2 𝐽23 3 = −𝛽1 [ 𝐾𝑟𝑒1𝛼1 − 𝑟𝑞2𝐸2 − 𝑑1𝑟 − 𝐾𝑒1𝛼1𝑞1𝐸1 𝐾𝑒1𝛼1 2 + 𝑟𝑔1 ] 𝐽31 3 = 0 𝐽32 3 = 0 𝐽33 3 = 𝑒2𝛼2 [ 𝐾(𝑞2𝐸2𝛼1 + 𝛼1𝑑1 + 𝑟𝑔1 − 𝑞1𝑔1𝐸1) 𝐾𝑒1𝛼1 2 + 𝑟𝑔1 ] −𝛽2 [ 𝐾𝑟𝑒1𝛼1 − 𝑟𝑞2𝐸2 − 𝑑1𝑟 − 𝐾𝑒1𝛼1𝑞1𝐸1 𝐾𝑒1𝛼1 2 + 𝑟𝑔1 ] − 𝑑2 − 𝑞3𝐸3 Characteristics equation matrix 𝐽(𝑇3) is 𝜆3 + 𝐵1𝜆 2 + 𝐵2𝜆 + 𝐵3 = 0 (9) where, 𝐵1 = −(𝐽11 3 + 𝐽22 3 + 𝐽33 3 ) 𝐵2 = 𝐽11 3 𝐽22 3 + 𝐽11 3 𝐽33 3 + 𝐽22 3 𝐽33 3 − 𝐽12 3 𝐽21 3 𝐵3 = 𝐽12 3 𝐽21 3 𝐽33 3 − 𝐽11 3 𝐽22 3 𝐽33 3 . To ensure the stability of model system with equilibrium point 𝑇3, the point must qualify of the Routh-Hurtwiz criteria. The equation (9) have negative roots if 𝐽11 3 + 𝐽22 3 + 𝐽33 3 < 0, 𝐽11 3 𝐽22 3 + 𝐽11 3 𝐽33 3 + 𝐽22 3 𝐽33 3 > 𝐽12 3 𝐽21 3 , 𝐽12 3 𝐽21 3 𝐽33 3 > 𝐽11 3 𝐽22 3 𝐽33 3 , and 𝐵1𝐵2 > 𝐵3. Theorem 4 Equilibrium point 𝑇4 local stable if 𝐽11 4 + 𝐽22 4 + 𝐽33 4 < 0, 𝐽11 4 𝐽22 4 + 𝐽11 4 𝐽33 4 + 𝐽22 4 𝐽33 4 > 𝐽12 4 𝐽21 4 + 𝐽13 4 𝐽31 4 + 𝐽23 4 𝐽32 4 , 𝐽11 4 𝐽23 4 𝐽32 4 + 𝐽12 4 𝐽21 4 𝐽33 4 + 𝐽13 4 𝐽31 4 𝐽22 4 > 𝐽11 4 𝐽22 4 𝐽33 4 + 𝐽12 4 𝐽23 4 𝐽31 4 + 𝐽13 4 𝐽32 4 𝐽21 4 , and 𝐶1𝐶2 > 𝐶3. Proof. The result of substitution equilibrium point T4 to Jacobian Matrix (3) 𝐽(𝑇4) = [ 𝐽11 4 𝐽12 4 𝐽13 4 𝐽21 4 𝐽22 4 𝐽23 4 𝐽31 4 𝐽32 4 𝐽33 4 ] (10) where M. Ikbal, Riskawati Dynamics of Predator-Prey Model Interaction with Harvesting Effort 99 𝐽11 4 = 𝑟 − 2𝑟𝑃∗ 𝐾 − 𝛼1𝐻1 ∗ − 𝛼2𝐻2 ∗ − 𝑞1𝐸1 𝐽12 4 = −𝛼1𝑃 ∗ 𝐽13 4 = −𝛼2𝑃 ∗ 𝐽21 4 = 𝑒1𝛼1𝐻1 ∗ 𝐽22 4 = 𝑒1𝛼1𝑃 ∗ − 2𝑔1𝐻1 ∗ − 𝛽1𝐻2 ∗ − 𝑑1 − 𝑞2𝐸2 𝐽23 4 = −𝛽1𝐻1 ∗ 𝐽31 4 = 𝑒2𝛼2𝐻2 ∗ 𝐽32 4 = −𝛽2𝐻2 ∗ 𝐽33 4 = 𝑒2𝛼2𝑃 ∗ − 2𝑔2𝐻2 − 𝛽2𝐻1 ∗ − 𝑑2 − 𝑞3𝐸3 Characteristics equation matrix J(T4) is 𝜆3 + 𝐶1𝜆 2 + 𝐶2𝜆 + 𝐶3 = 0 (11) where, 𝐶1 = −(𝐽11 4 + 𝐽22 4 + 𝐽33 4 ) 𝐶2 = 𝐽11 4 𝐽22 4 + 𝐽11 4 𝐽33 4 + 𝐽22 4 𝐽33 4 − 𝐽12 4 𝐽21 4 − 𝐽13 4 𝐽31 4 − 𝐽23 4 𝐽32 4 𝐶3 = 𝐽11 4 𝐽23 4 𝐽32 4 + 𝐽12 4 𝐽21 4 𝐽33 4 + 𝐽13 4 𝐽31 4 𝐽22 4 − 𝐽11 4 𝐽22 4 𝐽33 4 − 𝐽12 4 𝐽23 4 𝐽31 4 − 𝐽13 4 𝐽32 4 𝐽21 4 . To ensure the stability of model system with equilibrium point E4, the point must qualify of the Routh-Hurtwiz criteria. The equation (11) have negative roots if 𝐽11 4 + 𝐽22 4 + 𝐽33 4 < 0, 𝐽11 4 𝐽22 4 + 𝐽11 4 𝐽33 4 + 𝐽22 4 𝐽33 4 > 𝐽12 4 𝐽21 4 + 𝐽13 4 𝐽31 4 + 𝐽23 4 𝐽32 4 , 𝐽11 4 𝐽23 4 𝐽32 4 + 𝐽12 4 𝐽21 4 𝐽33 4 + 𝐽13 4 𝐽31 4 𝐽22 4 > 𝐽11 4 𝐽22 4 𝐽33 4 + 𝐽12 4 𝐽23 4 𝐽31 4 + 𝐽13 4 𝐽32 4 𝐽21 4 , and 𝐶1𝐶2 > 𝐶3. 4. Numerical Simulation In this section we simulated the model with some parameter values. The parameter values was adopted from literature [3], [6], [8], [9], [11], [15]. We try simulated the model with some condition. The first condition with parameter with E1 = 0.28773096, E2 = 0.24093140 and E3 = 0.13141604. The second condition with parameter E1 = 0.2, E2 = 0.4 and E3 = 0.3. The third condition with E1 = 0.2, E2 = 0.3 and E3 = 0.4. To see the system is in a stable state, a numerical simulation is performed with parameter estimates according to the following Table 1. Table 1. Parameter values Parameter Values r 1 K 100 α1 0.21 α2 0.274 e1 1 e2 1 g 1 0.1 g 2 0.2 β 1 0.05 β 2 0.06 d1 0.05 d2 0.06 𝑞1 1 𝑞2 1 𝑞3 1 𝐸1 0.5 Jurnal Matematika MANTIK Volume 6, No. 2, October 2020, pp. 93-103 100 (a) (c) (b) (b) Parameter Values 𝐸2 0.5 𝐸3 0.5 With the parameter values in Table 1, the simulation results are given nonnegative equilibrium points: T0 = (0, 0, 0) T1 = (50, 0, 0), T2 = (3.288183092, 0, 1.704810836), T3 = (3.669623060, 2.206208426 , 0), T4 = (3.149229952, 0.4190122146, 1.388741370). Figure 1. Numerical simulation model. (a) prey population density; (b) first predator population density; (c) second predator population density M. Ikbal, Riskawati Dynamics of Predator-Prey Model Interaction with Harvesting Effort 101 Figure 2. Time series of the model with E1 = 0.5, E2 = 0.5 and E3 = 0.5 Figure 3. Time series of the model with E1 = 0.2, E2 = 0.4 and E3 = 0.3 Figure 4. Time series of the model with E1 = 0.2, E2 = 0.3 and E3 = 0.4 Jurnal Matematika MANTIK Volume 6, No. 2, October 2020, pp. 93-103 102 Figure 1 shows the existence of each population. With a combination of the parameters presented, the food supply of predators is more than the predators themselves. the system is expected to last a long time with a combination of these parameters. Figure 2 shows the simulation results with given parameters with the same parameter (E1 = 0.5, E2 = 0.5 and E3 = 0.5) of harvest effort for all populations. In this simulation, the system is stable and lasts a long time with the same harvesting effort conditions. To demonstrate system dynamics, we simulate models with varying combinations of harvest effort parameters. By changing the number of harvesting business parameters to E1 = 0.2, E2 = 0.4 and E3 = 0.3, the simulation plot will look different. Figure 3 shows that the number of the first predator population is reduced compared to other populations because the number of harvesting efforts is also the highest. With different parameter numbers (E1 = 0.2, E2 = 0.3 and E3 = 0.4), Figure 4 shows the same symptoms, the second predator population is the smallest because the number of harvesting efforts is the largest. Apart from the stability of the system, the policy in harvesting in this study is a matter of focus. In this study, harvesting also determines the stability of a system, if the harvesting of a population changes and does not match the harvest rates in other populations, it will cause system stability disturbances. It can be seen clearly in the simulation that produces dynamic graphs in Figures 2, 3, and 4, the population density of one predator is sometimes more than other predators, and vice versa. The results of this study are intended to show in general that harvesting efforts have a large enough impact on population sustainability and also have an impact on the system. Harvesting effort can interfere with the growth and activity of predators. We look visually in the image with the selection of different harvesting efforts. Population that is harvested in large numbers to other populations will decrease in population. This research can be developed by considering other factors that make a system dynamic. For example, further researchers can optimize harvesting efforts so that economic benefits can be clearly measured. 5. Conclusions In this section we will make conclusions of this research. This study focuses on the dynamics of prey-predator populations with harvesting effort for all of population. There are 5 non-negative equilibrium point of the system. The interior point is 𝑻𝟒. The equilibrium point 𝑻𝟒 stabel if 𝑱𝟏𝟏 𝟒 + 𝑱𝟐𝟐 𝟒 + 𝑱𝟑𝟑 𝟒 < 𝟎, 𝑱𝟏𝟏 𝟒 𝑱𝟐𝟐 𝟒 + 𝑱𝟏𝟏 𝟒 𝑱𝟑𝟑 𝟒 + 𝑱𝟐𝟐 𝟒 𝑱𝟑𝟑 𝟒 > 𝑱𝟏𝟐 𝟒 𝑱𝟐𝟏 𝟒 + 𝑱𝟏𝟑 𝟒 𝑱𝟑𝟏 𝟒 + 𝑱𝟐𝟑 𝟒 𝑱𝟑𝟐 𝟒 , 𝑱𝟏𝟏 𝟒 𝑱𝟐𝟑 𝟒 𝑱𝟑𝟐 𝟒 + 𝑱𝟏𝟐 𝟒 𝑱𝟐𝟏 𝟒 𝑱𝟑𝟑 𝟒 + 𝑱𝟏𝟑 𝟒 𝑱𝟑𝟏 𝟒 𝑱𝟐𝟐 𝟒 > 𝑱𝟏𝟏 𝟒 𝑱𝟐𝟐 𝟒 𝑱𝟑𝟑 𝟒 + 𝑱𝟏𝟐 𝟒 𝑱𝟐𝟑 𝟒 𝑱𝟑𝟏 𝟒 + 𝑱𝟏𝟑 𝟒 𝑱𝟑𝟐 𝟒 𝑱𝟐𝟏 𝟒 , and 𝑪𝟏𝑪𝟐 > 𝑪𝟑. Harvesting effort have impact to system. Harvesting efforts will reduce the population size so that it can affect the stability of the dynamics of the prey-predator population system. The system will be stable and exist if we have control quantity of the harvesting effort. The system will be stable and exist if we have control’s the harvesting effort. 6. Acknowledgment Thanks to friends who helped this research and the Directorate of Research and Community Service Deputy for Strengthening Research and Development, The Ministry of Research and Technology/ National Research and Innovation Agency in Indonesia which has funded this research. References [1] J. Alebraheem and Y. Abu-Hasan, “Persistence of predators in a two predators-one M. Ikbal, Riskawati Dynamics of Predator-Prey Model Interaction with Harvesting Effort 103 prey model with non-periodic solution,” Appl. Math. Sci., vol. 6, no. 17–20, 2012. [2] K. pada Das, “A Mathematical Study of a Predator-Prey Dynamics with Disease in Predator,” ISRN Appl. Math., vol. 2011, 2011, doi: 10.5402/2011/807486. [3] C. Li, Y. Zhang, J. Xu, and Y. Zhou, “Global Dynamics of a Prey-Predator Model with Antipredator Behavior and Two Predators,” Discret. Dyn. Nat. Soc., vol. 2019, 2019, doi: 10.1155/2019/3586508. [4] R. P. Gupta and P. Chandra, “Dynamical properties of a prey-predator-scavenger model with quadratic harvesting,” Commun. Nonlinear Sci. Numer. Simul., vol. 49, 2017, doi: 10.1016/j.cnsns.2017.01.026. [5] T. K. Kar and H. Matsuda, “Global dynamics and controllability of a harvested prey- predator system with Holling type III functional response,” Nonlinear Anal. Hybrid Syst., vol. 1, no. 1, 2007, doi: 10.1016/j.nahs.2006.03.002. [6] B. Mukhopadhyay and R. Bhattacharyya, “Effects of harvesting and predator interference in a model of two-predators competing for a single prey,” Appl. Math. Model., vol. 40, no. 4, 2016, doi: 10.1016/j.apm.2015.10.018. [7] J. N. Ndam, J. P. Chollom, and T. G. Kassem, “A Mathematical Model of Three- Species Interactions in an Aquatic Habitat,” ISRN Appl. Math., vol. 2012, 2012, doi: 10.5402/2012/391547. [8] R. K. Upadhyay and S. N. Raw, “Complex dynamics of a three species food-chain model with Holling type IV functional response,” Nonlinear Anal. Model. Control, vol. 16, no. 3, 2011, doi: 10.15388/na.16.3.14098. [9] N. Ali, M. Haque, E. Venturino, and S. Chakravarty, “Dynamics of a three species ratio-dependent food chain model with intra-specific competition within the top predator,” Comput. Biol. Med., vol. 85, 2017, doi: 10.1016/j.compbiomed.2017.04.007. [10] A. Rojas-Palma and E. González-Olivares, “Optimal harvesting in a predator-prey model with Allee effect and sigmoid functional response,” Appl. Math. Model., vol. 36, no. 5, 2012, doi: 10.1016/j.apm.2011.07.081. [11] A. Chatterjee and S. Pal, “Interspecies competition between prey and two different predators with Holling IV functional response in diffusive system,” Comput. Math. with Appl., vol. 71, no. 2, 2016, doi: 10.1016/j.camwa.2015.12.022. [12] R. D. Parshad, E. Quansah, K. Black, R. K. Upadhyay, S. K. Tiwari, and N. Kumari, “Long time dynamics of a three-species food chain model with Allee effect in the top predator,” Comput. Math. with Appl., vol. 71, no. 2, 2016, doi: 10.1016/j.camwa.2015.12.015. [13] C. Ji, D. Jiang, and N. Shi, “A note on a predator-prey model with modified Leslie- Gower and Holling-type II schemes with stochastic perturbation,” J. Math. Anal. Appl., vol. 377, no. 1, 2011, doi: 10.1016/j.jmaa.2010.11.008. [14] X. Y. Meng, N. N. Qin, and H. F. Huo, “Dynamics analysis of a predator–prey system with harvesting prey and disease in prey species,” J. Biol. Dyn., vol. 12, no. 1, 2018, doi: 10.1080/17513758.2018.1454515. [15] T. K. Kar and S. K. Chattopadhyay, “A dynamic reaction model of a prey-predator system with stage-structure for predator,” Mod. Appl. Sci., vol. 4, no. 5, 2010, doi: 10.5539/mas.v4n5p183.