SQU Journal for Science, 2021, 26(1), 40-57 DOI:10.24200/squjs.vol26iss1pp40-57 Sultan Qaboos University 40 Predator-Prey Model with Refuge, Fear and Z- Control Ibrahim M. Elmojtaba, Kawkab Al-Amri and Qamar J. A. Khan* Department of Mathematics, College of Science, Sultan Qaboos University, P.O. Box 36, PC 123, Muscat, Sultanate of Oman. *E-mail: qjalil@squ.edu.om. ABSTRACT: In this paper, we consider a predator-prey model incorporating fear and refuge. Our results show that the predator-free equilibrium is globally asymptotically stable if the ratio between the death rate of predators and the conversion rate of prey into predator is greater than the value of prey in refuge at equilibrium. We also show that the co-existence equilibrium points are locally asymptotically stable if the value of the prey outside refuge is greater than half of the carrying capacity. Numerical simulations show that when the intensity of fear increases, the fraction of the prey inside refuge increases; however, it has no effect on the fraction of the prey outside refuge, in the long run. It is shown that the intensity of fear harms predator population size. Numerical simulations show that the application of Z- control will force the system to reach any desired state within a limited time, whether the desired state is a constant state or a periodic state. Our results show that when the refuge size is taken to be a non-constant function of the prey outside refuge, the systems change their dynamics. Namely, when it is a linear function or an exponential function, the system always reaches the predator-free equilibrium. However, when it is taken as a logistic equation, the system reaches the co-existence equilibrium after long term oscillations. Keywords: Predator-prey; Refuge; Fear; Z-control and Adaptive control. Z نموذج المفترس والفريسة مع ملجأ وخوف وتحكم خان ليل أحمدجال،كوكب العامري و قمر ىإبراهيم مجتب المفترس يكون مستقًرا بشكل يتضمن الخوف والملجأ. تظهر نتائجنا أن التوازن الخالي من الذي في هذه الورقة ، نعتبر نموذج المفترس والفريسة :صلخمال نظهر أيًضا أن مقارب عالميًا إذا كانت النسبة بين معدل وفاة المفترس ومعدل تحويل الفريسة إلى مفترس أكبر من قيمة الفريسة في الملجأ عند التوازن. القدرة االستيعابية. تظهر المحاكاة العددية أنه عندما نقطة توازن التعايش مستقرة محليًا بشكل مقارب إذا كانت قيمة الفريسة خارج الملجأ أكبر من نصف تبين أن شدة الخوف لها يليزداد جزء الفريسة داخل الملجأ ، ولكن ليس لها أي تأثير على جزء الفريسة خارج الملجأ ، على المدى الطو ،تزداد شدة الخوف سيجبر النظام على الوصول إلى أي حالة مرغوبة خالل فترة زمنية Zتطبيق التحكم تأثير سلبي على عدد الحيوانات المفترسة. تظهر المحاكاة العددية أن خارج الملجأ ، فإن محدودة ، سواء كانت حالة الرغبة حالة ثابتة أو حالة دورية. تظهر نتائجنا أنه عندما يتم اعتبار حجم الملجأ وظيفة غير ثابتة للفريسة ون دالة خطية أو دالة أسية ، فإن النظام يصل دائًما إلى التوازن الخالي من الحيوانات المفترسة. ومع ذلك ، عندما يتم . أي عندما تكديناميكياتهااألنظمة تغير أخذها كمعادلة لوجستية ، يصل النظام إلى توازن التعايش بعد التذبذبات طويلة المدى. التحكم التكيفي. ،التحكم -Z،يخاف أ،. لج،فريس ،المفترس:مفتاحيةالكلمات ال PREDATOR-PREY MODEL WITH REFUGE, FEAR AND Z-CONTROL 41 1. Introduction redator-induced stress in prey animals has not been much studied by physiologists or population ecologists. Canon [1] studied the concept of predator induced stress in a prey population. Population ecologists have focused their studies on predator induced stress effects on the birth rate of free-living prey populations and found that it affects the demography process of prey animals seriously. By experiments in laboratories and field studies it has been demonstrated that mere exposure of prey animals to predators affects the birth rate of prey, and this phenomenon of behavioral change of demography is called the “Ecology of Fear”, “Degree of Fear” or “Cost of Fear” by ecologists. Only a little work has been done on the subject of stress induced in a prey population due to predators. The concept of stress was limited to humans and it was thought that stress in a prey population is transitory and that it is not lifelong. Ecologists in the 1990’s [2-4] verified experimentally that fear of predators is a more powerful cause of demographical change in a prey population than direct killing, shortage of food, or parasitic infection. Fear of predators in prey persists even in the absence of predators and has long-lasting effects on their production of prey species. Sapolsky [5] explained the concept of stress in zebras due to fear of lions. Zanetteet al. [6] reported, after experimental verification, that the sparrow reduces offspring production by 40% just with intimidation by predators where direct killing is stopped by some means. Zanetteet al.[6], Eggers et al. [7], and Travers et al. [8] verified that female sparrows lay fewer eggs and, due to incubation disruption, fail to hatch eggs. Due to fear, they bring less food to their nests and, as a result, a greater proportion of their nestlings starve to death. Creel et al. [9], Creeland Christianson [10], and Creel et al. [11] reported that in the National Parks, USA, due to intimidation by wolves, elk pregnancy rates decline. Recently, Moncluset al. [12] examined an association between predator risk and birth rate of prey. Field studies have demonstrated that playback of predator sounds can affect the emotions of prey. Remage- Healey et al. [13] demonstrated that a playback sound of dolphins affects the emotions of gulf toad fish. Mateo [14] found that playback of the call of predators’ alerts squirrels and that they communicate predator risk to each other. The playback sound of predators increases the glucocorticoid level in prey and hence increases the fear or stress in the prey population. Wanget al. [15] studied how fear of a predator reduces the reproduction of prey animals and found that it could destabilize the system. A refuge is an area, such as island, where wild animals obtain protection from predation. In this protected territory the chances of being hunted by predator are reduced and these areas reduce the chances of extinction of prey species due to predatory killing. It is a natural phenomenon of prey species in an ecosystem to seek protection from predation (Cowlishaw [16], Sih [17]). Refuge habitats are of different types, such as burrows, trees, cliff faces, or dense vegetation (Clarke et al. [18], Dill and Houtman [19], Berger [20], Cassine[21]). Coral reefs provide refuge for prey fish (Friedlander and Martine [22], Sandinet al. [23]). Various control systems are used to prevent a species from drastic oscillations and avoid extinctions in an ecological system. The literature records various techniques for control in the multi-species Lotka-Volterra System. These controls are called adaptive control, back idea of control, impulsive control and applications of control theory to Lyapunov functions. One can refer to [24,25,26], where Z-control obtains the desired steady state quickly and prevents high amplitude oscillations. The Z-type control method is an error-based dynamic method, and in this method, it is certain that error function converges to zero. The error between desired outputs and actual outputs go to zero exponentially. There are two ways to apply Z-control in a predator-prey system. The first method is called direct control, where both prey and predator are controlled simultaneously to bring the population to a desired level. The second method is called indirect control, where either prey species or predator species is controlled through immigration, emigration or culling. The second species automatically comes to the desired level exponentially. Our model is similar to t h a t o f Wang et al. [27], and we examine this model by introducing predator fear. Prey only come out of refuge when they feel t h e r e i s less predation; otherwise, they go back to the protected area. We use two control measures to bring the model population of prey and predator to a desired level, and thus we can save it from becoming extinct. 2. Model formulation We studied a model where a prey species lives in two different habitats. One is called t h e refuge habitat where the prey species is saved from predation. It is assumed that all resources required for the growth of the prey species are available inside the refuge habitat and that their population grows logistically. It is also assumed that t h e prey species is fully protected from predation inside t h e refuge habitat. When pressure of predation fear i s released, then t he prey species moves to a second habitat outside the refuge and in this habitat t h e prey species can be killed by predators under the law of mass action. As predation fear increases in the prey species due to the presence of predators, the prey species migrates to the refuge habitat. In the absence of a prey species, predators die exponentially, because predators can only consume prey outside the refuge habitat. This predator- prey interaction is modelled by the following diagram and system of differential equations: P IBRAHIM M. ELMOJTABA ET AL. 42 Figure 1. Compartmental representation of the model. 𝑑𝑥1 𝑑𝑡 = 𝑎𝑥1 (1 − 𝑥1 𝑘 ) − 𝛼𝑥1 1+ 𝑒𝑦 + 𝛽𝑥2 𝑑𝑥2 𝑑𝑡 = 𝛼𝑥1 1+ 𝑒𝑦 − 𝛽𝑥2 − 𝑏𝑥2𝑦 𝑑𝑦 𝑑𝑡 = 𝑐𝑥2𝑦−𝑑𝑦 (1) All variables and parameters in model (1) are positive and defined below: 𝑥1 Prey density in the refuge habitat. 𝑥2 Prey density outside the refuge habitat. 𝑦 Abundance of predator species. 𝑑 Death rate of the predator. 𝑘 Carrying capacity of the prey in the refuge habitat. 𝑏 Feeding rate of the predator on the prey outside the refuge habitat. 𝑐 Conversion rate of prey to predator. 𝛼 Migration rate from the refuge habitat. 𝛽 Immigration rate into the refuge habitat. 𝑎 Birth rate of prey species inside refuge. 𝑒 The fear parameter. 3. Mathematical analysis of the model 3.1 Positivity of solutions Model (1) describes the dynamics of animal populations and therefore it is very important to prove that all quantities will remain positive for all time. We want to prove that all solutions of the model with positive initial data will remain positive for all time t>0. We can easily verify that 𝑑𝑥1 𝑑𝑡 | 𝑥1=0 = 𝛽𝑥2 ≥ 0 PREDATOR-PREY MODEL WITH REFUGE, FEAR AND Z-CONTROL 43 𝑑𝑥2 𝑑𝑡 | 𝑥2=0 = 𝛼𝑥1 1 + 𝑒𝑦 ≥ 0 𝑑𝑦 𝑑𝑡 | 𝑦=0 = 0 ≥ 0 . (2) Hence, all solutions will remain positive for all time. 3.2 Boundedness Proposition 1The trajectories of system (1), are bounded. Proof. Let 𝑤 = 𝑥1 +𝑥2 +𝑦. Take the time derivative along the solution of model (1) 𝑑𝑤 𝑑𝑡 = 𝑑𝑥1 𝑑𝑡 + 𝑑𝑥2 𝑑𝑡 + 𝑑𝑦 𝑑𝑡 = 𝑎𝑥1 (1− 𝑥1 𝑘 )− 𝑏𝑥2𝑦 + 𝑐𝑥2𝑦 − 𝑑𝑦 For any positive constant 𝑞we have: 𝑑𝑤 𝑑𝑡 +𝑞𝑤 ≤ 𝑎𝑥1 − (𝑏 − 𝑐)𝑥2𝑦 + 𝑞𝑥1 + 𝑞𝑥2 +𝑞𝑦 Where 𝑥1 ≤ 𝑘and 𝑥2 ≤ 𝑘. So 𝑑𝑤 𝑑𝑡 + 𝑞𝑤 ≤ (𝑎 + 2𝑞)𝑘 − [(𝑏 −𝑐) 𝑑 𝑐 − 𝑞]𝑦 Because𝑥2 > 𝑑 𝑐 . Now iff 𝑏𝑑 𝑐 > 𝑞 +𝑑, then 𝑑𝑤 𝑑𝑡 + 𝑞𝑤 ≤ (𝑎 +2𝑞)𝑘 Let (𝑎 +2𝑞)𝑘 = 𝐿.Therefore, we have 𝑤 ≤ 𝐿 𝑞 + 𝐴𝑒−𝑞𝑡 From which we can deduce that lim 𝑡→∞ sup𝑤 ≤ 𝐿 𝑞 Independently of the initial conditions. This completes the proof. Corollary 1: If 𝑏𝑑 𝑐 > 𝑞 +𝑑 > 0, then the region ∅ = {0 ≤ 𝑥1,𝑥2 ≤ 𝑘,0 ≤ 𝑦,𝑥1 +𝑥2 +𝑦 ≤ 𝐿 𝑞 } is an invariant region for model 1. Proof. This is a direct conclusion of Proposition 1 3.3 Equilibrium analysis Let𝑥1 ̅̅̅̅ ,𝑥2̅̅̅and �̅� be the equilibrium values of 𝑥1,𝑥2and𝑦.We find three biologically meaningful equilibrium points (i)𝐸0̅̅ ̅ = (0,0,0). The extinction of all populations, this equilibrium always exists. IBRAHIM M. ELMOJTABA ET AL. 44 (ii) 𝐸1̅̅ ̅ = (𝑘, 𝛼𝑘 𝛽 ,0).The prey species survive inside and outside the refuge habitat and the predator goes to extinction. (iii) 𝐸2̅̅ ̅ = (𝑥1̅̅̅ ,𝑥2̅̅̅ , �̅�). All populations survive. Note that at this equilibrium, and using the third equation of system (1), we get: 𝑥2̅̅̅ = 𝑑 𝑐 (3) From the second equation of system (1) at equilibrium we have 𝑥1̅̅̅ = 𝑥2̅̅̅ (1+ 𝑒�̅�)(𝛽 + 𝑏�̅�) 𝛼 . (4) Using the first equation of system (1) at equilibrium and substituting (4) we get 𝑓(�̅�) = 𝑎𝑥2̅̅̅ 𝑘𝛼 𝑒2𝑏2�̅�4 + 2𝑎𝑥2̅̅̅ 𝑒𝑏 𝑘𝛼 (𝑏 + 𝑒𝛽)�̅�3 + [ 𝑎𝑥2̅̅̅𝑏 2 𝑘𝛼 + 𝑎𝑥2̅̅̅𝑒 2𝛽2 𝑘𝛼 + 𝑎𝑒( 4𝑥2̅̅̅𝛽 2 𝑘𝛼 − 1)]�̅�2 +[𝑎𝑏( 2𝑥2̅̅̅𝛽 𝑘𝛼 − 1) + 𝑎𝑒𝛽( 2𝑥2̅̅̅𝛽 𝑘𝛼 − 1) +2𝑏]𝑦 ̅+𝑎𝛽( 𝑥2̅̅̅𝛽 𝑘𝛼 − 1) = 0 (5) If 1 2 < 𝑥2̅̅̅̅ 𝛽 𝑘𝛼 < 1 , then there will be only one positive root of �̅� of (5), because 𝑓(0) = 𝑎𝛽( 𝑥2̅̅̅̅ 𝛽 𝑘𝛼 − 1) < 0 , and lim�̅�→∞ 𝑓(�̅�) = ∞ , and then by the continuity of 𝑓(�̅�) and zero-point theorem,𝑓(�̅�) = 0 has one positive solution, so there will be a unique positive coexistence equilibrium. 3.4 Stability analysis In this subsection, we examine the stability of the system about the equilibrium points found in the previous subsection. (a) Stability analysis of the equilibrium (i): Consider a small perturbation about the equilibrium𝑥1 = 𝑥1̅̅̅+𝑢, 𝑥2 = 𝑥2̅̅̅+𝑣 and 𝑦 = �̅� +𝑤.Substitutingthese into the system (1), and neglecting products of small quantities, we obtain the stability matrix: ( 𝑎 − 𝛼 𝛽 0 𝛼 −𝛽 0 0 0 −𝑑 ) The corresponding characteristic equation is: −(𝜆 + 𝑑)(𝜆2 +𝜆(2+𝛽 −𝑎)−𝑎𝛽) = 0 (6) One of the values of 𝜆 is positive, so (0,0,0) is unstable and hence the both populations will never be extinct. (b) Stability analysis of the equilibrium (ii): Theorem 1. The predator free equilibrium𝐸1̅̅̅̅ = (𝑘, 𝛼𝑘 𝛽 ,0)of system (1) is locallyasymptotically stable if 𝑐𝛼𝑘 𝛽 < 𝑑, and unstable if 𝑐𝛼𝑘 𝛽 > 𝑑 .In fact, we can prove that 𝐸1̅̅̅̅ is globally asymptotically stable if 𝑐𝛼𝑘 𝛽 < 𝑑. PREDATOR-PREY MODEL WITH REFUGE, FEAR AND Z-CONTROL 45 Proof. Using the above mentioned, we obtain the stability matrix ( −𝑎 − 𝛼 𝛽 𝛼𝑏𝑘 𝛼 −𝛽 −𝛼𝑏𝑘 𝛽 −𝛼𝑏𝑘 0 0 𝑐𝛼𝑘 𝛽 − 𝑑 ) and the corresponding characteristic equation is ( 𝑐𝛼𝑘 𝛽 − 𝑑 − 𝜆)[𝜆2 +(𝑎 + 𝛼 +𝛽)𝜆+𝑎𝛽] (7) All roots will be negative if 𝑐𝛼𝑘 𝛽 < 𝑑, and hence this equilibrium will be stable. If 𝑔(𝑡) is a continuous and bounded function, then we define: 𝑔∞ ≜ lim 𝑡→∞ 𝑆𝑢𝑝 𝑔(𝑡), 𝑔∞ ≜ lim 𝑡→∞ 𝐼𝑛𝑓 𝑔(𝑡) For a system (1) with initial conditions 𝑥1 = 𝑥1(𝑡),𝑥2 = 𝑥2(𝑡)and = 𝑦(𝑡) , we have 0 ≤ 𝑥1∞ ≤ 𝑥1 ∞ ≤ ∞, 0 ≤ 𝑥2∞ ≤ 𝑥2 ∞ ≤ ∞, 0 ≤ 𝑦∞ ≤ 𝑦 ∞ ≤ ∞ . Using fluctuation lemma [28], we can say that there is a sequence {𝑡𝑛}, and when 𝑡𝑛 → ∞ we have 𝑥1(𝑡𝑛) → 𝑥1 ∞and 𝑑𝑥1(𝑡𝑛) 𝑡𝑛 → 0 as 𝑛 → ∞. On substituting 𝑡𝑛 into the third equation of system (1), We have 𝑑𝑦(𝑡𝑛) 𝑡𝑛 = (𝑐𝑥2(𝑡𝑛)−𝑑)𝑦(𝑡𝑛) Taking the limit on both sides lim 𝑛→∞ 𝑑𝑦(𝑡𝑛) 𝑡𝑛 = (𝑐 lim 𝑛→∞ 𝑥2(𝑡𝑛)−𝑑) lim 𝑛→∞ 𝑦(𝑡𝑛) Which gives 0 = (𝑐𝑥2 ∞ − 𝑑)𝑦∞ Therefore 𝑦∞ = 0or 𝑥2 ∞ = 𝑑 𝑐 . Adding the first and second equations of system (1), we obtain 𝑑𝑥1(𝑡𝑛) 𝑡𝑛 + 𝑑𝑥2(𝑡𝑛) 𝑡𝑛 ≤ 𝑎𝑥1(𝑡𝑛)(1− 𝑥1(𝑡𝑛) 𝑘 ) Taking the limit on both sides lim 𝑛→∞ 𝑑𝑥1(𝑡𝑛) 𝑡𝑛 + 𝑑𝑥2(𝑡𝑛) 𝑡𝑛 ≤ lim 𝑛→∞ 𝑎𝑥1(𝑡𝑛)(1− 𝑥1(𝑡𝑛) 𝑘 ) IBRAHIM M. ELMOJTABA ET AL. 46 0 ≤ 𝑎𝑥1 ∞ (1− 𝑥1 ∞ 𝑘 ) Therefore either 𝑥1 ∞ = 0or 0 < 𝑥1 ∞ ≤ 𝑘. According to the limit Theorem [29],we get lim𝑡→∞ 𝑥1(𝑡) = 𝑘. The second equation of system (1) yields 𝑑𝑥2(𝑡𝑛) 𝑡𝑛 ≤ 𝛼𝑥1(𝑡𝑛)−𝛽𝑥2(𝑡𝑛) Taking the limit on both sides lim 𝑛→∞ 𝑑𝑥2(𝑡𝑛) 𝑡𝑛 ≤ 𝛼 lim 𝑛→∞ 𝑥1(𝑡𝑛) −𝛽 lim 𝑛→∞ 𝑥2(𝑡𝑛) Which gives 𝑥2 ∞ ≤ 𝛼𝑘 𝛽 Therefore, and using the limit Theorem we have lim 𝑡→∞ 𝑥2 = 𝛼𝑘 𝛽 Again, using the second equation of system (1), we have 𝑑𝑥2(𝑡𝑛) 𝑡𝑛 ≤ 𝛼𝑥1(𝑡𝑛)−𝛽𝑥2(𝑡𝑛) Similarly, 0 < 𝛼𝑥1 ∞ −𝛽𝑥2 ∞; i.e. 𝛼𝑥1 ∞ > 𝛽𝑥2 ∞. If 𝑥1 ∞ = 0, then 𝑥2 ∞ < 0, which is not possible. Therefore, we take,𝑥1 ∞ = 𝑘. If we consider,𝑥2 ∞ = 𝑑 𝑐 , then 𝑑𝑥2(𝑡𝑛) 𝑡𝑛 ≤ 𝛼𝑥1(𝑡𝑛)−𝛽𝑥2(𝑡𝑛) Taking the limit on both sides 0 ≤ 𝛼𝑘−𝑘 𝑑 𝑐 Hence 𝑑 ≤ 𝑐𝛼𝑘 𝛽 This inequality is not true also because for stability of 𝐸1̅̅ ̅, we need 𝑑 > 𝑐𝛼𝑘 𝛽 . Therefore, the predator free equilibrium is globally asymptomatically stable if 𝑘 > 𝛽𝑑 𝑐𝛼 ; i.e. maximum prey population inside the refuge is greater than 𝛽𝑑 𝑐𝛼 . (a) Stability analysis of the equilibrium(iii): Theorem 2. The equilibrium point 𝐸2̅̅ ̅is locally asymptotically stable if𝑥1̅̅̅ > 𝑘 2 . Proof. The stability matrix of the system (1) around the equilibrium point 𝐸2̅̅ ̅ is ( 𝑝1 𝑝2 𝑝3 𝑞1 𝑞2 𝑞3 0 𝑐�̅� 0 ) PREDATOR-PREY MODEL WITH REFUGE, FEAR AND Z-CONTROL 47 The corresponding characteristic equation is 𝜆3 + (−𝑝1 − 𝑞2)𝜆 2 +(𝑝1𝑞2 −𝑝2𝑞1 −𝑐�̅�𝑞3)𝜆+(𝑐�̅�𝑞3𝑝1 −𝑐�̅�𝑝3𝑞1) = 0 (8) Where 𝑝1 = − 𝛼𝑥1̅̅̅ 𝑘 − 𝛽𝑥2̅̅̅ 𝑥1̅̅̅ < 0 𝑝2 = 𝛽 > 0 𝑝3 = 𝛼𝑒𝑥1̅̅̅ (1 + 𝑒�̅�)2 > 0 𝑞1 = 𝛼 1 + 𝑒�̅� > 0 (9) 𝑞2 = −𝛽 − 𝑏�̅� < 0 𝑞3 = − 𝛼𝑒𝑥1̅̅̅ (1 +𝑒�̅�)2 − 𝑏𝑥2̅̅̅ < 0 The equation (8) can be written as 𝜆3 + 𝑎1𝜆 2 +𝑎2𝜆+𝑎3 = 0 (10) Where 𝑎1 = (−𝑝1 −𝑞2),𝑎2 = (𝑝1𝑞2 −𝑝2𝑞1 −𝑐�̅�𝑞3)and 𝑎3 = 𝑐𝑦( 𝑞3𝑝1 −𝑝3𝑞1).The Routh-Hurwitz criteria for the third order system is given by: 𝑎1 > 0, 𝑎3 > 0and 𝑎1 𝑎2 > 𝑎3. Here 𝑎1 = (−𝑝1 −𝑞2) > 0 𝑎3 = 𝑐𝑦( 𝑞3𝑝1 −𝑝3𝑞1) = (−𝑎+ 2𝑎𝑥1̅̅̅ 𝑘 )( 𝛼𝑒𝑥1̅̅̅ (1 + 𝑒𝑦)2 + 𝑏𝑥2̅̅̅)+ 𝑏𝛼𝑥2̅̅̅ 1+ 𝑒𝑦 > 0 if 𝑥1̅̅̅ > 𝑘 2 . Now to show that 𝑎1 𝑎2 > 𝑎3; i.e. −(𝑝1 +𝑞2)(𝑝1𝑞2 −𝑝2𝑞1 −𝑐�̅�𝑞3) > 𝑐𝑦( 𝑞3𝑝1 −𝑝3𝑞1) Which, after simplification, gives (𝑝1𝑞2 −𝑝2𝑞1)(−𝑝1 −𝑞2) > 𝑐𝑦(𝑞2𝑞3 +𝑝3𝑞1) > 0 Now 𝑝1𝑞2 −𝑝2𝑞1 = ( 𝛼𝑥1̅̅̅ 𝑘 + 𝛽𝑥2̅̅̅ 𝑥1̅̅̅ )(𝛽 +𝑏𝑦)− 𝛽𝛼 1 +𝑒𝑦 = 𝛼𝑥1̅̅̅ 𝑘 (𝛽 + 𝑏�̅�) IBRAHIM M. ELMOJTABA ET AL. 48 So, clearly,𝑝1𝑞2 −𝑝2𝑞1 > 0. Note that −𝑝1 −𝑞2 > 0; therefore, 𝑎1 𝑎2 > 𝑎3. Hence the co-existing equilibrium𝐸2̅̅ ̅ = (𝑥1,̅̅̅̅ 𝑥2̅̅̅, �̅�) will be asymptotically stable if 𝑥1̅̅̅ > 𝑘 2 . Hence the co-existing equilibrium will be stable if therefuge prey population at equilibrium is higher than the carrying capacity of the prey in the refuge habitat. 4. Z-control To achieve predator population and prey population inside and outside the refuge to a desire level, direct Z- control strategy is used. The direct Z-control are functions that a r e incorporated in each equation of the system (1). This system is then described as follows 𝑑𝑥1 𝑑𝑡 = 𝑎𝑥1 (1 − 𝑥1 𝑘 )− 𝛼𝑥1 1+ 𝑒𝑦 + 𝛽𝑥2 − 𝑢1(𝑡)𝑥1 𝑑𝑥2 𝑑𝑡 = 𝛼𝑥1 1+𝑒𝑦 −𝛽𝑥2 −𝑏𝑥2𝑦 − 𝑢2(𝑡)𝑥2 (11) 𝑑𝑦 𝑑𝑡 = 𝑐𝑥2𝑦 − 𝑑𝑦 − 𝑢3(𝑡)𝑦 Then we define the error functions as 𝑥1 − 𝑥1𝑑 = 𝑒1 = 𝑢1, 𝑥2 −𝑥2𝑑 = 𝑒2 = 𝑢2, 𝑦 − 𝑦𝑑 = 𝑒3 = 𝑢3, where 𝑥1𝑑,𝑥2𝑑 and 𝑦𝑑 are desiredstates of prey inside the refuge, prey outside the refuge and the predator population respectively. These functions decay exponentially with time, i.e. 𝑒1,𝑒2 and 𝑒3 tends to zero. For achieving our purpose we adopt �̇�1 = −𝜆1𝑒1, �̇�2 = −𝜆2𝑒2and �̇�3 = −𝜆3𝑒3, with𝜆1,𝜆2,𝜆3 > 0.Now we have �̇�1 − �̇�1𝑑 = −𝜆1(𝑥1 −𝑥1𝑑) �̇�2 − �̇�2𝑑 = −𝜆2(𝑥2 − 𝑥2𝑑) (12) �̇� − 𝑦�̇� = −𝜆3(𝑦 −𝑦𝑑) Finally, with (11) and (12) we get the following three control functions 𝑢1 = 1 𝑥1 [𝑎𝑥1 (1− 𝑥1 𝑘 )− 𝛼𝑥1 1 + 𝑒𝑦 +𝛽𝑥2 −�̇�1𝑑 +𝜆1(𝑥1 − 𝑥1𝑑)] 𝑢2 = 1 𝑥2 [ 𝛼𝑥1 1 +𝑒𝑦 − 𝛽𝑥2 − 𝑏𝑥2𝑦 − �̇�2𝑑 +𝜆2(𝑥2 − 𝑥2𝑑)] (13) 𝑢2 = 1 𝑦 [𝑐𝑥2𝑦 − 𝑑𝑦 −𝑦�̇� +𝜆3(𝑦 −𝑦𝑑)] 5. Adaptive non-linear control We start by first non-dimensionalizing the system (1) by using the following transformations: 𝑥1 𝑘 = 𝑋1, 𝑥2 𝑘 = 𝑋2,𝑒𝑦 = 𝑌, 𝛼 𝑎 = 𝛼1, 𝛽 𝑎 = 𝛽1, 𝑏 𝑎 = 𝑏1,𝑎𝑡 = 𝜏, 𝑐𝑘 𝑎𝑒 = 𝑐1, 𝑑 𝑎 = 𝑑1.Now the system (1) takes the non-dimensional form 𝑑𝑋1 𝑑𝜏 = 𝑋1(1 − 𝑋1)− 𝛼1𝑋1 1 + 𝑌 + 𝛽1𝑋2 𝑑𝑋2 𝑑𝜏 = 𝛼1𝑋1 1+𝑌 − 𝛽1𝑋2 −𝑏1𝑋2𝑌 (14) PREDATOR-PREY MODEL WITH REFUGE, FEAR AND Z-CONTROL 49 𝑑𝑌 𝑑𝜏 = 𝑐1𝑋2𝑌 −𝑑1𝑌 We use non-linear feedback control for the system (14). This system can be represented as 𝑑𝑋1 𝑑𝜏 = 𝑋1(1 − 𝑋1)− 𝛼1𝑋1 1 +𝑌 + 𝛽1𝑋2 + 𝑢1 𝑑𝑋2 𝑑𝜏 = 𝛼1𝑋1 1+𝑌 − 𝛽1𝑋2 −𝑏1𝑋2𝑌 + 𝑢2 (15) 𝑑𝑌 𝑑𝜏 = 𝑐1𝑋2𝑌 −𝑑1𝑌 + 𝑢3 Where 𝑢1,𝑢2 and 𝑢3are adaptive nonlinear feedback control functions which will be the functions of 𝑥1,𝑥2 and y. If these feedback functions stabilize the system, then in a n infinitely long time state the variables converge to zero. Let 𝑒𝛼1 = 𝛼1 − �̂�1,𝑒𝛽1 = 𝛽1 − �̂�1,𝑒𝑏1 = 𝑏1 − �̂�1,𝑒𝑐1 = 𝑐1 − �̂�1,𝑒𝑑1 = 𝑑1 − �̂�1be unknown estimators, which give �̇�𝛼1 = −�̂�1̇, �̇�𝛽1 = −�̂�1 ̇ , �̇�𝑏1 = −�̂�1 ̇ , �̇�𝑐1 = −�̂�1̇, �̇�𝑑1 = −�̂�1 ̇ . To prove the global stability, we choose Lyapunov function as 𝑉(𝑋1,𝑋2,𝑌) = 1 2 𝑋1 2 + 1 2 𝑋2 2 + 1 2 𝑌2 + 1 2 𝑒𝛼1 2 + 1 2 𝑒𝛽1 2 + 1 2 𝑒𝑏1 2 + 1 2 𝑒𝑐1 2 + 1 2 𝑒𝑑1 2 Which has the derivative �̇� = 𝑋1 [𝑋1 −𝑋1 2 − 𝛼𝑋1 1+ 𝑌 + 𝛽𝑋2 + 𝑢1]+ 𝑋2 [ 𝛼𝑋1 1 +𝑌 − 𝛽𝑋2 − 𝑏1𝑋2𝑌 + 𝑢2]+ 𝑌[𝑐1𝑋2𝑌 −𝑑1𝑌 + 𝑢3]+ 𝑒𝛼1�̇�𝛼1 + 𝑒𝛽1�̇�𝛽1 + 𝑒𝑏1�̇�𝑏1 + 𝑒𝑐1�̇�𝑐1 +𝑒𝑑1�̇�𝑑1 Choosing adaptive non-linear controls 𝑢1 = −2𝑋1 + 𝑋1 2 + �̂�1𝑋1 1 +𝑌 − �̂�1𝑋2 𝑢2 = − �̂�1𝑋1 1+𝑌 +�̂�1𝑋2 + �̂�1 ̇ 𝑋2𝑌 − 𝑋2 (16) 𝑢3 = −�̂�1𝑋2𝑌 + �̂�1𝑌 − 𝑌 Using (16) in �̇�, we have �̇� = (−𝑋1 2 − 𝑋2 2 −𝑌2)+ 𝑒𝛼1 𝑋1𝑋2 1 + 𝑌 −𝑒𝛽1𝑋2 2 −𝑒𝑏1𝑋2𝑌 + 𝑒𝑐1𝑋2𝑌 − 𝑒𝑑1𝑌 2 −𝑒𝛼1�̂�1̇ −𝑒𝛽1�̂�1 ̇ − 𝑒𝑏1�̂�1 ̇ − 𝑒𝑐1�̂�1̇ − 𝑒𝑑1�̂�1 ̇ (17) = (−𝑋1 2 − 𝑋2 2 − 𝑌2)+ 𝑒𝛼1 (− 𝑋1 2 1 +𝑌 + 𝑋1𝑋2 1 +𝑌 − �̂�1̇) +𝑒𝛽1 (𝑋1𝑋2 − 𝑋2 2 − �̂�1 ̇ ) +𝑒𝑏1 (−𝑋2 2𝑌 − �̂�1 ̇ )+ 𝑒𝑐1(𝑋2𝑌 2 − �̂�1̇)+ 𝑒𝑑1 (−𝑌 2 − �̂�1 ̇ ) Considering parameter estimators �̂�1̇ = − 𝑋1 2 1 + 𝑌 + 𝑋1𝑋2 1 + 𝑌 +𝑒𝛼1 �̂�1 ̇ = 𝑋1𝑋2 − 𝑋2 2 + 𝑒𝛽1 (18) �̂�1 ̇ = −𝑋2 2𝑌 + 𝑒𝑏1 �̂�1̇ = 𝑋2𝑌 2 +𝑒𝑐1 IBRAHIM M. ELMOJTABA ET AL. 50 �̂�1 ̇ = −𝑌2 + 𝑒𝑑1 Using dynamics of unknown estimators (18) in (17) we will find �̇� = (−𝑋1 2 − 𝑋2 2 − 𝑌2)− 𝑒𝛼1 2 −𝑒𝛽1 2 − 𝑒𝑏1 2 − 𝑒𝑐1 2 − 𝑒𝑑1 2 Clearly the system will be globally stable, because �̇� < 0. 6. Numerical simulation In this section, we performed several numerical simulations for the system (1) to confirm our theoretical results and to acquire more knowledge about its dynamics and general behavior. The parameter values used are listed in the following table, some of them might be changed in order to study their effect: Table 1. Parameter values used for simulations. 6.1 Effect of adaptive control To study the effect of adaptive control on the system, we look at the stable co-existence equilibrium point. (Figure 2) shows that without the adaptive control the system took a long time to converge to this stable equilibrium point. However, with the use of adaptive control it is clear that the time to reach the equilibrium point is very short, and the system almost instantly started to reach this stable equilibrium point, as seen from (Figure 3). 6.2 Effect of the intensity of fear Parameter value 𝑎 0.07 𝛼 0.035 𝛽 0.0119 𝑘 0.8 𝑏 0.0112 𝑐 0.04 𝑑 0.07 𝑒 50 Figure 3.The effect of adaptive control on the convergence of the co-existence equilibrium point. The parameters are 𝑎 = 0.1,𝛼 = 0.035,𝛽 = 0.00119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.45,𝑑 = 0.07 and 𝑒 = 20. Figure 2. The convergence of the co-existence equilibrium point. The parameters are 𝑎 = 0.1,𝛼 = 0.035,𝛽 = 0.00119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.45,𝑑 = 0.07 and 𝑒 = 20. PREDATOR-PREY MODEL WITH REFUGE, FEAR AND Z-CONTROL 51 To study the effect of the intensity of fear, we simulate our model with parameter values taken from Table1, and the intensity of fear taken between 0.1 to 100. It is clear that when the intensity of fear increases, the fraction of the prey in the refuge increases. However, it has no impact on the fraction of prey outside the refuge, as shown from Figures 4-5. On the other hand, when the intensity of fear increases, the fraction of predators decreases as shown in ( Figure6), which dictates that fear is not in the favor of the predator. 6.3 Application of Z-type control with constant desired state Figures 7 - 9 show that with the help of z-type control all three populations, t h e prey in the refuge, the prey out of the refuge and the predators reaches the desired states as indicated. Clearly the time needed to reach the desired states is very short, and this is due to the power of z-type control, which takes the output to the desired state rapidly. Figures 1 0 - 1 2 show the control profile of all three populations, where ( Figure13) shows the error profiles of all three control profiles. Figure 4.The effect of the fear intensity on the prey in the refuge. The parameters are 𝑎 = 0.07,𝛼 = 0.035,𝛽 = 0.0119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.04 and 𝑑 = 0.07. Figure 5.The effect of fear intensity on the prey out the refuge. The parameters are𝑎 = 0.07,𝛼 = 0.035,𝛽 = 0.0119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.04 and 𝑑 = 0.07. Figure 6. The parameters are 𝑎 = 0.07,𝛼 = 0.035,𝛽 = 0.0119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.04 and 𝑑 = 0.07. IBRAHIM M. ELMOJTABA ET AL. 52 Figure 12. . The parameters are 𝑎 = 0.07,𝛼 = 0.035,𝛽 = 0.0119, 𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.04,𝑑 = 0.07and 𝑒 = 5. Figure 8. 𝑥2𝑑 = 1. The parameters are 𝑎 = 0.07,𝛼 = 0.035,𝛽 = 0.0119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.04,𝑑 = 0.07 and 𝑒 = 50. Figure 9.The convergence of the predator to its constant desired states 𝑦𝑑 = 0.5. The parameters are 𝑎 = 0.07,𝛼 = 0.035,𝛽 = 0.0119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.04,𝑑 = 0.07 and 𝑒 = 50. Figure 10. The parameters are 𝑎 = 0.07,𝛼 0.035,𝛽 = 0.0119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.04,𝑑 = 0.07 and 𝑒 = 50. Figure 11. The parameters are 𝑎 = 0.07,𝛼 0.035,𝛽 = 0.0119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.04,𝑑 = 0.07 and 𝑒 = 50. Figure 7.The convergence of the prey in the refuge to its constant desired states 𝑥1𝑑 = 1.5. The parameters are 𝑎 = 0.07,𝛼 = 0.035,𝛽 = 0.0119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.04,𝑑 = 0.07 and 𝑒 = 50. PREDATOR-PREY MODEL WITH REFUGE, FEAR AND Z-CONTROL 53 6.4 Application of Z-type control with periodic desired state Figures 14-16 show that Z-type control could be used to achieve a periodic desired state. From a biological point of view, it is very important to be able to reach a stable limit cycle instead of a constant equilibrium point, as periodic solutions (i.e. limit cycles) are of great interest for ecosystems and more generally for conservation biology. For the purpose of simulation, we consider the desired state to be of the form 𝑥1𝑑 = 𝑟1 + 𝜔1 cos 𝜋𝑡 100 𝑥2𝑑 = 𝑟2 + 𝜔2 cos 𝜋𝑡 200 𝑦𝑑 = 𝑟3 + 𝜔3 cos 𝜋𝑡 100 Figure 13. . The parameters are 𝑎 = 0.07,𝛼 = 0.035, 𝛽 = 0.0119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.04,𝑑 = 0.07and 𝑒 = 50. Figure 14.The convergence of the prey in the refuge to its periodic desired states. The parameters are 𝑎 = 0.07,𝛼 = 0.035,𝛽 = 0.0119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.04,𝑑 = 0.07,𝑒 = 50,𝑟1 = 1, 𝑟2 = 0.8, 𝑟3 = 0.8 ,𝜔1 = 0.5,𝜔2 = 0.25,𝜔3 = 0.8. Figure 15.The convergence of the prey out the refuge to its periodic desired states.The parameters are 𝑎 = 0.07,𝛼 = 0.035,𝛽 = 0.0119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.04,𝑑 = 0.07,𝑒 = 50,𝑟1 = 1, 𝑟2 = 0.8, 𝑟3 = 0.8 ,𝜔1 = 0.5,𝜔2 = 0.25,𝜔3 = 0.8. IBRAHIM M. ELMOJTABA ET AL. 54 Figure 17. . The parameters are 𝑎 = 0.07,𝛼 = 0.035,𝛽 = 0.0119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.04,𝑑 = 0.07,𝑒 = 50,𝑟1 = 1, 𝑟2 = 0.8, 𝑟3 = 0.8 ,𝜔1 = 0.5,𝜔2 = 0.25,𝜔3 = 0.8. Figure 16. . The parameters are 𝑎 = 0.07,𝛼 = 0.035,𝛽 = 0.0119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.04,𝑑 = 0.07,𝑒 = 50,𝑟1 = 1, 𝑟2 = 0.8,𝑟3 = 0.8 ,𝜔1 = 0.5,𝜔2 = 0.25,𝜔3 = 0.8. Figure 18. . The parameters are 𝑎 = 0.07,𝛼 = 0.035,𝛽 = 0.0119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.04,𝑑 = 0.07,𝑒 = 50,𝑟1 = 1, 𝑟2 = 0.8, 𝑟3 = 0.8 ,𝜔1 = 0.5,𝜔2 = 0.25,𝜔3 = 0.8. Figure 19. . The parameters are 𝑎 = 0.07,𝛼 = 0.035,𝛽 = 0.0119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.04,𝑑 = 0.07,𝑒 = 50,𝑟1 = 1,𝑟2 = 0.8,𝑟3 = 0.8 ,𝜔1 = 0.5,𝜔2 = 0.25,𝜔3 = 0.8. Figure 20. . The parameters are 𝑎 = 0.07,𝛼 = 0.035,𝛽 = 0.0119,𝑘 = 0.8,𝑏 = 0.0112,𝑐 = 0.04,𝑑 = 0.07,𝑒 = 5,𝑟1 = 1, 𝑟2 = 0.8, 𝑟3 = 0.8 ,𝜔1 = 0.5,𝜔2 = 0.25,𝜔3 = 0.8. PREDATOR-PREY MODEL WITH REFUGE, FEAR AND Z-CONTROL 55 6.5 Effect of different forms of 𝜷 It is more realistic than not to assume the refuge size (i.e. 𝛽) is not constant. Figures 21-23 show the fractions of prey in the refuge, the prey outside the refuge and the predator, where we took different forms of 𝛽. It is clear that when 𝛽 is taken as a linear function (i.e. 𝛽(𝑡) = (𝑎 +𝛽0)𝑥2) or as an exponential function of the prey outside the refuge (i.e. 𝛽(𝑡) = 𝑎𝑥2 exp(𝛽0) ), then both the prey outside the refuge and in the refuge reach their stable equilibrium and the predator goes to extinction. However, when 𝛽 is taken as a logistic function of the prey outside the refuge (i.e. 𝛽(𝑡) = 𝑎𝛽0𝑥2 (1 − 𝑥2 𝑘 )), then all three populations co-exist together after initial small oscillations. Note that 𝛽0 is a positive constant and represents the base-line value for 𝛽. 7. Conclusion All animals are threatened by predators and face the risk of predation. Prey populations change their behavior due to their fear of predators. In this paper we have studied the dynamics of predator-prey interaction where prey reside in two habitats, namely refuge and out of refuge. In refuge, the prey is safe from predatory killing and has sufficient resources to survive, and the population in the refuge grows logistically. Out of refuge, predators interact with the prey and may kill them. Prey live under the fear of predation, but when predator fear is diluted, prey come out of their refuge. On increasing predation fear, prey take shelter in the refuge. We obtain three biologically feasible equilibria and discuss their stability. The equilibrium free from prey and predator population will always exist and it is unstable i.e. prey, and predator populations will never be extinct. The equilibrium having zero predator population and non-zero prey population will always exist, and it will be globally stable if the maximum prey population inside the refuge is greater Figure 23. IBRAHIM M. ELMOJTABA ET AL. 56 than 𝛽𝑑 𝑐𝛼 . Otherwise it will be unstable. The co-existing equilibrium will be asymptotically stable if 𝑥1̅̅̅ is bigger than half of the carrying capacity, and otherwise will be unstable. To bring the population to the desired level and to protect it from extinction, we use z-control, where the population reaches the desired level in a short time. We performed a simulation where the desired state was limit cycles instead of an equilibrium point. Numerically it is shown that as the intensity of fear increases, the population of the prey in refuge increases, while the population of the predator decreases i.e. fear is not in the favor of predator populations. To make our study more ecologically realistic, we took different forms of refuge size (𝛽) i.e. linear, exponential, and logistic instead of constant. We observed that when refuge size 𝛽 is linear or exponential, the prey out of the refuge and in the refuge tend to attend their stable equilibrium and predators go to extinction. If we consider refuge size as a logistic function of the prey out of refuge, then after a little oscillation all three populations co-exist. The adaptive control inputs for asymptotic stability are obtained as non-linear feedback. We examined the stability of the system with and without control and noted that the system with control approaches stability faster than the system without control. Conflict of interest The authors declare no conflict of interest. Acknowledgement The authors would like to acknowledge anonymous reviewers, whose comments helped to increase the clarity and readability of the paper. References 1. Canon, W.B., Bodily changes in pain, hunger, fear and rage. D. Appleton and company, New York, NY,1915. 2. Peckarsky, B. L., Cowan, C.A., Penton, M.A. and Anderson, C., Sublethal consequences of stream-dwelling predatory stoneflies on mayfly growth and fecundity. Ecology, 1993, 74, 1836- 1846. 3. Krebs, C.J., Bouten, S., Boonstra, R., Sinclair, A.R.E., Smith, J.N.M., Dale, M.R.T. , Martin, K. and Turkington, R., Impact of food and predation on the snowshoe hare cycle, Science, 1995, 269, 1112-1115. 4. Schmitz, O.J., Beckerman, A.P., and O’Brien, K.M., Behaviorally mediated trophic cascades. Effects of predation risk on food web interactions. Ecology, 1997, 78, 1388-1399. 5. Spolsky, R.M., Why zebras don’t get ulcers, 3rdedn., Henry Holt and Company, New York, NY,2004. 6. Zanette, L.Y., White, A.F., Allen, M.C. and Clinchy, M., Perceived predation risk reduces the number of offspring songbirds produce peryear, Science, 2011, 334, 1398-1401. 7. Eggers, S., Griesser, M., Nystrand, M. and Ekman, J., Predation risk induces changes in nest-site selection and clutch size in the siberian jay. Proceeding of the Royal Society of London. B., Biological Sciences, 2006, 273, 701-706. 8. Travers, M., Clinchy, M., Zanette, L., Boonstra, R. and William, T.D., Indirect predator effects on clutch and the cost egg production. Ecology letters, 2010, 13, 980-988. 9. Creel, S., Christianson, D., Liley, S. and Winnie, J.A.Jr, Predation risk affects reproductive physiology and demography of elk. Science, 2007, 315,960. 10. Creel, S., Christianson, D., Relationships between direct predation and risk effects. Trends in Ecology and Evolution, 2008, 23, 194-201. 11. Creel, S., Christianson, D., and Winnie, J.A.Jr, A survey of the effects of wolf predation risk on pregnancy reates and calf recruitment in elk. Ecological Applications, 2012, 21, 2847-2853. 12. Monclu̇s, R.,Tiulim, J. and Blumstein, D.T., Older mothers fol low conservative strategies under predator pressure: the adaptive role of maternal glucocorticoids in yellow-bellied marmots. Hormones and Behavior, 2011, 60, 660-665. 13. Ramage-Healy, L, Nowacek, D.P., and Bass, A.H., Dolphin foraging sounds suppress calling and elevate stress hormone levels in a prey species, the Gulf toadfish, Journal of Experimental Biology, 2006, 209, 4441- 4451. 14. Mateo, J.M., Ecological and hormonal correlates of antipredator behavior in adult belding’s ground squirrels (Spermophilusbelding). Behavioral Ecology and Sociobiology, 2007, 62, 37-49. 15. Wang, X., Zanette, L., Zou, X., Modelling the fear effect in predator-prey interactions, Mathematical Biology, 2016, 73, 1179-1204. 16. Cowlishow, G., Refuge use and predation risk in a desert baboon population, Animal Behaviour, 1997, 54, 241-253. 17. Sih, A. To hide or not to hide? Refuge use in a fluctuating environment, Trends in Ecology and Evolution, PREDATOR-PREY MODEL WITH REFUGE, FEAR AND Z-CONTROL 57 1997, 12, 375-376. 18. Clarke, M.F., Da Silva, K.B., Lair, H., Pocklinton, R., Kramer, D.L., Mc laughlin, R.L., Site familiarity affects escape behavior of the eastern chipmunk, Tamius Striatus, Oikos, 1993, 66, 533-537. 19. Dill, L.M., Houtman, R., The influence of distance to refuge on flight initiation distance in the grey squirrel (Sciurus carolinensis), Canadian Journal of Zoology, 1989, 67,232-235. 20. Berger, J., Pregnancy incentives, predation constraints and habitat shifts experimental and field evidence for wild big horn sheep, Animal Behavior, 1991, 41, 61-77. 21. Cassine, M.H., Foraging under predation risk in the wild guinea pig caviaaperea. Oikos, 1991, 62, 20-24. 22. Friedlander, A.M., Martini, E.E., Contrasts in density, size, and biomass of reef fishes between the northwestern and the main hawaiian islands: the effects of fishing down apex predators, Marine Ecology Progress Series, 2002, 230, 253-264. 23. Sandin,S.A.,Smith,J.E.,DeMartine,E.E.,Dinsdale,E.A.,Donner,S.D.,Fiedlander, A.M., Konotchick, T., Malay, M., Maragos, J. E., Obura, D., Pantos, O., Paulay, G., Richie, M., Rohwer, F., Schroeder, R.E., Walsh, S., Jackson, J.B.C., Knowlton, N., Sala, E., Baselines and degradation of coral reef in the northern line islands, PLoS ONE 3, e1548,2008. 24. Zhange, Y., Yan, X., Liao, B., Zhang, Y., Ding, Y., Z-type control of populations for Lotka-Volterra model with exponential convergence, Mathematical Biosciences, 2016, 272,15-23. 25. Lacitignola, D., Diele, F., Marangi, C., Provenzale, A., On the dynamics of a generalized predator-prey system with Z-type control, Mathematical Biosciences, 2016, 280, 10-23. 26. Nadim, S.S., Samanta, S., Pal, N., Elmojtaba, I.M., Mukhopadhyay, I., Chattopadhyay, J., Impact of predator signals on the stability of a predator-prey system: AZ-Control Approach, Differential Equations and Dynamical Systems. https://doi.org/10.1007/s12591-018-0430-x, 2018. 27. Wang, H., Thanarajah, S. and Gaudreau, P., Refuge-mediated predator-prey dynamics and biomass pyramids, Mathematical Biosciences. 2018, 298, 29-45. 28. Hirsch, W.M., Hanisch, H., Gabriel, J.P., Differential equation models of some parasitic infections: Methods for the study of asymptotic behavior, Communication on pure and applied mathematics, 1985, 38, 733-753. 29. Castillo-Chavez, C., Thieme, H.R., Asymptotically autonomous epidemic models, Applied Mathematical Science, 1994. Received 8 December 2020 Accepted 28 February 2021 https://doi.org/10.1007/s12591-018-0430-x