Effect of prey refuge on a harvested predator–prey model with generalized functional response

Commun Nonlinear Sci Numer Simulat 16 (2011) 1052–1059

Contents lists available at ScienceDirect

Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns

Effect of prey refuge on a harvested predator–prey model with generalized functional response q Youde Tao a,b, Xia Wang a,*, Xinyu Song a a b

College of Mathematics and Information Science, Xinyang Normal University, Xinyang 464000, Henan, PR China Beijing Institute of Information Control, Beijing 100037, PR China

a r t i c l e

i n f o

Article history: Received 10 September 2009 Received in revised form 21 May 2010 Accepted 24 May 2010 Available online 1 June 2010 Keywords: Predator–prey model Harvesting Prey refuge Stability Equilibrium

a b s t r a c t A predator–prey model with generalized response function incorporating a prey refuge and independent harvesting in each species are studied by using the analytical approach. A constant proportion of prey using refuges is considered. We will evaluate the effects with regard to the local stability of equilibria, the equilibrium density values and the long-term dynamics of the interacting populations. Some numerical simulations are carried out. Ó 2010 Elsevier B.V. All rights reserved.

1. Introduction The study of the consequences of hiding behavior of prey on the dynamics of predator–prey interactions can be recognized as a major issue in applied mathematics and theoretical ecology [1–5]. Some of the empirical and theoretical work have investigated the effects of prey refuges and drawn a conclusion that the refuges used by prey have a stabilizing effect on the considered interactions and prey extinction can be prevented by the addition of refuges [6–11]. Economic progress and ecological balance always have conflicting interests. Catering to the necessities and comforts of human beings invariably robs the ecological structure of the nature. Therefore, reasonable harvesting policies is indisputably one of the major and interesting problems from ecological and economical point of view. The exploitation of biological resources and harvest of population species are commonly practiced in fishery, forestry and wildlife management. Concerning the conservation for the long-term benefits of humanity, there is a wide-range of interest in the use of bioeconomic modelling to gain insight in the scientific management of renewable resources like fisheries and forestries. The problem of predator–prey interactions under constant rate of harvesting of either species or both species simultaneously have been studied by some authors. For example, Brauer and Soudak [12–15] studied a class of predator–prey models under constant rate of harvesting and under constant quota of harvesting of both species simultaneously. They showed how to classify the possibilities of the quantitative behavior of the solutions to locate the set of initial values in which the trajectories of the solutions approach to either an asymptotic stable equilibrium or an asymptotically stable limit cycle. q Supported by the National Natural Science Foundation of China (No. 10771179), the Scientific and Technological Project of Henan Province (No. 092102210070), the National Science Foundation of the Education Department of Henan Province (No. 2010B110021) and the Young Backbone Teacher Foundation of Xinyang Normal University. * Corresponding author. E-mail address: [email protected] (X. Wang).

1007-5704/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2010.05.026

Y. Tao et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 1052–1059

1053

For consistency with the previous work in this field, we consider the harvested predator–prey interactions incorporating refuges used by prey as a new ingredient. In this paper, we propose a generalized predator–prey model incorporating refuges used by prey and investigate its both stabilizing and destabilizing effects on the dynamics of the considered model. 2. The model and boundedness Recently, Dai and Tang [16] studied the following predator–prey model in which two ecological interacting species are harvested independently with constant rates:

(

_ Þ  puðxÞyðtÞ  u; xðtÞ ¼ rxðtÞð1  xðtÞ K _ yðtÞ ¼ yðtÞðd þ quðxÞÞ  h:

ð2:1Þ

They showed that system (2.1) possesses very complicated dynamics [16]. Motivated by the papers [8–16], we consider the following predator–prey model with a generalized functional response incorporating a constant proportion prey refuge and harvesting:

(

  _ xðtÞ ¼ rxðtÞ 1  xðtÞ  puðð1  mÞxÞyðtÞ  q1 E1 xðtÞ; K _ yðtÞ ¼ ðquðð1  mÞxÞ  dÞyðtÞ  q2 E2 yðtÞ;

ð2:2Þ

where x(t) and y(t) denote the prey and predator population at any time t respectively, and r, K, p, q, d, q1, q2 are all positive constants and have its biological meanings accordingly. r represents the intrinsic growth rate of the prey, K is the carrying capacity of the prey in the absence of predator and harvesting, p > 0 is the maximal per capita consumption rate of predator, q > 0 is the efficiency with which predators convert consumed prey into new predator, d > 0 is the death rate of the predator, E1 P 0, E2 P 0 denote the harvesting efforts for the prey and predator, respectively. q1E1x and q2E2y represent the catch of the respective species, where q1, q2 represent the catchability coefficients of the prey and predator, respectively. The model incorporates a refuge protecting mx of the prey, where m 2 [0, 1) is a constant. This leaves (1  m)x of the prey available to the predator. The term u(x) denotes the functional response of the predator and satisfies the following assumption:

uð0Þ ¼ 0; u0 ðxÞ > 0; x > 0:

ð2:3Þ

The solutions of model (2.2) represent the densities of the interacting populations and have their own realistically ecological meanings, that is to say, they must be positive and bounded. Therefore, we have the following theorem: Theorem 2.1. All the solutions of system (2.2) which start in R2þ are positive and uniformly bounded. Proof. The positively of the solutions of system (2.2) can be easily obtained according to their ecological meanings. In the following, we only show the boundedness of system (2.2). Denote

p VðtÞ ¼ xðtÞ þ yðtÞ; q then we have

  xðtÞ p _  q1 E1 xðtÞ  ðd þ q2 E2 ÞyðtÞ: VðtÞj ¼ rxðtÞ 1  ð2:2Þ K q

Now, for each a > 0, we get

K p _ VðtÞ þ aVðtÞ 6 ða þ r  q1 E1 Þ2  ðd þ q2 E2  aÞ: 4r q We can choose a > d + q2E2, then the right-hand side is positive. Thus, the right-hand side is bounded for all ðx; yÞ 2 R2þ . _ Therefore, we choose a k > 0 such that VðtÞ þ aVðtÞ < k. Applying the theory of differential inequality [17], we obtain

0 < Vðx; yÞ <

k

a

ð1  eat Þ þ Vðxð0Þ; yð0ÞÞeat !

k

a

ðt ! 1Þ:

So, we have that all the solutions of system (2.2) which start in R2þ are confined to the region X, where

X ¼ fðx; yÞ 2 R2þ jVðx; yÞ 6

k

a

þ e;

for

e > 0g:

This completes the proof. h

3. Main results For simplicity, we take the following change of variables:

xðtÞ ¼

xðtÞ ; 1m

yðtÞ ¼

yðtÞ : 1m

ð3:1Þ

1054

Y. Tao et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 1052–1059

ðtÞ as x(t), y(t)): And then system (2.2) takes the following form (still denote  xðtÞ; y

(

xðtÞ _ Þ  puðxÞyðtÞ  q1 E1 xðtÞ; xðtÞ ¼ rxðtÞð1  ð1mÞK

ð3:2Þ

_ yðtÞ ¼ ðquðxÞ  dÞyðtÞ  q2 E2 yðtÞ: 3.1. The existence of equilibria of (3.2)

We now study the existence of equilibrium of (3.2). Particularly, we are interested in the positive equilibrium of the system. By a simple calculation, we can list all possible equilibria: (i) The trivial equilibrium P0(0, 0); (ii) Equilibrium in the absence of predator P1(x1, y1), where

x1 ¼ ð1  mÞKð1 

q1 E1 Þ; r

y1 ¼ 0;

if 0 < q1E1 < r; (iii) The positive equilibrium P*(x*, y*), where

  d þ q2 E2 x ¼ u1 ; q

y ¼

qu1



dþq2 E2 q

2

pðd þ q2 E2 Þ

4r  q1 E1 

r u1



dþq2 E2 q

ð1  mÞK

3 5:

From the expression of y*, it is clear that a positive equilibrium P*(x*, y*) exists for the system (3.2) only if the harvesting rate satisfies:

0

0 < q1 E1 < r @1 

u1



dþq2 E2 q

ð1  mÞK

1

A;

ð3:3Þ

that is x* < x1. 3.2. Dynamic behavior of (3.2) In this subsection we discuss the stability properties of the equilibria P0, P1 and P*. The jacobian matrix of the system (3.2) about the equilibrium P0(0, 0) is given by:

 JðP0 Þ ¼

r  q1 E1

0

0

d  q2 E2

 :

Obviously, the eigenvalues for the steady state P0 are r  q1E1 and d  q2E2, where d  q2E2 is always negative. Now we see that r  q1E1 is positive or negative according as the prey biotechnical productivity (BTP) qr1 is greater or less than the effort E1. Thus if the prey (BTP) exceeds the effort, then P0 is a saddle point. Otherwise P0 is a locally asymptotically stable node. The Jacobian matrix of the system (3.2) about the equilibrium P1(x1, 0) is given by

 JðP1 Þ ¼

ðr  q1 E1 Þ

puðx1 Þ

0

quðx1 Þ  d  q2 E2

 :

ð3:4Þ

Eigenvalues of matrix (3.4) are k1 = (r  q1E1) < 0 (if x1 > 0 exists) and k2 = qu(x1)  d  q2E2. Thus, P1(x1, 0) is locally asymptotically stable if and only if

   q E1 qu ð1  mÞK 1  1  d  q2 E2 < 0: r That is





u ð1  mÞK 1 

q1 E1 r

 <

d þ q2 E2 ¼ uðx Þ: q

0

Since u (x) > 0 for x > 0, we have

0 r @1 

u1



dþq2 E2 q

ð1  mÞK

1 A < q1 E1 < r:

ð3:5Þ

Y. Tao et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 1052–1059

1055

Next, we consider the Jacobian matrix for the equilibrium P*(x*, y*):

JðP Þ ¼ where



a11

a12

a21

0

 ð3:6Þ

;

  ru1 dþqq2 E2 r  q1 E1 a11 ¼ A1  A2 ; d þ q2 E2 ðd þ q2 E2 Þð1  mÞK     d þ q2 E2 d þ q2 E2 Þu1 ; A1 ¼ d þ q2 E2  qu0 ðu1 q q      d þ q2 E2 d þ q2 E2 ; u1 A2 ¼ 2ðd þ q2 E2 Þ  qu0 u1 q q pðd þ q2 E2 Þ < 0; q     2  3 q2 u0 u1 dþqq2 E2 u1 dþqq2 E2 ru1 dþqq2 E2 4r  q1 E1  5 > 0: ¼ qu0 ðx Þy ¼ pðd þ q2 E2 Þ ð1  mÞK

a12 ¼ puðx Þ  a21

Clearly, detJ(P*) = a12a21 > 0, hence, the sign of the eigenvalues depends only on trJ(P*) = a11. Therefore, the equilibrium P* is locally asymptotically stable if and only if a11 < 0. Otherwise, the equilibrium P* is asymptotically unstable. According to A2 = d + q2E2 + A1 >A1, we can easily find the case A1 > 0, A2 < 0 is impossible, so we onsider three cases as follows: Case A. If A1 > 0, then A2 > 0. The positive equilibrium P*(x*, y*) of system (3.2) is locally asymptotically stable if and only if

 1 1 0 u1 dþqq2 E2 A 2 A < q 1 E 1 < r @1  A: 0 < r @1  ð1  mÞK A1 ð1  mÞK 0

u1



dþq2 E2 q



Otherwise, this equilibrium P* is asymptotically unstable. Case B. If A1 < 0 and A2 > 0, then the positive equilibrium P* of system (3.2) is always locally asymptotically stable, because the condition a11 < 0 holds in perpetuity. Case C. If A2 < 0, then A1 < 0. Hence the positive equilibrium P* of system (3.2) is locally asymptotically stable if and only if

0





1 A 2 A: 0 < q 1 E 1 < r @1  ð1  mÞK A1

u1

dþq2 E2 q

Therefore, we obtain the following theorems: Theorem 3.1. Assuming A1 > 0, we have (1) If

0 0 < q 1 E 1 < r @1 





1 A2 A ; ð1  mÞK A1

u1

dþq2 E2 q

then the prey and predator populations oscillate around the unique positive equilibrium point. (2) If

0 r @1 

u1



dþq2 E2 q



1

0

u1



dþq2 E2 q

1

A2 A A; < q 1 E 1 < r @1  ð1  mÞK A1 ð1  mÞK

then the two interacting populations tend to reach a locally asymptotically stable equilibrium at the first quadrant. (3) If

0 r @1 

u1



dþq2 E2 q

ð1  mÞK

1 A < q1 E1 < r;

then the predators go extinct while prey population reaches its maximum environmental carrying capacity.

1056

Y. Tao et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 1052–1059

Theorem 3.2. Assuming A2 < 0, we obtain: (1) If

0





1 A 2A ; 0 < q1 E1 < r @1  ð1  mÞK A1

u1

dþq2 E2 q

then the prey and predator populations tend to reach a locally asymptotically stable equilibrium at the first quadrant. (2) If

 1 1 0 1 dþq2 E2 u q A 2A A; < q1 E1 < r @1  r @1  ð1  mÞK A1 ð1  mÞK 0

u1



dþq2 E2 q



then the two interacting populations oscillate around the unique equilibrium P*(x*,y*) at the first quadrant. (3) If

0 r @1 

u1



dþq2 E2 q

1 A < q1 E1 < r;

ð1  mÞK

then the predators go extinct while prey population reaches its maximum environmental carrying capacity. Corollary 3.1. Assuming A1 > 0, we have (1) If

0
u1



dþq2 E2 q



A2 ; Kð1  q1rE1 Þ A1

then the prey and predator populations oscillate around the unique positive equilibrium P*(x*, y*) of system (3.2). (2) If

1

u1



dþq2 E2 q



u1



dþq2 E2 q



A2 ;
then the two interacting populations tend to reach a locally asymptotically stable equilibrium at the first quadrant. (3) If

1

u1



dþq2 E2 q



Kð1  q1rE1 Þ

< m < 1;

then the predators go extinct while the prey population reaches its maximum environmental carrying capacity. Corollary 3.2. Assuming A2 < 0, we obtain (1) If

0
u1



dþq2 E2 q



A2 ; Kð1  q1rE1 Þ A1

then the prey and predator populations tend to reach a locally asymptotically stable equilibrium P*(x*, y*) at the first quadrant. (2) If

  u1 dþqq2 E2 A2 1
u1



dþq2 E2 q



then the two interacting populations oscillate around the unique equilibrium P*(x*, y*) at the first quadrant.

Y. Tao et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 1052–1059

(3) If

1

u1



dþq2 E2 q

1057



Kð1  q1rE1 Þ

< m < 1;

then the predators go extinct while prey population reaches its maximum environmental carrying capacity.

4. Numerical simulation conclusions Now we illustrate our theoretical results with a numerical study. For the purpose of simulating the system, we let x uðxÞ ¼ 1þax , and p = 0.1, q = 0.001, r = 2.0, k = 600, a = 0.002, d = 0.00046, q1 = 0.2, q2 = 0.02, E1 = 1.0 in appropriate units.

Fig. 1. (1) Phase space trajectories corresponding to different initial levels, which shows that P*(19.168, 17.909) is a global attractor. (2) Both the prey and predator population converge to their equilibrium values.

Fig. 2. (1) Time series of prey in (3.2) with initial value (50, 30). (2) Time series of predator in (3.2) with initial value (50, 30).

1058

Y. Tao et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 1052–1059

Fig. 3. (1) Phase diagram of the limit cycle surrounding P*(10.683, 18.017). (2) There exist Hopf-bifurcating small amplitude periodic solutions.

Fig. 4. (1) Time series of prey in (3.2) with initial value (50, 50). (2) Time series of predator in (3.2) with initial value (50, 50).

If we take m = 0.15, E2 = 0.9, then the conditions of Corollaries 3.1 and 3.2 hold, P*(19.168, 17.909) exists and stable (see Figs. 1 and 2); If we take m = 0.01, E2 = 0.5, then the conditions of Corollaries 3.1 and 3.2 hold, system (3.2) admits exactly one limit cycle, we see that the corresponding equilibrium point P*(10.683, 18.017) is unstable. The phase diagram, as shown in Fig. 3 is a limit cycle (see Figs. 3 and 4). 5. Conclusions In this paper, we consider a harvested predator–prey model with a generalized functional responses incorporating a prey refuge. When a constant number of preys in refuges, we also find a clear stabilizing effect or an increase in stability of the interior equilibrium point of the considered model. From Theorem 3.1 and Corollary 3.1, the stability of the positive equilibrium P* is enhanced with the fraction of prey feeding in refuges when some conditions hold. Our results referring to increase in the equilibrium density of prey population and decrease in that of predators. Moreover, from Theorem 3.2 and

Y. Tao et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 1052–1059

1059

Corollary 3.2, the stability of the positive equilibrium P* is enhanced with the fraction of prey feeding in refuges when other conditions hold. Figs. 1–4 also show that the dependence of the dynamic behavior of system (3.2) on the prey refuge m. When m is small, there is a stable limit cycle surrounding the unique positive equilibrium (see Figs. 3 and 4), and when m is large enough, the limit cycle is broken and both the prey and predator population converge to their equilibrium values respectively (see Figs. 1 and 2), which means that if we change the value of m, it is possible to prevent the cyclic behavior of the predator–prey system and to drive it to a required stable state. Our results and numerical simulation also indicate that dynamic behavior of the model very much depends on the prey refuge parameter m. Hence, it is possible to control the system in such a way that the system approaches a required state, using prey refuge m as controls. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]

Holling CS. Some characteristics of simple types of predation and parasitism. Can Entomol 1959;91:385–98. Hassel MP, May RM. Stability in insect host-parasite models. J Anim Ecol 1973;42:693–726. Maynard Smith J. Models in ecology. Cambridge: Cambridge University; 1974. Hassell MP. The dynamics of arthropod predator–prey systems. Princeton, NJ: Princeton University; 1978. Hoy MA. Almonds (California). In: Helle W, Sabelis MW, editors. Spider mites: their biology, natural enemies and control. World crop pests, vol. 1B. Amsterdam: Elsevier; 1985. Sih A. Prey refuges and predator–prey stability. Theor Popul Biol 1987;31:1–12. Taylor RJ. Predation. New York: Chapman and Hall; 1984. Gonzalez-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–46. Krivan V. Effects of optimal antipredator behavior of prey on predator–prey dynamics: the role of refuges. Theor Popul Biol 1998;53:131–42. Kar TK. Stability analysis of a prey–predator model incorporating a prey refuge. Commun Nonlinear Sci Numer Simul 2005;10:681–91. Huang YJ, Chen FD, Li Z. Stability analysis of a prey–predator model with Holling type III response function incorporating a prey refuge. Appl Math Comput 2006;182:672–83. Brauer F, Soudack AC. Stability regions and transition phenomena for harvested predator–prey systems. J Math Biol 1979;7:319–37. Brauer F, Soudack AC. Stability regions in predator–prey systems with constant prey harvesting. J Math Biol 1979;8:55–71. Brauer F, Soudack AC. Constant-rate stocking of predator–prey systems. J Math Biol 1981;11:1–14. Brauer F, Soudack AC. Coexistence properties of some predator–prey systems under constant rate harvesting. J Math Biol 1981;12:101–14. Dai G, Tang M. Coexistence region and global dynamics of a harvested predator–prey system. SIAM J Appl Math 1998;58(1):193–210. Birkoff G, Rota GC. Ordinary differential equations. Boston: Ginn; 1982.